A selection of H{\alpha} emitters at z = 2.1-2.5 using the Ks-band photometry of ZFOURGE

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.


INTRODUCTION
Since galaxies increase their masses through star formation, it is important to investigate properties of starforming 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 ∼ 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 ∼ 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 (∼ 0.3 dex) can be explained by a "selfregulated" 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 ) < 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α/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α/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α/UV 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α 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α 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α emission lines are sensitive to a diverse range of galaxies at z ∼ 2. Therefore we can obtain less-biased sample of star-forming galaxies with robust SFR measurements by using Hα emission line as a tracer.
There have been various studies to probe star-forming galaxies at z > 1 with Hα emission line by both spectroscopy and narrow-band imaging. For example, Reddy et al. (2018) measured Hα 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 narrowband imaging data of HiZELS (Geach et al. 2008;Sobral et al. 2009) to derive Hα 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 ) 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 broadband photometry, which allows the identification of lowmass galaxies over relatively large survey volumes. Indeed, several works have identified Hα emission line at z > 4 from color excesses in broad-band Spitzer/IRAC photometry (e.g., Smit et al. 2014Smit et al. , 2015Stark et al. 2013;Schaerer & de Barros 2009;Faisst et al. 2016Faisst et al. , 2019. At lower redshift, Saito et al. (2020) and Onodera et al. (2020) have extracted Hα 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α fluxes from composite SEDs consisting of stacked broad-band photometry SEDs.
In this work, we estimate Hα emission line fluxes of galaxies at z ∼ 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α fluxes extracted from broad-band photometry are compared to those from narrow-band imaging and spectroscopy. In Section 3, we investigate properties of Hα emitters at z = 2.1-2.5. Especially we discuss their Hα 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 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.

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), andCDFS (Giacconi et al. 2002). The survey area is ∼ 407 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 µ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 µ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 µ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 s -band imaging whose 5σ limiting magnitude is ≥ 25.5 across the fields, resulting in 80% completeness for galaxies with log (M/M ) ∼ 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 1 , J 2 , J 3 ) and two (H s , H l ), respectively. These filters efficiently sample the Balmer break (3636 A) and 4000Åbreak of galaxies at 1.5 z 3.5 and improve photo-z accuracy in this redshift range. Straatman et al. (2016) have derived photo-zs from SED fitting by EAZY (Brammer et al. 2008) and compared them with publicly available spectroscopic redshifts. They have found median σ z ∼ 0.013-0.022 in the three fields, where σ z ≡ 1.48 × |z phot − z spec |/(1 + z spec ).

