Broadband Selection, Spectroscopic Identification, and Physical Properties of a Population of Extreme Emission-line Galaxies at 3 < z < 3.7^*

We carried out a spectroscopic follow-up observation for a subset of the Ks-band excess EELG candidates. Spectroscopic targets were selected primarily to maximize the number of observed EELG candidates with 3 < z_phot < 3.7. We also supplementarily considered criteria rJL ≡ (J − [3.6]) − 1.4(r − J) > 0, which was proposed to select objects at 2.5 ≲ z ≲ 4 (Daddi et al. 2004), and iJHK colors to remove low-redshift (z ≲ 0.2) ELGs (Tadaki et al. 2013).

We used the Multi-Object InfraRed Camera and Spectrograph (MOIRCS; Ichikawa et al. 2006; Suzuki et al. 2008) at the Subaru Telescope. In 2015, the two detectors in MOIRCS were upgraded from Hawaii-2 to Hawaii-2RG (Fabricius et al. 2016; Walawender et al. 2016). We used the HK500 grism covering 1.3–2.3 μm with a 08 slit, which provides the spectral resolution of R ≃ 350. After the detector upgrade, the lower readout noise allows us to shorten the exposure time to reach the background-limited performance, enabling us to better capture the temporal variation of the sky background. Although the throughput in the HK band improved only by ≲5%, the lower readout noise results in a deeper limiting flux than the previous detectors by cleaner sky subtraction, which is crucial for the relatively low resolution of the HK500 grism used in this study.

