© The Authors 2024

1. Introduction

Over the past decades, several studies have discovered a population of galaxies that undergo very strong star formation events. They are called extreme emission line galaxies (EELGs) and are typically characterized by high equivalent widths (EWs) of hundreds to thousands of Angstroms which are driven by elevated specific star formation rates (sSFRs) up to 10−100 Gyr⁻¹ in galaxies with stellar masses ≲10⁹ M_⊙, with subsolar metallicities and little dust (e.g., van der Wel et al. 2011; Maseda et al. 2014; Amorín et al. 2015; Forrest et al. 2017). These large emission line EWs are most commonly seen in the [OIII]λλ4959,5007 (hereafter [OIII]) and Hα lines in the rest-frame optical. The nebular emission lines are often driven by the ionizing photons produced in massive and short-lived O and B stars – or active galactic nuclei (AGN)–, whereas the surrounding rest-optical continuum includes the contribution from longer-lived and less massive stars (Eldridge & Stanway 2022, for a review).

Such extreme systems are rare in the local Universe, but the number density of EELGs increases with increasing redshift (e.g., Smit et al. 2014) and they are expected to be common in the Epoch of Reionization (EoR). Recent results indicate that the population of galaxies in the EoR shows a mean EW([OIII]+Hβ) = 780 Å (e.g., Endsley et al. 2023a). Examples of EELGs in the local universe are the so-called HII galaxies (Terlevich et al. 1991) or blueberry galaxies (Yang et al. 2017), and Green Pea galaxies (Cardamone et al. 2009; Amorín et al. 2010, 2012; Izotov et al. 2011). Besides allowing us to study extreme conditions in galaxies in the early universe, such systems are also essential to understanding the process of reionization. Such starburst galaxies are currently presumed to be the main drivers of hydrogen reionization over 5.5 < z < 15 (e.g., Robertson et al. 2015; Finkelstein et al. 2019; Yung et al. 2020a,b), which serves to motivate many of the searches for high-z EELGs as well. To understand the properties of EELGs, large and representative samples of EELGs must be assembled.

The study of EELGs in the local Universe up to intermediate redshifts has been a significant field in recent years because they resemble the properties of galaxies in the EoR with high EWs of UV and optical lines (Smit et al. 2015; Hutchison et al. 2019). With the advent of the JWST era (Gardner et al. 2006, 2023), the discoveries of young galaxies are rapidly moving deep into the EoR, where such systems can be studied now directly up to z ∼ 9 with NIRCam and NIRSpec (e.g., Withers et al. 2023; Simmonds et al. 2024; Endsley et al. 2023b; Davis et al. 2023). The first step is to find and characterize such systems to conduct detailed spectroscopic studies.

Many of the searches for EELGs select samples using either narrow- (e.g., Sobral et al. 2013; Iglesias-Páramo et al. 2022; Lumbreras-Calle et al. 2022), broad- (e.g., van der Wel et al. 2011; Onodera et al. 2020; Kojima et al. 2020; Chen et al. 2024; Davis et al. 2023) or medium-band photometry (e.g., Cohn et al. 2018; Withers et al. 2023; Simmonds et al. 2024), and slitless spectroscopy (e.g., Maseda et al. 2018; Kashino et al. 2023). Using photometry, EELGs are selected based on the enhancing effect the emission line fluxes have in some specific filter, depending on its redshift. Unlike the commonly used Lyman break technique, broad-band color selection does not require high S/N in the rest-frame UV-optical continuum. Faint continuum EELGs are difficult to detect through their continuum emission alone but can be selected via their color excess due to the presence of emission lines. This population of very faint-continuum galaxies may play an essential role in reionization if found in great numbers (e.g., Endsley et al. 2023a), yet will not be detected with many of the commonly used selection criteria.

Deep imaging with new JWST capabilities opens a new window to identify such young metal-poor starbursts. Our sample originates from the Cosmic Evolution Early Release Science survey (CEERS, Finkelstein et al. 2023, 2024) which is a JWST Early Release Science program that obtained imaging and spectroscopy of the Extended Groth Strip (EGS, 14h19m00s +52°48′00″, Davis et al. 2007) field with three instruments and five coordinated parallel observing modes. We use the CEERS NIRCam and NIRSpec data to look at EELGs over a redshift range 4 ≲ z < 9 to study their physical properties in detail to determine what physical conditions are necessary for a galaxy to be an EELG. We focus on the [OIII]λ5007 line, but we also consider cases where Hα may also have extreme EWs.

This paper is organized as follows: in Sect. 2, we present the observational data and the templates of EELGs we consider in this work. Particularly in Sect. 2.5, we present our NIRCam synthetic observations to estimate the color excess to identify EELG candidates. In Sect. 2.6, we present the color-color diagrams used and the final sample of EELG candidates. Regarding the analysis of the data, in Sect. 3.1 we present the photometric fluxes and EWs of emission lines. In Sect. 3.2 we present the physical properties of our sample based on SED fitting and in Sect. 3.3 we present the physical sizes and the star formation rate (SFR) surface density of the galaxies in our sample. In Sect. 3.4 we present the NIRSpec spectroscopy of a subsample of 47 galaxies in our candidate list. In Sect. 4, we present our results and discussions on the ionization properties of the EELG candidates and the comparison with the control sample of non-EELG galaxies. Finally, in Sect. 5, we present our conclusions.

Throughout this paper, we adopt a Λ-dominated flat universe with Ω_Λ = 0.7, Ω_M = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹. All magnitudes are quoted in the AB system. Equivalent widths are quoted in the rest frame and are positive for emission lines. We consider log(O/H)_⊙ = 8.69 (Asplund et al. 2009).

2. Data and sample selection

2.1. The CEERS survey

The complete CEERS program involves imaging with the NIRCam short and long-wavelength channels in ten pointings, observed as coordinated parallels to primary observations with the NIRSpec and the MidInfrared Instrument (MIRI; Wright et al. 2015). In this paper, we use photometric data from NIRCam and NIRSpec spectra. The CEERS observations were split into two epochs. The first epoch was executed in June 2022, and the second in December 2022. An additional third epoch of NIRSpec prism spectra (pointings 11 and 12) was taken in February 2023, due to the MSA short that affected prism 9 + 10 in December 2022. In total, ten pointings were taken with NIRCam imaging, including seven filters per pointing (F115W, F150W, F200W, F277W, F356W, F410M, and F444W, see Fig. 1), reaching 5σ depths of ∼28.5–29.2 AB mag for point sources for a total area of ∼97 arcmin².

Fig. 1. CEERS NIRCam filters. The black lines are the estimated observed wavelength at a given redshift for [OII]λ3727 (dotted-dashed), [OIII]λ5007 (dashed) and Hα (solid). The filter response functions for NIRCam filters are arbitrarily scaled.

2.2. Photometric catalogs

We use version v0.51.2 of the CEERS Photometric Catalogs (Finkelstein et al., in prep.). The catalog contains 101808 sources. The NIRCam images used are publicly available, and we refer the reader to Bagley et al. (2023) for a complete description of the data reduction. We also make use of the HST data from CANDELS (Grogin et al. 2011; Koekemoer et al. 2011). For the CEERS pointings 1, 2, 3, and 6, the images are available in the Data Release 0.5¹, while for the pointings 4, 5, 7, 8, 9, and 10, the images are available in the Data Release 0.6².

A complete description of the photometric catalog will be presented in Finkelstein et al. (in prep.). Briefly, the photometry was performed with SExtractor (v2.25.0; Bertin & Arnouts 1996) with F277W and F356W as the detection image. The fiducial fluxes were measured in small Kron apertures corrected by large-scale flux, following the methodology in Finkelstein et al. (2023).

Photometric redshifts are estimated with EAZY (Brammer et al. 2008), following the methodology in Finkelstein et al. (2023). The code fits nonnegative linear combinations of user-supplied templates to derive probability distribution functions for the redshift based on the fit quality to the observed photometry for each detected source. The redshift is allowed to vary between 0 − 20. For this catalog, new templates from Larson et al. (2023a) were included to improve the goodness of the fit in high-z sources. This new set of templates includes models with low metallicity (5% solar), young stellar populations (10⁶, 10^6.5, 10⁷ yr), and high ionization parameters (log U = −2). In this paper, we use the best-fitting photometric redshift defined as za in EAZY.

In this paper, we use CEERS NIRCam photometry and synthetic observations of templates of local EELGs to find such systems at 4 ≲ z < 9. In this redshift range, bright emission lines such as [OIII]λ5007 and Hα may fall in the filters F277W, F356W, F410M, and F444W (see Fig. 1). We note that at 6.6 ≲ z < 9, Hα redshifts beyond the bandpass of F444W, and then in this redshift range the selection is based only on [OIII]+Hβ emission. The templates of local EELGs are described in the following Sect. 2.4.

