© A. Lumbreras-Calle et al. 2022

1. Introduction

Galaxies undergoing strong events of star formation and presenting compact morphologies have been identified and analyzed since the middle of the 20th century, with pioneering works by Haro (1956), Zwicky (1966), and Markarian (1967). One of the most fruitful developments in the detection of this type of galaxy was the use of an objective prism with which strong emission lines were identified (Markarian 1967). Early analyses of this class of galaxies reported their usually compact morphologies and blue colors, with their spectra resembling those of galactic HII regions, thus the label “HII galaxies” (Sargent & Searle 1970; Melnick et al. 1985). Low metallicities were also identified, as well as a recent enhancement of star formation (Searle & Sargent 1972).

With the advent of modern CCD detectors and wide-field surveys, new windows into the analysis of these objects were opened. An observational alternative for the identification of emission line galaxies (ELGs) is the use of multiband surveys; some examples are CADIS (Meisenheimer et al. 1998), COMBO-17 (Wolf et al. 2001, 2003), HiZELS (Geach et al. 2008; Sobral et al. 2015), ALHAMBRA (Moles et al. 2008), MUSYC (Cardamone et al. 2010), SHARDS (Pérez-González et al. 2013), COSMOS medium band (Taniguchi et al. 2015), PAU (Benítez et al. 2009), SC4K (Sobral et al. 2018), J-PLUS (Cenarro et al. 2019), S-PLUS (Mendes de Oliveira et al. 2019), and J-PAS (Benitez et al. 2014). Using a set of several narrow-band or medium-band filters, sometimes also broadband filters, these surveys easily detect emission lines by identifying an excess of flux in one of the bands (e.g., Hippelein et al. 2003; Maier et al. 2003; Gronwall et al. 2007; Sobral et al. 2013, 2015; Cava et al. 2015; Lumbreras-Calle et al. 2019a; Spinoso et al. 2020; Vilella-Rojo et al. 2021). While usually providing smaller fields of view and less wavelength coverage than broadband or objective prism surveys, they can reach deeper magnitudes than the former and higher spectral resolution than the latter. These surveys have yielded large samples of ELGs and sometimes simultaneously identified different emission lines and covered a wide redshift range.

Focusing on higher redshifts, the search for ELGs has been performed using deep, broadband surveys, especially those with data from the Hubble Space Telescope (HST), such as GOODS, COSMOS, or CANDELS (Kakazu et al. 2007; van der Wel et al. 2011; Maseda et al. 2018). They found an increase in the density of ELGs at higher redshifts, in coherence with the evolution of the cosmic star formation history (SFH; Madau & Dickinson 2014). At very high redshift (z > 6), typical low-mass galaxies are expected to be extreme emission line galaxies (EELGs; see, e.g., Labbé et al. 2013; Smit et al. 2014; Endsley et al. 2021), which play a crucial role in the reionization of the Universe (Ouchi et al. 2009; Bouwens et al. 2015). From an observational point of view, space-based grism spectroscopy in some of the HST deep fields provides an even greater insight into the nature of EELGs (Pirzkal et al. 2013). This approach has only been feasible in small fields, which is enough for higher redshifts, but insufficient for local Universe studies. Upcoming survey telescopes such as Euclid, the Nancy Grace Roman Telescope, and the Chinese Space Station Telescope will extend this grism analysis to fields covering thousands of square degrees. It is therefore vital to gather samples of EELGs that can be used as a reference for this type of survey, testing their limitations and providing a readily available reference sample.

The precise definition of EELGs varies in the literature. It is very often described as a threshold in the rest-frame equivalent width (EW) of the [OIII]5007 line, but values range from 100 Å in Amorín et al. (2015) and Pérez-Montero et al. (2020) to 300 Å in Jiang et al. (2019), with several surveys providing most objects above 500 Å (van der Wel et al. 2011; Yang et al. 2017; Maseda et al. 2018).

Recently, several projects have aimed at detecting samples of EELGs at low redshift using photometric surveys, such as Cardamone et al. (2009), Yang et al. (2017), and Senchyna & Stark (2019). They have demonstrated the power of broadband surveys to identify this type of galaxy with relatively high purity, while covering much wider areas than previous surveys, allowing the study of rare populations. Studies using narrow-band filters have also proven useful in the detection of ELGs over relatively modest fields (Kellar et al. 2012; Salzer et al. 2020) or brighter galaxies in very wide fields (Cook et al. 2019).

In addition to photometric surveys, the extensive database provided by the SDSS spectra has been thoroughly explored in the search of EELGs (Izotov et al. 2011) or extremely metal-poor galaxies (Sánchez Almeida et al. 2016). The physical information for individual galaxies that is accessible through this database is unparalleled and shows, for instance, how these local samples can serve as analogs for the very high redshift systems that cause the reionization of the Universe (Izotov et al. 2021). Nevertheless, a drawback of this data set is its poorly defined selection function, with galaxies often observed in surveys targeting other types of sources (such as quasars). As a result, no meaningful study regarding number densities, for example, can be performed, except for the very bright end, where the completeness is high (Strauss et al. 2002). A very deep spectroscopic survey (with exposure times of ∼10 h) that is free from these selection biases is the MUSE Hubble ultra-deep field (Bacon et al. 2017). With the MUSE integral field unit, a large sample of ELGs was detected, some of which were extremely faint and were not detected in the HST data. Nevertheless, this work covered a small area (∼10 arcmin²) that is insufficient for studies of uncommon sources at low redshift.

In this context, the Javalambre Photometric Local Universe Survey (J-PLUS; Cenarro et al. 2019) combines multiband observation with a very wide field, making it a unique tool for the identification of very rare emission line objects. In addition, as an imaging survey, it allows us to analyze the emission in the whole galaxy, not only within the the aperture in which the spectra are extracted, which is a limitation found in fiber-fed spectroscopy surveys (such as SDSS, GAMA, or DESI). This photometric survey has already produced relevant results in the detection and analysis of extragalactic systems with emission lines in their spectra. In the local Universe, Logroño-García et al. (2019) and Vilella-Rojo et al. (2021) have analyzed Hα emitters and measured the local star formation main sequence accurately. In addition, Spinoso et al. (2020) performed a search for high-redshift quasars with Lyman α emission in J-PLUS. They reported an unprecedented number of very high luminosity sources.

In the present work, we select low-redshift EELGs in J-PLUS by identifying objects showing intense emission in the J0515 filter, as illustrated in Fig. 1. We cover a much wider area than previous multiband surveys, reach deeper magnitudes than wide-field spectroscopic surveys, and are able to identify emission lines with much higher precision than broadband surveys. We aim at providing a new sample that will uncover many unclassified EELGs in the Local Universe, which can complement previous photometric and spectroscopic studies and serve as targets for follow-up observations.

Fig. 1. Illustrative example of an extreme [OIII] emitter, the J-PLUS source 67834-5013 (RA = 239.1019, Dec = 48.1127, zspec = 0.050). Left panel: color composite of the galaxy, obtained from the gri J-PLUS images. The sky location of the source is shown as a black dot, and the white ellipse indicates the three effective radius contours for the source. Right panel: J-PLUS 12-band PSFCOR photometry of the galaxy. The squares show the five SDSS-like filters (ugriz), and circles indicate the seven medium-band filters (J0378, J0395, J0410, J0430, J0515, J0660, and J0861). The solid line shows the spectra from SDSS with a downgraded resolution of R ∼ 180, normalized to the flux in the filter J0660. The location of the most prominent emission lines is marked: [OIII] (traced by J0515 and g), Hα (traced by r), and [OII] (traced by J0395).

The paper is structured as follows: In Sect. 2 we describe the J-PLUS survey and database, as well as our procedure for selecting the candidate sample and removing contaminants. In Sect. 3 we described the spectral energy distribution (SED) fitting analysis and its results, along with the final definition of the EELG sample and the visual morphology of the galaxies. Finally, we compare these results with the available spectra for the sample. We discuss the results in Sect. 4, testing the sample selection and the number densities of EELGs we find, and reviewing the main physical properties of the EELG sample. Finally, in Sect. 5 we summarize our work and present the conclusions.

Throughout this paper, all quoted magnitudes are in the AB system, all logarithms used are in base 10, and coordinates refer to the 2000 equinox. All mentions of the [OIII] line, unless specified otherwise, refer to the [OIII]5007+4959 doublet. The median is used to show the typical value of a magnitude, while the upper and lower limits refer to the 16th and 84th percentiles of the distributions. We assumed a ΛCDM cosmology with with xsΩ_Λ = 0.7, Ω_M = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹. In the SED fitting we use the Salpeter (1955) initial mass function (IMF).

2. Database and sample selection

2.1. J-PLUS second data release

J-PLUS¹ is being conducted at the Observatorio Astrofísico de Javalambre (OAJ; Cenarro et al. 2014) using the 83 cm Javalambre Auxiliary Survey Telescope (JAST80) and T80Cam, a panoramic camera of 9.2k × 9.2k pixels that provides a 2deg² field of view, with a pixel scale of 0.55 arcsec pix⁻¹ (Marín-Franch et al. 2015). The J-PLUS filter system is composed of 12 passbands (Table 1) spanning the full optical range (3500 − 10 000 Å), with seven medium-band and five broadband filters. The J-PLUS observational strategy, image reduction, and main scientific goals are presented in Cenarro et al. (2019). The J-PLUS photometric calibration is described in López-Sanjuan et al. (2019) and López-Sanjuan et al. (2021).

This work is based on the second data release (DR2) of J-PLUS. It covers 2176 square degrees (1941 deg² after masking), with 1088 individual images of 2 deg². The survey reaches limiting magnitudes ranging from 21.8 mag in r to 20.5 mag in z, with 21.0 mag in the J0515 filter (considering 5σ detection).

In order to create the galaxy sample we describe in this work, we ran queries using the Astronomical Data Query Language (ADQL) interface in the J-PLUS database. In the queries, we used the PSFCOR photometry of the catalog, which is designed to better capture the colors of the objects rather than their total flux (Molino et al. 2019). Briefly, PSFCOR photometry is measured considering elliptical apertures with a semimajor axis equal to the Kron radius (Kron 1980) in the reference band (r), half the size of the AUTO aperture. Then, the photometry in the rest of the bands is measured, correcting for the different point spread functions (PSFs), to produce accurate colors (see Molino et al. 2019 and Hernán-Caballero et al. 2021 for more details). This photometry, while optimal for comparing fluxes in different bands, underestimates the total flux of galaxies with respect to the AUTO photometry by 0.5 mag on average (González Delgado et al. 2021). Therefore, in the rest of the analysis of this work (unless specified otherwise), we rescaled the PSFCOR photometry to the AUTO photometry using the ratio of AUTO and PSFCOR fluxes in the r band. In addition, we apply a galactic extinction correction as provided in the J-PLUS database.

2.2. Sample selection

2.2.1. Selection of intense J0515 emitters

The first step in the process is to query the J-PLUS database for objects that show a large excess of emission in the J0515 filter because all extreme [OIII]4959+5007 emitters between redshift 0.007 and 0.06 are expected to show strong emission in this band (Fig. 1). The precise definition of “strong emission” is discussed in Sect. 3.1.1, but given the width of the J0515 filter, we aimed at selecting objects with EW(J0515) > 200 Å in order to secure clear detections. First, in order to compute the excess of emission corresponding to the line flux, we estimated the continuum underneath it. To perform a simple query, we assumed that the continuum flux (in erg s⁻¹ cm⁻²) at the wavelength of J0515 is roughly equivalent to the flux of the r band. This obviously depends on the slope of the SED and on the relative intensity of the Hα line, but it is enough to draw a first sample of emitters, from which we later select the most extreme ones. A possible bias against extreme Hα emitters introduced by this assumptions is addressed in Sect. 3.2. Specifically, we imposed the following condition:

$$\begin{array}{c}\hfill \frac{F(J0515)-F(r)}{F(r)}>1,\end{array}$$(1)

which, when we assume that r traces the continuum, would imply that the emission lines that lie within J0515 have an EW ≳ 200 Å. A precise estimation of the EW is performed in Sect. 3.1.2, when accurate estimates of the continuum have been computed. In addition to this cut, we imposed other constraints to ensure the quality of the sample. We selected objects with r < 22 mag and a signal-to-noise ratio (S/N) higher than 3 in the r, g, and J0515 bands. We also removed objects with a FLAGS parameter indicating poor photometric quality for those bands (saturation, proximity to a bright star or the edge of an image, etc.), except when the flag indicated a close companion or deblending. We preferred to draw a broad sample in order to apply further cuts only later, motivated by the properties and limitations of the sample. The full ADQL query is reproduced in Appendix A.

Nevertheless, in order to analyze the selected objects, we separated them into four groups, considering their CLASS_STAR value (> = 0.5 and < 0.5) and their deblending flag (either no deblending or some deblending). This allowed us to better identify different types of spurious detections that only appear in some of the groups, and to develop techniques to address them specifically. This first rough selection yielded a sample of 30 336 objects with flux excess in the J0515 filter.

2.2.2. Cross-match with SDSS spectra

Several types of sources can show excess in the J0515 filter, such as ELGs, quasi-stellar objects (QSOs), foreground stars, or any object with spurious photometric measurements. In addition, ELGs and QSOs at different redshifts can show emission in this filter, but due to different emission lines. Because we are only interested in low-redshift ELGs (z < 0.06, where the [OIII] line causes the excess emission), we purged the sample from all other sources. To test the accuracy of the methods for removing contaminants, we cross-matched the excess-flux sample with the 16th data release of the SDSS spectroscopic database (Ahumada et al. 2020). This survey covers a large fraction (∼80%) of the J-PLUS DR2 footprint (more than any other spectroscopic survey) and provides a spectral classification for the sources. We considered a 3 arcsec maximum separation in the cross-match and identified 2561 objects.