Flux excess in the K s -band
Hα emission lines (rest-frame 6563Å) of galaxies at 2.1 < z < 2.5 fall into the K s -band of the ZFOURGE, while [OIII] (rest-frame 4959 & 5007Å) and Hβ (restframe 4861Å) into the two medium-bands (H s and H l ), respectively as shown in Figure 1.
To derive the emission line fluxes, we define "flux excess" (in units of 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: where f obs and f 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.

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) 1 . They have calculated the emission-line strengths following a method of Inoue (2011), where strengths of the UVoptical emission lines from Lyα to those at ∼ 1 µm relative to Hβ are calculated varying ionization parameters, metallicities, and hydrogen number densities using CLOUDY (Ferland et al. 1998). Hβ luminosity is estimated to be: where f esc and N lyc represent the escape fraction and the production rate of the Lyman continuum photons in a galaxy, respectively. While N lyc is calculated for each Bruzual & Charlot (2003) model template based on each age, there is no consensus on f esc and its redshift dependency. Salmon et al. (2015) have provided calculations for two extreme cases, f esc = 0 and f esc = 1. We adopt f 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 ∼ 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 ∼ 0.8 dex.

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 esc and N 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 lyc is doubled in the SED fitting. Each panel shows flux difference due to N lyc as a function of continuum flux (left) and a histogram (right). We find a typical difference of 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α fluxes. In our sample, a typical difference between the Hα fluxes from the flux excesses versus those derived from direct SED fitting is ∼ 40%.
2.5. Selection of Hα emitter at z = 2.1-2.5 We select our initial candidates of Hα emitters at z = 2.1-2.5 by requiring the flux excesses in the K s -band to be f excess ≥ 2σ Ks , where σ Ks is photometric error in the band. This selection results in 2005 Hα emitters (699, 632, and 674 in COSMOS, UDS, and CDFS, respectively).
Hα fluxes (in units of erg s −1 cm −2 ) of our sample are derived from the flux excesses as follows: where r Hα is a ratio of Hα to total strengths of all emission lines in the Ks-band, f excess is the flux excess, and ∆λ is a bandwidth of the Ks filter. The other major emission lines expected to fall in the Ks-band are [NII] and [SII]. The constant r Hα 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 Ks-band. A median value of r Hα is ∼0.79 over a metallicity range from 0.20Z to Z . It increases up to r Hα ∼ 0.9 for an extreme case of Z = 0.02Z , however, the impact on the final f Hα estimate is still insignificant. These values are consistent with the line ratios of star-forming galaxies at z ∼ 2 measured by Kewley et al. (2016) from spectroscopy. Uncertainties of the Hα fluxes are estimated using photometric errors in the Ks-band. The Hα 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 ) > 9.0.
2.6. Hα 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α 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α from galaxies at z = 2.2.
To check the consistency of our method with a conventional narrow-band emitter selection, we derive Hα fluxes from (K s − N B209) color excesses and compare them with those derived from the K 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α 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α fluxes derived using the narrow-band color excesses and the K s -band excesses is shown in Figure 4. We find that 46/63 (73%) of the Hα emitters have consistent values within a factor of 2 in both methods. The typical scatter around the oneto-one relation is 0.16 dex. Although some outliers from the one-to-one relation exist, there is no clear trend with stellar mass.

Long-slit spectroscopy
Some of the ZFOURGE galaxies have been observed by a spectroscopic survey, ZFIRE . We find that 93 Hα emitters in our sample have ZFIRE spectroscopic observations. Figure 5 shows a   comparison of the Hα fluxes between our method and ZFIRE, where a correction for different aperture sizes is applied. We find that 62/93 (67%) of the Hα 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 ∼ 2.7.

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

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α fluxes between the three methods already reported in the literature.
The scatter between the ZFIRE and 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 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α flux derivation compared to other methods.
In addition, we have also applied our method to H medium bands where pared to spectroscopic measurements of 3D-HST (??) they show good agreement with a typical scatter of < 0.17 dex (N. Chen et al., in preparation).

Corrections for dust extinction
When converting the Hα 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: where SFR(IR) and 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 µm scaled to match an observed 24 µm flux, while the UV luminosity is estimated from flux density at 1600 A 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 Hα 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 ∼ 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 ffactor is ∼ 0.44, galaxies with high Hα equivalent widths ( 100Å) have f -factors close to unity. As the redshifts of our Hα 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 µm band (S/N < 3), we use relation between visual extinction derived from SED fitting by FAST (A V,FAST ) and IRX extinction obtained from stacked MIPS 24µm fluxes (A V,IRX ) as outlined in the following. First, we divide our sample into A V,FAST bins regardless of the S/N of the MIPS detections. Median of 24 µm and 1600 A fluxes are calculated in each bin to obtain A 1600 via Equation 4, which is then converted to A V,IRX assuming the extinction curve of Calzetti et al. (2000). Figure 7 shows a correlation betweenA V,IRX and A V,FAST , where the best-fit relation is given by We use this relation to estimate visual extinctions for galaxies without IR detection, which is then converted to A Hα using the curve of Calzetti et al. (2000) again. We find that there is an offset between A V,IRX and A V,FAST . For galaxies with SFR(Hα) > 10, a typical value of A Hα based on the IRX is ∼ 0.21 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-β 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. We derive intrinsic Hα luminosities corrected for dust extinction and convert them to Hα SFRs using the calibration of Kennicutt (1998) with a correction to the Chabrier (2003) IMF: SFR(Hα) = 7.9 × 10 −42 L Hα × 10 −0.24 , where d L is the luminosity distance corresponding to the redshift of a galaxy.
3. PROPERTIES OF Hα 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α emitters on the restframe U V J diagram (e.g. Wuyts et al. 2007). The restframe UVJ data is provided in the ZFOURGE catalogs. It is clear that most of the Hα 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α emitters in the quiescent region. We note that a few low-mass galaxies with M * < 10 9 M (shown as bluer color in Figure 8) classified as quiescent have poor SED fitting so that their stellar masses and Hα excesses are Quiescent Star-forming Figure 8. Locations of the Hα emitters on the restframe U V J 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.
likely not reliable. However, most others have well-fitted SEDs and indeed a more conservative excess selection thresholds only reduces the Hα emitters in the quiescent region from 37 to 22 with 3σ Ks excess threshold. One possibility is that our emitter selection is sensitive enough to detect weak Hα emission lines from galaxies in quenching phases. Indeed, the quiescent Hα 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 U V J data. In any case, the existence of the quiescent Hα emitters does not affect the results below because they only represent 1.8% of the galaxies in Hα luminosity range of log L Hα < 43.2.