The observation was carried out during the first half nights on 2017 April 9 and 17 under clear sky conditions with 05–08 seeing. Three masks for three pointings (each has a diameter of ≃6') were used to observe in total 23 EELG candidates. We used a modified AB dithering pattern with three dithering widths of 28, 30, and 32, to increase the signal-to-noise ratio (S/N) and to avoid bad pixels (Kriek et al. 2015), and 180 s on-source exposure per dithering position. Each mask has been integrated ≃40–110 minutes. The detailed observing log is shown in Table 1. We observed A0V-type stars at the beginning of the nights for flux calibration and telluric correction using the identical instrument setup to the science exposures.

Data reduction was carried out using the MCSMDP¹⁰ pipeline (Yoshikawa et al. 2010) and custom scripts. Since the original MCSMDP was developed for the previous detectors, we made modifications to properly handle the new Hawaii-2RG detectors. The data were flat-fielded by using dome-flat frames that were taken under the identical instrument setup to the science frames. Bad pixels and cosmic rays were identified by using bad-pixel maps created for the new detectors and by using the pair of images in the dithering, respectively, and linearly interpolated using adjacent pixels along the spatial direction. Background sky was subtracted by taking a difference between A and B positions. The optical distortion was then corrected by applying a polynomial function determined by the MOIRCS imaging reduction pipeline MCSRED.¹¹

Individual two-dimensional (2D) spectra were extracted from each frame, and the wavelength calibration was carried out by using OH-airglow lines (Rousselot et al. 2000). The typical uncertainty of the wavelength calibration is ≃3 Å, which is about 40% of the pixel scale. Based on the wavelength solution, the extracted 2D spectra were rectified so that sky lines are aligned along the spatial direction and the secondary sky subtraction procedure was made in order to remove residual background by fitting a linear function at each wavelength pixel.

Standard-star frames of A0V-type stars were processed in the same manner as the science object frames. The 1D spectra of the standard stars were extracted using the apall task in the Image Reduction and Analysis Facility (IRAF; Tody 1986, 1993), and the system total response curve was derived by comparing the observed spectra with a theoretical template (Kurucz 1979). Along with the telluric correction, we have carried out absolute flux calibration for the standard spectra by scaling them to Two Micron All Sky Survey magnitudes correcting for the slit loss. Flux-calibrated 2D frames of each object are then median stacked after 3σ clipping, with offsets corresponding to the dithering widths. One-dimensional (1D) noise spectra were calculated using the 1σ dispersion of the counts along the spatial direction at each wavelength pixel in the nondetected slits for each mask after the stacking procedure.

1D science and corresponding noise spectra were extracted by adopting the optimal extraction algorithm (Horne 1986). Here we assume a Gaussian spatial profile at the strongest detected emission line for weighting and a straight trace for extraction. Before the extraction, we confirmed that the trace is straight along the dispersion direction by using bright stars and galaxies put in the science masks as filler targets.

Figures 3 and 4, respectively, show 2D and 1D spectra of emission-line-detected EELG candidates. Twenty-one out of 23 observed EELG candidates show detected emission lines, and the origins of K-band excesses are confirmed as [O iii] and Hα at z ≃ 3.3 and z ≃ 2.2 for 19 and 2 objects, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As seen in Figures 3 and 4, 21 objects show significant detection of multiple emission lines, which makes the emission-line identification and rough redshift estimate straightforward. Then, we measured the emission-line properties in two fitting processes as explained below.

First, we only fit the primary emission line. The primary emission line here is the strongest emission line, i.e., [O iii]λ5007 and Hα for z ≃ 3.3 and z ≃ 2.2, respectively. We assume a Gaussian profile with a flat continuum. Therefore, there are four free parameters, namely, the redshift, line width σ, total line flux, and constant continuum flux. We used ±500 Å from the line center for the fitting procedure. When the primary line is Hα, adjacent [N ii]λλ6548,6583 lines are fit simultaneously, as they are not fully resolved. In this case, all three lines are assumed to have the same redshift, σ, and continuum flux. Also, [N ii]λ6548 flux is assumed to be one-third of [N ii]λ6583 flux.

For the fitting, we used emcee, a Python implementation of an affine-invariant Markov Chain Monte Carlo (MCMC) algorithm (Goodman & Weare 2010; Foreman-Mackey et al. 2013) through lmfit (Newville et al. 2019). Priors are assumed to be top hat in all free parameters with only minimal boundaries, namely, positivity of the σ and line flux. We used 100 walkers and 1000 steps per walker, of which 200 steps were used as a burn-in process and discarded from the posterior sampling. Since the S/Ns are typically high (≫10) for the primary emission lines, the convergence is quickly achieved and all parameters for line profiles are well constrained.

In the second fitting process, we fit the rest of the strong emission lines, namely, [O ii]λλ3726,3729, [Ne iii]λ3869, Hβ, and [S ii]λλ6717,6731, within the observed wavelength range. For z ≃ 2.2 objects, [O iii]λλ4959,5007 were also fit in this process. Here the redshift and line width σ were fixed to those determined in the first fitting process, but allowing them to be different within ±0.5σ and ±0.1σ from those of the primary line, respectively. We also used ±500 Å from the line centers for the fitting, but continuum was assumed to be a linear function across the spectrum. Because nonprimary lines are fainter than the primary and include nondetections, we used 100 walkers with 10,000 steps per walker and 2000 burn-in steps to sample the posterior distributions. These numbers are large enough to distinguish between converged and nonconverged parameters.

For each parameter, the best-fit parameter and the corresponding 1σ error are defined as the median and half the interval between 16th and 84th percentiles, respectively, of the posterior distribution. To judge whether the line is detected or not, we compared the best-fit line flux with the 3σ noise computed from 1D noise spectra integrated with a Gaussian weight corresponding to the primary's line profile. While continuum fluxes are not well constrained in general for our sample, as they are barely detected as seen in Figures 3 and 4, this does not affect the emission-line flux measurements. Table 2 shows the measured emission-line properties. Here we only list emission-line fluxes used in the following discussion, though we also fit [Ne iii]λ3869 and [S ii]λλ6717,6731 together in the second fitting process. As a reference, there are 5 out of 19 and 1 of 2 objects at z ≃ 3.3 and z = 2.2 with detected [Ne iii]λ3869 and [S ii]λ6717 at >3σ significance, respectively.

In the following sections, we focus on 19 spectroscopically confirmed EELGs at 3.1 < z < 3.6 unless explicitly mentioned. Their positions are highlighted in Figure 1, and the distributions of Δm_Ks and spectroscopic redshifts are shown in Figure 2. The median spectroscopic redshift of them is 3.3.

We carried out an SED fitting to broadband photometry to obtain primarily stellar masses. Subaru/HSC Y-band data are also used in addition to the photometric bands used to derive Δm_Ks. Photometry is converted to the total magnitudes from 3'' aperture magnitudes following Appendix A.2 in Laigle et al. (2016). Upper limits are assigned when the fluxes are below 3σ significance.

SED fitting was performed by using the 2018.0 version of Code for Investigating GALaxy Emission (CIGALE;¹² Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019). As SED templates, we used composite stellar population models generated from the simple stellar population models of Bruzual & Charlot (2003, hereafter BC03) with a Chabrier initial mass function (IMF; Chabrier 2003) with lower and upper mass cutoffs of 0.1 and 100 M_⊙, respectively. Stellar population age ranges in log age/yr = 7–9.6 with steps of 0.1 dex. The upper limit of the age is assumed not to exceed the age of the universe at z = 3. Metallicities are allowed to have Z = 0.004, 0.008, and 0.02.

We consider delayed-τ models, $\mathrm{SFR}\propto t\exp (-t/\tau )$, for the star formation history (SFH) with $\mathrm{log}\tau /\mathrm{yr}$ = 8–10 with steps of 0.1 dex. Recent studies of SFH of EELGs at low and high redshifts suggest that they are in a starburst phase. Telles & Melnick (2018) studied SFH of local H ii galaxies with a three-burst SFH and found that the bulk of stellar mass of them has been formed by the past star formation episode, while they are currently in a maximum starburst phase. At higher redshifts similar to our sample, Cohn et al. (2018) found that EELGs at 2.5 < z < 4 show evidence of a starburst in the most recent 50 Myr with rising SFH in the past 1 Gyr. These results suggests that delayed-τ models are a more representative SFH for EELGs than a constant and exponentially declining SFH. Note that in the later analysis we will only use the stellar mass from the SED fitting and the stellar mass is a robust parameter against the assumption of SFH.

For galaxies with intense emission lines like the EELGs studied here, inclusion of emission lines has critical importance (e.g., Schaerer et al. 2013; Stark et al. 2013; Onodera et al. 2016) to extract physical parameters from SED fitting. For example, changes in stellar masses can be ≳0.5 dex when the contribution of Hβ + [O iii]λλ4959,5007 fluxes is ≳50% in the Ks band (Onodera et al. 2016). Instead of subtracting emission-line contributions from broadband photometry (Onodera et al. 2016), we include nebular emissions as supported in CIGALE. We assumed an ionization parameter of $\mathrm{log}U=-2.5$, which is a typical value for SFGs on the star-forming main sequence (MS) at z ∼ 3.3 (Onodera et al. 2016). LyC photons are assumed to be entirely absorbed by neutral hydrogen, i.e., zero LyC escape fraction (f_esc) and no LyC absorption by dust. As we will discuss later, some of our EELGs at z ≃ 3.3 are likely to have nonzero f_esc. However, the photometric data from the COSMOS2015 catalog do not allow us to put meaningful constraints on the f_esc with CIGALE. Indeed, the output parameters used in the following discussion are consistent well within the corresponding errors even if we set f_esc as a free parameter. Therefore, we simply adopt the fitting results with f_esc = 0. Nebular contributions (i.e., emission lines and continuum) are calculated following Inoue (2011) in CIGALE. Although the line width is not important for broadband data, we set the FWHM of emission lines as 300 km s⁻¹.

In CIGALE dust attenuation is implemented as a modified starburst law, k_λ, which is based on the Calzetti curve (Calzetti et al. 2000), ${k}_{\lambda }^{\mathrm{starburst}}$, with flexibilities to add a UV bump and alter the overall slope as follows:

$\begin{eqnarray}&&{k}_{\lambda }={D}_{\lambda }+{k}_{\lambda }^{\mathrm{starburst}}{\left(\displaystyle \frac{\lambda }{550\mathrm{nm}}\right)}^{\delta }\displaystyle \frac{E{\left(B-V\right)}_{\delta =0}}{E{\left(B-V\right)}_{\delta }},\end{eqnarray} \tag{ 3 }$

where D_λ is the Drude profile to express a UV bump and δ modifies the slope (Noll et al. 2009). We assumed the average SMC Bar extinction curve from Gordon et al. (2003), which can be approximated by setting D_λ = 0 and δ = −0.62. The nebular extinction E(B − V)_neb is allowed to vary between 0 and 0.8 with steps of 0.01 in the fitting. Because nebular emission lines are explicitly taken into account in the fit, one needs to make further assumptions on the attenuation curve and the relation between attenuation affecting stellar continuum and nebular emission. Here we also used the same SMC Bar curve (Gordon et al. 2003) for nebular emission lines and assumed E(B − V)_neb = 3.06 E(B − V)_star. The latter relationship has been recently obtained by Theios et al. (2019) for the case of the SMC attenuation curve using a sample of SFGs at 2 < z < 2.7 (Steidel et al. 2014; Strom et al. 2017).

At z ≳ 3, absorption by the intergalactic medium (IGM) in the observed optical bands cannot be ignored (e.g., Madau 1995; Inoue & Iwata 2008). When shifting the template SEDs to the observed redshift, CIGALE applies the IGM absorption using the prescription of Meiksin (2006).

The best-fit values and their corresponding standard deviations are computed based on probability distribution functions implemented as the pdf_analysis module (see Noll et al. 2009 and Boquien et al. 2019 for the details). The observed and best-fit SEDs are shown in Figure 5, and the resulting stellar masses, UV spectral slopes (see Section 4.2), attenuation parameters, SFRs, and stellar ages for EELGs at z > 3 are listed in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Stellar Mass, Attenuation, Star Formation Rate, and Stellar Age

ID	$\mathrm{log}{M}_{\star }$	βUV,CIGALE	βUV,UV	E(B − V)star,CIGALE	E(B − V)star,UV	log SFRCIGALE	log SFRUV	log AgeCIGALE
	(M⊙)			(mag)	(mag)	(M⊙ yr−1)	(M⊙ yr−1)	(yr)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
312786	9.58 ± 0.15	−1.78 ± 0.17	−1.84 ± 0.43	0.06 ± 0.02	0.07 ± 0.04	1.22 ± 0.26	1.24 ± 0.09	8.79 ± 0.27
314763	9.55 ± 0.20	−1.17 ± 0.14	−1.24 ± 0.22	0.14 ± 0.02	0.12 ± 0.02	1.95 ± 0.25	1.70 ± 0.05	8.14 ± 0.39
324921	9.83 ± 0.14	−1.28 ± 0.15	−1.41 ± 0.25	0.11 ± 0.02	0.11 ± 0.02	1.66 ± 0.18	1.60 ± 0.06	8.60 ± 0.27
337592	8.29 ± 0.19	−2.13 ± 0.13	−2.16 ± 0.54	0.05 ± 0.02	0.04 ± 0.05	1.41 ± 0.15	0.96 ± 0.12	7.34 ± 0.64
340544	8.95 ± 0.31	−1.89 ± 0.13	−1.77 ± 0.33	0.07 ± 0.02	0.07 ± 0.03	1.74 ± 0.22	1.38 ± 0.07	7.88 ± 0.72
350488	10.36 ± 0.08	−0.20 ± 0.37	−1.16 ± 0.73	0.18 ± 0.05	0.13 ± 0.06	1.23 ± 0.24	1.10 ± 0.16	9.12 ± 0.16
796059	9.99 ± 0.10	−1.64 ± 0.13	−1.89 ± 0.26	0.07 ± 0.02	0.06 ± 0.02	1.36 ± 0.10	1.39 ± 0.06	8.97 ± 0.20
797658	9.02 ± 0.24	−1.46 ± 0.19	−1.56 ± 0.51	0.11 ± 0.03	0.09 ± 0.05	1.47 ± 0.31	1.17 ± 0.11	8.20 ± 0.52
798507	8.88 ± 0.14	−1.71 ± 0.10	−1.58 ± 0.25	0.10 ± 0.01	0.09 ± 0.02	1.90 ± 0.15	1.53 ± 0.06	7.43 ± 0.43
809877	8.46 ± 0.07	−2.26 ± 0.07	−2.31 ± 0.27	0.04 ± 0.01	0.03 ± 0.02	1.66 ± 0.09	1.19 ± 0.06	7.15 ± 0.24
810838	9.91 ± 0.13	−1.37 ± 0.18	−1.52 ± 0.35	0.10 ± 0.02	0.10 ± 0.03	1.45 ± 0.16	1.46 ± 0.08	8.84 ± 0.23
818357	8.75 ± 0.32	−1.30 ± 0.26	−1.40 ± 0.85	0.13 ± 0.03	0.11 ± 0.08	1.56 ± 0.25	1.07 ± 0.19	7.98 ± 0.87
870659	8.53 ± 0.25	−1.66 ± 0.17	−2.01 ± 0.76	0.10 ± 0.02	0.05 ± 0.07	1.54 ± 0.19	0.94 ± 0.16	7.61 ± 0.92
871455	8.87 ± 0.26	−2.14 ± 0.11	−2.06 ± 0.30	0.05 ± 0.02	0.05 ± 0.03	1.55 ± 0.21	1.27 ± 0.07	7.91 ± 0.54
891150	9.82 ± 0.16	−1.39 ± 0.14	−1.34 ± 0.21	0.10 ± 0.02	0.11 ± 0.02	1.79 ± 0.18	1.77 ± 0.05	8.48 ± 0.29
893160	8.30 ± 0.17	−2.27 ± 0.10	−2.16 ± 0.48	0.04 ± 0.01	0.04 ± 0.04	1.32 ± 0.16	1.00 ± 0.10	7.45 ± 0.49
894641	9.03 ± 0.27	−1.86 ± 0.23	−1.24 ± 0.76	0.06 ± 0.03	0.12 ± 0.07	1.20 ± 0.39	1.23 ± 0.16	8.55 ± 0.46
919252	8.72 ± 0.15	−2.06 ± 0.10	−2.17 ± 0.27	0.06 ± 0.01	0.04 ± 0.02	1.75 ± 0.15	1.30 ± 0.06	7.42 ± 0.46
919704	10.57 ± 0.06	−1.39 ± 0.18	−1.03 ± 0.57	0.04 ± 0.03	0.14 ± 0.05	0.79 ± 0.21	1.32 ± 0.13	9.16 ± 0.11
Low-mass composite	8.74 ± 0.33	−1.89 ± 0.37	−2.01 ± 0.35	0.06 ± 0.03	0.05 ± 0.03	1.55 ± 0.21	1.19 ± 0.18	7.61 ± 0.42
High-mass composite	9.87 ± 0.31	−1.38 ± 0.23	−1.37 ± 0.27	0.10 ± 0.05	0.11 ± 0.02	1.40 ± 0.32	1.42 ± 0.27	8.82 ± 0.31

Note. Column (1): object ID. Column (2): stellar mass from SED fitting. Column (3): UV β slope from the best-fit SED. Column (4): UV β slope from broadband photometry. Column (5): E(B − V) for stellar continuum from SED fitting. Column (6): E(B − V) for stellar continuum from UV β slope. Column (7): SFR from the best-fit SED. Column (8): attenuation-corrected SFR based on UV luminosity. Column (9): stellar population age from the best-fit SED.

Download table as:  ASCIITypeset image

As seen in Figure 5, most of the EELGs show a flat or power-law continuum in the optical to near-IR bands corresponding to young stellar populations with an excess flux in the K band due to the intense [O iii]λ5007 line. There are a couple objects, namely, 350488 and 919704, showing redder SEDs than the others. These two objects indeed have the best-fit age of ≳1 Gyr, indicating a significant contribution of old stellar populations.

It turned out that CIGALE underestimates [O iii]λ5007 fluxes and EWs of the EELG sample on average by a factor of two in our fitting run. To check this effect on the estimate of the SED properties, especially stellar masses, we ran CIGALE without Ks-band data while keeping other inputs identical and found that stellar masses from the no-Ks run are larger than those from our fiducial run up to 0.2 dex with a median of 0.04 dex. The difference tends to be larger at lower stellar masses ($\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$) with a median of 0.08 dex. These differences are, however, smaller than the uncertainties of stellar masses, and the results and discussion presented later will remain unchanged.

The underestimated best-fit predictions of [O iii] fluxes and EWs could be an indication of higher ionization parameters in EELGs at z ≃ 3.3 than $\mathrm{log}U=-2.5$ assumed in the CIGALE run. Although $\mathrm{log}U\gt -2.0$ seems unrealistically high for SFGs (e.g., Yeh & Matzner 2012; Shirazi et al. 2014; Strom et al. 2018), we carried out a CIGALE run with a fixed $\mathrm{log}U=-1$, which is a theoretical maximum (Yeh & Matzner 2012). It is found, however, that the predicted [O iii]λ5007 properties are still underestimated, though the discrepancies become smaller. Therefore, the issue cannot be solved solely by the higher ionization parameter.

The lower limit of the stellar population age in the SED fitting could also limit the [O iii]λ5007 properties. To check this possibility, we run CIGALE by adding 1–10 Myr stellar populations. We found that the best-fit EW([O iii]λ5007) does not always become larger by including younger stellar populations, but a clear improvement is seen for those with highest EW([O iii]λ5007), and the best-fit ages of <10 Myr are often derived for such cases. Similar to other CIGALE runs discussed above, there is only little change in the best-fit stellar mass (0.01 dex in median). Trying to find an exact match of the EWs and fluxes of the [O iii]λ5007 emission line from our EELG sample is not a focus of this study, but it appears that higher ionization parameters and younger stellar population ages are likely to play key roles in intensifying the emission-line strength (see also Section 5.6; De Barros et al. 2019).

Measurements of the amount of dust attenuation and SFR are carried out by using the UV part of the SED, namely, the UV spectral slope β_UV defined as f_λ ∝ ${\lambda }^{{\beta }_{\mathrm{UV}}}$ and attenuation-corrected far-UV (FUV) luminosity, respectively. We derived β_UV in the same way as Onodera et al. (2016). Briefly, we fit the total broadband photometry in r, ip, zpp, Y, and J bands with a linear function in the magnitude–$\mathrm{log}\lambda$ space to obtain FUV (λ ≃ 1530 Å) and near-UV (λ ≃ 2300 Å) magnitudes to compute β_UV following the prescription presented by Pannella et al. (2015). In the fitting, we put an additional weight of $1+\left|\lambda -{\lambda }_{\mathrm{FUV}}\right|/{\lambda }_{\mathrm{FUV}}$ in order to give more weights to filters closer to the rest-frame FUV band (Nordon et al. 2013). This procedure was repeated 5000 times by artificially perturbing the photometry with the corresponding errors. The median values and standard deviations based on the median absolute deviation (σ_MAD) of the resulting β_UV and rest-frame FUV magnitude distributions are considered as the values and associated 1σ errors, respectively.

The β_UV was then transformed into E(B − V)_star by using a recent calibration by Reddy et al. (2018a) for the SMC extinction curve (Gordon et al. 2003),

$\begin{eqnarray}&&{\beta }_{\mathrm{UV}}=-2.616+11.259\,E{(B-V)}_{\mathrm{star}}.\end{eqnarray} \tag{ 4 }$

The intrinsic slope of β_UV = −2.616 is derived by assuming a simple stellar population of the Binary Population and Spectral Synthesis (BPASS) models (Eldridge et al. 2017) with a constant star formation with an age of 100 Myr and metallicity of Z = 0.002, with a two-power-law IMF with α = 2.35 for M_⋆ > 0.5 M_⊙ and α = 1.30 for 0.1 < M_⋆/M_⊙ < 0.5. The adopted relationship in Equation (4) provides E(B − V)_star closer to that of the original calibration by Calzetti et al. (2000) with discrepancies within ≃0.1 mag, while using the relations derived by the same authors (Reddy et al. 2015) for the Calzetti et al. (2000) and Reddy et al. (2015) attenuation curves provides systematically higher E(B − V)_star than that from the SMC curve by 0.03–0.2 mag.

By using the derived E(B − V)_star, SFRs were computed from attenuation-corrected rest-frame FUV (1500 Å) luminosities (L₁₅₀₀) using the calibration by Theios et al. (2019),

$\begin{eqnarray}&&\mathrm{logSFR}/({M}_{\odot }\,{\mathrm{yr}}^{-1})=\mathrm{log}\nu {L}_{\nu }(1500\,\mathring{\rm A} )/(\mathrm{erg}\,{{\rm{s}}}^{-1})-43.46.\end{eqnarray} \tag{ 5 }$

This relation assumes the same intrinsic stellar population model as used to derive E(B − V)_star from β_UV, and it is slightly different from that of Kennicutt & Evans (2012), where the constant term is 43.35.

The attenuation properties and SFRs derived here are shown in Table 3. These properties are compared with those derived from SED fitting in Figure 6. Both estimates agree well for β_UV, E(B − V)_star, and L_ν(1500 Å), with a few outliers that are caused mainly by large uncertainties in the rest-frame UV photometric data. SFRs derived by CIGALE SED fitting are systematically higher than those based on UV luminosities, with a median difference of 0.28 dex. The systematic difference is likely due to different assumptions in the intrinsic SED and dust attenuation. In the SED fitting, we used Bruzual & Charlot (2003) models, while BPASS models (Eldridge et al. 2017) are used in the calibration by Reddy et al. (2018a) in the UV-based estimate. Moreover, the relation between E(B − V)_star and E(B − V)_neb applied for the CIGALE run is derived assuming BPASS models and a revised version of the SMC extinction curve (Reddy et al. 2016a; Theios et al. 2019). Reconciling the two estimates is beyond the scope of this paper, and we will use only stellar masses from the SED fitting and use UV-based estimates for the attenuation and SFR parameters in the following analysis.

Zoom In Zoom Out Reset image size

Because our EELG sample does not show continuum in the MOIRCS spectra in general (Figure 3), emission-line EWs cannot be measured directly on the spectra. Furthermore, since Ks broadband fluxes are dominated by [O iii] + Hβ up to ≃100%, estimated continuum levels based on the Ks-band photometry can be comparable to the photometric error. Therefore, we use the best-fit stellar SED derived in Section 4.1 assuming no error in the continuum. The model continuum is derived as the median flux at 4940–5028 Å in the rest-frame wavelength. The derived EW([O iii]λ5007) is shown in Table 4. It should be noted that the assumption of the error-free continuum is fairly strong. If we assume the same errors as the total Ks-band fluxes for the continuum, the resulting uncertainties of EW([O iii]λ5007) increase from ≃0.01–0.02 dex to ≃0.04 dex, which does not change any of our results and conclusions presented in this study.

Table 4.  Ionization Properties

ID	EW([O iii] λ5007)	O32	R23	$\mathrm{log}{\xi }_{\mathrm{ion},0}$
	(Å)			(erg−1 Hz)
(1)	(2)	(3)	(4)	(5)
312786	223 ± 11	>4.5	>10.5	<25.0
314763	217 ± 9	1.9 ± 0.7	6.5 ± 1.7	25.30 ± 0.07
324921	250 ± 9	1.5 ± 0.6	9.6 ± 2.8	25.33 ± 0.08
337592	2051 ± 82	>7.2	4.8 ± 3.1	25.86 ± 0.15
340544	1391 ± 57	2.9 ± 1.3	9.4 ± 3.4	25.63 ± 0.09
350488	233 ± 11	>4.3	>8.5	<25.7
796059	407 ± 16	2.7 ± 1.2	>6.5	<25.7
797658	1123 ± 33	⋯	⋯	<26.0
798507	557 ± 30	2.5 ± 1.0	8.1 ± 2.8	25.30 ± 0.11
809877	1853 ± 64	⋯	⋯	<25.8
810838	262 ± 14	1.5 ± 0.7	5.2 ± 2.0	25.70 ± 0.11
818357	670 ± 60	>3.0	>4.6	<25.7
870659	570 ± 43	>4.2	8.8 ± 8.1	25.21 ± 0.23
871455	472 ± 29	>5.8	>9.0	<25.0
891150	153 ± 13	1.3 ± 0.4	>11.2	<25.0
893160	1171 ± 58	>5.3	5.0 ± 2.8	25.54 ± 0.13
894641	572 ± 22	>2.6	7.7 ± 6.8	25.32 ± 0.20
919252	860 ± 27	3.3 ± 1.3	14.7 ± 5.0	25.15 ± 0.10
919704	133 ± 5	>2.5	5.7 ± 3.8	25.58 ± 0.16
Low-mass composite	791 ± 360	5.7 ± 2.5	10.8 ± 3.7	25.32 ± 0.18
High-mass composite	189 ± 30	2.7 ± 0.5	9.6 ± 1.5	25.43 ± 0.27

Note. Column (1): object ID. Column (2): rest-frame EW of [O iii]λ5007. Column (3): O₃₂ as defined in Equation (8). Column (4): R₂₃ as defined in Equation (7). Column (5): ξ_ion,0 defined in Equation (9) assuming f_esc = 0.

Download table as:  ASCIITypeset image

In Figure 7, we compare EW([O iii]λ5007) with EW([O iii] + Hβ) estimated from Δm_Ks and that predicted with a conversion derived by Reddy et al. (2018b) for a sample of SFGs at z ≃ 3.4. The former is derived by using Equation (1), while the latter is estimated as

$\begin{eqnarray}&&\mathrm{EW}([{\rm{O}}\,{\rm\small{III}}]+{\rm{H}}\beta )=c\,\mathrm{EW}([{\rm{O}}\,{\rm\small{III}}]).\end{eqnarray} \tag{ 6 }$

The correction factor c is defined as c = 0.581–0.0074x + 6.322 × 10⁻³x², where $x=\mathrm{log}{M}_{\star }/{M}_{\odot }$ (Reddy et al. 2018b). Although about half of our sample has lower stellar mass than the range of $9.0\lt \mathrm{log}M/{M}_{\odot }\lt 11.5$ for which the conversion is derived, the difference between the two estimates agrees well with a median of 0.06 dex and a standard deviation of 0.13 dex, which is smaller than the ≃0.2–0.3 dex scatter of the sample distribution presented in Reddy et al. (2018b).

Zoom In Zoom Out Reset image size

We derived R₂₃ and O₃₂ indices defined as

$\begin{eqnarray}&&{R}_{23}\equiv \displaystyle \frac{[{\rm{O}}\,{\rm\small{II}}]\lambda 3727+[{\rm{O}}\,{\rm\small{III}}]\lambda \lambda 4959,5007}{{\rm{H}}\beta }\end{eqnarray} \tag{ 7 }$

and

$\begin{eqnarray}&&{O}_{32}\equiv \displaystyle \frac{[{\rm{O}}\,{\rm\small{III}}]\lambda \lambda 4959,5007}{[{\rm{O}}\,{\rm\small{II}}]\lambda 3727},\end{eqnarray} \tag{ 8 }$

respectively. All emission-line fluxes are corrected for dust extinction using β_UV-based E(B − V)_neb (see Section 4.2) with a conversion E(B − V)_neb = 3.06 E(B − V)_star for the average SMC Bar extinction curve (Gordon et al. 2003; Theios et al. 2019).

We assumed [O iii]λ5007/[O iii]λ4959 = 3 to estimate the total [O iii] flux. When [O ii]λ3727 emission lines are not detected, we assign fluxes corresponding to the 3σ upper limit to derive the lower limit of O₃₂, while undetected [O ii]λ3727 fluxes are assumed to be zero for R₂₃ because they are essentially negligible compared to those of [O iii]λ5007. These upper limits are also corrected for dust extinction as described above.

Detected Hβ fluxes are also corrected for stellar absorption by using the best-fit SED derived in Section 4.1. The correction factors are ≲10% for less massive galaxies with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lt 10$, while one object with $\mathrm{log}{M}_{\star }/{M}_{\odot }\simeq 10.5$ has a correction factor of ≃40%. The latter value is slightly larger than the prediction using the correction factor–stellar mass relation derived by Zahid et al. (2014) for a sample of SFGs at z ≃ 1.6.

The derived O₃₂ and R₂₃ values are listed in Table 4.

The LyC photon production efficiency, ξ_ion, is defined as the ratio of the production rate of LyC photons, N(H⁰), to the attenuation-corrected UV continuum luminosity density, L_UV,corr:

$\begin{eqnarray}&&{\xi }_{\mathrm{ion}}=\displaystyle \frac{N({{\rm{H}}}^{0})}{{L}_{\mathrm{UV},\mathrm{corr}}}\left({{\rm{s}}}^{-1}/(\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1})\right).\end{eqnarray} \tag{ 9 }$

