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

C. Marconcini, F. d'Eugenio, R. Maiolino, S. Arribas, A. Bunker, S. Carniani, S. Charlot, M. Perna, B. Rodríguez del Pino, H. Übler, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

C. Marconcini, F. d'Eugenio, R. Maiolino, S. Arribas, A. Bunker, et al. GA-NIFS: the interplay between merger, star formation, and chemical enrichment in MACS1149-JD1 at z = 9.11 with JWST/NIRSpec. Monthly Notices of the Royal Astronomical Society, 2024, 533, pp.2488-2501. 10.1093/mnras/stae1971. insu-04851524

HAL Id: insu-04851524

https://insu.hal.science/insu-04851524v1

Submitted on 21 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

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

C. Marconcini<sup>1,2,3*</sup> F. D'Eugenio<sup>1,4</sup> R. Maiolino,<sup>1,4,5</sup> S. Arribas,<sup>6</sup> A. Bunker,<sup>7</sup> S. Carniani<sup>8</sup>

S. Charlot<sup>9</sup> M. Perna,<sup>6</sup> B. Rodríguez Del Pino<sup>6</sup>,<sup>6</sup> H. Übler<sup>1,4</sup> C. J. Willott,<sup>10</sup> T. Böker,<sup>11</sup> G. Cresci,<sup>3</sup>
M. Curti,<sup>12</sup> G. C. Jones,<sup>7</sup> I. Lamperti<sup>2,3</sup> E. Parlanti<sup>8</sup> and G. Venturi<sup>1,8</sup>

Accepted 2024 August 13. Received 2024 August 13; in original form 2024 July 12

## ABSTRACT