Deriving equivalent width
A total rest-frame equivalent width of all emission lines in the K s -band is calculated as follows: where W is the bandwidth of the K s -band filter. The Hα equivalent width is obtained by: using the r Hα ratio of Hα 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 As the K s photometry errors, σ Ks , do not vary significantly among the galaxies and f cont correlates with stellar masses, EW lim (Hα) is expected to show strong correlation with the stellar mass. We have estimated EW lim (Hα) for all the Hα emitters as shown in Figure 9.
In the data we find a correlation with stellar mass of the form: log EW lim (Hα) = −0.46 log (M/M ) + 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 ) > 9.0, which corresponds to the 80% mass completeness limit of ZFOURGE at z ∼ 2. A caveat on this conclusion is that we may still have missed emitters with relatively low equivalent widths.
In Figure 9 we also fit to only galaxies with log (M/M ) > 9.0 and find: log EW (Hα) = −0.24 log (M/M ) + 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 ) < 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α emitters with extremely high Hα equivalent widths (> 500Å). These values are much higher than the extrapolation of the relation at log (M/M ) > 9.0 to lower masses. Such galaxies are expected to have high sSFRs and we will further investigate them later. . Hα 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 ) > 9 is shown by the black line extrapolated toward the low-mass end. The blue line and intervals show the typical Hα equivalent width and its scatter reported by Reddy et al. (2018) for galaxies with log (M/M ) 9.3-10.7. Our method seems sensitive enough to detect most of the Hα emitters with log (M/M ) > 9.0, but only extreme ones can be detected below this limit. Figure 10 shows a correlation between the Hα SFR and estimated dust attenuation, which was derived two methods as described in Section 2.7. Clearly the galaxies with higher SFRs (i.e., more luminous in Hα) 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α SFRs and the amounts of the attenuation obtained from IRX relation in the sample is given by A Hα,IRX = 0.95 log SFR(Hα) − 0.51.

Attenuation by dust
A similar relation is obtained for UV+IR SFRs using only the galaxies detected in the Spitzer/MIPS 24 µm band with S/N>3. The above relation is for the A Hα derived from the IRX. The difference is consistent with  where A V,FAST shows systematically smaller values than A V,IRX .