The UV luminosity is derived at rest frame 1500 Å using broadband photometry (Section 4.2). For N(H⁰), we adopt the relation to the intrinsic Hα luminosity, L(Hα), by Leitherer & Heckman (1995),

$\begin{eqnarray}&&L({\rm{H}}\alpha )\left(\mathrm{erg}\,{{\rm{s}}}^{-1}\right)=1.36\times {10}^{-12}\,N({{\rm{H}}}^{0})\left({{\rm{s}}}^{-1}\right).\end{eqnarray} \tag{ 10 }$

Then, we assume L(Hα)/L(Hβ) = 2.86 assuming case B recombination with an electron density of 100 cm⁻³ and temperature of 10⁴ K (Osterbrock & Ferland 2006) to compute ξ_ion for our EELG sample at z ≃ 3.3. The Hβ fluxes are corrected for dust attenuation and stellar absorption as described in Section 4.4. One should keep in mind that the relation above is derived for an ionization-bounded H ii region with no LyC escape.

LyC escape fraction is suggested to be correlated with some of the galaxy properties such as O₃₂ (e.g., Nakajima & Ouchi 2014; Faisst 2016; Izotov et al. 2018a; Nakajima et al. 2020), E(B − V) (e.g., Reddy et al. 2016b), and Lyα profile (e.g., Verhamme et al. 2015; Dijkstra et al. 2016; Izotov et al. 2020). Although we have measurements of O₃₂ and E(B − V) in hand, we do not correct f_esc here, because the correlations between f_esc and these parameters show large scatters especially at large O₃₂ and low E(B − V) (e.g., Faisst 2016; Reddy et al. 2016b; Nakajima et al. 2020). Instead, we use ξ_ion,0 for the case of zero LyC escape fraction (Bouwens et al. 2016; Nakajima et al. 2016) to explicitly indicate the assumption of no LyC escape from the H ii region in Equation (10). The intrinsic ξ_ion then can be derived by ξ_ion = ξ_ion,0/(1 − f_esc). We will investigate the relationship between ξ_ion,0 and other galaxy properties in Section 5.4.3. The derived ξ_ion,0 are shown in Table 4.

