A MUSE/VLT spatially resolved study of the emission structure of Green Pea galaxies

Green Pea galaxies are remarkable for their intense star formation and serve as a window into the early universe. In our study, we used integral field spectroscopy to examine 24 of these galaxies in the optical spectrum. We focused on the interaction between their ionized interstellar medium and the star formation processes within them. Our research generated spatial maps of emission lines and other properties like ionization structures and chemical conditions. These maps showed that areas with higher levels of excitation are usually located where starbursts are occurring. Continuum maps displayed more intricate structures than emission line maps and hinted at low brightness ionized gas in the galaxies' outer regions. We also analyzed integrated spectra from selected areas within these galaxies to derive physical properties like electron densities and temperatures. In some galaxies, we were able to determine metallicity levels. Our observations revealed the presence of high-ionizing lines in three galaxies, two of which had extremely high rates of star formation. Our findings provide valuable insights into the properties and star-forming processes in Green Pea galaxies, contributing to our broader understanding of galactic evolution in the early universe.


Introduction
The cosmic dawn (6 ≲ z ≲ 10) marks a major phase transition of the Universe, during which the "first light" [metal-free stars (or the so-called PopIII-stars) and the subsequent formation of numerous low-mass, extremely metal-poor galaxies] appeared, putting an end to the dark ages.The details of the reionization history reflects the nature of these first sources, which is just starting to be constrained with the James Webb Space Telescope (JWST) (e.g., Curtis-Lake et al. 2022;Robertson et al. 2022;Naidu et al. 2022).However, exactly when and how the Universe was reionized remains one of the most important questions of modern astrophysics, and that will be explored during the next decades.
An immediately accessible approach to better understand the first sources is to identify galaxies at lower redshifts with properties similar to galaxies in the very early Universe (e.g., Schaerer et al. 2022;Chen et al. 2023).Among these local analogs, we find the objects known as GPs, which are a subset of galaxies inside the extreme emission line galaxies (EELGs) set (e.g., Pérez-Montero et al. 2021;Breda et al. 2022;Iglesias-Páramo et al. 2022).GPs are compact starbursts commonly found at redshift z ∈ (0.112, 0.360), corresponding to 2.5 − 4.3 Gyr ago.The upper size limit of these galaxies is 5 kpc in Hubble Space Telescope (HST) images (Yang et al. 2017a).The JWST has already shown the similarity between primeval galaxies (z ∼ 8) and GPs (Rhoads et al. 2022).The three high-z galaxies presented in Rhoads et al. (2022) are all strong line emitters, with spectra reminiscent of nearby GPs.They also show compact morphologies typical in GPs.So, without any doubt, GPs are excellent local analogs of high-redshift galaxies.
On average, a GP galaxy has a stellar mass of M ⋆ = 1 × 10 9 M ⊙ and a star formation rate (SFR) of 10 M ⊙ /yr (Cardamone et al. 2009;Izotov et al. 2011), and thus its mass-doubling timescale easily goes below 100 Myr.In this way, GPs present a very strong starburst similar to in high-redshift galaxies (Lofthouse et al. 2017).Furthermore, different studies suggest hot massive stars as the main excitation source in GPs, resulting in a highly ionized ISM (e.g., Jaskot & Oey 2013).

