© The Authors 2025

1. Introduction

Recent studies suggest that, at high redshift, the population of star-forming galaxies is characterized by an increased fraction of extreme-emission-line galaxies (EELGs), which have large equivalent widths (EWs) in their rest-optical emission lines (typically defined as EW ≳ 750 Å for [O III]λ5007 or Hα). This fraction increases from < 5% at z < 2 to ≳50% at z ≥ 6 (from spectroscopic studies, e.g., Boyett et al. 2022a,b, 2024; Matthee et al. 2023, while photometric studies find lower yet still elevated factions of ≳20%; e.g., Endsley et al. 2023, 2024). This redshift evolution is justified by an increase in young galaxies with bursty star-formation histories (SFHs) and is consistent with observations showing that high-z galaxies have, on average, bluer UV slopes, β_UV (indicating young stellar populations with low dust-attenuation values; e.g., Bouwens et al. 2014; Bhatawdekar & Conselice 2021; Nanayakkara et al. 2023; Cullen et al. 2023; caputi2024), and larger ionizing photon production efficiencies, ξ_ion (indicating bursty systems characterized by hard ionizing spectra; e.g., Tang et al. 2019, 2023; Simmonds et al. 2024; Harshan et al. 2024; Caputi et al. 2024). When focusing at redshifts of z ≥ 6, galaxies with lower rest-UV luminosities show weaker line emission, especially in [O III]+Hβ, possibly because of lower metallicities (Endsley et al. 2024). Alternatively, this EW trend may be influenced by an anti-correlation between UV luminosity and the Lyman-continuum (LyC) escape fraction. The latter explanation suggests that low-luminosity galaxies may be major contributors to cosmic reionization, in line with two other pieces of observational evidence: (i) low-mass galaxies are more efficient producers of ionizing radiation (e.g., Castellano et al. 2023; Simmonds et al. 2024; Atek et al. 2024; Harshan et al. 2024) and (ii) galaxies with extremely blue UV slopes (β_UV ≲ −2.8), requiring escape fractions of > 50%¹, are predominantly low-mass systems (e.g., Chisholm et al. 2022; Topping et al. 2024).

Despite many recent studies highlighting the redshift evolution of galaxy properties, boosted by the advent of the James Webb Space Telescope (JWST), little is known about their sub-galactic-scale structures. It is well established that galaxies are on average smaller at higher redshifts, partly following the same UV size–stellar-mass relation observed at low redshifts (e.g., Shibuya et al. 2015; Neufeld et al. 2022; Morishita et al. 2024; Sun et al. 2024; Ono et al. 2024). The most massive galaxies (log(M_⋆/M_⊙) > 9) seem to show no redshift evolution (in their mass–size relation), as opposed to the strong evolution observed for lower-mass systems (Langeroodi & Hjorth 2023). In general, galaxies at high z are very compact, with sub-kiloparsec effective radii (R_eff) already at z > 4, and R_eff ≲ 500 pc at z ≥ 6, on average (e.g., Langeroodi & Hjorth 2023; Morishita et al. 2024). Despite their compactness, the study of gravitationally lensed samples, which allows the characterization of ∼100 pc scales even at z ≥ 1, reveals that all star-forming galaxies host small clusters and clumps, which in many cases dominate their recent star-forming episodes. As a notable example, in the z = 2.38 system dubbed the Sunburst arc, a significant fraction (40 − 60%) of the rest-UV emission is located in 13 individual star clusters, with radii of 3 − 20 pc (Vanzella et al. 2022), confirming a trend observed in nearby galaxies, where galaxies with higher integrated star-formation-rate surface densities, Σ_SFR, have a larger fraction of their star formation taking place in bound star clusters (e.g., Adamo et al. 2020).

Large samples of lensed galaxies confirm the presence of compact clumps and clusters hosted by galaxies at any redshift within 1 ≤ z ≲ 8 (e.g., Meštrić et al. 2022; Claeyssens et al. 2023, 2024), and up to z = 10 (Adamo et al. 2024); such stellar systems are denser than is typically observed in local galaxies (Messa et al. 2024). Observations of molecular gas at ∼100 pc scales for a very limited sample of galaxies (at z = 1 − 3.5) suggest that the progenitor gas clouds of high-z stellar clumps are also denser than their local counterparts, and that the star formation process itself, converting dense gas into star complexes, is more efficient at higher redshifts (Dessauges-Zavadsky et al. 2019, 2023; Zanella et al. 2024). This trend seems to be driven by the gas-rich nature of high-z galaxies (e.g., Wisnioski et al. 2015; Girard et al. 2018; Dessauges-Zavadsky et al. 2020), as also suggested by models and simulations (e.g., Mandelker et al. 2017; Fisher et al. 2017; Renaud et al. 2021; Fensch & Bournaud 2021; Dekel et al. 2022; Garcia et al. 2023; Renaud et al. 2024; Sugimura et al. 2024). While this is a robust explanation for most of the observed stellar clumps at cosmic noon, z ≲ 3 (Zanella et al. 2019), whether or not it is also valid at higher redshifts, where galaxies are expected to be more affected by major mergers (e.g., Hopkins et al. 2010; Rodriguez-Gomez et al. 2015; Calura et al. 2022; Nakazato et al. 2024; Ceverino et al. 2024), remains unknown.

At low redshifts, the presence of massive and dense young star clusters is associated with observations of strong feedback from their massive stars (e.g., Sirressi et al. 2024), which can be sufficiently powerful to affect the host by creating strong galactic-scale outflows, as seen in M82 (Ohyama et al. 2002; Förster Schreiber et al. 2003) and NGC 253, (Westmoquette et al. 2011) and/or ionized channels that facilitate the escape of ionizing radiation. This is particularly true in dwarf starburst galaxies, such as ESO338-IG04 (Bik et al. 2015, 2018), Haro 11 (Sirressi et al. 2022; Komarova et al. 2024), and SBS0335-052E (Wofford et al. 2021), which are closer analogs of high-z systems, and has been observed directly in the Sunburst arc at z = 2.38 (Rivera-Thorsen et al. 2017, 2019; Pascale et al. 2023). Such observations suggest that compact massive clusters could be the source of the extreme line emission and ionizing production efficiency observed in high-z EELGs. This was also suggested by a few other direct examples, such as the Sunrise arc at z = 6, and the Firefly Sparkle at z = 8.3, both showing star clusters with radii of R_eff < 10 pc associated with large equivalent widths in Hα and [OIII]+Hβ (Vanzella et al. 2023a; Mowla et al. 2024). However, studies at small intrinsic scales and high z are, in most cases, beyond the capabilities of current facilities. In this paper, we exploit the combination of extreme magnification (μ > 15) from gravitational lensing and new observations from JWST to unveil the properties of a Lyman-α emitter at z = 6.145 down to scales of ∼10 pc, and try to identify the powering source of its intense nebular emission. The target system, along with the data reduction, is presented in Sect. 2. Using the new data, the system is analyzed in Sects. 3 and 4, while the main results are discussed in Sect. 5. Finally, in Sect. 6 we present our conclusions and outline future prospects of this work.

Throughout this paper, we adopt a flat Λ-CDM cosmology with H₀ = 68 km s⁻¹ Mpc⁻¹ and Ω_M = 0.31 (Planck Collaboration XVI 2014), the Kroupa (2001) initial mass function, and a solar metallicity of Z_⊙ = 0.018 (Asplund et al. 2009). All quoted magnitudes are on the AB system.

2. Target, observations, and data reduction

2.1. The target system at z = 6.145

The target system of the current study was initially identified as a Lyman-α arc at z = 6.145 in the galaxy cluster field MACS J0416.1–2403 (hereafter MACS0416, Fig. 1). The arc extends for ∼45″ on the sky and is composed by three multiple images (system 2 of Caminha et al. 2017). In the most magnified image, with μ > 15 and covering a region of ∼6 kpc on the source plane (Bergamini et al. 2023; see also Sect. 3.1), Vanzella et al. (2017, hereafter V17) identified 3 rest-UV sources in the Hubble Space Telescope (HST) data at the position of the peak of the emission, separated by < 2″ (two of them have photometry from the Astrodeep photometric catalog of Castellano et al. 2016). Following the nomenclature of Vanzella et al. (2019, hereafter V19), we refer to the three stellar systems as D1, T1 and UT1 (meaning Dwarf 1, Tiny 1 and Ultra-Tiny 1, respectively). V17 also identified a Lyα-emitting knot in the southern part of the arc (named T2, initially without detectable HST counterpart) and another lensed Lyα emission line of an object (named D2) at the same redshift as the ‘main’ arc but with an angular separation corresponding to ∼25 kpc in the source plane, implying it is a distinct system. Also, D2 has a well-detected UV-optical counterpart in the Hubble Frontier Field (HFF) deep photometry. Finally, V19 recognized a pair of sources, T3 and T4, at the edge of the Lyα arc and showing the same colors and dropout signature as D1 and T1, with a z_phot ≈ 6 in the catalog of Castellano et al. (2016), thus probably part of the same environment. The positions of the Lyα emission and of its UV/optical counterparts are shown in the map of Fig. 1. Searching for neutral gas in this system with the Atacama Large Millimeter/submillimeter Array (ALMA), Calura et al. (2021) report a 4σ tentative detection of [CII] emission; the low-luminosity observed in [CII] is possibly due to low-density gas and/or a strong radiation field caused by intense stellar feedback from the stellar sources. Deep observations with the Multi Unit Spectroscopic Explorer (MUSE) on the Very Large Telescope (VLT, Prog.ID 0100.A-0763(A), PI: E.Vanzella) reveal the presence of two other Lyα halos at the same redshift, separated from D1-T1 by ≳50 kpc in the source plane; these additional sources will be presented in a separate publication (Messa et al., in prep.). We name this entire system, consisting of several “islands” at the redshift of z = 6.145 within a ∼100 kpc radius area, the “Cosmic Archipelago”.

Fig. 1. Color image combining HST and JWST observations from the 2023-146 release, which is available at the following link, with over-plotted: (i) the Lyα emission from VLT-MUSE of the Cosmic Archipelago systems, at z = 6.145, as yellow contours; and (ii) the four NIRSpec-IFU pointings of program GO 1908, across the main arc and also covering the compact Lyα halo D2, as green squares. The pointing analyzed in detail in the current publication covers the D1 and T1 systems (darker green square). (iii) The other pointing of GO 1908, denoted with a cyan square and covering the metal-poor system LAP1 at z = 6.63, is presented in a separate publication (Vanzella et al. 2023b).

A first analysis of the physical properties of the D1, T1 and D2 systems, based on HST data, is presented in V17 and V19. Briefly, all systems appear to be young and compact. The light profile of D2 is almost unresolved, leading to an upper limit for the intrinsic effective radius R_eff < 100 pc, after accounting for de-lensing. The large magnification of the T1 system (μ ∼ 25) allows us to constrain an intrinsic effective radius of R_eff < 16 ± 7 pc, only slightly larger than the typical size of individual young star clusters and globular clusters observed in the local Universe (e.g., Brown & Gnedin 2021). A similarly large magnification in D1 allows us to split the source into a compact core (< 13 pc) within a larger (220 pc) structure. While the properties of the D1 system seem to be robustly constrained, with M_⋆ = (2 ± 1)⋅10⁷ M_⊙ (de-lensed value) and an age of 1 Myr (but with 3σ uncertainties spanning up to masses of 15 ⋅ 10⁸ M_⊙ and to ages > 700 Myr), the large photometric uncertainties of T1 and D2 imply large 1σ intervals on the derived physical properties, M_⋆ = 10⁶ − 10⁸ M_⊙ and age = 1 − 700 Myr for both T1 and D2. For what concerns the other sources, only magnitude values or limits were given for UT1, mag_UV = 32.1 and D2, mag_UV > 32.5 (de-lensed values, V17).

The lack of information (or large uncertainties) on the physical properties of the systems that compose the Cosmic Archipelago was filled with new photometric and spectroscopic information in rest-UV/optical from JWST observations. In particular, the compact rest-UV sources, identified in previous HST data, were covered with 4 NIRSpec integral field unit (IFU) pointings (GO 1908; Sect. 2.2), as shown in Fig. 1, while the entire cluster MACS0416 was imaged with NIRCam (Sect. 2.2). In the current work, we consider the data of one pointing only (named as “D1T1”, highlighted in Fig. 1). The analysis of the pointings relative to the T2, T3-T4 and D2 sources, as well as an overview analysis of the entire Cosmic Archipelago system, will be given in forthcoming publications.

2.2. NIRSpec IFU