In order to investigate average spectroscopic properties of EELGs at z ≃ 3.3, we create composite Subaru/MOIRCS spectra. We split the sample in the stellar mass at $\mathrm{log}{M}_{\star }/{M}_{\odot }\simeq 9.2$. There are 11 and 8 objects in the low-mass and high-mass bins, respectively. Before the stacking, each 1D spectrum is shifted and linearly interpolated to the rest-frame wavelength grid of 1 Å interval, while the corresponding noise spectrum is interpolated to the same wavelength grid in quadrature. Then, composite spectra are constructed by taking an average at each wavelength weighted by the inverse variance. The associated noise spectra are derived by the standard error propagation from the individual noise spectra. Emission-line fluxes are measured by fitting Gaussian components for emission lines and linear continuum similar to that applied to individual objects (Section 3.2). The resulting stacked spectra are shown in Figure 8, and emission-line fluxes are shown in Table 2.

Zoom In Zoom Out Reset image size

SED properties for the composite spectra are derived by taking the median and σ_MAD of the sample in each bin. EWs of [O iii]λ5007, emission-line ratios, and ξ_ion,0 of the composite spectra are derived by using the derived SED properties as described above and emission-line fluxes measured on the composite spectra. These parameters are appended in Tables 3 and 4.

The COSMOS2015 catalog provides photometric redshifts with a precision of σ_Δz/(1 + z) = 0.021, with a catastrophic failure of 13.2% for SFGs at z > 3 (Laigle et al. 2016). Although nebular emission-line contributions of Lyα, [O ii]λ3727, Hβ, [O iii]λλ4959,5007, and Hα have already been taken into account when deriving the photometric redshifts, their emission-line ratios are fixed to representative values for normal SFGs. For example, while they fixed O₃₂ = 0.36, our EELG sample at z ≃ 3.3 shows a significantly higher O₃₂ with the median of ≳3. Therefore, it is worth checking whether the photometric redshift and physical parameters in the COSMOS2015 catalog work for our EELG sample or not.

In the left panel of Figure 9, we compare the COSMOS2015 photometric redshifts with spectroscopic ones for our identified EELG sample at z > 3. Using the same metric as Laigle et al. (2016), we compute ${\sigma }_{{{\rm{\Delta }}}_{z/(1+z)}}=0.025$ with a systematic offset of −0.01 and no catastrophic failure. Two of our spectroscopic samples are at z_phot ≃ 2.2, and their spectroscopic redshifts show almost exact agreement. There are two objects without emission-line detections. Since their COSMOS2015 photometric redshifts are 2.476 and 3.774, respectively, the primary emission lines (Hα and [O iii]λ5007, respectively) may fall outside of the spectral coverage of MOIRCS. Here we conclude that the photometric redshift in COSMOS2015 for EELGs at z ≃ 3.3 is as precise as those for other galaxies at the similar redshift despite the fact that the emission-line ratios assumed in the COSMOS2015 catalog are significantly different from those in our sample. Note, however, that the spectroscopic follow-up was carried out for a sample with visually reliable SED shapes and no attempt was done to search for EELGs at z ≃ 3.3 with z_phot ≲ 2 and z_phot ≳ 4.

Zoom In Zoom Out Reset image size

Stellar masses are compared in the middle panel of Figure 9. COSMOS2015 stellar masses are systematically larger than those derived by our SED fitting, especially evident at lower stellar masses of $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$. We do not find any indications that the difference between two stellar mass estimates correlates with overestimates of photo-z. Rather, this can be understood by a realistic treatment of the emission-line strengths and ratios especially of [O iii] to the Ks-band fluxes in CIGALE using the recipe of Inoue (2011) compared to the empirical line ratios in the COSMOS2015 SED fitting. Therefore, in order to estimate more appropriate stellar masses, one may not be able to rely on the catalog values for low-mass ($\mathrm{log}{M}_{\star ,\mathrm{COSMOS}2015}/{M}_{\odot }\lesssim 9.5$) galaxies at z ≳ 3 and needs to perform a dedicated SED fitting depending on the subpopulation of interest.

