Mengtao Tang

Daniel P Stark

Jacopo Chevallard

Stéphane Charlot

MMT/MMIRS spectroscopy of z = 1.3 -2.4 extreme [OIII] emitters: Implications for galaxies in the reionization-era

Keywords: cosmology, observations -galaxies, evolution -galaxies, formationgalaxies, high-redshift

. We describe how this relationship helps explain the anomalous success rate in identifying Lyα emission in z 7 galaxies with strong [O III]+H β emission. We probe the impact of the intense radiation field on the ISM using O32 and Ne3O2, two ionization-sensitive indices. Both are found to scale with the [O III] EW, revealing extreme ionization conditions not commonly seen in older and more massive galaxies. In the most intense line emitters, the indices have very large average values (O32 = 9.1, Ne3O2 = 0.5) that have been shown to be linked to ionizing photon escape. We discuss implications for the nature of galaxies most likely to have O32 values associated with significant LyC escape. Finally we consider the optimal strategy for JWST spectroscopic investigations of galaxies at z 10 where the strongest rest-frame optical lines are no longer visible with NIRSpec.

INTRODUCTION

Over the past decade, deep imaging surveys with the Hubble Space Telescope (HST) have uncovered thousands of colorselected galaxies at z > 6 (e.g. McLure et al. 2013;Bouwens et al. 2015a;Finkelstein et al. 2015;Livermore et al. 2017;Atek et al. 2018;Oesch et al. 2018), revealing an abundant population of low mass systems that likely contribute greatly to the reionization of intergalactic hydrogen (Bouwens et al. 2015b;Robertson et al. 2015;Stanway et al. 2016;Dayal & Ferrara 2018). These galaxies have been shown to be compact with low stellar masses, blue ultraviolet continuum tangmtasua@email.arizona.edu slopes, and large specific star formation rates (see Stark 2016 for a review).

The James Webb Space Telescope (JWST) will eventually build on this physical picture, providing spectroscopic constraints on the stellar populations and gas conditions within the reionization era. Dedicated efforts throughout the last decade have delivered our first glimpse of what we are likely to see with JWST. Ground-based spectroscopic surveys have revealed a rapidly declining Lyα emitter fraction at z > 6 (e.g. Caruana et al. 2014;Pentericci et al. 2014;Schenker et al. 2014), suggesting that the intergalactic medium (IGM) is likely partially neutral at z 7 -8. While prominent Lyα emission is very rare in early galaxies, the rest-frame optical emission lines are found to be very strong relative to the underlying continuum, as expected for galaxies with large specific star formation rates (sSFR). The average [O III]+Hβ equivalent width (EW) inferred from broadband flux excesses in stacked spectral energy distributions (SEDs) at z 7 -8 is found to be 670 Å (Labbé et al. 2013, hereafter L13), while individual objects have been located with combined [O III]+Hβ equivalent widths reaching up to 1000 -2000 Å (Smit et al. 2014(Smit et al. , 2015;;Roberts-Borsani et al. 2016). While direct spectral measurements with JWST will ultimately be required to confirm the z ∼ 7 equivalent width distribution, these measurements suggest large [O III] EWs (> 300 Å) are likely to be very common.

These results indicate that the majority of reionizationera galaxies are caught in an extreme emission line phase, implying a large ratio of the flux density in the ionizing UV continuum (λ < 912 Å) to that in the optical (λ 4000-6000 Å). There are a variety of ways a stellar population can produce such spectra. A likely pathway follows a recent upturn in the star formation rate (SFR), as would occur for systems undergoing rapidly rising star formation histories (which are expected at z 6; i.e. Behroozi et al. 2013) or punctuated bursts of star formation that greatly exceed earlier activity. In both cases, the stellar population will be dominated by recently formed stars, leading to an enhanced contribution from O stars relative to the longer-lived main-sequence stars that dominate the non-ionizing UV continuum and optical continuum. In this framework, the galaxies with the largest optical line equivalent widths (i.e., those approaching EW = 1000 -2000 Å) will be dominated by extremely hot massive stellar populations, resulting in an intense extreme UV (EUV) radiation field.

The interstellar gas reprocesses this radiation field, powering the observed nebular emission lines. Uncertainties in the form of the ionizing spectrum thus have a significant effect on the interpretation of the nebular spectrum. While the implementation of new physics in population synthesis models (i.e., stellar multiplicity, rotation) continues to improve predictions of the EUV radiation field (e.g. Levesque et al. 2012;Eldridge et al. 2017;Götberg et al. 2018), the shape of the ionizing spectrum between 1 and 4 Ryd remains very poorly constrained for a given age and metallicity. This is perhaps especially true among the highest equivalent width line emitters for which the EUV spectrum is extremely sensitive to the nature of the hottest and most massive stars present in the galaxy. While large spectroscopic surveys have begun to constrain the ionizing spectrum of typical galaxies at z 2 -3 (e.g. Sanders et al. 2016;Steidel et al. 2016;Strom et al. 2017), much less work has focused on the extreme emission line population.

As the first nebular spectra of z 6 galaxies emerge, we are already beginning to see indications that the ionizing spectra of reionization-era sources differ greatly from typical galaxies at lower redshifts. Perhaps most surprising has been the detection of strong nebular C IV emission in two of the first galaxies observed at z 6 (Stark et al. 2015b;Mainali et al. 2017;Schmidt et al. 2017), implying a very hard EUV ionizing spectrum capable of triply ionizing carbon in the ISM. The presence of C IV in a high redshift galaxy typically indicates the power law spectrum of an active galactic nucleus (AGN). However, Mainali et al. (2017) have demonstrated that the observed line ratios in one of the z > 6 C IV emitters requires a strong break at the He + -ionizing edge, suggesting a metal poor (< 0.1 Z ) stellar ionizing spec-trum is likely responsible for powering the nebular emission. The precise form of the EUV spectrum produced by massive stars in this metallicity regime is very uncertain, making it difficult to reliably link the observed nebular emission to a unique physical picture. Similar challenges have emerged following detections of intense C III] emission in two galaxies at z 6 (Stark et al. 2015a(Stark et al. , 2017)). Whereas some have argued that the UV line emission can be powered by moderately metal poor, young stellar populations (Stark et al. 2017), others suggest an AGN contribution is necessary to explain the large C III] equivalent width (Nakajima et al. 2018). Until more is known about the EUV radiation field of very young stellar populations that dominate in extreme emission line galaxies (EELGs), it will be difficult to reliably interpret the the z 6 nebular line spectra that will be common in the JWST era.

The uncertainties in the EUV spectra of early galaxies also plague efforts to use Lyα as a probe of reionization. The success of the Lyα reionization test hinges on our ability to isolate the impact of the IGM on the evolving Lyα equivalent width distribution, requiring careful attention to variations in the intrinsic production of Lyα within the galaxy population. The intrinsic Lyα luminosity of a given stellar population (per SFR) depends sensitively on the production efficiency of Lyman continuum (LyC) radiation (e.g. Wilkins et al. 2016), often parameterized as ξ ion = N(H 0 )/L UV , the ratio of production rate of hydrogen ionizing photons (N(H 0 )) and the continuum luminosity of non-ionizing UV photons (L UV ). By definition, the extreme emission line galaxies that dominate in the reionization era should be extraordinarily efficient at producing hydrogen ionizing radiation. Efforts to measure ξ ion from both high quality spectra (e.g. Chevallard et al. 2018;Shivaei et al. 2018) and narrowband imaging (Matthee et al. 2017) have ramped up in recent years. Chevallard et al. (2018) demonstrated a correlation exists between ξ ion and [O III] EW among ten nearby galaxies, such that the most extreme line emitters are the most efficient in their ionizing photon production. The first hints at high redshift suggest a similar picture (Matthee et al. 2017;Shivaei et al. 2018), but the LyC production efficiency of extreme emission line galaxies remains poorly constrained. As a result, it is not known the extent to which ξ ion varies over the full range of equivalent widths present in reionization era galaxies (EW [OIII]+Hβ 300 -2000 Å), making it challenging to accurately account for how variations in the radiation field are likely to impact Lyα visibility.

The limitations of our physical picture of Lyα production in the reionization era has recently become apparent following the discovery of four massive z 7 -9 galaxies with extremely large [O III]+Hβ equivalent widths (EW [OIII]+Hβ 1000-2000 Å;Roberts-Borsani et al. 2016, hereafter RB16). Spectroscopic follow-up has revealed Lyα in each of the four galaxies from RB16 (e.g. Oesch et al. 2015;Zitrin et al. 2015;Stark et al. 2017), corresponding to a factor of ten higher success rate than was found in earlier studies (e.g. Schenker et al. 2014;Pentericci et al. 2014). The association between Lyα and extreme EW [O III] emission is also seen in other z > 7 Lyα emitters with robust IRAC photometry (Ono et al. 2012;Finkelstein et al. 2013). Remarkably these results suggest significant differences between the Lyα properties of the extreme line emitters that are typical at z 7 (EW [OIII]+Hβ 670 Å) and massive systems with slightly larger equivalent widths (EW [OIII]+Hβ 1000 -2000 Å). This could imply that the transmission of Lyα is enhanced in the higher equivalent width sources, as might be expected if these massive galaxies trace larger-than-average ionized patches of the IGM. Or it could be explained if the production efficiency of Lyα is much larger in the higher EW sources. In this case, efforts to model the evolving Lyα EW distribution would have to control for the large intrinsic variations in the Lyα luminosity per SFR that accompany sources of slightly different [O III] EW. Unfortunately without robust measurements of how ξ ion scales with [O III] EW at high redshift, it is impossible to know which picture is correct.

The shortcomings described above motivate the need for a comprehensive investigation of the spectral properties of high redshift extreme emission line galaxies. Prior to the launch of JWST, this is most easily accomplished at z 2 where both rest-UV and rest-frame optical nebular lines are visible from the ground. In this paper, we present the first results from a near-infrared spectroscopic survey targeting rest-frame optical nebular emission lines in over 200 extreme emission line galaxies at z = 1.3 -2.4. Our sample is selected using a combination of HST grism spectra (Brammer et al. 2012;Skelton et al. 2014;Momcheva et al. 2016) and broadband photometry. Over the past 2.5 years, we have obtained 100 hours of MMT and Keck spectroscopy following up this population in the near-infrared. The program is ongoing, and we ultimately aim to obtain useful constraints on the strongest rest-frame optical lines (

[O II], [Ne III], Hβ, [O III],
Hα, [N II]) for each galaxy in our sample. This program complements the KBSS (Steidel et al. 2014) and MOSDEF (Kriek et al. 2015) programs, including galaxies with lower masses and larger specific star formation rates than are common in these surveys.

We focus on three key topics in this paper. First we use our spectra to characterize how the production efficiency of LyC (and hence Lyα) photons scales with the [O III] and Hα EW. Using this information, we consider how variations in Lyα production efficiency impact the visibility of Lyα at z > 7. Second, we investigate how the intense radiation field of extreme emission line galaxies impacts the conditions of the ISM, characterizing how rest-frame optical line ratios that are sensitive to the ionization state of the gas vary with [O III] EW. We discuss the findings in the context of results showing an association between ionization-sensitive line ratios and the escape fraction of ionization radiation. Third we use our emission line measurements to help optimize spectroscopic observations of the highest redshift galaxies (z 10) that JWST will target. At these redshifts, the strongest rest-frame optical lines ([O III], Hα) are no longer visible with NIRSpec. We consider whether fainter rest-frame optical lines (i.e., [O II], [Ne III]) will provide viable alternatives and make predictions for the range of equivalent widths that are likely to be present in early galaxies.

The organization of this paper is as follows. We describe the sample selection and near-infrared spectroscopic observations in §2. We introduce our photoionization modeling procedure in §3 before discussing the rest-frame optical spectroscopic properties of the EELGs in §4. We discuss the implications in §5, and summarize our results in §6. We adopt a Λ-dominated, flat Universe with Ω Λ = 0.7, Ω M = 0.3 and H 0 = 70 km s -1 Mpc -1 . All magnitudes in this paper are quoted in the AB system (Oke & Gunn 1983), and all equivalent widths are quoted in the rest-frame. We assume a Chabrier (2003) initial mass function throughout the paper.

OBSERVATIONS AND ANALYSIS

We have initiated a large spectroscopic survey of z 2 galaxies with extremely large equivalent width optical emission lines. A complete description of the survey and the galaxies targeted will be presented in a catalog paper following the end of the survey. Here we provide a detailed summary of the survey and resulting emission line sample, describing the pre-selection of targets ( §2.1) and our campaign to follow-up these sources with spectrographs on MMT and Keck ( §2.2-2.4). We discuss the analysis of the spectra and our current redshift catalog in §2.5.

Pre-Selection of Extreme Emission Line Galaxies

The first step in our survey is to identify a robust sample of extreme emission line galaxies for spectroscopic followup. Since we aim to study how the radiation field and gas conditions vary within the reionization era population, we must select galaxies at that span the full range of equivalent widths expected at z 6 (EW [OIII]+Hβ 300 -2000 Å; e.g., Stark 2016). This pre-selection is most efficiently performed using publicly-available HST WFC3 G141 grism spectra in the five CANDELS (Grogin et al. 2011;Koekemoer et al. 2011) fields. The G141 grism covers between 1.1 µm and 1.7 µm, enabling identification of Hα at 0.7 < z < 1.6 and [O III]λλ4959, 5007 (blended at the grism resolution) at 1.2 < z < 2.4. Four of the fields (AEGIS, COSMOS, GOODS-S, UDS) have been targeted by the 3D-HST program (van Dokkum et al. 2011;Brammer et al. 2012;Skelton et al. 2014;Momcheva et al. 2016), and the fifth field (GOODS-N) was observed in GO-11600 (PI:Weiner). The latest 3D-HST public release brings together all the G141 data, providing photometry (Skelton et al. 2014) and spectral measurements (Momcheva et al. 2016) of galaxies in all five CANDELS fields.

Our primary selection is focused on systems in two redshift windows (1.37 ≤ z ≤ 1.70, 2.09 ≤ z ≤ 2.48) where the full set of strong rest-frame optical lines ([O II], Hβ, [O III], Hα) can be observed with ground-based spectrographs. At 1.37 ≤ z ≤ 1.70, we can detect Hβ and [O III] in the J band, and Hα in the H band. [O II] is situated in the Y band at 1.55 ≤ z ≤ 1.70 and in the optical at z < 1.55. At 2.09 ≤ z ≤ 2.48, we can detect [O II] in the J band, Hβ and [O III] in the H band, and Hα in the K band. We also identify secondary targets at 1.70 < z < 1.90 and 2.02 < z < 2.09, where a subset of the lines are visible using ground-based spectrographs. We will discuss the relative priorities we assign to galaxies when introducing the near-infrared spectroscopic follow-up in §2.2.