We present JWST/NIRSpec integral-field spectroscopy observations of the z ~ 9.11 lensed galaxy MACS1149-JD1, as part of the GA-NIFS programme. The data were obtained with both the G395H grating (R ~ 2700) and the prism (R ~ 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 ionized-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_e ~ 1.0 × 10³ cm⁻<fnref n="3" />), temperature (T_e ~ 1.6 × 10⁴ 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σ evidence of a spatially resolved anticorrelation 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 β, ionized gas velocity dispersion, and gas-phase metallicity.

**Key words:** galaxies: abundances – galaxies: high-redshift – galaxies: ISM – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

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 ~ 10 (Curti et al. 2023; Curtis-Lake et al. 2023; Harikane et al. 2023; Isobe et al. 2023; Larson et al. 2023; Robertson et al. 2023; Hsiao et al. 2023a; Abdurro’uf et al. 2024; Carniani et al. 2024; Sanders et al. 2024; Vikaeus 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; Claeysens 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; Förster Schreiber et al. 2011; Schreiber et al. 2015; Ferreira et al. 2020; Le Fèvre et al. 2020; Meng & Gnedin 2021). 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

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Kavli Institute for Cosmology, University of Cambridge, Madingley Road, CB3 0HE Cambridge, UK</footnote>

<footnote n="2"><sup>2</sup>Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Firenze, Italy</footnote>

<footnote n="3"><sup>3</sup>INAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="4"><sup>4</sup>Cavendish Laboratory, University of Cambridge, 19 J. J. Thomson Ave., CB3 0HE Cambridge, UK</footnote>

<footnote n="5"><sup>5</sup>Department of Physics and Astronomy, University College London, Gower Street, WC1E 6BT London, UK</footnote>

<footnote n="6"><sup>6</sup>Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850 Madrid, Spain</footnote>

<footnote n="7"><sup>7</sup>Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, OX13RH Oxford, UK</footnote>

<footnote n="8"><sup>8</sup>Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="9"><sup>9</sup>Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="10"><sup>10</sup>NRC Herzberg, 5071 West Saanich Rd, Victoria, V9E 2E7 BC, Canada</footnote>

<footnote n="11"><sup>11</sup>European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, 21218 MD, USA</footnote>

<footnote n="12"><sup>12</sup>European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany</footnote>

<footnote>*E-mail: cosimo.marconcini@unifi.it</footnote>
of galaxy evolution, shaping galaxy morphology and kinematics (Kaviraj 2014; Zanella et al. 2015; Jackson et al. 2022; Boylan-Kolchin 2023; Huško, Lacey & Baugh 2023; Zhang et al. 2023). Additionally, it is crucial to characterize other drivers of galaxy evolution such as the SF activity, the gas-phase distribution, and the multiphase 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 ionized 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 (μ = 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 ~ 15) and of a younger component. From SED fitting, Laporte et al. (2021) estimated a best-fitting stellar population of age 512 Myr. Moreover, they inferred past star-formation histories (SFH) and showed that more than 93 per cent of the observed stellar mass had already formed at z ~ 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 disc. Consistent with Hashimoto et al. (2018), Tokuoka et al. (2022) found that a ~300 Myr old stellar population of 2.3 × 10⁹ M☉ 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 α map, with a star formation rate (SFR) of 3.2–5.3 M☉ yr⁻<fnref n="1" />. Moreover, they identified two main clumps with indication of a recent stellar burst (≥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 characterize both the ionized gas kinematics and ISM properties such as the electron density and temperature (Kewley, Nicholls & Sutherland 2019; Abdurro'uf et al. 2024; Mazzolari et al. 2024; Morishita 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 ~ 2700) and low-resolution prism (R ~ 100) observations. This paper is organized as follows. In Section 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 Section 4, we discuss the impact of our results in the broad context of galaxy assembling processes at high-redshift. Finally, in Section 5, we summarize the main results. In this paper, we adopt a ΛCDM cosmology: H₀: 67.4 km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" />, Ω_m = 0.315, and Ω_Λ = 0.685.

# 2 OBSERVATIONS AND DATA REDUCTION

## 2.1 NIRSpec IFU observations

MACS1149-JD1 (hereafter JD1, RA 11h 49m 33.580s Dec. +22° 24′45.70″, J2000) was observed as part of the Galaxy Assembly with NIRSpec Integral Field Spectroscopy (GA-NIFS.<fnref n="1" />; 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 2023 June 5, using the NIRSpec IFS (Böker et al. 2022; Jakobsen et al. 2022), covering a 3.1 arcsec × 3.2 arcsec region, with a native spaxel size of 0.1 arcsec (Böker et al. 2022; Rigby et al. 2023). The observations consist of two configurations: high-spectral resolution with the G395H grating (R ~ 2700, covering 2.87–5.14 µm) and low-spectral resolution with the prism (R ~ 100, covering 0.6–5.3 µm). Both observations used a medium (i.e. ~0.5 arcsec) 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 18 207 s (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 s (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 multiband properties of this system as discussed in Section 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 μ 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 μ = 10 ±0.7 and corrected for the lensing magnification all the appropriate properties, thus assuming that all flux-dependent quantities are in units of (10/μ) (Hashimoto et al. 2018; Stiavelli et al. 2023; Álvarez-Márquez et al. 2024; Bradač et al. 2024).

## 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 arcsec, we adopted the *drizzle* method, with particular attention to use an official patch

<footnote n="1">¹GA-NIFS website: https://ga-nifs.github.io</footnote>
<figure><img src="image_4.png" /><figcaption>**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. 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.</figcaption></figure>

to account for a known issue affecting the calibration pipeline.<fnref n="2" /> 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 data cube, 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 arcsec) 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. ~450 pc) and 0.15" (i.e. ~670 pc), respectively. The two apertures correspond to the putative location of the two

<footnote n="2">²https://github.com/spacetelescope/jwst/pull/7306</footnote>
<figure><img src="image_5.png" /><figcaption>**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-fitting 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 ζ + He I λ3888 (left), and the H β and [O III] complex (right). Shaded grey areas are masked spectral channels.</figcaption></figure>

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 I λ3889], 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.

## 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-fitting 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-fitting 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

<figure><img src="image_6.png" /><figcaption>**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-fitting 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-fitting model are in black and red, respectively. The two Gaussian components are represented in blue and light blue solid lines.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**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⁻¹.</figcaption></figure>

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-fitting 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$ minimization, likely indicating the presence of outflow or merger. Similarly, Álvarez-Márquez et al. (2024) analysed MIRI/MRS observations tracing the H α kinematics and measured a difference in velocity dispersion of ≈60 km s⁻¹ among the two clumps, supporting the hypothesis of outflowing gas powered by the UV-bright clump JD1-S (Tokuoka et al. 2022, see also section 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 ~90 km s⁻¹. The zero-order moment map in Fig. 3 shows the [O III] λ5007 emission of ionized 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⁻¹.

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 ≈10° 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 Section 4.3.

From the spatially resolved [O III] λ5007 emission-line fit, we measured an average velocity dispersion of 62 ± 7 km s⁻¹, 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_{-52}^{+157} \sqrt{10/\mu}$ pc and Sérsic index of 1 ± 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.4}^{+0.5} \sqrt{10/\mu} \times 10^9 M_{\odot}$. This value is consistent with previous H α-based estimates of 2.4 ± 0.5 × 10⁹ M$_{\odot}$ from Álvarez-Márquez et al. (2024) and with resolved [O III] λ88 µm-based estimates of

<figure><img src="image_8.png" /><figcaption>**Figure 5.** AGN-diagnostic diagram of Mazzolari et al. (2024) for the excitation source in JD1 based on the [O III] λ4363/H γ versus [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).</figcaption></figure>

$0.7–3.7 \times 10^9 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–2.6 \sqrt{10/\mu} \times 10^9 M_{\odot}$, which does not change our conclusions.

## 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 ionization at high-z and low metallicity by using the standard BPT diagram (Feltre, Charlot & Gutkin 2016; Nakajima & Maiolino 2022; Kocevski et al. 2023; Maiolino et al. 2023; Scholtz et al. 2023; Übler et al. 2023). We therefore explored the possible source of ionization 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 γ versus [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 versus [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 ionization 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} \text{ erg s}^{-1} \text{ cm}^{-2}$, a $3\sigma$ upper limit flux on [N II] of $2.1 \times 10^{-18} \text{ erg s}^{-1} \text{ cm}^{-2}$ from
<figure><img src="image_9.png" /><figcaption>**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 ≤ 3 are masked.</figcaption></figure>

Á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 ionization nature for JD1, with the location of this system being located on the demarcation line between SF and AGN ionization. Therefore, from the results of both diagnostics we cannot conclusively state the nature of the main ionization source in JD1.

In order to trace regions of high ionization, we also computed the [O III] λ5007/H β line ratio map shown in Fig. 6. Due to the difference in ionization 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 characterized by high excitation, with [O III] λ5007/H β ~ 8, at variance with JD1-S which is characterized by an average [O III] λ5007/H β ~ 6.

## 3.4 Electron density and temperature

The high-spectral resolution G395H data include 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 \ge 9$. In particular, first, we used the PYNEBPYTHON package (Luridiana, Morisset & Shaw 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 per cent) and an electron density of $600 \text{ 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_2$), due to the absence of temperature diagnostics, we used the relation $T_2 = -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_2$ 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 Section 3.6, together with an estimate via prism continuum fitting. The average electron temperature and density of JD1-S are $1.4_{-0.7}^{+0.6} \times 10^4$ K and $670_{-240}^{+100} \text{ 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^4$ K, we estimated an electron density from integrated spectrum of $730 \text{ 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 Abdurro'uf et al. 2024). 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 Section 4 for a more detailed discussion).

## 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, Gunn & White 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, parametrized 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, 30, 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
<figure><img src="image_10.png" /><figcaption>**Figure 7.** ISM physical properties of JD1. From left to right: Electron density, temperature, gas-phase metallicity derived with the direct method (Te method) and gas-phase metallicity derived with PROSPECTOR performing SED fitting (see Sections 3.4 and 3.5 for details). Spaxels at S/N ≤ 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.</figcaption></figure>

<table><caption>**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.</caption><thead><tr><th></th><th>Parameter<br>(1)</th><th>Free<br>(2)</th><th>Description<br>(3)</th><th>Prior<br>(4)</th><th>Posterior<sub>JD1</sub><br>(5)</th></tr></thead><tbody><tr><td rowspan="8">Free parameters</td><td>log M<sub>*</sub>[M<sub>☉</sub>]</td><td>Y</td><td>total stellar mass formed</td><td>U(7, 13)</td><td>8.47<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>log Z[Z<sub>☉</sub>]</td><td>Y</td><td>stellar metallicity</td><td>U(-2, 0.19)</td><td>-1.62<sup>+0.10</sup><sub>-0.10</sub></td></tr><tr><td>log SFR ratios</td><td>Y</td><td>ratio of the log SFR between adjacent bins of the SFH</td><td>T(0, 0.3, 2)</td><td>—</td></tr><tr><td>τ<sub>V</sub></td><td>Y</td><td>optical depth of the diffuse dust</td><td>G(0.3, 1; 0, 2)</td><td>0.63<sup>+0.08</sup><sub>-0.10</sub></td></tr><tr><td>μ</td><td>Y</td><td>ratio between the optical depth of the birth clouds and τ<sub>V</sub></td><td>U(-1.0, 0.4)</td><td>0.02<sup>+0.05</sup><sub>-0.02</sub></td></tr><tr><td>σ<sub>gas</sub> [km s<sup>-1</sup>]</td><td>Y</td><td>intrinsic velocity dispersion of the star-forming gas<sup>‡</sup></td><td>U(0, 300)</td><td>64<sup>+44</sup><sub>-38</sub></td></tr><tr><td>log Z<sub>gas</sub>[Z<sub>☉</sub>]</td><td>Y</td><td>metallicity of the star-forming gas</td><td>U(-2, 0.19)</td><td>-0.79<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>log U</td><td>Y</td><td>ionization parameter of the star-forming gas</td><td>U(-4, -1)</td><td>-1.02<sup>+0.01</sup><sub>-0.02</sub></td></tr><tr><td rowspan="2">Other</td><td>log SFR<sub>10</sub>[M<sub>☉</sub> yr<sup>-1</sup>]</td><td>N</td><td>star-formation rate averaged over the last 10 Myr</td><td>—</td><td>1.46<sup>+0.04</sup><sub>-0.05</sub></td></tr><tr><td>log SFR<sub>100</sub>[M<sub>☉</sub> yr<sup>-1</sup>]</td><td>N</td><td>star-formation rate averaged over the last 100 Myr</td><td>—</td><td>0.47<sup>+0.04</sup><sub>-0.05</sub></td></tr></tbody></table>

Note. (1) Parameter name and units (where applicable). (2) Only parameters marked with ‘Y’ are optimized 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*, *τ*<sub>V</sub>, and *μ*, see Tacchella et al. (2022, their equations 4 and 5). (4) Parameter prior probability distribution; *N*(*μ*, *σ*) is the normal distribution with mean *μ* and dispersion *σ*; *U*(*a*, *b*) is the uniform distribution between *a* and *b*; *T*(*μ*, *σ*, *ν*) is the Student’s *t* distribution with mean *μ*, dispersion *σ* and *ν* degrees of freedom; *G*(*μ*, *σ*, *a*, *b*) is the normal distribution with mean *μ* and dispersion *σ*, truncated between *a* and *b*. (5) Median and 16th–84th percentile range of the marginalized posterior distribution for JD1; 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 (~500 km s<sup>-1</sup> at the wavelength of [O III] λ5007); indeed, the posterior probability distributions for JD1 is both consistent with 0 km s<sup>-1</sup>.

of the SFR surface density within the last 10 (SFR<sub>10</sub>) and 100 Myr (SFR<sub>100</sub>).

From the integrated spectrum of JD1, we found a total stellar mass budget of log(M<sub>*</sub>/M<sub>☉</sub>)<sub>JD1</sub> = 7.47<sup>+0.05</sup><sub>-0.05</sub> log(10/μ), 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(M<sub>*</sub>/M<sub>☉</sub>) = 7.88 log(10/μ)). 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 Section 3.2) and assuming a rotating disc, we estimated a (de-lensed) maximal dynamical mass of M<sub>dyn</sub> = 1.2<sup>+0.5</sup><sub>-0.4</sub>√10/μ × 10<sup>9</sup> M<sub>☉</sub>. 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-virialized motions (we discuss this finding in Section 4) or large amounts of dark matter (de Graaff et al. 2024) or a combination of the aforementioned effects.

From the SED fitting of the integrated spectrum we derived a SFR<sub>10</sub> = 2.8<sup>+0.1</sup><sub>-0.1</sub> M<sub>☉</sub> yr<sup>-1</sup> and SFR<sub>100</sub> = 0.29<sup>+0.1</sup><sub>-0.1</sub> M<sub>☉</sub> yr<sup>-1</sup>, 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 SFR<sub>10</sub> = 9.6<sup>+1.9</sup><sub>-1.7</sub> Gyr<sup>-1</sup> and SFR<sub>100</sub> = 0.97<sup>+0.2</sup><sub>-0.2</sub> Gyr<sup>-1</sup>. The sSFR maps show consistent values between JD1-N and JD1-S in the last 100 Myr. Within the last
<figure><img src="image_11.png" /><figcaption>Figure 8. JD1 properties derived from SED fitting of the low-resolution data cube with prospector (see Section 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 and 100 Myr. Black contours are arbitrary flux levels from the prism data cube collapsed over the NIRCam F200W filter.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>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.</figcaption></figure>

10 Myr instead, we observe higher sSFR over JD1-S compared to JD1-N. None the less, comparing the sSFR within the last 10 and 100 Myr for JD1-N, we conclude that the majority of its mass formed very recently (Stiavelli et al. 2023; Álvarez-Márquez et al. 2024; 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 the 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 \text{ erg s}^{-1} \text{ cm}^{-2}$, we computed the SFR as:

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

where $f_{H\alpha/H\beta} = 2.85$ for $T_e = 1.4 \times 10^4$ K (see Section 3.4). Correcting for the lensing factor, we obtain SFR(H β) = $2.5 \pm 0.3$ ($10/\mu$) $M_{\odot} \text{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

<figure><img src="image_13.png" /><figcaption>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.</figcaption></figure>

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 Section 4.3, JD1-N may be characterized by an extreme escape fraction of ionizing photons. The effect of an elevated escape fraction over JD1-N translates in a loss of ionizing 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 β$_{UV}$ slope, measured between 1200 to 1500 Å rest-frame (where β$_{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 β$_{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
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Continuum slope map derived from prism data cube fit in the spectral rest-frame range 1200 to 1500 Å (see Section 3.5). Black contours are arbitrary flux levels from the prism data cube collapsed over the NIRCam F200W filter.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Star formation histories of single spaxels (purple) and integrated JD1-N (green) and JD1-S (maroon) clumps. Error bars for the main clumps represent the 84^th^ and 16^th^ percentiles.</figcaption></figure>

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 α 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 ± 1.1 × 10⁻¹⁸ erg s⁻¹ cm⁻², 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 ~ 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 characterized 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).

## 3.6 Gas-phase metallicity

In Section 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 ionization 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 ± 0.1 and 7.71 ± 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 and 7.85±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 Section 3.5). The metallicity derived via prism fitting of the integrated spectrum of JD1 is log(O/H) = 7.89±0.05, consistent with the value of 7.90±0.04 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<sub>☉</sub>-Z<sub>gas</sub>-SFR that does not appear to evolve with redshift (Tremonti et al. 2004; Mannucci et al. 2010; Cresci, Mannucci & Curti 2019; Maiolino & Mannucci 2019). Recently,
<figure><img src="image_16.png" /><figcaption>**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 ≤ 3 are masked.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**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 Section 3.5 and Figs 7 and 8).</figcaption></figure>

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$, and anticorrelates with $\Sigma_{SFR}$. We employed the resolved maps of SFR$_{100}$ density and the gas-phase metallicity derived from PROSPECTOR to test such observed anticorrelation at $z \ge 9$. Fig. 14 shows the observed resolved anticorrelation between the gas phase metallicity and the SFR$_{100}$ 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$_{100}$ 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 Section 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 and 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).

# 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 summarize 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 ~ 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 ~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] $\lambda$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 ionized 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] $\lambda$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.

Interestingly, the right panel in Fig. 7 shows a peak of the gas-phase metallicity over JD1-N (12 + log (O/H) ~ 7.95), consistent with the integrated value obtained from the direct method. The higher metallicity and [O III] $\lambda$5007/H β line ratio support the idea of an 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] $\lambda$5007/H β can be produced by high levels of sSFR (see Fig. 9), which in turn enhances the ionization parameter and thus indicates high-SF efficiency (Dickey et al. 2016; Kewley et al. 2016; Sanders et al. 2016a; 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 and 9). This scenario is also supported by VLT/X-Shooter
<figure><img src="image_18.png" /><figcaption>**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 $\alpha(1+z)$ and $\alpha(1+z)^2$ 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. 2012a, b; 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 (Lamperti et al. 2024; Rodríguez Del Pino et al. 2024).</figcaption></figure>

observations of a diffuse Ly $\alpha$ bubble, supposedly inflated by UV emission of the older population in the JD1-S (Hashimoto et al. 2018).

## 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 (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 (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é, Finlator & Oppenheimer 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 Davies et al. 2021), 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, Ouchi & Harikane 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 (1+z)^{1.4-2}$ relation. On the other hand, recent works motivate the $n_e \propto (1+z)^1$ 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 anticorrelation 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, Kewley & Richard 2013), we get the dotted relation: $n_e \propto (1+z)^1$.

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

ALMA [O III] $\lambda$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] $\lambda$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; Bsemi Stiavelli et al. 2023; Bradač et al. 2024). In
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Left: Ratio between the [O III] λ5007 EW and the continuum between 1400 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.</figcaption></figure>

particular, Bradač et al. (2024) identified three unresolved star-
forming clumps with an underlying smooth galaxy component. From
our analysis, the high-resolution data do 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 ≈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
ionizing 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 (β<sub>UV</sub>) 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 β) and the UV slope β<sub>UV</sub> 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<sub>esc</sub> 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.

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 ionized 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 summarized as follows:

(i) The emission line multi-Gaussian fit of the high-resolution G395H data cube revealed the peculiar ionized 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 ~−200 km s⁻¹ in correspondence of JD1-N. The spatially resolved gas velocity dispersion map shows a maximum of 90 km s⁻¹ perpendicular to the direction of the blue tail, ~0.8 kpc NW from JD1-S.

(ii) We computed the electron density and temperature from spatially resolved emission line ratios of [O II] λλ3726,3729 doublet and [O III] λ5007/[O II] λ4363 (see Fig. 7). From integrated spectra of JD1-N and JD1-S, we found electron densities of 730 and 670 cm⁻³, 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.

(iii) We use recently proposed emission line ratio diagrams tailored to investigate the ionization source in high-z systems employing [O III] λ5007/[O III] λ4363 versus [O III] λ4363/H γ, [O III] λ5007/[O II] λ3727 versus [O III] λ4363/H γ, and Ne III λ3869/[O II] λ3727 versus [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 ionization from SF.

(iv) 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 characterized 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 ionized gas velocity dispersion (Fig. 3).

(v) 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 anticorrelation 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.

(vi) 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.4}^{+0.5} \sqrt{10/\mu} \times 10^9 M_{\odot}$ and a stellar-to-total mass fraction of 8 per cent. 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.

## ACKNOWLEDGEMENTS

CM, RM, and FDE acknowledge support by the Science and Technology Facilities Council (STFC), by the ERC Advanced Grant 695671 'QUENCH', and by the UKRI Frontier Research Grant RISEandFALL. CM also acknowledges the support of the INAF Large Grant 2022 'The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities' and of the grant PRIN-MUR 2020ACSP5K.002 financed by European Union-Next Generation EU. RM is further supported by a research professorship from the Royal Society, SA, BRdP, and MP acknowledge support from the research project PID2021-127718NB-I00 of the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI). HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship, GCJ and AJB acknowledge funding from the 'FirstGalaxies' Advanced Grant from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 789056). SC, EP, and GV acknowledge support by European Union's HE ERC Starting grant no. 101040227-WINGS. GC acknowledges the support of the INAF Large Grant 2022 'The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities'.