2.3. Spectroscopic data

The CEERS survey also includes six NIRSpec pointings, numbered p4, p5, p7, p8, p9, and p10. Each of these pointings has observations with the three NIRSpec medium resolution (G140M, G235M, and G395M) gratings and with the low-resolution Prism. Two more fields, p11 and p12, were observed with the prism in February 2023 because the prism observations p9 and p10 were severely impacted by a short circuit. The grating set covers from 0.97−5.10 μm with a resolving power R = λ/Δλ of ∼1000, while the prism covers from 0.60−5.30 μm with a resolving power of 30 < R < 300, depending on the wavelength. More objects can be observed simultaneously with the prism thanks to the shorter length (in pixels) of the prism spectra.

We adopt the NIRSpec data produced by the CEERS collaboration using the STScI JWST Calibration Pipeline³ (Bushouse et al. 2022). Specifically, we use the JWST pipeline to perform the standard reductions, including the removal of dark current and bias, flat-fielding, background, photometry, wavelength, and slitloss correction for each exposure. We also perform additional reductions to remove the 1/f noise and the snowballs. The 2D spectra of each target are then rectified and combined to generate the final 2D spectra. The details of the data reduction are presented in Arrabal Haro et al. (2023, and in prep.).

2.4. Empirical templates of EELGs

To create synthetic NIRCam observations of EELGs, we use the templates of the Automated Spectroscopic K-means-based (ASK) classes presented in Sánchez Almeida et al. (2010). According to this scheme, the ∼one-million SDSS-DR7 galaxies with an apparent magnitude brighter than 17.8 can be classified in only 28 ASK classes based exclusively on the features and shape of their rest-frame and normalized by the u magnitude of the optical spectrum.

The classification details are described in Sánchez Almeida et al. (2010) but very briefly; to classify a galaxy, the spectrum is considered as a multidimensional vector, and it is assumed that the vectors are clustered around several cluster centers, known as classes. Then, the algorithm works iteratively to assign each rest-frame spectrum to the nearest class and then the class template as the average over all the class members. In the end, the algorithm finds the number of classes and their templates and assigns to each galaxy spectrum one of the classes according to a certain probability. All the galaxies in a class have very similar spectra; their average spectrum is the template spectrum of the class. We note that each template has a rest-frame wavelength coverage from ∼3000 to ∼9000 Å.

These ASK classes are labeled according to their u − g color, from the reddest, ASK 0, to the bluest, ASK 27. Most (∼99%) galaxies in the SDSS-DR7 were classified into only 17 ASK major classes, with 11 additional minor classes including the remaining ∼1%. Sánchez Almeida et al. (2012) show that these rare classes correspond to metal-poor starbursts, and strong emission lines dominate their spectra. These minor classes, in particular, the classes ASK 15, ASK 17, ASK 20, ASK 21, ASK 25, ASK 26 and ASK 27 show EW([OIII]λ5007) and EW(Hα) > 170 Å in both lines (in Table 1 we list our measurements using LiMe⁴ which is a library that provides a set of tools to fit lines in astronomical spectra). The complete sample of local EELGs in the minor ASK classes is presented in Pérez-Montero et al. (2021) and a more detailed analysis of the sample and templates will be presented in a future paper (Amorín et al., in prep.). We also include in Table 1 the gas-phase metallicity associated with each template derived using the O32 calibration from Bian et al. (2018), which is the methodology that will be used in the following Section 4.3. We note that the templates are representative of metal-poor systems with metallicities ranging from 0.25−0.54 Z_⊙. We use the templates of these minor classes corresponding to EELG to estimate the expected colors of EELGs in the CEERS NIRCam photometry at z > 4.

2.5. Synthetic NIRCam observations

A simple and effective way of searching for EELGs using broad and medium band photometry is by using color selections, that is identifying the extreme colors caused by extreme emission lines producing excess up to ∼2 mag in neighboring filters (e.g., Withers et al. 2023).

Each spectrum of the 28 ASK classes was taken to a common wavelength range from 3030 to 9000 Å. We note that the ASK 15 has the shortest spectral range from 3747 to 8800 Å. In this case, we extrapolated the spectrum considering the continuum slope. Each spectrum was then redshifted to z = 2 − 10 (with Δz = 0.05 step sizes). Synthetic NIRCam observations for the seven filters observed in CEERS were created for each redshift to search for broad and medium band color excesses driven by strong [OIII]+Hβ and Hα emission. The mean density flux of each template in each filter for each redshift was calculated by $f_{\nu }\propto \frac{\int \lambda T{f}_{\lambda }\mathrm{d}\lambda }{\int \frac{T}{\lambda }\mathrm{d}\lambda }$, where T is the transmission curve of each filter, and f_λ is the normalized template spectrum. In this way, the color in AB mag is given by color₁₂ = −2.5[log(f_ν1)−log(f_ν2)]. Given that the templates have a limited spectral range, our synthetic NIRCam observations are limited only to the filters where the templates fall completely in the filters. Due to this, we only estimated the magnitudes in the filters where the template is within 10% of the filter transmission. This is illustrated in Fig. 2 for z = 6, where the magnitudes for ASK templates are estimated only in magenta filters. A video for the complete z range is available in this online movie.

Fig. 2. Illustration of the ASK templates at z ∼ 6 in the NIRCam filters. We illustrate which filters are considered to estimate the NIRCam color due to the limited wavelength coverage of the templates. In red are the NIRCam filters where the templates are outside the filter, while in purple are the filters where we estimate the NIRCam colors. Three ASK templates are illustrated as examples (ASK0, ASK10, and ASK15). A movie associated to this figure is available online.

In Fig. 3, we show the evolution with redshift of the colors used in this project, based on the synthetic NIRCam observations. With these results, we define color cuts based on the colors of the EELGs templates, which are separated from the other ASK classes in redshift windows. We note that depending on the redshift of an EELG and the filters it is observed with, its colors can drastically change by up to ±1 mag (Fig. 3), and the effect is stronger with increasing EWs of emission lines. In that sense, an EELG may appear as an extremely blue or red source as strong emission lines move in and out of different filter transmission curves (e.g., Papaderos & Östlin 2012; Papaderos et al. 2023).

Fig. 3. Redshift evolution of NIRCam colors of ASK templates. The major ASK, that is non-EELGs, classes are in gray-scale colors, according to the legend. The minor ASK classes, that is EELGs, are in colors according to the legend. The magenta and cyan-shaded regions are the z ranges and color cuts used to select EELGs.

2.6. Color-color diagrams

We define the color cuts based on the regions of redshift evolution of colors shown in Section 2.5. We use these color cuts to select candidates based on their position in color-color diagrams. While it is possible to identify EELGs by targeting emission from a single emission line complex, our strategy requires strong emission in both [OIII]+Hβ and Hα. In this paper, we focus on the most extreme EELG candidates, and because of that, we consider only the templates ASK 15, 17, 20, and 21, which show EW([OIII]+Hβ) > 680 Å.

We focus on two redshift ranges (see color-color diagrams in Fig. 4) which are selected due to the fact that EELGs can be distinguished from non-EELGs since bright emission lines enter or exit in a particular filter in the combination of colors. Additionally, we consider only sources with S/N > 3 in the photometric bands used in the color criteria. To remove objects with potential low-redshift solutions, we only considered sources whose ∫P(z > 3)dz > 0.5, ensuring that at least 50% of the posterior probability is above z = 3.

Fig. 4. Color diagrams to select EELGs at each redshift range. The color cuts are marked in magenta (top panel) and cyan (bottom panel) shaded regions in each panel according to Eqs. (1) and (2), respectively. The selected EELG candidates are in red circles, while parent samples at similar redshifts are in black circles. The sources with red crosses were discarded due to contaminated photometry as explained in the main text.

First, we consider sources whose photometric redshift is z = 3.78 − 5.3 (top panel in Fig. 4). In this redshift range, [OIII] falls in F277W, and Hα falls in F356W or F410M. The numeric thresholds for the color cuts were chosen to include the colored more extreme EELG template curves in the color-z space while excluding the grey curves for non-EELGs. We keep some regions where grey ordinary galaxies encroach upon the colored selection zones to allow a wider range in redshift since we are considering photometric redshifts. Based on our synthetic NIRCam color, we define the following color criteria:

$$\begin{array}{c}\hfill \{\begin{array}{c}F277W-F410M+\sigma (F277W-F410M)<0.35\\ F356W-F410M+\sigma (F356W-F410M)<-0.14\\ \vee \\ F356W-F410M-\sigma (F356W-F410M)>0.49\end{array}.\end{array}$$(1)

In the color criteria, σ refers to the observed color uncertainty. There are 4153 sources in this redshift range with S/N > 3 in the three filters. 634 sources (∼15%) satisfy the above conditions. An example of a candidate in this z range is shown in Fig. 5.

Fig. 5. Images and SED of CEERS1455, an example of an EELG candidate at zphot = 4.66. In the top panels, images (4″ × 4″) of the galaxy in filters (from left to right) RGB image, HST/F606W, HST/F814W, JWST/F277W, JWST/F356W, JWST/F410M, and JWST/F444W. The RGB image is made with Red = F410M, Green = F356W, and Blue = F277W. On the bottom panel, SED of the galaxy. The blue (black) squares are the NIRCam (HST) photometric points. NIRCam (HST) filters are in rainbow-scale (gray-scale) colors, according to legend. The filter transmission is plotted in arbitrary units depending on the observed flux density. The black solid line is the SED model and the inset text summarizes its physical properties (see description in Sect. 3.2). The vertical dashed lines represent the position of Lyα (black), [OIII] (blue), and Hα (red), according to their zphot.

We then consider sources whose photometric redshift is in the range z = 5.62 − 9 (bottom panel in Fig. 4). In this range, [OIII] may fall in F356W, F410M, or F444W and Hα in F410M or F444W. We note that at 6.6 ≲ z < 9, Hα redshifts beyond the bandpass of F444W, and then in this redshift range the selection is based only on [OIII]+Hβ emission. Based on our synthetic NIRCam color, we define the following color criteria:

$$\begin{array}{c}\hfill \{\begin{array}{c}F356W-F410M+\sigma (F356W-F410M)<0.1\\ F410M-F444W-\sigma (F410M-F444W)>0.28\\ \vee \\ F356W-F410M-\sigma (F356W-F410M)>0.6\\ F410M-F444W+\sigma (F410M-F444W)<-0.27\end{array}.\end{array}$$(2)

There are 2441 sources in this redshift range with S/N > 3 in the three filters. 420 sources (∼17%) satisfy the above conditions. An example of a candidate in this z range is shown in Fig. 6. We note that within this redshift range, for galaxies at z > 6.5, the condition ∫P(z > 6.5)dz > 0.5 is satisfied which implies at least 50% of the posterior probability is above z = 6.5.

Fig. 6. Images and SED of CEERS76140, an example of an EELG candidate at zphot = 7.48. Symbols as in Fig. 5.

With our method, we do not identify EELGs in the small interval between the two ranges, that is 5.3 < z < 5.62, because [OIII] starts to enter the F356W filter and we can not distinguish both populations using this filter. We also note that we are considering the photometric redshifts in the selection of the candidates. Based on recent results (e.g., Davis et al. 2023), the fraction of catastrophic redshifts (Δz > 3) can be as high as 55% along galaxies with signatures of extreme emission at similar redshifts. For this reason, we might be losing candidates in the selection if they have catastrophic redshifts. We restrict our selection using the photometric redshift to have the more secure candidates in the sample. We also note that a complete validation of the photometric redshift is in preparation within the collaboration by comparing them with the entire sample with spectroscopic redshifts. For now, some comparisons have been performed in reduced samples (e.g., Arrabal Haro et al. 2023; Chworowsky et al. 2023) finding good agreement between photometric and spectroscopic redshifts with differences lower than 0.5. In Sect. 3.4, we show a comparison with our sample of EELGs.

Our total sample of EELG candidates consists of 1054 galaxies. We highlight that by limiting our selection to galaxies with S/N > 3 in their photometry, we are limiting our sample to galaxies with < 29−30 mag in F200W and cover the magnitudes of their parent samples at similar redshifts and limited with the same S/N criteria in selecting filters (see Fig. 7). From these parent samples, we also define the control sample as the 1572 galaxies at the same redshift range and S/N that do not satisfy any selection color cuts in Eqs. (1) and (2). In Fig. 4, the control sample is selected from black small circles- galaxies that do not reside in the hashed areas. To clean the control sample from galaxies with contaminated photometry, we additionally require that the Kron radius is larger than 1.6 pixels to avoid including hot pixels as sources. We selected this threshold after visual inspection of galaxies with Kron radius < 1.6 pix which were mostly saturated pixels. We aim to compare the physical properties of this control sample with the selected EELG candidates.

Fig. 7. Distribution of the magnitude F200W for the sample of EELGs (unfilled histograms) and the parent samples at similar redshifts (filled histograms) in the two redshift ranges considered for the selection in blue and red, respectively. They are normalized to the total number of galaxies in each parent sample.

We clean the sample of EELGs by performing a visual inspection of the images of the candidates to remove galaxies whose photometry might be contaminated. We do not perform a clean process for possible AGN in the sample, we only clean for sources with contaminated photometry. We remove 46 sources with saturated pixels and 8 sources in the spikes of nearby stars. Our clean sample contains 1000 EELG candidates⁵. We flag (flag = 1) a total of 60 candidates that are close to the edge of the detector, close to bright sources, and sources with low surface brightness. We keep these sources, which represent ∼6% of the sample of EELGs, in the analysis. The final redshift distributions of the EELG candidates and the control sample are shown in Fig. 8. Each bin represents the percentage of the parent sample that is selected in the EELG and control sample. We note that there is a peak of the number of galaxies in the control sample at redshift z ∼ 7. In this redshift range, [OIII] starts to fall outside the F356W filter and inside the F410M and F444W filters. For this reason the colors tend to be zero and sources do not satisfy the selection criteria and most of the galaxies at this redshift range are in the control sample.

Fig. 8. Distribution of photometric redshift of the final sample of EELG candidates (in blue) and the control sample (in red). The gray-shaded region is the gap between the redshift range for selection. The normalization is made for the total of galaxies in the parent sample at similar redshift and S/N.

We will also perform a further cleaning of this sample of EELGs based on the photometric EWs described in the following Sect. 3.1. We validate our selection method by comparing it with available spectroscopy in the next Sect. 3.4.

The RGB (red = F410M, green = F356W, blue = F277W) images of the brightest F410M galaxies in our sample of EELGs are shown in Fig. 9. The colors in RGB images are associated with redshift ranges where bright emission lines fall within filters. For example, blueish galaxies are at z ∼ 3.78 − 4.82 where [OIII] falls in F277W and Hα in F356W. In this redshift range, in particular galaxies at z ∼ 3.8 may appear as greenish due to the low transmission of F277W compared with F356W for bright emission lines and/or Hα being brighter than [OIII]. Purpleish galaxies are at z ∼ 4.82 − 5.31 where [OIII] falls in F277W and Hα in F410M. Green galaxies are at z ∼ 5.62 − 6.63 where [OIII] falls in F356W and Hα in F444W. Redish galaxies are at z ∼ 6.63 − 7.68 where [OIII] falls in F410M. For higher redshift galaxies (z ∼ 7.68 − 9) [OIII] falls only in F444W but [OII] falls in F356W and then they show light green colors.

Fig. 9. RGB (F410M–F356W–F277W) images of the 100th brightest F410M EELG candidates in our sample. Each image has a size of 4″ × 4″.

We highlight that in the sample there are compact isolated sources (see for example CEERS29490, CEERS20040, CEERS64644 in Fig. 9) but also nearby sources that are likely to be interacting (see for example CEERS53557, CEERS14830, CEERS54272 in Fig. 9). Some galaxies in the sample look clumpy or chains of clumps (see for instance CEERS50196 or CEERS65253) dominated by one or more high surface brightness clumps, similar to EELGs at lower and intermediate redshifts (e.g., Amorín et al. 2015; Calabrò et al. 2017). There are cases with clear major companions that indicate they are currently experiencing a major merger (for instance CEERS65781 or CEERS62998 in Fig. 9).

3. Data analysis

3.1. Photometric fluxes

We estimate the fluxes of [OIII]+Hβ and Hα based on the observed photometry, assuming the flux excess with respect to the continuum is due to the emission line that falls in the filter. We estimate the continuum by assuming that it follows a power-law f_λ ∝ λ^β which matches the EELG templates. We consider all the NIRCam photometry observed in CEERS, excluding the bands where [OIII]+Hβ and Hα fall. We fit a linear model in log-log to estimate the β slope using LMFIT (Newville et al. 2016). We extrapolate the resulting model at the pivot wavelength of the filter to estimate the continuum of the line. The difference between the observed flux density and the estimated continuum is the flux of the line as F(line) = (f_λ − C_λ)Δλ, where Δλ is the filter width and C_λ is the continuum level.

Based on the line fluxes, we estimate the rest-frame EW = $\frac{F(\mathrm{line})}{C_{\lambda }(1+z_{\mathrm{phot}})}$. We note that these quantities might be overestimated since we are assuming all the flux excess is due to the bright emission line, but we are not considering other fainter lines that may contribute. For this reason, these quantities should be considered upper limits. We report our estimations in Table 2.

We perform this estimation for the sample of EELG candidates and the control sample. We note that in our sample of EELGs we find 264 galaxies with photometric EW([OIII]+Hβ) < 680 Å, which implies that our successful rate to detect EELGs with EW([OIII]+Hβ) > 680 Å using the selection criteria in Eqs. (1) and (2) is about 75%. In the cases where EW([OIII]+Hβ) < 680 Å, we note that the selection would improve if we added another condition in the selection, for example, the color F277W–F356W to probe the intense emission compared with the continuum for some redshift ranges. Or as already suggested in Davis et al. (2023), the inclusion of an additional observational medium band filter such as F300M would benefit to sample the continuum. In the following sections, we only consider the 736 EELG candidates with photometric EW([OIII]+Hβ) > 680 Å in the sample of EELGs.

The distributions can be found in Fig. 10. As expected, the sample of EELGs shows a higher fraction of galaxies with high fluxes of [OIII]+Hβ and Hα compared to the control sample. While the fraction of faint emitters (≲10⁻¹⁸ erg s⁻¹ cm⁻²) is larger in the control sample. We also note that the continuum, which is the projection of the linear fit, is similar in both samples and no differences are found in the ranges. Regarding the EWs, the EELG sample shows a larger fraction of galaxies with larger values compared to the control sample.

Fig. 10. Photometric fluxes of [OIII]+Hβ (left panels) and Hα (right panels) for the sample of EELG candidates (in blue) and the control sample (in red). From top to bottom: Integrated flux of the line in units of erg s−1 cm−2, continuum in the line wavelength in units of erg s−1 cm−2 Å−1, and rest-frame EW in units of Å.

We highlight that the galaxies in the control sample with high EW([OIII]+Hβ) > 680 Å show low EW(Hα) < 460 Å, which is below the threshold of EWs for our selection color cuts. For this reason, these intense [OIII]+Hβ emitters are not included in our sample of EELGs, since our selection criteria require that galaxies at z ≲ 6.75 (where Hα falls still within the F444W filter) are intense emitters of both [OIII]+Hβ and Hα above the thresholds defined by the color cuts. For this reason, there is an overlap between the sample of EELGs and the control sample regarding the high EWs.

3.2. Physical parameters

To obtain the physical properties of the sample of EELGs and the control sample, we perform SED fitting with only JWST photometry to have homogeneous data. We consider the entire set of filters used in the CEERS surveys, which are F115W, F150W, F200W, F277W, F356W, F410M, and F444W. We used BAGPIPES (Carnall et al. 2018) to estimate the physical parameters with the Bruzual & Charlot (2003) stellar population models. We fixed the redshift to the photometric redshift. We consider a delayed exponential τ-model for the SFH, where τ is the timescale of the decrease of the SFH. In the model, we consider the age ranging from 1 Myr to the age of the Universe at the observed redshift. We allow the τ parameter to vary between 0.1 to 10 Gyr. We allow the metallicity to vary up to 0.5 Z_⊙ freely. For the dust component, we consider the Calzetti et al. (2000) attenuation curve and let the A_V parameter vary between 0 − 2 mag. We also include a nebular component in the model, and we let the ionization parameter freely vary between −3 and −2. Some examples of the SED model are shown in Figs. 5 and 6. The distribution of the main obtained SED parameters is shown in Fig. 11, which includes stellar mass, SFR, sSFR, age, log U, stellar metallicity, absolute attenuation A_V, and τ.

Fig. 11. Distribution of the physical parameters based on the SED fitting for the sample of EELG candidates (in blue) and the control sample (in red). The vertical dotted (H−z sample) and dotted-dashed (M−z sample) lines are the mean values of each parameter.

We report the obtained parameters in Table 3. We find that our sample of EELGs shows stellar masses between 10^6.83 M_⊙ and 10^9.65 M_⊙ with a mean value of 10^7.84 M_⊙. They are actively forming stars with SFRs between 0.41 to 900 M_⊙ yr⁻¹, with a mean value of 13.64 M_⊙ yr⁻¹. This implies they show very high sSFR above 10^−8.47 yr⁻¹ and up to 10^−6.1 yr⁻¹, with a mean value of 10^−7.03 yr⁻¹. We also find that these galaxies are young, with times after the onset of star formation of roughly 1.7−792 Myr, with a mean value of 45 Myr. Regarding the ages, we also tested an exponential τ-model as SFH and we found a difference of 0.28 dex in the ages between both models which is a factor ∼ × 2 larger ages with the delayed exponential model.

Regarding the ISM properties of the EELGs, we find they show high ionization parameters with a mean value of log U = −2.17, with a saturation at log U = −2 which indicates harder ionizing spectra than those considered in the models are needed to model these sources. They show subsolar metallicities ranging from 0.008 to 0.44 Z_⊙ and a mean value of ∼0.25 Z_⊙. They also show low dust extinction, with A_V values ranging from ∼0 to 1.87 mag, but with a mean value as low as 0.34 mag. We note that we are considering a delayed τ-model for the SFH, but the results of the SED fitting indicate τ values of ∼5.1 Gyr on average, which indicates that they are consistent with a constant SFH.

As shown in Fig. 11, where we compare the sample of EELGs with the control sample, we note that both samples show similar ranges of stellar metallicity and absolute attenuation with similar mean values. Regarding other parameters, we find that the EELG candidates show lower stellar masses and higher SFRs and sSFRs compared to the control sample. Comparing the mean values, we find that the sSFRs of the EELGs are on average ∼0.8 dex higher than the mean value of the control sample. EELGs also show younger ages compared with the control sample which shows a mean value of 405 Myr. The control sample also shows a lower value of log U with a mean value of −2.48. To statistically verify these differences we performed a Kolmogórov–Smirnov (K-S) test and we found that for all parameters in Fig. 11 the two samples are significantly different, that is they show a p-value < 0.05. In the parameters where these differences are more significant, that is lowest p-values (∼10⁻²¹¹ − 10⁻¹¹³), are the sSFR, ages, and ionization parameter. While in the absolute attenuation, τ-parameter, and stellar metallicity these differences are less significant with higher p-values (∼10⁻²⁰ − 10⁻¹⁷) compared with the rest of the parameters.

We also split our samples of EELGs and control in two redshift bins: M−z for galaxies at z ≤ 6.5 and H−z for galaxies at z > 6.5 (see the vertical lines in Fig. 11 for their mean values) to check if there are differences between both populations. Regarding stellar mass, we find that both populations of EELGs show similar mean values. Regarding the SFR, we find that the M−z EELG galaxies show slightly lower (∼0.18 dex) SFRs than H−z EELGs. This implies that the M−z EELGs show slightly lower (∼0.2 dex) sSFRs than H−z EELGs. Similar trends are also observed in the control sample. Regarding the ages, the H−z EELGs show a mean value of 25 Myr while the M−z EELGs show a value of 51 Myr. Regarding the control sample, the M−z sample shows a mean age of 486 Myr while the H−z sample shows 235 Myr, which is expected due to the age of the universe. Even though we find differences in the mean values of these parameters when split by redshift, we note that the differences with the control sample are more significant than the differences with redshift.

In summary, EELGs are younger, more star-forming, and have higher ionization parameters compared to the control sample.

3.3. Physical sizes and SFR surface density

To measure the sizes of the EELG candidates and the control sample, we use the Galfit catalogs (v0.52) of the collaboration (McGrath et al., in prep.). Galfit (Peng et al. 2002, 2010) was run for sources with F356W < 28.5 mag using background-subtracted mosaics. Most of the EELG candidates are in the catalog (84% of the sample) independent of the quality of the fit. The photometric catalogs were used for making first guesses on source location, magnitude, size, position angle, and axis ratio. The Kron radius was used to determine an appropriate image thumbnail region for Galfit to fit. All galaxies within 3 mag of the primary source, but no fainter than 27 mag that fell within the thumbnail region were fit simultaneously. Galaxies that were not fit simultaneously were masked during the fitting process using the segmentation map. This procedure was performed for both F200W and F356W images. Here we only consider the F200W results since they are tracing mostly the rest-frame < 1750−4680 Å (depending on redshift) and therefore they represent well the size of the young stellar populations. In this way, we minimize the contribution from the ionized gas with bright emission lines. As a caveat, we note that the nebular continuum might be significant in EELGs and could contribute to the sizes that we use in this paper. We also note that the filter F200W may have a contribution of emission lines such as [OII]λ3727 for the galaxies at the lowest redshifts considered (3.8 < z < 5) which might lead to overestimations of their rest-UV sizes.

We only consider 360 sources with good fit quality (so-called Flag = 0). In Fig. 12, we show that the galaxies lie within 3σ of the size-mass relation at z ∼ 5 reported in Morishita et al. (2024) for rest-UV sizes. The subsample of EELGs shows effective radii ranging from 80 pc to 3 kpc, with a mean value of 0.57 kpc. Regarding the control sample, we find that 74% of them are in the Galfit catalog, but only 526 of them have good fit quality. This subsample is also displayed in Fig. 12. They show a slightly higher mean effective radii of 0.94 kpc.

Fig. 12. Size-mass relation for the EELG (in blue) and the control (in red) samples. The grey solid line is the size-mass relation at z = 5 from Morishita et al. (2024) and the shaded region the 3σ-scatter.

Based on the F200W sizes and the SFRs from the SED model, we estimate the projected SFR surface density as ${\mathrm{\Sigma }}_{\mathrm{SFR}}=\frac{{\mathrm{SFR}}_{\mathrm{SED}}}{2\pi r_{\mathrm{eff}}^2}$, where r_eff is the effective radius as the size in the F200W image. We find that the galaxies in our sample of EELGs show a wide range of Σ_SFR values ranging from 0.14 to 930 M_⊙ yr⁻¹ kpc⁻², with a mean value of 28.8 M_⊙ yr⁻¹ kpc⁻². Regarding the control sample, it shows a comparable range of values but a lower mean value of 7.6 M_⊙ yr⁻¹ kpc⁻².

3.4. Comparison with NIRSpec spectra

A total of 39 EELG candidates have NIRSpec spectra (see Table A.1). Two galaxies have spectra but in one of them (CEERS48859), the [OIII] would fall in the detector gap according to the z_phot, and no line measurements are performed in this galaxy. In the other galaxy (CEERS26436) [OIII] also falls in a detector gap but [OII]λλ3727,3729 (hereafter [OII]) and Hα are measured in this galaxy. From the sample of EELG candidates with NIRSpec spectra, 32 of them have prism spectra while 15 of them have medium-resolution grating spectra. Here, we use the data products of the collaboration DR as described in Sect. 2.3. We first estimate the spectroscopic redshift with [OIII] using LiMe (Fernández et al. 2023, 2024, for the technical details). We also measure the fluxes and EWs of Hα, and [OII]⁶. We report the obtained fluxes in Table A.1. Regarding the control sample, there are 44 galaxies with NIRSpec spectra. [OIII] is detected in 27 of those spectra, 17 with the prism and 10 with the medium-resolution gratings. In the other 17 galaxies (38%) with spectra, no bright emission lines are detected or the lines fall in detector gaps according to their z_phot.

In Fig. 13, we compare the photometric and the spectroscopic redshifts of the sample of EELGs. We find that photometric redshifts are on average Δz = 0.2 higher than the spectroscopic redshifts (and up to Δz = 0.8 higher). The larger differences are in the cases where the photo-z code confuses a galaxy with Hα falling in the red edge of filter F444W with a galaxy with Hα falling in the red edge of filter F410M or [OIII] falling in the blue edge of F356W. These degeneracies are due to the fact that strong emission lines fall on the edge of observed filters producing the largest differences in redshift in our sample.

Fig. 13. Comparison between photometric and spectroscopic parameters for the subsample of EELG candidates with NIRSpec spectrum. Photometric vs spectroscopic redshift with medium resolution (MR) grating (circles) and prism (triangles). The magenta dashed line is a difference of 0.5 in redshift. The black solid line is the 1:1 relation.

We also compare the line flux [OIII]+Hβ from the spectra and the photometric method (see the top panel in Fig. 14). The prism spectra show line flux within a factor of 0.1 and 1.77, with a mean value of 0.9, and most of the galaxies are within a factor of 3, compared to the photometric-based fluxes. Similar values are found in the medium-resolution spectra where the mean difference is a factor of 0.96, with values ranging from a factor of 0.27 to 1.72. We perform a K-S test and we find that the differences are significant between the spectroscopic fluxes and the ones estimated in Sect. 3.1. The disagreement could be due to a mismatch of the photometric calibration or residual slit loss correction. But the fact they have [OIII] in emission indicates that the main contributor to the flux excess in the observed fluxes is in fact dominated by strong emission lines. We also compare the EW([OIII]+Hβ) from photometry and the spectra (see the bottom panel in Fig. 14) and we find they also agree within a factor of 3, especially in the cases where the spectra are taking from the prism configuration. In the cases with MR spectra, the differences may be larger due to the lack of detecting continuum which makes it more uncertain to quantify the EW of a given line. We omitted the errorbars for data points from MR spectra (green circles) because they are huge and unrealistic given that the continuum is not detected.

Fig. 14. Comparison between photometric and spectroscopic parameters for the subsample of EELG candidates with NIRSpec spectrum. Top (bottom) panel: [OIII]+Hβ flux ([OIII]+Hβ EW) based on photometry vs spectra with grating (circles) and prism (triangles). In both panels, the black dashed lines are a difference by a factor of 3. The black solid line is the 1:1 relation.

In order to probe the mean properties of the sample, we performed two stacked spectra. One includes all galaxies with prism spectra and the other includes the galaxies with medium-resolution spectra. To perform the stacking, we calculate the median of all spectra in each bin of rest-frame wavelength. The results of this stacked analysis are presented in the following sections along with measurements in individual galaxies.

4. Results and discussion

4.1. Young starbursts

We use the Mass-Excitation (MEx, Juneau et al. 2014) diagram to study the ionization source in the sample of EELGs. We note that MEx is usually not preferred to other line ratio diagnostics given the added biases of using stellar masses (e.g., Coil et al. 2015; Cleri et al. 2023a), but we keep our analysis based on the MEx diagnostic because we can apply it to a large amount of galaxies in the sample. Since we can not resolve [OIII] and Hβ from the photometry, we focus our analysis on the subsample with spectroscopy. In Fig. 15, we show the position of our sample of EELGs in the classical MEx diagram. We find that their position is consistent with the location of metal-poor starbursts with high log([OIII]λ5007/Hβ)≳0.4 − 1. This region is populated by local analogs of these high-z systems. For example, in Fig. 15, we include the sample of EELGs up to z ∼ 1 from Amorín et al. (2015) to compare with our sample. We notice that the sample of low-z EELGs are in the same locus with slightly higher stellar masses. We also show that the stacked spectra show similar [OIII]/Hβ ∼ 7 ratios in both cases using prism or medium-resolution. This indicates the high ionization conditions in this subsample of EELGs. Regarding the control sample with confirmed spectra, we note in Fig. 15 that they show similar values of log([OIII]/Hβ), even though there are more galaxies where Hβ is an upper limit.

Fig. 15. Mass-Excitation (MEx) Diagram. The EELG candidates with NIRSpec spectra are in green symbols (circles for medium-resolution (MR) and triangles for prism). The stacks in this work are marked by star symbols. Galaxies in the control sample with NIRSpec spectra are shown with magenta circles (medium-resolution mode) and triangles (prism mode) symbols. The black dashed line indicates the demarcation between star-forming galaxies and AGN, according to Juneau et al. (2014). The black dashed-dotted line is the demarcation at z ∼ 2 (Coil et al. 2015). The black circles are EELGs up to z ∼ 1 (Amorín et al. 2015).

The excitation properties of low-mass compact systems are consistent with being dominated by young starbursts. However, we cannot rule out other non-thermal sources (such as AGN) for the entire sample of photometric-selected EELGs using only photometry (e.g., Backhaus et al. 2024; Barro et al. 2024; Cleri et al. 2023b; Davis et al. 2023; Larson et al. 2023b).

In Fig. 16 we show the position of the sample of EELGs and the control sample in the stellar mass-SFR plane. The contours are the 75th, 50th, and 25th percentiles, respectively. We find that the galaxies in the control sample are distributed along the so-called main sequence (e.g., Speagle et al. 2014) with smaller scatter compared to the observed scatter in the EELGs candidates. The sample of EELGs shows larger SFR for a given stellar mass compared with the control sample, which implies higher sSFR, as previously shown in Fig. 11. Some galaxies may be part of the main-sequence but most of the sample shows higher SFR which implies they are likely strong starbursts. This bimodality between main-sequence galaxies and starbursts has been observed in other works at z ∼ 3 − 6.5 (e.g., Rinaldi et al. 2022). Given that our sample covers a wide range of redshifts, we assume no evolution for z > 3 in the normalization of the main-sequence but an increase of the scatter at z > 4 due to more bursty SFH (e.g., Cole et al. 2023). Individual galaxies with NIRSpec spectra are also shown in Fig. 16 and they show a wide range in stellar mass (∼2 dex) and in SFR (∼2 dex). We note that given the wide range of parameters, this subsample is representative of all galaxies in the entire sample as will also be seen in the following Fig. 17.

Fig. 16. SFR–M⋆ diagram for our sample of EELGs (in blue scale contours) and control sample (in red scale contours). The contours are the percentiles 75th, 50th, and 25th respectively. The solid line is the main sequence from Iyer et al. (2018) at z ∼ 6 that covers a mass range from 106 M⊙. The green squares are the subsample of EELGs with NIRSpec spectra.

Fig. 17. Relation between EW([OIII]+Hβ) and sSFR (top panel) and age (bottom panel) for the sample of EELGs and control sample (same as in Fig. 16). The black squares are a sample of EELGs at z = 1.3 − 2.4 (Tang et al. 2019). The EWs of the sample with NIRSpec spectra are from photometry measurements.

EELGs also show systematically higher sSFR than the control sample, as can be seen in Fig. 17 (top panel). We also found an increase of sSFR with EW([OIII]+Hβ). This increase starts for sSFR > 10⁻⁸ yr. This trend is not observed in the control sample where the scatter is larger for a given EW. Similar results are found with the age of the starburst (see the bottom panel in Fig. 17) where the youngest bursts show the largest EW([OIII]+Hβ). We compare these results with a sample of [OIII] emitters at z ∼ 1.3 − 2.4 (Tang et al. 2019) and we find our results are consistent with the trend observed of intermediate-z analogs. A similar trend of EW with sSFR is found in Papovich et al. (2022) for ∼200 star-forming galaxies at z ∼ 1.1 − 2.3 using HST/WFC3 IR grisms. Similarly, the decreasing trend of age with EWs has also been found in other works of emission-line galaxies at z ∼ 7 − 8 Rinaldi et al. (2023).

4.2. Ionizing photon production efficiency ξ_ion

Ionizing photon production efficiency ξ_ion refers to the ratio of ionizing photons emitted by stars within a galaxy to the total number of photons produced by those stars. This efficiency is a crucial parameter to estimate the contribution of different types of galaxies to the reionization process. Galaxies with high ionizing photon production efficiencies are likely to play a more significant role in driving reionization. Scaling relations have been proposed to estimate such parameters based on observables, for example, the ξ_ion–EW([OIII]) relation proposed in Chevallard et al. (2018).

In Fig. 18, we show the scaling relations between EW([OIII]+Hβ) and the stellar mass. A clear negative correlation between emission line EW and stellar mass has been found in the literature at different redshifts (e.g., Tang et al. 2019; Lumbreras-Calle et al. 2022; Matthee et al. 2023). This indicates that galaxies with lower masses tend to have stronger recent star formation events relative to their mass (higher sSFR). However, we note that given the high average value of [OIII]/Hβ we measured, then EW([OIII]+Hβ) is dominated by [OIII], and the strength of [OIII] can depend on other excitation effects (such as metallicity and the hardness of the ionizing spectrum), and not only on stellar mass or sSFR. As shown in Fig. 18, we find that the spectroscopically-confirmed EELGs in our sample follow the trend observed at lower redshift (z ∼ 3, black dashed line in Fig. 18). This trend seems to be followed also by spectroscopically-confirmed galaxies in the control sample, with few cases with higher stellar mass and an overlap region with the EELGs region. This is consistent with other works based on stacking of z > 5.5 galaxies (red squares, Matthee et al. 2023). Overall, galaxies with low stellar masses and high sSFR tend to show higher EWs and then are more likely to show extreme ISM conditions reflected in their extreme EWs.

Fig. 18. Relation between stellar mass and EW([OIII]+Hβ). The contours are the same as in Fig. 16. In green symbols, the subsample of EELG candidates with NIRSpec spectroscopy, while the control sample is in magenta symbols. The blue star is the stack of the EELGs with prism spectra. The black dashed line is the relation at z ∼ 3 and the observed scatter (Llerena et al. 2023). The red squares are stacks of z > 5.5 galaxies (Matthee et al. 2023).

Based on photometry, our EELGs show a mean EW([OIII]+Hβ) = 1231 Å (σ = 421 Å) with values as high as 2932 Å. The stack from prism spectra shows a typical EW([OIII]+Hβ) of 1258 Å. The highest value measured in a single galaxy with NIRSpec/prism spectra is EW([OIII]+Hβ) = 2553 Å. For Hα, the mean EW is 808 Å (σ = 392 Å), which is consistent with the EW(Hα) = 793 Å of the stack. We use the relation in Tang et al. (2019) depending on EW([OIII]λ5007) and we extrapolate the relation to cover the EWs observed in the control sample and in the EELG sample. This relation is consistent with recent observations of galaxies at z = 3 − 5.7 (Boyett et al. 2024). We obtain ξ_ion for both samples. The values are displayed in Fig. 19. We find that roughly 65% of the EELGs (483 out of 736 candidates) show ξ_ion higher than the canonical values for a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013). This subsample of EELGs, in particular the galaxies with lower stellar masses (due to the mild trend we find), is efficiently producing ionizing photons and are ideal laboratories to probe the escape of ionizing photons. Actually, our sample may include strong candidates of reionization galaxies with LyC escape as will be discussed in Sect. 4.4. As expected from the lower EWs of the control sample, they show ξ_ion below the canonical values, and only 3% of the control sample show values within the range of canonical values. The sample of EELGs shows a mean value of log ξ_ion[Hz/erg] = 25.35 while the control sample shows a mean value of log ξ_ion[Hz/erg] = 24.90.

Fig. 19. Relation between stellar mass and ξion for the sample of EELGs and control sample (same as in Fig. 16). The gray shaded region is the canonical ξion value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013), assuming an escape fraction of 0.20. The horizontal dashed lines are the canonical ξion values for an escape fraction of 0.10 and 0.05, respectively.