Hα luminosity function
Before discussing the luminosity function itself, we evaluate completeness of the sample with respect to Hα 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 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α 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 s -band fluxes by adding the excesses to the stellar continuum fluxes. We then choose galaxies which have (1) K s -band flux larger than 5σ, and (2) the excess fluxes larger than 1σ, where σ is a photometric error, to be selected as Hα emitters. By calculating fractions of the emitters above the detec-tion limit in each Hα luminosity bin, we find that our sample is almost (∼ 96%) complete for galaxies with L Hα > 10 42.25 erg s −1 , where the Hα luminosities are corrected for attenuation by dust.
Number density of Hα 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 ∼ 5.51×10 5 Mpc 3 . Then the number density in a bin is obtained as follows: where L c is a central value of the bin. We count the number of emitters and divide it by the total volume (∆V ) and the bin width (∆ 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α luminosity function, which is given as where L * , α, and φ * are a characteristic luminosity, a faint-end slope, and a normalization, respectively. Figure 11 shows our Hα luminosity function with that of Sobral et al. (2013). The best-fit parameters are listed in Table 1. 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α luminosity function and that of Sobral et al. 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σ fitting uncertainties of our result are shown by red dashed curves.
(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α 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α profiles with AO-assisted imaging data. These results indicate at a fixed aperture size one potentially might miss Hα flux at large radii.
The flux measurements contained in the ZFOURGE catalog were originally performed with a 1.2 aperture on PSF-matched images and then corrected to total fluxes using growth curves from the K s -band (Straatman et al. 2016). In the HiZELS, Sobral et al. (2009) used a fixed 2 aperture for all galaxies so they might have missed some fractions of the Hα flux at radii > 2 .
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 where n is a Sérsic index, R e a half-light radius, and I e a surface intensity at that radius, I(R e ). Then b n is determined to satisfy the following equation: where Γ(n) and γ(n, x) are complete and incomplete gamma functions, respectively. We then create Hα 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: where β = 2.5. A FWHM of the seeing is given using α and β by In the case of the HiZELS, FWHM of the seeing is ∼ 1 (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α fluxes as a function of radius shown in Figure 12. A radius corresponding to the 2 aperture at z = 2.23 is ∼ 8 kpc and represented by the grey dashed vertical line in Figure 12. We can clearly see that at least ∼ 20% of the Hα 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 ∼ 40% of the Hα flux. We should note, however, that the measurements of radial profiles of the Hα fluxes in Nelson et al. (2016) are truncated at a radii of ∼ 10 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α flux.
A potential caveat of this analysis is the assumption we make in about the Hα sizes of our z =2.1-2.5 sample. Cumulative fraction [%] 9.0 < logM < 9.5 9.5 < logM < 10.0 10.0 < logM < 10.5 10.5 < logM < 11.0 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 ∼ 2 these average sizes are not appropriate, which will have implications for the analysis of missing Hα flux. Indeed, Wilman et al. (2020) have found Hα seems to follow continuum sizes, which means our Hα sizes we adopt for this analysis are likely overestimated. This would act to reduce the amount of missing Hα 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 ∼ 10%, which means more than 20% of Hα flux could have been missed even at z > 2. For comprehensive understanding on the missing flux in Hα, more observations are necessary to reveal structural properties of galaxies at z > 2. 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α luminosity and the estimated attenuation by dust. Sobral et al. (2013) adopted A Hα = 1 mag for all the galaxies, which is not supported by our data. They therefore have likely underestimated intrinsic luminosities of the brightest Hα 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α 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 brightend discrepancy between our luminosity function and that derived in Sobral et al. (2013) is caused by intrinsically extended Hα profiles and the larger A Hα in bright Hα emitters. This implies there are larger numbers of intrinsically bright Hα 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.

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 ∼ 13% of the total Hα 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 Hα = 0-10 45 erg s −1 , while the other is of L Hα = 10 41.6 -10 45 erg s −1 . These intervals are same as those used in Sobral et al. (2013), where L Hα = 10 41.6 erg s −1 corresponds to 0.01L * .
Given that the Hα emission line is a direct SFR tracer and that Hα emitters cover diverse populations of star- forming galaxies (Oteo et al. 2015), we expect a consistent Hα 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α emitter selection reflects the underlying star-forming galaxy population and that our correction for the dust extinction properly works to recover total Hα 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 s -band. The SFRD agreement is interesting given we find a larger number of bright Hα 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.

Star formation main sequence
We also investigate the relation between stellar mass and star formation rate of our Hα emitter sample. In the following analysis, all the AGNs identified by either X-ray, IR, or radio (38/2005) Figure 15. Hα 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 ) < 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α emitters with well-constrained SEDs classified by the U V J 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α emitter. Figure 15 shows a stellar masses-Hα SFRs plot of our sample with SFMS of previous studies (Speagle et al. 2014;Whitaker et al. 2014;Behroozi et al. 2013). The black line represents our best-fit SFMS, which is given by log SFR(Hα) = 0.66 log (M/M ) − 5.32.
In the fitting, we only used the galaxies in the range of 9 < log (M/M ) < 11 because the incompleteness becomes significant below this range. At log (M/M ) > 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α 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 con-trast, Hα selection is thought to reflect a more diverse galaxy population (Oteo et al. 2015).

Low-mass galaxies with high sSFR
At log (M/M ) > 9, we find that the scatter around the SFMS is ∼ 0.3 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 ) < 9, however, there are many lowmass galaxies with much higher SFRs than predicted from an extrapolation of the SFMS. We should note, of course, that only strong Hα emitters can be detected in the low-mass end due to the sensitivity limit of the K s -band of ZFOURGE. Nevertheless it is still interesting to investigate them further because the high sS-FRs 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 lowmass galaxies. It might be possible that the apparently high sSFRs are only caused by increased photometric errors in their Hα 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 ) = 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 ) < 9. Then we add fluctuations to the simulated SFRs based on the photometric errors in the K s -band (grey circles in Figure 16). 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 ( 0.5 dex) compared to the simulated scatter from photometric errors at log (M/M ) < 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.