## DATA AVAILABILITY

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

## REFERENCES

* Abdurro'uf et al., 2024, preprint (arXiv:2404.16201)

* Álvarez-Márquez J. et al., 2024, A&A, 686, A85

* Arribas S. et al., 2023, A&A, 688, 28

* Baker W. M. et al., 2023, MNRAS, 519, 1149,

* Berg D. A. et al., 2012, ApJ, 754, 98,

* Berg D. A. et al., 2022, ApJS, 261, 31,

* Böker T. et al., 2022, A&A, 661, A82

* Boylan-Kolchin M., 2023, Nat. Astron., 7, 731,

* Bradač M. et al., 2024, ApJ, 961, L21,

* Byler N., Dalcanton J. J., Conroy C., Johnson B. D., 2017, ApJ, 840, 44

* Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storch-Bergmann T., 2000, ApJ, 533, 682

* Cappellari M., 2023, MNRAS, 526, 3273

* Cappellari M., Copin Y., 2003, MNRAS, 342, 345,

* Cappellari M., Emsellem E., 2004, PASP, 116, 138,

* Cappellari M. et al., 2006, MNRAS, 366, 1126,

* Cappellari M. et al., 2013, MNRAS, 432, 1709,

* Carniani S. et al., 2024, preprint (arXiv:2405.18485)

