The H$\alpha$ and [O III] $\lambda 5007$ Luminosity Functions of $1.2

Euclid and the Roman Space Telescope (Roman) will soon use grism spectroscopy to detect millions of galaxies via H$\alpha$ and [O III] $\lambda 5007$ emission. To better constrain the expected galaxy counts from these instruments, we use a vetted sample of 4,239 emission-line galaxies from the 3D-HST survey to measure the H$\alpha$ and [O III] $\lambda 5007$ luminosity functions between $1.16<z<1.90$; this sample is $\sim 4$ times larger than previous studies at this redshift. We find very good agreement with previous measurements for H$\alpha$, but for [O III], we predict a higher number of intermediate-luminosity galaxies than previous works. We find that for both lines, the characteristic luminosity, $\mathcal{L}_*$, increases monotonically with redshift, and use the H$\alpha$ luminosity function to calculate the epoch's cosmic star formation rate density. We find that H$\alpha$-visible galaxies account for $\sim 81\%$ of the epoch's total star formation rate, and this value changes very little over the $1.16<z<1.56$ redshift range. Finally, we derive the surface density of galaxies as a function of limiting flux and find that previous predictions for galaxy counts for the Euclid Wide Survey are unchanged, but there may be more [O III] galaxies in the Roman High Latitude Survey than previously estimated.