2.2.3. Removal of contaminants via J-PLUS photometry

When performing a visual inspection of the samples with a blending flag, we identified several objects that were classified as galaxies, but in fact were stars measured with large apertures. These stars are included in our sample because of a spurious deficit in r emission, which implies a high J0515 to r ratio. Therefore, removing candidates with low r-band flux compared to redder filters (such as J0861 and z) allowed us to eliminate these detections and keep all other candidates. Additionally, this removed several QSOs with red (r − z) and (r − J0861) colors. Specifically, we rejected the 4491 objects that fulfilled (r − z) < 1 and (r − J0861) < 1. The selection was confirmed by verifying that none of the 148 objects with SDSS spectra removed in this step was a low-redshift galaxy (and also by inspecting a subsample of their images and SEDs).

Another subsample of contaminants are QSOs at z ∼ 1.7 with a prominent C III]1909 emission line, which at that redshift creates an excess of flux in the J0515 filter, mimicking a low-redshift ELG. In these cases, the object often also shows excess in the J0430 filter, corresponding to the strong C IV 1549 line. These emission lines have been used to detect AGNs in narrow-band surveys, for example, in Stroe et al. (2017a,b). We therefore removed the 255 objects with (J0515 − J0430) > 0. Out of these, 46 candidates have SDSS spectra, and none of them is a low-redshift galaxy.

The sample with J0515 excess after these cleaning steps includes 25 590 sources.

2.2.4. Removal of contaminants: WISE photometry and object sizes

In an early analysis, it became clear that J-PLUS photometry alone was not enough to distinguish precisely between different types of sources, especially at fainter magnitudes. Following the method performed by Spinoso et al. (2020), we cross-matched our excess flux selection with the infrared (IR) data from the Wide-Field Infrared Survey Explorer (WISE) satellite. This mission performed a very wide survey in four IR bands, from 3.4 to 22 μm. We used the unWISE catalog (Schlafly et al. 2019), which only has data for the bluest bands (W1 and W2), but provides a higher depth and spatial resolution. Out of the 25 590 sources with excess flux at this stage, 19 922 are present in the unWISE catalog. We removed 71 stars from this sample that were identified using Gaia DR2 data available in the J-PLUS tables.

After testing several alternatives, we found that the most precise way of separating the low-redshift [OIII] emitters from the other types of objects using unWISE data was the g − W2 versus r − W1 diagram (Fig. 2). The separation in this diagram of low-redshift galaxies and quasars can be understood regarding the different components that dominate the IR light emitted by these sources. In galaxies, the W1 and W2 filters cover a local minimum in their SEDs at a longer wavelength than the peak of emission from low-mass stars, but at a shorter wavelength than the peak of dust emission (see, e.g., Boquien et al. 2019). This results in blue optical W1/W2 colors. In low-redshift quasars, the W1/W2 bands are dominated by dust emission that is much brighter than their optical emission (dominated by the accretion disk), resulting in very red optical W1/W2 colors. Even for higher-redshift quasars, where the W1/W2 range is affected by a combination of disk and dust emission, the optical W1/W2 colors remain slightly red or close to zero (see, e.g., Hernán-Caballero et al. 2016), which is different from low-redshift galaxies. We plot all sources with available W1 and W2 photometry in this diagram and indicate the source type in different colors for those with available SDSS spectroscopy. It becomes clear that the low-z [OIII] emitters are clustered in a very specific region of Fig. 2, while the vast majority of stars, QSOs, and higher-redshift galaxies are excluded. We define a set of limits in Fig. 2 to select a sample of candidate EELGs, preserving a very high purity and completeness (a full account is provided in Sect. 2.3). The selected objects are located in the area that fulfills the following set of equations (blue lines in Fig. 2):

$$\begin{array}{cc}\hfill g-W2& >1.1·(r-W1)-0.8\hfill \\ \hfill g-W2& <-1.3·(r-W1)-0.19\hfill \\ \hfill g-W2& <1.1·(r-W1)+1.6\hfill \\ \hfill g-W2& >-1.3·(r-W1)-4.5.\hfill \end{array}$$(2)

Fig. 2. Color-color diagram used to separate low redshift ELGs from other types of sources. The g and r data are taken from J-PLUS DR2, and the W1 and W2 from the unWISE catalog. All candidates from the sample with J0515 excess flux are shown in dark grey dots. Sources with SDSS spectra are shown in colors according to their nature: red squares for QSOs, blue circles for galaxies at z < 0.06, brown triangles for galaxies at z > = 0.06, and green diamonds for stars. Filled symbols (or outlines only) are used if the sources are brighter (fainter) than r = 20 mag. The objects selected as candidates to be EELGs are enclosed by the blue lines, while those to the left of the dashed red line are visually inspected to select only the extended ones as EELGs. Top left corner: the typical errors (for sources with r < 20 mag) are shown.

Nevertheless, a large fraction of objects remains without W2 photometry. Figure 2 shows that the r − W1 color alone acts as a good discriminant on its own. When we plot r − W1 as a function of r magnitude (Fig. 3), the separation between AGNs and galaxies becomes less clear for high values of r. In addition, there is a lack of spectroscopic data at r > 20 and (r − W1) < 0, which casts doubts on the nature of the objects with these characteristics.

Fig. 3. Color-magnitude diagram used to separate low redshift ELGs from other types of sources, for those objects where the W2 flux was not available (grey dots). All sources with SDSS spectra are shown in colors. The color code for those is the same as the one used in Fig. 2, but in this case objects with available W2 measurements are shown in a lighter shade. The objects selected as candidates to be EELGs are enclosed by the blue lines, while those to the left of the dashed red line are visually inspected to select only the extended ones as EELGs. Top right corner: typical errors (for sources with r < 20 mag) are shown.

We again selected the area in Fig. 3 in which the low-redshift galaxies are clustered and avoided as many contaminants as possible. This results in conditions that must be fulfilled by objects without W2 detection as

$$\begin{array}{cc}\hfill (r-W1)& <0\hfill \\ \hfill (r-W1)& >-1.75\hfill \\ \hfill r& <20\hfill \\ \hfill r& >16.\hfill \end{array}$$(3)

After applying both sets of criteria (from Eq. (2) and (3)), we are left with 1447 candidate EELGs. If we had applied the selection in Eq. (3) to all sources with unWISE data, the final sample would be slightly different: it would be missing 16 objects and would includ 78 new ones. Based on the very clear separation and gap between spectroscopically confirmed QSOs and low redshift galaxies in Fig. 2, we consider that including the g − W2 color adds valuable information despite the small change. In view of the uncertainties in the color and magnitude values, a small fraction of the selected galaxies could fall outside the selection areas (and vice versa). This effect is small in the color cuts, with only ∼10% of the selected galaxies being located less than 1σ away from the dividing line between low-z and quasar-dominated regions in Fig. 2. Only 4% of the selected galaxies fulfil this condition for the r − W1 = 0 line in Fig. 3. More galaxies might swap from the selected to the rejected areas for those with magnitudes close to r = 20 in Fig. 3. This number reaches values close to 20% of the selected sample. Nevertheless, because we do not expect a swift physical change at r = 20 mag, we consider that the contamination induced by these uncertainties is relatively small.

Figures 2 and 3 also show an area with a small population of sources with very negative r − W1 color (marked in the figures with a dashed red line at (r − W1) < 1.75, and in Fig. 3 also with a dashed line at r = 20 mag). Some of the objects are spectroscopically confirmed to be low-redshift galaxies, but most are stars. In order to avoid missing a small but extreme population, we inspected the morphology and SEDs of these galaxies (10 selected in Fig. 2 and 19 selected in Fig. 3). This subsample is heavily contaminated by stars, but still shows some promising candidates. Therefore, we added to the ELG candidate sample the five objects in this subsample that have an extended morphology (1452 candidates up to this point).

There are nevertheless 5668 objects that remain in the parent sample and do not have any counterpart in the unWISE catalog. We removed 60 stars that were identified with a cross-match to Gaia DR2. For this subsample without WISE data, we used a method of removing contaminants considering only the r-band magnitude and the apparent effective radius of the object (R_eff) as measured with the R_EFF parameter in the J-PLUS database. This value is defined as the radius that encloses half of the total flux in the object in the r band, considering the run of the SExtractor software (Bertin & Arnouts 1996) performed in the J-PLUS data reduction pipeline (Cenarro et al. 2019). As shown in Fig. 4, at R_eff > 2.8 arcsec and r < 20 mag, the vast majority of sources with spectra are low-redshift ELGs: all five sources with spectra but without an unWISE counterpart that fulfil these conditions are in fact low-redshift ELGs. Using R_eff to separate contaminants is less efficient than using WISE photometry (therefore we do not use it in the remaining sample), but we still recover 92 additional candidates. At this stage, we have 1544 candidate EELGs.

Fig. 4. Diagram showing the effective radius as a function of r magnitude, used to select candidates to be EELGs for those objects with excess flux in J0515 and no data in the unWISE catalog (grey dots). All sources with SDSS spectra are shown in colors. The color code for those is the same as the one used in Fig. 2, but in this case objects with W1 data are shown in a lighter shade.

While this final selection biases our sample toward less compact EELGs, the effect is very small. Only 45 objects in the parent sample have no WISE data, r < 20 and R_eff < 2.8. Considering the selection of targets with WISE data, we would only expect about half of these objects to be low-redshift galaxies. Therefore, missing these galaxies would mean a very small reduction in the sample of candidates (∼1.4%), while the contamination rate would increase by ∼30 %. Therefore no attempt was made to include objects without WISE detection and R_eff < 2.8.

2.2.5. Cosmic ray removal

In some cases, one of the (usually three) individual frames that were combined to form the final J0515 image of an object was affected by a cosmic ray. This increased its flux significantly. In a few of these instances, the data reduction process failed to remove the contaminated frame from the coadded frames, and the object appears in our selection because it apparently shows a strong excess in J0515. In order to remove these contaminants, we downloaded all frames that contribute to the J0515 image of each candidate, measured the object flux in each, and determined if any individual frame deviated more than 20% from the median flux. In this way, we removed 21 objects, one of which is a spectroscopically confirmed star.

2.2.6. Photometric correction for deblended objects

In our sample, we chose to include objects with a blending flag to reach the highest possible completeness. This allowed us to accurately analyze objects with relatively close companions without being contaminated by their emission because the masks for each object were separated. Nevertheless, this implies that for some extended galaxies, we only selected a small star-forming region and the remaining galaxy was considered as a separate object. This would bias our selection because we would identify as “galaxies” what in fact are simply star-forming regions. In order to avoid this, we recomputed the photometry of all 375 selected objects with a deblending flag. We used SExtractor in dual mode, with a detection in the r band, tailoring the parameters of the code to make them appropriate for the sample of galaxies we analyzed. In this case, we did not perform any correction like the one applied to create the PSFCOR photometry, but because the galaxies we studied in this step are extended a correction like this would play a very minor role. After visual inspection, 171 galaxies were confirmed to be better represented by our SExtractor run rather than the original J-PLUS photometry, and we kept these data during the remaining analysis. During this stage, 30 galaxies were removed from the sample: 20 because they were spurious objects (almost all, spikes from bright stars), one was a repeated object, and in 9 cases, neither the original photometry nor the new SExtractor run were able to properly measure them. For more details, see Appendix B.

2.3. Summary of the sample selection

With all these considerations, we built a sample of 1493 sources that are candidate extreme [OIII] emitters at z < 0.06 by selecting objects with a high J0515 to r ratio. We removed the vast majority of contaminants from other types of sources or higher-redshifts galaxies, using J-PLUS and WISE photometry, and SDSS spectra. Out of 85 objects with available SDSS spectra in the candidate sample, only 3 are not galaxies at z < 0.06: one is a star (likely selected due to an error in the photometric measurement), and 2 are galaxies at z = 0.072 and z = 0.123, with a relatively high uncertainty in the J0515 flux. With these three interlopers, the purity of our candidate selection is ∼96%. Out of the 2560 objects with spectra in the excess flux sample, 89 are galaxies at 0.006 < z < 0.056 with r < 20 mag. Eighty-two of them were selected in the candidate sample, which translates into a completeness of ∼92% in our selection compared to the excess flux sample.

3. Results

To further confirm and characterize the nature of the sample of 1493 candidate EELGs at z < 0.06, we performed SED fitting using the CIGALE code (Noll et al. 2009; Boquien et al. 2016), which yields physical properties for the galaxies. We used the synthetic spectra to select only systems with the highest EW values. Additionally, the galaxy images were inspected to classify their morphologies, and the line fluxes estimated with CIGALE were compared to the available SDSS spectra.

3.1. SED fitting

To analyze the physical properties of the sample of galaxies candidate EELGs in detail, we performed SED fitting with CIGALE. It is a fast and flexible software implementing theoretical models for the different galaxy components. CIGALE creates a large grid of composite stellar populations based on single stellar population models and a variety of SFHs. It also includes models for dust extinction (both for gas and stars), dust emission, and most importantly, nebular emission (both lines and continuum). The resulting grid of models was fit to the photometric data, and the galaxy properties were estimated by analyzing the posterior likelihood distribution, producing a best-fit model, and a Bayesian estimate for each parameter.