MUSE data
We studied a sample of GPs observed with MUSE (Bacon et al. 2010) at the Very Large Telescope (VLT; ESO Paranal Observatory, Chile).MUSE is a panoramic IFS, which, operating in its wide field mode (WFM), provides a field of view (FoV) of 1 ′ × 1 ′ with a spatial sampling of 0.2 ′′ and a FWHM spatial resolution of 0.3 ′′ − 0.4 ′′ .The data were obtained in nominal mode (wavelength range λ4750Å−λ9350Å) with a spectral sampling of about 1.07 Åpix −1 and an average resolving power of R ∼ 3000.The selection criteria corresponding to the observations consist in all the GPs presented in Cardamone et al. (2009) at a declination < +20, so there is visibility from the Paranal ob-servatory.The program ID corresponding to the observations of these galaxies is 0102.B − 0480(A) (PI: Hayes, Matthew).
We retrieved the fully reduced data cubes from the ESO archive.The data reduction was performed with MUSE Instrument Pipeline v. 1.6.1 with default parameters (Weilbacher et al. 2020), which consists of the standard procedures of bias subtraction, flat fielding, sky subtraction, wavelength calibration, and flux calibration.
The sample consists of 24 GPs.It is an unbiased representative set of all GP galaxies: stellar masses, redshifts, metallicities, SFR, [Oiii]5007 EWs, and line ratios of our galaxies spanning the ranges typical of GPs (e.g., Cardamone et al. 2009;Amorín et al. 2010).The names, positions, redshifts, and information about the observations of the galaxies are in Table 1.A histogram of the redshift of the GPs in our sample can be seen in Fig. 1.From a first inspection of the data cubes, our GPs appear compact and almost point-like except for a few of them.In order to ensure that we can spatially resolve these galaxies, the FWHM of all the sources within the FoV has been measured in white light 1 .We checked the objects in the FoV that are stars according to SDSS DR16 (Ahumada et al. 2020).Then, we defined the stellar-like FWHM (FWHM ⋆ ) as the median stellar FWHM, and, in case there are no SDSS stars in the FoV, we defined the FWHM ⋆ as the FWHM of the least extended object within the FoV.We considered a GP to be extended if its FWHM (FWHM GP ) follows this condition: After applying this criterion, 12 GPs are considered to be extended out of 24.The FWHM ⋆ and the FWHM GP along with the number of stars in the FoV can be seen in Table 2.
The previous analysis allows us to determine which galaxies present a resolved core.Nevertheless, thanks to the very high VLT/MUSE sensitivity, we can also check the extended nature of these galaxies in low surface brightness.In order to prove it, an analogous analysis of the previous one is realized but using instead the full width at 1 10 of the maximum (FW 1 10 M) and the full width at 1 100 of the maximum (FW 1 100 M).These parameters trace the extension of the objects at lower surface brightness.We consider that a GP is extended in the low surface brightness regions if one of the following conditions is satisfied: After applying this criterion, 11 GPs are considered to be resolved in the low surface brightness regions out of 24.The FW 1 10 M and FW 1 100 M of the stellar-like sources and the GPs can be seen in Table 2.
Due to the redshift distribution in our sample of GPs, the spectral coverage of MUSE provides us with most of the emission lines in the optical wavelengths.In particular, for all GPs we have all the emission lines from Hδ to [S ii]λ6731Å.For 13 of them we can reach the [Oii]λ3727Å line.
Rayleigh scattering of the atmosphere would spread the light from a source depending on wavelength, and the extension of the emission line maps at substantially different wavelengths (e.g., Hα vs Hβ and [Oiii]5007Å vs. [Oii]3727Å, 3729Å)) might be affected.Aiming to correct by this effect, a detailed analysis was performed sampling the FWHM for all the sources (within the FoV of each GP) along the entire wavelength range, showing that blue images in the cubes have a lower spatial resolution than red images.We derived a decrease in the FWHM (as wavelength increases) of the sources ranging from 2.25 × 10 −5 arcsec Å −1 to 1.15 × 10 −4 arcsec Å −1 .This information allowed us to apply the appropriate gaussian kernel to red images to effectively compare them with the blue images.
The MUSE pipeline corrects the data for the effect of atmospheric differential refraction.As a sanity check, we measured the center of all the objects detected in the cubes in the spectral ranges of [4900,5100]Å and [8950,9150]Å, finding a mean difference of only 0.07 ′′ (Weilbacher et al. 2020).

SDSS spectra
We retrieved SDSS-DR16 integrated spectra for all the GPs (Ahumada et al. 2020).The diameter of the SDSS fiber is 3".The wavelength coverage is 3800-9200Å.The spectral resolution is 1500 at 3800Å and 2500 at 9000Å.The pixel spacing log wavelength is 10 −4 dex and the exposure times are in the range of 2800 s − 8500 s.

Flux measurements and spatially resolved structure
In this section, we describe the methodology followed to measure the flux of the emission lines in the galaxies.Additionally, we present the spatially resolved maps generated from the spaxel-to-spaxel measurements of fluxes and continuum along with the so-called BPT diagrams (Baldwin et al. 1981;Kewley et al. 2006).

