GA-NIFS: the interplay between merger, star formation and chemical enrichment in MACS1149-JD1 at z=9.11 with JWST/NIRSpec

Abstract

We present JWST/NIRSpec integral-field spectroscopy observations of the z $\sim$ 9.11 lensed galaxy MACS1149-JD1, as part of the GA-NIFS programme. The data was obtained with both the G395H grating (R$\sim$ 2700) and the prism (R$\sim$ 100). This target shows a main elongated UV-bright clump and a secondary component detected in continuum emission at a projected distance of 2 kpc. The R2700 data trace the ionised-gas morpho-kinematics in between the two components, showing an elongated emission mainly traced by [O iii]5007. We spatially resolve [O ii]3726,3729, [O iii]4959,5007, and [O iii]4363, which enable us to map the electron density (n${}_{\rm e}\sim 1.0\times 10^{3}$ cm^-3), temperature (T${}_{\rm e}\sim 1.6\times 10^{4}$ K), and direct-method gas-phase metallicity (-1.2 to -0.7 dex solar). A spatially resolved full-spectrum modelling of the prism indicates a north-south gas metallicity and stellar age gradient between the two components. We found 3-$\sigma$ evidence of a spatially resolved anti-correlation of the gas-phase metallicity and the star formation rate density, which is likely driven by gas inflows, enhancing the star formation in JD1. We employ high-z sensitive diagnostic diagrams to rule out the presence of a strong AGN in the main component. These findings show the unambiguous presence of two distinct stellar populations, with the majority of the mass ascribed to an old star formation burst, as suggested by previous works. We disfavour the possibility of a rotating-disc nature for MACS1149-JD1; we favour a merger event that has led to a recent burst of star formation in two separate regions, as supported by high values of [O iii]5007/H, ionised gas velocity dispersion, and gas-phase metallicity.

1 Introduction

James Webb Space Telescope (JWST) observations of high-redshift galaxies are providing unique constraints on the galaxy and interstellar medium (ISM) properties in the very early Universe, up to z $\sim$ 10 (Harikane et al., 2023; Curti et al., 2023; Larson et al., 2023; Isobe et al., 2023; Hsiao et al., 2023a; Robertson et al., 2023; Curtis-Lake et al., 2023; Sanders et al., 2024; Vikaeus et al., 2024; Carniani et al., 2024; Abdurro’uf et al., 2024). At such early epochs, JWST reveals that galaxies are fainter and smaller compared to the local Universe (Suess et al., 2022; Costantin et al., 2023; Huertas-Company et al., 2024). As a consequence, studies aided by gravitational lensing by massive galaxy clusters are crucial to magnify the intrinsic light and sizes of distant, early galaxies, opening a new window on the investigation of their properties (Vanzella et al., 2017, 2022; Meštrić et al., 2022; Mowla et al., 2022; Claeyssens et al., 2023; Hsiao et al., 2023a, b). In the first few billion years of cosmic history galaxies are observed to form giant star-forming clumps, typically resulting in more irregular shapes with respect to present day analogues (Delgado-Serrano et al., 2010; Le Fèvre et al., 2020; Schreiber et al., 2015; Meng & Gnedin, 2021; Förster Schreiber et al., 2011; Ferreira et al., 2020). Observations and simulations show that gas accretion via streams onto massive clumps and during mergers triggers high rates of star formation (SF) and is the preferential driver of galaxy evolution, shaping galaxy morphology and kinematics (Kaviraj, 2014; Zanella et al., 2015; Jackson et al., 2022; Boylan-Kolchin, 2023; Zhang et al., 2023; Huško et al., 2023). Additionally, it is crucial to characterise other drivers of galaxy evolution such as the SF activity, the gas-phase distribution, and the multi-phase gas kinematics to constrain how local galaxies built up at early epochs. The JWST finally allows the exploration of these aspects with superb sensitivity and high angular resolution in the Near-Infrared (NIR).

The goal of this paper is to discuss the ionised gas and stellar population properties in the lensed galaxy MACS1149-JD1 at z = 9.11 employing JWST/NIRSpec Integral Field Spectroscopy (IFS) observations at high and low spectral resolution. JD1 is highly magnified ($\mu$ = 10; Hashimoto et al., 2018; Stiavelli et al., 2023; Bradač et al., 2024) by the foreground cluster of galaxies MACS1149 at z = 0.544 (Zheng et al., 2012). This high-redshift source was first detected through a photometric excess at 4.5 µm with Spitzer/IRAC (Zheng et al., 2012), and it was later confirmed spectroscopically with ALMA (Atacama Large Millimeter Array; Hashimoto et al., 2018; Tokuoka et al., 2022). Hashimoto et al. (2018) presented tentative evidence of Ly emission via VLT/X-shooter observations and [O iii]88µm emission from ALMA coincident with the brighter component. To explain the observed Spectral Energy Distribution (SED) and photometric excess, they suggested the presence of a Balmer break due to a dominant old stellar component (formation redshift z $\sim$ 15) and of a younger component. From SED fitting, Laporte et al. (2021) estimated a best-fit stellar population of age 512 Myr. Moreover, they inferred past star-formation histories (SFH) and showed that more than 93 percent of the observed stellar mass had already formed at z $\sim$ 10. Tokuoka et al. (2022) performed a detailed morpho-kinematic analysis of the ALMA data, discarding the presence of multiple clumps and suggesting the presence of a smoothly rotating gas disk. Consistent with Hashimoto et al. (2018), Tokuoka et al. (2022) found that a $\sim$ 300 Myr old stellar population of 2.3 $\times$ 10⁹ $\mathrm{M_{\odot}}$ is well suited to reproduce the observed features.

The advent of JWST enabled us to drastically revise the properties of this object. Stiavelli et al. (2023) employed JWST NIRCam and NIRSpec MSA observations to first reveal that JD1 is made up of three spatially distinct components; one in the North (hereafter, JD1-N), and two in the South (hereafter collectively JD1-S). They also showed that JD1 is mostly dust free, with a robust metal enrichment (12+log(O/H) = 7.90), and shows at most only a very weak Balmer break. Bradač et al. (2024) presented JWST NIRCam and NIRISS observations of the system; in addition to the three unresolved clumps, they infer the presence of an underlying, edge-on disc component, containing the bulk of the stellar mass and made up of an older stellar population compared to the clumps. Taking advantage of JWST/MIRI MRS observations, Álvarez-Márquez et al. (2024) traced ongoing SF via H$\alpha$ map, with a star formation rate (SFR) of 3.2 – 5.3 $\mathrm{M_{\odot}}$ yr^-1. Moreover, they identified two main clumps with indication of a recent stellar burst ($\geq$ 5 Myr) coincident with JD1-N and a more distributed star-forming region in JD1-S.

To date, JD1 is among the highest-redshift sources for which it is possible to carry out a detailed analysis of the ISM properties via oxygen auroral-lines. In particular, thanks to the NIRSpec spectral coverage, including among many others the [O iii]5007, [O iii]4363 and [O ii]3726,3729 emission lines, we can characterise both the ionised gas kinematics and ISM properties such as the electron density and temperature (Kewley et al., 2019; Morishita et al., 2024; Abdurro’uf et al., 2024; Mazzolari et al., 2024).

In this work, we present spatially resolved gas properties in JD1, using new JWST/NIRSpec Integral Field Unit (IFU) high-resolution grating (R $\sim$ 2700) and low-resolution prism (R $\sim$ 100) observations. This paper is organised as follows. In Sec. 2 we present the NIRSpec observations, together with the data reduction. Section 3 is devoted to the description of emission line analysis and the following results on the ISM properties and gas kinematic. In Sec. 4 we discuss the impact of our results in the broad context of galaxy assembling processes at high-redshift. Finally, in Sec. 5 we summarise the main results. In this paper, we adopt a $\Lambda$CDM cosmology: H₀: 67.4 km s^-1 Mpc^-1, $\Omega_{\mathrm{m}}$ = 0.315, and $\Omega_{\Lambda}$ = 0.685.