As a first pass at our pre-selection, we use the 3D-HST v.4.1.5 grism catalogs (see Momcheva et al. 2016) to select extreme line emitters in all five CANDELS fields. These catalogs include extracted grism spectra for objects with WFC3 F140W magnitude < 26.0. We aim to select galaxies with optical line equivalent widths similar to those seen in the reionization era. Stacked broadband SEDs imply typical values of [O III]+Hβ EW = 670 Å at z 8 (L13). To ensure our sample spans the full range of line strengths seen at z 8, we consider all objects with [O III]+Hβ EW > 340 Å, roughly a factor of two below the average value quoted above. We thus select all galaxies in the grism catalogs with [O III]λλ4959, 5007 EW > 300 Å.1 Here we have assumed a flux ratio of [O III]λ5007/Hβ = 6, consistent with measurements of extreme emission line galaxies at z 2 (Maseda et al. 2014). For the sake of simplicity, we use the same equivalent width threshold for Hα. The precise value we choose does not impact our results, as we will sample a wide range of equivalent widths in our follow-up survey. For reference, if we crossmatch the MOSDEF spectroscopic sample to the 3D-HST grism catalogs (considering only galaxies in same redshift range as in our EELG sample), we find that the median [O III]λλ4959, 5007 EW is 154 Å for sources with S/N > 3 grism line detections. Our sample thus largely probes a distinct [O III] EW regime.

The efficiency of our follow-up spectroscopy program depends sensitively on how accurately the equivalent widths can be measured from the grism spectra, requiring robust measurements of both the line flux and the underlying continuum. To ensure the line detections are robust, we only include sources with significant (S/N > 3) [O III] or Hα flux and equivalent width measurements. However we find that the equivalent width S/N threshold excises many extreme line emitters from our sample owing to the low S/N continuum in many of the grism spectra. For these fainter continuum sources, we can more reliably identify extreme line emitters using the line flux from the grism catalogs and the continuum from the HST broadband imaging data. We estimate the equivalent width errors for this subset by adding the uncertainties of emission line fluxes and broadband photometry in quadrature. If the rest-frame [O III]λλ4959, 5007 or Hα equivalent width calculated in this manner is above 300 Å with S/N > 3, we include the source in our catalog for follow-up spectroscopy.

Finally we visually inspect the grism spectra and photometry of all galaxies that satisfy these criteria, removing sources with grism spectra that are unreliable owing to contamination from overlapping spectra or those that appear to have misidentified emission lines. We are left with a sample of 1587 galaxies with grism measurements that imply [O III]λλ4959,5007 or Hα rest-frame EW greater than 300 Å across the five CANDELS fields. Examples of the broadband SEDs of the extreme emission line galaxies in our sample are shown in Figure 1. In all cases, the measured flux in the specific broadband filters that sample the strong nebular emission lines are well above the flux of the underlying continuum.

The grism observations cover ∼ 75% of the area imaged by the CANDELS survey. While our spectroscopic follow-up observations are centered in these regions, often the edges of the slitmasks extend to areas of the CANDELS fields that lack grism spectra. In these sub-regions, we use simple color cuts to identify extreme emission line galax-ies based on the characteristic shape of their broadband SEDs. In the redshift range 1.6 < z < 1.8, galaxies with large equivalent width [O III]+Hβ emission show a large flux excess in the J 125 -band (Figure 1). Previous studies have successfully demonstrated that extreme [O III]+Hβ emitters can be easily selected by identifying galaxies in which the J 125 -band flux is significantly greater than that in the I and H 160 -bands (van der Wel et al. 2011;Maseda et al. 2014). We closely follow these selections, identifying galaxies in the Skelton et al. (2014) photometric catalogs that satisfy the following criteria: I -J 125 > 0.3 + σ(I -J 125 ) and J 125 -H 160 < -0.3-σ(J 125 -H 160 ), where σ is the color uncertainty. The required J-band flux excess is chosen to match the equivalent width threshold of our grism sample. Whereas the grism selection identifies strong [O III]λλ4959,5007 or Hα emission, the photometric selection picks out sources based on their combined [O III]+Hβ emission. Under the same assumptions described above, our color cuts correspond to a selection of galaxies at 1.6 < z < 1.8 with [O III]+Hβ EW > 340 Å. We have verified (in fields where grism spectra are available) that this selection identifies the same extreme emission line galaxies selected via the grism catalogs.

Following Maseda et al. (2014), we also select [O III]+Hβ emitter candidates at 2.0 < z < 2.4 via an H 160 -band flux excess technique. Because of the limited K-band sensitivity, the color cuts identify sources in which H 160 is significantly in excess of J 125 . To minimize contamination from Balmer break sources, Maseda et al. (2014) introduce a requirement that galaxies additionally have blue UV continuum (in V -I) slopes. We adopt similar color cuts, again making minor adjustments to the required color excess for the sake of consistency with our grism selection. In particular, we identify sources with V 606 -I < 0.1 -σ(V 606 -I) and J 125 -H 160 > 0.36 + σ(J 125 -H 160 ). We also require that the detections in V 606 , I, J 125 and H 160 have S/N > 5. The two color section techniques provide an additional 612 targets, resulting in a total sample of 2199 galaxies in our input catalog.

The final step is to remove galaxies that likely host AGNs from the sample. While recent studies of EELG at z 7 suggests that AGN may be present in several of the sources for which we are obtaining spectra (e.g. Tilvi et al. 2016;Laporte et al. 2017;Mainali et al. 2018), our goal is to focus first on the range of spectral properties in galaxies dominated by star formation. We use the deep Chandra X-ray imaging in AEGIS (Nandra et al. 2015), COSMOS (Civano et al. 2016), GOODS-N (Alexander et al. 2003;Xue et al. 2016), and GOODS-S (Xue et al. 2011) fields, as well as XMM X-ray imaging in UDS field (Ueda et al. 2008) to identify X-ray sources in the sample. We match the coordinates of our EELG candidates to the X-ray source catalogs using a 1. 0 search radius. There are 26 sources in our sample found to have X-ray counterparts within 1. 0. We remove these 26 sources, leaving 2173 targets in the final grism spectroscopy and photometry-selected EELG candidate sample.

MMT/MMIRS Spectroscopy