We analyze the data from the Cycle 1 GO program 1908 (PI: E. Vanzella). The observations consist of five pointings taken with NIRSpec in the IFU mode. These data were obtained applying a small-cycling dithering pattern with eight dithers, using the high-resolution G395H/F290LP grism/filter combination. For each dither we used 22 groups, with one integration per group, for a total integration time on source of 3.6 h.

We processed the data with the version 1.14.0 of the JWST Pipeline (Bushouse et al. 2023), developed by the Space Telescope Science Institute (STScI), and the Calibration Reference Data System (CRDS) context jwst_1230.pmap, indicating the reference files used to calibrate the data. The data reduction consists of three stages, each of them including several steps. Stage 1 applies detector-related corrections to the raw files (i.e., uncalibrated ramps), such as bias and dark subtraction, linearity correction and cosmic rays flagging. The final product of this stage are 2D count rate exposures (rate files), which are the input data for Stage 2. During this stage, important steps such as wavelength, flat-field, and flux calibrations are applied. The final products of this stage are calibrated exposures (cal files), which are then processed in Stage 3 to build the final datacube. We note that, at the end of Stage 2, also calibrated datacubes, corresponding to the eight dithers, are generated. Stage 3 combines the eight calibrated exposures accounting for the sub-pixel spatial shifts due to dithering and creates the final datacube. A further flagging of cosmic rays is also applied in this stage on the eight calibrated exposures.

We added further steps in order to improve the data reduction. In particular, before running Stage 2, we removed the residual 1/f noise presented in the rate files by subtracting the median of each spectral column, which was estimated after applying a sigma-clipping. A similar procedure was applied in other works (Perna et al. 2023; Rauscher 2024; Loiacono et al. 2024). At the end of Stage 2, we excluded all the saturated pixels and the ones with a bad flat-field solution from the build of the final datacube by updating their quality flag (“DQ” extension of the cal files). We also removed, at the end of Stage 3, possible residual cosmic rays by filtering out all the spikes persisting for less than four channels (corresponding to < 150 km s⁻¹ at 4.0 μm), which is the typical width of these features. We carefully inspected the voxels corresponding to the science targets before and after applying this step to check that none of the emission lines is affected by the algorithm. These last customized step, combined with the outlier detection implemented in the pipeline, turned out to be ineffective in removing part the of the spikes in the datacube (see Appendix A). Therefore, instead of using the datacube produced at the end of Stage 3 for the analysis, we developed a different approach to remove these features.

We used the eight datacubes generated at the end of Stage 2 for producing the final combined and clean cube. In particular, two main steps were implemented: (1) defects and spikes associated with the detector coordinates were identified and masked by calculating the “persistence” parameter (R) as the ratio between the median signal on the spaxel at the given physical coordinates and its one-sigma median deviation (both derived from the eight partial cubes) and (2) we combined the eight cubes on the WCS coordinates after aligning them according to the dithering pattern used in the observations (small cycling). In this step a 3-sigma clipping procedure was applied to clean from outliers among the cubes. In this step, along with the 3-sigma clipped average of each pixel that contributes to the final scientific datacube, the median deviation is stored in a separate datacube which collects the pixel errors from the combination of partial cubes. This error cube is then used to extract one-dimensional 1σ error spectra by propagating in quadrature the pixel errors from the datacube within the same mask used for spectral extraction². The alignment of the eight cubes was performed in the discrete grid of pixel without applying any rebinning, implying an accuracy better than 1/2 pixel (with the adopted small cycling dithering pattern, the average accuracy is 0.24 pixel). As aforementioned, the R quantity accounts for the defects associated with the detector, while the other cleaning by sigma-clipping erases possible outliers when combining the aligned cubes. We explored different values of R, from 2 (aggressive) to 10 (conservative) and find an optimal compromise with R = 5. However, the choice of R depends on the brightness and size of the targeted sources, along with the adopted dithering scheme. In fact, sources which spatially extend more than the dithering pattern amplitude can generate a smooth (persistent) positive signal on the detector coordinates before aligning the cubes. We verified that this signal was not recognized as possible spurious pixels (by high “R”). Moreover, to avoid this behavior in the construction of our masked pixels, a background was subtracted for each slice after calculating a moving median over large windows of 4 × 4 spaxels (we note that this is not the final background of the cube, it is needed to construct the mask). This process produces a mask for each slice on the detector coordinates, in which all pixels with R > 5 are flagged as nan.

The background was initially calculated as a single scalar median value per slice after masking the known sources (identified on NIRCam images); this method leaves a second-order residual background, spatially varying within the IFU (see Appendix A). The latter is calculated over adjacent 30 slices (≃200 Å) blueward and redward the actual slice (masking only the slices containing the known lines) and subtracted from the cube. This method can be applied in this specific case being the continuum of the targets too faint and not detected in our cubes.

The variation of the point spread function (PSF) across the wavelength range covered by the observations could impact the line ratios derived in this work, especially when distant lines are compared (e.g., Hα and Hβ; see also Venturi et al. 2024). To account for this possible bias, we implement a correction to the cube based on the PSF variation described by Eq. (2) of D’Eugenio et al. (2024); this is one of the few attempts to characterize the PSF variations in NIRSpec, based on the observation of a relatively bright star. We compare, in Appendix A, the main line ratios derived with and without this PSF correction; the relative difference among the two is below ∼7%.

The resulting combined cube shows fluxes compatible with the ones obtained using the standard reduction pipeline (though the latter shows several systematics and artifacts along the entire wavelength). We also verified the flux calibration by comparing the extracted line fluxes from our targets (the most prominent lines Hβ+[OIII]4959, 5008 and Hα) with the NIRCam photometric excesses which arise from the same areas. The photometric and spectroscopic line fluxes extracted for D1 (see Fig. 2 and Appendices C, D) are in agreement, with the photometric inferred ones lying within 20% from the spectroscopic measurements. We show in Fig. 2 a map of the Hα and OIII emission, along with complete spectra (and 1σ error) extracted for the D1 and T1 regions.

Fig. 2. Imaging and spectroscopy of the D1T1 system. Top panels: Cutouts showing the observed systems in the sum of NIRCam SW (F115W, F150W, F200W, left) and LW (F277W, F356W, F410M, F444W, center) filters and a NIRSpec-IFU map with the sum from Hα, Hβ, and [OIII] line emission (right). White ellipses mark the apertures used to extract photometry of the three main regions (labeled in the left panel). The white dashed aperture is used only to compare NIRCam to NIRSpec-derived quantities for the T1 region (see Sect. 3.4). The right panel also contains the pixel masks (blue and red contours) used to extract the spectra shown in the bottom panels. Middle panel: Spectra for D1 and T1 regions (in blue and red, respectively). No stellar or nebular continuum is detected in the spectra. The gray shaded area marks the wavelength range falling in the NIRSpec detectors’ gap. The insets in the bottom panels are zoom-ins around the detected lines and display also the uncertainty of the spectra (as wide shaded lines). The best-fit center of the lines is displayed as vertical lines and reveals a small shift in redshift between the two regions (corresponding to a velocity shift of ∼42 km/s; see also Fig. 5, Sect. 3.4 and Appendix D).

2.3. NIRCam imaging

The NIRCam products combine all observations on MACS0416 from the Prime Extragalactic Areas for Reionization and Lensing Science program (PEARLS, PID 1176; Windhorst et al. 2023) and the CAnadian NIRISS Unbiased Cluster Survey (CANUCS, PID 1208; Willott et al. 2022). Both programs utilize eight passbands, namely, F090W, F115W, F150W, and F200W in the short wavelength (SW) channel and F277W, F356W, F410M, and F444W in the long wavelength (LW) channel, for a total of exposure times between 15 100 and 17 700 s per filter. Because of the redshift of our target, we do not consider the observations in F090W throughout this work. Although the main systems are detected in F090W, at their redshift the filter encompass the LyC region, with intergalactic medium (IGM) absorption, and the resulting photometry would therefore be difficult to model. The cluster was centered on the module B of NIRCam in both programs.

The reduction of these data followed the same procedures as described in Yan et al. (2023). We retrieved the data from the Mikulski Archive for Space Telescopes (MAST). Our reduction started from the so-called Stage 1 “uncal” products, which are single exposures after the Level 1b processing through the default JWST pipeline (Bushouse et al. 2023). Our further reduction was based on the JWST pipeline version 1.9.4 in the calibration context of jwst_1063.pmap. The astrometry of each single exposure was registered to the public HFF products. All the individual exposures in each band were then projected onto the same grid and combined. For ease of photometry incorporating the HST HFF products, these NIRCam stacks were made at the pixel scale of 0.04″ to match that of the HST images. In addition, we create stacks for the SW filters at a pixel scale of 0.02″ in order to leverage the angular resolution of the instrument and study the smallest substructures of the system. The stack of the observations from the SW and from the LW filters are shown in Fig. 2.

To perform multiband aperture photometry (Sect. 3.2) we create images PSF-matched to F444W, the filter with the broadest PSF; details of the matching technique are given in Merlin et al. (2022). The PSF models are observationally-derived from stars in the FoV of the observations.

3. Analysis of the main regions

NIRCam imaging reveals three main regions, (D1, T1, and UT1 following the naming of V17 and V19), characterized by substructures, the most evident being the bright compact “cores” of D1 and T1 (D1_core and T1_core), and local peaks in the rest-UV emission (see Fig. 2); we name the latter by adding lowercase letters to the name of the main region they belong (e.g., D1a; see Fig. 3). In the first part of the analysis, we characterize photometrically and spectroscopically, the three main systems as a single entity each; their subregions are studied by analysing individually the properties of cores and peaks in Sect. 4. A complete NIRCam view of all main and subregions discussed in this section is given in Appendix B.

Fig. 3. Positions of sub-regions in the image and source planes. Left panel: Labels of the main rest-UV peaks and subregions of the D1-T1-UT1 systems on the stack of the SW filter observations. The inset shows a stacking of all filters where the three main regions (D1, T1 and UT1) are delimited by the red contours, while larger contours enclosing also the low surface-brightness emission are shown in green. Right panel: Mapping of the main and subregion positions on the source plane at z = 6.145 obtained from the best-fitting lens model by Bergamini et al. (2023). The compact cores (D1core, T1core and UT1a) are marked by white filled circles, while the other peaks of UV emission are marked by white “X” symbols. A 1 kpc (0.17″) reference scale is given in the top left corner. The FWHM of the F200W is also given as reference. The location of T1core coincides with the peak of the Lyα emission introduced in Sect. 1. The bridge region refers to the region of faint line emission discussed in Sect. 3.5.

3.1. Lens model and source plane reconstruction

In order to derive the intrinsic properties of our targets, we need to rely on a lensing model. In this work, we refer to the parametric lens model developed by Bergamini et al. (2023). This benefits from deep integral field VLT/MUSE observations (Vanzella et al. 2021) to obtain the multiple image and cluster member samples and to measure galaxy stellar velocity dispersions. In detail, the model is based on 237 spectroscopically confirmed multiple images, spanning the redshift range 0.9 < z < 6.6, and is characterized by a global precision of 0.43″ in reproducing the observed positions of these images. The subhalo mass component of the model counts 213 cluster member galaxies, whose total mass is accurately determined by exploiting the additional lensing-independent kinematic information obtained from the measured central stellar velocity dispersion values of 64 galaxies. Using this model, we estimate for the D1, T1, and UT1 systems the magnification factors on the model predicted positions and the associated 68% confidence level intervals, computed from 500 different realizations of the model obtained by randomly extracting samples of free-parameter values from the final Markov chain Monte Carlo (MCMC) chain.

The best-fitting lens model is also used to de-lense onto the source plane at z = 6.145 (see Sect. 3.4 for the redshift measurement) the observed positions of the substructures composing the system D1-T1 (Fig. 3). The source-plane view of the system reveals how the three “macro” components of this region are close to each other, with separations of ∼700 pc between the cores of D1 and T1 and between T1 and UT1. A similar distance exists between D1_core and D1a or D1d, the two “extremities” of the D1 region; this is caused by D1 system extending in the direction transversal to the shear of the gravitational lens. The T1 system is well-fitted, in the image plane, by a Sersic profile with an effective (i.e., half-light) radius R_eff = 0.05″ along the tangential direction of the magnification (while it is unresolved in the transversal direction); this observed radius corresponds to R_eff = 18.6 pc in the source plane, considering the tangential magnification of the lens model at the position of T1 (μ_tan = 15.8). Overall, all the subregions of D1, T1 and UT1 are enclosed in a region with 1 kpc radius. This radius is slightly larger than the average value found for galaxies with similar magnitudes (R_eff ∼ 300 pc for galaxies with M_UV ≈ −18 mag at 5 < z < 7, Morishita et al. 2024), yet is well within the observed scatter. We cannot easily discern if the D1, T1 and UT1 systems are subregions of a common (unobserved) galactic structure or instead are different structures in the course of merging to form a single one, as predicted by cosmological simulations (Calura et al. 2022). The possible nature of the system(s) is further discussed in a forthcoming publication (Messa et al., in prep.).