3.1.1. CIGALE parameters

We used a simple model of two stellar populations with an exponentially declining SFH: an old population selected to represent the underlying galaxy, and a young one to reproduce the strong starburst causing the extreme emission lines. This type of modeling has been used in the literature, especially for galaxies with recent events of star formation (e.g., Nilsson et al. 2011; Catalán-Torrecilla et al. 2015; López-Sanjuan et al. 2017; Lumbreras-Calle et al. 2019a; Arrabal Haro et al. 2020). The very high EW of the emission lines detected in the galaxies in our sample implies a very short timescale for the event, therefore a very young population is more appropriate than a continuous star formation model. The existence of an old, underlying stellar population that accounts for the majority of the stellar mass has been extensively demonstrated in the literature. Targeting systems similar to ours (e.g., blue compact dwarf or green pea galaxies) and using different methods (resolved stellar populations, spectroscopic analysis, and 2D analysis of multiwavelength photometry), many authors have identified this host component (Östlin 2000; Cairós et al. 2003; Amorín et al. 2007, 2009, 2012; Janowiecki & Salzer 2014; Telles & Melnick 2018; Clarke et al. 2021). Telles & Melnick (2018) and Lopes et al. (2021) have shown that three separate populations produce better fits, but for the purposes of our work and in order to reduce degeneracies, we used only two populations. In order to further simplify the models, we fixed the τ value of the exponential to 50 Myr for the old population and 1 Myr for the young one.

We used Bruzual & Charlot (2003) stellar population models, while the nebular emission used by CIGALE is predicted based on the photoionization models by Inoue (2011). We assumed the Calzetti et al. (2000) extinction law for dust extinction.

In Table 2 we summarize the main free parameters we used to fit the galaxy SEDs, and we present a representative example in Fig. 5. For a more detailed account of the SED fitting process, including the values used for all the different parameters, see Appendix C.

Fig. 5. Example of the J-PLUS SED and CIGALE fit of an EELG. Red triangles correspond to the observed photometry, grey dots to the CIGALE synthetic photometry, and blue squares to the CIGALE synthetic photometry only considering the stellar and nebular continua (not considering the emission lines). The thick lines represent the CIGALE synthetic spectrum, considering different components: the grey line takes into account the full model, while the blue does not take into account the emission lines. The red line represents the stellar continuum, while the brown and green only consider the old and young stellar populations, respectively. The thin grey lines in the background represent the transmission profiles of all 12 J-PLUS filters.

3.1.2. SED results

In Fig. 6 and Table 3 we show the distribution of the main stellar and nebular parameters that CIGALE delivered for our sample of 1493 EELG candidates. We chose to use the best-fit parameters instead of the Bayesian estimates in the CIGALE output for consistency because CIGALE does not provide Bayesian estimates for individual components in the synthetic spectra, and we rely on them to compute EWs. In addition, the Bayesian-estimated photometric redshifts are less accurate than the best-fitting ones (see Sect. 3.4).

Fig. 6. Histograms showing the best-fitting values for the J-PLUS galaxies, derived using the SED fitting software CIGALE. We show in gray the results for the 1493 galaxies in the candidate sample and in red the 466 galaxies in the EELG sample (EW[OIII] > 300 Å). From left to right and top to bottom: we show the age of the young population, the stellar metallicity, the ionization parameter, the color excess, the mass ratio of young and old populations, and the total stellar mass. The gaps observed in some of the histograms are due to the sampling in the parameters (see Table 2 and Appendix C).

It is important to note the limits on the accuracy of these SED fits. They are mainly driven by the goal to obtain accurate estimates of the flux of the most intense emission lines as well as the underlying continuum flux. Additionally, we relied on the stellar mass estimates for comparison purposes, and on other parameters (burst age, extinction, and metallicity) for the broad properties of the sample, given the assumptions in the CIGALE run. The results of the fits cannot be analyzed in extreme detail, especially in some sparsely sampled, highly degenerated parameters such as metallicity, log(U), and escape fraction. For an additional discussion of this topic, see Appendix C.

The age of the star formation burst is very young in view of the best-fit parameters; almost no galaxies have a burst older than 8 Myr. This is consistent with the extreme EW values we measured because that parameter decreases very rapidly in the first few million years of a star formation burst (see, e.g., Leitherer et al. 1999).

The metallicity and ionization parameter of the sample are low. This is typical for low-mass galaxies with strong bursts of star formation, especially at high redshift (Khostovan et al. 2015; Tang et al. 2021; Matthee et al. 2021). Nevertheless, the selection process for our sample, imposing a very high EW of [OIII], may prevent the selection of some extremely metal-poor galaxies (see Sect. 4.1.3 for more details). The typical extinction values that CIGALE derives for the main sample are low, which is consistent with previous results obtained for samples of low-mass star-forming galaxies (Garn & Best 2010; Duarte Puertas et al. 2017; Lumbreras-Calle et al. 2019a).

Even considering the strong burst of star formation in these galaxies, the mass ratio of young to old population is typically low (log$(M_{\mathrm{you}}/M_{\mathrm{old})}=-2.00\begin{array}{c}+0.70\\ -0.60\end{array}$). This is not surprising, however, because the mass-luminosity ratios vary strongly between old and very young populations. Therefore, while the old population dominates the total mass of the galaxy, a relatively low-mass) young population can have strong effects on the integrated photometry of a low-mass galaxy. The total stellar masses of the galaxies in the sample are not constrained by input parameters, and they span most of the range of what is usually considered dwarf galaxies, with log(M_⋆/M_⊙)∼6.5 − 9.5.

Two of the main parameters in the analysis of ELGs are the flux and EW of the emission lines (in our case, especially the [OIII] line, which we used for the sample selection). We focused on the J0515 filter, in which most of the line flux comes from the [OIII]4959+[OIII]5007 lines. To compute them, we first convolved the synthetic stellar and nebular continuum derived by CIGALE (blue line in Fig. 5) through the transmission curve of the J0515 filter to obtain an estimate of the continuum in this filter (F_cont., blue squares in Fig. 5). Then, using the measured flux in that filter (F_J0515, red triangles in Fig. 5) and the classical formula (assuming a flat continuum within the filter, and considering fluxes in units of erg/s cm⁻²/Å), we obtained the EW as

$$\begin{array}{c}\hfill \mathrm{EW}[\mathrm{\AA }]=\frac{F_{J0515}-F_{\mathrm{cont}.}}{F_{\mathrm{cont}.}}\mathrm{\Delta }=\frac{F_{\mathrm{em}}}{F_{\mathrm{cont}.}}\mathrm{\Delta },\end{array}$$(4)

where Δ is the width of the filter, 200 Å, and ΔF_em is the line flux. In Fig. 7 we show the histogram of the EW for the J0515 filter with the galaxies with S/N > 3 in the J0515 EW in red.

Fig. 7. Histogram of EW([OIII]), measured using the J-PLUS J0515 photometry and the CIGALE fits for the continuum. The filled white histogram represents the whole sample of candidates, and the filled red histogram shows only the galaxies with S/N > 3 in the EW([OIII]) value. The blue vertical line represents the minimum EW threshold for selecting EELGs in this work, 300 Å.

3.1.3. Comparison between CIGALE and J-PLUS fluxes

In the redshift range we covered, 0.0075 < z < 0.06, the J0515 filter for most of the sample is affected by the [OIII]5007 and [OIII]4959 emission lines. According to the photometric redshifts obtained from CIGALE, in a small fraction of galaxies (∼18%, at the higher redshift end), the [OIII]5007 line lies in a low transmission region of the J0515 filter (or even entirely out of it). The Hβ line enters the filter wavelength range in almost half of the sample, while at the very low redshift end, some galaxies (∼5%) lack the contribution of the [OIII]4959 line. Nevertheless, the limited photometric redshift precision (see Sect. 3.4) prevents us from providing a detailed account of this distribution of emission lines, galaxy by galaxy. We can still study the whole sample by comparing the emission line flux in the whole filter, estimated using the J0515 photometry and the CIGALE-derived continuum value, with the fluxes computed by CIGALE for each emission line. The same can be done for the r broadband filter, which is significantly contaminated by the Hα emission line. In Figs. 8 and 9 we show the logarithm of line fluxes derived directly from the J-PLUS photometry as a function of the CIGALE-estimated individual line fluxes. The integrated fluxes are clearly dominated by the brightest emission lines in several filters. In the r filter (Fig. 9), taking the Hα line alone into account is enough to recover a very good one-to-one relation, with a small offset (∼ − 0.04 dex) and low scatter (1σ ∼ 0.13). If we include the [NII]6584 flux in the analysis, the offset is greatly reduced (down to −0.008 dex), with a similarly low scatter (1σ ∼ 0.12 dex). Nevertheless, given the small impact of [NII]6584, we consider that the line emission in the r filter is mostly dominated by Hα. For the J0515 filter (Fig. 8), taking the [OIII]4959 and [OIII]5007 lines for all galaxies into account provides a reasonably good fit, with a negligible offset (∼ − 0.007 dex) and and very low scatter (∼0.096 dex). Nevertheless, in the highest redshift range of our sample (z ≥ 0.05), the [OIII]5007 line falls in a very low transmission region of the J0515 filter. This translates into the parallel subsample, ∼0.35 dex below the one-to-one relation, clearly noticeable in Fig. 8, and only ∼7.5% of the sample are affected. The contamination by the Hβ line for galaxies at 0.038 < z < 0.045 results in a small offset in the opposite direction. Because of the uncertainty in the photometric redshift estimates, we cannot correct for this offsets. In Fig. 8, a deviation from the one-to-one line can be seen at low CIGALE fluxes. The very few sources with large deviations are caused by issues in the modeling process (e.g., estimating a high Hβ emission, or failing to reproduce the observed J0515 flux entirely). The lack of sources with large scatter below the one-to-one line is caused by the limit in measured J-PLUS J0515 line flux, at around −14.1 log(erg s⁻¹ cm⁻²). The limit is a byproduct of the EW and S/N limits, as well as of the r-magnitude cut. This results in a scatter only above the identity line because the sources below it were rejected. Nevertheless, this effect is small, and it only affects a few sources at the low S/N limit.

Fig. 8. Comparison between the emission line flux measured in the J0515 filter and the [OIII]5007+4959 flux estimated by the CIGALE SED fit for the EELG sample. The black line represents the one-to-one relation.

In conclusion, considering the generally good agreement shown in Figs. 8 and 9, we estimate that the J0515 EW and line fluxes correspond to the sum of the [OIII]4959 and [OIII]5007 emission lines, and the r line flux and EW correspond to Hα (see Table 4). In addition, for the 18 galaxies with photometric redshift z_phot < 0.017, the Hα line falls within the J0660 filter wavelength range, and the precision of the Hα measurement is much higher, with an offset of ∼0.002 dex and a scatter of ∼0.02 dex.

Fig. 9. Comparison between the emission line flux measured in the r filter and the Hα flux estimated by the CIGALE SED fit for the EELG sample. The black line represents the one-to-one relation.

3.2. Selection of the EELG sample

Our original selection of candidates was performed using a very rough estimate of the stellar continuum below the [OIII] line (simply the r-band flux), and this results in an extended [OIII] EW distribution (see Fig. 7). We decided to perform an additional cut to obtain a clearly defined EELG sample, with a strict limit in EW. We considered (only for galaxies with S/N > 3 in [OIII] EW) the variation in the logarithm of the number of galaxies in each EW bin with EW given the expected decreasing exponential relation (see Fig. 10). This exponential relation is clear (with R² = 0.998), for example, when selecting galaxies with EW([OIII]) > 10 Å and S/N > 3 in the SDSS spectroscopic database (the GalSpecLine table from DR8, Aihara et al. 2011). In our case, for a given lower limit of EW, the number of galaxies does not increase as fast as an exponential: that would be our completeness limit. In order to compute it accurately, we performed linear fits with EW bins of 0.1 dex in log(Å), limiting the analysis in the high EW range to the bins showing more than ten galaxies, which translates into log(EW) < 3.05 log(Å). For the low EW limit, we performed three separate linear fits, covering the bins from 2.45, 2.55, and 2.65 log(EW). For log(EW) values lower than 2.5 (EW ∼ 300 Å), the measured value drops below that of a linear fit to the bins above it (Fig. 10), with the fits limited to 2.55 and 2.65 showing comparable parameters. Therefore we placed our threshold at EW = 300 Å in [OIII] to keep a complete sample, obtaining a total of 466 EELGs. Because several emission lines can fall within the J0515 filter wavelength range depending on the redshift of the source (see the discussion in Sect. 3.1.3), the 300 Å limit may result in moderately different cuts for different sources. Nevertheless, considering the good agreement presented in Fig. 8, we are confident that for most of the sample, it represents the [OIII]4959+5007 EW accurately. Finally, in the higher EW limit, completeness is limited by the r < 20 threshold defined in Sect. 2 because fainter (and lower mass) galaxies tend to show higher EWs (see Sect. 4.2.1). This magnitude cut is nevertheless necessary in order to provide a clearly defined sample without significant contamination of high-z AGNs (as seen in Fig. 3).

Fig. 10. Diagram used to asses the completeness of the EELG sample. We plot the number of galaxies per bin of log(EW[OIII]) as a function of log(EW[OIII]). We overplot the linear fits to the data considering three different limits in the lower end, from 2.45 log(Å) (black line), 2.55 (red), and 2.65 (green). The points considered in at least one of the fits are shown in a darker shade of gray.