We use the catalog of extreme equivalent width line emitters described in §2.1 as input for our ground-based spectroscopic follow-up program, the majority of which is conducted at the MMT using the MMT and Magellan Infrared Spectrograph (MMIRS; [START_REF] Mcleod | [END_REF]. MMIRS is a wide-field nearinfrared imager and multi-object spectrograph (MOS) with a field of view of 4 × 6 .9. We used the "xfitmask 2 " software to design our MMIRS slit masks, using a slit length of 7. 0 and a slit width of 1. 0 for science targets. Data were obtained over six observing runs between the 2015B and 2018A semesters. Thus far, we have collected 81 hours of on-source integration with MMIRS, allowing us to target 313 galaxies on 17 separate masks, with 58 galaxies observed on two masks.

Our selection function is defined primarily by the redshift and equivalent width of the targets. We first prioritize galaxies at 1.37 ≤ z ≤ 1.70 and 2.09 ≤ z ≤ 2.48, where the full set of strong rest-frame optical emission lines can be obtained from ground-based observations. For galaxies at 1.37 ≤ z ≤ 1.70, we give the highest priority to those at 1.55 ≤ z ≤ 1.70, where all strong rest-frame optical lines can be probed by MMIRS observations. Galaxies at z < 1.55 will need additional follow-up with red-sensitive optical spectrographs to detect the [O II] emission line. We next adjust the target priority in each redshift interval based on the grism equivalent width or the flux excess implied by photometry. Sources with the highest equivalent widths are extraordinarily rare, with only 15 -20 (3) galaxies per CANDELS 2 http://hopper.si.edu/wiki/mmti/MMTI/MMIRS/ObsManual/ MMIRS+Mask+Making field with [O III]λ5007 EW > 1000 (2000) Å. To ensure that our sample contains an adequate number (i.e., 5 -10) of the largest equivalent width galaxies, we increase the priority of objects with the largest equivalent widths in each redshift interval. This is the minimum number required for us to characterize the average line ratios in this [O III] EW range. Finally, we prioritize grism-selected targets over those identified by photometric excesses when both are available in the same region of a slitmask.

The details of the MMIRS masks obtained to date are summarized in Table 1. Spectra were taken with the J grism + zJ filter, H3000 grism + H filter, and K3000 grism + Kspec filter, providing wavelength coverage of 0.95-1.50 µm, 1.50 -1.79 µm and 1.95 -2.45 µm, respectively. The 1. 0 slit widths result in resolving power of R = 960, 1200, and 1200 for J + zJ, H3000 + H, and K3000 + Kspec grism and filter sets, respectively. The total integration time in each filter ranges from 1 hr to 8 hr, with the specific value chosen depending on the predicted brightness of the set of emission lines we are targeting. The average seeing was between 0. 6 and 1. 6 (see Table 1). The masks all contained an isolated star to monitor the throughput and seeing during the observations and to compute the absolute flux calibration. For each mask and filter combination, we also observed A0V stars at a similar airmass to derive the response spectrum and correct for telluric absorption.

We reduced our MMIRS data using the publicly avail-able data reduction pipeline3 developed by the instrument team. The MMIRS data reduction pipeline is implemented in IDL and described in Chilingarian et al. (2015). Individual two-dimensional (2D) spectra are extracted from the original frames and flat-fielded after spectral tracing and optical distortion mapping. The wavelength calibration uses atmosphere airglow OH lines, or the internal argon arc lamp if the OH based computation fails. For sky subtraction, the pipeline creates a night sky spectrum model using science frames and the precise pixel-to-wavelength mapping, and subtracts the night sky emission from the 2D spectra. Finally, the pipeline provides telluric correction by first computing the atmosphere transmission curve (including spectral response function of the detector) from the observed telluric star spectra, and then correcting the transmission using a grid of atmospheric transmission models to account for the airmass difference between observations of the science target and the telluric star. The pipeline outputs the sky-subtracted and telluric-corrected 2D spectra of every slit on the mask.

One-dimensional (1D) spectra are then extracted from the fully-reduced 2D spectra. In cases where an emission line or continuum is detected with confidence, we use an optimal extraction procedure (Horne 1986). In all other cases we use a simple boxcar extraction, with the extraction aperture matched to the object spatial profile. We have verified that the extraction method does not significantly impact on the line flux measurements. The absolute flux calibration is then applied to the extracted 1D spectra using the slit stars. Slit loss corrections are performed following the similar procedures in Kriek et al. (2015). We first extract a postage stamp of each galaxy from the HST F160W image (Skelton et al. 2014), then smooth the postage stamp and fit the smoothed image with a 2D Gaussian profile. We compute the fraction of the light within the slit to that of the total Gaussian profile. Each spectrum is multiplied by the ratio of the in-slit light fraction measured for the slit star to that for each galaxy.

Keck/MOSFIRE Spectroscopy

We also obtained spectra from the Multi-object Spectrometer for Infrared Exploration (MOSFIRE; [START_REF] Mclean | Ground-based and Airborne Instrumentation for Astronomy IV[END_REF] on the Keck I telescope. We placed galaxies on masks using the same selection function described in §2.2. We used the MAGMA4 software to design the slit masks, and adopted a slit width of 0. 7. Data were obtained on observing runs in 2015 November (one mask in the COSMOS field, one mask in the GOODS-S field, and one mask in the UDS field), and 2016 April (two masks in the AEGIS field). We have summarized the details of the Keck/MOSFIRE observations in Table 2. Spectra were taken in the Y, J, and H bands, providing wavelength coverage of 0.972 -1.125 µm, 1.153 -1.352 µm, and 1.468 -1.804 µm. The slit width of 0. 7 results in resolving power of R = 3388, 3318, and 3660 for Y, J, and H, respectively. The total integration time for each mask in each band ranges from 0.5 hr to 4.3 hr. The average seeing during the observations was between 0. 46 and 0. 94. The masks all contained 1 -2 isolated stars to monitor the throughput and seeing. We use these stars to derive the absolute flux calibration. We also observed A0V stars in Y, J, and H bands to derive the response spectrum and correct for telluric absorption. In total, we observed 94 targets on 5 separate masks with MOSFIRE.

Our MOSFIRE data were reduced using the publicly available data reduction pipeline5 (DRP), which is implemented in PYTHON. The MOSFIRE DRP first generates a flat-fielded image and traces the slit edge for each mask. Using a median-combined image of all science exposures, wavelength solutions are fit interactively for the central pixel in each slit using the night sky OH emission lines, and propagated spatially along the slit. After wavelength calibration and background subtraction, the two A-B and B-A stacks are shifted, combined and rectified to produce the final 2D spectra. We derive the telluric correction spectra using longslit observations of a telluric standard star. The absolute flux scaling factors are computed by comparing the count rates of slit star spectra with the flux densities in the 3D-HST photometric catalogs. Slit loss corrections were applied following the procedures in Kriek et al. (2015) as described in §2.2.

MMT/Red Channel Spectroscopy

For the subset of our extreme emission line sample at z < 1.55, we require optical spectra to detect the [O II] doublet. As a first step toward obtaining such spectra, we observed two galaxies at z < 1.55 from our spectroscopic sample using the Red Channel Spectrograph (Schmidt et al. 1989) on the MMT telescope on 2018 January 09 (UT). The objects have extremely large [O III]λλ4959, 5007 equivalent widths (EW 2500 Å for UDS-29267; EW 650 Å for UDS-14019) and bright continuum fluxes (V 606 = 23.8 for both galaxies), making them ideal targets for securing [O II] detections in feasible integrations with MMT Red Channel.

The observations are summarized in Table 3. We used the 1200 lines mm -1 grating centered at 9290 Å with a UV-36 order blocking filter, providing spectral coverage of 805 Å. The observations were conducted in long slit mode, and we used a 1. 0 slit width, providing a spectral resolution of 1.8 Å. We observed object UDS-29267 for 2 hours and UDS-14019 for 1.67 hours, with individual exposures of 20 minutes. Thin clouds were present during observations of both sources. The average seeing was 1. 1 during the Red Channel observations of UDS-29267 and 1. 3 for UDS-14019.

The Red Channel long slit spectra were reduced using standard IRAF routines. Wavelength calibration was performed with HeAr/Ne arcs. We corrected the atmosphere transmission and instrument response by using the standard star spectrum. A slit star was also put on each slit, allowing us to perform an absolute flux calibration. Slit loss corrections were performed following the methods described in §2.2.

Emission Line Measurements

Examples of one-dimensional spectra are shown in Figure 2. Redshift confirmation is performed by visual inspection of the MMT and Keck spectra. We require two emission line detections for a robust redshift measurement. Currently we have confirmed redshifts of 240 of the galaxies from our grism input catalog of extreme line emitters (see §2.1). The success rate for redshift confirmation is 92% for those objects with complete near-infrared spectral coverage. In the remaining 8% of objects for which we fail to measure a red- shift, the spectra are either very low S/N or the emission lines are contaminated by sky line residuals. The photometric flux excess selection contributes an additional 18 galaxies, resulting in a total sample of 258 extreme emission line galaxies with follow-up near-infrared spectra. This includes 227 galaxies in the redshift range 1.3 < z < 2.4; the remaining 31 galaxies are z 1 Hα emitters. Finally we measure redshifts for an additional 35 galaxies (included as filler targets) with equivalent widths below our selection threshold.

The redshift distribution of the 227 galaxies at z = 1.3 -2.4 is shown in Figure 3. The distribution peaks in the 1.55 ≤ z ≤ 1.70 redshift bin, with 85 of the 227 targets falling in this window. Roughly 29% of the sample (66/227) is at z > 2, and an additional 21% (48/227) is at 1.37 < z < 1.55. As described above, this latter subset will require optical spectroscopic follow-up for measurement of [O II] and [Ne III]. Our ultimate goal is to provide spectra coverage between [O II] and Hα for the majority of extreme line emitters in this sample. In our current sample, we have obtained complete coverage for 53 of the 227 galaxies. The subset with spectral coverage between Hβ and Hα is larger (73), with 64 of the 73 sources having robust detections (S/N > 3) of both Hα and Hβ. The completeness of our survey and the redshift distribution will change as we acquire more near-infrared and optical spectra in the future.

Line fluxes are determined from fits to the extracted 1D spectra. We first fit the [O III]λ5007 emission line (or Hα in the case that [O III] is not available) with a single Gaussian function. The central wavelength from the fit is then used to compute the redshift of the object, which we in turn use to identify the other emission lines. In cases where the lines are well-measured (i.e., S/N > 7), we derive the flux using a Gaussian fit to the line profile; otherwise we derive the line flux using direct integration. In cases where nearby lines are partially resolved by MMIRS or MOSFIRE, we use multiple Gaussian functions to fit the data. If the flux is measured with S/N < 2, we consider the line undetected and derive a 2σ upper limit. We correct the Hα and Hβ fluxes for Balmer absorption using the best-fitting stellar population synthesis models (described in §3). We find that median correction for Balmer lines is less than 10% of the measured emission line fluxes.

The [O III]λ5007 emission line fluxes for the extreme line emitters range between 1.1×10 -17 erg cm -2 s -1 and 6.3×10 -16 erg cm -2 s -1 with a median of 7.5 × 10 -17 erg cm -2 s -1 . The range of [O III] fluxes is very similar to those of galaxies in the MOSDEF survey (Kriek et al. 2015); while the equivalent widths in our sample are larger, the continuum magnitudes tend to be fainter than typical MOSDEF sources. The Hβ fluxes of our galaxies are fainter than [O III], with a median of 1.9×10 -17 erg cm -2 s -1 . The emission line fluxes we derive from our ground-based spectra generally agree with the WFC3 grism spectra measurements. We find a median offset of 7% and a scatter of 36%, in agreement with the comparison between MOSFIRE spectra and 3D-HST spectra reported by the MOSDEF survey (Kriek et al. 2015).

We estimate the nebular attenuation by comparing the measured Balmer decrement (i.e., the observed Hα/Hβ intensity ratio) with the value expected by Case B recombination in the case of zero dust reddening, Hα/Hβ = 2.86 (assuming electron temperature T e = 10, 000 K; [START_REF] Osterbrock | Astrophysics of gaseous nebulae and active galactic nuclei[END_REF]. To facilitate comparison with recent nearinfrared spectroscopic studies of z ∼ 2 galaxies (e.g. Reddy et al. 2015;Steidel et al. 2016;Shivaei et al. 2018), we assume the Galactic extinction curve in Cardelli et al. (1989) to compute the dust reddening toward H II regions. We will discuss the extinction values implied by this analysis in §4.1.

To investigate how the average spectral properties vary with optical line equivalent width, we create composite spectra by stacking galaxies in four bins of [O III] 4). To minimize the contribution of sky lines from individual spectra, we stack spectra using inverse-variance weighted luminosities in each wavelength bin (i.e., weighted by 1/σ 2 i , where σ 2 i is the variation of the ith individual spectrum). To check that the resulting composites are not dominated by a few individual spectra, we also create composites using uniform weights and verify that the resulting spectra have similar line measurements as the weighted composite.

The composite spectra for all the four [O III]λ5007 EW bins are shown in Figure 4. For these spectra, we have included the subset of our current sample with complete coverage between [O II] and Hα. There are 7, 26, 13, and 14 galaxies included in the 0 -225 Å, 225 -450 Å, 450 -800 Å, and 800 -2500 Å stacks shown in this figure, respectively. We measure the luminosities of the strong rest-frame optical lines ([O II], [Ne III], Hβ, [O III], and Hα) in each of the four composites. In particular, we are interested in deriving the dependence of the Balmer decrement, and the O32 and Ne3O2 indices on the [O III]λ5007 equivalent width. Each line ratio requires slightly different spectral coverage for a robust measurement. For example, the Balmer decrement requires spectra that span between Hβ and Hα, whereas the Accordingly, we can include more objects in the stacks used to investigate the dependence of the Balmer decrement on the [O III]λ5007 EW. We will discuss the resulting line ratios in §4. The detailed spectral coverage requirements and the resulting number of objects included in the composite are summarized in Table 4.

PHOTOIONIZATION MODELING

We fit the broadband fluxes of galaxies in our spectroscopic sample using the Bayesian galaxy SED modeling and interpreting tool BEAGLE (for BayEsian Analysis of GaLaxy sEds, version 0.10.5; Chevallard & Charlot 2016), which incorporates in a consistent way the production of radiation from stars and its transfer through the interstellar and intergalactic media. Broadband photometry is obtained from 3D-HST using the Skelton et al. (2014) catalog, and we utilize multi-wavelength data covering 0.3 -2.5 µm. We also test the impact of adding Spitzer/IRAC constraints for the subset of galaxies that are not strongly confused. For each object, we remove fluxes in filters that lie blueward of Lyα to avoid introducing uncertain contributions from Lyα emission and Lyα forest absorption. We also simultaneously fit the available strong rest-frame optical emission line fluxes explore the impact of C/O variations on the derived gas properties in a follow-up paper focused more closely on the photoionization model results. To account for the effect of dust attenuation, we first assume the Calzetti et al. (2000) extinction curve. We also fit galaxies assuming the Small Magellanic Cloud (SMC) extinction curve in Pei (1992). Finally, we adopt the prescription of Inoue et al. (2014) to include the absorption of IGM.

We assume constant star formation history for model galaxies in BEAGLE, and parameterize the maximum age of stars in a model galaxy in the range from 1 Myr to the age of the universe at the given redshift. We fix the redshift of each object to the spectroscopic redshift measured from the MMIRS or MOSFIRE spectra. We adopt a standard Chabrier (2003) initial mass function and assume that all stars in a given galaxy have the same metallicity, in the range -2.2 ≤ log (Z/Z ) ≤ 0.25, where BEAGLE uses the solar metallicity value Z = 0.01524 from Caffau et al. (2011). The interstellar metallicity is assumed to be the same as the stellar metallicity (Z ISM = Z ) for each object. The ionization parameter is allowed to freely vary in the range -4.0 ≤ log U S ≤ -1.0, and the dust-to-metal mass ratio is allowed to span the range ξ d = 0.1-0.5. We adopt an exponential distribution prior on the V-band dust attenuation optical depths, fixing the fraction of attenuation optical depth arising from dust in the diffuse ISM to µ = 0.4 (see Chevallard & Charlot 2016). With the above parameterization, we use the BEAGLE tool to fit the broadband SEDs and available emission line constraints for the z = 1.3 -2.4 galaxies in our sample. Emission line fluxes and broadband fluxes are put on the same absolute scale using the aperture correction procedures described in §2.2. We obtain the output posterior probability distributions of the free parameters described above, and those of other derived physical parameters such as the ionizing photon production efficiency inferred from model (ξ ion ). We use the posterior median value as the bestfitting value of each parameter. In Figure 1, we overlay the best-fitting BEAGLE models on the broadband SEDs. It is clear that the models generally do a good job recovering the shape of the continuum and the large flux excesses caused by nebular emission.

The distribution of stellar masses implied by the BEA-GLE models is shown in the left panel of Figure 5. Here we include the 64 spectroscopically-confirmed extreme [O III] emitters with spectral coverage between Hβ and Hα, and with significant detections of Hβ and Hα (S/N > 3). The median stellar mass in this sample is 4.9×10 8 M , well below the typical stellar masses ( 10 10 M ) found for z 2.3 galaxies in the KBSS and MOSDEF surveys (e.g. Strom et al. 2017;Sanders et al. 2018). We find that galaxies in our sample with the largest [O III] equivalent widths have the lowest stellar masses (see the upper left panel of Figure 6). Among the subset of sources with EW [OIII]λ5007 in excess of 800 Å, the median stellar mass is just 4.5 × 10 7 M . This increases to 4.7 × 10 8 M for sources with 450 Å < EW [OIII]λ5007 < 800 Å and 1.4 × 10 9 M for 225 Å < EW [OIII]λ5007 < 450 Å. This variation is to be expected. The near-infrared filters which are sensitive to the stellar mass are heavily contaminated by nebular emission (lines and continuum) in the highest equivalent width systems, implying a smaller contribution from stellar continuum for a given near-infrared magnitude. This effect is further amplified by the dependence of the stellar mass to light ratio on the age of the stellar population. We note that our results do not change significantly when we add Spitzer/IRAC constraints. In particular, we find that the derived stellar masses are always within 0.1 dex of the values derived without IRAC constraints.

In the middle and right panels of Figure 5, we show the best-fitting model SFR and sSFR of the same 64 line emitters described above. The median of the SFR distribution is 7.8 M yr -1 . In contrast to the stellar mass, we do not find strong variations in the SFR with [O III] EW. The median sSFR is 17 Gyr -1 , well above the average values ( 2 Gyr -1 ) for galaxies in the KBSS and MOSDEF surveys (e.g. Strom et al. 2017;Sanders et al. 2018). As expected, the sSFR increases with [O III] EW within our sample. The median sSFR ranges from 7.9 Gyr -1 (225 Å < EW [OIII]λ5007 < 450 Å) to 21 Gyr -1 (450 Å < EW [OIII]λ5007 < 800 Å) to 125 Gyr -1 (800 Å < EW [OIII]λ5007 < 2500 Å). The strong variation in sSFR seen in Figure 6 can equivalently be described as a trend in the luminosity weighted age of the stellar population. Under our assumed constant star formation history, the BEAGLE models predict that the median of the maximum stellar age [O III] emitters with spectral coverage between Hβ and Hα). Stellar masses, SFRs, and sSFRs are derived from stellar population synthesis modeling using BEAGLE. The median stellar mass is 4.9 × 10 8 M , the median of SFR is 7.8 M yr -1 , and the median of sSFR is 17 Gyr -1 .

parameter is 130, 50, and 8 Myr for the three [O III]λ5007 equivalent width bins described above. These very young values only refer to the burst that is currently dominating the observed SED and do not negate the presence of faint older stars from earlier activity.

We re-calculate equivalent widths using the newlyobtained near-infrared spectra together with the underlying continuum predicted by BEAGLE. The models provide an improved determination of the continuum, accounting for the emission line contamination in the near-infrared broadband filters. The ground-based spectra provide higher S/N detections of the fainter lines, while also allowing us to separate doublets that are blended in the grism spectra. We calculate the rest-frame equivalent widths of [O II], Hβ, [O III]λ5007, and Hα. We will use these measurements when investigating trends with rest-frame optical equivalent widths. In Table 5, we report the [O III]λ5007 and Hα equivalent widths of the 64 extreme line emitters (EW [OIII]λ5007 > 225 Å) at z = 1.3 -2.4 with robust (S/N >3) Hα and Hβ detections. We will release the full optical line properties of our sample in a catalog paper after the survey is completed.

The BEAGLE models allow us to characterize the production efficiency of hydrogen ionizing photons in the EELGs. There are various definitions of this quantity in the literature. We follow the nomenclature used in our earlier work (Chevallard et al. 2018) which we briefly review below. First, we define ξ ion as the hydrogen ionizing photon production rate (N(H 0 )) per unit UV stellar continuum luminosity (L UV ). The quantity L UV is the luminosity produced by the stellar population prior to its transmission through the gas and dust in the galaxy, and as such it does not include nebular continuum emission or absorption from the ISM. As is commonplace, we evaluate the UV luminosity at a rest-frame wavelength of 1500 Å. BEAGLE includes a determination of ξ ion in its output parameter file for each source. Secondly, we define ξ ion , as the hydrogen ionizing photon production rate per unit observed UV luminosity (L UV ), again evaluated at 1500 Å. L UV includes emission from the stellar population and nebular continuum and is not corrected for the attenuation provided by the ISM. Finally we define ξ HII ion as the hydrogen ionizing photon production rate per unit L HII UV , the observed UV luminosity at 1500 Å (including nebular and stellar continuum) corrected for dust attenuation from the diffuse ISM. This is the most commonly used definition of the hydrogen ionizing production efficiency in the literature, and we will focus primarily on this quantity in the following section.

In deriving ξ HII ion , we follow a similar procedure adopted in several recent studies (Matthee et al. 2017;Shivaei et al. 2018). We compute the Hα luminosity from our spectra and apply a correction for dust attenuation using the measured Balmer decrement (see §2.5). We then calculate the ionizing photon production rate from the Hα luminosity [START_REF] Osterbrock | Astrophysics of gaseous nebulae and active galactic nuclei[END_REF]:

L(Hα) [erg s -1 ] = 1.36 × 10 -12 N(H 0 ) [s -1 ].
(1)

This assumes radiation-bounded nebula, with negligible escape of ionizing radiation. The dust correction through the Balmer decrement traces dust only outside the H II regions, and does not account for the absorption of ionizing photons before they ionize hydrogen (Petrosian et al. 1972;Mathis 1986;Charlot & Fall 2000). Thus, the N(H 0 ) computed in Eq. ( 1) is somewhat lower than the true production rate of ionizing photons emitted by stars. However, in §4.1 we demonstrate that the extreme [O III] emitters in our sample have very little dust, indicating that the effect of absorption of ionizing photons inside the H II regions should be very small. We infer the observed UV continuum luminosity from the best-fitting BEAGLE model using a flat 100 Å filter centered at 1500 Å. We then apply a dust correction using the reddening inferred from BEAGLE assuming first a Calzetti and then a SMC extinction law. Finally, the ionizing photon production efficiency ξ HII ion is computed as follows:

ξ HII ion [erg -1 Hz] = N(H 0 ) [s -1 ] L HII UV [erg s -1 Hz -1 ]
.

(2)

In the following section, we will investigate how the ionizing production efficiency varies with the nebular line equivalent widths. Here we attempt to build some basic physical intuition about this relationship. Since the [ the stellar population 6 , the variation of ξ ion will mirror the time evolution of N(H 0 ) and L UV , the former powered by O 6 The Hα and [O III] EWs are additionally regulated by the stellar metallicity, and the [O III] EW will also be affected by the gas properties. Here we consider only the impact of age on the Hα EW at fixed metallicity with the goal of building a basic understanding of the dependence of ξ ion on the optical line equivalent widths.

stars and the latter by O to early B stars. For the very young stellar populations probed by the most extreme line emitting galaxies in our sample (EW Hα > 1000 Å), the ratio of O to B stars will be maximized, resulting in very efficient ionizing photon production. Among more moderate line emitters (EW Hα = 400 -600 Å), a larger population of B stars will have emerged, boosting L UV relative to N(H 0 ). The ionizing production efficiency will thus be reduced for these systems.

Finally, at yet lower equivalent width (EW Hα = 50 -200 Å), the O and B star populations will be closer to equilibrium, with the number of newly-formed stars nearly balanced by those exiting the main sequence. As a result, both N(H 0 ) and L UV will not vary significantly with age, and the production efficiency of ionizing photons should thus begin to plateau to a near-constant value for galaxies with EW Hα < 200 Å. While this basic physical picture guides our expectations, the precise dependence of ξ ion on the [O III] and Hα equivalent width depends on unknown physics (i.e., binary stars, rotation) governing the formation of hot stars. Empirical constraints on ξ ion are thus critical for assessing the ionizing output of early galaxies.

REST-FRAME OPTICAL SPECTROSCOPIC PROPERTIES OF EXTREME [O III]

EMITTERS AT z = 1.3 -2.4

We now investigate what the rest-frame optical spectra reveal about the radiation field and gas conditions of galaxies with large equivalent width rest-frame optical nebular line emission. We will focus primarily on empirical trends that can be extracted from line ratios, leaving a more detailed analysis of the photoionization models to a later paper. In §4.1, we use measurements of the Balmer decrements to constrain the nebular attenuation. We then characterize the radiation field and gas conditions through measurement of the ionizing photon production efficiency ( §4.2) and standard ionization-sensitive emission line ratios ( §4.3).

Balmer Decrement Measurements

The physical interpretation of our spectra requires robust determination of the nebular attenuation between [O II] and Hα. As introduced in §2.5, this is most commonly done through comparison of the Balmer decrement, I(Hα)/I(Hβ), to the ratio expected in absence of dust. As noted in §3, the Balmer decrement corrects for dust outside the H II regions but does not account for the possible absorption of LyC photons by dust before they ionize hydrogen. The first statistical measurements of the Balmer decrement distribution at high redshift have begun to emerge from the KBSS and MOS-DEF surveys (e.g. Reddy et al. 2015;Steidel et al. 2016). These investigations provide an important baseline for comparison to the extreme emission line galaxies in our sample, so we briefly review the results below. To enable comparison with investigations of nebular attenuation in the KBSS and MOSDEF galaxies, here we adopt a Cardelli et al. (1989) extinction curve.

The average nebular properties of the KBSS survey are well described by the KBSS-LM1 composite spectrum, a weighted average of galaxies at z 2.11 -2.57 with median stellar mass of 6×10 9 M and specific star formation rate of 3.5 Gyr -1 (Steidel et al. 2016) The Balmer decrement derived from the composite is I(Hα)/I(Hβ) = 3.61±0.07 (Steidel et al. 2016), similar to that found in individual KBSS galaxies (Strom et al. 2017). For the MOSDEF survey, the average Balmer decrement measured from a composite of 213 galaxies with significant Hβ detections is I(Hα)/I(Hβ) = 4.1 (Reddy et al. 2015). For a Cardelli et al. (1989) I(Hα)/I(Hβ) = 2.86, the KBSS and MOSDEF composites imply typical color excesses of E(B-V) gas = 0.22 and 0.36, respectively. We note that Steidel et al. (2016) adopt a slightly larger intrinsic Balmer decrement, I(Hα)/I(Hβ) = 2.89, that is consistent with their photoionization model fits. As a result, they report a slightly different color excess, E(B-V) gas = 0.21, for the KBSS-LM1 composite.

In our current sample, we can measure Balmer decrements in 64 extreme emission line galaxies (EW [OIII]λ5007 > 225 Å) with significant detections (S/N > 3) of Hα and Hβ (Table 4). The measurements suggest that extreme line emitters suffer less nebular attenuation than most galaxies in the KBSS and MOSDEF surveys. In Figure 7, we show the distribution of Balmer decrements for EELGs in our sample with EW [OIII]λ5007 > 450 Å. This threshold is chosen to correspond roughly to the average [O III]+Hβ equivalent width at z 7 -8 (L13), assuming an [O III]/Hβ ratio that is characteristic of EELGs (see §2.1). The average Balmer decrement in this histogram is I(Hα)/I(Hβ) = 3.09, implying a color excess of just E(B -V) gas = 0.08 for the Cardelli et al. (1989) extinction curve. For reference we also measure β, the UV continuum slope, for this same subset of galaxies. We compute the UV slope from broadband fluxes by fitting a power-law ( f λ ∝ λ β ) at rest-frame wavelengths 1250 -2600 Å (the same wavelength range used in Calzetti et al. 1994; filters covering Lyα emission line are not included). The resulting UV slopes are very blue, with a median of β = -2.04 for galaxies with EW [OIII]λ5007 > 450 Å, implying that the stars are also minimally reddened by dust.

The level of nebular attenuation varies with the [O III] EW, with the most extreme line emitters having the least dust. In our sample, the median Balmer decrement for the 16 galaxies with EW [OIII]λ5007 = 450-800 Å is I(Hα)/I(Hβ) = 3.12. When we consider the 21 more extreme line emitters with EW [OIII]λ5007 = 800 -2500 Å, we find a Balmer decrement of just I(Hα)/I(Hβ) = 2.95. This implies color excesses of E(B -V) gas = 0.09 and 0.03 for the two respective EW bins. The composite spectra reveal a similar picture. The measured Balmer decrement in the stacks decreases from I(Hα)/I(Hβ) = 3.12 ± 0.04 (EW [OIII]λ5007 = 450 -800 Å) to I(Hα)/I(Hβ) = 2.90 ± 0.05 (EW [OIII]λ5007 = 800 -2500 Å). This indicates that among the most extreme line emitters, the Balmer decrements are very close to the intrinsic case B recombination value (for T e = 10 4 K gas), suggesting that the emission lines in these galaxies face little to no attenuation from dust.

The results described above clearly demonstrate that the EELG population has a different distribution of Balmer decrements than the typical KBSS and MOSDEF galaxies. This result is not surprising given trends between the Balmer decrement and galaxy properties found previously (Reddy et al. 2015), likely reflecting both the low stellar mass (and hence moderately low metallicities) and large specific star formation rates (and hence young stellar populations) of the galaxies in our sample (Figure 6). Importantly this implies that very small adjustments are required to correct the observed fluxes for reddening.

An important consequence is that uncertainties in the high redshift attenuation law should not significantly affect our interpretation of the nebular line spectra of EELGs. Among more massive star forming galaxies at z 2, this is not always the case. Shivaei et al. (2018) demonstrated that the ξ HII ion inferred assuming a Calzetti attenuation law is systematically 0.3 dex lower than that derived for an SMC attenuation law, making it difficult to robustly determine ξ HII ion . While our measurements of ξ HII ion also depend on the assumed dust law, the variation is typically only 0.1 dex when considering the SMC and Calzetti attenuation laws. For the most extreme line emitting sources (EW [OIII]λ5007 > 800 Å), the dust content is low enough that the average offset in ξ HII ion is just 0.05 dex.

The Ionizing Photon Production Efficiency

As a first step in our investigation of the radiation field of extreme emission line galaxies, we characterize the production efficiency of hydrogen ionizing photons in our spectroscopic sample. Recent efforts have quantified ξ HII ion in more massive star forming galaxies at z 2 (e.g Matthee et al. 2017;Shivaei et al. 2018), revealing typical values of log [ξ HII ion (erg -1 Hz)] = 25.06 for a Calzetti UV attenuation law (Shivaei et al. 2018). These systems tend to have [O III]λ5007 equivalent widths of 115 Å (see §2.1), well below those considered in this paper, as expected from older stellar populations. It has been shown that ξ ion scales with the [O III] equivalent width among local star forming galaxies (Chevallard et al. 2018), suggesting that the higher equivalent width systems which become common in the reionization era may be more efficient ionizing agents. Here we seek to investigate whether a similar relationship between the ionizing production efficiency and optical lines holds at z 1 -2. In the following, we will first describe the distri-bution of ξ HII ion values, using the dust corrections discussed in §3. Then we will compare our measurements of ξ ion to the relation found locally, investigating whether there is any evidence for strong redshift evolution in ξ ion values at fixed [O III To estimate the typical ionizing efficiency for galaxies with this [O III] EW, we group those systems in our sample with EW [OIII]λ5007 = 300-600 Å. Assuming a Calzetti attenuation law for the UV, the median ionizing production efficiency of this sub-sample is log [ξ HII ion (erg -1 Hz)] = 25.22, ∼ 1.5× greater than found in more typical galaxies at z 2 (Shivaei et al. 2018). For an SMC extinction law for the UV, the ionizing production efficiency is log [ξ HII ion (erg -1 Hz)] = 25.32, slightly larger than what we found above for the Calzetti law. This suggests that galaxies which are common at z 6 should make a greater contribution to reionization than was originally thought, as has been noted by several recent investigations of IRAC excesses and UV metal lines at z 5 (Rasappu et al. 2016;Stark et al. 2015bStark et al. , 2017)).

Attention is now focused on galaxies at z 7 with even larger [O III] equivalent widths (EW [OIII]λ5007 = 800 -2500 Å). While such sources may not be the norm at z 7, they are readily detectable in the CANDELS fields (e.g. Roberts-Borsani et al. 2016) and appear to have an enhanced Lyα visibility at redshifts where the IGM is thought to be significantly neutral (e.g. Oesch et al. 2015;Zitrin et al. 2015;Stark et al. 2017). One possible explanation is that these systems are situated in the largest ionized bubbles, thereby allowing enhanced transmission of their Lyα radiation. Alternatively, in addition to sitting in ionized bubbles, the galaxies may be producing more Lyα than other galaxies of similar far-UV luminosities. This would follow if the production efficiency of Lyman-alpha is larger in systems with the most prominent [O III] emission.

Our results provide some insight into the situation. Since Lyα is produced by reprocessed hydrogen ionizing photons, the Lyα production efficiency (L Lyα /L HII UV ) should be directly correlated with ξ HII ion . Sources that are efficient at producing ionizing radiation should also be efficient in powering Lyα emission. In Figure 8, we see that the production efficiency of ionizing radiation steadily increases between EW [OIII]λ5007 = 450 Å and 2500 Å. Among those sources with EW [OIII]λ5007 > 800 Å, the production efficiency has a median of log [ξ HII ion (erg -1 Hz)] = 25.58 (Calzetti) and 25.62 (SMC). This value is very similar to estimates of ξ HII ion in the RB16 sample of z 7 galaxies with similarly strong [O III] (see red squares in Figure 8). Our results suggest that galaxies with [O III] EWs similar to the RB16 galaxies may produce two times more ionizing photons than EW for the z = 1.3 -2.4 galaxies in our sample (blue circles). The ξ HII ion is computed based on the spectroscopic measurements of Hα and Hβ. A linear relation is fitted between ξ HII ion and log EW [OIII]λ5007 (black dashed line, with grey area shows the 95 percent credible region). The average EW of z 7 -8 galaxies (EW [OIII]λ5007 450 Å; L13) is highlighted by yellow dashed line. The orange area shows EWs of z > 7 galaxies in RB16, and red squares show the individual sources in RB16, in which the ξ HII ion are inferred from photoionization models (Stark et al. 2017). Right Panel: Relation between ξ HII ion and Hα EW for the z ∼ 2 galaxies in our sample (blue circles). A linear relation is fitted between ξ HII ion and log EW Hα (black dashed line, with grey area shows the 95 percent credible region).

typical reionization-era sources (the latter of which have EW [OIII]λ5007 450 Å) with similar non-ionizing UV luminosities. The larger-than-average Lyα equivalent widths that we are now seeing in the RB16 population should thus partially reflect an enhanced production of ionizing radiation. In §5.2, we will explore in more detail whether variations in internal galaxy properties can explain the anomalous visibility of Lyα that is now being found in strong [O III] emitters at z 7.

The relationship between the production efficiency of ionizing photons and the [O III]λ5007 EW shown in Figure 8 can be described by a simple scaling law. We fit a linear relation between the two quantities for those galaxies in our sample with 225 Å < EW [OIII]λ5007 < 2500 Å. The best-fit relation for the Calzetti UV attenuation law is

log ξ HII ion = (0.76 ± 0.05) × log (EW [O III]λ5007) + (23.27 ± 0.15). (3)
For an SMC attenuation law, we find a very similar relationship:

log ξ HII ion = (0.73 ± 0.08) × log (EW [O III]λ5007) + (23.45 ± 0.23). (4)
The derived fitting function is overlaid on the data in Figure 8. We emphasize that the scaling laws is only valid for the large equivalent widths stated above. Indeed, as we motivated at the end of §3, we expect ξ HII ion may deviate from the relationship at lower equivalent widths. But for sources with EW [OIII]λ5007 > 225 Å, these fitting functions can be used to predict the distribution of ξ HII ion given an observed distribution of [O III] equivalent widths. As suggested in Chevallard et al. (2018), and provided that this relationship does not evolve strongly with redshift (as we will show below for 0 < z < 2), this can be used in conjunction with IRAC flux excesses to estimate the production efficiency of ionizing photons in the reionization era.

We can also derive a relationship between Hα and ξ HII ion using the results in Figure 8. Because the Hα strength is less sensitive to the gas physical conditions than [O III], we expect it to have a more universal relationship with the ionizing production efficiency. Applying a similar fitting procedure and assuming a Calzetti attenuation law for the UV, we derive the following relationship:

log ξ HII ion = (0.85 ± 0.06) × log (EW Hα) + (23.03 ± 0.15). (5)
For the SMC law, we again find a very similar relationship:

log ξ HII ion = (0.87 ± 0.07) × log (EW Hα) + (23.06 ± 0.19). (6)
To assess the magnitude of the redshift evolution in the relationship between the production efficiency and the restframe optical line equivalent widths, we compare our results to those derived for local galaxies in Chevallard et al. (2018). The local galaxy relation is derived using ξ ion , the ratio of the intrinsic production rate of ionizing photons to the stellar continuum luminosity in the non-ionizing far-UV, with both quantities derived from BEAGLE model fits to the data (see definitions in §3). In cases where the nebular continuum contributes significantly to the observed UV luminosity, we expect ξ ion to deviate from ξ HII ion , providing a more realistic description of the efficiency of a stellar population at producing ionizing radiation. We find that ξ ion is indeed system-atically larger than ξ HII ion in our sample, but the differences are relatively small. Among galaxies with [O III]λ5007 equivalent widths between 450 and 800 Å (800 and 2500 Å), the difference is 0.04 (0.08) dex. This suggests that the relations between ξ HII ion and optical line equivalent widths that we derived above should closely approximate those inferred using ξ ion . Nevertheless for the sake of consistency, we consider trends with ξ ion when comparing the production efficiencies derived at z = 1.3 -2.4 with those found locally.

In Figure 9, we show the dependence of ξ ion on the Hβ and [O III] EW for our high redshift sample (blue circles) and the nearby galaxies (green squares) investigated in Chevallard et al. (2018). The local sample is comprised of ten systems with deep UV and optical spectra (Senchyna et al. 2017), each selected to have extreme radiation fields based on the presence of He II emission. The relationship between ξ ion and Hβ EW appears largely similar in the two redshift samples. We can quantify this by comparing the linear fit to the log ξ ionlog EW Hβ relationship in the range EW Hβ = 80 -400 Å. The best-fit slope and intercept for the z 0 sample (1.23 ± 0.25 and 22.89 ± 0.56) are both consistent (within 1σ) with those derived for the z 2 galaxies (1.09 ± 0.07 and 23.23 ± 0.16). Given the dependence of [O III] on ionized gas conditions (i.e., metallicity, ionization parameter), we expect the relationship between [O III] EW and ξ ion may potentially show more variation with redshift. However, the current data reveal little evidence to this effect. The best fit intercept and slope at z 0 (1.14 ± 0.16 and 22.28 ± 0.44) remain broadly similar (within 2σ) to those derived at z 2 (0.88±0.06 and 22.98±0.16). While larger samples may eventually reveal some mild differences, our results suggest that the redshift evolution in the relationship between ξ ion and the Hβ and [O III] EW is not likely to be strong.

The Physical Conditions of the Nebular Gas

The results of §4.2 show that extreme emission line galaxies are considerably more efficient at producing ionizing radiation than the more massive star-forming galaxies which are typical at z ∼ 2. The intense radiation field of these system may impact the ionization state of the ISM, potentially aiding the escape of Lyα and LyC radiation. In this subsection, we investigate the ISM of extreme emission line galaxies by quantifying the dependence of The O32 index is one of the most commonly used probes of nebular gas ionization state and is often employed as an empirical proxy for the ionization parameter (the ratio of the number density of incident hydrogen-ionizing photons to the number density of hydrogen atoms in the H II region; Penston et al. 1990). The average O32 ratios of the massive star forming galaxies in the KBSS (O32 = 2.0; Steidel et al. 2016) and MOSDEF (median O32 = 1.3; Sanders et al. 2016) surveys are ∼ 4 -7× higher than those of z ∼ 0 galaxies with similar stellar masses (median O32 = 0.3; Abazajian et al. 2009;Sanders et al. 2016). While the origin of the higher ionization parameters remains a matter of some debate, it has been argued to reflect a shift toward lower metallicities at fixed mass in high redshift galaxies (Sanders et al. 2016).

The extreme emission line galaxies targeted in our sur-vey have lower stellar masses and higher specific star formation rates than most galaxies from the KBSS and MOSDEF surveys (see §3). Both factors could lead to different ionization conditions. Our current sample contains O32 measurements with dust corrections via the Balmer decrement for 44 galaxies (Table 4) with large equivalent width We see a comparable O32 trend in the individual galaxy measurements. For galaxies with EW [OIII]λ5007 = 300 -600 Å (similar to the average EW expected in the reionization era), we find a median [O III]/[O II] ratio of O32 = 3.5. Among the galaxies in our sample with [O III] The Ne3O2 index is another proxy for the ionization parameter (e.g. Levesque & Richardson 2014), providing an independent probe of the ionization state of the ISM. Owing to the short wavelength baseline between [Ne III] and [O II], the flux ratio is largely insensitive to the effect of reddening. At the highest redshifts probed by JWST (z 10), [Ne III] and [O II] will be among the brightest lines visible in the NIRSpec bandpass (see §5.3 for more discussion), motivating considerable interest in the effectiveness of Ne3O2 in proving nebular gas conditions. Among the more massive star forming galaxies probed by the KBSS-LM1 composite spectrum, [Ne III] tends to be much fainter than [O II], with Ne3O2 = 0.15 (Steidel et al. 2016). Based on the dependence of O32 on [O III] EW described above, we expect to see larger Ne3O2 values among the extreme emission line population.

In our spectroscopic sample, we have obtained measurements of Ne3O2 for 26 galaxies (Table 4) with large equivalent width [O III] emission (EW [OIII]λ5007 > 225 Å). In the bottom panel of Figure 10, we plot the dependence of Ne3O2 on the [O III]λ5007 EW, showing individual galaxies (blue circles) and measurements from the composite spectra (red squares). The results show that the Ne3O2 index increases with the [O III] EW (albeit with scatter), reaching values that are much greater than found among more massive galaxies in the KBSS survey. In our composite spectra, we measure Ne3O2 = 0.17 ± 0.06, 0.20 ± 0.02, 0.26 ± 0.03, and 0.51 ± 0.04 for stacks including galaxies with EW ters in our sample, the [Ne III] line is just as strong as the individual components [O II] doublet. In §5.3, we will consider the feasibility of detecting both features at z 10 with JWST.

Both the O32 and Ne3O2 measurements suggest a picture in which the ISM of galaxies with prominent optical line emission is characterized by extreme ionization conditions. A large number of factors can modulate O32 (e.g. Nakajima & Ouchi 2014). We will discuss some of these in §5.1, but a detailed investigation of the physical origin of the trend between O32 and [O III] EW will be considered in a followup paper focused on the photoionization models described in §3.

DISCUSSION

In the previous section, we showed that the ionizing photon production efficiency and the ionization state of the nebular gas scales with optical line equivalent width. Here we consider implications for the escape of ionizing radiation (and its potential association with large O32), the Lyα visibility test, and observability of extreme line emitting galaxies with JWST at z 10. Over the past few years, attention has focused on a possible connection between large values of O32 and the escape of ionizing radiation (e.g. Izotov et al. 2016Izotov et al. , 2017Izotov et al. , 2018;;Naidu et al. 2018;Fletcher et al. 2018). It has been argued that large [O III]/[O II] ratios may arise in density-bounded nebula in which ionizing radiation is able to escape (Guseva et al. 2004;Jaskot & Oey 2013;Nakajima & Ouchi 2014), although as we note below this is not a unique interpretation. An investigation of five compact z 0.3 Green Peas with relatively large [O III]/[O II] ratios (O32 = 6.4 -8.9 with a median of O32 = 6.5) revealed substantial ionizing photon escape fractions ( f esc = 0.06-0.13) in all five systems (Izotov et al. 2016). More recent work has shown even larger escape fractions ( f esc > 0.3) in several galaxies with O32 > 10, while also demonstrating that there is considerable scatter in LyC escape fractions at large O32 (e.g. Izotov et al. 2018). Large O32 may be a necessary but not sufficient condition for the detection of LyC emission in low redshift samples (Izotov et al. 2017). Work at higher redshift indicates a similar picture, with the small number of known LyC leakers typically exhibiting large [O III]/ [O II] ratios (e.g. Vanzella et al. 2016;de Barros et al. 2016;Fletcher et al. 2018). High redshift galaxies with moderate [O III]/[O II] ratios (O32 > 3) tend not to show LyC detections (Naidu et al. 2018), implying relatively low escape fractions ( f esc < 0.1). This is consistent with the local results indicating the largest escape fractions are generally found in galaxies with the largest O32 values. We emphasize that the relationship between O32 and LyC leakage does not uniquely point to density-bounded nebula, as stellar ionizing spectra are capable of producing the majority of O32 values seen in these studies (e.g. Stasińska et al. 2015). But it does suggest that the ISM of galaxies with large O32 ( 6.5) is conducive to the escape of ionizing radiation.

While we do not have direct constraints on the escape fraction of the galaxies in our sample, the results in §4 show that not all extreme emission line galaxies power the O32 ratios that appear to be associated with large escape fractions at lower redshifts. The extreme [O III]/[O II] ratios (O32 > 10) associated with many LyC leakers ( f esc 0.3) become the norm for galaxies with extremely strong [O III] emission (EW [OIII]λ5007 > 1100 Å), implying very large specific star formation rates ( 100 Gyr -1 ) and young stellar populations (∼ 3 -10 Myr for constant star formation; see §3). The O32 values that appear to be associated with more moderate LyC escape fractions (O32 6.5, f esc 0.1) map to galaxies with EW [OIII]λ5007 760 Å, corresponding to (constant star formation) ages of 10 -30 Myr. There are several important caveats to this picture. First the relation between O32 and f esc remains poorly sampled and needs additional constraints at low and high redshift. Second, significant LyC leakage has recently been reported by Shapley et al. (2016) in a high redshift galaxy that appears to be very different from what we have described above, with a comparatively old age ( 1 Gyr for constant star formation) and relatively weak [O III]+Hβ emission (EW [OIII]+Hβ < 256 Å). Clearly young stellar populations are not the only route toward LyC escape. Finally, the absolute ages we quote above are subject to standard systematic uncertainties associated with our assumed star formation history and population synthesis model (see §3) and may be somewhat larger if binary stars contribute significantly to ionizing production. Nevertheless the data suggest a relative trend whereby the conditions for LyC escape appear to be optimized in galaxies with extremely young stellar populations and intense optical line emission. This is consistent with a picture whereby LyC escape can reach large values for a short window following a burst of star formation (and after the birth cloud has been cleared), as has been predicted in many simulations (e.g. Kimm & Cen 2014;Wise et al. 2014;Ma et al. 2015;[START_REF] Osterbrock | Astrophysics of gaseous nebulae and active galactic nuclei[END_REF].

The potential relationship between LyC leakage and rest-frame optical line properties (in particular, O32 and [O III] EW) has implications for our understanding of ionizing photon escape in the reionization era. Over the past few years, attention has focused on the intense [O III] emission that is thought to become common at z 7 -8 (Labbé et al. 2013;Smit et al. 2014Smit et al. , 2015)). While extreme emission line galaxies may be the norm at these redshifts, the average equivalent width (EW [OIII]λ5007 450 Å; see §4.1) is considerably lower than that typically seen in LyC leakers. If the mapping between [O III] EW and O32 (valid at 0.3 z 2.4; Figure 10) is similar at z 6, it would imply typical [O III]/[O II] ratios of O32 3.6. As noted above, high redshift galaxies with similar optical line properties tend to have relatively low escape fractions (Naidu et al. 2018). It is of course possible that the galaxies at z 6 have different ISM conditions at fixed [O III] EW. Spectroscopic measurements with JWST will soon clarify the [O III] EW and O32 distribution at z 6, allowing constraints to be placed on the fraction of the population with rest-frame optical line properties that are correlated with significant LyC escape.

Implications for the Lyα visibility test

The Lyα visibility test provides one of the few probes of the IGM ionization state at 7 z 8. The rapid decline in the population at z 7 is consistent with a substantially neutral IGM at 7 z 8 (e.g. Mason et al. 2018a;Weinberger et al. 2018). The 100% success rate in identifying Lyα in the RB16 sample of [O III] emitters (Oesch et al. 2015;Zitrin et al. 2015;Stark et al. 2017) stands in sharp contrast to the weak Lyα that is typical of the reionization era, likely indicating that Lyα visibility is enhanced in massive galaxies with intense optical line emission. Specifically, the fraction of these sources with Lyα EW > 25 Å (x Lyα = 0.50 ± 0.29) is a factor of five greater than the general population (Stark et al. 2017). It has been suggested that these systems may trace overdense regions that produce extremely large ionized bubbles, boosting transmission of Lyα in a mostly neutral IGM. While the galaxies must lie in ionized regions for Lyα to escape, we suggest that their atypical visibility may instead be related to properties internal to the galaxies. The trend between ξ HII ion and [O III] EW (Figure 8) suggests the RB16 galaxies (EW [OIII]λ5007 > 800 Å) are likely to be twice as efficient at producing Lyα as typical galaxies (EW [OIII]λ5007 450 Å) in the reionization era. In addition to enhanced production, the transmission through the galaxy is likely to be increased due to the highly ionized state of the ISM in the most extreme line emitters (see Figure 10).

To test the impact of ξ HII ion and ISM variations on the Lyα EW distribution, we consider the Lyα properties of lower redshift samples with intense optical line emission. Stark et al. (2017) characterized the fraction of large EW (> 25 Å) Lyα emission in a sample galaxies at 3.8 < z < 5.0 with intense Hα emission. The Hα EWs are determined using Spitzer/IRAC color excesses. For the sake of comparison with the results at z 7, Stark et al. (2017) isolate luminous galaxies (M UV < -21) with blue UV slopes (β < -1.8) and an Hα threshold (EW Hα+[NII]+[SII] > 600 Å) that is is chosen to match the [O III] EW cut of the RB16 sample. The Lyα emitter fraction in this sub-sample is found to be x Lyα = 0.53 ± 0.17, nearly a factor of five greater than that of the general population of z 4 -5 galaxies of the same luminosity (x Lyα = 0.12; Stark et al. 2010). Because the IGM is highly ionized at these redshifts, the increased Lyα output that is observed should reflect factors internal to the galaxies.

A similar test can be conducted with Green Peas, the population of EELGs at z 0.3. The Lyα properties of 43 Green Peas have recently been reported by Yang et al. (2017). The [O III] emission lines in their sample are very strong, with a median EW of EW [OIII]λ5007 900 Å, comparable to the RB16 galaxies. Large EW Lyα is shown to be common in the HST/COS spectra, with 28 of 43 galaxies having Lyα EW > 25 Å. This implies a Lyα emitter fraction (x Lyα = 0.65) that is similar to that measured in the extreme line emitters at 3.8 < z < 5.0 and at z 7. The low redshift samples described above confirm that the conditions are conducive to strong Lyα emission in galaxies with intense optical line emission. The factor of five enhancement in x Lyα seen in the RB16 galaxies is fully consistent with that measured in EELGs at z 4 -5, indicating that IGM variations (i.e., larger than average ionized bubbles) are not required to explain the large Lyα success rate in this sample of galaxies. It has also been argued that the visibility is likely further enhanced in massive galaxies owing to their larger Lyα velocity offsets (e.g., Stark et al. 2017;Mason et al. 2018a). Taken together, these results suggest that massive galaxies with extremely large [O III] EW are the systems most likely to have detectable Lyα at very high redshifts where the IGM is still substantially neutral. While less massive systems will provide the most sensitive probe of the IGM in the middle and late stages of reionization (e.g., Mason et al. 2018a;Weinberger et al. 2018), their Lyα emission will be mostly extinguished at earlier epochs. Because of their prodigious Lyα output, massive extreme emission line galaxies may offer a unique path toward studying the IGM at the earliest stages of reionization (see Mason et al. 2018b).

Implications for JWST observations of z 10 galaxies

The discovery that reionization era systems tend to have strong [O III] emission is important not only for what it tells us about the nature of early galaxies, but also for what it implies about the feasibility of future spectroscopic surveys with JWST. If large equivalent line emission is commonplace, it will be straightforward to build large redshift samples down to faint continuum flux levels. But at the highest redshifts that JWST will hope to probe (z 10), the prospects for spectroscopy are not so straightforward. The Over the past few years, there has been considerable efforts invested in the UV lines, revealing equivalent widths that approach 20 Å for 10 -20% of the population (Mainali et al. 2018). In this section, we turn our attention to the EW distribution of the fainter rest-frame optical lines in the extreme emission line population and consider which features are likely to provide the best route to spectroscopic confirmation at z 9.

We first consider the feasibility of detecting -18 , 1.1×10 -18 , and 4.5 × 10 -19 erg cm -2 s -1 for a z 10 galaxy with apparent magnitudes of H = 26, 27, and 28, respectively. Here we assume a flat continuum in f ν to extrapolate to the underlying continuum near [O II]. We compute the detectability of [O II] assuming each component of the doublet has a line width typical of the EELGs in our sample (FWHM = 190 km s -1 ). For R 1000 NIRSpec MSA observations with the G395M/F290LP disperser-filter combination, the NIRSpec exposure time calculator predicts the doublet can be detected with S/N = 5 in 27 minutes (assuming underlying continuum of H = 26), 2.7 hours (H = 27), and 16 hrs (H = 28). This is considerably more efficient than targeting the C III] doublet with NIRSpec. Assuming the most optimistic value for its equivalent width (EW CIII] 20 Å), likely present in only a subset of early galaxies (Mainali et al. 2018), we predict exposure times that are 6× larger than what is required for detection of [O II].