Hα/UV ratio
If a galaxy has undergone an instantaneous starburst within the past 100 Myr, its luminosity ratio between Hα 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α  Figure 15 but individual observed galaxies are colored according to luminosity ratios between Hα 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 ) < 9 suggest that high photometric scatter is not the likely causing the high SFRs at low masses. luminosity compared to the observed UV. As a consequence, when a short starburst occurs, the Hα/UV ratio is boosted (relative to a constant SFH) during first 10 Myr. Subsequently the ratio is returns to a lower value 10-100 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α/UV ratios. Figure 17 shows a comparison of the Hα luminosities and the UV luminosities of our Hα emitters, where both have been corrected for the dust extinction. Grey line represents the expected relation when assuming a constant SFH (Faisst et al. 2019) and shaded area shows a ±0.3 dex region about the relation. Most of the galaxies with log (M/M ) < 9.0 have very high Hα/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α/UV ratios. These results suggest that the low-mass Hα emitters with high sSFRs found here likely recently experienced short period starbursts.

SFHs of the starburst galaxies
Here we bring the results together and explore the possibility that there are two starbust modes: one for lowmass 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α/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α/UV ratio in each stellar mass bin, as shown in Figure 18. It is clear that Hα/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α/UV ratios of the starburst galaxies, as Hα/UV ratios can be overestimated if UV-faint galaxies at fixed Hα luminosities and/or Hα-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αbright galaxies at fixed UV luminosities is plausible, so we explore it further in Figure 19. The Figure shows the Hα/UV ratios of the starburst galaxies in various stellar mass bins. As the starburst galaxies are, by definition, selected to have high Hα 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 im- pact on the Hα/UV ratios. Therefore the impact of the bias toward luminous Hα emitters on the Hα/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α/UV ratio implies a bursty SFH, suggesting that a starburst has occurred within past ∼ 10 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 lowmass galaxies, causing high Hα/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 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.

SUMMARY
We have identified 2005 Hα emitters at 2.1 < z < 2.5 from the ZFOURGE catalog. We take advantage of the deep ZFOURGE K s photometry and wide wavelength coverage and derive Hα flux estimates from observed K s photometry. We specifically measure the K 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. . 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α SFR to be classified as the starburst galaxies in our definition while dasheddotted horizontal lines is our Hα 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.
We have corrected Hα 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α depends on both Hα 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α 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α 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 aperture used in Sobral et al. (2013) would have missed ∼ 40% of the Hα fluxes for the luminous-end galaxies. Moreover, they have underestimated intrinsic Hα luminosities of the most luminous galaxies by at most 1 mag as they assumed A Hα = 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α 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α-luminous galaxies is small. Larger surveys and spectroscopic follow-ups will allow us to further investigate this rare Hα-bright population.
The SFMS from our sample agrees well with previous results above the stellar mass completeness limit, log (M/M ) > 9. Interestingly, we find there are many Hα 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α/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 shortertimescale 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α 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β+[OIII] of galaxies in the same redshift range as Hα emitters in this work, respectively.