4.3. Gas-phase metallicity

For the galaxies (20 out of 39) with a detected Hα and Hβ (S/N > 3), we estimated the nebular attenuation assuming Hα/Hβ = 2.79 under case B approximation for T_e = 15 000 K and n_e = 100 cm⁻³ (Pérez-Montero 2017) and considering the Cardelli law (Cardelli et al. 1989). For galaxies with prism, we find a mean value of 0.29 mag (Hα/Hβ = 3.79), while for medium-resolution we find a mean value of 0.21 mag (Hα/Hβ = 3.31). We use these mean values to estimate the nebular attenuation in galaxies with unobserved Hα or undetected Hβ. We corrected the observed fluxes reported in Table A.1 assuming the Reddy et al. (2015) law with R_V = 2.505.

We used the O32 = [OIII]λ5007/[OII]λλ3727,29 calibration from Bian et al. (2018) to estimate the gas-phase metallicity in individual galaxies and in the stacks. In Fig. 20, we show the gas-phase metallicity for the sample of EELGs with NIRSpec spectra. We estimate gas-metallicity only in the subsample of galaxies where [OII] is detected or where we were able to put an upper limit. The galaxies with detected [OII] in prism spectra show a mean redshift of z ∼ 5.39 and a mean gas-phase metallicity of log(O/H) + 12 = 8.07 (∼24% solar). On the other hand, the galaxies with detected [OII] in medium-resolution spectra show a mean redshift of z ∼ 5.96 and a mean gas-phase metallicity of log(O/H) + 12 = 8.00 (∼20% solar). Regarding the stacks, we find consistent values of log(O/H) + 12 = 7.85 (±0.3) and 8.01 (±0.3) for prism and grating stacks, respectively. These values are not different from the typical average metallicities measured in EELG at low and intermediate redshifts (e.g., Amorín et al. 2010, 2014; Maseda et al. 2014; Calabrò et al. 2017; Pérez-Montero et al. 2021; Tang et al. 2021). We also note that the obtained values are consistent with the gas-phases metallicities of the templates ASK 17, 20, 21 with values 0.25−0.35 Z_⊙. We also estimate the metallicities using the calibration for EELGs presented in Pérez-Montero et al. (2021) and we find a mean difference of 0.03 towards lower metallicities using this O32 calibration, which means no significant differences are using both calibrations. We also explore the same method used in Sanders et al. (2021) based on multiple diagnostics to estimate the metallicity minimizing their Eq. (6). We used the diagnostics [OIII]λ5007/Hβ, [OII]λ3727/Hβ, O32 and R23 = $\frac{[\mathrm{OIII}]\lambda \lambda 4959,5007+[\mathrm{OII}]\lambda 3727}{\mathrm{H}\beta }$. We find a mean difference of 0.1 dex compared with using only the O32 calibration from Bian et al. (2018), which would reduce our metallicities to ∼15−20%. Similar differences of up to 0.1 dex are found using the low metallicity branch of the R23 calibration from Papovich et al. (2022) which is tested in HST grism spectra with galaxies at z ∼ 1.1 − 2.3.