If the ISM conditions at very high redshift are similar to what we have observed at z 0 -2, we expect [Ne III] will be 0.5× the flux of the combined [O II] doublet in the most extreme line emitters. [Ne III] will be narrower than [O II] in the R 1000 grism, allowing it to be detected in just 1 hr (H = 26) and 6 hrs (H = 27) for sources with EW [NeIII] = 50 Å. As recently pointed out in Shapley et al. (2017), deep NIRSpec integrations targeting z 10 galaxies should easily be able to recover both lines in sources with H 27, providing a unique path toward constraining the ISM conditions in the highest redshift galaxies JWST will probe.

The Hβ line will not be detectable by NIRSpec at z 10, as it shifts out of the bandpass at roughly the same redshift as [O III]. But it will serve as an important diagnostic for reionization-era systems. Since Hβ is likely to be among the faintest rest-frame optical lines targeted, it is valuable to know the range of Hβ equivalent widths expected in the extreme emission line population. In the bottom panel of Figure 11 While much attention has recently focused on the rest-UV metal lines, we show here that the fainter rest-frame optical emission lines (i.e., [O II], [Ne III], Hβ) are likely to be provide a more efficient route toward spectroscopic confirmation of EELGs at 9 z 12. For sources that are bright enough (H 26), the UV nebular lines will complement the picture provided by the rest-frame optical lines, providing unique constraints on the nature of the ionizing sources (e.g. Feltre et al. 2016;Jaskot & Ravindranath 2016;Mainali et al. 2017;Senchyna et al. 2017;Byler et al. 2018) and the abundance of carbon in the early universe (e.g. Berg et al. 2016). But initial JWST/NIRSpec efforts aimed at building a statistical sample of redshifts for z 10 galaxies are likely to be more successful targeting the observed wavelengths around [O II] and [Ne III].