3.2. Photometric properties

We derive photometry from large elliptical apertures (see Fig. 2); the same apertures are also considered for spectral extraction in the IFU data, using the spaxel masks shown in the top-right panel of Fig. 2. Those apertures and masks are chosen to enclose at least 90% of rest-UV-optical and line emission. Larger apertures would add more noise, especially in the IFU data. The T1 region resides close to the border of the NIRSpec field of view and the photometric aperture extends beyond the IFU coverage (solid ellipse in Fig. 2). However, we make use of a smaller aperture (dashed ellipse in the top panels of Fig. 2) tracing the mask used in the IFU when we need to compare line measurements between NIRCam and NIRSpec (Sect. 3.4).

We perform aperture photometry on images PSF-matched to the F444W filter (see Sect. 2.3). Sky background is estimated for each source as a sigma-clipped median in a region surrounding its aperture and is subtracted from the aperture flux. We perform aperture correction, and assuming PSF-like sources we derive −0.13, −0.25 and −0.55 mag for D1, T1 and UT1, respectively³. The final fluxes are corrected for galactic reddening (which is small, ranging from 0.033 mag to 0.005 mag from F115W to F444W, respectively). The resulting photometry is shown in Fig. 4. We have converted the observed magnitudes in F115W into absolute rest-frame FUV magnitudes and, in turn, into star formation rate (SFR) values, following the prescription of Kennicutt & Evans (2012). By considering the reference magnification value for each of the systems (see Appendix C), we reach the intrinsic (i.e., de-lensed) values SFR_UV: 0.35 ± 0.02, 0.12 ± 0.01 and 0.03 ± 0.01 M_⊙ yr⁻¹ for D1, T1 and UT1, respectively⁴. Those are typical values observed for compact (≲100 pc) star-forming regions at similar redshifts (e.g., Messa et al. 2024), and are larger than the average SFRs observed for their local counterparts in main-sequence galaxies (SFR_UV < 0.01 M_⊙ yr⁻¹ at z ≲ 1, e.g., Kennicutt et al. 2003). On the other hand, they are similar to the rates inferred in compact regions within the most active local blue compact dwarf galaxies (e.g., Calzetti et al. 1997; Annibali et al. 2003; Messa et al. 2019).

Fig. 4. Photometry of the main regions, D1, T1, and UT1 (black circles), obtained via aperture photometry (apertures shown in Fig. 2). The best-fit UV beta slopes are shown as blue lines (and relative blue-shaded uncertainties) and reported at the top of the panels. The results of the broadband-SED fitting are over-plotted (red empty squares and red solid line: spectrum of the maximum-likelihood model, with red-shaded uncertainty spectrum) and reported at the bottom of the panels.

All three sources have steep rest-UV continuum slopes (β_UV = −2.4 ± 0.1, −2.7 ± 0.1 and −2.5 ± 0.2 for D1, T1 and UT1 respectively), bluer than the average slope of galaxies at similar redshifts (β_UV ≈ −2, e.g., Castellano et al. 2012; Bouwens et al. 2014; Jiang et al. 2020; Topping et al. 2024; Nanayakkara et al. 2023; Weibel et al. 2024) and hinting at the presence of young stellar populations (age ≲ 10 Myr) with sub-solar (Z < 0.5 Z_⊙) metallicity (e.g., Bolamperti et al. 2023). Further evidence of the young ages of D1 and T1 systems is given by the emission lines observed in the spectra of the sources (Sect. 3.4), reflected also in the photometric jumps observed in the F356W (from [OIII]+Hβ) and F444W (from Hα) filters. On the other hand, the lack of strong lines in UT1 suggests an age ≳10 Myr for that system (see also the next section).

3.3. Spectral energy distribution fitting

By fitting the broadband⁵ spectral energy distribution (SED) we derive the main physical properties of the sources (age, stellar mass, extinction, metallicity, logU), which are reported in Appendix C. The SED fitting is performed via the publicly available code Bagpipes (Carnall et al. 2018). We use as reference the Binary Population and Spectral Synthesis (BPASS) stellar models (v2.2.1, Eldridge et al. 2017; Stanway & Eldridge 2018), but we also test the “standard” Bruzual & Charlot (2003) models⁶ (Appendix E). We consider a Kroupa (2001) IMF with an upper-mass limit at 300 M_⊙; the stellar models also include Cloudy (Ferland et al. 2013) output to account for nebular emission. Due to the limited number of bands available we fix the star formation history (SFH) to an exponentially declining model (described by the decline timescale τ). Extinction and metallicity are left as free-parameters; the indications coming from NIRSpec line analysis (Balmer decrement and metallicity index) are compared to the SED results a posteriori (see Sect. 3.4). We assign a flat uninformative prior to extinction (A_V), log(M_⋆) and logU, while for age, τ and metallicity a flat prior in logarithmic space is used.

The spectra of the maximum-likelihood models (and associated uncertainties from posterior distributions; see Appendix C) are over-plotted to the photometric data in Fig. 4. T1 is the youngest system with a (mass-weighted) age of only $1_{-0}^{+3}$ Myr and a present-age SFR = 17.3 M_⊙ yr⁻¹, while the best-fit ages for D1 and UT1 are both > 30 Myr. However, while the star formation in D1 is characterized by a slow decline (τ ≈ 500 Myr), resulting in a present-age SFR = 0.34 M_⊙ yr⁻¹ (consistent with SFR_UV), the SFR in UT1 declines quickly (τ = 12 Myr) resulting in SFR ≈ 0 M_⊙ yr⁻¹ at present-age. The intrinsic (i.e., de-lensed) masses of the systems span from 10⁶ M_⊙ in T1 to 2.5 ⋅ 10⁷ M_⊙ in D1. Best-fit extinctions (A_V ≤ 0.08 mag) and metallicities (Z < 0.3 Z_⊙) are low and consistent with the values derived from the line ratios, within uncertainties (see Sect. 3.4).

We note that for T1 the best-fit ionization parameter is (within uncertainties) in the range log U = −1.8 to −1.1, larger than in star forming (and starburst) galaxies, where typical values are log U ≤ −2, (e.g., Yeh & Matzner 2012). High values of the ionization parameter, sometimes observed in dense (individual) HII regions (e.g., Snijders et al. 2007; Indebetouw et al. 2009), are associated with large equivalent widths for the nebular emission lines (e.g., Simmonds et al. 2024), which is the case also for T1 (see the next section), and are typical of systems with large SFR densities (e.g., Reddy et al. 2023a,b). In addition, despite the very young age, lack of extinction, low-metallicity and high ionization of the system, the best-fit model struggles to reproduce the steep UV slope of T1; we refer to Sect. 5.1 for a discussion of this result.

3.4. Spectroscopic properties

In parallel to NIRCam photometry and SED-fitting, we create a “master” spectrum of D1 and T1 regions from the IFU data, by summing all the spaxels within the selected apertures (Fig. 2). In the resulting spectra no stellar continuum is detected, but we identify, both for D1 and for T1, Hα, Hβ, [O III]5007 and [OIII]4959 lines with ≥5σ confidence⁷ (bottom panels of Fig. 2). The Hγ line is detected in the D1 mask, while it has a very low signal-to-noise ratio in T1 (S/N < 3). The EW of the main lines are reported in Appendix D. The line emission map in Fig. 2 (right) shows a low-luminosity bridge between these two regions, slightly displaced from the position of UT1; we therefore consider the latter as lacking emission lines (as also suggested by photometry, Fig. 2). The faint line emission from this bridge region is shown and discussed in Sect. 3.5. From the same map we note how, in D1, the rest-UV peak (from the NIRCam imaging) is spatially displaced from the peak of nebular emission; this difference is further investigated in Sect. 4.

We start by testing the consistency of the brightest lines’ fluxes between NIRCam and NIRSpec, assuming in both cases that the stellar continuum below the lines (un-detected in NIRSpec) is traced by the medium band photometry in F410M. In the case of D1 we recover F_Hα = (3.8 ± 0.2)⋅10⁻¹⁸ erg s⁻¹ cm⁻² for the Hα line, while putting together Hβ with [OIII]4959,5007 lines we recover a total flux F_{Hβ, [OIII]} = (8.8 ± 0.3)⋅10⁻¹⁸ erg s⁻¹ cm⁻²; these values are consistent within uncertainties with the values derived from the flux excesses in the NIRCam F444W and F356W filters, F_F444W = (3.4 ± 0.7)⋅10⁻¹⁸ and F_F356W = (9.4 ± 0.7)⋅10⁻¹⁸ erg s⁻¹ cm⁻², respectively. These line fluxes convert to (rest-frame) equivalent widths EW(Hα) = 455 ± 26 Å and EW(Hβ,[O III]) = 605 ± 24 Å. These EWs are in line with the average values found in z ∼ 6 galaxies (mostly ranging between 300–1000 Å and 200–900 Å, for Hα and [O III]+Hβ, respectively, e.g., Endsley et al. 2023; Boyett et al. 2024). For what concerns T1, the line fluxes derived from the IFU mask (F_Hα = (3.0 ± 0.2)⋅10⁻¹⁸ erg s⁻¹ cm⁻² and F_{Hβ, [OIII]} = (6.7 ± 0.3)⋅10⁻¹⁸ erg s⁻¹ cm⁻²) are lower than what is estimated from apertures covering the same region in NIRCam (F_F444W = (4.2 ± 0.6)×10⁻¹⁸ erg s⁻¹ cm⁻² and F_F356W = (8.2 ± 0.6)×10⁻¹⁸ erg s⁻¹ cm⁻²). This large difference is not driven by the difference in PSF between instruments; while the IFU fluxes are not aperture-corrected the mask used to extract the flux is considerably larger than the FWHM of NIRSpec (∼0.2″, i.e., 2 pixels). We consider the proximity of the T1 region to the border of the IFU detector (where noise and defects are more prominent than in the rest of the detector) as the possible main cause of this discrepancy. For this reason, we consider the NIRCam-derived fluxes more robust to derive equivalent widths in this system. The NIRCam photometry in F444W and F356W leads to rest-frame EW(Hα) = 2180 ± 630 Å and EW(Hβ, [OIII]) = 2810 ± 860 Å, making T1 standing out as a system with extreme equivalent widths when compared to “average” galaxies at similar redshift. The location of T1 coincides with the peak of the Lyα halo covering the entire region (V19, see also Fig. 1).

The Hα luminosities (L_Hα) can be used as tracers of the ionizing photon production (Q_H⁰) and, if compared to rest-UV luminosities, give the ionizing photon production efficiency, ξ_ion, for example, following Bouwens et al. (2016), Emami et al. (2020):

$$\begin{array}{c}\hfill {\xi }_{\mathrm{ion}}=\frac{Q_{{\mathrm{H}}^0}}{L_{\nu }(\mathrm{UV})};Q_{{\mathrm{H}}^0}[{\mathrm{s}}^{-1}]=\frac{L_{\mathrm{H}\alpha }[{\mathrm{erg}\mathrm{s}}^{-1}]}{1.36\times {10}^{-12}},\end{array}$$(1)

where L_ν(UV) is the UV-continuum luminosity density (per unit frequency) around 1500 Å. We use the photometry in F115W (rest-frame pivotal wavelength 1615 Å) to trace the UV luminosity and we derive log(ξ_ion/erg⁻¹ Hz) = 25.2 ± 0.1 and 25.7 ± 0.1 for D1 and T1, respectively; these values are in line with the efficiencies measured in extreme emission line galaxies (EELGs), z = 7 − 9 (e.g., Tang et al. 2019, 2023) and in the most metal-deficient compact star-forming local galaxies (Izotov et al. 2024). Two factors that may affect these measures are (i) the presence of extinction, making L_ν(UV) fainter, thus increasing ξ_ion, and (ii) the escape of ionizing radiation, f_esc > 0, in which case Q_H⁰, and consequently ξ_ion would be underestimated); in our calculation, Eq. (1) would become:

$$\begin{array}{c}\hfill {\xi }_{\mathrm{ion}}=\frac{Q_{{\mathrm{H}}^0,\mathrm{obs}.}{10}^{-0.4·E(B-V)·k_{{\mathrm{H}}_{\alpha }}}}{(1-f_{\mathrm{esc}})}·\frac{{10}^{0.4·E(B-V)·k_{\mathrm{UV}}}}{L_{\nu }{(\mathrm{UV})}_{\mathrm{obs}.}}.\end{array}$$(2)