Fig. 20. Mass-Metallicity relation for EELGs. Individual galaxies of the EELG sample (control sample) are in green (magenta) symbols. Results from the stacks are represented with red (grating) and blue (prism) stars. Red squares are stacks at z ∼ 5 − 7 from Matthee et al. (2023) and red circles are stacks at z ∼ 3 − 10 from Curti et al. (2024). The black solid line is the MZR at z ∼ 3.3 from Sanders et al. (2021) and its extrapolation to lower stellar masses in the dashed line. The red dashed line is the predicted MZR at z = 5 − 12 from FIRE-2 simulation (Marszewski et al. 2024).

We note that these values are consistent with the mass-metallicity relation at z ∼ 3.3 reported in Sanders et al. (2021), and comparable with those predicted in FIRE-2 simulation at z ∼ 5 − 12 (Marszewski et al. 2024). Compared with other recent works, our results are comparable to those of Matthee et al. (2023) at z ∼ 5 − 7 and Curti et al. (2024) at z ∼ 3 − 10 in the similar stellar mass range. The metallicities for the galaxies in the control sample tend to also follow the MZR as can be seen in Fig. 20, with slightly higher values in the high mass end compared to the EELG sample.

Overall, we find that our sample includes galaxies with subsolar metallicities at around ∼20−25% solar, consistent with the mass-metallicity relation, which suggests they do not show different metallicities compared with the general population at similar redshifts.