SUMMARY

We report results from a large MMT and Keck spectroscopic survey targeting strong rest-frame optical emission lines in extreme emission line galaxies at high redshift. We have confirmed redshifts of 227 galaxies at z = 1.3 -2.4 with rest-frame [O III]λ5007 equivalent widths spanning 225 Å < EW [OIII]λ5007 < 2500 Å. This range is chosen to include the average [O III] EW of reionization-era galaxies (EW [OIII]λ5007 450 Å; L13) and the more extreme values (EW [OIII]λ5007 > 800 Å) that characterize most of the known Lyα emitters at z > 7. By obtaining spectra of similarly intense line emitters at slightly lower redshifts, we hope to gain insight into the range of stellar populations and gas conditions present in the reionization era as well as what distinguishes the Lyα emitting galaxies from the rest of the population at z 7.

We fit the broadband continuum SED and emission line fluxes using the BEAGLE tool (Chevallard & Charlot 2016). The stellar masses implied by the models extend down to 10 7 M with a median of 4.9 × 10 8 M , both well below typical values in the KBSS and MOSDEF surveys. The sSFR distribution extends up to 300 Gyr -1 with a median of 17 Gyr -1 , consistent with very young stellar populations expected after a substantial burst or recent upturn of star formation. We summarize our findings below:

1. The Balmer decrement decreases with increasing [O III] and Hα EW, reaching a value close to that expected for case B recombination (I(Hα)/I(Hβ) = 2.90 ± 0.05) in our composite spectrum of galaxies with EW [OIII]λ5007 = 800 -2500 Å. The UV continuum slopes are uniformly blue, with median values of β = -2.04 for galaxies with EW [OIII]λ5007 > 450 Å. This suggests that both the hot stars and nebular gas are minimally obscured by dust in the most extreme line emitters. 2. We find that the ionizing production efficiency scales with the [O III] and Hα EW at high redshift, reaching the largest values (log [ξ HII ion (erg -1 Hz)] 25.6) in the most intense emitters (EW [OIII]λ5007 > 800 Å). We derive simple scaling laws between the ionizing production efficiency and the [O III] and Hα equivalent widths. Comparison to a recent investigation of nearby galaxies (Chevallard et al. 2018) reveals that there is not strong evolution in the link between ξ ion and the optical lines over 0 z 2. The trend between the ionizing efficiency and [O III] EW indicates that the extreme line emitters which become common in the reionization era are likely to be very efficient ionizing agents.

3. We discuss the impact of the ξ ion -[O III] EW relationship on the visibility of Lyα in the reionization era. A large ionizing efficiency results in a larger than average intrinsic Lyα EW, boosting the visibility relative to systems with lower [O III] EW. The detectability is further enhanced by the ionization conditions in the ISM and CGM and the large Lyα velocity offsets associated with massive galaxies. We suggest that the anomalous Lyα success rate in massive galaxies with intense [O III] is likely driven by factors internal to the galaxies and does not necessarily require these objects to sit in the very largest ionized bubbles at these early epochs. Because of their atypical Lyα visibility, massive extreme emission line galaxies may provide a unique probe of the earliest stages of reionization (see Mason et al. 2018b) where substantial neutral fractions (and smaller velocity offsets) make lower mass galaxies very difficult to detect in Lyα (Mason et al. 2018a;Weinberger et al. 2018).

4. We show that two ionization-sensitive line ratios (O32 and Ne3O2) increase with the [O III] EW, revealing extreme ionization conditions in the ISM of the most intense line emitters. In a stacked spectrum of galaxies with EW [OIII]λ5007 = 800 -2500 Å, the ionization-sensitive ratios have very large values (O32 = 9.1 and Ne3O2 = 0.51) that are commonly associated with LyC leakers at 0 z 2. We find that very similar O32 and Ne3O2 relationships are also present in Green Pea population, indicating very little redshift evolution over 0.3 z 2.4 and suggesting a close connection between the [O III] EW and the ionization state of the ISM in the extreme line emitter population.

5. We discuss the O32 -[O III] EW relationship in the context of emerging evidence that large O32 values tend to accompany LyC leakers. Our results demonstrate that the O32 values that appear associated with moderate ( f esc 0.1) escape fractions (O32 6.5) are most commonly found in EELGs with EW [OIII]λ5007 > 760 Å. The extreme O32 (O32 > 10) values linked with larger escape fractions ( f esc 0.3) become the norm for galaxies with extremely strong [O III] emission (EW [OIII]λ5007 > 1100 Å). These [O III] EW values are well above the threshold commonly adopted for EELG selections, corresponding to extremely young stellar populations ( 3 -30 Myr for constant star formation history), consistent with a picture whereby LyC leakage is enhanced during a short window after a burst of star formation.