2 Observations and data reduction

2.1 NIRSpec IFU Observations

MACS1149-JD1 (hereafter JD1, R.A. 11^h 49^m 33.580^s DEC +22^∘ 24^′45.70^′′, J2000) was observed as part of the Galaxy Assembly with NIRSpec Integral Field Spectroscopy (GA-NIFS, PIs: S. Arribas, R. Maiolino) Guaranteed Time Observations (GTO) programme, included in Program ID $\#$1262 (PI: N. Lützgendorf). The observations were carried out on June 5 2023, using the NIRSpec IFS (Jakobsen et al., 2022; Böker et al., 2022), covering a 3.1^′′ $\times$ 3.2^′′ region, with a native spaxel size of 0.1^′′ (Böker et al., 2022; Rigby et al., 2023). The observations consist of two configurations: high spectral resolution with the G395H grating (R$\sim$ 2700, covering 2.87–5.14 µm) and low spectral resolution with the prism (R$\sim$100, covering 0.6–5.3 µm). Both observations used a medium (i.e. $\sim$ 0.5^′′) cycling pattern of eight dithers. For G395H we used 31 groups per integration and one integration per exposure, with the NRSIRS2 readout mode (Rauscher et al., 2017), giving a total integration time of 18207 seconds (5.05 h). For the prism we used 33 groups per integration and one integration per exposure with the NRSIRS2RAPID readout, giving a total integration time of 3968 seconds (1.1 h).

JD1 resides behind the MACS J1149.6+2223 cluster and was observed by ground- and space-based instruments providing insightful information on the multi-band properties of this system as discussed in Sec. 1. Similarly to the magnification analysis performed in Bradač et al. (2024), we estimated the 2D spatial magnification using different publicly available lens models. Overall, all models show negligible spatial variation of the magnification factor $\mu$ over the spatial scale of JD1, despite on average different lens models can provide very different different median magnification (see Bradač et al. (2024)). Therefore, in this work we adopted a $\mu$ = 10 $\pm$0.7 and corrected for the lensing magnification all the appropriate properties, thus assuming that all flux-dependent quantities are in units of (10/$\mu$) (Hashimoto et al., 2018; Álvarez-Márquez et al., 2024; Stiavelli et al., 2023; Bradač et al., 2024).

Figure 1: JWST NIRSpec collapsed images and spectra before continuum subtraction. Left: G395H (top) and prism (bottom) data cubes collapsed over the total [O iii]5007 emission line and 1.75 - 2.23 µm, respectively. Middle: Zoom-in centred on JD1 of the data cubes shown in the left panels. Yellow filled circles represent the fiducial PSF size. Right: Integrated spectra from JD1-N and JD1-S for the G395H (top) and prism (bottom) data cubes extracted from the apertures shown in middle panels. Labels indicate the identified emission lines. North is up and East is to the left.

2.2 Data Reduction

We performed the data reduction using the well tested JWST calibration pipeline version 1.8.2 with CRDS context jwst1068.pmap. The steps and changes we performed with respect to the standard pipeline in order to improve the final data quality are discussed in Perna et al. (2023). To combine each integration and create the final data cube with spaxel size of 0.05^′′, we adopted the $drizzle$ method, with particular attention to use an official patch to account for a known issue affecting the calibration pipeline. For the G395H data we performed no background subtraction as we will account for its spectral contribution during the emission line fitting. For the low-resolution prism instead, we performed the background subtraction on the detector image, for each of the eight observations separately. From a first datacube, we created a source mask including all contiguous regions with detected [O iii]5007; we enlarged this mask by adding padding of width 0.1 arcsec. We then deprojected this on-sky source mask back onto each of the eight 2-d detector images, using the function blot in the jwst pipeline v1.12.5. For each 2-d ‘cal’ image produced by the pipeline stage 2, we fit a simple linear slope across each of the slices and at each pixel along the dispersion direction, using the 2-d mask to exclude pixels where the source is present. We then smooth the resulting background in the dispersion direction using 7-pixel median filtering. The resulting background images are subtracted from the cal images. We then re-run the pipeline stage 3 on the background-subtracted 2-d images.

3 Analysis and Results

3.1 Emission line fitting

The G395H high-resolution data cube was analysed by means of a set of custom Python scripts with the goal of subtracting the continuum and then fit the observed emission lines (e.g. see Marasco et al., 2020; Tozzi et al., 2021). First, we performed Voronoi binning (Cappellari & Copin, 2003) on the continuum level requiring a minimum S/N in each spectral channel of 5. Then we used the Penalized Pixel-Fitting software (pPXF; Cappellari & Emsellem, 2004; Cappellari, 2023) to fit the continuum in each bin with a second order polynomial while simultaneously fitting the emission lines with one Gaussian component. During this procedure, we masked bad pixels throughout all the data cubes. Then, we subtracted the continuum spaxel by spaxel and matched the spatial resolution within the Field of View (FoV) by convolving each spectral slice of the data cube with a $\sigma$ = 1 pixel ( i.e., 0.05 ^′′) Gaussian kernel. Finally, we fit the emission lines of the continuum-subtracted cube with one or two Gaussians while imposing the same velocity and widths of each component for all the lines and leaving the intensities as a free parameter. For the emission line doublet [O iii]4959,5007 we fixed an intrinsic ratio of 1/3 between the two lines Osterbrock & Ferland (2006). To decide the number of Gaussian components in each spaxel we performed a reduced $\chi^{2}$ analysis with the purpose of identifying any potential region where two components might be needed. In particular, we compared the residuals of the fitting in the wavelength range of the [O iii]4959,5007 doublet using a Kolgomorov-Smirnov test to check if the residuals with an additional Gaussian component are statistically different from those of a model with a single component (see Marasco et al., 2020, for details). As a result, we adopted two Gaussian components only in a few spaxels where a broad line wing is detected. Finally, we obtained an emission-line model cube for each emission line.

The left panels in Fig.  1 show images of the G395H data cube obtained collapsing over the total [O iii]5007 emission and the prism data cube collapsed in the wavelength range 1.75 - 2.23 µm (i.e. 1700 - 2200 Å rest-frame), corresponding to F200W band in NIRCam. While the G395H image is dominated by the nebular emission, the prism one mainly traces the continuum emission with somewhat better angular resolution due to the sharper JWST PSF at these shorter wavelengths. The right panels in Fig. 1 show the integrated spectra of the high- and low-resolution (background subtracted) data cubes extracted from the apertures marked with black circles in the central panels. North and south apertures have radii of 0.1^′′ (i.e. $\sim$ 450 pc) and 0.15^′′ (i.e. $\sim$ 670 pc), respectively. The two apertures correspond to the putative location of the two components identified by Hashimoto et al. (2018) (see also Stiavelli et al., 2023). Besides the known components of JD1, we did not detect any other source in the 3-by-3-arcsec field of view of NIRSpec/IFS. Therefore, in this work, we focused on a subcube centred on the main galaxy and shown in the left and middle columns of Fig. 1. Overall, from the high resolution data cube we fitted the following emission lines: H, H, H, H, H, [O ii]3726, [O ii]3729, He i3889, [Ne iii]3869, [Ne iii]3967, [O iii]4363, [O iii]4959, and [O iii]5007. The NIRSpec G395H spectrum extracted from a pixel in JD1-S is shown in Fig. 2 with the single Gaussian model used to fit the aforementioned emission lines.

Figure 2: JWST/NIRSpec G395H spectrum extracted from a single spaxel of JD1-S. The top panel shows the total spectrum focusing on the 3700 - 5050 Å spectral region, which contains the emission lines of interest. The solid black and red lines show the observed spectrum and best-fit model, respectively. The grey dashed lines show the location of the fitted emission lines. The bottom two panels show a zoom on the regions containing the [O ii] and [Ne iii] doublets, H, H +HeI$\lambda$3888 (left), and the H and [O iii] complex (right). Shaded grey areas are masked spectral channels.