* Charlot S., Fall S. M., 2000, ApJ, 539, 718,

* Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102

* Christensen L. et al., 2012a, MNRAS, 427, 1953,

* Christensen L. et al., 2012b, MNRAS, 427, 1973,

* Claeysens A., Adamo A., Richard J., Mahler G., Messa M., Dessauges-Zavadsky M., 2023, MNRAS, 520, 2180

* Conroy C., Gunn J. E., 2010, ApJ, 712, 833,

* Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486,

* Costantin L. et al., 2023, ApJ, 946, 71,

* Cresci G., Mannucci F., Maiolino R., Marconi A., Gnerucci A., Magrini L., 2010, Nature, 467, 811,

* Cresci G., Mannucci F., Curti M., 2019, A&A, 627, A42

* Curti M. et al., 2023, MNRAS, 518, 425

* Curtis-Lake E. et al., 2023, Nat. Astron., 7, 622

* D'Eugenio F. et al., 2023, preprint (arXiv:2308.06317)

* Davé R., Finlator K., Oppenheimer B. D., 2011, MNRAS, 416, 1354

* Davies R. L. et al., 2021, ApJ, 909, 78

* de Graaff A. et al., 2024, A&A, 684, A87

* Dekel A. et al., 2009, Nature, 457, 451

* Delgado-Serrano R., Hammer F., Yang Y. B., Puech M., Flores H., Rodrigues M., 2010, A&A, 509, A78