6. We show that [O II] is likely to be one of the efficient spectroscopic probes at very high redshifts (z 10) where the strongest optical lines are shifted out of the NIRSpec bandpass. In spite of the large O32 values in the strongest line emitters, the median [O II] doublet EW remains fairly large (EW [OII] 100 Å) throughout the EELG population. We show that [O II] is likely to be detected in considerably less exposure time than the strongest of the UV nebular lines (i.e., C III], O III]). [Ne III] should also be readily detectable if Ne3O2 ratios are as large as indicated by our survey.

  λ5007 EW: 0 -225 Å, 225 -450 Å, 450 -800 Å, and 800 -2500 Å. We first shift individual spectra to the rest frame using the redshifts measured from [O III]λ5007 emission lines. Each spectrum is then interpolated to a common rest-frame wavelength scale of 1 Å and normalized by the individual Hα luminosity (when investigating the average [Ne III]/[O II] ratio in §4.3, we normalize the individual spectrum by its [O III] luminosity since the Hα may not be observed in sources used for [Ne III]/[O II] analysis; see also Table

Figure 2 .Figure 3 .

 23 Figure 2. Examples of 1D MMT and Keck spectra of two extreme [O III] emitters at z = 1.3 -2.4. Wavelengths are shifted to the rest-frame. The upper panels show the spectrum of a galaxy with similar equivalent width as the z ∼ 7 galaxies in the RB16 sample ([O III]λ5007 EW = 814 Å). The lower panels show the spectrum of a galaxy with similar equivalent width as the average value of z ∼ 7 galaxies in L13 ([O III]λ5007 EW = 437 Å). The 1σ error spectra are also shown as grey regions. The locations of strong rest-frame optical emission lines are highlighted by black dotted lines. Zoom-in spectra of relatively faint [O II] and [Ne III] emission lines are shown in the left panels. Emission line spectrum of the EW = 814 Å galaxy shows higher [O III]/[O II] and [Ne III]/[O II] ratios than that of the galaxy with EW = 437 Å.

Figure 5 .

 5 Figure 5. Stellar mass, star formation rate, and specific star formation rate distributions for z = 1.3 -2.4 extreme [O III] emitters ([O III]λ5007 EW > 225 Å) with robust Hα and Hβ measurements (S/N > 3) in our sample (64 of total 73 extreme[O III] emitters with spectral coverage between Hβ and Hα). Stellar masses, SFRs, and sSFRs are derived from stellar population synthesis modeling using BEAGLE. The median stellar mass is 4.9 × 10 8 M , the median of SFR is 7.8 M yr -1 , and the median of sSFR is 17 Gyr -1 .

Figure 8 .

 8 Figure 8. Left Panel: Relation between the ionizing photon production efficiency ξ HII ion and [O III]λ5007EW for the z = 1.3 -2.4 galaxies in our sample (blue circles). The ξ HII ion is computed based on the spectroscopic measurements of Hα and Hβ. A linear relation is fitted between ξ HII ion and log EW [OIII]λ5007 (black dashed line, with grey area shows the 95 percent credible region). The average EW of z 7 -8 galaxies (EW [OIII]λ5007 450 Å; L13) is highlighted by yellow dashed line. The orange area shows EWs of z > 7 galaxies in RB16, and red squares show the individual sources in RB16, in which the ξ HII ion are inferred from photoionization models(Stark et al. 2017). Right Panel: Relation between ξ HII ion and Hα EW for the z ∼ 2 galaxies in our sample (blue circles). A linear relation is fitted between ξ HII ion and log EW Hα (black dashed line, with grey area shows the 95 percent credible region).

  [O III]/[O II] (O32) and [Ne III]/[O II] (Ne3O2), two ionization-sensitive emission line ratios, on the [O III] and Hα EW.

  [O III] emission (EW[OIII]λ5007 > 225 Å). In the top panel of Figure10, we show the dependence of O32 on the [O III]λ5007 EW for individual galaxies (blue circles) and for our composite spectra (red squares). The results clearly demonstrate that O32 increases with the [O III]λ5007 EW, indicating more extreme ionizing conditions in the high EW galaxies. The composite spectra reveal [O III]/[O II] ratios of O32 = 1.4 ± 0.2, 3.6 ± 0.2, 4.8 ± 0.3, and 9.1 ± 0.5 for stacks including galaxies with EW [OIII]λ5007 = 0 -225 Å, 225 -450 Å, 450 -800 Å, and 800 -2500 Å, respectively.

  equivalent widths comparable to the RB16 z 7 Lyα emitters (EW [OIII]λ5007 = 800 -2500 Å), we find evidence for even more extreme ISM conditions, with 8 of the 11 sources with robust [O III]/[O II] measurements having O32 > 6. Not surprisingly given these trends, the two highest [O III]/[O II] ratios observed in our survey (O32 = 17 ± 4 and 22 ± 3) are found in two of the most extreme line emitters (EW [OIII]λ5007 = 1216±121 Å and 1893 ± 59 Å, respectively). The highest O32 values appear to become commonplace in the extremely young systems ( 10 Myr for constant star formation, §3) that power intense optical line emission.

Figure 9 .

 9 Figure 9. BEAGLE model-derived ξ ion as a function of Hβ EW (left panel) and [O III]λ5007 EW (right panel) at z = 1.3 -2.4 (blue circles; this work) and z ∼ 0 (green squares; Chevallard et al. 2018). The ξ ion is computed based on the intrinsic UV stellar luminosity (assuming a SMC extinction curve in order to be consistent with Chevallard et al. 2018), and only objects with similar EWs in the two samples (Hβ EW = 80 -400 Å in the left panel; [O III]λ5007 EW = 400 -2500 Å in the right panel) are shown.

5. 1

 1 Implications for the ISM Conditions and the Escape of Ionizing Radiation in Reionization-Era Galaxies In §4.3, we demonstrated that [O III]/[O II] ratio increases with [O III] EW, reaching very large value (O32 > 10) for sources with EW [OIII]λ5007 > 1000 Å. The physical factors regulating O32 have been discussed in detail elsewhere (e.g. Kewley et al. 2013; Nakajima & Ouchi 2014; Sanders et al. 2016). These studies show that the [O III]/[O II] ratio depends strongly on the ionization parameter of the gas with a weaker dependence on its metallicity. Assuming ionization balance, the ionization parameter in turn depends on the ionization photon production rate, the gas density, and the volume filling factor of ionized gas. The hardness of the ionizing spectrum also impacts the [O III]/[O II] ratio, increasing O32 at fixed ionization parameter (Sanders et al. 2016). While a comprehensive investigation of the trend seen in Figure 10 is beyond the scope of this paper, the boosted O32 values likely reflect the combined effect of a large ionization parameter and hard ionizing spectrum, the latter expected in the very young stellar populations present in the largest [O III] EW galaxies. Extreme optical line emitters are not the only sources at high redshift known to have large O32 ratios. Spectroscopic investigations have shown large O32 values in high redshift galaxies with large EW Lyα emission. Nakajima et al. (2016) measure [O III]/[O II] ratios in the range O32 = 6.0 -11.5 for four Lyα emitting galaxies (EW Lyα = 26 -61 Å) at z 3.1, comparable to the galaxies with EW [OIII]λ5007 > 800 Å in our sample. Similarly Erb et al. (2016) report [O III]/[O II] ratios of 2.4 to 12.2 (median O32 = 4.9) for 11 Lyα emitters at z 2 with EW Lyα > 20 Å. The large EW Lyα emitters with low masses likely probe a very similar population. Strong [O III] emitters at z 1 have previously been shown to have [O III]/[O II] ratios (O32 3 -100) similar to those we report for our sample (Kakazu et al. 2007). Among local galaxy samples, the Green Peas also exhibit large [O III]/[O II] ratios and large [O III] equivalent widths (Cardamone et al. 2009). In Figure 10, we overlay the O32 and Ne3O2 values of the Green Peas on our high redshift relationships (green open circles). The scaling between the ionization-sensitive ratios and the [O III] EW appears very similar to what we have found at high redshift.

  Figure 10. Relation between [O III]λ5007 EW and ionizationsensitive line ratios ([O III]/[O II], upper panel; [Ne III]/[O II], lower panel) of the z = 1.3 -2.4 galaxies in our sample (blue circles). Line ratios of objects without robust [O II] (cyan triangles in the upper panel) or [Ne III] detections (cyan upside-down triangles in the lower panel) are shown as 2σ lower or upper limits. The average line ratios derived from composite spectra of galaxies with different EWs are shown by red squares. The average EW of z 7 -8 galaxies (L13) is highlighted by yellow dashed line, and the orange area shows EWs of z > 7 galaxies in the RB16 sample. For comparison, [O III]λ5007 EWs and line ratios of Green Peas (Cardamone et al. 2009; Hawley 2012) are shown by green open circles.

  [O III] line shifts out of the NIRSpec bandpass above z 9, leaving only faint nebular lines in the rest-frame optical ([O II], [Ne III])) and rest-UV (C III], O III], C IV).

  Figure 11. Left Panel: [O II]λλ3727, 3729 EW as a function of [O III]λ5007 EW for the z = 1.3 -2.4 galaxies in our sample. Sources with > 2σ detections of [O II] are shown by blue circles, while cyan upside-down triangles show the 2σ upper limits. Red squares present the median [O II] EWs of galaxies in different [O III]λ5007 EW bins. Right Panel: Hβ EW as a function of [O III]λ5007 EW for the z = 1.3 -2.4 galaxies in our sample. Sources with > 2σ detections of Hβ are shown by blue circles, while cyan upsidedown triangles show the 2σ upper limits. Red squares present the median Hβ EWs of galaxies in different [O III]λ5007 EW bins. For comparison, EWs of Green Peas (Cardamone et al. 2009; Hawley 2012) are shown by green open circles.

  , we plot the dependence of the Hβ EW on the [O III]λ5007 EW in our spectroscopic sample. The Hβ EW scales with the [O III] EW, with values of EW Hβ = 53, 98, and 198 Å for [O III]λ5007 EW bins of EW [OIII]λ5007 = 225-450, 450-800, and 800 -2500 Å, respectively. As can be seen in Figure 11, the relationship between the Hβ EW and the [O III] EW is very similar for the Green Peas (at EW [OIII]λ5007 > 350 Å), suggesting very little evolution in the trend over 0.3 z 2.4. If the evolution remains slow at z 3, it indicates that Hβ should be detectable in typical (H 27) reionization-era galaxies with deep ( 10 hr) NIRSpec integrations.

  Examples of the broadband SEDs of the extreme[O III] line emitting galaxies at z = 1.3 -2.4 in our sample. The upper panels show three of the sources with very large EWs in our sample ([O III]λ5007 EW > 1000 Å) that similar to the z > 7 galaxies in RB16. The lower panels show three of the sources with equivalent widths slightly above the average EW of z > 7 galaxies ([O III]λ5007 EW 450 Å; L13). Observed broadband photometry is shown as solid black circles. The best-fitting SED models (inferred from BEAGLE, see §3)

		22							24						24
		23
									25						25
	mag (AB)	24
		25		[OIII]	Hα		z=1.5190	26		[OIII]	Hα	z=1.6549	26	[OIII]	Hα	z=1.6772
				[OIII]/[OII]=21.6			UDS-29267			[OIII]/[OII]=17.1		UDS-13027		[OIII]/[OII]=9.2	GOODS-N-36886
				log ξ HII ion =25.80			EW[OIII]λ5007=1893	Å			log ξ HII ion =25.62		EW[OIII]λ5007=1216	Å		log ξ HII ion =25.59	EW[OIII]λ5007=1072	Å
		26	0	1		2	3		27	0	1	2	3		27	0	1	2	3
		22							24						23
		23							25
	mag (AB)														24
		24			[OIII]	Hα	z=2.2953	26		[OIII]	Hα	z=1.6847		[OIII]	Hα	z=1.4612
				[OIII]/[OII]=4.1			AEGIS-22931			[OIII]/[OII]=5.0		GOODS-N-29675		[OIII]/[OII]=4.8	UDS-14019
				log ξ HII ion =25.34			EW[OIII]λ5007=645	Å			log ξ HII ion =25.36		EW[OIII]λ5007=561	Å		log ξ HII ion =25.21	EW[OIII]λ5007=490	Å
		25	0	1		2	3		27	0	1	2	3		25	0	1	2	3
				λobs (µm)							λobs (µm)				λobs (µm)
	Figure 1. are plotted by solid blue lines, and synthetic photometry is shown as open green squares. Strong rest-frame optical emission lines, [O
	III]λ5007 and Hα, are highlighted by dashed black lines. Spectroscopic redshift, object ID, and rest-frame [O III]λ5007 EW are listed in
	each panel. The LyC photon production efficiency (ξ HII ion ) and the ionization-sensitive line ratio [O III]/[O II] of each object are also listed.

Table 1 .

 1 Summary

	Mask Name Number of Target	R.A.		Decl.	P.A.		Grism	Filter	Exposure Time Average Seeing
		(hh:mm:ss)	(dd:mm:ss)	(deg)				(seconds)	( )
	(1)	(2)	(3)		(4)	(5)		(6)	(7)	(8)	(9)
	udsel1	2:17:49.000	-5:13:02.00	-11.00	J		zJ	14400	1.2
	udsel2	2:17:26.000	-5:13:20.00 -116.00	J		zJ	14400	0.8
	udsel3	2:17:26.000	-5:12:13.00	88.00		J		zJ	28800	0.8
	udsel3	2:17:26.000	-5:12:13.00	88.00		H3000	H	7200	0.8
	egsel3	14:19:25.973	52:50:08.43	177.00		J		zJ	21000	0.6
	egsel3	14:19:25.973	52:50:08.43	177.00		H3000	H	6000	1.0
	egsel4	14:19:58.000	52:55:04.01	179.00		J		zJ	10200	0.8
	egsel4	14:19:58.000	52:55:04.01	179.00		H3000	H	2400	1.0
	gdnel3	12:37:14.173	62:19:29.10	-162.00	J		zJ	14400	0.8
	gdnel3	12:37:14.173	62:19:29.10	-162.00	H3000	H	12000	1.2
	gdnel4	12:36:31.204	62:14:51.55	-125.00	J		zJ	12000	0.8
	gdnel4	12:36:31.204	62:14:51.55	-125.00	H3000	H	7200	0.8
	udsel5	2:17:17.088	-5:14:53.13	-86.00	J		zJ	7200	1.2
	udsel5	2:17:17.088	-5:14:53.13	-86.00	H3000	H	14400	1.2
	udsel6	2:17:49.000	-5:13:02.00	106.00		J		zJ	7200	1.2
	udsel6	2:17:49.000	-5:13:02.00	106.00		H3000	H	14400	1.6
	egsel1	14:19:37.358	52:50:22.97	-179.00	J		zJ	14400	1.1
	egsel1	14:19:37.358	52:50:22.97	-179.00 K3000	Kspec	3600	0.7
	egsel2	14:20:07.295	52:54:13.13	-134.00	J		zJ	3600	0.6
	egsel2	14:20:07.295	52:54:13.13	-134.00 K3000	Kspec	3600	1.6
	gdnel5	12:37:08.882	62:19:36.39	-162.00	H3000	H	3600	0.8
	gdnel5	12:37:08.882	62:19:36.39	-162.00 K3000	Kspec	3600	0.7
	gdnel6	12:37:22.384	62:18:02.04	-143.00	J		zJ	7200	1.0
	egsel5	14:20:11.119	52:56:50.92	120.00		J		zJ	7200	0.8
	udsel7	2:17:15.656	-5:14:04.59	93.00		J		zJ	10800	0.8
	udsel7	2:17:15.656	-5:14:04.59	93.00		K3000	Kspec	3600	0.8
	udsel8	2:17:36.777	-5:11:40.87	-77.00	J		zJ	10800	0.8
	udsel8	2:17:36.777	-5:11:40.87	-77.00	K3000	Kspec	3600	0.8
	egsel6	14:19:15.000	52:46:29.00	-146.00	J		zJ	13200	1.0
	egsel6	14:19:15.000	52:46:29.00	-146.00	H3000	H	7200	1.0
	egsel6	14:19:15.000	52:46:29.00	-146.00 K3000	Kspec	3600	1.0
	Mask Name Number of Targets	R.A.		Decl.	P.A.	Band Exposure Time Average Seeing
			(hh:mm:ss)	(dd:mm:ss)	(deg)		(seconds)	( )
	(1)	(2)	(3)		(4)			(5)	(6)	(7)	(8)
	cos y4	12	10:00:23.43	2:20:27.06	167.0	Y	8640	0.74
	gs y5	16	3:32:14.63	-27:44:14.42 158.0	Y	11880	0.94
	uds yjh2	27	2:17:07.02	-5:10:29.00	49.0	Y	2160	0.52
	uds yjh2	27	2:17:07.02	-5:10:29.00	49.0	J	1680	0.54
	uds yjh2	27	2:17:07.02	-5:10:29.00	49.0	H	1920	0.46
	egs z8.6	20	14:20:02.22	52:54:49.08	292.0	H	15360	0.72
	egs z7.5	19	14:20:28.82	53:00:55.53	107.0	J	7080	0.90

of MMT/MMIRS observations. Column (1): mask name; Column (2): number of science targets on each mask, alignment stars and slit stars are not included; Column (3): right ascension of the mask center; Column (4): declination of the mask center; Column (5): position angle of the mask; Column (6): grism of the mask observed; Column (7): filter of the mask observed; Column (8): Total exposure time of the mask in each grism + filter set; Column (9): average seeing during the observation, in full width at half maximum (FWHM).

Table 2 .

 2 Summary

of Keck/MOSFIRE observations. Column (1): mask name; Column (2): number of science targets on each mask, alignment stars and slit stars are not included; Column (3): right ascension of the mask center; Column (4): declination of the mask center; Column (5): position angle of the mask; Column (6): band of the mask observed; Column (7): Total exposure time of the mask in each grating; Column (8): average seeing during the observation, in full width at half maximum.

Table 3 .

 3 Summary

	Target Name	R.A.	Decl.	P.A.	z spec	Configuration Exposure Time Average Seeing
		(hh:mm:ss)	(dd:mm:ss)	(deg)		line mm -1	(seconds)	( )
	UDS-29267	02:17:25.322 -05:10:40.40 -168.99 1.5190	1200	7200	1.1
	UDS-14019	02:17:28.554 -05:13:44.86 +114.81 1.4616	1200	6000	1.3

of MMT/Red Channel observation. Column (1): target name; Column (2): right ascension of the target; Column (3): declination of the target; Column (4): position angle of the slit; Column (5): spectroscopic redshift of the target; Column (6): configuration used; Column (7): Total exposure time on the target; Column (8): average seeing during the observation, in full width at half maximum.

  ([O II]λλ3727, 3729, Hβ, [O III]λλ4959, 5007, and Hα). The version of BEAGLE used in this work adopts the recent photoionization models of star-forming galaxies of GutkinFigure 4. MMT and Keck composite spectra of extreme [O III] emitters at z = 1.3 -2.4, grouped by sources with different [O III]λ5007 EWs. The spectra are ranked by [O III]λ5007 EW ranges, from the most extreme line emitters ([O III]λ5007 EW = 800 -2500 Å) on the top to the objects with similar EWs as the typical galaxies at z ∼ 2 ([O III]λ5007 EW = 0 -225 Å) in the bottom. The 1σ error spectra are also shown as grey regions. The locations of strong rest-frame optical emission lines are highlighted by black dotted lines. Zoom-in spectra of relatively faint [O II] and [Ne III] emission lines are shown in the left panels. The composite spectra of galaxies with higher [O III]λ5007 EWs show higher [O III]/[O II] and [Ne III]/[O II] ratios than those of galaxies with lower EWs.

		5
	λ (arbitrary unit)	1 2 3 4	[OIII]λ5007 EW = 800 -2500 Å	[OII]	3700 3730 0.0 0.2 0.4 [OII]	3850 3880 [NeIII]	[NeIII]	Hβ	[OIII]	[OIII]	[NII]	Hα	[NII]
	f
		0
		-1	3650	3700	3750	3800	3850	4850	4900	4950	5000	5050	6550	6600
		5									rest-frame wavelength ( Å)
	λ (arbitrary unit)	1 2 3 4	[OIII]λ5007 EW = 450 -800 Å	[OII]	3700 3730 0.0 0.2 0.4 0.6 0.8 [OII]	3850 3880 [NeIII]	[NeIII]	Hβ	[OIII]	[OIII]	[NII]	Hα	[NII]
	f
		0
		-1	3650	3700	3750	3800	3850	4850	4900	4950	5000	5050	6550	6600
		5				1.2					rest-frame wavelength ( Å)
	λ (arbitrary unit)	1 2 3 4	[OIII]λ5007 EW = 225 -450 Å	[OII]	3700 3730 0.0 0.4 0.8 [OII]	3850 3880 [NeIII]	[NeIII]	Hβ	[OIII]	[OIII]	[NII]	Hα	[NII]
	f
		0
		-1	3650	3700	3750	3800	3850	4850	4900	4950	5000	5050	6550	6600
		5									rest-frame wavelength ( Å)
	λ (arbitrary unit)	1 2 3 4	[OIII]λ5007 EW = 0 -225 Å	[OII]	3700 3730 0.0 0.4 0.8 1.2 1.6 [OII]	3850 3880 [NeIII]	[NeIII]	Hβ	[OIII]	[OIII]	[NII]	Hα	[NII]
	f
		0
		-1	3650	3700	3750	3800	3850	4850	4900	4950	5000	5050	6550	6600
											rest-frame wavelength ( Å)