Emission line and continuum maps
For each spaxel, all the emission lines were fit to a gaussian profile, which gives us the flux of each line.The errors in the fluxes were calculated using the bootstrap method (Efron & Tibshirani 1985).We combined the line fluxes with the position of the fibers in the sky to create the maps of emission lines presented in this paper.
As an example, the emission line maps of the galaxy GP06 are shown in Fig. 2.Among the galaxy sample, GP06 is the one that presents the most complex structure in these maps, where even in dim lines like [S ii]6716Å and [S ii]6731Å a small bump is seen in the south-western position with respect to the central  burst.Furthermore, this galaxy presents the brightest nebular HeII4686Å emission (also detected in GP20 and GP15) in our sample.This line indicates the presence of very high energy photons (E > 4 Ry) that strongly ionize the gas and raise the electron temperature (e.g., Kehrig et al. 2015).The emission line maps corresponding to the rest of the galaxies are presented in Appendix A. All emission lines in our set of GPs peak in the center of the galaxy.This indicates that the total brightness of the galaxies is dominated by a compact region where the ionization sources are present.This region is slightly resolved for only a few GPs (e.g., GP13; see Fig.To make the continuum maps, we selected a rest frame spectral range common to all galaxies to integrate the continuum. Since each galaxy has a different redshift, they present a different rest frame spectral range.So, we can take the intersection of all these to define our rest frame spectral range for defining the continuum.This is indeed between 4154Å and 7026Å, very close to the spectral range of the human eye.To make these maps we removed the lines in the spectrum of each spaxel.Then, we integrated the spectrum in our selected rest-frame spectral range.In Fig. 3, we present the continuum map corresponding to GP06.The continuum maps of all GPs are in Fig. C.1.As seen in these maps, all galaxies present a richer structure than in the emission line maps.One of the reasons behind this is the fact that we are collecting more light, since the rest frame spectral range is 2872Å wide.In order to collect the same amount of light from a line, its EW must be similar, while their values are up to 2000Å.Continuum maps show bumps attached to the central part of the galaxy in many cases (e.g., GP03, GO06, GP13, GP15).Such disturbed morphologies of the continuum could indicate a spread in the underlying stellar population due to recent mergers (Lofthouse et al. 2017).
Even if we are collecting more light in the continuum than in the emission lines, the galaxies show a more compact appearance in the continuum maps than the [Oiii]5007 emission line maps.This possibly indicates that we are tracing the ionized gas in the outer parts of the galaxies, since the underlying stellar population is less extended than the ionized gas.In particular, GP06 presents two regions in the north east and north west (with respect to the Hα peak), where the [Oiii]5007 emission line extends much further (up to 10kpc) than the continuum (Figs. 2 & 3).In these galactic regions, there is almost no extinction (Hα/ Hβ ∼ 2.86), and the gas has higher excitation (log([Oiii]λ5007/Hβ) ∼ 0.7) in comparison to the other outer regions where log([Oiii]λ5007/Hβ) ∼ 0.45 (see Section 3, Figs. 4 and 6).Similar results were found for the galaxy SBS 0335-052E where [Oiii]5007 and Hα filaments with low Hα/ Hβ ratio were discovered (Herenz et al. 2023).These observations  point to the fact that these regions can be channels where LyC photons can escape.In particular, for the case of GP06 it is found that this galaxy presents an optically thick, neutral outflow along the line of sight (LOS) (Jaskot & Oey 2014), making the escape of high-energy photons difficult in this particular direction (with f esc (Lyα) = 0.05 as found in Jaskot et al. (2017)).However, this galaxy displays one of the highest ratios of [Oiii] to [Oii] in this study ([Oiii] / [Oii] = 6.56).This suggests that there could be potential escape paths for ionizing photons in other directions (Izotov et al. 2022).

Emission line ratio maps
Here, we present maps of some of the relevant line ratios to study the ionization structure of the gas in our set of GPs.These line ratios are corrected for reddening using the corresponding c(Hβ) for each spaxel; c(Hβ) was computed from the ratio of the measured-to-theoretical Hα/Hβ assuming the reddening law of Cardelli et al. (1989), and case B recombination with the electron temperature T e = 10 4 K and electron density n e = 100 cm −3, which give an intrinsic value of Hα/Hβ= 2.86.In Fig. 4, we show the uncorrected Hα/Hβ map for GP06.The Hα/Hβ map is a tracer of dusty regions, and the zones with higher extinction are likely to present dust that acts like a wall for high-energy photons (Weingartner et al. 2006).Hence, if there is an escape of LyC photons, the preferred direction would be the one that presents lower extinction.The emission line ratio maps for all GPs are provided in Appendix B.
The MUSE spectral coverage allows us to observe the [Oii]λλ3727, 3729 lines for 13 of our GPs due to their redshifts, so we can analyze the map of the [Oiii]λ5007 / [Oii]λλ3727, 3729 ([Oiii] / [Oii]) line ratio (see Appendix B), which traces the ionization of the gas, in these galaxies.As an example, in Fig. 5  High values of [Oiii]/Hβ correspond to the areas of ionized gas with relatively higher excitation.The [Oiii]/Hβ maps do not show significant spatial variations (with a maximum difference of < 0.1 dex) for any of the GPs except for GP06, indicating low spatial gradients in gas excitation.For this galaxy, the spatial variation of [Oiii]/Hβ goes up to 0.4 dex, with a maximum  [S ii]/Hα maps trace the opacity of the column of gas (i.e., gas in the line of sight in each spaxel) (Pellegrini et al. 2012).[Oi]/Hα peaks in the ionization front (edge of ionizationbounded regions).Low values on both coefficients indicate thin columns of gas (and as well, high gas excitation) where LyC photons are likely to escape (Paswan et al. 2022).We can see that most GPs present a blister-type morphology, which is characterized by an optically thin galactic center, and as we go to the outer parts (around 5-10kpc) the medium tends to become optically thick.(see Appendix B).An extremely low [S ii]/Hα value could indicate a hole from which photons can escape (Wang et al. 2021), being GP06, GP15, and GP20 the most extreme cases again (see Figs. B.6,B.15,and B.20).
The last two maps presented here are [S ii]λ6716/[S ii]λ6731 and [Oiii]λ4363/[Oiii]λ5007. Examples of maps corresponding to these line ratios are displayed in Fig. 7.These emission line ratios give us information about the electron density and electron temperature, respectively (e.g., Pérez-Montero 2017).Higher [S ii]λ6716/[S ii]λ6731 values correspond to lower electron densities, whereas higher [Oiii]λ4363/[Oiii]λ5007 values correspond to higher electron temperatures.The maps corresponding to these two line ratios do not present clear radial variances, while they show a wide variety of morphologies despite the low spatial extension of these dim lines (see Appendix B).In the particular case of GP06, the values of the [S ii]λ6716/[S ii]λ6731 ratio are higher where the gas shows lower excitation (as traced by [Oiii]/Hβ and [S ii]/Hα) and higher extinction regions (traced by Hα/Hβ).This indicates that for this galaxy the electron density tends to a decrease in the regions with higher dust content and lower gas excitation.Moreover, the highest values for [Oiii]λ4363/[Oiii]λ5007 (line-ratio indicator of the electron temperature) are found in the galaxies with HeII4686Å emission (see Figs. B.6,B.15 and B.20), which reinforces the existence of a harder ionizing radiation field in these objects (see, e.g., Kehrig et al. 2016).
We present the radial profiles of the emission line ratios previously mentioned.By plotting the emission line ratio as a function of radial distance, we effectively visualize the radial variations in the line ratios across the galaxy.This representation facilitates the identification of radial trends and provides a comprehensive understanding of how these line ratios evolve as we move away from the galactic center.
To obtain the radial profiles, we utilized a technique that integrates circular crowns centered around the peak of Hα emission.The code calculates the average value of the emission line ratio within each crown and plots it against the radial distance from the center.The radial profiles of Hα/Hβ are shown in Fig. 8.The rest of the radial profiles are in Appendix D.
One notable observation is that the overall radial tendencies of the emission line ratio profiles are similar across the entire set of galaxies.However, there are substantial variations in the absolute values of the ratios between different galaxies, typically on the order of 0.5 dex.In contrast, within each individual galaxy, the radial change in the emission line ratios is relatively small, typically less than 0.1 dex.This small radial variation within each galaxy could potentially be attributed to the low spatial resolution of our observations.Despite the low radial change within each galaxy, the mean radial profiles of all galaxies present clear tendencies for the various emission line ratios studied.These profiles provide valuable insights into the ionization state within the galaxies.
We observe a consistent decrease in the [Oiii] / [Oii] and [Oiii]5007/Hβ ratios and an increasing trend in the [S ii]/Hα, [Oi]/Hα, and [Nii]/Hα ratios as we move away from the galac-Fig.6: BPT maps corresponding to GP06.tic center.This implies that the highest level of ionization and density-bounded tracers are predominantly concentrated in the central regions of the galaxies.Still, the mean variation in the [Oiii]5007/Hβ and [Nii]/Hα ratios are really low (0.03 dex and 0.015 dex, respectively), indicating that the radial changes in these emission line ratios are nearly imperceptible.
Furthermore, the radial tendency in the Hα/Hβ ratio indicates a decreasing trend as we move toward the outer parts of the galaxies.This suggests a relatively higher level of extinction or enhanced dust attenuation toward the central regions, resulting in a higher Hα/Hβ ratio compared to the outskirts.
In contrast, the radial changes in the [S ii]λ6716/[S ii]λ6731 and [Oiii]λ4363/[Oiii]λ5007 ratios are small (on the order of the previously mentioned [Nii]/Hα ratio) and do not exhibit a clear radial trend.These ratios may be less sensitive to radial variations or they are influenced by other factors such as excitation conditions or observational uncertainties.
Overall, the mean radial profiles of the emission line ratios consistently reveal the spatial variations in the ionization state within the galaxies.The observed trends support the notion that the central regions of the galaxies exhibit higher ionization levels and greater density-bounded tracers, while the outer parts experience lower ionization conditions.

BPT diagrams
The BPT diagrams for all our GPs are shown in Fig. 9 on a spaxel-by-spaxel basis; there, we can see the line ratios [Oiii]5007/Hβ versus [Nii]/Hα, [S ii]/Hα, and [Oi]/Hαplotted.Small dots correspond to measurements from individual spaxels for each galaxy.Big green dots show the line ratios derived from the integrated spectra (see Section 4), these values are shown in Table F.2.
Overall, our GPs fall in the general locus of SF objects according to the spectral classification scheme proposed by Baldwin et al. (1981) and Kewley et al. (2006).This suggests that photoionization from hot massive stars is the dominant excitation mechanism within these galaxies.In particular, GPs are lo-Fig.7: Maps corresponding to GP06.Fig. 8: Radial profile of Hα/Hβ.The left image shows the profiles of all GPs with available information about the lines in gray.In black it is represented the mean of all galaxies and extends up to a radius where only 30% of the galaxies present values at larger radii.On the right is a zoomed-in image of the mean profile.The same display applies to all line ratios presented in the appendix.cated in the top left part of the diagram where the most extreme galaxies reside (i.e., lower metallicity and higher excitation of the ionized gas).Previous studies confirm this trend for GPs (e.g., Cardamone et al. 2009 Furthermore, there is no presence of the Fe[X]λ6374Å line in any GP.This line is a tracer of black hole (BH) activity.Luminosities of this line on the order of 10 36 − 10 39 erg/s correspond to the presence of BHs with a mass of ∼ 10 5 M ⊙ (Molina et al. 2021).Nevertheless, due to the distance of GPs, it is not possible to measure lines with luminosities << 10 40 erg/s.Such limitations lead us to conclude that GPs do not present actively accreting BHs with a mass >> 10 5 M ⊙ .
In addition, we used a method for estimating BH masses involving the use of a BH mass -stellar mass relation-which states that it is nearly independent of redshift (Shankar et al. 2020).Nevertheless, our range of stellar masses, which covers from 10 8.3 M ⊙ to 10 10 M ⊙ in the low end, is well below the predictive capacity of the relation.Thus, we are only using it for the GPs with higher masses (∼ 10 10 M ⊙ ); for these galaxies, the expected BH mass is no greater than 10 6.7 M ⊙ .
This estimation is clearly above 10 5 M ⊙ solar masses, which suggests that if these galaxies do contain BHs, these BHs have a low mass compared to the total stellar mass of the galaxy.It is important to note, however, that with the use of the Shankar et al. (2020) relation we cannot accurately predict BH masses in galaxies with M ⊙ < 10 9.8 .That is why it is so difficult to retrieve any information from the BHs in the low stellar mass GPs.Further observations and analysis are needed to confirm the presence and measure the exact masses of any BHs in these galaxies.
In the BPT diagrams, the distance from the center of the galaxy as a parameter is also represented.We can see a common tendency for all galaxies.As we go closer to the center of the galaxy, it shows a higher [Oiii]5007/Hβ ratio (i.e., higher Fig. 9: BPT diagram.Small dots in the color scale (from yellow to dark purple) correspond to spaxels in each galaxy.Big green dots correspond to the integrated spectra of the galaxy.Next to the green dots we can see the name of each galaxy.The color of the small dots corresponds to the distance to the center of the galaxy.The closer to the center the more yellow they are, and farther away from the center they become darker.The lines that delimitate each region are taken from Kewley et al. (2001) and Kauffmann et al. (2003).
excitation close to the center) and lower [S ii]/Hα and [Oi]/Hα ratios (i.e., the center of the galaxies are optically thinner, and neutral oxygen is more abundant in the outer part).Regarding the [Nii]/Hα ratio, it is almost independent of distance, possibly indicating that metallicity gradients are low.All these results are in agreement with the ones presented in the previous section (Section 3.2).

Properties of GPs from integrated spectra
We also took advantage of our IFS data to produce the 1D spectra of selected galaxy regions.To do so, for each GP, we added the flux in all spaxels for which the Hα flux measurements present a signal-to-sigma sky ratio greater than three (S /σ sky > 3).As an example, the integrated spectrum of the galaxy GP06 can be seen in Fig. 10.The integrated spectra of all GPs are shown in figure E.1.We did not find in any of the spectra features of WR stars.
We used the integrated spectrum of each GP to identify and measure the most relevant emission lines.The same procedure as the one used in Section 3 was used to calculate emission line fluxes, errors, and extinction correction.Tables in F lists all the data of the emission lines.To check the wellness in the flux measurement, the MUSE integrated spectra and the SDSS spectra were compared, retrieving that, on average, the ratio between the fluxes per line is 1.003 ± 0.038.Additionally, the spectral coverage of SDSS extends to shorter wavelengths than MUSE, enabling us to obtain the flux of the [Oii] line for the 11 GPs that do not present this line in MUSE.We used the MUSE fluxes in all other cases.

SFR
By means of the Kennicutt relation (Kennicutt Jr 1998), we derived the SFR through the Hα luminosity (Hα luminosity was retrieved using luminosity distances from the UCLA cosmology calculator 2 ) in the following way: S FR(M ⊙ /yr) = 1.26 × 10 41 L Hα (erg/s).In Fig. 11, we display sSFR versus stellar mass.Here, GPs occupy a different region in the diagram than present-day galaxies (galaxies at z<0.05) (Catalán-Torrecilla et al. 2015).GPs have lower masses and much higher sSFRs.In fact, GPs share the same space in this diagram with high-redshift galaxies z = 1.1 − 4. The mass-doubling timescale of GPs is on average 2 dex lower than that of present-day galaxies, reaching down to ∼ 15 Myr for GP20, GP22, GP18, and GP15.
The depletion timescale, which indicates the starburst duration, is on the order of the mass-doubling timescale (which is the inverse of the sSFR) only if the mass of available hydrogen for creating new stars (M HII ) is on the order of the M ⋆ .Nevertheless, if M HII > M ⋆ the starburst duration is larger than the mass-doubling timescale.If we consider that M HII ≃ M ⋆ , this could indicate that GPs are short-lived events, since they will not be able to support such an incredibly high SFR for a long time.The mechanisms that would stop the SFR are mainly the exhaustion of the gas that fuels star formation and the stellar feedback via supernovae (Amorín et al. 2012).Observationally, we have very little information on whether GPs will immediately quench or not.However, if they continue forming stars, they would quite rapidly build stellar mass and increase the stellar luminosity in the blue, and the EW of [Oiii]5007Å will decrease.Thus, they would no longer be selected as GPs, given the effective criterion of having high EW in [Oiii]5007Å.Using this kind of argument one could state that they are the analogs of the early phases of galaxies that reionized the Universe.

Electron density, temperature, and abundances of the ionized gas
In this section, we discuss how we derived chemical abundances of the ionized gas and electron temperatures (T e ) and densities (n e ).The emission line [Oiii]4363Å is essential to calculate electron temperature and hence abundances.There are six GPs (out of the 24 GPs analyzed in this work) in which this line can be detected (i.e., the line is above the three-sigma detection limit).For this subset of galaxies, we derived the T e values of In Fig. 12, we can see the electron density, electron temperature, and metallicity for the six GPs with the [Oiii]4363Å line measured.The same properties are also represented for 35 galaxies selected from the NASA-Sloan Atlas 3 with W(λ5007) > 1000Å and 37 galaxies from the COS Legacy Archive Spectroscopic SurveY (Berg et al. 2022) galaxies, which are the closest local analogs (z<0.18) of high-redshift galaxies in the epoch of reionization.The derivation of gas parameters in these galaxies was done by Peng (2021).The metallicity of our set of GPs is low and in agreement with previous results (i.e., 12 + log(O/H) = 7.6 − 8.4) (e.g., Amorín et al. 2010).For this subset of GPs with [Oiii]4363Å emission, the electron temperature ranges from 11500 K to 15500 K and the electron density ranges from 30 cm −3 to 400 cm −3 .Fig. 12: Electron density, electron temperature, and metallicity representation.Big labeled dots correspond to the six GPs where the direct method can be well applied.The rest of the small dots correspond to all galaxies selected from the NASA-Sloan Atlas with W(λ5007) > 1000Å and 37 COS Legacy Archive Spectroscopic SurveY galaxies.For a given stellar mass, GPs show lower oxygen abundances than the bulk of the SDSS galaxies.This is in accordance with the result from Amorín et al. (2010), where GPs follow a relation between mass and metallicity parallel to the one defined by the SDSS SFGs, but it is offset by ∼ 0.3 dex to lower metallicities.The nitrogen over oxygen (N/O) ratio (Fig. 14) of GPs follows the tendency of SDSS galaxies, except for the low-mass end, where GPs present a higher N/O (see GP22, GP18 and GP20 in Fig. 14).This ratio ranges from log(N/O) = −1.5, −0.85.
Despite the lack of metals in GP galaxies, the metal ratio (N/O) is mostly conserved.One scenario that can explain this would be a massive and recent inflow of metal-poor gas (basically neutral hydrogen clouds) from the HI halo/reservoir of the galaxy.This accretion could dilute the O/H keeping the N/O unaltered (Köppen & Hensler 2005).The fact that their N/O is in most cases the expected for the stellar masses of these galaxies supports this scenario (Amorín et al. 2010).
The proposed scenario of a massive inflow of pristine gas that both sustains low metallicity and increases the SFR is plausible.However, this raises a complex issue regarding the presence of such a process at low redshifts.The origins of this unenriched gas remain an area of inquiry.Nonetheless, empirical data provides insight with evidence from extremely metal-poor galaxies proximal to us (z=0.03) that exhibit substantial amounts of neutral hydrogen, resembling a halo (e.g., Lequeux & Viallefond 1980;Herenz et al. 2023).This suggests that similar circum-stances could exist for GPs, thereby supporting the proposed scenario.
Chemical evolution model predictions (e.g., Mollá et al. 2006;Vincenzo et al. 2016) suggest that the transition between primary and secondary N dominance in the N/O versus O/H plane depends on the SF history, particularly on the star formation efficiency (SFE).A "bursty" galaxy, that is, one having experienced a very recent starburst, will see a quicker increase in N/O than a galaxy with a smoother and longer SF history.Another factor that could elevate the N/O at low metallicity is the initial mass function (IMF), where a higher fraction of massive stars could enhance primary N production at low metallicity.However, we observe that the N/O -stellar mass relation generally holds.For galaxies deviating from this relation (see GP22, GP18, and GP20 in Fig. 14), adjusting the IMF could be a viable solution.

Summary and conclusions
We have presented physical and chemical properties of 24 GPs using MUSE/VLT data cubes.These galaxies are one of the best local analogs of high-redshift galaxies.Thus, their study is fundamental in order to understand the very first epoch of the formation and assembling of galaxies, and in particular the reionization.
For our set of GPs, we wanted to confirm the spatial extension of these sources.To do so, a study regarding the extension of all the sources within the FoV of the MUSE data cube was carried out.We established a criterion to determine whether a GP is resolved based on the comparison between the FWHM, FW 1 10 M, and FW 1 100 M of a stellar-like object and the GP itself.Retrieving seven spatially extended GPs in the core and the low surface brightness region, five GPs extended in the core, and four GPs extended in the low surface brightness region (see Table 2).
We compare the emission line maps and the continuum maps.The only four emission line maps with no circular symmetry are the ones corresponding to GP06, GP07, GP13, and, marginally, GP01.Continuum maps present a much richer structure for all GPs that are proxies for the stellar underlying population.[Oiii]5007 maps are the most extended ones, tracing the low surface brightness regions of ionized gas.
Regarding the ionization structure, Hα/Hβ maps trace dusty regions and zones with low extinction where photons can travel without being absorbed; these maps present very different morphologies for the different GPs.The ionization parameter (as traced by [Oiii]/[Oii]) tends to peak in the center of the galaxies, indicating that the highest ionization is near the star-forming region.GP20, GP06, GP15, and GP22 present the strongest ionization (see Table F The BPT diagrams indicate hot massive stars are confirmed as the main source of ionizing photons.In particular, our GPs are located in the top left part of the diagram, where the most extreme galaxies reside (i.e., those with the lowest metallicity and highest excitation of the ionized gas).The absence of the [FeX]λ6374 in all the spectra discards BHs with masses >> 10 5 M ⊙ contributing to ionizing the gas.
We also produced the integrated spectrum of each GP by integrating the flux of a region defined by the Hα map.The SFRs derived from the luminosity of the Hα line indicate bursts of star-formation with mass-doubling timescales 2 dex lower than common star-forming galaxies.Our study of the ionized gas properties using emission lines indicates low gas metallicity (i.e., 12 + log(O/H) = 7.6 − 8.4), high electron temperatures ranging from 11500 K to 15500 K, and electron densities ranging from 30 cm −3 to 530 cm −3 .The nitrogen over oxygen ratio versus stellar mass of the GPs (see Fig. 14) generally follows the tendency of SDSS galaxies and ranges between log(N/O) = −1.5 and −0.85, whereas these are clearly above the SDSS sequence in the log(N/O) versus 12 + log(O/H) diagram (see Fig. 15), which possibly indicates the presence of an inflow of pristine gas into the galaxies.
We detected the nebular HeIIλ4686 line in the galaxies GP06, GP15 (this particular galaxy being a confirmed LyC leaker (Izotov et al. 2016)), and GP20, indicating the presence of very hard ionizing photons (E > 4Ry).We checked that none of these GPs show WR features in their spectra, which suggests that WR stars are not the main HeII excitation contributors (e.g., Senchyna et al. 2017;Kehrig et al. 2015Kehrig et al. , 2018)).We also note that two of these GPs present among the highest sSFR (> 8 × 10 8 yr −1 ), suggesting that, besides a low metallicity, high sSFR can be a dominant factor to determine the HeII emitting nature of a galaxy (Kehrig et al. 2020;Pérez-Montero et al. 2020).A detailed analysis of the origin of the HeII ionization, which keeps challenging up-to-date stellar models (see e.g., Eldridge & Stanway 2022), is beyond the scope of this work and will be investigated in future work.
A.13)  The only four GPs that present non-circular symmetry in the low surface brightness regions are GP06, GP13, GP07, and marginally GP01 (see Figs.2, A.13, A.7, and A.1).The most extended structures in the maps are the ones corresponding to the [Oiii]5007Å line and, in lesser terms, to the Hα line.The high intensity of these lines allows us to trace the low surface brightness regions corresponding to the ionized gas in the outer parts of the galaxies.Dimmer emission lines (e.g., Hβ, [Nii]6584Å, [S ii]6716Å, [S ii]6731Å, [Oi]6300Å, and [Oiii]4363Å) present less extension and more circular symmetry.

Fig. 2 :
Fig. 2: GP06 emission line maps.The black line in the bottom left corresponds to a distance of 10 kpc.The circle in the bottom right is the representation of the seeing.The peak in the Hα emission is represented by the black point in the middle.The green and black contour represents the 3 σ sky level of the continuum map corresponding to the same galaxy.All maps presented in this work have the previously mentioned features except for the continuum maps, which do not present the continuum contour.Black contours are regions with different signal = k × σ sky levels (with k = 10, 100), and all spaxels represented in all maps are above 3 σ sky .The red contour indicates the FWHM of the map.The same contours are shown in each emission line map and continuum map in this study.

Fig. 3 :
Fig. 3: Continuum map of GP06.Integration is done between 4154Å and 7026 in rest frame.
we show the map of this ratio for GP13, which is the most extended galaxy in our sample.High values of [Oiii] / [Oii] trace galactic regions with high ionization.These regions tend to be near the center of the galaxy where the Hα peaks and the main ionizing sources are also located.The highest values of [Oiii] / [Oii] are found in GP22 and GP15 (see Figs. B.22 and B.15), reaching a value of 6.5.The ratios used in the BPT diagrams are [Oiii]λ5007/Hβ ([Oiii]/Hβ), [Nii]λ6584/Hα ([Nii]/Hα), [S ii]λλ6716, 6731/Hα ([S ii]/Hα), and [Oi]λ6300/Hα ([Oi]/Hα).Examples of maps corresponding to these line ratios are displayed in Fig. 6.

2
See https://astro.ucla.edu/wright/CosmoCalc.html.Most of the stellar masses (M ⋆ ) were reproduced from Izotov et al. (2011); nevertheless, Izotov et al. (2011) did not calculate the mass for all the galaxies in our sample, and so the Cardamone et al. (2009) value was used instead.Izotov et al. (2011) recalculated the masses of the GPs and obtained systematically lower values.Their values are lower because in fitting the SED, they subtracted the contribution from gaseous continuum emission.No errors were listed in the original studies.The corresponding specific SFRs were derived as follows: sS FR = S FR/M ⋆ .

Fig. 10 :
Fig. 10: Integrated spectrum of GP06.The bottom image is a zoomed-in view of the y-scale of the top image. F.2.

Fig. 11 :
Fig. 11: Stellar mass versus sSFR.Green points are the GPs presented in this work, light green are GPs with HeII emission.The green star represents the median of the GPs, corresponding to a stellar mass of 2.6 × 10 9 M ⊙ , an sSFR of 12 Gyr −1 , and thus a mass-doubling timescale of 47 Myr.Black points represent a set of local galaxies (z<0.03)from the Califa survey (Catalán-Torrecilla et al. 2015) (Torrecilla15).The black star represents the median of this sample, it corresponds to a stellar mass of 3.2 × 10 10 M ⊙ , a sSFR of 0.12 Gyr −1 and thus a mass-doubling timescale of 8.3 Gyr, which is on the order of the age of the Universe.We also show the sSFR-mass relation at a variety of redshifts from Tasca et al. (2015) (dash-dot line, T15), Karim et al. (2011) (dashed line, K11), and Whitaker et al. (2012) (dotted line, W12).
[Oiii]  using the[Oiii]λ4363/[Oiii]λ4959, 5007 line ratio and the values of T e corresponding to[Oii]  from the empirical relation between[Oii]  and[Oiii]  electron temperatures given byCampbell et al. (1986).We obtained the electron densities, from the [S ii]λ6716/[S ii]λ6731 line ratio.The oxygen ionic abundance ratios, O + /H + and O 2+ /H + , were derived from the [Oii]λ3727 and [Oiii]λλ4959, 5007 lines, respectively, using the corresponding electron temperatures.The total oxygen abundance is assumed to be O/H = O + /H + + O 2+ /H + .The nitrogen ionic abundance ratio, N + /H + , was calculated using the [Nii]λ6584 emission line and assuming T e [Nii] ∼ T e [Oii]; the N/O abundance ratio was computed under the assumption that N/O = N + /O + , based on the similarity of the ionization potentials of the ions involved.All this was computed by implementing the Pyneb code(Luridiana et al. 2015).We calculated the final errors in the derived quantities by error propagation and taking into account errors in flux measurements.Furthermore, for the subset of galaxies that do not present the [Oiii]4363Å line, we used the HII-CHI-Mistry code (Pérez-Montero 2014), which calculates oxygen and nitrogen over oxygen abundances (and errors) without this line.The values corresponding to all these gas properties are shown in TableF.2.

Figures
Figures 13, 14, & 15 display the oxygen abundance versus stellar mass, nitrogen over oxygen versus stellar mass, and nitrogen over oxygen versus oxygen abundance, respectively.A comparison is made between ∼ 200000 star-forming galaxies (SFGs) from SDSS (Duarte Puertas et al. 2022) and the GPs.For a given stellar mass, GPs show lower oxygen abundances than the bulk of the SDSS galaxies.This is in accordance with the result fromAmorín et al. (2010), where GPs follow a relation between mass and metallicity parallel to the one defined by the SDSS SFGs, but it is offset by ∼ 0.3 dex to lower metallicities.The nitrogen over oxygen (N/O) ratio (Fig.14) of GPs follows the tendency of SDSS galaxies, except for the low-mass end, where GPs present a higher N/O (see GP22, GP18 and GP20 in Fig.14).This ratio ranges from log(N/O) = −1.5, −0.85.Despite the lack of metals in GP galaxies, the metal ratio (N/O) is mostly conserved.One scenario that can explain this would be a massive and recent inflow of metal-poor gas (basically neutral hydrogen clouds) from the HI halo/reservoir of the galaxy.This accretion could dilute the O/H keeping the N/O unaltered(Köppen & Hensler 2005).The fact that their N/O is in most cases the expected for the stellar masses of these galaxies supports this scenario(Amorín et al. 2010).The proposed scenario of a massive inflow of pristine gas that both sustains low metallicity and increases the SFR is plausible.However, this raises a complex issue regarding the presence of such a process at low redshifts.The origins of this unenriched gas remain an area of inquiry.Nonetheless, empirical data provides insight with evidence from extremely metal-poor galaxies proximal to us (z=0.03) that exhibit substantial amounts of neutral hydrogen, resembling a halo (e.g.,Lequeux & Viallefond 1980;Herenz et al. 2023).This suggests that similar circum-

Fig. 14 :
Fig. 14: Nitrogen over oxygen versus stellar mass.GPs are the green points.Black points correspond to ∼ 200000 galaxies from Duarte Puertas et al. (2022).
.2 and Figs.B.15 and B.22).The [Oiii]/Hβ, [Nii]/Hα, [S ii]/Hα, and [Oi]/Hα ratios are studied in maps and in the BPT diagrams.These indicate a tendency for higher excitation of the gas in the center of the galaxy (i.e., higher ratios of [Oiii]/Hβ and lower ratios of [S ii]/Hα and [Oi]/Hα close to the Hα peak).Still, the [Oiii]/Hβ and [Nii]/Hα ratios do not present much spatial variation (with a maximum difference of 0.14 dex in all GPs, except for GP06 reaching over 0.4 dex), indicating uniformity in the gas excitation and metallicity.[S ii]/Hα and [Oi]/Hα trace the boundaries of the ionized gas, and they present their lower values close to the center of the galaxies, suggesting a blister-type morphology (e.g., GP06, GP08, GP10, GP13, and GP23).

Fig
Fig. B.7: Line ratio maps for GP07.Spaxels > 20 kpc to the west from the Hα peak are not reliable due to sky contamination.

Table 1 :
Name, position, redshift, and observation data.The names shown in the first column are the ones adopted in this work.SDSS name refers to the name of each galaxy in SDSS.Exposure time, seeing and date are taken from the ESO archive.Positions and redshifts are taken from the cube header.Each seeing is referring to the corresponding night of observation.

Table 2 :
Extension of stellar-like sources and GPs.Column (1): Name of the galaxy; bold text indicates that the galaxy satisfies the extension criteria.Column (2): Number of SDSS stars in the FoV.Column (3): FWHM of stellar-like sources.