4.4. Conditions for the escape of ionizing photons

In this section, we explore indirect tracers of the escape of ionizing photons to understand if our sample of young starburst is made of strong LyC leakers. We use indirect tracers since in galaxies in the EoR, as the ones included in this paper, a direct probe of the LyC is not possible due to the high fraction of neutral hydrogen and then indirect tracers are required (e.g., Mascia et al. 2024).

We first explore the relationship with the star formation surface density (e.g., Flury et al. 2022). As displayed in Fig. 21, we find that our EELG candidates tend to show higher EWs while increasing the Σ_SFR. Regarding the spectroscopically-confirmed EELGs, we find that they show high Σ_SFR at high EW. Regarding the control sample, we show that they have Σ_SFR ≲ 10 M_⊙ yr⁻¹ kpc⁻², while the EELGs show up to larger values. We note that Σ_SFR ≳ 10 M_⊙ yr⁻¹ kpc⁻² is the range in which the fraction of detected LyC leakers is > 40% according to observations of local galaxies (Flury et al. 2022). This indicates that galaxies with higher EWs and concentrated SFR have the conditions that should facilitate the escape of LyC photons in agreement with the $f_{\mathrm{esc}}\propto {\mathrm{\Sigma }}_{\mathrm{SFR}}^{0.4}$ relation described in Naidu et al. (2020).