The value we selected is higher than some limits found in the literature for starburst galaxies or EELGs (100 Å in Amorín et al. 2015, 80 Å in Hinojosa-Goñi et al. 2016), but it is lower than most of the galaxies in other surveys (Cardamone et al. 2009; Yang et al. 2017), especially those at high redshift (van der Wel et al. 2011; Atek et al. 2011). We consider that the intermediate value of 300 Å is appropriate for the width of the filter we used to select the sample. It provides a large sample, showing undoubtedly strong emitters, allowing both a statistical analysis and low contamination. A higher threshold would have removed interesting objects and restricted the sample, while a lower threshold could provide a high rate of contaminants unless we removed the lower brightness objects, which would be large fraction or the sample. By comparing with the SDSS spectroscopic database (using again the GalSpecLine table from DR8), we see that with a 300 Å threshold, we keep only the top 2% of the EW([OIII]) distribution (considering only galaxies with EW([OIII]) > 10 Å), which further confirms that the emission properties of our sample are extreme.

We also considered the effect of the initial selection (Eq. (1)) in the Hα flux distribution of the final sample. Since the Hα lines falls within the wavelength range of the r filter and Eq. (1) selects objects with high J0515 to r ratio, we could be biasing the sample, removing objects with relatively high [OIII] EW yet low [OIII]5007/Hα ratio. To test if this effect was present, we artificially increased the CIGALE-estimated Hα fluxes in the selected galaxies (those with EW([OIII]) > 300 Å) and check if they still fulfilled Eq. (1). A 25% increase in the Hα flux was used, since that produced a significant offset in the [OIII]5007/Hα ratio as a function of EW([OIII]), unobserved in the spectroscopic data. Using that 25% increase, only 16 objects (3.2%) were rejected, therefore we are confident that no significant bias is added to our sample against galaxies with low [OIII]5007/Hα ratio due to Eq. (1).

This sample selection results in 466 EELGs. In order to determine how novel the J-PLUS database and the present work are, we investigated the number of objects in our sample that were not previously identified as EELGs. To do this, we first queried the NASA Extragalactic Database (NED) to select cataloged objects near the positions of the galaxies in our sample. Only 53 have a reported spectroscopic redshift, and we considered that these galaxies have been already classified as EELGs. In addition, 14 galaxies without spectroscopic redshift in NED have been referenced in at least one publication, and 2 of them were described as EELGs or attributed a similar category. Finally, we confirmed that all objects with spectroscopic confirmation in this EELG sample are indeed low-redshift galaxies, which translated into a theoretical purity of 100%. In conclusion, we here identify 411 galaxies (88% of the sample) as EELGs that were previously unknown to belong to this class. This highlights that the J-PLUS survey is well suited to finding new extreme emitters because it reaches fainter magnitudes than spectroscopic surveys that cover wide areas of the sky.

3.3. Morphology of the EELG sample

We performed a simple visual classification of the morphology of the extreme [OIII] emitters. The EELGs in the sample were classified into three categories: compact (c), semicompact (s), and extended (e). Compact galaxies show a circular shape, without further discernible structure. We classified them into the semicompact category when their shape is dominated by a circular and bright clump, but some fainter structure is revealed: a tail, a halo, and so on. We considered that extended galaxies show a complex morphology, have no bright clump, or clump or clumps do not dominate the light profile in the galaxy. As an example, we show the RGB postage stamps of four galaxies in Fig. 11.

Fig. 11. RGB composites (made using images taken with z, r, and g filters) of four galaxies in the EELG sample. Bottom right corner of each panel, a letter indicates the type of morphology: compact (c), semicompact (s), or extended (e). Bottom left corner of each panel: 10 arcsecond line for scale. All images have been obtained using the Legacy Survey viewer with data from the BASS and MzLS surveys (Zou et al. 2019; Dey et al. 2019).

For this analysis, we took the images from the Legacy Survey database², which were obtained in the framework of the BASS and MzLS surveys (Zou et al. 2019; Dey et al. 2019). These surveys provide images in the g, r, and z bands, with median depths of 23.65, 23.08, and 22.60 mag respectively (considering 5σ detections of point sources), and typical full width at half maximum (FWHM) values of their PSFs of 1.61″, 1.47″, and 1.01″(Dey et al. 2019). These images reach low surface brightness detection limits in r of 27.9 mag arcsec² for 3σ detection of a 100 arcsec² feature, which is fainter than the SDSS and Panstarrs PS1 surveys, and comparable to SDSS Stripe-82 (Hood et al. 2018). Therefore, even if fainter features than those visually noticeable in the RGB composite images may exist (and change the morphological classification), we consider that our analysis provides results that are sufficiently accurate for this work.

The compact class is the most common and covers 43% of the EELG sample, followed by the semicompact class (38%). The remaining galaxies are classified as extended. This is consistent with previous studies of this topic (e.g., Izotov et al. 2011; Yang et al. 2017), which focused on compact or semicompact galaxies. Nevertheless, we show that a significant amount (19%) of the EELGs in our sample are more properly classified as extended systems, suggesting possible extensions of previous EELG searches.

3.4. Comparison with spectra

We assessed the accuracy of the properties derived using J-PLUS photometry and the CIGALE SED analysis by discussing the results obtained for the 82 sources in the candidate sample with SDSS spectra. For this analysis, we removed the misclassified objects (one star and two high-redshift galaxies).

First, we measured the accuracy of the photometric redshifts determined with CIGALE, which are compared in Fig. 12 with the spectroscopic redshifts. The agreement is very good. We computed a parameter to quantify the quality of the agreement, the σNMAD, defined as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{NMAD}}=1.48\times \mathrm{median}|\mathrm{\Delta }z-\mathrm{median}(\mathrm{\Delta }z)|,\end{array}$$(5)

Fig. 12. CIGALE photometric redshift as a function of SDSS-derived spectroscopic redshift for the galaxies in our sample of candidates with available spectra (excluding the three misclassified objects). The line represents the one-to-one relation. In the top left corner of the figure, we show the minimum step in photometric redshift considered with CIGALE (0.0025) as error bars.

where Δz = (z_best − z_spec)/(1+z_spec). We obtained a value of σNMAD ∼ 0.003, similar to the mini-JPAS results, with 56 narrow-band filters (Hernán-Caballero et al. 2021). This high accuracy is likely due in part to the limited redshift range covered, to the extreme intensity of the emission lines, but most importantly, to the fact that these lines ([OII]3727, [OIII], and Hα) lie inside or near the wavelength range of narrow or medium-band filters, which produces very clear photometric features. For more than half of this subsample (56%), the difference between photometric and spectroscopic redshift is smaller than the minimum step used in the photometric redshift determination (0.0025). This high accuracy diminishes the possible uncertainties in our determinations of absolute parameters, such as luminosities and stellar masses. We used the best-fitting CIGALE models instead of the Bayesian estimates, in part because the agreement with spectroscopic redshift values drops when using the latter. In consequence, the values discussed for the SED fitting (ages, masses, etc.) correspond to the best-fit models.

We also compared the fluxes and EWs estimated in the J-PLUS photometry with those measured using the SDSS spectra. We used our own code to measure the emission line fluxes and EWs in the spectra, a revised version of the procedure used in Lumbreras-Calle et al. (2019a). Briefly, we computed the continuum below each emission line by masking the line and performing a linear fit to the remaining spectrum in a 200 Å aperture around the line. We estimated the values and uncertainties of the linear fits by performing bootstrap simulations. We computed the flux of each emission line by simply adding the measured flux over the spectral window considered and subtracting the continuum value (and then dividing by the continuum value to obtain the EW).

It is important to recall that the photometry we used is the PSFCOR J-PLUS photometry, re-scaled to the AUTO aperture value in r. In order to perform an accurate comparison, we need to use the 3ARCSEC aperture in the J-PLUS catalog, which should match the SDSS spectroscopic measurement (which is performed with a 3 arcsec wide fiber in most of the sample) more closely. Nevertheless, some sources of discrepancy still exist: different seeing conditions, small differences between the sky position of the objects, or a possible offset between absolute calibrations. We compared the 3ARCSEC r fluxes in J-PLUS with the synthetic r fluxes computed from the spectra, and found very good agreement (a 1σ = 0.12 dex scatter) with a small but noticeable offset (0.08 dex). Therefore, we rescaled the 3ARCSEC photometry to the SDSS flux using this r -band offset.

For a precise comparison, line by line, we need to use the CIGALE model output, which provides us with fluxes for several emission lines ([NII]6584, Hα, Hβ, [OIII]5007, [OIII]4959, and [OII]3727). In Figs. 13 and 14 we show the comparison between spectra and SED models for the [OIII]5007 and Hα lines, respectively. For [OIII]5007, the agreement between the measurements is very good, with almost no bias (a median difference of −0.001 dex) and very low scatter, a 1σ value of ∼0.18. We can also compare our result with others in the literature. In the Census of the Local Universe (CLU) preliminary fields, Cook et al. (2019) showed the comparison between photometric and spectroscopic Hα fluxes in their Fig. 10. Considering only their 5σ detections, we estimate a 1σ scatter of ∼ 0.25. This value is higher than ours, which is striking considering that their filters are notably narrower than ours (from 76 to 92 Å compared to 200 Å). For a more similar comparison, using J-PLUS data, we examined the results in Logroño-García et al. (2019). Considering Hα fluxes measured in the J0660 filter, these authors reached more accurate values than we did, with a 1 σ scatter of ∼0.11. This result is expected because the J0660 is narrower than J0515 (∼140 versus ∼200 Å), and because they took extreme care in matching the apertures of the photometric and spectroscopic measurements, which included integral field unit data.

Fig. 13. [OIII]5007 flux estimated using CIGALE SED fitting of J-PLUS data as a function of [OIII]5007 flux measured in SDSS spectra. Same sample as in Fig. 12. The black line represents the one-to-one relation.

Fig. 14. Hα flux estimated using CIGALE SED fitting of J-PLUS data as a function of Hα flux measured in SDSS spectra. Same sample as in Fig. 12. The black line represents the one-to-one relation.

The comparison for the Hα line is also very successful, with a similar offset (∼0.04 dex) and scatter (∼0.18 dec.) even if the line flux is estimated using the broadband r filter. This value is likely low because the data in the J0660 filter very accurately trace the continuum near the Hα line for galaxies at z > 0.017. In Figs. 13 and 14, a small offset can be seen at low spectroscopic fluxes, which is caused by a similar effect as described in Sect. 3.1.3 for Fig. 8.

For the few galaxies in our sample with available spectra and z < 0.017, CIGALE used data from the J0660 filter to estimate the Hα flux and obtained an even better agreement (with an offset of ∼ − 0.01 and a scatter of ∼0.12). Other emission lines show a worse agreement between spectra and CIGALE models, some simply because they are less intense (e.g., the [NII]6584 line, which in addition is very close to the much more intense Hα line). Some others are not only faint, but also located in regions in which the spectroscopic and photometric analysis is more complex. For example, the [OII]3727 line lies outside the wavelength range of most of the spectra (and at the edge of some), and it can fall in the gap between two blue narrow-band J-PLUS filters (J0378 and J0395), which are also harder to calibrate than the redder filters. Nevertheless, it is worth noting that making direct measurements on the J-PLUS data, the [OIII]/[OII] ratio of the sample of EELGs reaches extremely high values compared to typical star-forming galaxies. This indicates very hard ionizing radiation and is similar to very high redshift galaxies. We will explore this result in depth in an upcoming paper, in which we analyze the spectra of some of the galaxies in the current sample.

The previous comparison between line fluxes, even if very successful for [OIII]5007 and Hα, relies heavily on the SED models. For a less model-dependent version, we compared the EW measured directly in the J-PLUS photometry and the synthetic EW values for each filter, obtained by convolving the emission lines in the SDSS spectra with the filter transmission curves. We should note that some modeling is still required because we used the CIGALE fits in order to estimate the continuum in each filter. We show the results for J0515 and r in Figs. 15 and 16, respectively. In this case, the agreement is even better, with similarly small offsets (0.02 for J0515 and 0.04 for r), but lower scatter values (0.12 and 0.17, respectively). This result strengthens our confidence in the accuracy of J-PLUS photometry and in our SED and EW analysis.

Fig. 15. J0515 EW directly measured on J-PLUS data as a function of J0515 EW estimated convolving SDSS spectra with the filter transmission. Same sample as in Fig. 12. The black line represents the one-to-one relation.

Fig. 16. r EW directly measured on J-PLUS data as a function of r EW estimated convolving SDSS spectra with the filter transmission. Same sample as in Fig. 12. The black line represents the one-to-one relation.

4. Discussion

In this section we compare our sample with different works in the literature regarding their broadband colors, number density values, and the ratio of [OIII] and the Hα and Hβ emission lines. In addition, we discuss the physical properties derived for our EELG sample, considering EW, star formation rate (SFR), mass, and other properties in the context of the literature at different redshifts.

Because the photometric and spectroscopic redshifts and the photometric and spectroscopic line fluxes and EWs agree well, we are able to discuss the results of this work without performing follow-up spectroscopic observations of the whole sample of galaxies. This step has been necessary in previous works, however, that used only broadband selection (Yang et al. 2017; Senchyna & Stark 2019; Kojima et al. 2020) or narrow-band selection like Hα dots (Kellar et al. 2012; Salzer et al. 2020) and the CLU survey (Cook et al. 2019). Even if the physical information we can obtain is limited compared to what can be derived from spectroscopic observations, it is still significant, considering that the data span 2176 sq. deg. down to r = 20. In addition, this photometric analysis allows for a more efficient spectroscopic follow-up that targets galaxies with specific physical properties. These characteristics will improve significantly in the upcoming J-PAS survey, with deeper observations and higher spectral resolution.