* Dickey C. M. et al., 2016, ApJ, 828, L11

* Falcón-Barroso J., Sánchez-Blázquez P., Vazdekis A., Ricciardelli E., Cardiel N., Cenarro A. J., Gorgas J., Peletier R. F., 2011, A&A, 532, A95

* Feltre A., Charlot S., Gutkin J., 2016, MNRAS, 456, 3354

* Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761

* Ferreira L., Conzelice C. J., Duncan K., Cheng T.-Y., Griffiths A., Whitney A., 2020, ApJ, 895, 115

* Fletcher T. J., Tang M., Robertson B. E., Nakajima K., Ellis R. S., Stark D. P., Inoue A., 2019, ApJ, 878, 87

* Förster Schreiber N. M. et al., 2011, ApJ, 739, 45

* Garcia F. A. B., Ricotti M., Sugimura K., Park J., 2023, MNRAS, 522, 2495

* Gburek T. et al., 2019, ApJ, 887, 168

* Giménez-Arteaga C. et al., 2023, ApJ, 948, 126

* Harikane Y. et al., 2023, ApJS, 265, 5

* Hashimoto T. et al., 2018, Nature, 557, 392

* Hsiao T. Y.-Y. et al., 2023a, preprint (arXiv:2305.03042)

* Hsiao T. Y.-Y. et al., 2023b, ApJ, 949, L34