We explore O32 as an alternative indicator of LyC leaking (e.g., Izotov et al. 2020). In Fig. 22 (top panel), we show the relation of O32 with EW([OIII]λ5007). We show there is an increase of EW with O32. We find that 9 out of 12 (75%) spectra with detected [OII] and [OIII] in their spectra show O32 > 5. The galaxies with EW([OIII]λ5007) ≳ 400 Å are the ones that show O32 > 5 which indicates strong LyC leaker candidates based on local galaxies (Flury et al. 2022). We note that a value of O32 = 5 corresponds to a metallicity of 0.27 Z_⊙ using the Bian et al. (2018) calibration, which is approximately the mean metallicity of the sample of EELGs. The trend we find is consistent with analogs at lower redshifts (Tang et al. 2019). The stack with galaxies with prism spectra show O32 = 14.6 and EW([OIII]) = 874 Å (blue star in Fig. 22), which indicates that is a population of galaxies with conditions of escape of ionizing photons and deeper observations are needed in individual galaxies to probe their nature. Regarding the control sample, we note a lower fraction of strong candidates to be strong leakers.

Fig. 21. Relation between EW([OIII]+Hβ) and ΣSFR for the sample of EELGs and control sample (same as in Fig. 16). The green squares are the subsample of EELGs with NIRSpec spectra. The red dashed line is the lower EW limit of the ASK templates in this paper.

Fig. 22. Relation between O32 and EW([OIII]λ5007) (top panel) and ΣSFR (bottom panel). We include lower limits due to the non-detection of [OII]. In green symbols, the sample of EELGs with NIRSpec spectra, while in magenta symbols, the control sample with NIRSpec spectra. The vertical dashed line is the separation between weak and strong LyC leaker according to Flury et al. (2022). On the top panel, black squares are a comparison sample at z = 1.3 − 2.4 (Tang et al. 2019), and the blue dashed-dotted line is the relation found in Boyett et al. (2024) for a sample of galaxies at z = 3 − 9.5.

As expected from previous results, we find an increase of Σ_SFR with O32 (bottom panel in Fig. 22). This indicates that indeed the EELGs with higher concentrations of star formation provide the feedback necessary to clear LyC escape paths in the ISM. Similarly, we note a lower fraction of galaxies in the control sample with a high concentration of star formation compared to the sample of EELGs. We find also a lower fraction of strong candidates of strong leakers in the control sample.

5. Conclusions

We use NIRCam photometry and empirical templates to select EELGs. We demonstrate that NIRCam can identify a large sample of previously unknown EELGs in a wide range of redshifts showing unique properties similar to the rare metal-poor local starburst. The proposed color selections can effectively identify galaxies with EW([OIII]+Hβ) > 680 Å at all redshifts targeted in this work. We use the broad-band filters F277W, F356W, and F444W and the medium-band filter F410M to select EELGs at 4 ≲ z ≤ 9. We selected a sample of 1000 candidates. 47 of them have already NIRSpec spectra. We define a control sample to compare the properties with the sample of EELGs. We use BAGPIPES to estimate their physical properties considering a delayed τ-model for the SFH and a nebular component to model the intense emission lines. We also use the F200W filter to estimate the physical size of the young stellar populations. We find:

• Our sample of EELGs shows a mean stellar mass of 10^7.84 M_⊙ with high sSFRs with a mean value of 10^−7.03 yr⁻¹. They are young with a mean value of the time after the onset of star formation of 45 Myr.

• Compared with the control sample, they show slightly lower stellar masses but similar stellar metallicities and dust attenuation. The larger differences are in the sSFR, ages, and ionization parameter. This suggests they may have similar underlying stellar populations but the young starburst may be the reason for the extreme emission lines.

• Based on their confirmed emission lines, we find that they are in the locus of metal-poor starbursts with high log([OIII]/Hβ) ≳ 0.4−1 which indicates that star-formation may be the dominant source of ionization in these galaxies. The stacking analysis confirms this result with a small sample of NIRSpec spectra. The starburst nature of these galaxies is also observed in their position above the main-sequence.

• Based on the photometric fluxes, we find a mean rest-frame EW([OIII]+Hβ) of 1231 Å for our sample of EELGs. The EW([OIII]+Hβ) shows an increase with sSFR and a decrease with age and stellar mass. In the control sample, the scatter is larger in these relations.

• We use the EWs to estimate the ionizing photon production efficiency and we found that roughly 65% (483 out of 736 candidates) of the sample of EELGs show ξ_ion values higher than the canonical values, which implies they are efficiently producing ionizing photons and are ideal laboratories to probe the escape of ionizing photons. The sample of EELGs shows a mean value of log ξ_ion[Hz/erg] = 25.35.

• We find sub-solar gas-phase metallicities for the sample of EELGs based on the O32 calibration with a mean value of 20−25% solar. They follow the MZR at z ∼ 3.3 which suggests they do not show different gas-phase metallicities than the general population of galaxies.

• We find an increase of EW with Σ_SFR. The sample of EELGs can reach Σ_SFR > 10 M_⊙ yr⁻¹ kpc⁻² which indicate they are strong candidates of LyC leakers. Another indirect indicator is the high values of O32 > 5 (corresponding to gas-phase metallicities < 0.27 Z_⊙) that can be reached for some galaxies in the sample. This indicates that they may have the conditions to facilitate the escape of ionizing photons.

Data availability

Movie associated to Fig. 2 is available at https://www.aanda.org

The full Tables 2 and 3 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/691/A59

Acknowledgments

We thank the anonymous referee for the detailed review and useful suggestions that helped to improve this paper. We wish to thank all our colleagues in the CEERS collaboration for their hard work and valuable contributions to this project. We wish to thank to Jorge Sánchez-Almeida for providing us with the ASK templates used in this paper. MLl acknowledges support from the PRIN 2022 MUR project 2022CB3PJ3 – First Light And Galaxy aSsembly (FLAGS) funded by the European Union – Next Generation EU. RA acknowledges the support of ANID FONDECYT Regular Grant 1202007. RA acknowledges financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. This work has made extensive use of Python packages astropy (Astropy Collaboration 2018), numpy (Harris et al. 2020), Matplotlib (Hunter 2007) and LiMe (Fernández et al. 2024).