4.1. Testing the sample selection

4.1.1. Comparison with broadband color-color selection

Several works over the past years have used broadband colors from large photometric surveys (most notably SDSS) in order to select galaxies with strong emission lines and/or extremely metal-poor gas (Cardamone et al. 2009; Yang et al. 2017; Senchyna & Stark 2019; Kojima et al. 2020). For the selection of extreme [OIII] emitters with broadband data, it is necessary to define regions of extreme g − r color (or r − i at higher redshift) in order to avoid selecting typical galaxies. Line ratios and redshift affect these colors, preventing some galaxies with very high [OIII] emission from showing extreme colors. This is presented in Fig. 17, where many extreme [OIII] emitters in our sample are indistinguishable from the main SDSS galaxy population using only the g, r, and i bands. Our sample of EELGs covers the region of the Yang et al. (2017) blueberry galaxies (since they present extreme [OIII] emission at our redshift range), as expected, but it also covers the same color space as extremely metal-poor galaxies both from observations and models (Kojima et al. 2020). This effect is mainly due to the different Hα/[OIII] line ratios, given how dominated the broadband fluxes are by the line emission in this type of galaxies. If this ratio is high, the g − r color is not extreme, while if the ratio is low, the g − r color is highly negative. This ratio, even if high, does not prevent us from identifying these extreme emitters because we study the flux in the J0515 medium-band filter. This again proves the added value of the medium-band filters in the J-PLUS survey.

Fig. 17. Broadband color-color diagram, created using SDSS data. We plot our EELG sample in blue dots, alongside the blueberry galaxies in Yang et al. (2017) (black triangles) and samples of extremely metal poor galaxies (XMP) from Kojima et al. (2020) (red squares) and the other literature works (green diamonds), as reported in (Kojima et al. 2020). The blue line represent the Yang et al. (2017) sample selection limits, while the green lines follow the evolutionary tracks of models in Kojima et al. (2020). The small black dots represent the density contours of typical SDSS galaxies (Kojima et al. 2020).

4.1.2. Number density of EELGs

The detection of EELGs presented in this work follows clear and reproducible procedures, without a preselection of targets (other than a magnitude limit). We can therefore use our sample to estimate the number density of this class of galaxies in the redshift range we probed, and compare it with other works in the literature. In order to roughly estimate uncertainties, we considered two sources of errors: Poisson noise (associated with the discreet nature of galaxies) and cosmic variance. We approximated the Poisson noise by $\sqrt{N}$ (where N is the number of galaxies) except when N < 100, where we used the asymmetric prescriptions in Gehrels (1986). We estimated the error associated with cosmic variance using the code presented in Driver & Robotham (2010). We added in quadrature both errors to obtain an uncertainty estimation. The results of this analysis are shown in Table 5.

First, the number density of EELGs in our J-PLUS sample was computed. We limited the analysis of our sample to the redshift range in which the two [OIII] lines lie within the J0515 filter, 0.016 < z < 0.048, in order to ensure completeness. We find 394 EELGs in this range, which translates into a number density of (2.45 ± 0.27)×10⁻⁴ Mpc⁻³. This is about one EELG every 4000 Mpc³.

In order to compare the number density of EELGs we measure with other works, we imposed limitations in magnitude (or mass), [OIII] EW, and redshift range. This was done for the literature samples and our sample in order to match them appropriately. We started the comparisons with the GAMA survey (DR3 Baldry et al. 2018), which provides precise spectroscopic data while being still relatively deep and complete. They obtained fiber spectra of essentially all sources brighter than a certain magnitude threshold in the targeted fields. In order to perform a density comparison, we limited ourselves to three of the four fields with available data (excluding G02 because they used a different input catalog). To further homogenize the selection, we limited the comparison to galaxies brighter than Petrosian magnitude 19 in the r band, where GAMA and our sample are complete. Using the data in the GaussFitSimple table within the SpecLineSFR data management unit, we selected galaxies with EW([OIII]5007) > 225 Å in GAMA, which would correspond to ∼300 Å in our J-PLUS analysis ([OIII]4959+5007). This yields 13 GAMA galaxies over the 0.016 < z < 0.048 redshift range. A similar cut in Petrosian magnitude and EW in our sample yields 113 objects. When the relative areas and the volume of the Universe between these redshifts are taken into account, the density of GAMA sources is $(8.{77}_{-3.37}^{+3.96})\times {10}^{-5}$ Mpc⁻³ and (6.96 ± 0.96)×10⁻⁵ Mpc⁻³ for J-PLUS.

We also compared to the SDSS database, limiting our analysis in this section to the main legacy survey and the MPA-JHU catalog (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004; Salim et al. 2007), presented in the GalSpecLine table from DR8 (Aihara et al. 2011). The clear selection of this sample (essentially Petrosian magnitude brighter than 17.7 in r; Strauss et al. 2002) allows for an accurate comparison. We chose galaxies in the same redshift and EW ranges as with the GAMA survey. Their density in this case is $(1.{37}_{-0.23}^{+0.24})\times {10}^{-5}$ Mpc⁻³, compatible with the $(0.{99}_{-0.27}^{+0.33})\times {10}^{-5}$ Mpc⁻³ density in our J-PLUS data for the same cuts.

The results from the previous two comparisons show that the density of EELGs computed in this work using J-PLUS is compatible with the densities derived from magnitude-limited spectroscopic surveys, which can be considered as ground truth. This result shows the strength of our contribution, which is able to select essentially all EELGs at this redshift range down to a certain magnitude threshold. Moreover, our work has some advantages over spectroscopic analyses: we cover much wider areas than pencil-beam surveys such as GAMA, and we are complete down to deeper magnitudes than a wide-field spectroscopic survey such as the SDSS.

We can also compare our density values with broadband-selected samples such as the blueberry galaxies (Yang et al. 2017) and green peas (Cardamone et al. 2009). In this case, we restricted our sample to only the compact and semicompact categories to match their sample selection better. In both works, the authors selected only galaxies with very high EW, but without a clear threshold value. Therefore we limited our comparison to the galaxies with the highest EW values of these samples, in order to approach completeness. For Yang et al. (2017), we placed the EW limit roughly where the number of galaxies per EW bin decreases with increasing EW of the [OIII]5007 emission line (around 1200 Å, corresponding to 1600 Å in [OIII]5007+4959). In addition, we cut the comparison sample at r < 20 to simulate our limit in brightness, given the similar redshift range covered in both works. Estimating these limits for Cardamone et al. (2009) is more complex because of their EW distribution and higher redshift range (which results in higher masses). As an exercise, we place the limits in EW > 300 Å and log (M_⋆/M_⊙) > 9.5, obtaining only two galaxies in our sample, against 27 in their case. These thresholds result in a density of $(0.{56}_{-0.22}^{+0.31})\times {10}^{-6}$ Mpc⁻³ for Yang et al. (2017), and more than ten times higher in J-PLUS, $(8.{38}_{-0.24}^{+3.01})\times {10}^{-6}$ Mpc⁻³. The comparison between the Cardamone et al. (2009) sample and our work is even more striking, with values of $(0.{020}_{-0.004}^{+0.005})\times {10}^{-6}$ Mpc⁻³ vs. $(1.{24}_{-0.81}^{+1.64})\times {10}^{-6}$ Mpc⁻³, but more uncertain.

Other surveys use narrow-band imaging to identify ELGs. In principle, they should be more sensitive to low EW emission lines than J-PLUS because their filters are narrower than our J0515 filter, and therefore the contrast in brightness between narrow- and broadband should be higher. Nevertheless, because we limited our analysis only to the most extreme events of star formation, the strong contrast is enough to avoid missing a significant number of EELGs. The sky area and redshift range that narrow-band surveys cover is smaller than ours, which limits their ability to identify extreme, rare objects. To our knowledge, no narrow-band survey has targeted the [OIII] emission line in the local Universe, so we cannot perform any direct comparison. A similar example are the Hα dots identified in the ALFALFA Hα survey (Kellar et al. 2012; Salzer et al. 2020). While they selected galaxies based on Hα emission at our redshift range, they provided spectroscopic follow-up and thus [OIII] EW. We computed the density of Hα dots showing [OIII] EW higher than 300 Å and brighter than r = 20 located at z < 0.024 (beyond this, the transmission of their reddest Hα filter drops). The value we obtain for their number density (1.60${}_{-1.49}^{+1.65}\times {10}^{-3}$ Mpc⁻³) is higher than the corresponding value for our sample (0.70 ± 0.10 × 10⁻³ Mpc⁻³), again considering only compact or semicompact galaxies in our J-PLUS comparison. Nevertheless, because of the high uncertainty in the density value derived for the Hα dots survey (∼80%, driven mostly by cosmic variance), the density values are compatible within one standard deviation.

4.1.3. [OIII]5007/Hα and [OIII]5007/Hβ ratios

We intend to provide a sample of EELGs selected purely on magnitude and EW([OIII]), without direct bias on line ratios, because it is only based on the detection of flux excess in a medium-band filter. In contrast with this, other works based on broadband photometry detection (e.g., Yang et al. 2017 and Cardamone et al. 2009) are biased against low [OIII]5007/Hα systems (and therefore also against low [OIII]5007/Hβ systems). This is because they demand that the flux of the filter where the [OIII] doublet falls (g in Yang et al. 2017, r in Cardamone et al. 2009) must be significantly stronger than the filter where the Hα line flux contribution lies (r in Yang et al. 2017, i or z in Cardamone et al. 2009). According to stellar population synthesis and nebular photoionization models (e.g., Inoue 2011), the [OIII]5007/Hα ratio reaches a maximum value for gas metallicities around 12+log(O/H) ∼ 8.0, which is further confirmed by the typical and metallicities of Yang et al. (2017) and Cardamone et al. (2009). Therefore, selecting galaxies with high [OIII]5007/Hα ratios biases the samples against very low metallicities (Senchyna & Stark 2019; Kojima et al. 2020).

Because we only select extreme [OIII] emitters, we are certainly missing systems that have [OIII] EW below our threshold, but above it, we should recover all galaxies, regardless of their [OIII]5007/Hα ratio. In Fig. 18 we compare the distribution of the [OIII]5007/Hα and [OIII]5007/Hβ ratios for several samples of ELGs. We consider the [OIII]5007/Hβ ratios more reliable because the emission lines are much closer in wavelength, which minimizes the effect of dust extinction correction and/or flux calibration issues. Our work shows lower median values for these ratios than broadband selections such as Yang et al. (2017) or Cardamone et al. (2009), and we reach much lower values. Our selection shows slightly higher values than a pure spectroscopic selection (Amorín et al. 2015), especially in the [OIII]5007/Hα ratio. A stronger difference is present than in Hα -selected samples (the Hα dots, Kellar et al. 2012; Salzer et al. 2020) and extremely metal-poor galaxies (Hsyu et al. 2018; Guseva et al. 2017), especially if we include systems with low Hβ EW (a proxy for low sSFR). Nevertheless, the range of values that our survey covers overlaps with most of the range of the extremely metal-poor and strongly star-forming samples, in contrast with Yang et al. (2017) and Cardamone et al. (2009). In addition, we also plot the typical values for other types of systems, some with strong ionizing spectra: HeII emitters (Pérez-Montero et al. 2020) and Lyman α emitters (Matthee et al. 2021; Izotov et al. 2021). In these cases, their ratios are more similar to ours.

Fig. 18. Distribution of the ratios [OIII]5007/Hα (filled red rectangles and dots) and [OIII]5007/Hβ (filled white rectangles and gray dots). The thick black line represents the median value, and the box ranges from the first to the third quartile. The error bars represent the maximum and minimum value without outliers, which are plotted as dots (and defined as the values that lie beyond 1.5 times the interquartilic range). We plot in this figure results from our work and several samples from the literature (Yang et al. 2017; Cardamone et al. 2009; Amorín et al. 2015; Onodera et al. 2020; Hsyu et al. 2018; Guseva et al. 2017; Pérez-Montero et al. 2020; Matthee et al. 2021; Izotov et al. 2021), as well as the KISS survey (Wegner et al. 2003; Gronwall et al. 2004; Jangren et al. 2005; Salzer et al. 2005) and the Hα dots class (Kellar et al. 2012; Salzer et al. 2020). We selected the subsample with extremely low metallicity in Hsyu et al. (2018) (H18_XMP), as well as galaxies with EW(Hβ) > 100 Å in both Hsyu et al. (2018) XMPs (H18_XMP Hβ > 100 Å) and the KISS survey (KISS Hβ > 100 Å). Different samples are separated by vertical blue lines.

Finally, it is also important to keep in mind that several properties in the samples affect the [OIII]5007/Hα ratios, in addition to selection biases and metallicity. The SFR is also positively correlated with the ratio (e.g., when restricting the Hsyu et al. (2018) sample to objects with high Hβ EW). Stellar mass also plays a role (higher-mass galaxies show lower ratios), but it is less relevant than the other factors: the Cardamone et al. (2009) sample has a typical mass ∼ 100 times higher than those of Hsyu et al. (2018) or Guseva et al. (2017), but their ratios are higher. While these factors play a role in the results shown in this section, our selection method is clearly open to the selection of low [OIII]/Hα and Hβ galaxies, in contrast to others.