We use the ratios of the Balmer lines to derive average extinctions. The low Hα/Hβ ratios of 2.8 ± 0.2 (D1) and 2.8 ± 0.3 (T1) are consistent with no or very little extinction (A_V < 0.1 mag). These results are consistent with the outcomes of the SED fitting. The Hγ line is detected in D1; the ratio to the other Balmer lines (e.g., Hα/Hγ = 6.0) is still consistent with A_V ∼ 0 mag. The very low dust extinction is also supported by the presence of prominent Lyα emission associated with this system and along an arc-like shape (V19). The possible presence of escaping ionizing photons is much harder to assess and is discussed in Sect. 5.1. One possible reason for the extreme EW and ξ_ion values measured is that lensing is zooming into compact star-forming regions (especially in T1, where R_eff < 20 pc; see Sects. 3.1 and 4.2), which are supposedly the main driver of the ionizing photon production even in un-resolved galaxies. A deeper focus on the intrinsic small scales features of D1 and T1 systems is given in Sect. 4. On the other hand, efficiencies above log(ξ_ion) > 25.5, as found in T1, are hardly reached by “standard” stellar populations (e.g., Stanway & Eldridge 2023).

In order to estimate gas metallicity, we employ the indirect metallicity index based on the R3 = ([O III]λ5007/Hβ) strong-line method, widely used in the literature (e.g., Pagel et al. 1979; Maiolino & Mannucci 2019; Nakajima et al. 2023; Katz et al. 2023; Maseda et al. 2023; Curti et al. 2024; Sanders et al. 2023). For the D1 and T1 systems, we have R3 = 4.1 ± 0.3 and 3.8 ± 0.4. This strong-line method is dependent on the ionization parameter, which increases the scatter of the conversion factor (Izotov et al. 2021a). In order to account for this, we use the calibration given by Eq. (1) in Nakajima et al. (2022), estimated for samples at different EW(Hβ). More specifically, given that we observe EW(Hβ) = 89 ± 7 Å for D1 and EW(Hβ) = 351 ± 68 Å for T1, we use their “Small EW” calibration (< 100 Å) in the first case, and the “Large EW” calibration (> 200 Å) in the second. Under these assumptions, we derive $Z=0.{14}_{-0.06}^{+0.11}{\mathrm{Z}}_{\odot }$ and $Z=0.{05}_{-0.02}^{+0.02}{\mathrm{Z}}_{\odot }$ for D1 and T1, respectively. Those values are consistent with the ranges derived from the SED fitting in Sect. 3.3. Other lines commonly used in the literature to study the metallicity of galaxies at high-z are, among others, [OIII]λ4363 and [NII]λ6585; the latter is not detected in any of our spectra, and the upper limit we can derive is quite uninformative, implying a metallicity Z < 60% Z_⊙ (following the empirical calibrations of Curti et al. 2020; Nakajima et al. 2022, see Appendix F). We observe a faint signal (S/N ∼ 3) at the expected wavelength of [OIII]λ4363 in the spectrum of D1 (more details are given in Appendix F). If true, this tentative detection would imply that the ionized gas in D1 is characterized by either high temperatures (T_e ≫ 10⁴ K) or high densities (n_e > 10⁴ cm⁻³), or by a combination of both; deeper spectra would be necessary to confirm this detection.

The line profiles of T1 are slightly resolved, with velocity dispersion σ_[OIII]5007 = 20.0 ± 8.6 km s⁻¹ (after correcting for the spectral resolution⁸ at 5007 Å). Similarly, the de-convolved Hα line width results in a velocity dispersion σ_Hα = 30.8 ± 3.8 km s⁻¹, consistent within the uncertainties. These dispersions are similar to the ones observed in local star clusters (e.g., Bastian et al. 2006) with extreme stellar densities (Σ_M ∼ 10⁴ − 10⁵ M_⊙ pc⁻²) comparable to T1 (Σ_T1 ∼ 10⁴). On the other hand, the T1 system is very young (≲2 Myr) and shows intense star formation activity, and is therefore possibly a non-virialized system; the observed velocity dispersion could be explained by (radiation) feedback from star formation, as discussed in He et al. (2019, 2020).

The [O III]λ5007 line profile of D1 is unresolved; the velocity dispersion derived from the Hα line is σ_Hα = 27.9 ± 2.5 km s⁻¹, similar to what found for T1. However, D1 is a (morphologically) much more complex system, with a compact bright region which dominates the rest-UV brightness, but not the line emission (see also the discussion in Sect. 4.1).

Lastly, we observe a small shift in the peak wavelength of the brightest lines between D1 and T1 (bottom panels of Fig. 2), which converts to a velocity difference of 42 ± 3 km s⁻¹ (averaging the shift derived from Hα and from [O III]λ5007). The same velocity difference can be inferred from the velocity map in Fig. 5; the velocity within the D1 region is quite uniform, with a small gradient toward the north (i.e., the D1a subregion), while the velocity seems to change abruptly by ∼40 km s⁻¹ between the D1 and T1 systems. Assuming that D1 and T1 are part of a rotating disk, this velocity difference is lower than what observed in (gas) rotation curves of galaxies at z ≳ 4 (e.g., Jones et al. 2017; Rizzo et al. 2020, 2021; Lelli et al. 2021; Fujimoto et al. 2024), although the latter have much larger stellar masses (M_⋆ ≳ 10⁹ M_⊙) than what measured for D1-T1. Alternatively, D1 and T1 may be satellite (gravitationally interacting) systems, ≲1 kpc apart (see Sect. 3.1) and with a relative line-of-sight velocity ∼40 km s⁻¹.

Fig. 5. Velocity map of the system with ellipses marking the positions of the main regions (these are the same apertures as those used for the photometric analysis in Sect. 3.2 and shown over the NIRCam data in Fig. 2). The uncertainties on the velocity measurements (based on bootstrapping from the variance map of the IFU cube) are on average ∼1 km s−1, and ≤4 km s−1 in all pixels. The map is created combining the signal of all the four main lines observed (Angora et al., in prep., for details). The FWHM of the IFU PSF is given in the bottom-right corner.

3.5. Faint diffuse emitting regions

The IFU observations reveal a faint “bridge” of line emission between D1 and T1 and two regions of extremely faint emission on the north side of D1 (see Fig. 6, central panel). The bridge region between D1 and T1 has moderate line emission, with Hα and [O III]λ5007 lines detected with S/N > 5 (but Hβ and [O III]λ4959 are undetected, with S/N < 3; see Fig. 6, right panel). By assuming a Hα/Hβ ratio of 2.86⁹ we derive a metallicity index R3 = 2.3 ± 1.0, slightly lower than for the main D1 and T1 systems. Using the same method described in Sect. 3.4 we derive a metallicity of $Z=0.{03}_{-0.01}^{+0.02}{\mathrm{Z}}_{\odot }$ in the Large EW assumption ($Z=0.{05}_{-0.04}^{+0.14}{\mathrm{Z}}_{\odot }$ in the Small EW case). This bridge region is located close to the UT1 source (white crosses in Fig. 6), but the line emission map does not show a peak in correspondence with UT1. Indeed, NIRCam observations show rest-UV-optical emission coming from a region broader than UT1a and UT1b alone.

Fig. 6. Imaging and spectroscopy of the faint regions. Left panel: NIRCam sum of all filters (both from the SW and from the LW channels). The apertures used to derive the magnitudes (or upper limits) of the faint regions are shown as elliptical contours. The FoV of the IFU is shown as a white contour. Central panel: IFU observations showing the sum of Hα and [OIII] emission lines with the masks used to extract the spectra shown in the right panel. Two white “X” symbols mark the position of the UT1a and UT1b peaks. Right panel: Spectra from the two masks of the central panel, using their same color-coding. As done in Fig. 2, thick shaded lines are used to show the uncertainties of each spectrum.

The faint region close to D1d (which we name D1d_tail) is barely detected in NIRCam, with mag_SW = 29.3 ± 0.2 and mag_LW = 29.0 ± 0.1, corresponding to de-lensed magnitudes mag_SW = 32.3 and mag_LW = 32.0 (Fig. 6, left). Its spectrum reveals Hα and [O III]λ5007 lines detected at S/N ∼ 3; in this case, we derive a ratio R3 = 1.6 ± 0.8 indicating a low metallicity $Z=0.{02}_{-0.01}^{+0.02}{\mathrm{Z}}_{\odot }$ ($Z=0.{03}_{-0.02}^{+0.05}{\mathrm{Z}}_{\odot }$ in the Small EW assumption). The combination of Hα and rest-UV fluxes (following Eq. (1)) indicate a high ionization efficiency, log(ξ_ion/erg⁻¹Hz) = 25.5 ± 0.2.

Finally, NIRCam observations show a faint region southern of T1 (T1_tail); we note that it follows the shear direction of the lens model, extending toward another bright compact source at the same redshift of the D1-T1 system (e.g., V17). The unavailability of IFU coverage in that region prevents us from knowing if also this faint tail is line-emitting. We will discuss this region in the broader context of all the lensed z = 6.145 sources of the Cosmic Archipelago in a forthcoming publication (Messa et al., in preparation).

4. Zooming into the “micro” regions

We leverage the exquisite spatial resolution of the JWST, combined with gravitational lensing, to dig into the substructures of the three main regions just analyzed. Multiple peaks in rest-UV luminosity are particularly evident in the NIRCam SW filters (top-left panel in Fig. 2 and of Fig. 3); when isolating the filters containing line emission (e.g., F356W), also the RGB colors reveal substructures with different properties (e.g., the left panel of Fig. 7 reveals the regions with the most intense line emission in green). A similar conclusion can be derived by inspecting (photometrically-derived) maps of the line equivalent widths and of the rest-UV slope (central and right panel of Fig. 7), where we observe that the bright core of D1 is surrounded by structures with larger EWs, but shallower β slopes. Overall, both the line equivalent widths and the β_UV show variations on small scales: in the following subsections, we characterize these substructures in each of the three main regions.

Fig. 7. Mapping of the main spectral properties. Left panel: RGB composite of the LW filters (red: F410M, green: F356W, blue: F277W); the green channel, containing the emission in Hβ+[OIII] highlights the regions with the strongest line emission. The position of emission peaks in NIRCam SW (see also Figs. 3, 8, 9 and 10) are marked by red crosses (by “star” markers in the central and right panels); the bright core regions D1core and T1core discussed in Sects. 4.1 and 4.2 are marked by red asterisks (by plus plus symbols in the central and right panels). The (white) boxes outline the zoom-in regions shown in the other panels of the figure. Central panels: EW maps of [O III]+Hβ, as traced by the F356W filter (with F410M tracing the continuum); the data have been smoothed by a 2 px box kernel and low signal-to-noise pixel have been masked-out. The map has the same pixel-scale as the LW data (0.04 arcsec/px). Right panels:β-slope maps obtained by the pixel-by-pixel fitting of F115W, F150W and F200W. The data are in a 0.02 arcsec/px scale and have been smoothed by a 4 px box kernel. Low-signal pixel have been masked-out. Both in the central and left panels, a black line of 0.24″ is shown, for scale. The FWHM of the PSF in each panel is reported as a black/white circle in the bottom-right angle. Fig. 8. Zoom-in into the D1 region, showing the combined SW observations (F115W+F150W+F200W), the best-fit model of the core region, and its residuals. The apertures used for photometry are shown as white ellipses; for the D1core region, photometry is derived from a light-profile fitting (as described in the main text). For each of the apertures/regions, photometry is shown as black circles, with the SED best-fit model over-plotted as colored empty squares and relative spectrum. The best-fit of the UV slope is also shown and reported in the panels. Fig. 9. Zoom into the T1 system, shown as the sum of observations in the SW filters, the best-fit model of the core emission and its residuals. The results of photometry and of the broad-band SED fitting are shown in the right panel, both for T1core (blue lines) and for T1diff (orange lines), similarly to what done for the subregions of D1 in Fig. 8. Fig. 10. Same as Fig. 9 but for the UT1 system, with the best-fit model of the UT1a source (left panel) and the photometry + broad-band SED fitting results (right panel, green line for UT1a, red line for UT1b).

4.1. The subcomponents of D1