3.2 Gas kinematics

As a result of the spaxel-by-spaxel fit described in the previous section we obtain the one- or two-Gaussian best-fit model in each spaxel. From these models, we create moment maps for all emission lines of interest, which provide precious information on the projected gas kinematics. By comparing these maps and as a consequence of the fit kinematic constraints described in the previous section, we conclude that there is no statistically significant difference between the kinematics of different emission lines. For this reason, in the following we focus only on the [O iii]5007 emission line, which has the highest S/N. Fig. 3 shows the emission-line flux (moment 0), the flux-weighted line of sight (L.O.S.) velocity (moment 1), and the velocity dispersion maps (moment 2) for the [O iii]5007 emission line, computed from the best-fit emission line model cube. The zero velocity is set assuming z = 9.1130 as a reference. Black contours show arbitrary [O iii]5007 flux levels. The peak of the [O iii]5007 emission line originates from the main UV clump (JD1-S), with a blue-shifted tail elongated towards JD1-N. The right panel shows a horizontal central region with high velocity dispersion, which may be due to the beam smearing of the strong velocity gradient. However, the map shows also an intriguing excess of velocity dispersion towards the Western region. The bottom panel of Fig. 3 shows the spectrum extracted from a spaxel in such region of enhanced velocity dispersion, with the best-fit multi-Gaussian model in red. We found that over this region of high velocity dispersion the fits benefit from a second Gaussian component to reproduce the observed [O iii]5007 emission line profile based on $\chi^{2}$ minimisation, likely indicating the presence of outflow or merger. Similarly, Álvarez-Márquez et al. (2024) analysed MIRI/MRS observations tracing the H$\alpha$ kinematics and measured a difference in velocity dispersion of $\approx$ 60 km s^-1 among the two clumps, supporting the hypothesis of outflowing gas powered by the UV-bright clump JD1-S (Tokuoka et al., 2022, see also Sec. 4.3).

From the peak of integrated [O iii]5007 emission line over the two components we measured redshifts of z = 9.1130 $\pm$ 0.0002 and z = 9.1099 $\pm$ 0.0002, for JD1-S and JD1-N, respectively. Therefore, assuming the redshift of JD1-S as reference, JD1-N is blue-shifted by $\sim$ 90 km s^-1. The zero-order moment map in Fig. 3 shows the [O iii]5007 emission of ionised gas connecting the two components. The channel maps of [O iii]5007 in Fig. 4 highlight this blue-shifted elongated tail, showing a clear velocity gradient from the JD1-N to JD1-S, with average projected velocities between -150 and 0 km s^-1.

From these kinematic maps alone, we cannot conclusively distinguish between a rotating disc scenario or a merger. Indeed, in both cases we would expect to see the velocity dispersion peaking halfway between the blue-shifted and red-shifted components (Fig. 3). The kinematic axis of the putative disc seems tilted by $\approx 10^{\circ}$ relative to the major axis of the photometry (which runs approximately along the line connecting JD1-N to JD1-S; Bradač et al., 2024). The possible disc or merger nature of JD1 will be discussed more in-depth in Sec. 4.3.

From the spatially resolved [O iii]5007 emission-line fit we measured an average velocity dispersion of 62 $\pm$ 7 km s^-1, consistent with ALMA (Hashimoto et al., 2018; Tokuoka et al., 2022). To calculate a reference dynamical mass, we assumed the structure of JD1 to be consistent with a disc and used the stellar virial estimator of van der Wel et al. (2022), using our velocity dispersion, a (de-lensed) effective radius of 585${}^{+157}_{-52}$ $\sqrt{10/\mu}$ pc and Sérsic index of 1 $\pm$ 0.2 measured by Bradač et al. (2024) from NIRCam observations, and an axis ratio of 0.6. We obtained a dynamical mass of M_dyn = 1.2${}^{+0.5}_{-0.4}$ $\sqrt{10/\mu}$ $\times$ 10⁹$\mathrm{M_{\odot}}$. This value is consistent with previous H-based estimates of 2.4 $\pm$ 0.5 $\times$ 10⁹ $\mathrm{M_{\odot}}$ from Álvarez-Márquez et al. (2024) and with resolved [O iii]88µm-based estimates of 0.7 - 3.7 $\times$ 10⁹$\mathrm{M_{\odot}}$ of Tokuoka et al. (2022). Using different virial calibrations (e.g., Cappellari et al., 2006, 2013; Wisnioski et al., 2018) would give dynamical masses in the range $0.8\text{--}2.6\sqrt{10/\mu}\times 10^{9}$$\mathrm{M_{\odot}}$, which does not change our conclusions.

Figure 3: [O iii]5007 G395H moment maps and zoom of the [O iii] complex. Top panel: From left to right the integrated emission line flux, the flux-weighted LOS velocity and the velocity dispersion maps. All moment maps are calculated from the total best fit model. Moment maps are masked at S/N smaller than 3. Black contours are arbitrary [O iii]5007 levels. Bottom panel: Zoom of the [O iii]4959,5007 emission line doublet extracted from the spaxel marked with the blue star. Data and best fit model are in black and red, respectively. The two Gaussian components are represented in blue and light blue solid lines.

Figure 4: [O iii]5007 collapsed images at different velocity channels, assuming as zero-velocity the redshift z = 9.1130. From left to right the intervals are: (-150,-100), (-100,-50), (-50,50), (50,100) km s^-1.

3.3 Excitation Diagnostics and Line Ratios

We used the results of the emission line fitting of the high-resolution data cube to examine the excitation source in JD1. Recent works have highlighted the difficulty in separating SF and AGN ionisation at high z and low metallicity by using the standard BPT diagram (Feltre et al., 2016; Nakajima & Maiolino, 2022; Kocevski et al., 2023; Übler et al., 2023; Scholtz et al., 2023; Maiolino et al., 2023). We therefore explored the possible source of ionisation by employing the new diagnostic diagrams proposed by Mazzolari et al. (2024), which use the ratio of [O iii]4363/H to identify narrow-line Type II AGN at high redshift. Fig. 5 shows one of such diagnostic diagrams using the [O iii]4363/H vs [O iii]5007/[O iii]4363 line ratios, together with a map of the most likely excitation source in JD1 and the [O iii]4363/H line ratio. In this diagram, being above the demarcation line is a sufficient but not necessary condition for an object to be identified as an AGN, therefore spaxels which are below the demarcation line can still be potentially associated with AGN excitation. The diagram shows that no regions in JD1 show unambiguous evidence for AGN. Beside the one shown in Fig. 5 we also explored the other two diagrams discussed in Mazzolari et al. (2024) (employing [O iii]5007/[O iii]4363 and [Ne iii]3869/[O ii]3726,3729 emission line ratio vs [O iii]4363/H) and observed that all the selected spaxels in each diagram do not show unambiguous evidence for an AGN, favouring the SF ionisation as most likely scenario. For completeness, we estimated the position of JD1 in the standard BPT diagram adopting an H flux of 1.05 $\pm$ 0.07 $\times$ 10^-17 erg s^-1 cm^-2, a 3$\sigma$ upper limit flux on [N ii] of 2.1 $\times$ 10^-18 erg s^-1 cm^-2 from Álvarez-Márquez et al. (2024) and an average log([O iii]/H) = 0.8 ratio from our spatially resolved analysis. Similarly to the result obtained with the diagnostic diagram shown in Fig.  5, the standard BPT diagram report an ambiguous ionisation nature for JD1, with the location of this system being located on the demarcation line between SF and AGN ionisation. Therefore, from the results of both diagnostics we cannot exclusively state the nature of the main ionisation source in JD1.

In order to trace regions of high ionisation we also computed the [O iii]5007/H line ratio map shown in Fig. 6. Due to the difference in ionisation potential of [O iii]5007 and H this ratio is a good indicator of the hardness of the radiation field, and thus of instantaneous bursts of SF. We found that JD1-N is characterised by high excitation, with [O iii]5007/H $\sim$ 8, at variance with JD1-S which is characterised by an average [O iii]5007/H $\sim$ 6.