4.2. Physical properties of the EELG sample

4.2.1. [OIII] EW and stellar mass

A clear negative correlation between emission line EW and stellar mass has been found in the literature (e.g., Fumagalli et al. 2012; Sobral et al. 2014; Khostovan et al. 2016; Reddy et al. 2018; Lumbreras-Calle et al. in prep.). This indicates that galaxies with lower masses tend to have stronger recent star formation events relative to their mass (higher sSFR). The relation between the two variables was fit with a linear model (in logarithmic units) and depends on the emission line studied (both intercept and slope) and the redshift of the galaxies (mostly the intercept), although discrepancies are still present (i.e., between Khostovan et al. 2016 and Reddy et al. 2018, but see Khostovan et al. 2021). It is beyond the scope of the present paper to shed light onto the values of that linear relation because our selection process focuses only on a region of the EW – mass diagram, making a linear fit to our data meaningless in this context. Nevertheless, we can compare the values obtained to those available in the literature for the [OIII] line (see Fig. 19).

Fig. 19. [OIII] EW of the EELG (color dots) and candidate (gray dots) samples as a function of stellar mass. We plot the galaxies in the EELG sample according to their morphology in different colors: red corresponds to compact, green to semicompact, and blue to extended objects. In addition, we overplot linear relations from Khostovan et al. (2016) at different redshifts (black at z ∼ 0.84, blue at z ∼ 1.42, pink at z ∼ 2, and orange at z ∼ 3.24). The linear relations are plotted in as continuous lines in the mass ranges in which they were defined, and as dashed lines in their extensions to lower masses.

The clearest result is that even if our galaxies reside in the very low redshift Universe, they almost exclusively populate the regions of the diagram defined by the linear fits to galaxies at intermediate redshift (z > 0.84), with EW values similar to those typically seen at z ∼ 1.4 (Khostovan et al. 2016). Our sample likewise contains some galaxies with an EW that is similar to those at medium-high redshift: z ∼ 2.2 − 3.4 for Khostovan et al. (2016), z ∼ 3.4 for Reddy et al. (2018), and even z ∼ 8 (De Barros et al. 2019). The very high EW values, low [OII]/[OIII] ratios, and typically compact morphologies are similar to those found in very high-redshift galaxies (e.g., Onodera et al. 2020). According to detailed spectroscopic studies in the low-redshift universe (Izotov et al. 2021), this type of galaxies shares physical properties with those that form stars at very high redshift, likely leaking Lyman continuum radiation, and playing an important role in the reionization of the Universe. In consequence, our sample may be useful in improving the understanding of the very early Universe.

We split the sample into the three morphological categories described in Sect. 3.3, and show them in different colors in Fig. 19. Compact and semicompact galaxies have slightly lower masses, with typical values of $8.06\begin{array}{c}+0.65\\ -0.54\end{array}$ and $8.15\begin{array}{c}+0.58\\ -0.59\end{array}$ respectively, compared to $8.38\begin{array}{c}+0.56\\ -0.67\end{array}$ for extended galaxies. Driven mostly by this difference and the EW – mass anticorrelation, compact galaxies show lower typical EW values. Nevertheless, for the lower mass range, compact galaxies show clearly higher EW values than extended ones even at similar masses. This suggests that an additional physical effect favors higher EW values in compact galaxies over extended ones.

The linear relation shown in Khostovan et al. (2015) and Reddy et al. (2018) between EW([OIII]) and stellar mass appears to have a limit at EW([OIII]) ∼3000 Å (Reddy et al. 2018). This limit, reached for galaxies dominated by the star-forming burst, is likely to be physical in nature rather than purely observational because for galaxies with a given mass, there is no bias against selecting higher EW systems. Understanding the reason for this limit might provide insight into the star formation cycle in the most extreme environment, in galaxies at the highest redshifts and with the lowest masses and metallicities. Our sample provides several previously unknown nearby candidates around that limit, which can be targeted with follow-up spectroscopic observations to extract relevant conclusions for this topic.

4.2.2. SFR and main sequence

One of the most important properties of star-forming galaxies is their SFR. We used one of the most prominent methods to estimate this parameter: the flux in the Hα emission line. In particular, we used the flux estimated by the CIGALE SED fit, with the prescription in Kennicutt (1998) to transform Hα flux into SFR. We correct for extinction using the E(B − V) value from the SED fit and the Calzetti et al. (2000) extinction law with R_V = 4.05 for consistency with the CIGALE implementation (using R_V = 3.1 would only create an offset with a median value of ∼ − 0.04 dex in SFR). We used the Hα flux because it is the most direct SFR tracer that we can access with the J-PLUS data, and using other sources would entail other corrections that are beyond the scope of this paper. The Hα flux value, even if calculated using CIGALE (and thus not entirely independent of the stellar mass value we used), agrees very well with both photometric and spectroscopic estimates (see Figs. 9 and 14). Taking the SFR estimate by CIGALE directly would be much less independent and make it harder to compare with other works in the literature. In any case, we only used the SFR values to illustrate the extreme nature of our sample, not to provide any definitive physical conclusions.

A tight correlation has been proven to exist between the stellar mass of galaxies and their SFR; this is often called “star formation main sequence” (e.g., Brinchmann et al. 2004; Elbaz et al. 2007; Noeske et al. 2007; Speagle et al. 2014). This relation has been used to define both starburst galaxies (above the relation) and quiescent galaxies (below it). Our selection method is intentionally oriented to identify galaxies with very high SFR, and thus our data cannot be used to compute the main-sequence parameters. In Fig. 20 we plot the SFR of the galaxies in our sample as a function of their stellar mass, as well as the constant specific SFR (sSFR) lines. We overplot several main-sequence relations for different redshifts (Leslie et al. 2020; Elbaz et al. 2007), with the fits derived in Duarte Puertas et al. (2017) and Vilella-Rojo et al. (2021) as local comparisons. It is clear that our sample of EELGs lies well above the local main-sequence relation, with sSFR values more comparable to those of typical high-redshift star-forming galaxies. We computed the difference between the SFR values measured in our sample and those expected if the galaxies followed the Vilella-Rojo et al. (2021) relation exactly. The typical difference is very high, with our EELG sample showing SFR values $1.31\begin{array}{c}+0.55\\ -0.46\end{array}$ dex higher. In particular, the differences are higher at lower stellar masses, with M_⋆ < 10^7.5 M_⊙ galaxies having $1.71\begin{array}{c}+0.54\\ -0.48\end{array}$ dex difference and those with M_⋆ > 10^8.5 M_⊙ presenting $1.08\begin{array}{c}+0.34\\ -0.32\end{array}$ offset. This is partially due to selection effects because the detection of galaxies with low stellar masses and low sSFR is less likely. Nevertheless, the detection of high sSFR low-mass galaxies and the lack of high sSFR high-mass galaxies is not due to biases. This result is consistent with many works that have discovered that low-mass galaxies undergo bursts of star formation that are more intense and more frequent than their higher-mass counterparts, and therefore they are more likely to be observed in a starburst phase (Sparre et al. 2017; Guo et al. 2016; Khostovan et al. 2021). This effect is already present (albeit softer) in typical star-forming galaxies, hence the slope of the main sequence is lower than one.

Fig. 20. SFR as a function of stellar mass for the galaxies in the EELG sample (red dots) and the rest of the candidate sample (gray dots). Large black dots represent the median values for the EELG sample grouped into mass bins. Colored lines represent main-sequence fits from the literature: Vilella-Rojo et al. (2021) in blue, Duarte Puertas et al. (2017) in red, Elbaz et al. (2007) in orange, and Leslie et al. (2020) in black (z ∼ 1) and green (z ∼ 5). Dashed and continuous lines follow the prescription in Fig. 19. Dotted black lines represent constant sSFR.

These differences place the vast majority of our sample above the threshold for galaxies undergoing starbursts, which is usually set around 0.48–0.6 dex above the main sequence (Rodighiero et al. 2011; Elbaz et al. 2018; Bluck et al. 2020). They are in fact more similar to main-sequence galaxies at very high redshift, with some reaching values that are typical at z ∼ 5, considering the extrapolation to low-mass galaxies of the main sequence at that redshift.

4.2.3. Comparison with typical low-redshift star-forming galaxies

Some properties derived from the SED fitting process can be complex to evaluate based on the uncertainties in the models, the assumptions taken (most notably the SFH), and the relatively sparse wavelength coverage of our data. In order to provide a fair comparison with the literature, we performed the same SED fitting analysis, adopting the same assumptions on a different sample of galaxies. The absolute values obtained for ages, metallicities, and so on may differ from others due to our choices in the SED fitting process and not due to physical differences in the galaxies, but the relative difference in these magnitudes when the same method is applied provides more accurate insights.

In order to do this, we chose the sample from Vilella-Rojo et al. (2021) as a comparison. In that work, the authors used the J-PLUS first data release (DR1; Cenarro et al. 2019) to identify 805 local (z < 0.017) star-forming galaxies selected as Hα emitters. The comparison is therefore straightforward because the filters that were used are the same. CIGALE was run over the Vilella-Rojo et al. (2021) sample of Hα emitters, with the same parameters as for our analysis (Sect. 3.1.1). Even if their galaxies are similar to ours, they selected brighter (r < 18 mag) objects at a slightly lower redshift. Therefore, their typical masses are higher than ours, with median and 1σ limits of $log(M_{\star }/M_{\odot })=8.64\begin{array}{c}+1.12\\ -0.88\end{array}$ compared to $log(M_{\star }/M_{\odot })=8.15\begin{array}{c}+0.55\\ -0.59\end{array}$. In order to correct for this difference, we limited our analysis to their galaxies showing stellar masses below 10^8.9 M_⊙, which lowers their typical mass values to $log(M_{\star }/M_{\odot })=8.21\begin{array}{c}+0.44\\ -0.80\end{array}$, in agreement with our sample (see panel d in Fig. 21).

Fig. 21. Histograms comparing the properties of our EELG sample (blue) with a subsample of typical star-forming galaxies from Vilella-Rojo et al. (2021) (red). We show the best-fitting values of several magnitudes, obtained by fitting J-PLUS photometry to SED models using CIGALE. From left to right and top to bottom, they are the age of the young stellar population, the mass ratio of the young and the old populations, the ionization parameter, and the stellar masses.

The results of the comparison between the other SED parameters are shown in the remaining panels in Fig. 21 (we only plot the most relevant histograms for this analysis). The most striking difference appears in the age of the young stellar populations, with extremely low values for our extreme [OIII] emitters sample ($3.0\begin{array}{c}+2.7\\ -1.0\end{array}$ Myr) compared to more typical values for the Vilella-Rojo et al. (2021) sample (7$\begin{array}{c}+2\\ -1\end{array}$ Myr). This result was expected, considering the difference in the sample selection in their work (with a minimum EW in Hα of around 12 Å) and ours (with a minimum EW in [OIII] of 300 Å), resulting in 95% of our galaxies showing Hα EWs higher than 56 Å, and 90% above 135 Å. It is clear in all models of star formation that the EW of the emission lines decreases rapidly after the initial burst, in a few million years. Therefore, galaxies with higher EW values will tend to have younger bursts of star formation. Other parameters (fraction of young population, metallicity, ionization parameters) will also play a role, but the very high EW in our work necessarily implies very young ages.

In line with the EW differences stated, the mass fraction of the young population is slightly higher in our sample than in Vilella-Rojo et al. (2021), with $log(M_{\mathrm{young}}/M_{\mathrm{old}})=-2.0\begin{array}{c}+0.7\\ -0.3\end{array}$ and $log(M_{\mathrm{young}}/M_{\mathrm{old}})=-2.3\begin{array}{c}+1.0\\ -0.3\end{array}$, respectively. The metallicities (not shown in the figure) are indistinguishable between the samples, with a p-value of 0.27 in the Kolmogorov-Smirnoff test. The ionization parameter shows a difference that extends from $logU=-3.0\begin{array}{c}+1.0\\ -0.5\end{array}$ in their sample to $logU=-2.5\begin{array}{c}+1.0\\ -0.5\end{array}$ in ours. However, it must be kept in mind that mass fraction and ionization parameter are subject to degeneracies and are hard to accurately measure just with the photometric data we have.

While the results presented in this work are insufficient to explain why the EELGs are undergoing such extreme events, some inferences can be made. The very young ages of the bursts indicate that it is likely that these galaxies only show these properties for briefs periods of time, which explains in part why they are so uncommon. Many more galaxies may have undergone similar starburst phases in other moments of their SFH (see, e.g., Sánchez Almeida et al. 2008, 2018). The trigger of the extreme star formation events remains an open question, but some hypotheses have been presented in the literature. For decades, mergers have been considered likely to be causing a fraction of the events (Barnes & Hernquist 1991), and more recently, their importance has also been shown in dwarf galaxy starbursts both from simulations (Bekki 2008) and observations (Stierwalt et al. 2015). Disrupted morphologies have also been observed with increased frequency in EELGs (e.g., Calabrò et al. 2017). While in our work most galaxies are compact, there is a significant fraction that show extended, often complex and clumpy morphologies. In addition, while some may appear undisturbed, traces of a recent merger may be hidden due to their low surface brightness (Martínez-Delgado et al. 2012). Another key driver for starbursts is the infall of cold gas into the galaxy (Ceverino et al. 2010; Sánchez Almeida et al. 2015), which is in line with the typically low metallicity of these systems, especially if the gas is near pristine and falling from the cosmic web. This process also explains why the EELGs are more common at higher redshift, where cold gas was much more abundant. Another physical process that could trigger the starbursts is the star formation feedback (Sparre et al. 2017). It could be due to the removal of gas from the galaxy, that later rains down on it, causing the starburst, or by compressing the gas inside the galaxy, triggering further star formation (Tenorio-Tagle et al. 2005). While distinguishing the causes of the extreme events of star formation is beyond the scope of this paper, the sample presented here it can be useful in this regard in future analyses.