A compact bright clump that dominates the overall rest-UV emission of the D1 region (D1_core) was already investigated in V19. We fit its light profile in F115W following the methodology described in Messa et al. (2019, 2022). In brief, the source is fitted on the image plane and is modeled with a convolution of a 2D Gaussian function and the filter PSF, on top of a local background described by a one-degree polynomial. The source is barely resolved along the direction of the tangential shear of the lens model (while it is unresolved in the radial direction), with R_eff, tan = 17 ± 1 mas (Fig. 8). Assuming that the clump morphology is intrinsically circular, we use its tangential size and magnification (μ_tan = 13.1; see Appendix C) to derive its size, namely $R_{\mathrm{eff}}\equiv R_{\mathrm{eff},\tan }·{\mu }_{\mathit{tan}}^{-1}$, finding $R_{\mathrm{eff}}=7.6_{-0.3}^{+0.4}$ pc¹⁰. Giving its compact size, D1_core can be considered a single stellar cluster, as already pointed out by the work of V19, who gave an upper-limit < 13 pc using HST observations. By keeping the observed shape of the cluster fixed, and fitting its flux in all filters (as already done for clump studies in Messa et al. 2022, 2024; Claeyssens et al. 2023), we derive the photometry shown in Fig. 8.

The broadband SED is best fitted by a model with a short SFH (τ = 1 Myr), and an age of 12 Myr; the derived intrinsic mass, 4.8 ⋅ 10⁶ M_⊙, implies a surface density¹¹, ΣM = 7.7 · 10³ M_⊙ pc⁻², higher than the average value of stellar clusters in nearby galaxies (Brown & Gnedin 2021) and similar to other high-redshift compact star-forming regions (e.g., Messa et al. 2024). Combining the size and mass of D1_core we also derive (following Gieles & Portegies Zwart 2011) a crossing time of ∼1 Myr; this is shorter than the derived age, suggesting that the source could be a gravitationally bound system (see the discussion in e.g., Gieles & Portegies Zwart 2011). One of the most striking features of the SED of D1_core is the very steep rest-UV slope (β = −2.8 ± 0.1) combined with the absence of nebular emission lines¹² (Fig 8); the best-fit model is unable to reproduce these features, in particular the observed UV slope; we discuss possible reasons for this discrepancy in Sect. 5.1. We tested that even adding a second stellar population (allowing the simultaneous presence of a young and an old population), does not improve the fit. On the other hand, different stellar models could reproduce the absence of emission lines, and the ∼0.5 mag “break” at λ ∼ 4.4 μm (Fig. 8) with an older (∼60 Myr) stellar population, as tested in Appendix E.

A fainter emission is observed around D1_core; the stacking of the SW filters reveals four peaks (which we label with letters, a to d). Despite the flux in each peaks is detected with high S/N in all filters (Fig. 8), the morphological complexity of this region¹³ prevent a light-profile fitting analysis like the one performed for the core. We perform aperture photometry on the four apertures shown in Fig. 8. These apertures have radii ≥0.8″, which is larger than the PSFs (also in the reddest bands) with this setup. As shown in Sect. 3.1, these peaks are found at intrinsic relative distances between 0.2 and 1.6 kpc (Fig. 3). The results of photometry and relative SED fitting are shown in Fig. 8 and collected in Appendix C. We extract the spectra of these subregions from the IFU, using the masks shown in Appendix G: the main properties derived from the spectra are listed in Appendix D. Due to the proximity of D1b and D1c regions, and the faint line emission of the latter, we choose to use a mask that includes both (this merged region is named D1bc in Appendix D and in Appendix G).

All the regions considered are young, consistent with the presence of bright emission lines, as shown by the EW map derived from the NIRCam filters (Fig. 7); two of them have best-fit (mass-weighted) ages between 2 and 5 Myr, the others between 8 and 13 Myr. The (intrinsic) masses are all very similar, in the range 1−2 · 10⁶ M_⊙. D1a region has large equivalent widths (780 Å and 1830 Å for Hα and [OIII]+Hβ), but a shallow rest-UV slope, which could be attributed to a moderate level of extinction (best-fit A_V = 0.6 mag). However, such an extinction is inconsistent with the Balmer decrement we measure from the IFU spectrum extracted at the position of D1a, indicating no extinction. We note how the slope is strongly driven by a dearth of signal in the F115W filter only (Fig. 8), especially in few pixels inside the aperture (as can be also noted in the β map of Fig. 7); excluding that filter from the fit would return a steeper slope, β_UV = −2.26 instead of −1.75. In the case of D1c and D1d the best-fit models do not fully reproduce the observed SED and in particular their steep UV slopes; D1c, in particular, similarly to D1_core, is characterized by faint emission lines but a ${\beta }_{\mathrm{UV}}=-2.7_{-0.1}^{+0.2}$; these discrepancies are discussed in more details in Sect. 5.1. Finally, the metallicities of D1 subregions (0.07 − 0.014 Z_⊙), as inferred from the R3 line ratio, are in agreement with the value found for the entire region (see Appendix D), suggesting that metallicity is overall uniform within D1.

4.2. The core of T1

The structure of the T1 system can be split into a bright core (T1_core) and a fainter stretched “diffuse” emission (T1_diff, Fig. 9). We fit the F115W emission assuming two Gaussian sources; the compact core is barely resolved along the stretch direction (R_eff = 11 mas in the source plane), implying an intrinsic size $R_{\mathrm{eff}}=3.9_{-0.0}^{+0.2}$ pc¹⁴. Within this model, the small size of T1_core qualifies it as a stellar cluster, hosted in the more diffuse T1_diff region. As done for the core of D1, we keep the shape of the sources fixed, and fit their flux in all filters, recovering the photometry shown in Fig. 9. Due to the coarser resolution, the separation between the core and the diffuse contribution to the flux is hard to establish in the LW filters, leading to large photometric uncertainties. For the same region, we prefer to infer the photometry of the diffuse component from aperture photometry on the core-subtracted data (see Fig. 9). Both T1_core and T1_diff are characterized by steep UV slopes and bright emission lines; however, while the core has a steeper slope than the diffuse part, it also has smaller rest-frame equivalent widths; this spatial trend can be appreciated from the maps in Fig. 7, where we observe that the strongest emitting region is slightly displaced from the location of the core. This is confirmed by the EW derived from the IFU data by using a mask centered on the core (shown in Appendix G), showing that the core has a flux ∼50% lower than that inferred for the entire T1 mask (see Appendix D). Similarly, the ionization efficiency derived from the Hα flux in the core mask is relatively lower, log(ξ_ion/erg⁻¹Hz) = 25.1. On the other hand, the Hα/Hβ ratio and the R3 index remain similar to the values obtained for the entire T1 system.

The SED fitting analysis returns a very young age (≤2 Myr) both for T1_core and for T1_diff. The derived (intrinsic) mass of the core is $2.1_{-0.2}^{+8.8}·{10}^5{\mathrm{M}}_{\odot }$, implying a stellar mass surface density ΣM = 1.5 · 10³ M_⊙ pc⁻², consistent with the typical density of local star clusters (e.g., Brown & Gnedin 2021). As was the case for the T1 system overall, the best-fit model of the core cannot reproduce the β slope. The diffuse part of T1 retains most of the mass, $M_{\star }=9.8_{-1.8}^{+0.5}·{10}^5{\mathrm{M}}_{\odot }$.

4.3. A bound stellar cluster in UT1

Observations in the SW filters split UT1 into a diffuse NW component (UT1b) and a brighter and more compact SE source (UT1a); the latter is an unresolved source with R_eff < 9.5 mas (before de-lensing) and is therefore consistent with a star cluster of size R_eff < 3.8 pc¹⁵ (see Fig. 10). Aperture photometry and best-fit results are presented for both components in Fig. 10 and Appendix C. The broad-band SED fitting returns ages robustly older than 10 Myr and (intrinsic) masses 1.2−4.2 · 10⁶ M_⊙, despite large uncertainties on the best-fit values for both quantities (see Appendix C). As already pointed out in the previous section, the derived ages of UT1 (and of its subcomponents) indicate that this system experienced a burst of star-formation that is already over, unlike D1 and T1 systems. Combining the mass and size of the UT1a system, we derive a stellar mass density ${\mathrm{\Sigma }}_M>7.7_{-3.4}^{+29.7}·{10}^3{\mathrm{M}}_{\odot }$ pc⁻² and a crossing time of 1 Myr: the latter is considerably shorter than the age of the clump, suggesting that UT1a is a gravitationally bound dense stellar cluster. Despite its steep UV slope (${\beta }_{\mathrm{UV}}=-2.{68}_{-0.25}^{+0.18}$) combined with a lack of line emission, similar to what found in D1_core, the best-fit model reproduces the broadband SED, partly due to the larger photometric uncertainties.

5. Discussion

5.1. The cores of D1 and T1: Steep UV slopes and possible escape of ionizing radiation

The large equivalent widths found for T1 undoubtedly classify the system as very young. The analysis of line ratios (Sect. 3.4) also allows us to infer that the system is metal-poor, has very little or no extinction and is characterized by an efficient production of ionizing radiation A very steep rest-UV slope is also observed, which is typical of young and low-metallicity systems, although the latter is not fully reproduced by our best-fit model (Fig. 9), which tends to return a shallower slope. This is even more pronounced when focusing on T1_core, where the slope is more extreme (β_UV = −3.1). We point out that this discrepancy is not driven by the stellar models used, as it is observed also in the best-fit model obtained with nonbinary models (Appendix E). More in general, BPASS models have on average bluer UV slopes (by ∼0.1) than models not including binaries; however, all population synthesis models including nebular emission struggle to reproduce slopes steeper than β < −2.8 (e.g., Bouwens et al. 2010; Topping et al. 2022; Bolamperti et al. 2023). We also point out that our measures of the UV slope are not due to the adopted photometric approach and apertures, as suggested by the pixel-by-pixel map of Fig. 7. Galaxies with extremely blue slopes (β < −2.8), although rare, are observed especially at high redshifts, z ≳ 5; reproducing their SEDs requires models of density-bounded HII regions with non-zero ionizing photon escape (e.g., Topping et al. 2024). If this was the case for T1, the escape fraction cannot be 1, as that solution would be in contrast with the large nebular line equivalent widths observed. However, models with positive escape fractions display elevated OIII+Hβ EWs (> 2000 Å) even in case of f_esc ∼ 0.8 (Topping et al. 2022, 2024). Interestingly, we have already noted in Sect. 4.2 how T1_core is a region with bluer slope but weaker line emission (in terms of EWs and ξ_ion; see also Fig. 7 and Appendix D), compared to its surroundings, hinting at the possibility that the core is a region of ionizing radiation leakage. Starting from the properties of low-z leaking galaxies (Flury et al. 2022a,b), indirect tracers of the escape of ionizing photons have been proposed in the literature (e.g., Mascia et al. 2024, 2023); such indicators suggest an increase of the escape for compact galaxies with high Hβ equivalent widths (EW ≳ 200 Å) and steep UV slopes, consistently with the properties of the T1 region. In particular, slopes ≤ − 2.5 seem to be robust tracers of f_esc ≳ 0.1 (Chisholm et al. 2022). We remind however that these tracers have been calibrated and studied for entire galaxies at larger scales (R_eff ≳ 100 pc) and lower redshifts than our systems.

Another region where we note a discrepancy between the best-fit model and observations is D1_core, characterized by a steep UV slope (β_UV = −2.8); the absence of nebular lines “forces” the fit of the SED to ages of ∼10 Myr (or older; see Sect. 4.1 and Appendix E), though slightly under-estimating the slope (β_model = −2.6). Another possibility is that the source is younger (age < 10 Myr) and leaking (f_esc > 0). Recently, other high-z galaxies with strong UV but no line emission has been discovered by the JWST (e.g., Looser et al. 2024; Topping et al. 2024). Some possible scenarios used to describe this kind of SED are a bursty star formation (in which the source is observed just after the stop of its starburst phase, i.e., it is in an ’off’ mode), and/or the presence of spatially varying dust obscuration (e.g., Faisst & Morishita 2024); in our case, both from photometry and from spectroscopy we did not find any hint of moderate dust obscuration. The spatial variation in UV slope within the D1 region (Fig. 7, right panel) may suggest a varying dust effect. However, the pixels with shallower slopes (β > −2) are located in regions of low signal; as a consequence the slopes derived there are less robust and could also be naturally driven by the presence of older stellar components than in the bright clumps. Concerning the “bursty” scenario, this is contemplated by the best-fit model, described by a short burst with τ ∼ 1 Myr and an age ∼10 Myr; however, the discrepancy observed between the bets-fit result and the photometry in the rest-UV filters seem to disfavor it.

Other subregions analyzed show features similar to the cores; for example, the SED of the young subregion D1d is characterized by large EWs and a slope (β = −2.6) which is not fully reproduced by the best-fit model. On the other hand, D1c has an even steeper slope (β = −2.7) and fainter line emission, leading to a best model with age ∼10 Myr, similarly to D1_core; like in the latter, a younger age with f_esc > 0 would explain the UV slopes measured.