On the other hand, SFRs between UV-based and COSMOS2015 estimates agree well. At high SFR of log SFR ≳ 1.5, our UV-based SFRs seem to be slightly smaller than those in the COSMOS2015 catalog. As shown in Figure 6, UV-based SFRs tend to be lower than those derived in the CIGALE SED fitting. Therefore, if SFR_UV is replaced by SFR_CIGALE in Figure 9, points will move toward the right and SFR_CIGALE would be higher than SFR_COSMOS2015 by ≃0.3 dex.

Note also that the SED fittings in the COSMOS2015 catalog and this study use different assumptions on the set of intrinsic SED templates, dust attenuation recipes, and the explored ranges of age, metallicity, dust attenuation, and ionization parameters, which could introduce systematics for the estimate of physical parameters.

Intense [O iii] emission with a large EW and large [O iii]/Hβ ratio can originate from active galactic nuclei (AGNs) as well as SFGs (e.g., Baskin & Laor 2005; Kewley et al. 2006; Caccianiga & Severgnini 2011). Since our primary motivation in this study is to search for a population of EELGs at z ≳ 3 as an analog of SFGs in the EoR, we would like to check whether AGNs dominate the sample or not.

Among 240 objects in the parent sample, 142 objects have 3 < z_phot < 3.7 and none of them show counterparts in X-ray in the Chandra COSMOS Legacy Survey catalog (Civano et al. 2016; Marchesi et al. 2016) and in radio in the public VLA data.

There are four EELG candidates showing Spitzer/MIPS 24 μm detections ranging from 70 to 240 μJy. The 24 μm fluxes roughly correspond to SFR ≳ 1000 M_⊙ yr⁻¹ using the Dale & Helou (2002) templates to convert 24 μm fluxes to total IR luminosities and the Kennicutt (1998) relation to convert total IR luminosities to SFRs (e.g., Wuyts et al. 2008; Marchesini et al. 2010). If the 24 μm fluxes are due to star formation, they are well above the MS at z ≃ 3.3 (e.g., Speagle et al. 2014), and they would be detected in far-IR to submillimeter wavelength. However, none of the 140 EELG candidates at 3 < z_phot < 3.7 are detected in 100–500 μm in the Herschel/PACS PEP survey (Lutz et al. 2011) and Herschel/SPIRE HerMes survey (Oliver et al. 2012). Therefore, these MIPS 24 μm sources can be hosts of obscured AGNs, but they are rare among the selected EELG candidates, only ≲5%. Although the discussion above does not rule out the presence of low-luminosity or less obscured AGNs undetected in 24 μm, based on the available information for emission-line properties, it is unlikely that AGNs dominate the EELG population at z ≃ 3.3 as we discuss below.

The emission-line widths of [O iii]λ5007 of identified EELGs at z ≃ 3.3 are FWHM ≃ 400–500 km s⁻¹, which is actually narrower than the expected resolution of 960 km s⁻¹ of the MOIRCS HK500 grism. This could be due to the combination of better seeing size than the slit width and intrinsically compact nature of the objects. Since none of the members of our EELG sample at z ≃ 3.3 show a broad [O iii] line exceeding FWHM ≃ 1000 km s⁻¹ indicative of AGN-driven outflow (e.g., Holt et al. 2008; Alexander et al. 2010; Genzel et al. 2014), the line widths are rather consistent with star formation.

Emission-line ratios can also be used to diagnose the ionization mechanism through, e.g., the BPT diagram (Baldwin et al. 1981; Kewley et al. 2001; Kauffmann et al. 2003). However, for our EELGs at z ≃ 3.3, only Hβ and [O iii] lines are available. As an alternative to the BPT diagram, Juneau et al. (2011) proposed to use the mass–excitation (MEx) diagram by using [O iii]λ5007/Hβ and stellar mass. Figure 10 shows the MEx diagram for our EELGs at z ≃ 3.3 and Sloan Digital Sky Survey (SDSS) DR7 objects taken from the OSSY catalog (Oh et al. 2011), with the revised demarcation lines to classify star formation, composite, and AGN (Juneau et al. 2014). In the MEx diagram, most of the individual EELGs and the low-mass composite are in the region of either star-forming or composite galaxies, but a couple of objects, as well as the high-mass composite, are classified as Seyfert galaxies. These two objects are the most massive ones among the spectroscopically confirmed objects, and they show noticeably redder SEDs with old stellar populations of ≳1 Gyr as seen in Figure 5 and Table 3. Moreover, they show significantly lower SFRs relative to the MS of SFGs at z ≃ 3.3 (see Section 5.3 and Figure 11). Therefore, they could be obscured AGNs despite no indication in the other diagnostics. On the other hand, one MIPS 24 μm source is actually classified as an SFG. Note that the demarcation lines of Juneau et al. (2014) are, however, calibrated by using local galaxies from SDSS. Strom et al. (2017) showed that SFGs at z ≃ 2.3 from the KBSS-MOSFIRE survey distribute toward higher [O iii]/Hβ ratios up to 0.8 dex at a given stellar mass in the MEx diagram. Our sample of objects classified as AGNs by Juneau et al. (2014) lines are indeed consistent with the KBSS-MOSFIRE star-forming sample at z ≃ 2.3, as well as AGNs when considering the lower limits more conservatively (Strom et al. 2017).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, looking at various possible diagnostics of AGNs for our sample, we conclude that EELGs at z ≃ 3.3 are most likely to be dominated by star formation, though we cannot completely rule out the presence of AGNs at the massive end (M_⋆ ≳ 10^10.5 M_⊙). In the following analysis, we consider all the spectroscopically identified EELGs at z ≃ 3.3 as SFGs.

Figure 11 shows the relation between SFR and stellar mass for our EELG sample at z ≃ 3.3, together with normal SFGs at z ≃ 3.3 selected via UV continuum brightness or expected Hβ flux (Onodera et al. 2016) and Lyα emitters (LAEs) at z = 3–4 studied by Nakajima et al. (2016). Objects taken from Nakajima et al. (2016) and Onodera et al. (2016) distribute around the fiducial MS of SFGs at z ≃ 3.3 (Speagle et al. 2014). While our spectroscopically identified EELGs at z ≃ 3.3 are on the MS at $\mathrm{log}M/{M}_{\odot }\gtrsim 9.5$, those with $\mathrm{log}M/{M}_{\odot }\lesssim 9.0$ clearly show higher SFR than the MS by ≳0.5–1 dex.

The elevated SFR relative to the MS in less massive galaxies is also seen in EELGs at 1.3 < z < 2.3 identified by Maseda et al. (2014), as also shown in Figure 11 with light-green circles. They conclude that the star formation in low-mass galaxies producing intense [O iii] and Hα emission lines is undergoing a burst-like stage with a timescale of ∼50 Myr. The mass-weighted ages from the SED fitting with CIGALE for our EELG sample are indeed young, ≃10–100 Myr, at $\mathrm{log}M/{M}_{\odot }\lesssim 9.0$. Such young stellar population ages for low-mass galaxies are also reported by van der Wel et al. (2011) for a sample of starbursting EELGs at z ≃ 1.7 with $\mathrm{log}M/{M}_{\odot }\simeq 8$. Cohn et al. (2018) carried out an analysis of SFH for EELGs at 2.5 < z < 4 and found evidence of a recent starburst within 50 Myr, which is also a consistent picture with other lower-redshift EELGs, as well as this study. Therefore, we conclude that low-mass EELGs at z ≃ 3.3 are also likely to be in a starburst phase with a short timescale of ≲100 Myr.

Although our low-mass EELGs at z ≃ 3.3 show similar nature to lower-redshift counterparts in the SFR–M_⋆ diagram and stellar ages, one should keep in mind that our EELG selection at a fainter magnitude of Ks^obs ≳ 23 may fail to select objects located between the EW cut of Δm_Ks = 0.3 and the 3σ cut shown as the dashed–dotted and upper solid lines, respectively, in Figure 1. This means that the systematic trend of elevated SFR for low-mass EELGs at z ≃ 3.3 relative to the MS could be at least partly due to a selection bias. However, it is not straightforward to quantify the bias for [O iii] emitters, because translating [O iii] flux to SFR (or Balmer line fluxes) is a complex function of various physical parameters, such as the gas-phase metallicity, ionization parameter, hardness of UV radiation field, and electron density (Kewley et al. 2013), and also Ks^obs is a combination of emission lines and stellar continuum (i.e., stellar mass). The [O iii]/Hα ratios derived by Faisst et al. (2016) using local analogs show a broad distribution of log [O iii]/Hα ≃ −0.4 to 0.4, which is also consistent with the range measured for SFGs at z ≃ 2 from the literature (e.g., Steidel et al. 2014; Shapley et al. 2015; Suzuki et al. 2016). On the other hand, as we will show in Section 5.4.1, our spectroscopically identified low-mass EELGs follow the EW([O iii]λ5007)–M_⋆ relation defined for more massive normal SFGs at z ≃ 3.3 (Reddy et al. 2018b). Therefore, it is not likely that we miss a large amount of low-mass EELGs with significantly lower EW([O iii]λ5007) because of increasing photometric errors at the fainter magnitude.

5.4.1. [O iii]λ5007 Equivalent Width and SED Properties

Figure 12 compares EW([O iii]λ5007) with SED properties, namely, stellar mass, SFR, and sSFR for our EELGs at z ≃ 3.3 (both individual and composite measurements) and other SFGs from the literature, namely, local SFGs selected from SDSS (Oh et al. 2011), GPs at z ≃ 0.2 (Cardamone et al. 2009), and normal SFGs on the MS at z ≃ 3.3 (Onodera et al. 2016). We also compare them with the best-fit relationship derived by Reddy et al. (2018b) for SFGs at z ≃ 2.3 and 3.3.

Zoom In Zoom Out Reset image size