* Huertas-Company M. et al., 2024, A&A, 685, A48

* Huško F., Lacey C. G., Baugh C. M., 2023, MNRAS, 518, 5323

* Isobe Y., Ouchi M., Nakajima K., Harikane Y., Ono Y., Xu Y., Zhang Y., Umeda H., 2023, ApJ, 956, 139

* Izotov Y. I., Stasińska G., Meynet G., Guseva N. G., Thuan T. X., 2006, A&A, 448, 955

* Izotov Y. I., Worseck G., Schaerer D., Guseva N. G., Thuan T. X., Fricke Verhamme A., Orlitová I., 2018, MNRAS, 478, 4851

* Jackson R. A., Kaviraj S., Martin G., Devriendt J. E. G., Noakes-Kettel E. A., Silk J., Ogle P., Dubois Y., 2022, MNRAS, 511, 607

* Jakobsen P. et al., 2022, A&A, 661, A80

* Johnson B. D., Leja J., Conroy C., Speagle J. S., 2021, ApJS, 254, 22

* Kaasinen M., Bian F., Groves B., Kewley L. J., Gupta A., 2017, MNRAS, 465, 3220

* Kartaltepe J. S. et al., 2023, ApJ, 946, L15

* Kashino D. et al., 2017, ApJ, 835, 88