The hypothesis of leakage seems supported also by the Lyα emission of the system (presented in V19); when comparing it to the [O III]λ5007 line, taken as reference for the systemic redshift, we note that Lyα peaks at a redshift corresponding to Δv ∼ +100 km s⁻¹, both in D1 and T1 (Fig. 11). We assume that we are observing the “red” peak of the line, while the “blue” one is not observed due to intervening IGM absorption, which at z ∼ 6 can attenuate the line up to ∼80% (e.g., Laursen et al. 2011; Tang et al. 2024). In this scenario the observed peak would be further redshifted (Laursen et al. 2011); we conclude that a velocity shift of only ∼100 km s⁻¹ despite the IGM absorption and internal radiation transfer processes suggest low opacity for the system, facilitating the escape of ionizing radiation. For this reason, the Lyα peak separation has proved to be a strong indirect tracer of LyC escape (e.g., Verhamme et al. 2017; Flury et al. 2022a,b), also in case of low-mass (M_⋆ < 10⁸ M_⊙) galaxies (Izotov et al. 2021b).

Fig. 11. Comparison between the [[O III]λ5007] line emission of the D1 (solid blue) and T1 (solid red) regions (extracted from the spectra presented in Fig. 2) and the Lyα flux observed at the same position (from VLT/MUSE; see Vanzella et al. 2019, 2021), rescaled by a factor 0.5 in order to ease the comparison of the peak relative velocities. For both regions, the center of the [O III]λ5007 line has been used as the systemic redshift and therefore as reference for the x-axis. The Lyα line peaks at Δv ∼ +100 km s−1 in both regions.

5.2. Small-scales analysis in comparison to high-z emitters in the literature

In the absence of lensing magnification, the entire D1-T1-UT1 system would have been barely resolved, with radius ∼1 kpc (Fig. 3). By merging together the photometric and spectroscopic analyses of the main subregions (Sect. 3), we recover the following features: (i) an intrinsic UV magnitude of −17.8 mag; (ii) large equivalent widths, ∼650 Å in Hα and ∼900 Å in Hβ+OIII and (iii) a steep UV slope, β = −2.5, both indicating the presence of a young stellar population; line ratios indicating (iv) no extinction, and (v) low metallicity, Z ∼ 10% Z_⊙ (from [O III]λ5007/Hβ = 3.9); (vi) an intrinsic SED-derived stellar mass in range M ∼ (1−3) · 10⁷ M_⊙. These properties would place the galaxy as a typical (Ly-α emitting) system among the UV-faint objects (M_UV > −18) at its redshift (e.g., Bouwens et al. 2014; Shibuya et al. 2015; Bhatawdekar & Conselice 2021; Topping et al. 2022, 2024; Weibel et al. 2024), with large equivalent widths that would classify it among the EELGs (e.g., Tang et al. 2023). The large magnification of the Cosmic Archipelago system gives therefore the unique opportunity of a detailed view into the sub-galactic scales of a z = 6 Lyα emitter. First of all, the intense line emission can be resolved down to few young sources with EW(Hβ+O III) > 1000 Å and a hard ionizing field (log(ξ_ion/erg⁻¹Hz) ≥ 25.5), on scales approaching the ones of individual young massive clusters, down to a few pc (e.g., T1, D1a). Our analysis also indicates a reasonable possibility of ionizing radiation leaking at the same scales of individual clusters (e.g., T1_core, D1_core, D1c). No clear hint of this possibility would come from an un-resolved study of the region, in case of indirect tracers like EW(Hβ) or the UV slope (e.g., Mascia et al. 2024, 2023); the information given by the lack of emission lines combined to the steep UV slope in the core of D1 is easily washed out when considering the system on a larger scale (Sect. 3). The famous galaxy dubbed Sunburst arc at z = 2.4 shows that the escape of LyC emission comes from a very compact region, consistent to a massive cluster, with R_eff < 10 pc (Rivera-Thorsen et al. 2017, 2019; Vanzella et al. 2022), and co-spatial to the observation of a steep UV slope (β ≈ −3, Kim et al. 2023); similar results come from the spatial characterization of f_esc in the very nearby galaxy NGC4214, where the escape of ionizing radiation (as high as f_esc ∼ 40%) is observed in small sub-galactic regions but would have been undetectable from an integrated study (Choi et al. 2020). These studies demonstrate how an insight on the small scales of emitting galaxies could be fundamental to gain a deeper view on the escape of ionizing radiation from galaxies, and, in turn, on the study of reionization.

5.3. Possible biases of integrated studies

Both the sub-galactic scale variation of EW, β_UV, ξ_ion and the possible identification of leakers suggest that galactic-scale SED analyses could lead to biased properties (in terms of e.g., masses and SFRs) of the galaxy stellar population within the galaxy. A similar concern was raised recently by Giménez-Arteaga et al. (2023, 2024) and Bradač et al. (2024) via the sub-galactic scale analysis of gravitationally lensed galaxies at redshifts between 5 and 9; the authors demonstrated how the presence of an old, massive population is commonly missed by SED analysis of high-z UV-bright galaxies. In our analysis, spatially dissecting the D1-T1-UT1 systems into small scales, no obvious sign of older populations was found; all the subregions analyzed show overall similar colors and properties, in particular ages well below 100 Myr. We however recognize the possibility that an old population can be hidden behind the brightest UV regions; we will test this possibility in detail in a forthcoming publication (Messa et al. in prep.).

On the other hand, we observed, especially for D1, how the SED analysis of the entire region returns a moderate age (∼40 Myr), which is not found when breaking down the SED into smaller apertures (≲10 Myr). We attribute this discrepancy to the fact that the photometry of the entire region mixes up the small-scale variation easily observed from the maps of Fig. 7. In addition, the limited number of filters available for SED fitting required the assumption of a “simple” SFH model, which may not represent the complexity of the region, and also prevents breaking down the degeneracies associated with the parameters (e.g., the SFR decline τ and the age of the system), as also tested in Appendix E.

5.4. Low-metallicity regions

The emission line analysis (Sect. 3.4) indicates that the system has a sub solar metallicity, Z ≲ 0.15 Z_⊙; the R3 diagnostics used for this measure are calibrated via the study of local galaxies, and are associated with large uncertainties, especially in the metal-poor regime (≲ 0.1 Z_⊙, Nakajima et al. 2022). We improve the accuracy of the metallicity indicator by using the equivalent width of Hβ, tightly tracing the ionization states; the large EWs of T1 and D1a allow to robustly infer metallicities Z ∼ 0.05 Z_⊙ in those regions. Such low values are in line with the deficiency of [CII] emission and the non detection of dust continuum in the D1-T1 system (Calura et al. 2021).

Interestingly, two faint subregions (D1d_tail and the bridge region between D1 and T1) show even smaller [OIII]/Hβ ratios, indicating a more metal-poor environment (Z ∼ 0.02−0.05 Z_⊙; Sect. 3.5); these are among the most metal poor regions observed at any redshifts (e.g., Nishigaki et al. 2023; Vanzella et al. 2023b; Izotov et al. 2024; Curti et al. 2024). Those regions are intrinsically separated by hundreds of parsecs from the main UV-bright regions (Fig. 3) and, if not for the lensing, they would have been outshined by them in an unresolved view of the systems. Their faintness implies that the continuum regions associated with the emission has a low stellar mass. In the case of D1d_tail, assuming that it has a young age comparable to what found for D1d, the observed mag_SW = 29.3 ± 0.2 converts to an intrinsic mass M_⋆ ∼ 10⁵ M_⊙; yet, its measured ionizing photon production efficiency is large, log(ξ_ion/erg⁻¹Hz) = 25.5, indicating the presence of massive stars (e.g., Raiter et al. 2010; Stanway & Eldridge 2023; Schaerer et al. 2025). Another similar region, with low mass and metallicity but high ionization rate, is found near the T2 region of the same Cosmic Archipelago system; also in this case the low-Z region is found ∼200 pc away from the main UV source, in the source plane. The source is described in detail in Vanzella et al. (2024) and we refer to the latter for a detailed description of its implications.

5.5. Bound star clusters at high-z

Finally, the small-scale view of the galaxy reveals individual clusters, including a possibly gravitationally bound system with age of 13 Myr. We point out that the “boundness” mentioned in this work is based on the comparison between crossing time and age of the system, and therefore refers to its “natal” condition; from local studies we know that while ∼90% of the star formation in galaxies take places in clustered environment, only a small fraction is intrinsically bound, while the rest is dissolved in few Myr (e.g., Lada & Lada 2003). The long-term survival of the clusters observed in this work will be affected by the interaction with the host galaxy, via dynamical friction and tidal interactions; only a tiny fraction of clusters is expected to survive for cosmological times (e.g., Katz & Ricotti 2014; Reina-Campos et al. 2022, 2023).

Candidate bound clusters (sometimes referred to as “proto-globular clusters”) have currently been observed in other galaxies at similar or even higher redshifts, thanks to JWST observations, like in the Sunrise arc at z = 6 (Vanzella et al. 2023a), in the Firefly Sparkle galaxy at z = 8.3 (Mowla et al. 2024) and in the Cosmic Gems arc at z ∼ 10 (Adamo et al. 2024). Other compact (R_eff ≲ 10 pc) and bound clusters are found within systems of the Cosmic Archipelago; we leave their analysis and the overall comparison (in terms of mass and density) to other high-z star clusters to a forthcoming publication (Messa et al., in prep.).

6. Conclusions

The Cosmic Archipelago is a system of gravitationally lensed galaxies at z = 6.145 in the galaxy cluster field MACS J0416.1–2403, and was initially discovered as Lyα halos, but also showed compact rest-UV morphologies from HST data (Caminha et al. 2017; Vanzella et al. 2017, 2019). The entire galaxy cluster field has been covered with JWST-NIRCam imaging in eight filters by the Prime Extragalactic Areas for Reionization and Lensing Science program (PEARLS, PID 1176; Windhorst et al. 2023) and the CAnadian NIRISS Unbiased Cluster Survey (CANUCS, PID 1208; Willott et al. 2022). Part of the Cosmic Archipelago has also been observed with four pointings of the JWST NIRSpec-IFU, using the high-resolution G395H/F290LP grism and filter combination (GO 1908, PI: E. Vanzella). In this work, we exploit these new JWST observations to study the systems named D1, T1, and UT1, which are covered by one of the IFUs. The other IFU pointings, and the rest of the Cosmic Archipelago in general, will be analyzed and discussed in forthcoming publications (Vanzella et al. 2024, Messa et al., in prep., Bolamperti et al., in prep.).

The rest-UV and optical morphology of the target system is made of three main regions, D1, T1, and UT1, which are confined within an intrinsic radius of ≲1 kpc. This is slightly (approximately three times) larger than the average size of galaxies with similar UV magnitudes (M_UV ∼ −18 mag) at a redshift of ∼6 (Morishita et al. 2024), suggesting that they may be separate satellite systems, instead of subcomponents of a single galaxy. This conclusion is supported by a velocity difference of ≤40 km s⁻¹ between D1 and T1, which is smaller than the typical velocities observed in the rotation curves of galaxies at similar redshifts.

The large magnification of the system (μ ≥ 17) allows us to distinguish several subregions, namely a bright compact core within D1 (D1_core), one in T1 (T1_core), and other rest-UV peaks (D1a, b, c, and d within D1, and UT1a and UT1b within UT1). We performed photometry both for the main regions and for the subregions, which is fitted via Bagpipes to obtain ages, masses, extinctions, and SFRs. In parallel, we extracted spectra from the IFU, where we measured the flux of Hβ, [O III]λ4959, [O III]λ5007, and Hα lines, and the relative equivalent widths. Line ratios probe the attenuation and the metallicity of the systems.

T1 is the youngest system (1 Myr old), and also shows the largest line EWs (∼2800 Å but reaching values ≳ 3000 Å in the southern part; see the map in Fig. 7), the largest ionizing photon-production efficiency (log(ξ_ion/erg⁻¹Hz) = 25.7), and the steepest UV slope (β_UV = −2.7). These are extreme values, rarely observed at similar redshifts even in EELGs (e.g., Endsley et al. 2023; Boyett et al. 2024; Topping et al. 2024; Nanayakkara et al. 2023). In addition, the SED fitting also suggests a larger ionization parameter than that usually found in star forming and starburst galaxies (logU > −2); the value we find is typical of systems with large SFR densities, as is the case for T1. The latter is a small system with R_eff = 19 ± 1 pc, M_⋆ = (1.0 ± 0.1) · 10⁶ M_⊙, and $\mathrm{SFR}=(0.8_{-0.2}^{+0.0}){\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1}$ (implying Σ_M⋆ ∼ 300 M_⊙ pc⁻² and Σ_SFR ∼ 200 M_⊙ yr⁻¹ kpc⁻²), which are consistent with a compact HII region; its extreme ionization properties suggest the presence of massive stars.