INTRODUCTION
Since its observational discovery by Riess et al. (1998) and Perlmutter et al. (1999), dark energy has been at the forefront of astronomical research. While the ΛCDM paradigm has been extremely successful in predicting the properties of the cosmological microwave background and reproducing observables such as large scale structure and cosmic abundances, there are many open questions remaining, especially about the nature of dark matter and the evolution of dark energy (e.g., Bull et al. 2016;Amendola et al. 2018, and references therein). In order to better constrain cosmological models, we must continue accruing more accurate and precise observational probes.
Such efforts include the use of large-scale spectroscopic surveys to measure baryonic acoustic oscillations and redshift space distortions throughout cosmic time.
NOAO Extremely Wide Field Infrared Imager (NEW-FIRM) Hα Survey (Ly et al. 2011), achieve the same efficiency but are usually restricted to a very limited slice of redshift space, and thus require many different filters to survey large volumes.
Normal galaxies have no strong emission lines between Lyα at 1216Å and [O II] λ3727, so beyond z ∼ 1, redshift surveys are most efficiently performed at near-infrared wavelengths, with lines such [O II] λ3727, Hβ, [O III] λ5007, and Hα. While it is possible to photometrically-select z 1 objects and then refine their redshifts with follow-up spectroscopy (e.g., Davis et al. 2003;Steidel et al. 2004;Lilly et al. 2007;DESI Collaboration et al. 2016), emission line galaxy (ELG) surveys to probe large swaths of cosmic time are most easily performed from space. GRAPES (Pirzkal et al. 2004), which used the ACS G800L grism of the Hubble Space Telescope, as well as WISP (Atek et al. 2010) and 3D-HST (Brammer et al. 2012;Momcheva et al. 2016), which used the HST /WFC3 G102 and G141 grisms, represent some of the first efforts to create space-based ELG samples. In this study, we use the 3D-HST sample described in Nagaraj et al. (2021b) and Nagaraj et al. (2021a), hereafter referred to as Paper I and Paper II, to further explore the emission-line properties of 1.2 z 1.9 galaxies.
Two near-future missions, Euclid (Laureijs et al. 2011(Laureijs et al. , 2012 and Roman (Green et al. 2012;Dressler et al. 2012;Spergel et al. 2015), will identify millions of galaxies at redshifts 0.7 z 2.7 using their Hα and [O III] λ5007 emission. Given the similarity in the spectroscopic survey designs of 3D-HST, the Euclid Deep Survey, and the Roman High Latitude Survey (see Figure 5 in Paper I for a visual depiction of the survey limits and observed 3D-HST fluxes), we can use 3D-HST as a pathfinder for these missions. In particular, by evaluating the Hα and [O III] λ5007 luminosity functions and biases with respect to other galaxy samples and dark matter distributions, we can estimate how many galaxies these programs will find and how well they will be able to measure cosmological parameters.
Other investigations have been limited by sample size (e.g., Shim et al. 2009) or spectral resolution (e.g., Pirzkal et al. 2013, where [O III] and Hβ are a blended feature). Particularly notable is the work of Colbert et al. (2013) and the follow-up study by Mehta et al. (2015), which used data from the HST /WFC3 Infrared Spectroscopic Parallel (WISP) survey (Atek et al. 2010) to create a sample of approximately 1,000 ELGs between 0.3 < z < 2.3.
Recently, Nagaraj et al. (2021b) carefully vetted a sample of 4350 1.2 < z < 1.9 emission-line galaxies which were originally identified on the 3D-HST grism frames by Momcheva et al. (2016). Here we use a subsample of 3,187 sources to determine the luminosity function and equivalent width distribution of [O III] λ5007 and Hα in this redshift range. Since the depth and resolution of the 3D-HST data are similar to the grism surveys planned for Euclid and Roman, we can use our measurements to refine the predictions for these studies, and improve our measurement of the amount of intermediate-redshift star formation that is occurring in emission-line galaxies. These data will also allow us to examine the relationship between [O III] λ5007 emission and star-formation rate at an epoch intermediate between the local universe, where emission from oxygen is mostly from [O II] λ3727, and the redshifts studied by Bowman et al. (2021), where [O III] λ5007 dominates.
This paper is the third in series that analyses 3D-HST sources at redshifts 1.2 < z < 1.9. Paper I showed empirical relations between stellar mass and various observational and physical properties, such as absolute magnitude in a rest-frame optical filter. Paper II focused on the relationships among stars, gas, and dust in galaxies with available mid-and far-IR data.

DATA AND SELECTION EFFECTS
In this section we describe the fluxes used for the Hα and [O III] λ5007 luminosity function calculations, including the critical issue of incompleteness for lowluminosity objects.

Data
As the data used to compute the Hα and [O III] λ5007 luminosity functions have been described in Papers I and II, we give only a brief overview here. From an initial list of 9341 1.2 z 1.9 candidates identified on WFC3/G141 grism frames by the 3D-HST survey (GO-11600, 12177, 12328;Brammer et al. 2012;Momcheva et al. 2016), we carefully identified a clean sample of 4350 ELGs brighter than the catalog's F125W (J) + F140W (JH) + F160W (H) magnitude limit of m J+JH+H = 26. Since the 625 arcmin 2 covered by the 3D-HST survey coincides with regions of the Cosmic Assembly Nearinfrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011;Koekemoer et al. 2011), the fields have a wealth of photometric observations, and these data helped fuel the results of Paper I and Paper II.
Paper I describes the process we used to identify active galactic nuclei (AGN) masquerading as normal galaxies in our sample. First, we used X-ray matching from the deep surveys of the CANDELS fields (especially GOODS-S), to identify 72 AGN in our ELG sample. Specifically, we used a cross-correlation search radius of 1 to identify possible X-ray counterparts to our ELGs. Any galaxy with an X-ray luminosity greater than 10 42 erg s −1 in the 2 − 10 keV band was considered an AGN and eliminated from the analysis. Stacking of the remaining ELGs then found X-rays levels consistent with those expected from simple star formation, suggesting that the vast majority of objects remaining in our sample are normal galaxies with no obvious AGN activity.
While X-rays are very effective for AGN identification, they fail when the gas column densities are too high (N H 5 − 50 × 10 23 cm −2 ), which tends to correspond to highly dusty, Compton-thick systems (e.g., Brandt & Alexander 2015, and references therein). To find such objects, we used the IRAC AGN selection criteria described by Donley et al. (2012); the method identified 50 AGN candidates (with 11 being previously excluded via their X-ray luminosity). This whittled down our ELG sample to 4239 objects. Bowman et al. (2019) ran nearly the same procedure on 3D-HST galaxies at 1.9 ≤ z ≤ 2.35. After removing AGN, they were left with a sample of 1964 ELGs with trustworthy redshifts and clean spectra. These data are used in §4.3, where we combine the datasets to examine the evolution of the [O III] λ5007 luminosity function over the full redshift range from 1.16 ≤ z ≤ 2.35.

Completeness and Selection Effects
In this work, we use the Hα and [O III] fluxes derived from the 3D-HST grism spectra by Momcheva et al. (2016) A larger concern is our inability to separate Hα from the bracketing forbidden lines of [N II] λλ6548, 6584. While corrections for [N II] exist in the literature (e.g., the procedure and references outlined in Price et al. 2014), gas-phase metallicities defined through strong line indicators are subject to degeneracies, inaccuracies, and other issues (e.g., Kewley & Ellison 2008). Furthermore, the ionization balance between the first and second excited states of oxygen changes dramatically, from predominantly O + at low redshift to O ++ at z ∼ 2 (e.g., Kewley et al. 2019, and references therein). Since N + has an ionization potential only slightly less than O + , we can expect a similar trend in our data. This shift in ionization balance means that the [N II]/Hα ratio must also evolve with redshift.
Another concern associated with our understanding of the data involves the issue of dust attenuation. The wavelength dependence of attenuation has a variety of shapes in different galaxies (see Shivaei et al. 2020;Salim & Narayanan 2020, and references therein). An accurate modeling of this behavior is difficult, especially given its degeneracy with the other parameters of the complex stellar populations that make up a galaxy's SED. Even at the wavelength of Hα (6563Å), dust attenuation represents an uncertainty. Given the issues associated with [N II] and dust, we do not attempt to disentangle the flux of Hα from that of [N II]; instead, we present a luminosity function for the combined lines. We discuss this further in §4.2 where we calculate the star formation rate density between 1.2 z 1.9.
An extremely important issue for any luminosity function analysis is completeness. Following Bowman et al. (2021), we parameterize the completeness of the [O III] λ5007 and Hα detections as a function of flux, f , using a modified Fleming et al. (1995) function, In the equations, F c (f ) is the completeness, α F describes how quickly the completeness drops off as a function of flux, f 50 is the flux where the recovery fraction of objects is 50%, and f 10 is the flux at which the sample is 10% complete. As the Fleming function is completely described by f 50 and α F , f 10 is not an independent parameter but a quantity directly derived from Equation 1. An important point to note is that for G141 grism data, the behavior of the completeness curve is virtually independent of wavelength (e.g., Zeimann et al. 2014;Bowman et al. 2021). There are two ways to estimate the parameters α F and f 50 . The first method, which was used by Bowman et al. (2021) in their analysis of the [O III] λ5007 luminosity function of 1.90 < z < 2.35 grism-selected galaxies, is to use the data themselves and simultaneously fit the galaxy luminosity function ( §3) and each field's completeness parameters to the observed distribution of emission-line luminosities. Such a procedure is complex, since, even if α F is assumed to be the same across all five CANDELS fields, the calculation still involves at least 9 separate variables (a minimum of three for the luminosity function, one value of f 50 for each field, and α F ).
Alternatively, it is possible to fit the completeness curves separately from the luminosity function by assuming the intrinsic flux distribution of faint emission lines is a power law. This approach is reasonable, given the expected nature of the faint galaxy number counts, and was the method employed by Bowman et al. (2019) in their study of the physical properties of z ∼ 2 grismselected galaxies.
Our experiments show that both methods yield similar results for the galaxy luminosity function. Therefore, to avoid any degeneracies associated with high-dimensional fits and to simplify the analysis, we chose to decouple the question of completeness from the luminosity function calculation and adopt the values of f 50 and α F found by Bowman et al. (2021). These parameters, which were derived for 1.90 ≤ z ≤ 2.35 [O III] λ5007 emission on the same 3D-HST frames used here, were found using the simultaneous fitting technique described above. Since completeness should only be a function of line flux, and not depend on the specific line being observed, the Bowman et al. (2021) values should be equally applicable to our survey.
There is one potential caveat to this last assertion. As detailed by Momcheva et al. (2016), the 1σ sensitivity limit of the 3D-HST survey depends not only on the strength of an emission line, but also the angular size of the emitting source. Since the present study focuses on galaxies at lower redshift than those measured by Bowman et al. (2021), this difference has the potential to cause a shift in the survey's 50% completeness limit. However, based on values from the 3D-HST catalog (Momcheva et al. 2016), the mean size of 1.16 ≤ z ≤ 1.90 emission-line galaxies is only 1.15 times that of the Bowman et al. (2021) systems. In comparison, the size spread amongst the galaxies in the sample is a factor of ∼ 4. For this reason, we do not account for this size difference in our analysis. We list the completeness parameters in Table 1 and show the completeness curves in Figure 1. For the purposes of the luminosity functions presented in this paper, we only consider flux measurements above the 50% levels given in Table 1. The inclusion of objects much below this limit induces strong effects on the luminosity function by amplifying the uncertainties in the completeness curves. Conversely, if the flux limit is too strict, the sample of galaxies becomes too small for any reliable measurement of the function's low-luminosity end. Our 50% cutoff results in 2947 [O III] and 1892 Hα measurements being used for our luminosity function calculations. For the analysis of evolution in the [O III] λ5007 luminosity function, this number can be incremented using the 1.9 ≤ z ≤ 2.35 ELGs found by Bowman et al. (2021), which were identified in the same manner as the galaxies used here. This results in a sample of 4519 [O III] λ5007 ELGs above the 50% completeness limit. The maximum likelihood estimator (MLE) is commonly used in astronomy to measure and fit the parameterized variables of a luminosity function. Simplifications to the MLE integral have also led to computationally efficient procedures for creating discrete luminosity functions derived from techniques such as the 1/V max (Schmidt 1968(Schmidt , 1970Huchra & Sargent 1973;Avni & Bahcall 1980) and the C − (Lynden-Bell 1971) methods.  Table 1.
Here we do both, and derive both a digital representation of the number of galaxies versus emission-line (log) flux, and a convenient analytical representation of the luminosity function.
We assume Poisson statistics hold, and define the likelihood of a luminosity function (scripted as P to avoid confusion with luminosity) following the derivation by Ciardullo et al. (2013). For computational convenience, we use L ≡ log L rather than the linear luminosity; in both cases, the math is identical, although the units are different. With the equations in this section, we are able to both measure the galaxies' discrete luminosity function and fit the data to a parameterization of our choosing (i.e., the Schechter 1976, function).
From Ciardullo et al. (2013), the likelihood of observing any luminosity function, φ is given by (4) where N is the number of galaxy luminosities included in the sample. Note that there is a distinction between the true luminosity function φ(L, z), which we model as a Schechter (1976) function with parameters α, L * , and log φ * , i.e., and the observed luminosity function, φ , which is affected by incompleteness and measurement error.
Following Marshall et al. (1983) and Marshall (1985), we define a function Ω(L, z) that takes into account both the flux completeness and the limits of the survey area. If f is the flux, which can be thought of as a function of luminosity and redshift, and p(f,n) is the completeness function in a given directionn, then we have Equation 6 below, where d 2n is the integral over the unit sphere.
In our case, we consider p(f,n) constant over any given field used in the analysis. If Ω i is the effective survey area of field i, then Ω(L, z) can be approximated as Assuming that the true luminosity is isotropic, we can simply connect the observed and true luminosity functions through Going back to Equation 4, z 1 and z 2 represent the survey redshift limits. The minimum luminosity, L min (z), is somewhat arbitrary and is simply taken to be a luminosity lower than what is observable. In practice, the upper limit of the integral is also fixed at a value at, or above, that which no galaxies are expected to exist. Finally, dV /dz is the differential volume element. Assuming a spatially flat universe, this is simply where d A (z) represents the comoving angular diameter distance and H(z) is the Hubble parameter. Given this definition, the true luminosity function represents the number of galaxies per dex (log luminosity units) per comoving volume element. Moreover, given the nature of the 3D-HST survey, we assume that the flux completeness function is constant within each CAN-DELS field, so in a given field, p(f,n) ≡ p(f ). As described in §2, we use the modified Fleming et al. (1995) completeness curve F c (f ) for p(f ).
The likelihood analysis presented here has not included measurement errors in the luminosity (i.e., heteroscedasticity). We follow the convention adopted in Mehta et al. (2015), in which the luminosity errors are assumed to be normal with mean L i and standard deviation σ i . In that case, we can replace Equation 4 with The drawback of this approach is that it is computationally expensive, as the normal distribution has a possibly different mean and standard deviation for every flux measurement. The results given in §4 are therefore based on Equation 4. We find that including luminosity errors has little to no effect at the bright end of the [O III] luminosity function fit but does slightly lessen the low-luminosity slope α. In other words, the observed faint-end slope is affected by Eddington (1913) bias. Bowman et al. (2021) provides a more detailed analysis of the effects of photometric uncertainties on the shape of the luminosity function.
If we assume that the luminosity function remains unchanged over the redshift interval of interest, the MLE for discrete points becomes much simpler to compute. Let us parameterize the luminosity function as a sum of discrete Dirac-delta functions at M different luminosities, i.e., Then, using Equations 4 and 10, the likelihood can be expressed as Setting the derivative of the likelihood to zero, the MLE solution for the luminosity function φ i at L i becomes This is a slight generalization of the historic 1/V max method (Schmidt 1968(Schmidt , 1970Huchra & Sargent 1973;Avni & Bahcall 1980) for any given completeness curve p(f,n). We use this formalism, which we will call the V eff method, as a benchmark for our more sophisticated MCMC approaches.

Computing the Luminosity Function
Our simplest method of computing the luminosity function is the V eff approach introduced in §3.1. We calculate φ i (Equation 12) at every source luminosity in our sample and collect the results into luminosity (or log luminosity) bins. Wide bins allow for larger numbers of sources per bin and are thus more reliable, but do not convey as much information given their coarseness. We find that dividing the measurements into ∼ 50 bins of equal size in log luminosity space works quite well for balancing the number of sources per bin against the complexity of the results.
The errors on our fitted parameters are generated via a bootstrap analysis. We take the true V eff -method result using the original set of N (L i , φ i ) values, and then generate B bootstrap samples, in which the N values of L i and φ i are generated randomly with replacement. The sample variance is taken to be the error on the luminosity function measurement, where θ is the binned luminosity function measured at a particular interval, and B is typically set to B = 100. We include all fluxes down to the 50% completeness limit ( §2) in the V eff method. For the computation of the Schechter (1976) function parameters, our MCMC code uses uniform priors on α, log φ * , and L * with bounds [−3, 1], [−8, 5], and [40,45], respectively, while the completeness parameters are fixed at the values given in Table 1. We define the relative likelihood of a solution either through Equation 4 (no observational errors) or Equation 9 (with observational errors) and employ the emcee package (Foreman-Mackey et al. 2013) to explore the parameter space. We compute both a static (non-evolving) luminosity function and one that evolves over time.
To explore time evolution in the luminosity function, we use a method based on the technique described in Leja et al. (2020). We let both log φ * and L * be quadratic functions of redshift. However, rather than fitting the coefficients of the quadratic, whose priors are difficult to physically motivate , we use the values of log φ * and L * at three specific redshifts (z 1 = 1.20, z 2 = 1.76, and z 3 = 2.32 for [O III] λ5007 and z 1 = 1.18, z 2 = 1.36, and z 3 = 1.54 for Hα) to define the quadratic formulation of L * (z) and log φ * (z). As in our static luminosity-function calculation, the prior on L i * is uniform over the range [40,45] and the prior on log φ i * uniform on [−8, 5].
For the analysis, we fix the completeness parameters and also constrain α to be constant over time. Doing so avoids exacerbating degeneracies between α and the other two Schechter parameters. (In fact, in the case of [O III], we find that leaving α as a free parameter over the large redshift range 1.16 ≤ z ≤ 2.35 leads to failures in the MCMC fits. For the [O III] λ5007 line, we therefore we fix α = −1.5.) For the reader's convenience, in Table 2 we have listed all the parameters used in the non-evolving and redshiftvarying luminosity function calculations, as well as the priors applied in the MCMC code.
In all cases, we employ 100 walkers and 1000 steps, resulting in 100,000 MCMC realizations. In other words, 100 points in the parameter space are randomly selected as initial states, and the MCMC algorithm takes 1000 steps from each initial state to find regions of higher likelihood. There is no "best-fit" solution, but we find that given these generous numbers for walkers and steps, the solutions generally do converge, suggesting that the true best-fit is closely approached.

Cosmic Variance
One source of uncertainty in the normalization of our emission-line luminosity functions is cosmic variance. To estimate the expected amplitude of this effect, we use the cosmic variance calculator 1 of Trenti & Stiavelli (2008), which employs both the extended Press-Schechter formalism (Press & Schechter 1974) and numerical simulations to compute the expected variance in any pencil-beam region of the sky. Following Colbert et al. (2013) and Bowman et al. (2021), we compute the cosmic variance by taking the result for each of the five disconnected CANDELS fields and then dividing the average of these estimates by √ 5. The results of this calculation show that for a nonevolving luminosity function over the redshift range 1.16 ≤ z ≤ 1.56, the cosmic variance expected for our Hα-emission luminosity function is ∼ 6.7%, while that for [O III] λ5007 galaxies between 1.16 ≤ z ≤ 1.90, this number is ∼ 5.3%. For the redshift-varying case, the process of measuring the cosmic variance is not so straightforward, as we have modeled cosmic evolution using a quadratic equation, represented using the values of φ * and L * at three redshifts. However, if we divide the surveyed redshift range into three bins, we can estimate the effect of cosmic variance on each bin. We find that for both Hα and [O III] λ5007, the variance should slightly increase with redshift, with the uncertainties being roughly 12% and 9%, respectively.
We include the effects of cosmic variance in our calculations for the expected galaxy counts (Tables 3 and  4) as well as §4.4) and the star formation rate density 1 https://www.ph.unimelb.edu.au/∼mtrenti/cvc/ ( §4.2). Given the sizes of our galaxy samples and the volumes of space being surveyed ( 0.7 and 1.4 × 10 6 Mpc 3 for the Hα and [O III] λ5007 studies, respectively), the effects of cosmic variance should be small in the case of the static luminosity function, but non-negligible for our analysis of cosmic evolution.

RESULTS
In this section, we present the Hα + [N II] and [O III] λ5007 luminosity functions of 3D-HST ELGs, along with the Hα-based cosmic star formation rate density contained in the emission-line galaxies. Tables 3, 4, and 5 present the overall results for our sample.
In Tables 3 and 4, we list our best-fit static luminosity functions along with those of Shim et al. (2009), Colbert et al. (2013, Sobral et al. (2013), Khostovan et al. (2015), and Sobral et al. (2015). As a cautionary statement, the luminosity functions being given in the tables are not directly comparable, as detailed in the columns labeled "Notes". Moreover, because of the well-known degeneracies between the three Schechter (1976) parameters, our values of α, L, and φ i are not necessarily in agreement with those of the previous studies. Nevertheless, we find that the overall form of our Hα luminosity function is compatible with the luminosity functions derived by Colbert et al. (2013) and Sobral et al. (2013), but somewhat distinct from the literature measurements for [O III] λ5007.
In Table 5, we give the parameters of our best-fit Hα and [O III] λ5007 redshift-evolving luminosity functions. For the latter, we also extend the redshift range to z = 2.35 using the measurements of Bowman et al. (2021), since their galaxy sample was defined in exactly the same manner as our dataset.
We note that in a grism survey, the true survey area is difficult to calculate as contamination from overlapping spectra and edge effects reduce the number of objects included in analyses. Bowman et al. (2021) studied this censoring by masking out those regions of the 3D-HST survey where emission-line detections are compromised, and fitting a luminosity function using only those galaxies in the unmasked areas. They found that the effective survey area of 3D-HST is ∼ 85% of the total survey area; this is consistent with the correction applied by Colbert et al. (2013) and the estimate made by Ciardullo et al. (2014). In this work, we reduce the quoted area of the 3D-HST survey by 15% to approximate the effects of overlapping spectra and edge-losses.

Hα + [N II] Luminosity Function
As mentioned in §1, the low resolution of the G141 grism and the morphological broadening associated with    In Figure 2, the top plot shows the Hα + [N II] luminosity function with the three fitted parameters from Equation 5: α, L * , and log φ * . Fits from two hundred MCMC iterations are shown in red, and the median fit in displayed gray. The small spread in the solutions suggests a stable and well-characterized result.
The plots in the lower triangle show the marginal posterior distributions for the three parameters, along with 2-D contour plots of the MCMC chains. The strong correlations in the contour plots confirm that the three parameters are not independent. This is expected, as the Schechter parameters are not orthogonal variables. However, the ubiquity of the function in the literature makes continued efforts in such a parameterization worthy.
While a non-evolving luminosity function for Hα + [N II] λ6584 is valuable, the non-negligible redshift range of the data, 1.16 < z < 1.56 enables us to study the luminosity function's evolution. As mentioned in §3.2, we did this by employing an approach based on Leja et al. (2020), in which both L * and log φ * are assigned to be quadratic functions of redshift. Because of the dearth of galaxies at very low luminosities, we forced α to be the same across all redshifts; this avoids the issue of degeneracies between α, L * , and log φ * seen in Figure  2.
We find that the redshift-evolving luminosity function solution is quite similar to the static case, with a statistically indistinguishable low-luminosity slope α and similar L * and log φ * values. All relevant quantities are presented in Table 5. Figure 3 shows the best-fit luminosity functions with redshift indicated by the color. The use of quadratic functions for L * and log φ * allows for smooth, differentiable evolution of the luminosity function. The red points on top of the color curve show the V eff result for the entire redshift range. The points are in greater agreement with the redshift-constant luminosity function but are still within the bounds of the redshift-varying function showed here.
We also show how our results compare to the literature. We include comparisons to Shim et al. (2009), who found 80 Hα emitters at redshifts 0.7 < z < 1.9 through a Hubble-NICMOS grism survey; Colbert et al. (2013), who analyzed a sample of 517 Hα emitters at 0.9 < z < 1.5 found through the Hubble WISP program; and Sobral et al. (2013), who obtained a dataset of 515 Hα emitters at z = 1.47 through the HiZELS narrow-band imaging. To ensure an apples-to-apples comparison, the literature luminosity functions have been modified to reflect the inclusion of [N II] in our data. Shim et al. (2009) and Colbert et al. (2013) assume F Hα = 0.71F Hα+ [NII] , while Sobral et al. (2013) uses a formula for correcting [N II] based on equivalent width, which gives an average correction of 25% over all their data. We undo these corrections, along with the 1 mag internal extinction correction applied by Sobral et al. (2013). Figure 3 shows the results. Our sample of Hα-emitting galaxies is ∼ 4 times larger than that of any previous study. But from the figure, it is clear that the Hα luminosity functions of Colbert et al. (2013) and Sobral et al. (2013) are in good agreement with our work. The only serious discrepancy is with the curve produced by Shim et al. (2009), but since that measurement was based on only 80 objects, this difference is not a concern.
From both Table 5 and Figure 3, we see that L * increases with redshift. Given that we are highly complete at all redshifts above L > 41.8, this finding must reflect a physical difference as we go back in cosmic time: there are more high-luminosity ELGs at earlier epochs. From Figure 3, we also find fewer low-luminosity objects at higher redshifts, but this result is more subject to completeness issues, and our result may not be robust.
We explore the effects of incompleteness in Figure 4. For our main results, we have fixed the minimum completeness fraction at 50%, i.e., we have excluded from the analysis all objects with monochromatic fluxes fainter than the 50% limit shown in Figure 1. But in Figure 4, we perform an experiment in which we vary the minimum threshold from 1% to 80% and observe the effects on the best-fit luminosity function. As shown in the left panels of the figure, the best-fit value for α decreases (gets steeper), L * increases, and φ * decreases when the minimum completeness fraction increases.
Nevertheless, from the right panel, we see that the overall luminosity function does not vary significantly in the luminosity range we are able to observe. This is because the three Schechter parameters are correlated: various sets of values can lead to the same overall luminosity function. We see from this experiment that the decision of which flux measurements to include changes the parameterization of the luminosity function, but not the overall shape curve.  Colbert et al. (2013). The MCMC results are also compatible with our V eff data point (red triangles). In regard to evolution, we find that the knee of the luminosity function, L * , increases with redshift while the normalization factor log φ * decreases, although this may not be due to true redshift evolution (see text for more details). In any case, the evolution is not particularly strong in our [1.16, 1.56] redshift range.

Cosmic Star Formation Rate Density
The evolution of the star formation rate density (SFRD) of the universe is an important indicator of galaxy growth and evolution. As reviewed by Madau & Dickinson (2014), we have the general picture that star formation peaked around z ∼ 2, an era dubbed as "cosmic noon," and has been declining at ∼ 0.1 dex per Gyr ever since. However, while this outline is known, further measurements of the SFRD, especially for certain populations of galaxies such as ELGs, are useful for improving our knowledge of the evolution of star formation and quantifying how selection effects propagate into this understanding. In this section, we use the Hα luminosity function to calculate the SFRD between 1.16 ≤ z ≤ 1.56.

Conversion between Hα luminosity and Star Formation Rate
Hα luminosity is often used as a direct proxy for very recent star formation (under 10 Myr; e.g., Kennicutt & Evans 2012). As such, applying a conversion formula from Hα luminosity to star formation rate (SFR) is a common process. Typically, the conversion from  show that α decreases (gets steeper), L * increases, and φ * decreases as the completeness limit for the analysis increases. This experiment demonstrates the importance of completeness in determining the Schechter (1976) function parameters. The right panel shows the luminosity functions for the Schechter parameters in the left panels. The lines are colored according to the minimum completeness fraction. We can see that the effects of completeness on the shape of the luminosity function are much less pronounced than that variables used in the parameterization.
Nevertheless, the aforementioned relation is calibrated in the local universe, therefore raising the concern of its application to higher redshifts. For example, at lower metallicities, main sequence stars tend to be bluer and hotter, changing the amount of ionizing radiation emitted by massive stars, and thus the calibration.
Furthermore, the Hα-SFR calibration summarized by Kennicutt & Evans (2012) applies to dust-corrected Hα brightness, and at higher redshifts, the details of dust attenuation constitute a major source of uncertainty (e.g., Bouwens et al. 2012;Nagaraj et al. 2021a). This error propagates directly into the SFR conversion, and is then compounded by the fact that our Hα data is contaminated by [N II] λλ6548, 6584, while the SFR calibration of Kennicutt & Evans (2012) is for Hα only. Thus, the applicability of the local conversion of Hα luminosity to SFR is unclear.
We can address this issue directly by calculating our own relation between SFR and observed Hα + [N II] luminosity. In Paper I, we discussed the SED fitting procedure of the entire sample of 1.16 ≤ z ≤ 1.90 3D-HST emission-line galaxies. We used the Bayesian MCMC SED code MCSED  to estimate the physical properties of our ELG sample, assuming a Kroupa (2001) initial mass function (IMF), a single-valued (but free) stellar metallicity, a binned star formation history, with the SFR in each bin a free pa-rameter, and a dust attenuation parameterization from Noll et al. (2009) and Kriek & Conroy (2013).
MCSED uses simple stellar population (SSP) SEDs from the Flexible Stellar Population Synthesis (FSPS) library  with Padova isochrones (Bertelli et al. 1994;Girardi et al. 2000;Marigo et al. 2008). Nebular emission is treated via interpolation in tables of CLOUDY models ) computed by Byler et al. (2017). Gas-phase metallicity (Z gas ) is set equal to the stellar metallicity in MCSED, and we fixed the ionization parameter at log U = −2.5 based on the high [O III]/Hβ ratios observed for our sources.
Our SED fits are based primarily on the 3D-HST photometry collected by Skelton et al. (2014), which consists of 147 filter bands distributed over the five CANDELS fields and covering the wavelength range from 3,000Å to 80,000Å (observed frame). In Paper I, we merged these data with photometry from Swift and GALEX, which extended the wavelength coverage for over 400 sources to ∼ 2, 000Å (observed-frame). Also, mid-and far-IR measurements from Spitzer and Herschel were added for over 600 sources in Paper II, but to maintain consistency within the sample, we do not include these dust-sensitive wavelengths here. In addition, MCSED is able to employ emission line fluxes in its SED-solution, making it an ideal analysis tool for grism-based surveys. We included Hα, Hβ, and [O III] λ5007 fluxes for galaxies whenever available.
Our MCSED-based SFR estimates are fairly robust. We find that the mean uncertainty on our SFRs is 0.24 dex, with a standard deviation of 0.08 dex. Moreover, in Paper I, we examined how the basic assumptions underlying our SED fits affected the derived properties of the 3D-HST galaxies. For example, we found that changing the ionization parameter, modifying the weights assigned to the emission line fluxes, or fixing galaxy metallicity all led to statistically indistinguishable distributions for the galaxies' SFRs. In other words, ionization parameter, metallicity, and the inclusion of emission line fluxes do not noticeably affect the systematics of an SFR measurement. In addition, Bowman et al. (2020) showed that for a large sample of 3D-HST galaxies, the choice of dust attenuation curve does not strongly affect the SFR estimate.
On the other hand, we do find that the details of a galaxy's assumed star formation history (SFH) do affect our SFR estimates. Our fits are based on a "nonparametric" SFH, i.e., one in which the SFR of each epoch in a galaxy's history is fit independently of the other epochs. Such fits have been proven to reduce a bias in SFR measurements that is introduced by the use of parameterized SFHs (e.g., Conroy 2013; Leja et al. 2017Leja et al. , 2019Bowman et al. 2020), though the magnitude of this bias reduction depends on the prior used in each age bin (Leja et al. 2019). Our quoted SFRs represent the star formation rate during the most recent age bin, i.e., over the last 100 Myr of cosmic time.
Given that we are trying to calibrate the SFR-Hα relation, the lack of influence of the Hα flux on the SFR measurement suggests that correlations found between SFR and Hα are not artificially induced by the method of measuring SFR. Furthermore, the benefit of using MCSED SFRs is that we have a way of connecting the observed Hα + [N II] fluxes with SFRs that takes into account dust attenuation and contamination by [N II]. This reduces the bias and uncertainty introduced by applying single values for the dust and [N II] corrections.
To calculate the mean SFR-Hα relation, we use a procedure very similar to that employed by Bowman et al. (2021) in their analysis of the [O III] λ5007 luminosity function. We first divide the Hα + [N II] luminosities into 25 bins with equal numbers of objects each interval. In each bin, we adopt the mean linear SFR and linear Hα luminosity as representative values. We then take the logarithm of those values and fit a line. We chose this process because if we average the measurements in logarithmic space, we would underestimate the total SFR of the population. As suggested by Feigelson & Babu (1992), we use the orthogonal distance regression technique to fit the line, since there is non-negligible heteroscedasticity associated with the measurement errors for both luminosity and SFR.
To estimate errors in our solution, we bootstrap the Hα + [N II] luminosity measurements and repeat the procedure above 1000 times. This error estimation process is the same as that described in §3.2 for the V eff method.
We present the results for our sample in Figure 5. The data are plotted as blue dots, while the linear mean binned values are shown as amber diamonds. Finally, the best-fit line is shown in red, with the 1σ uncertainty displayed via the shaded region (nearly too small to be visible). The correlation between apparent Hα + [N II] luminosity and SFR is quite clear.
The best-fit relation is given by Equation 14. The line has a slope consistent with 1, suggesting that the relation between SFR and uncorrected Hα + [N II] luminosity is very close to linear. To be consistent with the SFRD compilation of Madau & Dickinson (2014), we have divided the (linear) SFR by 0.67 as done in their work; this converts an SFR based on the Kroupa (2001) IMF to one based on the Salpeter (1955) To be complete, we also calculate the SFRD using the Kennicutt & Evans (2012) Sobral et al. (2013) use a formula based on equivalent width and find a median correction of 25%. The 29% correction, which we now adopt, is equivalent to a −0.15 dex shift in log luminosity.
To test the appropriateness of this 29% correction (log [N II]/Hα = −0.54), we determined how [N II]/Hα should vary as a function of population age and metallicity using the CLOUDY  nebular emission tables created by Byler et al. (2017), under the assumption of log U = −2.5 (the same value used in our MCSED fits). We show the result in Figure 6. According to the CLOUDY lookup tables, the correction is valid only for most galaxies with log(Z/Z ) −0.2. For galaxies with log(Z/Z ) −0.2, the true relative strength of [N II] is lower than the correction. Figure 5. Calibration between SFR and observed (not corrected for dust) Hα + [N II] luminosity for the sample used in this paper. We show the data as blue dots, the logarithm of the mean linear SFR and luminosity in 25 bins as amber diamonds, and the best-fit line to the binned mean values in red, which has a slope nearly identical to one (a linear relationship in linear space). There is a clear correlation in the data. Figure 6. CLOUDY prediction for log [N II]/Hα as a function of metallicity and age when log U = −2.5. While log [N II]/Hα is nearly independent of age, it is highly dependent on metallicity. For galaxies with log(Z/Z ) −0.7, [N II] emission is less than 1% of the Hα emission.
Stellar metallicity is difficult to measure using SED fits to mostly broadband photometric data (e.g., Conroy 2013; Lower et al. 2020). Nevertheless, from our analysis (see Paper I) we find that 82% of the galaxies in our sample have log(Z/Z ) < −0.2; this is consistent with the relatively small stellar masses of the galaxies (median log M * /M ∼ 9.5). In other words, for the majority of Hα-emitting galaxies, our 29% correction over-estimates the contamination by [N II], and thus underestimates the population's SFRD. Our MCSED-calibration bypasses this issue by removing the need to correct for [N II]. Topping et al. (2021) do a stacking analysis on z ∼ 1.5 galaxies found in the Multi-Object Spectrometer for Infra-Red Exploration Deep Evolution Field (MOS-DEF) survey. They find that log [N II]/Hα is a strong function of stellar mass, which is correlated with metallicity. From their Figure 2, we notice that for galaxies with log(M * /M ) 10.5, the 29% correction for [N II] is an overestimate. We find that 81% of our galaxies have masses log(M * /M ) < 10.5, which is in perfect agreement with the aforementioned finding using CLOUDY lookup tables.
As for dust, based on results from Paper I, we find that the average differential extinction, E(B − V ), for our galaxy sample is 0.16 mag. To calculate A(Hα) for our sample, we use the conversions given by Reddy et al. (2020), and apply an average value of A(Hα) = 0.88 to the entire sample. The combination of the [N II] and dust corrections results in a multiplicative factor of 1.60.

Star Formation Rate Density Results
Given the relation between SFR and uncorrected Hα + [N II] luminosity (Equation 14), we can calculate the total SFRD using the luminosity function results in this paper. If φ(L) is the true luminosity function, SFR(L) is the relation given by Equation 14, and L is the uncorrected Hα + [N II] luminosity, then the total SFRD contained in 3D-HST ELGs brighter than some luminosity L min is Alternatively, we can use the local calibration between SFR and dust-corrected Hα luminosity (Equation 15) and apply the average dust and [N II] corrections as described in §4.2.1. The net effect of these two factors is included in the constant k = 1.6, and which simplifies to where Γ represents the incomplete gamma function. Figure 7 shows the evolution of the cosmic SFRD based on results compiled by Madau & Dickinson (2014) from both UV-based and IR-based measurements, along with the best-fit curve from their review paper. We also include SFRD measurements from Gruppioni et al. (2020) using sub-mm data. Our SFRD results for both the non-evolving and redshift-varying luminosity functions, based on our SED-based SFR calibration, are shown as a green star and red line, respectively. Our SFRD measurement derived from the non-evolving luminosity function with the Kennicutt & Evans (2012) relation is shown as an amber diamond. In order to make a fair comparison to the literature, we use L min = 0.03L * , as adopted by Madau & Dickinson (2014).
Our measurements derived via the SFR calibration shown in Figure 5 yield an SFRD that is larger by 0.12 dex than that found using the local Hα-SFR calibration given by Kennicutt & Evans (2012). This discrepancy is consistent with the conclusion of Section 4.2.1, i.e., that the application of a 29% correction for the contribution of [N II] to the Hα measurement leads to an underestimate in star-formation rate and the epoch's SFRD. We therefore believe that our SEDbased SFR calibration is more accurate. However, we do note that the uncertainty on the local-calibration-based SFRD is underestimated since the uncertainties on the dust and [N II] corrections are not propagated into the analysis. Thus, the two numbers may still be consistent.
We also observe that the cosmic SFRD calculated for our sample of ELGs is at least ∼ 0.09 dex below the value expected from best-fit curve of Madau & Dickinson (2014). This implies that Hα-selected galaxies contain 81% of the star formation in the z ∼ 1.4 universe. The difference between the cosmic SFRD at z ∼ 1.4, as determined by Madau & Dickinson (2014), and our Hα-based measurement suggests that not all of the epoch's star-formation is detectable via surveys for restframe optical emission lines. Such a result is easily explained if some star-forming galaxies are heavily obscured by dust; since emission-line gas is generally attenuated more than star-light, this is a reasonable hypothesis (e.g., Charlot & Fall 2000;Calzetti et al. 2000;Reddy et al. 2020). These dusty star forming galaxies have been shown to contribute significantly to the star formation rate around cosmic noon, i.e., 1 < z < 3 (e.g., Casey et al. 2013).
At higher redshifts (z 3), far-IR and sub-mm studies have shown that UV-based studies miss a significant portion of the star formation due to dust attenuation in extremely high-SFR galaxies (e.g., Gruppioni et al. 2015;Rowan-Robinson et al. 2016;Wang et al. 2019;Williams et al. 2019;Gruppioni et al. 2020;Loiacono et al. 2021;Khusanova et al. 2021). Nevertheless, the discrepancy between the UV-NIR and FIR-sub-mm methods of calculating star formation (e.g., see the extensive discussion by Katsianis et al. 2021) is much less pronounced at the redshifts of our sample. For example, in Figure 7, we include sub-mm results from Gruppioni et al. (2020) and find that for z = 1 and z = 2, the results are still consistent with the Madau & Dickinson (2014) best-fit curve.
Coming back to our results, we cannot discount the possibility that our SFRD is consistent with the Madau & Dickinson (2014) curve, as the discrepancy is still within the variance of the literature values. Thus, there is no conclusive evidence that the emission line surveys are missing a significant fraction of sources found in UVbased and/or IR/sub-mm-based surveys between redshifts 1.16 < z < 1.56.
While the redshift-evolving SFRD is not a flat curve, given the non-negligible uncertainties due to cosmic variance, we cannot dismiss the possibility of no evolution. In other words, our Hα emission-line measurements find no strong evidence for SFRD evolution over our 1.16 < z < 1.56 redshift range.

[O III] λ5007 Luminosity Function
Unlike Hα, [O III] λ5007 does not suffer from uncertainties due to blending. At the redshifts under consideration, the G141 grism has enough resolution and the sources are sufficiently small so that [O III]  In Figure 8, we show our fitted non-evolving [O III] λ5007 luminosity function (over the range 1.16 < z < 1.90) and place it in the context of the literature. This includes fits to 192 0.7 < z < 1.5 and 58 1.5 < z < 2.3 [O III] emitters by Colbert et al. (2013), 371 z = 1.42 [O III] + Hβ emitters from Khostovan et al. (2015), and 1343 z ∼ 1.4 [O III] + Hβ emitters from Sobral et al. (2015). For our comparison, the [O III] λλ4959, 5007 values of these works have been converted to [O III] λ5007 to match our measurements.
The studies by Khostovan et al. (2015) and Sobral et al. (2015) identify Hβ emitters as well as [O III] galaxies, since their narrow-band photometry is unable to distinguish the two object classes without spectroscopic follow-up. Thus, it is unclear how their samples compare to ours. Sobral et al. (2015) find that in their set  II] are not included in the error bar. This latter calibration produces an SFRD that is smaller by 0.12 dex than SED-based value; this is consistent with the overestimate of [N II] in low-mass (low-metallicity) galaxies. At z ∼ 1.4, the SFRD from ELGs is 0.09 dex lower than the Madau & Dickinson (2014) curve, but this difference is still within the variance seen in the literature. We are not able to find conclusive evidence of SFRD evolution between 1.16 ≤ z ≤ 1.56.
of sources with spectroscopic confirmation, Hβ emitters constitute around 16% of their z ∼ 1.4 sample. Interestingly, they notice that these Hβ emitters tend to have lower luminosities than the [O III]-identified emitters, and this may be reflected in the lower L * values implied by the plot. Still 5/6's of their (spectroscopically confirmed) sample is made up of [O III] emitters, so the effect of these contaminants should not be large. Khostovan et al. (2015) find a similar phenomenon in their sample.
As shown in Figure 8, at 41.5 log L 42.3, our luminosity function predicts more [O III] galaxies than any of the other studies, though at higher luminosities our results are well within the bounds of the literature. Moreover, our MCMC result is in good agreement with the V eff result (blue triangles).  Colbert et al. (2013), Khostovan et al. (2015), and Sobral et al. (2015). We predict higher counts of intermediate-luminosity galaxies than the other studies, but at the high-luminosity end, our data are well within the range of values previously derived. There is good agreement between the MCMC and V eff methods.
There are a few factors that may lead to our distinct result for [O III]. The first is that our redshift range is unique. No other study focuses specifically on 1.16 < z < 1.90 galaxies. As the luminosity function is an evolving quantity, the particular redshifts involved play a role in determining the measured parameters. Another important consideration is that our sample size is the largest to date, and the larger the sample size, the more accurate the luminosity function. In fact, as shown in Table 4, the value for log ∞ 0.03L * φ(L) dL derived by Sobral et al. (2015) is closest to our study, and it is also the measurement that is most consistent with our value.
Alternatively, we note that our analysis encounters problems at both the faint and bright ends of the luminosity function. At the faint end, we are subject to rapid loss of completeness while at the bright end, our measurements of the luminosity function suffer from cosmic variance as we are dealing with small numbers of galaxies. Coupled with degeneracies between Schechter parameters, these effects can lead to divergences between different studies. Extrapolations to lowluminosity galaxies, especially when α is fixed, are often uncertain and should be considered as such. Nevertheless, given our large sample size and the number of galaxy measurements at luminosities significantly less than L * , we believe that the differences between our results and those of previous studies are real and not an artifact of our analysis.
In Figure 9, we combine our [O III] λ5007 measurements with those of Bowman et al. (2021) to show the redshift evolution of the [O III] λ5007 luminosity function between 1.16 ≤ z ≤ 2.35. For this analysis, we fix α = −1.5, as it very difficult to fit α as a free parameter when the limiting luminosity depends so much on redshift. Like in Figure 2, we show the 1D and 2D cross sections of the MCMC chains in the lower left panels. Once again, we find that L * and log φ * are correlated, but only at the same redshift. In other words, L 1 * and log φ 1 * are highly correlated but L 1 * and log φ 2 * are not.
The data of Figure 9 show that the characteristic luminosity L * increases with redshift. Moreover, as redshift increases, the overall luminosity function increases at all but the lowest luminosities. In other words, there are many more [O III] λ5007-visible galaxies, and especially more [O III] λ5007-bright systems, at earlier epochs of cosmic history. This is consistent with the findings of Zeimann et al. (2014), Khostovan et al. (2015), and Bowman et al. (2019)

Number Counts
The precision to which one can measure cosmological parameters through galaxy surveys depends on the number of galaxies and the square of the bias of the observed galaxy population relative to dark matter. Surveys planned for Euclid and Roman will observe millions of Hα-visible galaxies at 0.9 z 1.8 and [O III]-visible galaxies at 1.5 z 2.7, and thus measure quantities such as the angular diameter distance and Hubble parameter at these distant epochs. Specifically, the Euclid Wide Survey (WS) will observe ∼ 15000 deg 2 of the sky down to a flux limit of ∼ 2×10 −16 erg cm −2 s −1 (Euclid Collaboration et al. 2022), while the Roman High Latitude Survey (HLS) will observe ∼ 2200 deg 2 down to a flux limit of ∼ 6 × 10 −17 erg cm −2 s −1 (Spergel et al. 2015). The Euclid Deep Survey will observe 50 deg 2 in three fields to a similar flux limit (Vavrek et al. 2016).
In this section, we use our measurements of the Hα and [O III] λ5007 luminosity functions to calculate the number of galaxies these surveys are likely to measure, corrected for completeness. One factor to note is that since we have removed AGN from the sample, the number counts we derive will be a slight underestimate.
We begin with the total number of galaxies with luminosities greater than L min (z) from redshifts z 0 to z 1 over area Ω, dL Ω(L, z)φ(L, z) Given the Schechter function with parameters α, L * , and φ * (that are assumed to be invariant over the effective survey area Ω 0 ), we can simplify Equation 21 to where Γ once again represents the incomplete gamma function.
In Figure 10, we show the total number of galaxies per square degree that are expected to have Hα at 1.2 ≤ z ≤ 1.6 (left) and [O III] λ5007 at 1.5 ≤ z ≤ 1.9 (right) above a given threshold. We calculate these values using Equation 22 with Ω 0 = 1 deg 2 , while translating L min to emission line flux using the luminosity distance D L via For the figure, we use the non-evolving luminosity functions (though the redshift varying luminosity functions yield similar results) and perform the same calculations for the luminosity functions given in the literature. In all cases, the values represent galaxy number counts, assuming 100% completeness. In addition, we include the direct measurements of counts obtained by Bagley et al. (2020) using samples of ELGs from the WFC3 Infrared Parallel Spectroscopic Survey (WISPS), 3D-HST, and A Grism H-Alpha SpecTroscopic survey (AGHAST). This work also made corrections for completeness.
In the left panel, we observe that our Hα counts are in excellent agreement with those of Colbert et al. (2013) at all values of the limiting flux. The counts also agree with those of Bagley et al. (2020) around the limiting flux of the Euclid WS, but less so at brighter fluxes. Given the consistency of our results with those of Colbert et al. (2013) and Bagley et al. (2020), there is no need to update the prediction of ∼ 3300 deg −2 for the Euclid WS at 100% completeness (Bagley et al. 2020;Euclid Collaboration et al. 2022) and 16.4 million 7σ Hα galaxy detections at 1.06 < z < 1.88 for the Roman HLS assuming 70% completeness (Spergel et al. 2015).
On the other hand, our result for [O III] (right side of Figure 10) is distinct from the literature. For limiting fluxes of 1 to 3 × 10 −16 erg cm −2 s −1 , our predicted galaxy counts agree with the results of Colbert et al. (2013) and Bagley et al. (2020). As these limits are similar to those of the Euclid WS (2 × 10 −16 erg cm −2 s −1 ) and the ≥ 7σ limit of 10 −16 erg cm −2 s −1 used by Spergel et al. (2015) for predicting Roman HLS counts, our results do not change the predictions for these surveys; at 70% completeness, the Roman HLS should detect 1.4 million [O III] 1.88 < z < 2.77 galax- In the lower left, we show the 1D and 2D cross sections of the parameter chains. Note that at each redshift, L * and log φ * are correlated, but the correlation is not strong across redshifts. The upper-right panel shows the evolution of the luminosity function. Our data show that L * increases with redshift, and, at most luminosities, the luminosity function is larger at higher redshift. In other words, [O III] λ5007-bright galaxies were much more common at z ∼ 2 than at z ∼ 1.
ies with > 7σ confidence (Spergel et al. 2015). However, as we push to lower flux limits, our luminosity function predicts higher counts. For example, for the Roman HLS nominal flux limit of 6.0 × 10 −17 erg cm −2 s −1 , we predict twice as many galaxies as Colbert et al. (2013) in the redshift range 1.5 < z < 1.9.

CONCLUSION
In Papers I and II, we used the G141 grism data from the Hubble 3D-HST Treasury program (GO-11600, 12177, 12328;Brammer et al. 2012;Momcheva et al. 2016) to identify a clean sample of 4350 normal (i.e., non-AGN) star-forming emission lines galaxies in the redshift range 1.16 ≤ z ≤ 1.90. In this work, we use the line fluxes of these galaxies, along with a similar set of fluxes measured by Bowman et al. (2021), to measure the galaxies' emission-line luminosity functions. These data include 1892 Hα emitting galaxies between 1.16 ≤ z ≤ 1.56 and 4519 [O III] emitters with 1.16 ≤ z ≤ 2.35, all with line fluxes above the 3D-HST 50% completeness limit. While there have been several previous efforts to calculate Hα and [O III] luminosity Figure 10. Predictions for total number of galaxies per square degree corrected for completeness, as a function of the limiting flux for Hα at 1.2 ≤ z ≤ 1.6 (left) and [O III] λ5007 at 1.5 ≤ z ≤ 1.9 (right). We use our non-evolving luminosity functions for these calculations and include cosmic variance in the error budget. Our Hα counts agree remarkably well with those from Colbert et al. (2013) at all limiting fluxes as well as with the direct counts from Bagley et al. (2020) around F lim ∼ 2 × 10 −16 erg cm −2 s −1 . This agreement suggests no changes to the Hα predictions for Euclid and Roman are necessary. While our [O III] λ5007 count predictions are distinct from the literature, they still do not greatly change the expected counts from Euclid. However, for the HLS flux limit of 6.0 × 10 −17 erg cm −2 s −1 , we predict twice as many galaxies as Colbert et al. (2013). functions in these redshift ranges (see §1 for a detailed list), none used the large sample sizes used here.
We employ a generalization of the classical 1/V max method to derive the emission-line luminosity functions for our entire sample of galaxies, and samples of galaxies broken down by redshift. We then use Markov Chain Monte Carlo (MCMC) Bayesian techniques to fit these data to the Schechter (1976) luminosity function with α held constant across redshift. We find very good agreement between our Hα results and those from the literature, and our [O III] luminosity function is also a good match to prior measurements for line luminosities brighter than log L = 42.3 (ergs s −1 ). However, at fainter [O III] luminosities (where completeness corrections might be an issue), we infer an excess of objects. These results are shown in Figures 2, 3, 8, and 9 and summarized in Tables 3 -5. We also compute the star formation rate density (SFRD) of the 1.16 ≤ z ≤ 1.56 epoch using the Hα luminosity function. We find that our SFRD is ∼ 19% smaller than the best-fit z ∼ 1.4 value found by Madau & Dickinson (2014) (Figure 7), though this discrepancy is within the variance found in the literature. If the difference is real, then one possible explanation is that not all z ∼ 1.4 star formation takes place in galaxies with observable Hα emission lines. In particular, surveys such as 3D-HST will miss heavily obscured galaxies where the emission lines are too extinguished to make it into the sample. We find no evidence for or against cosmic evolution of the SFRD between 1.16 < z < 1.56 Finally, we predict total galaxy counts per square degree as a function of the limiting flux ( Figure 10). For Hα + [N II] λ6584, our results are consistent with those from Colbert et al. (2013) at all limiting fluxes and Bagley et al. (2020) down to a limiting flux of ∼ 2×10 −16 erg cm −2 s −1 , suggesting the previous predictions for the Euclid (Euclid Collaboration et al. 2022) and Roman (Spergel et al. 2015) surveys are accurate. For [O III] λ5007, our numbers agree with previous estimates for the Euclid Wide Survey, but depending on where exactly we define the flux limit for the Roman High Latitude Survey, our data may imply a significantly larger number of detectable galaxies. For example, at 10 −16 erg cm −2 s −1 , our results are consistent with the analysis of Colbert et al. (2013), but at 6 × 10 −17 erg cm −2 s −1 , our number counts are higher by a factor of two. Roman may find more faint [O III] emitters than previously anticipated.
The ΛCDM paradigm has garnered many resounding successes in explaining observations of our universe at a variety of scales. However, there are still inconsistencies and unknowns leaving the cosmological model incomplete. To better constrain the model as well as alternate or additional theories, we need to continue honing our observations. Large galaxy surveys represent an important avenue to constrain cosmological parameters through the measurement of baryonic acoustic oscillations and redshift space distortions.
Galaxy surveys with precise redshifts will be especially useful for generating the necessary constraints, and IFU and slitless spectroscopy are the most efficient ways of performing these surveys. In the near future, Euclid and Roman will greatly enhance samples of emission line galaxies. The similarities between 3D-HST and these planned surveys make it a perfect pathfinder mission. Our measurements of the Hα and [O III] λ5007 luminosity functions with galaxy samples that are several times larger than any previous study of the z ∼ 1.5 redshift range help cement the predictions for the expected yield of the ongoing and future surveys.