* Kaviraj S., 2014, MNRAS, 440, 2944

* Kewley L. J. et al., 2016, ApJ, 819, 100

* Kewley L. J., Nicholls D. C., Sutherland R. S., 2019, ARA&A, 57, 511

* Kocovski D. D. et al., 2023, ApJ, 954, L4

* Kriek M., Conroy C., 2013, ApJ, 775, L16

* Lamperti I. et al., 2024, preprint (arXiv:2406.10348)

* Laporte N., Meyer R. A., Ellis R. S., Robertson B. E., Chisholm J., Roberts-Borsani G. W., 2021, MNRAS, 505, 3336
Larson R. L. et al., 2023, ApJ, 953, L29

Le Fèvre O. et al., 2020, A&A, 643, A1

Leja J., Carnall A. C., Johnson B. D., Conroy C., Speagle J. S., 2019, ApJ, 876, 3

Luridiana V., Morisset C., Shaw R. A., 2012, in Proc. IAU Symp. 283, Planetary Nebulae: An Eye to the Future, p. 422

Luridiana V., Morisset C., Shaw R. A., 2015, A&A, 573, A42

Maiolino R., Mannucci F., 2019, A&A Rev., 27, 3

Maiolino R. et al., 2023, preprint (arXiv:2308.01230)

Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010, MNRAS, 408, 2115