D1 has an older mass-weighted age of 42 Myr, yet the best fit returns a prolonged SFH that would explain the nebular emission lines observed. Interestingly, the individual fits of its compact subregions (selected in the rest-UV) return ages of ≤15 Myr in all cases, and very young ages (≤5 Myr) for some regions (namely D1a and D1b). In the same way, while the overall equivalent width of the system is EW(Hβ,[O III]) ∼ 600 Å, the spatial map of Fig. 7 reveals regions with EW > 1000 Å, at the same positions as the rest-UV peaks (e.g., EW(Hβ,[OIII]) ≥ 1500 Å in D1a). This result, in addition to what is found for T1, suggests that the large EWs characterizing most star-forming galaxies at z ≳ 6 (e.g., Boyett et al. 2024) can be powered by individual compact clumps (star clusters and/or HII regions) at small sub-galactic scales. Finally, UT1 has a similar mass-weighted age to D1, but a much shorter star-formation duration, which is consistent with the absence of emission lines.

For the compact regions D1_core, T1_core, and UT1a, we are able to measure an effective radius, which is in all cases R_eff < 10 pc. These systems have masses in the range of 0.2 − 4.8 ⋅ 10⁶ M_⊙, leading to Σ_M = (1−8) · 10³ M_⊙ pc⁻² and are therefore consistent with being individual star clusters. At 14 Myr old, UT1a is much older than its crossing time (∼1 Myr), implying that the cluster is gravitationally bound, and thus a proto-globular cluster. The same is true for D1_core; however the best-fit model is unable to reproduce the SED shape, in particular the steep UV slope (β_UV = −2.82 ± 0.05) combined with the lack of nebular lines. Few such SEDs have been observed in high-z galaxies (Topping et al. 2024) and are consistent with young systems leaking ionizing radiation, resulting in both a steepening of β_UV and a suppression of the nebular line emission. T1_core, on the other hand, shows strong line emission, undoubtedly characterising the system as young. However, due to an extreme β_UV = −3.1 ± 0.2 (rarely observed at any z), this is another case where the best-fit model cannot reproduce the observed SED, suggesting that this cluster may also be a LyC leaker; this scenario is supported by the fact that the EWs of T1_core are lower than those of the entire T1 region. The hypothesis that the cores of D1 and T1 are leakers is consistent with the Lyα emission extracted at their position (see Vanzella et al. 2019, 2021), which peaks at Δv ∼ +100 km s⁻¹ from the systemic redshift, therefore implying low opacity in those systems. We conclude that not only is most of the ionizing radiation produced at the scales of star clusters, but it can also escape at that same approximately parsec scale: being able to resolve small sub-galactic scales may therefore be fundamental in understanding the process of reionization and the nature of the sources driving it (e.g., Ricotti 2002).

Both D1 and T1 are low-metallicity systems, with $Z=0.{14}_{-0.06}^{+0.11}{\mathrm{Z}}_{\odot }$ and $Z=0.{05}_{-0.02}^{+0.02}{\mathrm{Z}}_{\odot }$, respectively, inferred from the indirect tracer [O III]λ5007/Hβ; similar values are also found in all their subregions. However, we note the presence of some regions with faint line emission that are characterized by lower R3 values, leading to smaller inferred metallicities. This is the case, in particular, for D1d_tail, a region separated by hundreds of parsecs from D1d and, more in general, from the main D1 region (Fig. 3), and characterized by Z ∼ 0.02 Z_⊙. The same region, while very faint in the rest-UV, has a measured log(ξ_ion/erg⁻¹Hz) = 25.5, indicating a system with a stellar mass of ≲ 10⁵ M_⊙ and the presence of massive stars. D1d_tail may therefore be a region where a very recent episode of star formation has not yet enriched its metal-poor gas. Another, similar region has been observed in the Cosmic Archipelago system, which is presented in a separate publication (Vanzella et al. 2024).

In conclusion, the sub-galactic view of this z ∼ 6 Lyα emitter, obtained using a combination of strong gravitational lensing and state-of-the-art JWST observations, revealed small and faint details that would be missed in studies of integrated galaxies yet are fundamental to understand the observed integrated properties of high-z galaxies (e.g., large EWs, steep UV slopes, large ξ_ion). In addition, our analysis reveals the formation of proto-GCs, whose properties can be directly studied at their formation epoch. While the results presented here are based on only one IFU cube, referring to one of the Cosmic Archipelago ∼1 kpc regions (D1-T1-UT1), three more IFU observations and their relative systems will be presented in a forthcoming publication (Messa et al., in prep.).

Acknowledgments

This research made use of Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2023). This work made use of v2.2.1 of the Binary Population and Spectral Synthesis (BPASS) models as described in Eldridge et al. (2017) and Stanway & Eldridge (2018). The research activities described in this paper have been co-funded by the European Union – NextGeneration EU within PRIN 2022 project n.20229YBSAN – Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch. M.M. acknowledges the financial support through grant PRIN-MIUR 2020SKSTHZ. FL acknowledges support from the INAF 2023 mini-grant “Exploiting the powerful capabilities of JWST/NIRSpec to unveil the distant Universe. AA acknowledges support by the Swedish research council Vetenskapsrådet (2021-05559). R.A.W., S.H.C., and R.A.J. acknowledge support from NASA JWST Interdisciplinary Scientist grants NAG5-12460, NNX14AN10G and 80NSSC18K0200 from GSFC. MB acknowledges support from the ERC Grant FIRSTLIGHT and from the Slovenian national research agency ARIS through grants N1-0238 and P1-0188.

References

Appendix A: Post-processing of the NIRSpec IFU cubes and background subtraction

The methodology used for NIRSpec data reduction are presented in detail in Sect. 2.2. In this appendix we present the direct comparison of line fluxes and ratios, measured in the main target regions, for the cubes obtained at different stages of the process; we also show some relevant intermediate products of the data reduction.

We discussed in the main text how, after stage 2 of the JWST-NIRSpec official pipeline, we perform a post-processing in order to remove spikes and artifacts still on the cube. An example of such features is shown in Fig. A.1 (left panel); the slice shown is at λ = 4.69 μm, close to the wavelength of the Hα emission at the redshift of the source (and a strong Hα emission can be observed for the T1 system; see e.g., Fig 2). For three sets of pixels close to the target, flagged as nan (indicated by white arrows in Fig. A.1), the pipeline had problems deriving robust values. In addition, spurious signal is detected in the same area (indicated by the blue arrows). The post-processing treatment applied to the cube correct for both Fig. A.1 (central panel), leaving a “clean” cube. The differences between the pipeline-output cube and the post-processed one, in line fluxes and ratios, measured in the main regions are on the order of ≲5% (Table A.1). The final background estimate is based on a moving median across 30 slices in wavelength. This methodology takes care of the spatial variations on the plane of the sky (Fig. A.1, right panel) which a scalar median value cannot account for; it also accounts for “oscillations” in the background, on ranges of ∼ 0.05 μxsm. These are secondary effects (Fig. A.2) that affect the line fluxes and ratios only up to ≲8% (Table A.1).

Finally, the comparison of line fluxes before and after the PSF correction of the cube is also shown in Table A.1. By construction, the PSF correction leaves the cube at the Hα wavelength unaltered, and applies a smoothing at lower wavelengths; for this reason the largest difference is observed for the Hα/Hβ ratio (∼7%), while for lines close in wavelength space the difference is minimal (e.g., ∼1% for the R3 ratio).

Fig. A.1. Left panel: A slice (corresponding to λ = 4.69 μm) of the final cube produced via the standard JWST-NIRSpec pipeline, showing Hα emission from the T1 region (see also Fig. 2) but also spurious signal blending with the real one (indicated by the blue arrows) and nan pixels (indicated by the white arrows). The same slice, for the cube post-processed by our custom pipeline, is shown in the central panel. Right panel: the residual background, after subtracting a scalar median one (i.e., the dark green background in Fig. A.2), is shown for the same slice; background spatial variations on the order of ∼10−21 erg s−1 cm−2 Å−1 are observed.

Fig. A.2. Left panel: Spectra extracted in the D1 mask (see Fig. 2) for the cube without background subtraction (light red line) and after the background subtraction (gray line), where the solid black line is the total measured background (based on a moving median across 30 slices) within the mask. The solid red line is the background obtained using one scalar value per slice. Right panel: the difference between the total background and the scalar one is shown as a dark green thick line, in a zoom-in at the wavelengths of Hβ, [O III]λ4959 and [O III]λ5007 emission lines. The residual signal produced using only the scalar background is also shown as a thin green histogram, for reference.

Appendix B: NIRCam maps with all main and subregions

Fig. B.1. Left panel: the three main regions (D1, T1, and UT1; Sect. 3), along with the two “faint” regions discussed in Sect. 3.5, on top of the sum of all SW and LW filters’ observations. Right panel: All the micro regions discussed in Sect. 4, on top of the sum of the SW filters’ observations. The ellipses marking the extent of the main regions are kept, in transparency. In both panels the FoV of the IFU is shown as a black contour. All the IFU masks used in this work are shown in Appendix G.

Appendix C: Table of photometric results

Appendix D: Table of spectroscopic results

Appendix E: Testing stellar models for the SED fit

We report, in Fig. E.1 the comparison of the best-fit models between BPASS, the reference stellar population synthesis model considered in the current work, and the Bruzual & Charlot (2003) models (B&C03), the default one of Bagpipes (Carnall et al. 2018). The comparison is given for the cases where, as discussed in the main text, the best-fit model was unable to reproduce the steep UV slopes observed (T1, T1_core, D1_core and D1c, all discussed in Sect. 5.1). The best fit of T1(total) is almost identical in the case of the two different models, both in terms of predicted magnitudes and of derived properties. The results are similar also in the case of T1_core; the B&C03 fit returns an older result (with an age of 7 Myr) that perform slightly worse in the LW filters (as also reflected by the larger ${\chi }_{\mathrm{red}}^2$). Much larger differences are observed when comparing models for D1_core and D1c; in both cases the B&C03 fit prefers prolonged SFRs (with τ ≈ 200 Myr for D1c) resulting in older ages and larger masses by up to a factor 10. In the case of D1_core the B&C03 fit is closer to the observed photometry in the SW filters but perform worse in the LW ones, while more comparable magnitudes are found in the case of D1c. Despite the large differences in the derived physical properties, the predicted magnitudes are extremely close to the observed ones in both cases, suggesting some degeneracy between the exponentially-declining SFH model (and in particular the τ parameter defining the decline of the SFR) and the age of the systems. We note that, in all cases, both models systematically underestimate the steep UV slopes observed, even if in some cases they lie within the observational uncertainties.

Fig. E.1. Photometry of the four regions with the steepest UV slopes observed (black markers), along with the best fit model using BPASS (blue diamonds and line) and Bruzual & Charlot (2003) (red squares and line) stellar synthesis population models. On the bottom of each panel the best-fit derived physical properties are given (mass-weighted age, extinction, stellar mass and reduced χ2).

Appendix F: Faint and undetected lines

We discuss here the upper limit of undetected lines tracing the properties of the ionized gas (e.g., electron density, n_e, and temperature, T_e), starting from the Oxygen line [OIII]λ4363. We find, in the spectrum of the D1 region, at the expected wavelength of the line, a signal with a S/N of 3.5 (Fig. F.1). When fitting this signal with a gaussian model, as done for all the other detected lines, its observed width (σ = 2.0 Å) is smaller than the spectral resolution at its wavelength (σ = 6.3 Å). Moreover, its fitted central wavelength deviates by Δz = 4 · 10⁻⁴ from the wavelength expected from the redshift derived from the two strongest lines (Hα and [O III]λ5007). Both these properties suggest that this signal may be spurious; on the other hand they could simply be explained by the low S/N. We also point out that a Δz = −4 · 10⁻⁴ difference in inferred redshift is measured when comparing Hα to Hγ lines in D1 (the latter with a S/N ∼ 7, i.e., brighter than [OIII]λ4363). In this work we consider this signal as a “tentative” detection. If considered real, the ratio between its flux (3.1 · 10⁻¹⁹ erg s⁻¹ cm⁻²) and the flux of [O III]λ5007, suggests that the gas in D1 is either very hot (T_e ≫ 10⁴ K) or very dense (n_e > 10⁴ cm⁻³), or a combination of both; the absence of detected lines sensitive to the gas density prevent to break this degeneracy. Assuming a temperature T_e ≈ 25.000 K and, as consequence, a density n_e ≈ 6 · 10⁴ cm⁻³, using the “direct method” (e.g., Izotov et al. 2006) we would derive a low-metallicity for the system (12+log(O/H)  ∼ 7.2, i.e., ≲ 5% Z_⊙) consistent with what discussed in Sect. 3.4. The metallicity would be higher in case the gas was cooler, T_e ≈ 12.000 K, but, according to the [OIII] line ratio, this would require extreme densities, n_e ∼ 10⁶ cm⁻³.