At a fixed stellar mass, SFGs at z ≃ 3.3 from this study, Onodera et al. (2016), and Reddy et al. (2018b) show elevated EW([O iii]λ5007) by ≃1.5 dex compared to local SFGs selected from SDSS (Oh et al. 2011) as shown in the left panel of Figure 12. This strong redshift evolution of the EW([O iii]λ5007)–M_⋆ relation for normal SFGs has been reported before (e.g., Khostovan et al. 2016), and the amount of the evolution in our study is similar to their result. The linear relation for a sample of SFGs at z ≃ 3.4 is derived at $\mathrm{log}{M}_{\star }/{M}_{\odot }\gtrsim 9.5$ by Reddy et al. (2018b). They argued that the evolution of the EW([O iii]λ5007)–M_⋆ relation for the bulk of the SFG population can be explained by the redshift evolution of the SFR–M_⋆ and mass–metallicity relation. Our EELG sample at z ≃ 3.3 with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$ appears to follow the same relation extrapolated to lower stellar masses, indicating that EELGs may be no longer a rare galaxy population of SFGs at z ≳ 3, especially at the low stellar masses. On the other hand, EELGs are a rare population in the local universe. Cardamone et al. (2009) estimated the spatial number density of GPs as 2 deg⁻² with r < 20.5 mag. Only looking at GPs, their EW([O iii]λ5007)–M_⋆ relation is similar to that of EELGs and normal SFGs at z ≃ 3.3.

On the other hand, our EELG sample at z ≃ 3.3 shows up to ≃3 dex and ≃1.5 dex higher EW([O iii]λ5007) at a fixed SFR compared to those of local and z ≃ 2.3 SFGs, respectively (middle panel of Figure 12), while the trend for GPs is quite similar to our EELG sample showing elevated EW([O iii]λ5007) at a given SFR. Since Reddy et al. (2018b) used SFR based on Hα luminosity, they did not present the relationship at z ≳ 3, where Hα falls in a longer wavelength than the K-band coverage. At z ≃ 3.3, the MS does not seem to evolve much since z ≃ 2.3. For example, the formula provided by Speagle et al. (2014) suggests a ≃0.1 dex increase of the normalization of the MS from z = 2.3 to z = 3.3. Suzuki et al. (2015) indeed found no difference in the normalization of the MS between narrowband-selected Hα emitters at z = 2.2 and 2.5 and [O iii] emitters z = 3.2 and 3.6, respectively. Note that most of their narrowband-selected objects are normal SFGs at the corresponding redshifts, as the narrowband method selects those with less extreme emission-line strengths than the broadband excess method we applied. This is also confirmed by the distribution of normal SFGs on the MS (Onodera et al. 2016), which is on average consistent with the best-fit relation at z ≃ 2.3 by Reddy et al. (2018b). In addition, the average EW([O iii]λ5007) for SFGs shows little evolution at z ≃ 2–3 (Khostovan et al. 2016). Therefore, the elevated EW([O iii]λ5007) of our EELG sample, as well as GPs, does not seem to be explained simply by the evolution of a normal SFG population between z ≃ 3.3 and 2.3.

Combining the above two parameters, the right panel of Figure 12 shows the relationship between EW([O iii]λ5007) and sSFR. EELGs at z ≃ 3.3 show elevated EW([O iii]λ5007) at a fixed sSFR relative to normal SFGs at z ≃ 3.3 (Onodera et al. 2016) and those at z ≃ 2.3 (Reddy et al. 2018b), while the slope of the relationship does not appear to be very different between two epochs. This trend can also be seen for GPs. Those with lower sSFR of ≲10⁻⁸ yr⁻¹ appear to show more elevation of EW([O iii]λ5007) at a given sSFR than our EELG sample. A more spectroscopic sample of EELGs at z ≃ 3.3 would be required for the detailed comparison between high-redshift EELGs and local counterparts at this sSFR range. Low-mass objects in our EELG sample with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lt 9$ tend to cluster around the upper right part of the figure as also seen from the low-mass composite point. Normal SFGs at z ≃ 3.3 presented by Onodera et al. (2016) also contain objects with elevated sSFR of log sSFR/yr⁻¹ ≳ − 8 and EW([O iii]λ5007) of ≳500 Å, similar to our EELG sample. As seen in the left panel of the figure, most of them are low-mass objects with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$.

Sanders et al. (2020b) studied low-mass galaxies with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$ at z = 1.5–3.5 with detected auroral [O iii]λ4363 emission lines to derive robust metallicities via the direct method. Their sample also shows SFR above the MS at z ≃ 2.3 and elevated EW([O iii]λ5007) similar to our EELGs at z ≃ 3.3. They concluded that the auroral-line sample is younger and more metal-poor than normal galaxies on the MS. EELGs at 1.3 < z < 2.4 presented by Tang et al. (2019) also occupy similar parameter space to the auroral-line sample by Sanders et al. (2020b). Therefore, our EELGs at z ≃ 3.3, especially low-mass ones, are also likely to be made of galaxies that are young and low metallicity and would be promising candidates for a systematic [O iii]λ4363 follow-up spectroscopic observation. As an example, assuming an electron density of n_e = 500 cm⁻³ and an electron temperature of T_e = 15,000 K, which are reasonable values for high-redshift, low-metallicity SFGs (e.g., Shirazi et al. 2014; Kojima et al. 2017; Sanders et al. 2020b), [O iii]λ4363 flux is derived as ≃2% of the [O iii]λ5007 flux by using PyNeb¹³ (Luridiana et al. 2015). Scaling from a typical [O iii]λ5007 flux of our EELG sample of ≃10⁻¹⁶ erg s⁻¹ cm⁻², their [O iii]λ4363 flux is estimated to be ≃2 × 10⁻¹⁸ erg s⁻¹ cm⁻². Our Subaru/MOIRCS setup detected the aforementioned [O iii]λ5007 flux with the S/N ≃ 20–30 with a ≃1.5 hr integration. Then, it will require ≳50–100 hr of telescope time to detect [O iii]λ4363 with 3σ significance from individual objects with the same instrument setup. Therefore, the follow-up observation of [O iii]λ4363 for the EELG sample with Subaru/MOIRCS appears to be impractically expensive, and such observations can be made more efficiently with telescopes with larger apertures on the ground and those in space.

5.4.2. ISM Ionization Properties

Figure 13 shows O₃₂ as functions of R₂₃ (left panel) and stellar mass (right panel). O₃₂ is an indicator of the ionization parameter (e.g., McGaugh 1991; Kewley & Dopita 2002), and R₂₃ is a strong-line gas-phase metallicity indicator (e.g., Pagel et al. 1979; Kobulnicky & Phillips 2003). Compared to local SFGs from SDSS DR7 (Oh et al. 2011), our EELG sample at z ≃ 3.3 occupies the upper part of the plot, which means that they are more metal-poor and have higher ionization parameters. The solid line in the left panel of Figure 13 is the locally defined relation of the two line ratios presented by Maiolino et al. (2008) as a function of gas-phase oxygen abundance, $12+\mathrm{log}({\rm{O}}/{\rm{H}})$. Small squares on the line of the local relationship indicate $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ = 7.5, 8.0, 8.5, and, 9.0. While the majority of local galaxies have $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ ≃ 8.5–9.0, our EELGs appear to show R₂₃ consistent with $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ = 7.5–8.5. EELGs at z ≃ 3.3 also show higher ionization parameters compared to those of local SFGs. The observed increase of ≃1 dex in O₃₂ from z = 0 to z ≃ 3.3 can be translated to an increase of ≳1.5 dex in the ionization parameter (e.g., Nakajima & Ouchi 2014).

Zoom In Zoom Out Reset image size

In the figure, we also show normal SFGs at z ≃ 3.3 (Onodera et al. 2016) and LAEs at z ≃ 3.1 (Nakajima et al. 2020). As shown in Nakajima & Ouchi (2014) and seen in the figure, LAEs at z = 2–3 tend to show higher O₃₂ indices than normal SFGs at z ≃ 2–3. The O₃₂ of our entire EELG sample at z ≃ 3.3 appears to be distributed between the two populations.

It is also noticeable that EELGs with higher O₃₂ values are more likely to show larger EW([O iii]λ5007). As shown in Figure 12, there is an anticorrelation between EW([O iii]λ5007) and galaxy stellar mass. Massive EELGs at z ≃ 3.3 with $\mathrm{log}{M}_{\star }/{M}_{\odot }\gtrsim 9.5$ typically have EW([O iii]λ5007) ≲ 300 Å, while lower-mass ones show EW([O iii]λ5007) ≳ 500 Å. From the right panel of Figure 13, it is clearer that at z ≳ 3 the distribution of the two line ratios considered here for massive EELGs is similar to that of normal SFGs and that for low-mass ones it is closer to that of LAEs (see also Suzuki et al. 2017). These trends are also seen for the composite measurements.

5.4.3. Ionizing Photon Production Efficiency

In the top left panel of Figure 14, we compare ξ_ion,0 as a function of β_UV. Our EELG sample at z ≃ 3.3 shows a flat relation between the two parameters with a median value of $\mathrm{log}{\xi }_{\mathrm{ion},0}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})=25.54$, though 8 out of 19 objects have upper limits for ξ_ion,0. This is also the case for the composite spectra. Emami et al. (2020) have found a similar trend, i.e., no β_UV dependence of ξ_ion,0 for a sample of low-mass (M_⋆ = 10^7.8–10^9.8 M_⊙) galaxies at z ≃ 2 (see also Matthee et al. 2017, for a sample of Hα emitters). They corrected dust attenuation by using the SMC curve for the UV luminosity from SED fitting and the Milky Way (MW) extinction curve (Cardelli et al. 1989) for the Hα luminosity based on the Balmer decrement. On the other hand, different trends have also been reported in the literature. Bouwens et al. (2016) reported a ≃0.3 dex increase of mean ξ_ion,0 for the bluest galaxies (β_UV < −2.3) for a sample of sub-L* SFGs at z = 4–5, while a constant value of $\mathrm{log}{\xi }_{\mathrm{ion},0}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})=25.34$ was measured at β_UV > −2.3. Their conclusion was not affected by the choice of the dust attenuation curve (Calzetti or SMC). Shivaei et al. (2018) found a similar trend to Bouwens et al. (2016) for a sample of SFGs at z ∼ 2 from the MOSDEF survey (Kriek et al. 2015) using the Calzetti curve for dust correction. When using the SMC curve, their relationship becomes almost flat with a minimum of ξ_ion,0 at β_UV ≃ −1.5 (Shivaei et al. 2018).