3.4 Electron density and Temperature

The high spectral resolution G395H data includes the spectrally resolved [O ii]3726,3729 doublet and [O iii]4343 emission lines, which allow us to compute the first spatially resolved estimate of the electron density and temperature at z$\geq$ 9. In particular, first we used the pyneb Python package (Luridiana et al., 2012, 2015) to compute the electron temperature T_e in the [O iii]-emitting region from the [O iii]5007/[O iii]4363 line ratio assuming no dust attenuation (Hashimoto et al., 2018; in any case, even a dust attenuation $A_{V}=1$ mag would only change T_e by 20 percent) and an electron density of 600 cm^-3 (which is a good assumption based on the value measured in high-z JWST galaxies; Isobe et al., 2023; Abdurro’uf et al., 2024). To compute the temperature of [O ii]-emitting regions (T₂), due to the absence of temperature diagnostics, we used the relation T₂ = -0.744 + T_e $\times$ (2.338 - 0.610$\times$ T_e) from Izotov et al. (2006) to derive it from the temperature of the [O iii]-emitting regions.

Using T₂ and the observed [O ii]3726/[O ii]3729 ratio, we then compute the electron density. Fig. 7 shows spatially resolved maps of the electron gas density (left panel) and temperature (middle panel). The derived gas-phase metallicity computed with the direct T_e method is discussed in the Sec. 3.6, together with an estimate via prism continuum fitting. The average electron temperature and density of JD1-S are 1.4${}^{+0.6}_{-0.7}$ $\times$ 10⁴ K and 670${}^{+100}_{-240}$ cm^-3, with the latter being consistent with the assumed density to compute T_e . As shown in Fig. 7, due to the low S/N of the [O ii] doublet over JD1-N we could not compute a spatially resolved estimate of the electron density. Therefore, we integrated the spectrum from the aperture of JD1-N (see Fig. 1) and assuming an electron temperature of 10⁴ K we estimated an electron density from integrated spectrum of 730 cm^-3, which is slightly higher than the integrated value of JD1-S but still consistent within the uncertainties. We computed the uncertainties on the electron density on the basis of a Monte Carlo technique. Our findings are consistent with the observed redshift evolution of the electron density reported by Isobe et al. (2023, see also , ). In particular, it is possible that, due to the more compact morphology of gas nebulae and the lower metallicity at high redshift, the electron density is higher than in the local Universe (see Sec. 4 for a more detailed discussion).

Figure 5: AGN-diagnostic diagram of Mazzolari et al. (2024) for the excitation source in JD1 based on the [O iii]4363/H vs [O iii]5007/[O ii]3727 line ratio. Dashed black line is the demarcation between pure AGN (above) and AGN or SF excitation source (below).

3.5 Continuum fitting - Stellar population properties

To provide a spatially resolved analysis of the stellar mass in JD1, we performed the continuum fitting of the prism data cube both spaxel by spaxel and by integrating the total spectrum using prospector (Johnson et al., 2021), a Bayesian SED modelling framework built around the stellar-population synthesis tool fsps (Conroy et al., 2009; Conroy & Gunn, 2010). We performed a SED modelling of the observed spectrum from 0.6 to 5.27 µm following the procedure described in Tacchella et al. (2023) and Pérez-González et al. (2023). Using the prism data cube we performed a spatial smoothing, to obtain approximately the same PSF at all wavelengths. We assumed a 2-d Gaussian kernel, with different full-width half-maximum values along and across slices, as a function of wavelength (D’Eugenio et al., 2023).

For the modelling, we configured fsps to use the MILES stellar atmospheres (Falcón-Barroso et al., 2011) and MIST isochrones (Choi et al., 2016). The nebular emission is modelled using pre-computed grids from cloudy (Ferland et al., 1998), as described in Byler et al. (2017); this approach takes into account possible stellar absorption at the wavelength of emission lines (see e.g. Pérez-González et al., 2003, 2008). Finally, we accounted for dust attenuation using a flexible dust attenuation law, consisting of a modified Calzetti law (Calzetti et al., 2000) with a variable power-law index and UV-bump strength (Noll et al., 2009; Kriek & Conroy, 2013). Stars younger than 10 Myr are further attenuated by an extra dust screen, parametrised as a simple power law (Charlot & Fall, 2000). No dust emission is included, because our reddest wavelengths are still in the rest-frame optical range, where dust emission is negligible. The star-formation history (SFH) uses 9 fixed time bins between z = 9.1130 and z = 20; the first three bins are at 10 Myr, 30 Myr, and 100 Myr, the remaining 6 bins are logarithmically spaced in time. We use a continuity prior to relate the log ratio of the SFRs between adjacent time bins (Leja et al., 2019). The model free parameters and their prior probability distributions are listed in Table 1. The left panel in Fig. 8 shows the spatially resolved stellar mass surface density estimated with this method, together with the spatially resolved maps of the SFR surface density within the last 10 (SFR₁₀) and 100 Myr (SFR₁₀₀).

From the integrated spectrum of JD1, we found a total stellar mass budget of $\log(\mathrm{M_{\star}}/\mathrm{M_{\odot}})_{\rm JD1}=7.47^{+0.05}_{-0.05}\log(10/\mu)$, which is on average an order of magnitude lower compared to previous estimates of the stellar mass in JD1 (Laporte et al., 2021; Stiavelli et al., 2023; Bradač et al., 2024). Such a significant discrepancy has to be ascribed to differences in the data. Our prism spectrum displays a clearly detected continuum, with evidence for a Balmer jump due to nebular emission (Fig.1). In contrast, other studies find evidence of a Balmer break, surmising an old stellar population which greatly contributes to the total stellar mass. In two cases, the Balmer break is seen only from photometry (Laporte et al., 2021 and Bradač et al., 2024, their fig. 2G), whereas in Stiavelli et al. (2023) the Balmer break is not clearly seen in the spectroscopy (their fig. 2). Also, our estimate of the total stellar mass is 0.4 dex lower than the value we obtain from modelling the spectrum in each spaxel and then adding up the resulting stellar masses (which gives $\log(\mathrm{M_{\star}}/\mathrm{M_{\odot}})=7.88\log(10/\mu)$). This discrepancy may be due to outshining of fainter, older components by younger stars in the integrated flux (Giménez-Arteaga et al., 2023). From the spatially resolved kinematic analysis (see Sec. 3.2) and assuming a rotating disc, we estimated a (de-lensed) maximal dynamical mass of M_dyn = 1.2${}^{+0.5}_{-0.4}$ $\sqrt{10/\mu}$ $\times$ 10⁹$\mathrm{M_{\odot}}$. This value is much higher than the total stellar mass in JD1 and is possibly due either to a large cold gas content, which is feeding the ongoing SF (e.g. Álvarez-Márquez et al., 2024), or, alternatively, to non-virialised motions (we discuss this finding in Sec. 4) or large amounts of dark matter (de Graaff et al., 2024) or a combination of the aforementioned effects.

Figure 6: [O iii]5007/H emission line ratio in JD1 in logarithmic scale. Black solid line are arbitrary [O iii]5007 levels. Spaxels with S/N $\leq$ 3 are masked.

Figure 7: ISM physical properties of JD1. From left to right: Electron density, temperature, gas-phase metallicity derived with the direct method (T_e method) and gas-phase metallicity derived with prospector performing SED fitting (see Sec. 3.4 and Sec. 3.5 for details). Spaxels at S/N $\leq$ 3 for the [O ii]3726,3729 doublet are masked. Black contours are arbitrary [O iii]5007 flux levels. The black cross marks the position of the [O iii]5007 line peak.