Another line commonly used in literature to study the metallicity of galaxies and star forming regions is [NII]λ6585; we do not observe any signal at the expected wavelength of the line in the spectrum of D1 (Fig. F.1, right panel). This non-detection is expected, given the low metallicity of the system. The 5σ upper limit we derive (F_NII < 4.8 · 10⁻¹⁹ erg s⁻¹ cm⁻²) is not informative: its ratio to the Hα flux, used as a metallicity proxy, implies Z < 60% Z_⊙ (following the calibrations of, e.g., Curti et al. 2020; Nakajima et al. 2022). Finally, we point out the neither [OIII]λ4363 nor [NII]λ6585 lines are observed in the spectrum of the T1 regions (where the lines detected are fainter than in D1).

Fig. F.1. Spectra of the D1 region covering the wavelengths nearby the Hγ(left panel) and the Hα(right panel) lines, including the best-fits of the observed lines. The left panel also includes a tentative detection of the [OIII]λ4363 line. An artificial 5σ signal at the expected wavelength of the [O III]λ4363 and of the [NII]λ6585 lines is shown as a blue-dashed line.

Appendix G: Spectra of the main subregions

Fig. G.1. Left panel: map of the Hα+[OIII] emission from NIRSpec-IFU; the masks used to extract the spectra of the subregions discussed in Sect. 4.1 and Sect. 4.2 are over-imposed. The relative spectra (covering the wavelengths around Hβ, [[O III]λλ4959, 5007] and Hα) are shown in the right panels. All the other IFU masks used in this work (and shown in Fig. 2 and 6) are included for comparison as white contours (with black labels).

All Tables

All Figures

Fig. 1. Color image combining HST and JWST observations from the 2023-146 release, which is available at the following link, with over-plotted: (i) the Lyα emission from VLT-MUSE of the Cosmic Archipelago systems, at z = 6.145, as yellow contours; and (ii) the four NIRSpec-IFU pointings of program GO 1908, across the main arc and also covering the compact Lyα halo D2, as green squares. The pointing analyzed in detail in the current publication covers the D1 and T1 systems (darker green square). (iii) The other pointing of GO 1908, denoted with a cyan square and covering the metal-poor system LAP1 at z = 6.63, is presented in a separate publication (Vanzella et al. 2023b).

In the text

Fig. 2. Imaging and spectroscopy of the D1T1 system. Top panels: Cutouts showing the observed systems in the sum of NIRCam SW (F115W, F150W, F200W, left) and LW (F277W, F356W, F410M, F444W, center) filters and a NIRSpec-IFU map with the sum from Hα, Hβ, and [OIII] line emission (right). White ellipses mark the apertures used to extract photometry of the three main regions (labeled in the left panel). The white dashed aperture is used only to compare NIRCam to NIRSpec-derived quantities for the T1 region (see Sect. 3.4). The right panel also contains the pixel masks (blue and red contours) used to extract the spectra shown in the bottom panels. Middle panel: Spectra for D1 and T1 regions (in blue and red, respectively). No stellar or nebular continuum is detected in the spectra. The gray shaded area marks the wavelength range falling in the NIRSpec detectors’ gap. The insets in the bottom panels are zoom-ins around the detected lines and display also the uncertainty of the spectra (as wide shaded lines). The best-fit center of the lines is displayed as vertical lines and reveals a small shift in redshift between the two regions (corresponding to a velocity shift of ∼42 km/s; see also Fig. 5, Sect. 3.4 and Appendix D).

In the text

Fig. 3. Positions of sub-regions in the image and source planes. Left panel: Labels of the main rest-UV peaks and subregions of the D1-T1-UT1 systems on the stack of the SW filter observations. The inset shows a stacking of all filters where the three main regions (D1, T1 and UT1) are delimited by the red contours, while larger contours enclosing also the low surface-brightness emission are shown in green. Right panel: Mapping of the main and subregion positions on the source plane at z = 6.145 obtained from the best-fitting lens model by Bergamini et al. (2023). The compact cores (D1core, T1core and UT1a) are marked by white filled circles, while the other peaks of UV emission are marked by white “X” symbols. A 1 kpc (0.17″) reference scale is given in the top left corner. The FWHM of the F200W is also given as reference. The location of T1core coincides with the peak of the Lyα emission introduced in Sect. 1. The bridge region refers to the region of faint line emission discussed in Sect. 3.5.

In the text

Fig. 4. Photometry of the main regions, D1, T1, and UT1 (black circles), obtained via aperture photometry (apertures shown in Fig. 2). The best-fit UV beta slopes are shown as blue lines (and relative blue-shaded uncertainties) and reported at the top of the panels. The results of the broadband-SED fitting are over-plotted (red empty squares and red solid line: spectrum of the maximum-likelihood model, with red-shaded uncertainty spectrum) and reported at the bottom of the panels.

In the text

Fig. 5. Velocity map of the system with ellipses marking the positions of the main regions (these are the same apertures as those used for the photometric analysis in Sect. 3.2 and shown over the NIRCam data in Fig. 2). The uncertainties on the velocity measurements (based on bootstrapping from the variance map of the IFU cube) are on average ∼1 km s−1, and ≤4 km s−1 in all pixels. The map is created combining the signal of all the four main lines observed (Angora et al., in prep., for details). The FWHM of the IFU PSF is given in the bottom-right corner.

In the text

Fig. 6. Imaging and spectroscopy of the faint regions. Left panel: NIRCam sum of all filters (both from the SW and from the LW channels). The apertures used to derive the magnitudes (or upper limits) of the faint regions are shown as elliptical contours. The FoV of the IFU is shown as a white contour. Central panel: IFU observations showing the sum of Hα and [OIII] emission lines with the masks used to extract the spectra shown in the right panel. Two white “X” symbols mark the position of the UT1a and UT1b peaks. Right panel: Spectra from the two masks of the central panel, using their same color-coding. As done in Fig. 2, thick shaded lines are used to show the uncertainties of each spectrum.

In the text

Fig. 7. Mapping of the main spectral properties. Left panel: RGB composite of the LW filters (red: F410M, green: F356W, blue: F277W); the green channel, containing the emission in Hβ+[OIII] highlights the regions with the strongest line emission. The position of emission peaks in NIRCam SW (see also Figs. 3, 8, 9 and 10) are marked by red crosses (by “star” markers in the central and right panels); the bright core regions D1core and T1core discussed in Sects. 4.1 and 4.2 are marked by red asterisks (by plus plus symbols in the central and right panels). The (white) boxes outline the zoom-in regions shown in the other panels of the figure. Central panels: EW maps of [O III]+Hβ, as traced by the F356W filter (with F410M tracing the continuum); the data have been smoothed by a 2 px box kernel and low signal-to-noise pixel have been masked-out. The map has the same pixel-scale as the LW data (0.04 arcsec/px). Right panels:β-slope maps obtained by the pixel-by-pixel fitting of F115W, F150W and F200W. The data are in a 0.02 arcsec/px scale and have been smoothed by a 4 px box kernel. Low-signal pixel have been masked-out. Both in the central and left panels, a black line of 0.24″ is shown, for scale. The FWHM of the PSF in each panel is reported as a black/white circle in the bottom-right angle. Fig. 8. Zoom-in into the D1 region, showing the combined SW observations (F115W+F150W+F200W), the best-fit model of the core region, and its residuals. The apertures used for photometry are shown as white ellipses; for the D1core region, photometry is derived from a light-profile fitting (as described in the main text). For each of the apertures/regions, photometry is shown as black circles, with the SED best-fit model over-plotted as colored empty squares and relative spectrum. The best-fit of the UV slope is also shown and reported in the panels. Fig. 9. Zoom into the T1 system, shown as the sum of observations in the SW filters, the best-fit model of the core emission and its residuals. The results of photometry and of the broad-band SED fitting are shown in the right panel, both for T1core (blue lines) and for T1diff (orange lines), similarly to what done for the subregions of D1 in Fig. 8. Fig. 10. Same as Fig. 9 but for the UT1 system, with the best-fit model of the UT1a source (left panel) and the photometry + broad-band SED fitting results (right panel, green line for UT1a, red line for UT1b).

In the text

Fig. 8. Zoom-in into the D1 region, showing the combined SW observations (F115W+F150W+F200W), the best-fit model of the core region, and its residuals. The apertures used for photometry are shown as white ellipses; for the D1core region, photometry is derived from a light-profile fitting (as described in the main text). For each of the apertures/regions, photometry is shown as black circles, with the SED best-fit model over-plotted as colored empty squares and relative spectrum. The best-fit of the UV slope is also shown and reported in the panels.

In the text

	Fig. 9. Zoom into the T1 system, shown as the sum of observations in the SW filters, the best-fit model of the core emission and its residuals. The results of photometry and of the broad-band SED fitting are shown in the right panel, both for T1core (blue lines) and for T1diff (orange lines), similarly to what done for the subregions of D1 in Fig. 8.
In the text

	Fig. 10. Same as Fig. 9 but for the UT1 system, with the best-fit model of the UT1a source (left panel) and the photometry + broad-band SED fitting results (right panel, green line for UT1a, red line for UT1b).
In the text

Fig. 11. Comparison between the [[O III]λ5007] line emission of the D1 (solid blue) and T1 (solid red) regions (extracted from the spectra presented in Fig. 2) and the Lyα flux observed at the same position (from VLT/MUSE; see Vanzella et al. 2019, 2021), rescaled by a factor 0.5 in order to ease the comparison of the peak relative velocities. For both regions, the center of the [O III]λ5007 line has been used as the systemic redshift and therefore as reference for the x-axis. The Lyα line peaks at Δv ∼ +100 km s−1 in both regions.

In the text

Fig. A.1. Left panel: A slice (corresponding to λ = 4.69 μm) of the final cube produced via the standard JWST-NIRSpec pipeline, showing Hα emission from the T1 region (see also Fig. 2) but also spurious signal blending with the real one (indicated by the blue arrows) and nan pixels (indicated by the white arrows). The same slice, for the cube post-processed by our custom pipeline, is shown in the central panel. Right panel: the residual background, after subtracting a scalar median one (i.e., the dark green background in Fig. A.2), is shown for the same slice; background spatial variations on the order of ∼10−21 erg s−1 cm−2 Å−1 are observed.

In the text

Fig. A.2. Left panel: Spectra extracted in the D1 mask (see Fig. 2) for the cube without background subtraction (light red line) and after the background subtraction (gray line), where the solid black line is the total measured background (based on a moving median across 30 slices) within the mask. The solid red line is the background obtained using one scalar value per slice. Right panel: the difference between the total background and the scalar one is shown as a dark green thick line, in a zoom-in at the wavelengths of Hβ, [O III]λ4959 and [O III]λ5007 emission lines. The residual signal produced using only the scalar background is also shown as a thin green histogram, for reference.

In the text

Fig. B.1. Left panel: the three main regions (D1, T1, and UT1; Sect. 3), along with the two “faint” regions discussed in Sect. 3.5, on top of the sum of all SW and LW filters’ observations. Right panel: All the micro regions discussed in Sect. 4, on top of the sum of the SW filters’ observations. The ellipses marking the extent of the main regions are kept, in transparency. In both panels the FoV of the IFU is shown as a black contour. All the IFU masks used in this work are shown in Appendix G.

In the text

	Fig. E.1. Photometry of the four regions with the steepest UV slopes observed (black markers), along with the best fit model using BPASS (blue diamonds and line) and Bruzual & Charlot (2003) (red squares and line) stellar synthesis population models. On the bottom of each panel the best-fit derived physical properties are given (mass-weighted age, extinction, stellar mass and reduced χ2).
In the text

	Fig. F.1. Spectra of the D1 region covering the wavelengths nearby the Hγ(left panel) and the Hα(right panel) lines, including the best-fits of the observed lines. The left panel also includes a tentative detection of the [OIII]λ4363 line. An artificial 5σ signal at the expected wavelength of the [O III]λ4363 and of the [NII]λ6585 lines is shown as a blue-dashed line.
In the text

Fig. G.1. Left panel: map of the Hα+[OIII] emission from NIRSpec-IFU; the masks used to extract the spectra of the subregions discussed in Sect. 4.1 and Sect. 4.2 are over-imposed. The relative spectra (covering the wavelengths around Hβ, [[O III]λλ4959, 5007] and Hα) are shown in the right panels. All the other IFU masks used in this work (and shown in Fig. 2 and 6) are included for comparison as white contours (with black labels).

In the text