Zoom In Zoom Out Reset image size

To check whether our result changes by using the Calzetti curve instead of the SMC curve, we derived ξ_ion,0 assuming the Calzetti curve for stellar continuum following the conversion from β_UV to E(B − V)_star presented in Reddy et al. (2018a) and the MW extinction curve for nebular emission lines with a conversion E(B − V)_neb,MW = 1.34 E(B − V)_star,Calz derived by Theios et al. (2019). Similar to Shivaei et al. (2018), ξ_ion,0 becomes 0.1–0.4 dex smaller with the Calzetti curve for dust correction than with the SMC curve, with a median difference of 0.28 dex. The difference appears to become larger for higher-β_UV objects as also reported by Shivaei et al. (2018), but there still seems to be no correlation between ξ_ion,0 and β_UV, likely to be partly due to the small number statistics with a large number of upper limits and a relatively large scatter of the distribution of the sample. In addition to the choice of the attenuation curves for dust correction, the conversion factor between the nebular and stellar E(B − V) could be a function of galaxy properties such as stellar mass, SFR, and hence sSFR as seen in the local universe (Koyama et al. 2019). Although it is not clear whether such a relation is applicable to galaxies at z ≃ 3.3 with intense [O iii] emission lines, this would potentially introduce an additional uncertainty in ξ_ion,0.

The main reason why our EELGs at z ≃ 3.3 show no dependence of ξ_ion,0 on β_UV could be their selection as the most intense [O iii] + Hβ emitting objects. Unlike our selection method, other studies mentioned above did not impose any cuts in emission-line strengths. In fact, objects with $\mathrm{log}{\xi }_{\mathrm{ion},0}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})\simeq 26$ exist at all β_UV up to β_UV ≃ −1 in the tail of the distribution in Shivaei et al. (2018).

The ξ_ion,0 of our EELG sample and those of SFGs at z ≃ 2 (Emami et al. 2020), LAEs at z ≃ 3.1 (Nakajima et al. 2020), SFGs at z ≃ 4–5 (Bouwens et al. 2016), and LAEs at z = 4.9 (Harikane et al. 2018) are shown as a function of UV absolute magnitude M_UV in the top right panel of Figure 14. LAEs at z = 4.9 from Harikane et al. (2018) follow the trend of the ξ_ion,0–M_UV relation of SFGs at z ≃ 4–5 (Bouwens et al. 2016; see also Lam et al. 2019). Our EELGs at z ≃ 3.3 show relatively bright M_UV, similar to brighter objects of Bouwens et al. (2016), but have on average higher ξ_ion,0 by up to ≃0.8 dex. This may be also an indication of the extreme nature of our sample based on the selection with significant K-band excesses.

Comparing to LAEs at z ≃ 3.1 from Nakajima et al. (2020), our EELGs show a similar range of ξ_ion,0 but brighter UV luminosities by 2–3 mag despite the similar stellar mass range. A part of the large difference in M_UV could be due to the assumption of a uniform E(B − V)_star = 0.01 based on the SED fitting by Fletcher et al. (2019). Note, however, that even a small increase of the dust attenuation parameter can result in a significant increase of the UV luminosity especially when using the SMC extinction curve. If we assume that our EELGs are extinction free, M_UV becomes fainter on average by ≃1.1 mag. Therefore, the difference in dust correction methods cannot fully explain the difference in UV luminosities between two samples and the reason why the LAE selection can pick up fainter objects with similar stellar masses. Another possible reason of the difference is that Fletcher et al. (2019) performed the SED fitting by using BPASS v2.1 stellar population models (Eldridge et al. 2017), while we adopted BC03 models in the CIGALE framework. Eldridge et al. (2017) pointed out that BC03 models are redder in UV at very young ages of <10 Myr and also redder in red−optical to near-IR wavelength at ages more than a few times 10² Myr. The choice of different population synthesis models can result in a systematic effect in estimating stellar masses. SFGs at z ≃ 2 by Emami et al. (2020) include objects with much fainter UV luminosities, M_UV ≳ −18, and they found no correlation between ξ_ion,0 and M_UV in contrast to Nakajima et al. (2020). Note that to derive the UV luminosities Emami et al. (2020) used BC03 templates and an SMC curve for dust correction that are different from those used in Nakajima et al. (2020), and that the Emami et al. sample shows on average smaller EW([O iii]λ5007) of ≃100 Å than that of Nakajima et al. (2020) of 600–1000 Å, which may be reasons for the discrepancy at the faintest UV magnitude bin.

A strong correlation between ξ_ion,0 and EW([O iii]) has been claimed for high-redshift EELGs (Emami et al. 2020; Tang et al. 2019; Nakajima et al. 2020), as well as nearby SFGs (Chevallard et al. 2018). The bottom left panel of Figure 14 shows the relation between ξ_ion,0 and EW([O iii]λ5007) for our EELG sample at z ≃ 3.3, along with that for LAEs at z ≃ 3.1 (Nakajima et al. 2020) and the best-fit relations derived for EELGs at 1.3 < z < 2.4 (Tang et al. 2019) and for local EELGs (Chevallard et al. 2018). Our EELG sample at z ≃ 3.3 follows a similar trend to the relationships in the literature. The correlation between ξ_ion,0 and EW([O iii]) appears to be the strongest among four parameters shown in Figure 14. On the other hand, the composite data do not show such correlation. This is probably because we adopted the inverse-variance weighting, which results in putting more weights on objects with the strongest, i.e., the highest-S/N, emission-line fluxes. As shown in the figure, our EELGs with EW([O iii]λ5007) ≳ 300 Å are exclusively low-mass galaxies with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$, which can also be clearly seen in the composite spectra. Moreover, those with EW([O iii] λ5007) ≳ 1000 Å show particularly high ξ_ion,0 of ≳25.5, suggesting that those with low stellar mass and high EW([O iii]) are very efficient in producing LyC photons.

In the bottom right panel of Figure 14, we show the relationship between ξ_ion,0 and O₃₂. EELGs at z ≃ 3.3 do not show an apparent correlation, unlike LAEs at z ≃ 3.1 (Nakajima et al. 2020) as shown with small red pentagons. Quantitatively, the Spearman's rank correlation coefficients are 0.05 and 0.48, and the corresponding two-sided p-values are 0.84 and 0.001, for our EELGs at z ≃ 3.3 and LAEs at z ≃ 3.1, respectively. Shivaei et al. (2018) also found a 0.3 dex enhancement of ξ_ion at large $\mathrm{log}{O}_{32}\gtrsim 0.5$, while ξ_ion is constant at $\mathrm{log}{O}_{32}\lesssim 0.5$. They also found that this enhancement can be observed also at low $\mathrm{log}{O}_{32}$ galaxies when the SMC extinction curve is used for dust correction instead of the Calzetti curve. The systematic effect due to dust correction could be one of the reasons why there is no correlation between ξ_ion,0 and O₃₂ for our EELG sample. Another reason would be again the biased selection of EELGs as the strongest [O iii] emitters. Compared to LAEs at z ≃ 3.1, our EELGs show similar ξ_ion,0 but systematically lower O₃₂ values by ∼0.2 dex, though the relatively large number of objects with only limits hampers a detailed comparison of the distribution of two populations. Because O₃₂ can be used as indicators of ionization parameters and gas-phase metallicity (e.g., Maiolino et al. 2008), it is expected that objects with larger O₃₂ are efficient hydrogen-ionizing photon producers and show elevated ξ_ion values. As seen before, most of the EELGs with larger $\mathrm{log}{O}_{32}$ of ≳0.4 are low-mass galaxies with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$ that follow the distribution of LAEs at a similar redshift well.

In summary, our EELG sample at z ≃ 3.3 contains efficient hydrogen-ionizing photon producers, in particular at lower stellar masses of $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$. The relationships between ξ_ion,0 and various observed properties (Figure 14) suggest that our EELGs have similar properties to luminous SFGs at z ≳ 4 (Bouwens et al. 2016) and EELGs at 1.4 ≲ z ≲ 2.5 (Tang et al. 2019). Comparison with LAEs at z ≃ 3.1 indicates that LAEs are systematically less luminous in UV and more highly ionized objects than the EELGs at z ≃ 3.3. Among the investigated parameters, ξ_ion,0 of our EELGs shows a marginal correlation with EW([O iii]), while other parameters do not have an impact on ξ_ion,0. This would be mainly due to our biased selection of objects with large EW([O iii] + Hβ). Another possible source of systematics is the choice of the attenuation curve for dust correction. If we use the Calzetti curve instead of the SMC curve, ξ_ion,0 are reduced by 0.1–0.4 dex. Although the change is more prominent in objects with larger β_UV, the conclusions of the comparison above are still unchanged. The fact that many of our EELG objects do not have detected Hβ and [O ii]λ3727 also dilutes correlations if any. A deeper spectroscopic follow-up observation (e.g., Nakajima et al. 2020) is required to make a more detailed comparison.