From the SED fitting of the integrated spectrum we derived a SFR${}_{10}=2.8^{+0.1}_{-0.1}$ $\mathrm{M_{\odot}}$ yr^-1 and SFR${}_{100}=0.29^{+0.1}_{-0.1}$ $\mathrm{M_{\odot}}$ yr^-1, integrating the SF history (SFH) of the last 10 and 100 Myr, respectively. Combining the result for the total stellar mass formed and the SFR, we computed the sSFR within the last 10 and 100 Myr. The resulting spatially resolved maps are shown in Fig. 9. From the integrated spectrum instead, we derived specific SFR of sSFR₁₀ = 9.6${}^{+1.9}_{-1.7}$ Gyr^-1 and sSFR₁₀₀ = 0.97 ${}^{+0.2}_{-0.2}$ Gyr^-1. The sSFR maps show consistent values between JD1-N and JD1-S in the last 100 Myr. Within the last 10 Myr instead, we observe higher sSFR over JD1-S compared to JD1-N. Nonetheless, comparing the sSFR within the last 10 and 100 Myr for JD1-N, we conclude that the majority of its mass formed very recently (Álvarez-Márquez et al., 2024; Stiavelli et al., 2023; Bradač et al., 2024), which could explain the observed enhancement of the gas-phase metallicity (in absence of significant inflows of low-metallicity gas; see next section and right panels of Fig. 7).

As a comparison with the SFR derived via SED fitting, we measured the H flux to provide a direct estimate of the SFR in JD1. We assumed a Case B recombination scenario (Osterbrock & Ferland, 2006) and based on the measured H flux of 1.47 $\pm$ 0.2 erg s^-1 cm^-2, we computed the SFR as:

$$SFR(\text{H\textbeta})=5.5\times 10^{-42}\ L_{\text{H\textbeta}}(\mathrm{erg\,s^{-1}})\times f_{\text{H\textalpha}/\text{H\textbeta}}\ \mathrm{M_{\odot}}\mathrm{yr}^{-1},$$ (1)