5. Summary and conclusions

We have used the J-PLUS DR2 to select a sample of 466 EELGs in the local Universe over 2176 square degrees, with EW[OIII] above 300 Å. Out of these, 411 (89%) were previously unclassified as EELGs.

To build the sample, we identified objects with excess emission in the unique ∼200Å J-PLUS J0515 filter compared to the r band. To remove contaminants (stars and high-redshift systems), we relied on color-color and color-magnitude diagrams using IR WISE observations, reaching 96% purity and 92% completeness considering the available SDSS spectra.

We performed SED fitting on the J-PLUS data using the CIGALE software. We found the following key properties of the sample of galaxies:

• Low stellar masses, with typical values of $log(M_{\star }/M_{\odot })=8.13\begin{array}{c}+0.61\\ -0.58\end{array}$.

• Very young ages, $3\begin{array}{c}+2.7\\ -2\end{array}$ Myr, the main driver of the very high EW.

• Moderately low metallicities, Z=$0.008\begin{array}{c}+0.012\\ -0.004\end{array}$, and high ionization parameters, $logU=-2.50\begin{array}{c}+1.00\\ -0.35\end{array}$, which also contribute to the exceptionally high EW.

• Low extinction, $E(B-V)=0.1\begin{array}{c}+0.2\\ -0.1\end{array}$, typical of low-mass galaxies with high star formation activity.

We have compared our sample with several results in the literature. The agreement for those galaxies that have spectroscopic data is very good, comparing both line fluxes and EW. We measured negligible offsets in the comparisons between spectroscopic and our photometric J-PLUS measurements in the Hα and [OIII]5007 lines, with very little scatter in fluxes (∼0.18 dex) and EWs (∼0.1 dex).

We also analyzed the efficiency of our selection process compared to other recent results. We find that the number density of EELGs we calculated agrees with magnitude-limited spectroscopic surveys such as the main SDSS legacy survey and GAMA, while reaching fainter galaxies than both and much wider areas than GAMA. We measured comparable densities as narrow-band surveys, with caveats due to selection criteria. We are, in contrast, much more efficient than searches made using broadband surveys, such as the blueberry galaxies (Yang et al. 2017) or the green peas (Cardamone et al. 2009). We find 20–50 times more EELGs per unit of volume when controlling for EW and magnitude or stellar mass. It is also likely that we are able to access lower metallicity systems than broadband surveys because we are not directly biased against high Hα/[OIII]5007 systems, which is a property of extremely metal-poor galaxies (Senchyna & Stark 2019). The [OIII]5007 EW and SFR as a function of mass diagrams place our sample of EELGs in the typical values for the high-redshift Universe (z ∼ 2 − 5) while being located at z < 0.06.

In conclusion, we have presented a sample of mostly previously unclassified EELGs, which despite residing in the local Universe, share characteristics (mass, sSFR, EW, and metallicity) that make them similar to those at high redshift. Therefore, spectroscopic follow-up of this sample may shed light onto the properties of the galaxies forming in the very early Universe.

Acknowledgments

We thank the anonymous referee for their helpful comments. This work is based on observations made with the JAST80 telescope at the Observatorio Astrofísico de Javalambre (OAJ), in Teruel, owned, managed, and operated by the Centro de Estudios de Física del Cosmos de Aragón. We acknowledge the OAJ Data Processing and Archiving Unit (UPAD, Cristóbal-Hornillos et al. 2012) for reducing and calibrating the OAJ data used in this work. Funding for the J-PLUS Project has been provided by the Governments of Spain and Aragón through the Fondo de Inversiones de Teruel; the Aragón Government through the Reseach Groups E96, E103, and E16_17R; the Spanish Ministry of Science, Innovation and Universities (MCIU/AEI/FEDER, UE) with grants PGC2018-097585-B-C21 and PGC2018-097585-B-C22; the Spanish Ministry of Economy and Competitiveness (MINECO) under AYA2015-66211-C2-1-P, AYA2015-66211-C2-2, AYA2012-30789, and ICTS-2009-14; and European FEDER funding (FCDD10-4E-867, FCDD13-4E-2685). The Brazilian agencies FINEP, FAPESP, and the National Observatory of Brazil have also contributed to this project. L.A.D.G. thanks for financial support from the State Agency for Research of the Spanish MCIU through the ‘Center of Excellence Severo Ochoa’ award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709). This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This project uses data from the Legacy Surveys. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID 2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID 2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID 2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant 11433005). This project uses data from the Sloan Digital Sky Survey survey. Funding for the SDSS IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This paper uses data from the GAMA survey. GAMA is a joint European-Australasian project based around a spectroscopic campaign using the Anglo-Australian Telescope. The GAMA input catalogue is based on data taken from the Sloan Digital Sky Survey and the UKIRT Infrared Deep Sky Survey. Complementary imaging of the GAMA regions is being obtained by a number of independent survey programmes including GALEX MIS, VST KiDS, VISTA VIKING, WISE, Herschel-ATLAS, GMRT and ASKAP providing UV to radio coverage. GAMA is funded by the STFC (UK), the ARC (Australia), the AAO, and the participating institutions. The GAMA website http://www.gama-survey.org/.

References

Appendix A: ADQL query

In this appendix we reproduce the ADQL query used in Section 2 to obtain the parent sample for this work. In addition to the main condition (significant emission in J0515 compared to r), described in Sect. 2, we include some other restrictions that we list below.

• Objects must have been detected with at least 3σ significance in g, r, and J0515.

• The flux in the J0515 and g filters cannot be compatible, considering the uncertainties. This is to ensure that there is significant emission in J0515, not just a very steep slope toward r.

• The r fluxes in both PSFCOR and AUTO photometry must be higher than a certain threshold. This is only to ensure appropriate detections and avoid numeric errors because the thresholds are much lower than the final imposed threshold (r < 20 mag).

• Objects with FLAGS other than 0, 1, or 2 are rejected. We only kept those with either no issues (0), close neighbors (1), deblended (2), or a combination.

• Objects with any MASK_FLAG are rejected.

SELECT flu.*, z.*, s.*, e.*
FROM jplus.FLambdaDualObj flu,
jplus.PhotoZLephare z,
jplus.StarGalClass s, jplus.MWExtinction e
WHERE flu.TILE_ID = z.TILE_ID
AND z.TILE_ID = s.TILE_ID
AND s.TILE_ID = e.TILE_ID
AND flu.NUMBER = z.NUMBER
AND z.NUMBER=s.NUMBER
AND s.NUMBER=e.NUMBER
AND flu.FLUX_AUTO[jplus::rSDSS]>44.8
AND flu.FLUX_PSFCOR[jplus::rSDSS]>17.8
AND flu.FLUX_RELERR_PSFCOR[jplus::rSDSS]
<0.333
AND flu.FLUX_RELERR_PSFCOR[jplus::gSDSS]
<0.333
AND flu.FLUX_RELERR_PSFCOR[jplus::J0515]
<0.333
AND (flu.FLUX_PSFCOR[jplus::J0515] -
flu.FLUX_PSFCOR[jplus::rSDSS])/
flu.FLUX_PSFCOR[jplus::rSDSS] > 1
AND (flu.FLUX_PSFCOR[jplus::J0515]
- flu.FLUX_RELERR_PSFCOR[jplus::J0515]
*flu.FLUX_PSFCOR[jplus::J0515])
>(flu.FLUX_PSFCOR[jplus::gSDSS]+
flu.FLUX_RELERR_PSFCOR[jplus::gSDSS]
*flu.FLUX_PSFCOR[jplus::gSDSS])
AND flu.FLAGS[jplus::rSDSS] <4
AND flu.MASK_FLAGS[jplus::rSDSS] = 0
AND flu.FLAGS[jplus::J0515] < 4
AND flu.MASK_FLAGS[jplus::J0515] = 0

Appendix B: Recalculation of J-PLUS photometry for extended sources with a blending flag

As mentioned in Section 2.2.6, some of the objects in our sample present photometric flags indicating that they have been deblended from other source. In some cases, this just means that there was a different object nearby (e.g., a star or a different galaxy), and the photometric data we used only belong to the galaxy in which we are interested. Nevertheless, in some cases, the nearby object is in fact the main galaxy body, and the object we detected is only a specific region of the galaxy. Keeping these regions as "galaxies" would bias our analysis toward lower masses and higher EWs. In order to avoid this, we ran SExtractor over the J-PLUS images.

B.1. Visual classification

We inspected the 371 objects with deblending flag in the main sample, selecting those where our object of interest is in fact part of a larger galaxy. While this can ultimately be a judgment call, we followed some rules to use our own SExtractor photometry of an object instead of the original J-PLUS one. These rules are listed below.

• If in the Legacy survey image (or in the PANSTARRS, if outside the Legacy footprint) there is a clear connection of blue emission between our object and an extended galaxy that is not covered by the J-PLUS aperture.

• If there is a large galaxy nearby, with several HII regions showing a similar color as our object, at similar distances.

Some cases were not considered as improper deblending:

• If the color of the object near our target is very red and no clear structure is found between them, we consider that they do not belong to the same physical object and keep the J-PLUS original photometry.

• If the object and our target have a similar size and color and there is no emission between them, we consider them as satellites and keep the original J-PLUS photometry.

In some cases, particularly involving mergers, the decision between considering a system one galaxy or several based upon the photometry is arbitrary, to a certain degree. We chose to keep the original photometry if the system was very disturbed and our object was physically separated from the larger galaxy or galaxies. We also preferred to keep the original photometry when the SExtractor apertures were contaminated significantly by other nearby objects that were masked out in the J-PLUS photometry, however. After this visual inspection, we identified 175 galaxies where the original J-PLUS photometry was appropriate, and 179 where additional SExtractor runs where necessary to capture all the flux from the galaxy. In addition, we identified 20 spurious objects that were removed from the sample. Most of them were spikes from bright stars, measured over large apertures. An example of a galaxy with deblending flag where we decided to keep the SExtractor photometry is shown in Fig. B.1.

Fig. B.1. Example of an improperly deblended source. Left panel: Legacy surveys RGB image of the galaxy, with the AUTO ellipse aperture used in the J-PLUS catalog overlaid on top. Right panel: J-PLUS r band image of the galaxy, with the contour of the aperture defined in the custom SExtractor run that covered the whole galaxy overlaid on top.

B.2. SExtractor runs

We ran SExtractor over the sample of 179 galaxies that were not accurately represented by the original J-PLUS photometric aperture, using the default SExtractor parameter set except for a few parameters shown in Table B.1 (Set 1). After the first run using Set 1, we visually inspected the result and considered that for 36 of these galaxies, neither the SExtractor run nor the original photometry were appropriate. In most of these cases, our target was a region of a larger galaxy, but both the original photometry and our first SExtractor run failed to identify them as the same object. We chose to run SExtractor with a second set of parameters (Set 2 in Table B.1). This was enough to recover accurate photometry for 27 out of the 36 galaxies. The remaining 9 galaxies were removed from the sample because they are so extended that they would not make the EW cut to join the EELG sample.

Appendix C: CIGALE SED fits and parameters

In this appendix we review the SED fitting process in more detail and describe in full the CIGALE input parameters we used. We reproduce here the full set of input parameter in the main SED fit presented in this work.

   [[sfh2exp]]

e-folding time of the main stellar
population model in Myr.
tau_main = 50

e-folding 1 of the late starburst
population model in Myr.
tau_burst = 1

Mass fraction of the late burst
population.
f_burst = 0.0005, 0.0025,
0.005, 0.01, 0.03, 0.05, 0.075

Age of the main stellar population in the
galaxy in Myr. The precision is 1 Myr.
age = 200, 500, 1000, 2000, 5000

Age of the late burst in Myr.
The precision is 1 Myr.
burst_age = 1, 2, 3, 4, 5, 7, 12

Value of SFR at t = 0 in M_sun/yr.
sfr_0 = 1.0

Normalise the SFH to produce one
solar mass.
normalise = True

   [[bc03]]

Initial mass function: 0 (Salpeter) or
1 (Chabrier).
imf = 0

Metalicity. Possible values are:
0.0001, 0.0004, 0.004, 0.008, 0.02, 0.05.
metallicity = 0.0001, 0.0004, 0.004,
0.008, 0.02

   [[nebular]]

Ionisation parameter
logU = -4.0, -3.5, -3.0, -2.5, -2.0, -1.5

Fraction of Lyman continuum photons
escaping the galaxy
f_esc = 0.0, 0.2

Fraction of Lyman continuum photons
        absorbed by dust
f_dust = 0.0

Line width in km/s
lines_width = 300.0

Include nebular emission.
emission = True


   [[dustatt_calzleit]]

E(B-V)*, the color excess of the stellar
continuum light for the young population.

E_BVs_young = 0, 0.1, 0.2, 0.3, 0.5

Reduction factor for the E(B-V)* of the old
population compared to the young one (<1).
E_BVs_old_factor = 0.44, 1

Central wavelength of the UV bump in nm.
uv_bump_wavelength = 217.5