Marasco A. et al., 2020, A&A, 644, A15

Mazzolari G. et al., 2024, preprint (arXiv:2404.10811)

Meng X., Gnedin O. Y., 2021, MNRAS, 502, 1433

Meštrić U. et al., 2022, MNRAS, 516, 3532

Morishita T. et al., 2024, ApJ, 971, 14,

Mowla L. et al., 2022, ApJ, 937, L35

Nakajima K., Maiolino R., 2022, MNRAS, 513, 5134

Nakajima K., Ellis R. S., Robertson B. E., Tang M., Stark D. P., 2020, ApJ, 889, 161

Noll S. et al., 2009, A&A, 499, 69

Ono Y. et al., 2023, ApJ, 951, 72

Ormerod K. et al., 2024, MNRAS, 527, 6110

Osterbrock D. E., Ferland G. J., 2006, Astrophysics of gaseous nebulae and active galactic nuclei, Sausalito, CA: University Science Books

Pérez-González P. G., Gil de Paz A., Zamorano J., Gallego J., Alonso-Herrero A., Aragón-Salamanca A., 2003, MNRAS, 338, 508

Pérez-González P. G. et al., 2008, ApJ, 675, 234

Pérez-González P. G. et al., 2023, ApJ, 946, L16

Perna M. et al., 2023, A&A, 679, A89

Rauscher B. J. et al., 2017, PASP, 129, 105003

Reddy N. A., Topping M. W., Sanders R. L., Shapley A. E., Brammer G., 2023, ApJ, 952, 167

Rigby J. et al., 2023, PASP, 135, 048001

Robertson B. E. et al., 2023, Nat. Astron., 7, 611

Rodríguez Del Pino B. et al., 2024, A&A, 684, A187

Sanders R. L. et al., 2016a, ApJ, 816, 23

Sanders R. L. et al., 2016b, ApJ, 825, L23

Sanders R. L., Shapley A. E., Topping M. W., Reddy N. A., Brammer G. B., 2024, ApJ, 962, 24

Scholtz J. et al., 2023, preprint (arXiv:2311.18731)

Schreiber C. et al., 2015, A&A, 575, A74

Shibuya T., Ouchi M., Harikane Y., 2015, ApJS, 219, 15

Steidel C. C. et al., 2014, ApJ, 795, 165

Stiavelli M. et al., 2023, ApJ, 957, L18

Suess K. A. et al., 2022, ApJ, 937, L33

Swinbank A. M. et al., 2019, MNRAS, 487, 381

Tacchella S. et al., 2022, ApJ, 926, 134

Tacchella S. et al., 2023, MNRAS, 522, 6236

Tokuoka T. et al., 2022, ApJ, 933, L19

Tozzi G. et al., 2021, aap, 648, A99

Tremonti C. A. et al., 2004, ApJ, 613, 898

Tripodi R. et al., 2024, preprint (arXiv:2403.08431)

Übler H. et al., 2023, A&A, 677, A145

van der Wel A. et al., 2022, ApJ, 936, 9

Vanzella E. et al., 2017, ApJ, 842, 47

Vanzella E. et al., 2022, ApJ, 940, L53

Vikaeus A. et al., 2024, MNRAS, 529, 1299

Wisnioski E. et al., 2018, ApJ, 855, 97

Yuan T. T., Kewley L. J., Richard J., 2013, ApJ, 763, 9

Zackrisson E. et al., 2017, ApJ, 836, 78

Zanella A. et al., 2015, Nature, 521, 54

Zhang H., Behroozi P., Volonteri M., Silk J., Fan X., Hopkins P. F., Yang J., Aird J., 2023, MNRAS, 518, 2123

Zheng W. et al., 2012, Nature, 489, 406

This paper has been typeset from a TeX/LA TeX file prepared by the author.