Our primary motivation to search for EELGs with intense [O iii] emission at z > 3 is to investigate a population of galaxies that can be considered as analogs that contributed to reionizing the universe in the EoR as close as possible to the epoch, but not too high redshift, where the universe is fully opaque to the LyC emission. Indeed, we have shown that we were able to select such galaxies at z ≃ 3.3 with EW([O iii]λ5007), reaching up to ≃2000 Å efficiently by the Ks-band excess method. Such large EW(Hβ + [O iii]) values are comparable to those inferred for SFGs at the EoR (e.g., Smit et al. 2014, 2015; Roberts-Borsani et al. 2016; Stark et al. 2017; De Barros et al. 2019; Endsley et al. 2020).

In addition to such large EW([O iii]), we have shown that the EELG sample at z ≃ 3.3 in this study is characterized by enhanced SFR (and sSFR) compared to the MS, large O₃₂ values, low level of dust attenuation, bright in UV, and high ξ_ion,0, suggesting that they are young, low-metallicity galaxies in a bursty phase of star formation with highly ionized ISM. These physical properties are more pronounced in the less massive ones, similar in many aspects to LAEs at a similar redshift and blue SFGs at even higher redshifts, as well as EELGs at lower redshifts. Some of these galaxies in the literature are found to show various degrees of LyC escape. Here we compare observational properties that we have in hand in this study and are considered as indicators of the LyC leakage.

Figure 15 shows the relation between O₃₂ and EW([O iii]λ5007) for our EELG sample at z ≃ 3.3 and confirmed LyC leakers at z ≃ 0.3 (Izotov et al. 2016, 2018b) and Ion2 at z = 3.2 (Vanzella et al. 2016, 2020). We also overplot LAEs at z = 3.1 presented by Nakajima et al. (2020), including those without LyC detection. LyC leakers are color-coded by the LyC escape fraction (f_esc). It is notable that there is a large diversity in f_esc among the LyC leakers from a few percent to >70%, and they tend to have large O₃₂ values of ≳3–5.

Zoom In Zoom Out Reset image size

Compared to LyC leakers shown in Figure 15, our EELGs follow a similar O₃₂–EW([O iii]λ5007) relationship to confirmed LyC leakers, though many of our EELGs have only lower limits in O₃₂ (i.e., nondetection in [O ii]λ3727). Our EELGs at z ≃ 3.3 show large O₃₂ values of ≳3 for about a half of the sample, and most of them are low-mass galaxies with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$ as also seen in the composite measurements. Also, as shown before, these low-mass EELGs show large EW([O iii]λ5007) ≳ 500 Å.

The comparison makes it tempting to classify these EELGs at z ≃ 3.3 as LyC leakers with large values of f_esc. Indeed, photoionization models suggest that O₃₂ can be used as an indicator of f_esc (e.g., Jaskot & Oey 2013; Nakajima & Ouchi 2014), and observational data also appear to follow the predictions (Faisst et al. 2016, and references therein). However, recent studies suggest that having large O₃₂ values is a necessary condition, not a sufficient condition, possibly due to geometrical effects along the line of sight (Izotov et al. 2018b; Naidu et al. 2018; Jaskot et al. 2019; Nakajima et al. 2020). Bassett et al. (2019) have made a detailed comparison of the relation between f_esc and O₃₂ with density-bounded models using the MAPPINGS V photoionization code (Allen et al. 2008; Sutherland et al. 2018). They found that O₃₂ shows the largest variation at f_esc < 0.1, while it seems to provide little constraint at larger f_esc than ≃0.1, especially for low-metallicity cases. The fact that many of the LAEs without LyC detections have comparably large O₃₂ to LyC leakers also supports these scenarios (Nakajima et al. 2020). Yet EELGs with very large O₃₂ of ≳10 are highly likely to contain LyC leakers, because such high O₃₂ ratios can emerge from a density-bounded nebula (i.e., fully ionized gas cloud) from which LyC emission can escape easily compared to an ionization-bounded nebula where the size of the nebula is determined by the Strömgren sphere (Nakajima & Ouchi 2014; Jaskot et al. 2019).

Additionally, Bassett et al. (2019) have also discussed that gas-rich galaxy mergers may play a role in breaking the O₃₂–f_esc relation by redistributing gas in the system. Such mergers of low-mass, metal-poor galaxies or inflow of metal-poor gas into the system are suggested to trigger bursty star formation by a sudden change in the gas accretion rate, which could result in breaking the equilibrium of self-regulation of star formation (Lilly et al. 2013; Onodera et al. 2016; Amorín et al. 2017). Alternatively, other mechanisms that are not necessarily classified as mergers, but minor interactions and violent disk instabilities, can also produce star-forming clumps and cause mixing of gas in galaxies (Ribeiro & Le Fèvre 2017). Such clumpy structures are often observed in low-mass metal-poor EELGs at z < 1 (e.g., Kunth et al. 1988; Lagos et al. 2014; Amorín et al. 2015; Calabrò et al. 2017). Among 11 EELGs in our sample with $\mathrm{log}{M}_{\star }/{M}_{\odot }\lesssim 9$, 5 of them show indications of ongoing merger or clumpy structure, while the rest show compact, point-like morphology in HST/ACS F814W images shown in Figure 16. The apparent diversity of morphology of our EELG sample at z ≃ 3.3 could partly reflect different merger stages. Note, however, that these merger- or clumpy-like structures can be due to chance projections along the line of sight. In order to fully characterize the nature of individual components of each object, deep high-resolution multiband imaging and integral field spectroscopy using currently available facilities, such as HST and ground-based 8–10 m class telescopes with adaptive optics, and future ones, such as the James Webb Space Telescope and ground-based 30 m class telescopes, will be required.

Zoom In Zoom Out Reset image size

Whether it is possible for SFGs in the EoR to provide sufficient hydrogen-ionizing photons (i.e., LyC photons) escaped from them is crucial to assess the dominant ionizing population of the reionization. Important parameters here are the hydrogen-ionizing photon production efficiency in the case of no LyC escape ξ_ion,0 and LyC escape fraction f_esc, as the product of the two parameters determines the amount of LyC photons outside of the system. As a canonical value, $\mathrm{log}{\xi }_{\mathrm{ion},0}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})$ = 25.2–25.3 has been widely used (e.g., Robertson et al. 2013). The value is suggested from BC03 stellar population models at about solar metallicity with a young stellar population with ages of ≃1–10 Myr for a single-burst SFH. It is also observationally supported by the analysis of β_UV slope of galaxies at z ≃ 7–9 (Dunlop et al. 2013; Robertson et al. 2013).

As seen in Figure 14 and Table 4, the majority of our EELGs at z ≃ 3.3 have $\mathrm{log}{\xi }_{\mathrm{ion},0}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})$ larger than the canonical value, though they include those with upper limits. This means that the majority of our EELGs are more efficiently producing LyC photons than expected by the BC03 models with the aforementioned parameters. Therefore, in order to explain the elevated values of ξ_ion,0, different properties of massive stars in galaxies would be needed, such as the IMF, SFH, inclusion of binary populations, and metallicity. Shivaei et al. (2018) have made a detailed investigation of the effect of these parameters using the BPASS and BC03 models. Their fiducial model consists of the Salpeter IMF (Salpeter 1955), with a lower and upper mass cut of 0.5 and 100 M_⊙, respectively, and a slope α = −2.35, constant SFH over 300 Myr, Z = 0.15 Z_⊙, and without binary populations, resulting in $\mathrm{log}{\xi }_{\mathrm{ion}}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})=25.23$. In their analysis, BC03 models with the same parameterization as the fiducial model give $\mathrm{log}{\xi }_{\mathrm{ion}}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})=25.06$, ∼0.1 dex lower than the canonical value mentioned above. Including the binary population, increasing the upper mass cut of IMF to 300 M_⊙, increasing the IMF slope to α = −2.00, and lowering the metallicity to 0.07 Z_⊙ result in an increase of the ξ_ion by ∼0.17, ∼0.16, ∼0.18, and 0.02 dex, respectively (see Table 1 in Shivaei et al. 2018). Since our EELGs at z ≃ 3.3 typically have $\mathrm{log}{\xi }_{\mathrm{ion},0}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})\simeq 25.5$ and as large as ≃26.0, a combination of more than one change from the fiducial model would be required to explain these values.

In addition to the effects related to their stellar populations as discussed above, the nonzero LyC escape fraction can scale the ξ_ion,0 by 1/(1 − f_esc). Statistically speaking, f_esc = 10%–20% is supposed to be necessary for galaxies to fully ionize the universe at z > 6 (e.g., Ouchi et al. 2009; Robertson et al. 2013), and in this case the intrinsic ξ_ion can be even larger by ≃0.05–0.1 dex. If we simply use the best-fit relation between O₃₂ and f_esc derived by Faisst (2016), f_esc = 0.05, 0.1, and 0.45 correspond to O₃₂ = 3, 5, and 10, respectively, which can be translated to the increase of the intrinsic ξ_ion of 0.02, 0.05, and 0.26 dex, respectively. Therefore, even a moderate escape of LyC from our EELGs can further increase the discrepancy. For those with large values of O₃₂ ≃ 10 and ξ_ion,0, the intrinsic ξ_ion can be as large as $\mathrm{log}{\xi }_{\mathrm{ion}}/({\mathrm{erg}}^{-1}\,\mathrm{Hz})\gtrsim 26$, which cannot be explained even by the most extreme set of parameters for BPASS models considered by Shivaei et al. (2018). As suggested by Figures 11 and 13, an even younger stellar population age with a burst-like SFH, lower metallicity, or a combination of them would be required to explain the most extreme objects in our sample. Note that, as discussed in the previous section, the LyC escape does not seem to be determined solely by the ionization properties inferred by, for example, O₃₂ and EW([O iii]), but the geometry of ionized gas is suggested to be crucial (e.g., Nakajima et al. 2020).