Width (FWHM) of the UV bump in nm.
uv_bump_width = 35.0

Amplitude of the UV bump.
For the Milky Way: 3.
uv_bump_amplitude = 0.0

Slope delta of the power law modifying
the attenuation curve.
powerlaw_slope = 0.0

It is important to note that the main focus of the SED fits is to compute accurate values for the continuum and the emission line fluxes, while also providing a stellar mass estimate. The specific results for the stellar population parameters are subject to our choices of free parameters and several physical degeneracies. Extracting high-level scientific conclusions from them is therefore very challenging. We performed several CIGALE runs on fractions of the main sample, with variations on the input parameters, to test the effect that different choices may have on the results.

The SFH we chose, with two single stellar populations, represents two almost instantaneous bursts of star formation. Other possibilities may also provide good fits to the data, for instance, including options for a more extended SFH in the old stellar populations (τ_old = 50, 500, and 1000 Myr). This would provide accurate fits as well and have almost no impact on the main results (negligible in log(U) or E(B-V), and just 0.1 dex in stellar mass, 0.02-0.04 dex in emission line fluxes, and 0.6 Myr in age of the young population). Even comparing the age of the old stellar population, we only see a small scatter (0.15 dex) and no offset. Attempting to fit the sample of galaxies only with a young population and entirely removing the old one results in generally successful fits with slightly higher typical reduced χ² values that reach from ${0.19}_{-0.10}^{+0.17}$ to ${0.26}_{-0.11}^{+0.20}$. These relatively successful fits are likely due to the fact that the optical flux generated in these galaxies is dominated by the nebular and stellar emission from the young population. To compensate for the lack of an old population in the models, the young populations show ages slightly older (by ∼ 3 Myr) and higher extinction levels (with typical E(B-V) of ∼ 0.3, which is 0.2 mag higher than in the two population models). This increased extinction estimation is significantly higher than the typical values measured in the available SDSS spectra of this sample (∼ 0.15 mag) or those in the literature for our mass range (∼ 0.1 mag, Garn & Best 2010; Duarte Puertas et al. 2017). Because of the large difference in mass-to-light ratios, the stellar masses would be about 1 dex lower if there were only a young population of stars. In conclusion, we consider it very unlikely that these galaxies are truly forming their first generation of stars, especially given the results shown in many detailed analyses of similar populations of galaxies (see Sect 3.1.1).

In our main run, we only considered two values (0 and 0.2) for the fraction of Lyman continuum photons escaping the galaxy (f_esc), and only one value (0) for the fraction of Lyman photons absorbed by dust (f_dust). We included some variation to provide more flexible fits, but we did not attempt to extract physical conclusions from these parameters. An increase in the set of possible values to f_esc = 0, 0.2, 0.5 and f_dust=0, 0.5 yields very similar results to the original SED fits: the median values of the difference between parameters remain virtually the same, with only small 1σ scatters: 0.08 dex in stellar mass, 0.02-0.04 dex in emission line fluxes, 1.5 Myr in age of the young population, 0.05 mag in E(B-V), and negligible in log(U).

Finally, we changed our assumption of a Salpeter (1955) IMF to a Chabrier (2003) IMF. In this case, we see the expected 0.23 dex offset in the stellar mass, along with a 0.16 dex scatter. The offsets in the rest of the parameters are negligible, with low scatter, as in the previous tests: 0.04 - 0.07 dex in emission line fluxes, 1.5 Myr in the age of the young population, 0.1 mag in E(B-V), and negligible in log(U).

All Tables

All Figures

Fig. 1. Illustrative example of an extreme [OIII] emitter, the J-PLUS source 67834-5013 (RA = 239.1019, Dec = 48.1127, zspec = 0.050). Left panel: color composite of the galaxy, obtained from the gri J-PLUS images. The sky location of the source is shown as a black dot, and the white ellipse indicates the three effective radius contours for the source. Right panel: J-PLUS 12-band PSFCOR photometry of the galaxy. The squares show the five SDSS-like filters (ugriz), and circles indicate the seven medium-band filters (J0378, J0395, J0410, J0430, J0515, J0660, and J0861). The solid line shows the spectra from SDSS with a downgraded resolution of R ∼ 180, normalized to the flux in the filter J0660. The location of the most prominent emission lines is marked: [OIII] (traced by J0515 and g), Hα (traced by r), and [OII] (traced by J0395).

In the text

Fig. 2. Color-color diagram used to separate low redshift ELGs from other types of sources. The g and r data are taken from J-PLUS DR2, and the W1 and W2 from the unWISE catalog. All candidates from the sample with J0515 excess flux are shown in dark grey dots. Sources with SDSS spectra are shown in colors according to their nature: red squares for QSOs, blue circles for galaxies at z < 0.06, brown triangles for galaxies at z > = 0.06, and green diamonds for stars. Filled symbols (or outlines only) are used if the sources are brighter (fainter) than r = 20 mag. The objects selected as candidates to be EELGs are enclosed by the blue lines, while those to the left of the dashed red line are visually inspected to select only the extended ones as EELGs. Top left corner: the typical errors (for sources with r < 20 mag) are shown.

In the text

Fig. 3. Color-magnitude diagram used to separate low redshift ELGs from other types of sources, for those objects where the W2 flux was not available (grey dots). All sources with SDSS spectra are shown in colors. The color code for those is the same as the one used in Fig. 2, but in this case objects with available W2 measurements are shown in a lighter shade. The objects selected as candidates to be EELGs are enclosed by the blue lines, while those to the left of the dashed red line are visually inspected to select only the extended ones as EELGs. Top right corner: typical errors (for sources with r < 20 mag) are shown.

In the text

	Fig. 4. Diagram showing the effective radius as a function of r magnitude, used to select candidates to be EELGs for those objects with excess flux in J0515 and no data in the unWISE catalog (grey dots). All sources with SDSS spectra are shown in colors. The color code for those is the same as the one used in Fig. 2, but in this case objects with W1 data are shown in a lighter shade.
In the text

Fig. 5. Example of the J-PLUS SED and CIGALE fit of an EELG. Red triangles correspond to the observed photometry, grey dots to the CIGALE synthetic photometry, and blue squares to the CIGALE synthetic photometry only considering the stellar and nebular continua (not considering the emission lines). The thick lines represent the CIGALE synthetic spectrum, considering different components: the grey line takes into account the full model, while the blue does not take into account the emission lines. The red line represents the stellar continuum, while the brown and green only consider the old and young stellar populations, respectively. The thin grey lines in the background represent the transmission profiles of all 12 J-PLUS filters.

In the text

Fig. 6. Histograms showing the best-fitting values for the J-PLUS galaxies, derived using the SED fitting software CIGALE. We show in gray the results for the 1493 galaxies in the candidate sample and in red the 466 galaxies in the EELG sample (EW[OIII] > 300 Å). From left to right and top to bottom: we show the age of the young population, the stellar metallicity, the ionization parameter, the color excess, the mass ratio of young and old populations, and the total stellar mass. The gaps observed in some of the histograms are due to the sampling in the parameters (see Table 2 and Appendix C).

In the text

	Fig. 7. Histogram of EW([OIII]), measured using the J-PLUS J0515 photometry and the CIGALE fits for the continuum. The filled white histogram represents the whole sample of candidates, and the filled red histogram shows only the galaxies with S/N > 3 in the EW([OIII]) value. The blue vertical line represents the minimum EW threshold for selecting EELGs in this work, 300 Å.
In the text

	Fig. 8. Comparison between the emission line flux measured in the J0515 filter and the [OIII]5007+4959 flux estimated by the CIGALE SED fit for the EELG sample. The black line represents the one-to-one relation.
In the text

	Fig. 9. Comparison between the emission line flux measured in the r filter and the Hα flux estimated by the CIGALE SED fit for the EELG sample. The black line represents the one-to-one relation.
In the text

	Fig. 10. Diagram used to asses the completeness of the EELG sample. We plot the number of galaxies per bin of log(EW[OIII]) as a function of log(EW[OIII]). We overplot the linear fits to the data considering three different limits in the lower end, from 2.45 log(Å) (black line), 2.55 (red), and 2.65 (green). The points considered in at least one of the fits are shown in a darker shade of gray.
In the text

Fig. 11. RGB composites (made using images taken with z, r, and g filters) of four galaxies in the EELG sample. Bottom right corner of each panel, a letter indicates the type of morphology: compact (c), semicompact (s), or extended (e). Bottom left corner of each panel: 10 arcsecond line for scale. All images have been obtained using the Legacy Survey viewer with data from the BASS and MzLS surveys (Zou et al. 2019; Dey et al. 2019).

In the text

	Fig. 12. CIGALE photometric redshift as a function of SDSS-derived spectroscopic redshift for the galaxies in our sample of candidates with available spectra (excluding the three misclassified objects). The line represents the one-to-one relation. In the top left corner of the figure, we show the minimum step in photometric redshift considered with CIGALE (0.0025) as error bars.
In the text

	Fig. 13. [OIII]5007 flux estimated using CIGALE SED fitting of J-PLUS data as a function of [OIII]5007 flux measured in SDSS spectra. Same sample as in Fig. 12. The black line represents the one-to-one relation.
In the text

	Fig. 14. Hα flux estimated using CIGALE SED fitting of J-PLUS data as a function of Hα flux measured in SDSS spectra. Same sample as in Fig. 12. The black line represents the one-to-one relation.
In the text

	Fig. 15. J0515 EW directly measured on J-PLUS data as a function of J0515 EW estimated convolving SDSS spectra with the filter transmission. Same sample as in Fig. 12. The black line represents the one-to-one relation.
In the text

	Fig. 16. r EW directly measured on J-PLUS data as a function of r EW estimated convolving SDSS spectra with the filter transmission. Same sample as in Fig. 12. The black line represents the one-to-one relation.
In the text

Fig. 17. Broadband color-color diagram, created using SDSS data. We plot our EELG sample in blue dots, alongside the blueberry galaxies in Yang et al. (2017) (black triangles) and samples of extremely metal poor galaxies (XMP) from Kojima et al. (2020) (red squares) and the other literature works (green diamonds), as reported in (Kojima et al. 2020). The blue line represent the Yang et al. (2017) sample selection limits, while the green lines follow the evolutionary tracks of models in Kojima et al. (2020). The small black dots represent the density contours of typical SDSS galaxies (Kojima et al. 2020).

In the text

Fig. 18. Distribution of the ratios [OIII]5007/Hα (filled red rectangles and dots) and [OIII]5007/Hβ (filled white rectangles and gray dots). The thick black line represents the median value, and the box ranges from the first to the third quartile. The error bars represent the maximum and minimum value without outliers, which are plotted as dots (and defined as the values that lie beyond 1.5 times the interquartilic range). We plot in this figure results from our work and several samples from the literature (Yang et al. 2017; Cardamone et al. 2009; Amorín et al. 2015; Onodera et al. 2020; Hsyu et al. 2018; Guseva et al. 2017; Pérez-Montero et al. 2020; Matthee et al. 2021; Izotov et al. 2021), as well as the KISS survey (Wegner et al. 2003; Gronwall et al. 2004; Jangren et al. 2005; Salzer et al. 2005) and the Hα dots class (Kellar et al. 2012; Salzer et al. 2020). We selected the subsample with extremely low metallicity in Hsyu et al. (2018) (H18_XMP), as well as galaxies with EW(Hβ) > 100 Å in both Hsyu et al. (2018) XMPs (H18_XMP Hβ > 100 Å) and the KISS survey (KISS Hβ > 100 Å). Different samples are separated by vertical blue lines.

In the text

Fig. 19. [OIII] EW of the EELG (color dots) and candidate (gray dots) samples as a function of stellar mass. We plot the galaxies in the EELG sample according to their morphology in different colors: red corresponds to compact, green to semicompact, and blue to extended objects. In addition, we overplot linear relations from Khostovan et al. (2016) at different redshifts (black at z ∼ 0.84, blue at z ∼ 1.42, pink at z ∼ 2, and orange at z ∼ 3.24). The linear relations are plotted in as continuous lines in the mass ranges in which they were defined, and as dashed lines in their extensions to lower masses.

In the text

Fig. 20. SFR as a function of stellar mass for the galaxies in the EELG sample (red dots) and the rest of the candidate sample (gray dots). Large black dots represent the median values for the EELG sample grouped into mass bins. Colored lines represent main-sequence fits from the literature: Vilella-Rojo et al. (2021) in blue, Duarte Puertas et al. (2017) in red, Elbaz et al. (2007) in orange, and Leslie et al. (2020) in black (z ∼ 1) and green (z ∼ 5). Dashed and continuous lines follow the prescription in Fig. 19. Dotted black lines represent constant sSFR.

In the text

Fig. 21. Histograms comparing the properties of our EELG sample (blue) with a subsample of typical star-forming galaxies from Vilella-Rojo et al. (2021) (red). We show the best-fitting values of several magnitudes, obtained by fitting J-PLUS photometry to SED models using CIGALE. From left to right and top to bottom, they are the age of the young stellar population, the mass ratio of the young and the old populations, the ionization parameter, and the stellar masses.

In the text

	Fig. B.1. Example of an improperly deblended source. Left panel: Legacy surveys RGB image of the galaxy, with the AUTO ellipse aperture used in the J-PLUS catalog overlaid on top. Right panel: J-PLUS r band image of the galaxy, with the contour of the aperture defined in the custom SExtractor run that covered the whole galaxy overlaid on top.
In the text