et al.

(2016)

, which describes the emission from stars and interstellar gas based on the combination of the latest version of

Bruzual & Charlot (2003)

stellar population synthesis model with the photoionization code CLOUDY

(Ferland et al. 2013)

. The main adjustable parameters of the photoionized gas are the interstellar metallicity, Z ISM , the typical ionization parameter of a newly ionized H II region, U S (which characterizes the ratio of ionizing-photon to gas densities at the edge of the Str omgren sphere), and the dustto-metal (mass) ratio, ξ d (which characterizes the depletion of metals on to dust grains). We consider models with C/O abundance ratio equal to the standard value in nearby galaxies (C/O) ≈ 0.44. The chosen C/O ratio does not significantly impact the results presented in this paper. We will

Table 4 .

 4 Statistics for sub-samples of extreme [O III] emitters ([O III]λ5007 EW > 225 Å) used for line ratio analysis in §4. Spectral coverage requirements and the resulting sub-sample sizes are listed. For each line ratio analysis, individual objects with robust line ratio (Balmer decrement Hα/Hβ, [O III]/[O II], or [Ne III]/[O II]) measurements have more strict requirements than objects used to make composite spectra. For example, the Balmer decrement requires spectra that span between Hβ and Hα (73 sources), while measuring robust Balmer decrement in an individual object also requires detecting Hβ and Hα with 3σ significance or higher (64 sources).

	Notes.
	1 Sub-sample used for Balmer decrement analysis in §4.1.

2 Sub-sample used for [O III]/[O II] ratio analysis in §4.3. 3 Sub-sample used for [Ne III]/[O II] ratio analysis in §4.3.

Table 5 .

 5 O III] and Hα EW are directly linked to the luminosity-weighted age of Summary

	Object ID	R.A. (hh:mm:ss)	Decl. (dd:mm:ss)	z spec	m F606W (AB)	EW [OIII]λ5007 ( Å)	EW Hα ( Å)
	AEGIS-29446	14:19:31.989 +52:53:22.789 1.5506 26.13 ± 0.10	1059 ± 152	662 ± 111
	AEGIS-26095	14:19:33.549 +52:52:52.797 1.5625 26.19 ± 0.11	340 ± 47	164 ± 26
	AEGIS-22931	14:19:33.787 +52:52:07.560 2.2953 24.05 ± 0.02	645 ± 23	790 ± 51
	AEGIS-23885	14:19:31.483 +52:51:57.961 1.5721 24.64 ± 0.03	608 ± 33	578 ± 38
	AEGIS-21396	14:19:32.014 +52:51:27.694 2.3206 25.35 ± 0.05	241 ± 13	458 ± 173
	AEGIS-31255	14:19:17.457 +52:51:09.896 1.5970 24.42 ± 0.03	532 ± 34	463 ± 36
	AEGIS-19479	14:19:31.179 +52:50:52.304 1.5969 24.79 ± 0.03	296 ± 17	366 ± 27
	AEGIS-20493	14:19:24.174 +52:49:48.020 1.6725 23.73 ± 0.01	376 ± 7	498 ± 13
	AEGIS-11745	14:19:33.812 +52:49:28.451 1.6069 26.15 ± 0.10	620 ± 95	728 ± 165
	AEGIS-11452	14:19:26.785 +52:48:04.364 1.6720 23.33 ± 0.02	464 ± 17	324 ± 14
	AEGIS-03127	14:20:07.698 +52:53:11.858 2.3089 24.51 ± 0.03	671 ± 56	686 ± 90
	AEGIS-14156	14:19:50.376 +52:52:59.100 1.6738 24.86 ± 0.03	249 ± 18	249 ± 38
	AEGIS-00361	14:20:21.976 +52:54:46.588 1.6678 23.78 ± 0.01	814 ± 53	1028 ± 57
	AEGIS-15032	14:19:26.935 +52:49:02.770 1.6132 26.21 ± 0.10	1121 ± 210	716 ± 170
	AEGIS-04337	14:19:35.878 +52:47:54.834 1.3985 26.80 ± 0.16	2296 ± 356	1876 ± 285
	AEGIS-04711	14:19:34.958 +52:47:50.219 2.1839 24.06 ± 0.01	1060 ± 25	671 ± 33
	AEGIS-16513	14:19:19.204 +52:48:00.052 1.4633 24.78 ± 0.03	361 ± 24	351 ± 28
	AEGIS-24361	14:19:08.525 +52:47:59.956 1.4714 25.86 ± 0.07	1471 ± 145	1192 ± 156
	AEGIS-08869	14:19:20.544 +52:46:21.299 1.5792 25.59 ± 0.06	909 ± 95	1017 ± 88
	AEGIS-15778	14:19:11.210 +52:46:23.414 2.1716 25.24 ± 0.04	1001 ± 42	891 ± 47
	AEGIS-01387	14:19:22.694 +52:44:43.040 2.1700 24.77 ± 0.03	426 ± 32	461 ± 37
	AEGIS-06264	14:19:17.443 +52:45:05.370 2.1879 24.72 ± 0.03	651 ± 46	456 ± 46
	AEGIS-07028	14:19:16.128 +52:45:03.763 2.2925 24.91 ± 0.04	1023 ± 67	783 ± 62
	AEGIS-20217	14:19:02.413 +52:45:56.017 1.6306 25.24 ± 0.05	341 ± 27	601 ± 55
	AEGIS-17118	14:19:03.204 +52:45:16.919 1.6055 24.02 ± 0.02	332 ± 19	376 ± 21
	GOODS-N-38085 12:37:18.165 +62:22:29.258 1.5214 24.80 ± 0.02	930 ± 79	913 ± 80
	GOODS-N-37906 12:37:16.752 +62:22:04.152 1.6847 24.65 ± 0.02	316 ± 23	379 ± 30
	GOODS-N-37876 12:37:12.919 +62:21:59.569 1.6148 25.53 ± 0.04	283 ± 23	357 ± 32
	GOODS-N-37878 12:37:06.160 +62:22:00.005 1.5508 25.86 ± 0.05	456 ± 46	305 ± 32
	GOODS-N-37296 12:37:24.594 +62:21:00.900 1.5136 24.47 ± 0.01	685 ± 38	741 ± 44
	GOODS-N-36886 12:37:31.259 +62:20:36.596 1.6772 25.34 ± 0.03	1072 ± 144	1114 ± 151
	GOODS-N-36583 12:37:33.171 +62:20:23.489 1.5495 25.27 ± 0.03	480 ± 50	699 ± 70
	GOODS-N-36852 12:37:23.247 +62:20:34.465 1.5970 25.94 ± 0.06	1678 ± 380	1209 ± 274
	GOODS-N-36684 12:37:22.580 +62:20:26.984 1.5965 25.44 ± 0.04	442 ± 44	392 ± 40
	GOODS-N-36273 12:37:01.427 +62:20:10.543 2.2422 26.05 ± 0.06	405 ± 40	541 ± 119
	GOODS-N-29675 12:37:07.081 +62:17:18.971 1.6847 25.29 ± 0.03	561 ± 43	938 ± 70
	GOODS-N-35204 12:37:06.813 +62:19:35.508 2.0976 23.89 ± 0.01	322 ± 12	403 ± 29
	UDS-06377	02:17:42.856 -05:15:19.134 1.6642 25.31 ± 0.05	943 ± 94	1169 ± 135
	UDS-12539	02:17:53.733 -05:14:03.196 1.6211 24.36 ± 0.02	882 ± 33	803 ± 31
	UDS-24003	02:17:54.729 -05:11:44.020 1.6205 24.55 ± 0.03	259 ± 11	266 ± 17
	UDS-07447	02:17:18.162 -05:15:06.275 1.5972 24.54 ± 0.02	666 ± 36	529 ± 26
	UDS-13027	02:17:16.355 -05:13:56.240 1.6549 25.28 ± 0.04	1216 ± 121	878 ± 145
	UDS-21873	02:17:17.096 -05:12:09.652 1.6552 24.59 ± 0.02	271 ± 19	387 ± 35
	UDS-15658	02:17:24.262 -05:13:25.612 1.6869 24.76 ± 0.04	307 ± 23	334 ± 34
	UDS-29267	02:17:25.322 -05:10:40.397 1.5190 23.77 ± 0.01	1893 ± 59	1691 ± 73
	UDS-26182	02:17:28.134 -05:11:17.272 1.4867 25.45 ± 0.04	1120 ± 98	1175 ± 137
	UDS-28931	02:17:31.478 -05:10:45.926 1.6000 23.33 ± 0.01	244 ± 5	308 ± 10
	UDS-30015	02:17:36.517 -05:10:31.256 1.6649 24.62 ± 0.02	917 ± 55	607 ± 48
	UDS-19818	02:17:02.612 -05:12:34.628 1.6241 25.16 ± 0.04	418 ± 35	470 ± 47
	UDS-29624	02:17:00.684 -05:10:34.903 1.6632 26.26 ± 0.11	881 ± 195	779 ± 208
	UDS-36954	02:17:14.900 -05:09:06.174 1.6578 25.41 ± 0.04	238 ± 20	347 ± 68
	UDS-14019	02:17:28.554 -05:13:44.857 1.4612 23.80 ± 0.01	490 ± 15	614 ± 17
	UDS-05122	02:17:22.253 -05:15:33.613 1.4628 24.67 ± 0.03	307 ± 16	291 ± 14
	UDS-05342	02:17:09.216 -05:15:31.990 1.5981 23.58 ± 0.01	437 ± 12	363 ± 8
	UDS-11693	02:17:03.893 -05:14:13.664 2.1854 24.52 ± 0.02	596 ± 33	592 ± 73
	UDS-15128	02:17:38.209 -05:13:32.092 2.3113 24.25 ± 0.02	371 ± 16	276 ± 97
	UDS-12435	02:17:38.609 -05:14:05.366 1.6112 24.06 ± 0.01	247 ± 9	173 ± 6
	UDS-11387	02:17:39.228 -05:14:17.689 1.4048 24.64 ± 0.02	602 ± 40	469 ± 34
	UDS-17713	02:17:42.450 -05:13:02.615 2.2484 24.09 ± 0.02	328 ± 16	249 ± 50

of the 64 extreme [O III] emitters (EW [OIII]λ5007 > 225 Å) with robust Balmer line detections (S/N > 3) at z = 1.3-2.4. Column (1): galaxy ID in 3D-HST v4 catalogs; Column (2): right ascension of the galaxy; Column (3): declination of the galaxy; Column (4): redshift measured from MMT or Keck spectra; Column (5): HST F606W magnitude (in AB magnitude); Column (7): rest-frame [O III]λ5007 equivalent width; Column (8): rest-frame Hα equivalent width.

Table 5 .

 5 Continued.

	5
	age /yr)
	log(t

Figure 6. Physical parameters (stellar mass, star formation rate, specific star formation rate, and stellar age) as functions of [O III]λ5007

EW for the 77 emission line galaxies with robust Hα and Hβ measurements (S/N > 3) at z = 1.3 -2.4 in our sample (64 extreme

[O III]

emitters with EW [OIII]λ5007 > 225 Å and 13 additional galaxies with lower EWs). All the parameters are derived from stellar population synthesis modeling using BEAGLE. The upper left, the upper right, the lower left, and the lower right panel shows the stellar mass, the star formation rate, the specific star formation rate, and the stellar age versus [O III]λ5007 EW, respectively.

  Galactic extinction curve and an intrinsic Balmer decrement,

		0.7
	Normalized Counts	0.2 0.3 0.4 0.5 0.6				[OIII]λ5007 EW Å = 800 -2500 [OIII]λ5007 EW Å = 450 -800 z = 2.1 -2.6; Steidel+2016 z = 1.3 -2.6; Reddy+2015
		0.1
		0	0.0	0.1	0.2	0.3	0.4
					E(B -V ) gas
	Figure 7. Probability distribution of nebular dust reddening for
	extreme [O

III] emitters with EW [OIII]λ5007 > 450 Å at z = 1.3-2.4.

The extinction E(B -V ) gas is inferred from Balmer decrement, assuming the

Cardelli et al. (1989)

extinction curve and intrinsic Hα/Hβ intensity ratio = 2.86. The extinction values inferred from composite spectra of extreme

[O III]

emitters with EW [OIII]λ5007 = 800-2500 Å and EW [OIII]λ5007 = 450-800 Å are plotted by red solid and dashed lines. For comparison, we also show the extinction values inferred from more massive, older systems at z ∼ 2 from KBSS

(Steidel et al. 2016)

as the orange dotted line and MOSDEF survey

(Reddy et al. 2015)

as the orange dashed-dotted line.

  ] EW.Our current sample contains 64 large equivalent width[O III] emitters (EW [OIII]λ5007 > 225 Å) with the spectroscopic measurements of Hα and Hβ necessary to infer ξ HII ion . In Figure 8, we show the implied ξ HII ion values as a function of the [O III]λ5007 and Hα equivalent widths. It is immediately clear that ξ HII ion scales with both optical lines at EW [OIII]λ5007 > 225 Å. We first consider galaxies with EW [OIII]λ5007 450 Å. For the [O III]/Hβ line ratios exhibited by extreme emission line galaxies (see §2.1), this [O III]λ5007 EW is comparable to that implied by IRAC flux excesses in composite SEDs of z 7 -8 galaxies (L13).

MNRAS 000, 1-24 (2019)

Adopting the theoretical flux ratio I (5007)/I (4959) = 3, this selection is equivalently stated as [O III]λ5007 EW

> 225 Å.

https://bitbucket.org/chil_sai/mmirs-pipeline

http://www2.keck.hawaii.edu/inst/mosfire/magma.html

https://www2.keck.hawaii.edu/inst/mosfire/drp.html MNRAS 000, 1-24 (2019)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000, 1-24 (2019)

ACKNOWLEDGEMENTS

We are grateful for enlightening conversations with Michael Maseda, Alice Shapley, Irene Shivaei, and David Sobral. DPS acknowledges support from the National Science Foundation through the grant AST-1410155. Observations reported here were obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution. We thank the MMT/MMIRS queue observers, ShiAnne Kattner and Chun Ly, who helped a lot collect the data. This work was partially supported by a NASA Keck PI Data Award, administered by the NASA Exoplanet Science Institute. Data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to the National Aeronautics and Space Administration through the agency's scientific partnership with the California Institute of Technology and the University of California. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors acknowledge the very significant cultural role that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.