References

Appendix A: Fluxes of the subsample of EELGs with NIRSpec spectra

All Tables

All Figures

	Fig. 1. CEERS NIRCam filters. The black lines are the estimated observed wavelength at a given redshift for [OII]λ3727 (dotted-dashed), [OIII]λ5007 (dashed) and Hα (solid). The filter response functions for NIRCam filters are arbitrarily scaled.
In the text

Fig. 2. Illustration of the ASK templates at z ∼ 6 in the NIRCam filters. We illustrate which filters are considered to estimate the NIRCam color due to the limited wavelength coverage of the templates. In red are the NIRCam filters where the templates are outside the filter, while in purple are the filters where we estimate the NIRCam colors. Three ASK templates are illustrated as examples (ASK0, ASK10, and ASK15). A movie associated to this figure is available online.

In the text

	Fig. 3. Redshift evolution of NIRCam colors of ASK templates. The major ASK, that is non-EELGs, classes are in gray-scale colors, according to the legend. The minor ASK classes, that is EELGs, are in colors according to the legend. The magenta and cyan-shaded regions are the z ranges and color cuts used to select EELGs.
In the text

Fig. 4. Color diagrams to select EELGs at each redshift range. The color cuts are marked in magenta (top panel) and cyan (bottom panel) shaded regions in each panel according to Eqs. (1) and (2), respectively. The selected EELG candidates are in red circles, while parent samples at similar redshifts are in black circles. The sources with red crosses were discarded due to contaminated photometry as explained in the main text.

In the text

Fig. 5. Images and SED of CEERS1455, an example of an EELG candidate at zphot = 4.66. In the top panels, images (4″ × 4″) of the galaxy in filters (from left to right) RGB image, HST/F606W, HST/F814W, JWST/F277W, JWST/F356W, JWST/F410M, and JWST/F444W. The RGB image is made with Red = F410M, Green = F356W, and Blue = F277W. On the bottom panel, SED of the galaxy. The blue (black) squares are the NIRCam (HST) photometric points. NIRCam (HST) filters are in rainbow-scale (gray-scale) colors, according to legend. The filter transmission is plotted in arbitrary units depending on the observed flux density. The black solid line is the SED model and the inset text summarizes its physical properties (see description in Sect. 3.2). The vertical dashed lines represent the position of Lyα (black), [OIII] (blue), and Hα (red), according to their zphot.

In the text

	Fig. 6. Images and SED of CEERS76140, an example of an EELG candidate at zphot = 7.48. Symbols as in Fig. 5.
In the text

	Fig. 7. Distribution of the magnitude F200W for the sample of EELGs (unfilled histograms) and the parent samples at similar redshifts (filled histograms) in the two redshift ranges considered for the selection in blue and red, respectively. They are normalized to the total number of galaxies in each parent sample.
In the text

	Fig. 8. Distribution of photometric redshift of the final sample of EELG candidates (in blue) and the control sample (in red). The gray-shaded region is the gap between the redshift range for selection. The normalization is made for the total of galaxies in the parent sample at similar redshift and S/N.
In the text

	Fig. 9. RGB (F410M–F356W–F277W) images of the 100th brightest F410M EELG candidates in our sample. Each image has a size of 4″ × 4″.
In the text

	Fig. 10. Photometric fluxes of [OIII]+Hβ (left panels) and Hα (right panels) for the sample of EELG candidates (in blue) and the control sample (in red). From top to bottom: Integrated flux of the line in units of erg s−1 cm−2, continuum in the line wavelength in units of erg s−1 cm−2 Å−1, and rest-frame EW in units of Å.
In the text

	Fig. 11. Distribution of the physical parameters based on the SED fitting for the sample of EELG candidates (in blue) and the control sample (in red). The vertical dotted (H−z sample) and dotted-dashed (M−z sample) lines are the mean values of each parameter.
In the text

	Fig. 12. Size-mass relation for the EELG (in blue) and the control (in red) samples. The grey solid line is the size-mass relation at z = 5 from Morishita et al. (2024) and the shaded region the 3σ-scatter.
In the text

	Fig. 13. Comparison between photometric and spectroscopic parameters for the subsample of EELG candidates with NIRSpec spectrum. Photometric vs spectroscopic redshift with medium resolution (MR) grating (circles) and prism (triangles). The magenta dashed line is a difference of 0.5 in redshift. The black solid line is the 1:1 relation.
In the text

	Fig. 14. Comparison between photometric and spectroscopic parameters for the subsample of EELG candidates with NIRSpec spectrum. Top (bottom) panel: [OIII]+Hβ flux ([OIII]+Hβ EW) based on photometry vs spectra with grating (circles) and prism (triangles). In both panels, the black dashed lines are a difference by a factor of 3. The black solid line is the 1:1 relation.
In the text

Fig. 15. Mass-Excitation (MEx) Diagram. The EELG candidates with NIRSpec spectra are in green symbols (circles for medium-resolution (MR) and triangles for prism). The stacks in this work are marked by star symbols. Galaxies in the control sample with NIRSpec spectra are shown with magenta circles (medium-resolution mode) and triangles (prism mode) symbols. The black dashed line indicates the demarcation between star-forming galaxies and AGN, according to Juneau et al. (2014). The black dashed-dotted line is the demarcation at z ∼ 2 (Coil et al. 2015). The black circles are EELGs up to z ∼ 1 (Amorín et al. 2015).

In the text

	Fig. 16. SFR–M⋆ diagram for our sample of EELGs (in blue scale contours) and control sample (in red scale contours). The contours are the percentiles 75th, 50th, and 25th respectively. The solid line is the main sequence from Iyer et al. (2018) at z ∼ 6 that covers a mass range from 106 M⊙. The green squares are the subsample of EELGs with NIRSpec spectra.
In the text

	Fig. 17. Relation between EW([OIII]+Hβ) and sSFR (top panel) and age (bottom panel) for the sample of EELGs and control sample (same as in Fig. 16). The black squares are a sample of EELGs at z = 1.3 − 2.4 (Tang et al. 2019). The EWs of the sample with NIRSpec spectra are from photometry measurements.
In the text

Fig. 18. Relation between stellar mass and EW([OIII]+Hβ). The contours are the same as in Fig. 16. In green symbols, the subsample of EELG candidates with NIRSpec spectroscopy, while the control sample is in magenta symbols. The blue star is the stack of the EELGs with prism spectra. The black dashed line is the relation at z ∼ 3 and the observed scatter (Llerena et al. 2023). The red squares are stacks of z > 5.5 galaxies (Matthee et al. 2023).

In the text

Fig. 19. Relation between stellar mass and ξion for the sample of EELGs and control sample (same as in Fig. 16). The gray shaded region is the canonical ξion value given by a simple stellar population at constant SFR over 100 Myr (Robertson et al. 2013), assuming an escape fraction of 0.20. The horizontal dashed lines are the canonical ξion values for an escape fraction of 0.10 and 0.05, respectively.

In the text

Fig. 20. Mass-Metallicity relation for EELGs. Individual galaxies of the EELG sample (control sample) are in green (magenta) symbols. Results from the stacks are represented with red (grating) and blue (prism) stars. Red squares are stacks at z ∼ 5 − 7 from Matthee et al. (2023) and red circles are stacks at z ∼ 3 − 10 from Curti et al. (2024). The black solid line is the MZR at z ∼ 3.3 from Sanders et al. (2021) and its extrapolation to lower stellar masses in the dashed line. The red dashed line is the predicted MZR at z = 5 − 12 from FIRE-2 simulation (Marszewski et al. 2024).

In the text

	Fig. 21. Relation between EW([OIII]+Hβ) and ΣSFR for the sample of EELGs and control sample (same as in Fig. 16). The green squares are the subsample of EELGs with NIRSpec spectra. The red dashed line is the lower EW limit of the ASK templates in this paper.
In the text

Fig. 22. Relation between O32 and EW([OIII]λ5007) (top panel) and ΣSFR (bottom panel). We include lower limits due to the non-detection of [OII]. In green symbols, the sample of EELGs with NIRSpec spectra, while in magenta symbols, the control sample with NIRSpec spectra. The vertical dashed line is the separation between weak and strong LyC leaker according to Flury et al. (2022). On the top panel, black squares are a comparison sample at z = 1.3 − 2.4 (Tang et al. 2019), and the blue dashed-dotted line is the relation found in Boyett et al. (2024) for a sample of galaxies at z = 3 − 9.5.

In the text