where $f_{\text{H\textalpha}/\text{H\textbeta}}$ = 2.85 for T_e = 1.4 $\times$ 10⁴ K (see Sec. 3.4. Correcting for the lensing factor we obtain $\rm SFR(\text{H\textbeta})$ = 2.5 $\pm$ 0.3 (10/$\mu$) $\mathrm{M_{\odot}}\mathrm{yr}^{-1}$, which is consistent with the value derived via SED fitting.

Fig. 10 shows maps of estimated rest-frame equivalent width (EW) of the H and [O iii]5007 emission lines. Both maps show a mild decrease of EW from the centre to JD1-S outwards. A similar trend of radially decreasing EW of H in high-redshift galaxies has already been reported (Tripodi et al., 2024). We also find higher EW values over JD1-S compared to JD1-N, as also shown by the sSFR maps in Fig. 9. Both findings points towards a higher and more recent SF in JD1-S, consistent with the result obtained via SED fitting. Moreover, as discussed in Sect. 4.3, JD1-N may be characterised by an extreme escape fraction of ionising photons. The effect of an elevated escape fraction over JD1-N translates in a loss of ionising photons that in turn reduces the intensity of emission lines with respect to the continuum, decreasing the observed EW. We note that the SED fitting algorithm that we used for the continuum analysis is not tailored to account for any fraction of escaping photons. Similarly to the effect on the EW, including a higher escape fraction in JD1-N would enhance the estimated resolved SFR and sSFR density maps, consistently with a scenario of a recently formed stellar population.

Fig. 11 shows a map of the $\beta_{UV}$ slope, measured between 1200 to 1500 Å rest-frame (where $\beta_{UV}$ is defined as $F_{\lambda}\propto\lambda^{\beta_{UV}}$), computed via continuum fitting of the prism data cube. Interestingly, we observed a smooth north–south gradient, highlighting different properties of the two main companions. The $\beta_{UV}$ has a maximum of –2.3 over JD1-S and then decreases to a steeper slope of –2.6 towards JD1-N. Given that we find no evidence of dust from the emission-line ratios, these differences in UV slope likely indicate an older and younger stellar population for JD1-N and JD1-S, respectively (Stiavelli et al., 2023; Bradač et al., 2024). At variance with previous works, we do not detect any Ly$\alpha$ emission in the low-resolution data cube, probably due to the low spectral resolution of our prism data (Hashimoto et al., 2018; Bradač et al., 2024). Similarly to Bradač et al. (2024), we injected in our low-resolution data the Ly emission line with an integrated (lensed) flux of 4.3 $\pm$ 1.1 $\times$ 10^-18 erg s^-1 cm^-2, as measured by Hashimoto et al. (2018). Based on the assumed flux estimate the line should be detected in our prism data at a S/N $\sim$ 9. However, we do not detect any emission line at the expected wavelength (see Fig.  1) and thus conclude that there is no evidence for Ly emission (Bradač et al., 2024, see also).

Fig. 12 shows the SFH derived for each spaxel and for the integrated JD1-N and JD1-S clumps. The SFH highlights that JD1-N is a very recently formed component, with the majority of the stellar mass formed within the last 10 Myr, thus characterised by a very recent peak of SF. Similarly, despite the flatter SFH of JD1-S we observed a strong peak of SF at the source redshift, followed by a flat SFH between z = 9.2 and 10.6, and then a significant decrease to higher redshift. At variance with Bradač et al. (2024), we do not detect any secondary peak of SF at high redshift for any clump in JD1. Despite this, overall our SFH is qualitatively consistent with previous analysis, showing a recent SF burst (Laporte et al., 2021; Bradač et al., 2024).

Table 1: Summary of the parameters, prior probabilities and posterior probabilities of the fiducial SED fitting prospector model (see also Fig. 8). The values reported here are not corrected for lensing magnification.

	Parameter	Free	Description	Prior	PosteriorJD1
	(1)	(2)	(3)	(4)	(5)
Free parameters	$\log\mathrm{M_{\star}}[\mathrm{M_{\odot}}]$	Y	total stellar mass formed	$\mathcal{U}(7,13)$	$8.47^{+0.05}_{-0.05}$
	$\log Z[\mathrm{Z_{\odot}}]$	Y	stellar metallicity	$\mathcal{U}(-2,0.19)$	$-1.62^{+0.10}_{-0.10}$
	$\log\mathrm{SFR}$ ratios	Y	ratio of the $\log\mathrm{SFR}$ between adjacent bins of the SFH	$\mathcal{T}(0,0.3,2)$	—
	$\tau_{V}$	Y	optical depth of the diffuse dust	$\mathcal{G}(0.3,1;0,2)$	$0.63^{+0.08}_{-0.10}$
	$\mu$	Y	ratio between the optical depth of the birth clouds and $\tau_{V}$	$\mathcal{U}(-1.0,0.4)$	$0.02^{+0.05}_{-0.02}$
	$\sigma_{\mathrm{gas}}\;[\mathrm{km\,s^{-1}}]$	Y	intrinsic velocity dispersion of the star-forming gas‡	$\mathcal{U}(0,300)$	$64^{+44}_{-38}$
	$\log Z_{\mathrm{gas}}[\mathrm{Z_{\odot}}]$	Y	metallicity of the star-forming gas	$\mathcal{U}(-2,0.19)$	$-0.79^{+0.05}_{-0.05}$
	$\log U$	Y	ionisation parameter of the star-forming gas	$\mathcal{U}(-4,-1)$	$-1.02^{+0.01}_{-0.02}$
Other	$\log SFR_{10}[\mathrm{M_{\odot}}\,\mathrm{yr^{-1}}]$	N	star-formation rate averaged over the last 10 Myr	—	$1.46^{+0.04}_{-0.05}$
	$\log SFR_{100}[\mathrm{M_{\odot}}\,\mathrm{yr^{-1}}]$	N	star-formation rate averaged over the last 100 Myr	—	$0.47^{+0.04}_{-0.05}$

(1) Parameter name and units (where applicable). (2) Only parameters marked with ‘Y’ are optimised by prospector; parameters marked with ‘N’ are either tied to other parameters (see Column 4), or are calculated after the fit from the posterior distribution (in this case, Column 4 is empty). (3) Parameter description. For the dust attenuation parameters $n$, $\tau_{V}$ and $\mu$ see Tacchella et al. (2022, their eq.s 4 and 5). (4) Parameter prior probability distribution; $\mathcal{N}(\mu,\sigma)$ is the normal distribution with mean $\mu$ and dispersion $\sigma$; $\mathcal{U}(a,b)$ is the uniform distribution between $a$ and $b$; $\mathcal{T}(\mu,\sigma,\nu)$ is the Student’s $t$ distribution with mean $\mu$, dispersion $\sigma$ and $\nu$ degrees of freedom; $\mathcal{G}(\mu,\sigma,a,b)$ is the normal distribution with mean $\mu$ and dispersion $\sigma$, truncated between $a$ and $b$. (5) and (6) Median and 16^th–84^th percentile range of the marginalised posterior distribution for the north and south clump, respectively; for some nuisance parameters we do not present the posterior statistics (e.g., log SFR ratios). ^‡ The velocity dispersion of the emission lines is a nuisance parameter, due to the low spectral resolution of the prism data ($\sim 500$ $\mathrm{km\,s^{-1}}$ at the wavelength of [O iii]5007); indeed, the posterior probability distributions for JD1 is both consistent with 0 $\mathrm{km\,s^{-1}}$.

Figure 8: JD1 properties derived from SED fitting of the low-resolution data cube with prospector (see Sec. 3.5) not corrected for the gravitational lensing effect. From left to right the surface density of the total stellar mass formed, the SFR within the last 10 Myr and 100 Myr. Black contours are arbitrary flux levels from the prism data cube collapsed over the NIRCam F200W filter.

3.6 Gas-phase metallicity

In Sec. 3.4 we derived spatially resolved estimates of the ISM electron density and temperature for JD1-S from [O ii]3726/[O ii]3729 and [O iii]4363/[O iii]5007 line ratios, respectively (see left panels in Fig. 7). Then, we used the getIonAbundance routine of pyneb, to compute the ionic abundances of the oxygen ions O²⁺ and O⁺ given their respective electron temperatures, gas densities, and the flux ratio of [O iii]5007/H and ([O ii]3726 +[O ii]3729)/H. Here, we approximate the oxygen abundance O/H as (O²⁺ + O⁺)/H, neglecting higher ionisation species of O and assume that it is representative of the total gas metallicity (for a detailed discussion on the underlying assumptions of this so called T_e method see Maiolino & Mannucci, 2019). The third panel in Fig. 7 shows the direct metallicity map of spaxels with resolved [O ii] doublet and [O iii]4363 emission line measurements. The average value of the gas-phase metallicity of JD1-N and JD1-S is computed from integrated spectra of the selected aperture shown in Fig. 1. We computed 12 + log(O/H) = 8.3$\pm$ 0.1 and 7.71 $\pm$ 0.25, for the North and South apertures, respectively, consistently on average with the resolved values (see Fig.  7).

The right panel of Fig. 7 shows the metallicity map as derived with prospector (i.e. using its grid of photionization modelling) from the nebular fitting of the prism cube. Interestingly, we derive consistent values for the gas metallicity derived via prospector photoionization models fitting of the prism spectra and from direct metallicity measurement using the grating spectra. As a more quantitative comparison, from the prism SED fitting we computed integrated gas metallicities for JD1-N and JD1-S of log(O/H) = 8.00${}^{+0.02}_{-0.013}$ and 7.85${}^{+0.03}_{-0.03}$, respectively, which are only slightly below (for JD1-N) and above (for JD1-S) the result from the direct method (see Fig.  7 and Sec. 3.5). The metallicity derived via prism fitting of the integrated spectrum of JD1 is log(O/H) = 7.89${}^{+0.05}_{-0.05}$, consistent with the value of 7.90${}^{+0.04}_{-0.05}$ reported by Stiavelli et al. (2023). Nevertheless, we observe a North-South gradient for the metallicity content from Fig. 7. JD1-N appears to be more metal enriched compared to JD1-S.

Then, we used pyneb to compute the Neon abundance and derived spatially resolved maps of Log(Ne/O) and Log(Ne/H), as shown in Fig. 13, from [Ne iii]3869, H and [O iii]5007 emission lines. Our measurements show higher abundances of Ne/H towards JD1-N. Similarly, the Ne/O abundance has a peak co-spatial to the Ne/H, without showing a second maximum towards JD1-N. Izotov et al. (2006) showed a correlation between log(Ne/O) and the gas-phase metallicity, claiming that in high-metallicity H ii regions the oxygen is depleted onto dust grains. Since JD1 is supposed to be almost dust free we should not expect a trend of log(Ne/O) with the gas-phase metallicity. Overall, our estimates of the Ne abundances are consistent with previous analysis (Stiavelli et al., 2023).

Many works have shown evidence of a three-dimensional relationship between M_⋆-Z_gas-SFR that does not appear to evolve with redshift (Tremonti et al., 2004; Mannucci et al., 2010; Maiolino & Mannucci, 2019; Cresci et al., 2019). Recently, Baker et al. (2023) have provided evidence of a spatially resolved version of this so-called Fundamental Metallicity Relation (FMR), showing that the local metallicity depends primarily on $\Sigma_{M_{\star}}$ and anti-correlates with $\Sigma_{SFR}$. We employed the resolved maps of SFR₁₀₀ density and the gas-phase metallicity derived from prospector to test such observed anti-correlation at z$\geq$ 9. Fig. 14 shows the observed resolved anti-correlation between the gas phase metallicity and the SFR₁₀₀ in JD1. As shown in Fig. 14, we performed a Spearman correlation test and found a correlation coefficient $\rho$ = -0.3 and p-value = 0.001, clearly indicating a significant, but mild, correlation between the SFR₁₀₀ density and the gas metallicity. Overall, our findings are qualitatively consistent with the effect of a stream of gas inflowing towards JD1-S from JD1-N (see Sec. 3.2). In this scenario, the inflowing gas dilutes the metal content of the ISM, decreasing the metallicity (see right panels in Fig. 7) and replenishing the gas reservoir of JD1-S, thus increasing the SFR (see Figs. 8, 10). Similarly, Arribas et al. (2023) found a prominent north-south metallicity gradient in a z $\sim$ 6.9 galaxy and interpreted it as the result of accretion of metal poor gas from the circumgalactic medium (CGM) into the galaxy (see also Dekel et al., 2009; Cresci et al., 2010; Rodríguez Del Pino et al., 2024). Similar results were also found for a sample of compact galaxies studied with NIRSpec/MSA (Tripodi et al., 2024).

Figure 9: Specific Star formation rate (sSFR) within the last 10 (Left) and 100 (Right) Myr of JD1 not corrected for the gravitational lensing effect. Black contours are arbitrary flux levels from the prism data cube collapsed over the NIRCam F200W filter, as in Fig.  8.

Figure 10: Equivalent width maps of H (left) and [O iii]5007 (right) in logarithmic scale. Black contours are arbitrary flux levels from the prism data cube collapsed over the NIRCam F200W filter, as in Fig.  8

Figure 11: Continuum slope map derived from prism data cube fit in the spectral rest-frame range 1200 to 1500 Å (see Sec. 3.5). Black contours are arbitrary flux levels from the prism data cube collapsed over the NIRCam F200W filter.

4 Discussion

In this paper we presented the spatially resolved properties of JD1, made possible thanks to the combination of high spectral resolution and broad spectral coverage of the G395H and prism NIRSpec/IFS data. Here we summarise and discuss the connection between the north and south companions, JD1-N and JD1-S, also by comparing our results with previous works.

4.1 Structure and interaction between the JD1 components

Fig. 1 clearly shows two different components in the collapsed continuum, separated by $\sim$ 2 kpc in projection (Hashimoto et al., 2018; Stiavelli et al., 2023; Álvarez-Márquez et al., 2024; Bradač et al., 2024). ALMA, HST, and JWST observations support the hypothesis of two different stellar populations in JD1 (Hashimoto et al., 2018; Stiavelli et al., 2023). In particular, the most credited scenario is that about half of the current stellar mass of JD1-S formed $\sim$130 Myr before the observed epoch. Then, the galaxy gas reservoir has been replenished by inflows that lead to an ongoing secondary burst of SF, which boosted the observed [O iii]88µm and Ly emission. Here we speculate that the secondary burst of SF was triggered by the encounter of JD1-S with JD1-N, which is transferring substantial amounts of gas to sustain the recent SF episode in JD1-S. From the kinematic point of view, this idea is consistent with projected gas kinematics in Fig. 3 and 4. Indeed, as clearly shown in the second panel in Fig.  4, we observed a collimated filament of ionised gas connecting the two components, with the highest velocity slightly offset from JD1-N and directed towards south. In this scenario we are witnessing the merger between JD1-N and the more massive JD1-S, with gas being dragged towards the gravitational centre. JD1-N has the highest blue-shifted velocity as it is above the plane of the sky with respect to the major galaxy, with the filament representing accreting gas towards South, still dragged by JD1-N, as supported by channel maps in Fig. 4, thus generating the extended blue-shifted tail.

The inflowing gas moving from North to South replenishes the gas reservoir of the system, possibly sustaining further SF. The panels in Fig. 8 show a peak of the stellar mass density and the SFR density within the last 100 Myr co-spatial with the enhanced velocity dispersion observed in Fig. 3. This peak is slightly offset from the peak of the [O iii]5007 emission line and is likely indicating that recent SF is powering outflowing gas, thus increasing the gas velocity dispersion or alternatively is due to beam smearing.

Figure 12: Star formation histories of single spaxels (purple) and integrated JD1-N (green) and JD1-S (maroon) clumps. Errorbars for the main clumps represent the 84^th and 16^th percentiles.

Figure 13: Neon abundance in JD1 with respect to Hydrogen (left) and Oxygen (right). Black solid lines are arbitrary [O iii]5007 levels. Spaxels with S/N $\leq$ 3 are masked.

Interestingly, the right panel in Fig. 7 shows a peak of the gas-phase metallicity over JD1-N (12 + log(O/H) $\sim$ 7.95), consistent with the integrated value obtained from the direct method. The higher metallicity and [O iii]5007/H line ratio support the idea of a SF burst over JD1-N which is enriching the ISM with newly forming metals. In particular, many works show that high values of [O iii]5007/H can be produced by high levels of sSFR (see Fig. 9), which in turn enhances the ionisation parameter and thus indicates high SF efficiency (Kewley et al., 2016; Sanders et al., 2016a; Dickey et al., 2016; Reddy et al., 2023). Finally, the low sSFR observed over the peak of the stellar mass within the last 100 Myr (see Fig.  9) indicates that most of the mass of the system must be older than 100 Myr, consistent with previous analysis (Hashimoto et al., 2018; Stiavelli et al., 2023; Álvarez-Márquez et al., 2024; Bradač et al., 2024). Our findings are consistent with a younger generation of stars distributed over JD1-N, and an older stellar population in JD1-S raised by the replenishing of gas dragged by this companion (see Fig. 8, 9). This scenario is also supported by VLT/X-Shooter observations of a diffuse Ly bubble, supposedly inflated by UV emission of the older population in the JD1-S (Hashimoto et al., 2018).

Figure 14: Distribution of the spatially resolved gas-phase metallicity and SFR density colour-coded with the surface density of the total stellar mass formed in JD1. The SFR density and stellar mass surface density are not corrected for magnification effects. All properties are derived with prospector (see Sec. 3.5 and Figs. 7, 8).

4.2 Resolved electron density at z = 9.11

Taking advantage of the high spatial and spectral resolution of the G395H data cube, we estimated the electron density from the integrated spectra of the two main companions in JD1 and provide the first spatially resolved electron density map at high redshift of the main galaxy. Fig. 15 shows a collection of electron density estimates from the local Universe up to redshift 8.68 compiled by Isobe et al. (2023) with the addition of electron density estimates at z $\sim$ 3.5 by the GA-NIFS project (Lamperti et al., 2024; Rodríguez Del Pino et al., 2024) and z$\sim$ 10.1 by Abdurro’uf et al. (2024). The black dashed and dotted lines represent the n_e $\propto$ $\times$ (1+z)^p relation, with p = 2 and 1, respectively. Our estimates for the electron density over the integrated JD1-N, JD1-S components, as well as the total integrated value for JD1 and the values derived from the spatially resolved analysis, represented as scattered blue points, are in agreement with the predicted trend. In particular, by using high-z JWST observations and extrapolating the trend of low-z galaxies up to z $\sim$ 9, Isobe et al. (2023) suggest that the electron density up to z = 9 is well fitted by n_e $\propto$ $\times$ (1+z)^p, with p = 1 - 2. A possible scenario to explain the observed trend (Isobe et al., 2023; Abdurro’uf et al., 2024) relates both to the decrease of diffused metals (Davé et al., 2011) into the ISM and to the redshift evolution of galaxy sizes with redshift (Kartaltepe et al., 2023; Ormerod et al., 2024). In particular, as suggested by Isobe et al. (2023, see also , ), the gas density of SF regions at high redshift are higher with respect to local analogues due to their higher compactness, due to the galaxy sizes scaling with redshift as R $\propto$ (1+z)^-n, with $n\sim$ 0.7 - 1 (Shibuya et al., 2015; Ono et al., 2023; Ormerod et al., 2024). Assuming a disc-like configuration with redshift-independent thickness, we would get n_e $\propto$ R^-2 which gives the n_e $\propto$ $\times$ (1+z)${}^{1.4\text{-}2}$ relation. On the other hand, recent works motivate the n_e $\propto$ $\times$ (1+z)¹ relation as mainly due to the metallicity evolution with redshift (Abdurro’uf et al., 2024). In particular, hydrodynamic simulations of high-z galaxies revealed a clear anti-correlation between the mean cloud density and the gas metallicity ($Z$) described by $\bar{n}\ \propto Z^{-1}$ (Garcia et al., 2023), as due to the high temperatures and low metallicity of high redshift systems (Kewley et al., 2019). Combining the mean gas density evolution with the observed and predicted gas metallicity trend with redshift: $Z\propto(1+z)^{-1}$ (Davé et al., 2011; Yuan et al., 2013), we get the dotted relation: n_e $\propto$ $\times$ (1+z)¹.

Figure 15: Electron density n_e as a function of redshift, adapted from Isobe et al. (2023) (see references therein). Dotted and dashed curves represent n_e following the $\propto$(1 + z) and $\propto$(1 + z)² relations, respectively. Blue, green and red crosses represent the integrated values for JD1-N, JD1-S and for the total integrated emission, respectively. Small orange circles represent the electron density of single spaxels as derived from our spatially resolved analysis. Blue squares and green hexagons show median n_e of SFMS galaxies at z $\sim$ 0–1 (Berg et al., 2012; Swinbank et al., 2019; Davies et al., 2021) and z $\sim$ 1–3 (Steidel et al., 2014; Sanders et al., 2016a; Kaasinen et al., 2017; Kashino et al., 2017; Davies et al., 2021), respectively. Cyan diamond and yellow pentagon show the average n_e from high specific SFR galaxies at z =0 (Berg et al., 2022) and z =1-3 (Christensen et al., 2012b, a; Sanders et al., 2016b; Gburek et al., 2019), respectively. Small blue circles and pink crosses show single JWST galaxies at z $\sim$ 4–6 and z $\sim$ 7–9, respectively (Isobe et al., 2023). The large purple circle and double circle denote the 50th and 16th–84th percentiles of the properties of the JWST galaxies at z $\sim$ 4–6 and z $\sim$ 7–9, respectively. Orange squares show n_e derived from single clumps and stacked spectra of a single JWST z$\sim$ 10.1 lensed galaxy from Abdurro’uf et al. (2024). The black and grey circles are electron density estimates at z $\sim$ 3.5 by GA-NIFS observations (Rodríguez Del Pino et al., 2024; Lamperti et al., 2024).

4.3 On the rotating disc nature and the third clump of JD1

ALMA [O iii]88µm observations of JD1 show a smooth gas distribution, with no resolved clumps within the system. Tokuoka et al. (2022) computed spatially resolved kinematics of the [O iii]88µm line that shows a smooth north–south velocity gradient, as we observed from [O iii]5007 emission line, and concluded that the best scenario for JD1 is a rotating disc galaxy. JWST NIRCam F150W and F200W observations bring evidence of two smaller clumps building up the main JD1-S clump JD1-Sa and JD1-Sb, see Fig.  16; Stiavelli et al., 2023; Bradač et al., 2024. In particular, Bradač et al. (2024) identified three unresolved star-forming clumps with an underlying smooth galaxy component. From our analysis, the high-resolution data does not allow us to resolve the JD1-N clump, which instead is well detected in the low-resolution data. Interestingly, we observe that the position of JD1-Sb coincides with the peak of the [O iii]5007 emission line, whereas JD1-Sa is shifted towards North by $\approx$ 450 pc, thus approximately located along the enhanced velocity dispersion region identified by [O iii]5007 emission line (see Fig. 3). We investigated the role of JD1-Sa and JD1-Sb in the big picture of the kinematic in JD1 to test whether the disc nature of the system can be confirmed or not. Assuming that JD1-Sb is the centre of the rotation is actually inconsistent with the observed kinematics, as there is no positional agreement with the peak of velocity dispersion shown from [O iii]5007 analysis (see Fig. 3). Additional insights come from the analysis of the escape fraction of ionizing photons. Recent works show evidence of high escape fraction of UV photons in bright [O iii]5007 line emitters (Izotov et al., 2018; Fletcher et al., 2019; Nakajima et al., 2020). Therefore, we compared the equivalent width of the [O iii]5007 emission line with the UV continuum to estimate the fraction of ionising photons produced by young stars that can escape the galaxy, without being absorbed by the diffuse ISM. We expect higher escape fractions in regions with lower emission line to continuum ratio and steeper continuum. The left panel in Fig. 16 shows the ratio of the [O iii]5007 EW and UV continuum emission between 1400 and 1500 Å rest frame. Overall, we observe deep minima of the ratio to be coincident with the location of JD1-N and JD1-Sa, and steeper UV slope ($\beta_{UV}$) with respect to JD1-Sb (see Fig. 11). The right panel in Fig. 16 shows the escape fraction map of JD1, derived interpolating the predicted distribution of the escape fraction as a function of EW(H$\beta$) and the UV slope $\beta_{UV}$ for simulated dust-free galaxies between z = 7 - 9 (Zackrisson et al., 2017). While the precise value of these estimates is model dependent, they show a clear trend of increasing f_esc away from JD1-Sb, including towards JD1-Sa and JD1-N. Under the hypothesis of a highly inclined disc galaxy, as suggested by the global morphology, it is unlikely that the gravitational centre is associated with higher levels of escape fraction with respect to the galaxy edges. Finally, the absence of dust (Hashimoto et al., 2018) brings additional evidence of a merger scenario instead of the rotating disc scenario; this is because an edge-on disc – particularly with the high gas fraction we inferred from the dynamical mass estimate – may have relatively high dust column densities, even at low metallicity. Instead, there is no evidence for dust from the region with high velocity dispersion, which would represent the centre of the disc and, therefore, the highest column densities.

Figure 16: Left: Ratio between the [O iii]5007 EW and the continuum between 1400 µm and 1500 µm rest frame. The black cross marks the peak of the [O iii]5007 emission line. The three photometric unresolved clumps are highlighted by white circles. Right: Escape fraction map derived from the relations in Zackrisson et al. (2017). Black contours are arbitrary flux levels from the prism data cube collapsed over the NIRCam F200W filter.

5 Conclusions

In this work, we used integral-field spectroscopy data from JWST NIRSpec/IFS to perform a detailed spatially resolved analysis of MACS1149-JD1, a bright, lensed galaxy at z=9.11. We combined the high spectral resolution of the G395H grating with the excellent sensitivity and broad spectral coverage of the prism to draw a comprehensive picture of this system. We studied the ionised gas kinematics and excitation properties through emission-line fitting, and we performed a spaxel-by-spaxel full-spectrum fitting analysis of the prism and obtained the stellar mass and gas metallicity. The main results are summarised as follows:

• The emission line multi-Gaussian fit of the high-resolution G395H data cube revealed the peculiar ionised gas kinematics of this system. Velocity channel maps in Fig. 4 show a smooth velocity gradient between JD1-S and JD1-N, with a peak of projected velocity of $\sim$ -200 km s^-1 in correspondence of JD1-N. The spatially resolved gas velocity dispersion map shows a maximum of 90 km s^-1 perpendicular to the direction of the blue tail, $\sim$ 0.8 kpc NW from JD1-S.

• We computed the electron density and temperature from spatially resolved emission line ratios of [O ii]3726,3729 doublet and [O iii]5007/[O iii]4363 (see Fig. 7). From integrated spectra of JD1-N and JD1-S we found electron densities of 730 cm^-3 and 670 cm^-3, respectively (see Fig. 15). From direct T_e measurements, we estimated the oxygen ion abundance and computed the gas-phase metallicity distribution. From integrated spectra of the two clumps we measured 12 + log(O/H) = 8.3 and 7.71, in JD1-N and JD1-S, respectively.

• We use recently proposed emission line ratio diagrams tailored to investigate the ionisation source in high-z systems employing [O iii]5007/[O iii]4363 vs [O iii]4363/H, [O iii]5007/[O ii]3727 vs [O iii]4363/H, and [Ne iii]3869/[O ii]3727 vs [O iii]4363/H emission line ratios. Both these diagrams and the standard BPT diagram show no clear AGN signatures, thus everything is well consistent with ionisation from SF.

• We fit the SED from the low-resolution data cube and provide spatially resolved estimates of the stellar mass budget, SFR and sSFR at different epochs, and gas-phase metallicity. We confirm the presence of two distinct populations, with JD1-N being characterised by a recent SF burst and the major contribution to the stellar mass to be ascribed to an older population distributed over JD1-S (see Fig. 12), co-spatial with the enhancement of the ionised gas velocity dispersion (Fig. 3).

• From SED fitting and direct-method analysis we measured a north-south gradient of gas-phase metallicity, with higher values over JD1-N. We provided first evidence of a resolved anti-correlation between the metallicity and SFR density at such high redshift. We motivate the observed inverse correlations as due to an inflow of gas over JD1-S from JD1-N, which is diluting the ISM metal content and boosting SF.

• As for the nature of JD1, we weigh both the merger and rotating disc hypotheses. The gas kinematics are broadly consistent with a disc, albeit with a tilted kinematic axis, and a peak of velocity dispersion offset from the morphological centre. In this case, JD1 would have M_dyn = 1.2${}^{+0.5}_{-0.4}$ $\sqrt{10/\mu}$ $\times$ 10⁹$\mathrm{M_{\odot}}$ and a stellar-to-total mass fraction of 8 percent. However, this scenario does not explain the smooth metallicity gradient (Fig. 7), the negligible dust attenuation, the low gas content, and the high escape fraction near the putative centre of the disc (Fig. 16). All of these observables can instead be naturally explained if JD1-S and JD1-N were two distinct systems interacting.

Our findings represent a step forward in the comprehension of the interplay between star formation processes, chemical enrichment and merger at high-redshift. Moreover, this work highlights the major impact of JWST/NIRSpec IFS observation in constraining the ISM properties and galaxy evolution processes in the early Universe.

Data Availability

The JWST/NIRSpec data used in this work has been obtained within the NIRSpec-IFU GTO programme (program ID 1262, Observation 7) and are publicly available since June 5, 2024. The data presented in this work will be shared upon reasonable request to the corresponding author.