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

# GA-NIFS: Multi-phase outflows in a star-forming galaxy
at $z \sim 5.5$

Eleonora Parlanti, Stefano Carniani, Giacomo Venturi, Rodrigo Herrera-Camus, Santiago Arribas, Andrew J Bunker, Stephane Charlot,

Francesco d'Eugenio, Roberto Maiolino, Michele Perna, et al.

## ► To cite this version:

Eleonora Parlanti, Stefano Carniani, Giacomo Venturi, Rodrigo Herrera-Camus, Santiago Arribas, et al.. GA-NIFS: Multi-phase outflows in a star-forming galaxy at $z \sim 5.5$. Astronomy & Astrophysics - A&A, 2025, 695, pp.A6. 10.1051/0004-6361/202451692. hal-04679965

HAL Id: hal-04679965

https://hal.science/hal-04679965v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>
# GA-NIFS: Multiphase analysis of a star-forming galaxy at $z \sim 5.5$

Eleonora Parlanti<sup>1,★</sup>●, Stefano Carniani<sup>1,●</sup>, Giacomo Venturi<sup>1,●</sup>, Rodrigo Herrera-Camus<sup>2</sup>, Santiago Arribas<sup>3,●</sup>, Andrew J. Bunker<sup>4</sup>, Stéphane Charlot<sup>5,●</sup>, Francesco D'Eugenio<sup>6,7</sup>, Roberto Maiolino<sup>6,7,8</sup>, Michele Perna<sup>3,●</sup>, Hannah Übler<sup>6,7</sup>, Torsten Böker<sup>9</sup>, Giovanni Cresci<sup>10,●</sup>, Mirko Curti<sup>11</sup>, Gareth C. Jones<sup>4,●</sup>, Isabella Lamperti<sup>12,10,3,●</sup>, Pablo G. Pérez-González<sup>3</sup>, Bruno Rodríguez Del Pino<sup>3</sup>, and Sandra Zamora<sup>1,●</sup>

1. Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy

2. Departamento de Astronomía, Universidad de Concepción, Barrio Universitario, Concepción, Chile

3. Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz E-28850, Madrid, Spain

4. University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX13RH, United Kingdom

5. Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France

6. Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

7. Cavendish Laboratory – Astrophysics Group, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK

8. Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

9. European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, MD 21218, USA

10. INAF – Osservatorio Astrofisco di Arcetri, largo E. Fermi 5, 50127 Firenze, Italy

11. European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany

12. Dipartimento di Fisica e Astronomia, Università di Firenze, Via G. Sansone 1, 50019 Sesto F.no (Firenze), Italy

Received 29 July 2024 / Accepted 1 February 2025

## ABSTRACT

In this study, we present a detailed multiphase analysis of HZ4, a main-sequence star-forming galaxy at $z \sim 5.5$, known for being a turbulent rotating disk and having a detection of a [CII] outflow in the ALMA observations. We exploited JWST/NIRSpec observations in the integral field spectroscopy mode with low- and high-spectral resolution which allow us, for the first time, to spatially resolve the rest-frame UV and optical emission of the galaxy to investigate the galaxy properties. In particular, the high-resolution dataset allowed us to study the kinematics of the ionized gas phase, and the conditions of the interstellar medium, such as the excitation mechanism, dust attenuation, and metallicity. The lower spectral-resolution observations allowed us to study the continuum emission and infer the stellar populations' ages and properties. Our findings suggest that HZ4 is a galaxy merger rather than a rotating disk as previously inferred from lower-resolution [CII] data. The merger is associated with an extended broad, blueshifted emission, potentially indicative of an outflow originating from a region of intense star formation and extending up to 4 kpc. In light of these new observations, we reanalyzed the ALMA data to compare the multiphase gas properties. If we interpret the broad components seen in [CII] and [OIII]$\lambda$5007 Å as outflows, the neutral and ionized components are co-spatial, and the mass loading factor of the ionized phase is significantly lower than that of the neutral phase, aligning with trends observed in multiphase systems at lower redshifts. Nonetheless, additional observations and larger statistical samples are essential to determine the role of mergers and outflows in the early Universe and to clarify the origin of the broad emission components observed in this system.

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

## 1. Introduction

Several theoretical models suggest that galactic outflows driven by star formation (SF) and active galactic nuclei (AGNs) are crucial to explain the lack of galaxies in both the high- and low-mass ends of the galaxy mass function (e.g., Fabian 2012; Naab & Ostriker 2017; Somerville & Davé 2015), the inefficiency of galaxies in turning baryons into stars (Dekel & Silk 1986; Behroozi et al. 2019), the metal enrichment of the circumgalactic medium (CGM) and the intergalactic medium (IGM) (Oppenheimer & Davé 2006), and the shape of the mass-metallicity relation (Chisholm et al. 2018). In particular, low-mass (<10 $M_{\odot}$) galaxies are believed to self-regulate the buildup of their stellar mass through the action of SF-driven outflows (Hopkins et al. 2012). The energy and momentum injected by young massive stars in the shape of stellar winds, supernova (SN) explosions, and radiation pressure can provide feed-

back, by heating and expelling the gas from the galaxy (e.g., Somerville & Davé 2015).

Local observations of low-mass galaxies have revealed that they also host accreting supermassive black holes at their centers (e.g., Sartori et al. 2015; Reines et al. 2020), indicating that SF activity can be regulated by AGN-driven outflows. Theoretical works and simulations have indeed pointed out that AGNs can boost the outflow temperature and velocity by two orders of magnitude (Silk 2017; Dashyan et al. 2018; Koudmani et al. 2019), which means that AGN-driven outflows may affect the CGM around low-mass galaxies by polluting it with metals and heating it. Koudmani et al. (2021) have also shown that AGN feedback in low-mass galaxies may be more relevant at $z > 2$, when AGN are expected to be more active. In conclusion, the primary mechanisms responsible for regulating star formation in low-mass galaxies are still unclear and debated. Theoretical models offer a wide range of solutions to this problem, and multiwavelength observations are needed to test them.

<footnote>* Corresponding author; eleonora.parlanti@sns.it</footnote>
In the last ten years, observations and simulations have shown that galactic outflows include gas at different temperatures, densities, and states (molecular, neutral atomic, and ionized) (e.g., Muratov et al. 2015; Janssen et al. 2016; Li et al. 2017; Nelson et al. 2019; Fluetsch et al. 2019, 2021). In particular, four components have been identified: (a) the hot (T ~ 10^{6-7} K) ionized outflowing gas, which emits in the X-rays; (b) the warm ionized component (T ~ 10^4 K), which is usually detected as blueshifted broad wings in the brightest optical emission lines, such as [OIII]λ5007 Å and Hα; (c) the cool neutral atomic part (T ~ 100 K) that is usually observed via NaID λλ5890,5896 Å in the local Universe and [CII]λ158 µm<fnref n="1" /> in high-z galaxies; and (d) the cold molecular component (T ~ 10 K), traced usually by the broad wings of the CO molecular transitions or by P-Cygni profiles of molecular transitions (e.g., OH 119 µm and H₂O). Aside from a few studies, the majority of the outflows are observed only through one tracer, making it difficult to estimate the global outflow properties such as the total outflow mass and mass loading factor, hence the effect on the host galaxy.

While at cosmic noon ($1 \le z \le 3$) outflows are usually identified from broad, often blueshifted components of optical emission lines (e.g., Coatman et al. 2019; Förster Schreiber et al. 2019; Villar Martín et al. 2020; Concas et al. 2022; Cresci et al. 2023), at high redshift ($z > 3$) outflow studies have focused on the far-infrared (FIR) emission lines. This is because optical emission lines are redshifted to near-infrared (NIR) wavelengths, which become hard to access from the ground. In particular, Atacama Large Millimetre Array (ALMA) observations have revealed outflow signatures in molecular transitions (P-Cygni profiles of OH 119µm and broad CO and H₂O line wings; Herrera-Camus et al. 2019; Jones et al. 2019; Butler et al. 2023) and the atomic carbon transition, [CII] (Gallerani et al. 2018; Ginolfi et al. 2020; Bischetti et al. 2019; Herrera-Camus et al. 2021; Solimano et al. 2024; Tripodi et al. 2023). However, the majority of these studies have mainly focused on AGN host galaxies, massive submillimeter galaxies, and lensed dusty star-forming galaxies (Jones et al. 2019; Spilker et al. 2020). Evidence for SF-driven neutral outflows in galaxies in the early Universe has been observed for the first time by stacking the [CII] emission of samples of galaxies at $z \sim 4-5$ (Gallerani et al. 2018; Ginolfi et al. 2020). So far, broad [CII] in the spectra of individual main-sequence galaxies has been reported only for two $z > 4$ sources (Herrera-Camus et al. 2021; Solimano et al. 2024).

With the advent of the James Webb Space Telescope (JWST), we finally have access to the rest-frame optical emission lines of high-redshift galaxies, which allows us to study the outflows of high-z galaxies with the same tools developed for low-z ones. This also gives us the possibility to simultaneously characterize up to two (or three, see above) outflow components in the same target. Early studies with JWST have shown that a non-negligible fraction of star-forming galaxies – about 20–30% – exhibit ionized outflows, even up to the cosmic dawn (Carniani et al. 2024; Zhang et al. 2024; Calabrò et al. 2024). Therefore, the synergy between ALMA and JWST promises to be a rich probe into the physical nature of high-z outflows.

In this work, we focus on HZ4 (RA = 09h58m28.5s, Dec = +02d03m06.7s), which is the most distant star-forming galaxy with evidence of neutral outflows through the [CII] emission line (Herrera-Camus et al. 2021). HZ4 was

identified as a Lyman Break Galaxy (LBG) at $z \sim 5.5$ in the Cosmic Evolution Survey (COSMOS) field (Scoville et al. 2007) and was spectroscopically confirmed with Keck Deep Extragalactic Imaging Multi-Object Spectrograph (DEIMOS) observations (Mallery et al. 2012) which detected Lyα emission (Cassata et al. 2020; Le Fèvre et al. 2020). It was then followed up with ALMA observations to target the [CII] and dust continuum emission with low- (1.01" × 0.85") and high- (0.39" × 0.34") angular resolution (Capak et al. 2015; Béthermin et al. 2020; Herrera-Camus et al. 2021; Jones et al. 2021). In the ALPINE survey (Le Fèvre et al. 2020; Béthermin et al. 2020), HZ4 is dubbed DEIMOS_COSMOS_494057. The high-resolution ALMA [CII] observations characterize the source as a prototypical high-z, turbulent, but rotating disk galaxy (Herrera-Camus et al. 2022; Parlanti et al. 2023). These spatially resolved observations also reveal the presence of a broader component in the spectrum of the [CII] line emission which was interpreted as a neutral outflow (Herrera-Camus et al. 2021). JWST observations HZ4 allow us to investigate the warm ionized gas phase and compare its properties with those of the cold gas.

This work is structured as follows: in Sect. 2 we describe the JWST and ALMA observations and data reduction, in Sect. 3 we describe the analysis of the integrated rest-frame optical spectrum, in Sects. 4 and 6 we analyze the JWST and ALMA data at a spaxel level. In Sect. 5, we present the analysis of the spectral energy distribution (SED) and its results. In Sect. 7, we discuss possible interpretations of the broad emission lines that we observe. In Sect. 8 we speculate that the broad components that we observe are tracing a multiphase outflow, hence we analyze and discuss the outflow properties and their impact on the galaxy. Finally, in Sect. 9 we draw our conclusions. Throughout this work, we adopt the cosmological parameters from Planck Collaboration XIII (2016): $H_0 = 67.7 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.307$, and $\Omega_\Lambda = 0.691$, giving $1'' = 6.09 \text{ kpc at } z = 5.54$.

## 2. Observations

### 2.1. JWST data

HZ4 was observed with the NIRSpec instrument (Jakobsen et al. 2022) onboard JWST as part of the GA-NIFS<fnref n="2" /> survey (program 1217 – Integral Field Spectroscopy in COSMOS, PI: Nora Luetzgendorf) on 23 April 2023. The galaxy was observed in the Integral field spectroscopy (IFS) mode (Böker et al. 2022) with both the G395H/F290LP and PRISM/CLEAR grating/filter combination with an eight-dither position of the medium cycling pattern. The total integration time was ~5 h for the high spectral resolution observations ($R \sim 1670–3600$, hereafter also referred as R2700), which cover the wavelength range between 2.87–5.27 µm, and ~1 h for the prism observations ($R \sim 30–400$, hereafter also referred as R100), which span the wavelength range between 0.6–5.3 µm.

We retrieved the raw data from the Mikulski Archive for Space Telescopes (MAST) archive, then we processed them with the JWST pipeline version 1.11.1 and Calibration Reference Data System (CRDS) context jwst_1097.pmap. We run the standard steps of the public JWST pipeline to produce the final cubes. In the pipeline process, we used some customized code to improve the quality of the final data (see Perna et al. 2023). In particular, during the first stage “calwebb detector 1” that accounts for detector level corrections, we also corrected

<footnote n="1">¹ The fraction of [CII] emission arising from ionized gas is only 1–30% of the total emission (Díaz-Santos et al. 2017; Cicone et al. 2018).</footnote>

<footnote n="2">² https://ga-nifs.github.io</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Integrated [OIII]λ5007 Å flux map of HZ4 from the JWST NIR-Spec/IFS R2700 cube. We show the 5, 25, and 75σ of [OIII] as green contours. Red contours illustrate the [CII] emission from the ALMA Briggs-weighted (robust = 0.5) cube, at 5, 15, and 25σ, respectively. HST/WFC3 F160W image is reported in black contours. The pink, white, and green ellipses represent the ALMA beam size, the HST F160W point spread function (PSF), and NIRSpec IFS PSF at 3.2 µm estimated by D'Eugenio et al. (2024), respectively.</figcaption></figure>

for the 1/f noise. At the end of the stage “calwebb_spec2” that calibrates the data, we visually inspected the calibrated exposures, and we manually masked all regions affected by cosmic rays, persistences, and failed open shutters in the MSA mask. Moreover, we removed the outliers by calculating the derivative of the count rate along the dispersion direction and removing all data where the measurement was above the 98th and 95th percentage of the distribution for R2700 and R100, respectively (D'Eugenio et al. 2024). Then, we applied the third stage, “calwebb_spec3” to combine each exposure and create the final datacubes with a spaxel size of 0.05” by using the “drizzle” weighting. Finally, we subtracted the background emission from each cube by removing the median spectrum computed in the target-free regions of the cubes.

## 2.2. ALMA data

The [CII] raw data were retrieved from the ALMA archive (2012.1.00523.S PI: Capak; 2017.1.00428.L, PI: Le Fevre; 2018.1.01605.S, PI: Herrera-Camus). The source was observed for ~20 min, ~30 min, and 4.7 h, respectively, in Band 7. We used CASA (McMullin et al. 2007) to combine these three observations and calibrate the uv visibilities by using the script included in the dataset. The measurement sets were then processed with CRISTAL's reduction pipeline (see Posses et al. 2024; Solimano et al. 2024; Villanueva et al. 2024; Herrera-Camus et al., in prep., for details). In brief, first, the continuum emission is subtracted in the visibility space using the task `uvcontsub`, the task `tclean` is run with a Briggs weighting (robust parameter = 0.5 and 2), `hogbom` deconvolver and a spaxel scale of 0.0445". The resulting datacubes have a beam size of 0.33" × 0.29", and 0.45" × 0.41", respectively. The cube with Briggs robust parameter = 2 (natural weighting) reaches a higher signal-to-noise ratio (S/N) at the expense of the angular resolution, while the cube with Briggs robust parameter = 0.5 is a compromise between the natural and the uniform (high angular resolution but lower S/N) weighting, providing us a good S/N and

a higher angular resolution with respect to the natural weighting cube. All cubes have a channel width of 40 km s⁻¹.

ALMA observations also allow us to detect the FIR dust continuum emission at ~160 µm in the galaxy rest frame. The continuum image was created with the task `tclean` and a Briggs weighting with a robust parameter = 2 (natural weighting) to maximize the S/N. The final continuum image has a beam size of 0.46" × 0.41" and a sensitivity of 6 µJy/beam.

## 2.3. Astrometric registration

We verified that the astrometric coordinates of JWST are aligned with those from ALMA. Since many other NIRSpec IFS observations showed an offset between the position indicated in the datacube header and the real position of the source on the sky (e.g., Arribas et al. 2024; Jones et al. 2024a; Lamperti et al. 2024), we assessed whether the astrometry of the JWST cubes of HZ4 was correct or not.

First, we retrieved the F160W images of HST/WFC3 (PID: 13641) from the MAST archive, with astrometry corrected a posteriori by using the Gaia eDR3 catalog with uncertainties on the order of 0.01". We then compared the image created by convolving the R100 cube with the F160W filter response and the image from the F160W filter from HST. We fitted a 2D Gaussian profile to find the coordinates of the center both in the HST and JWST images. The mismatch between the 2 images corresponds to a ~0.10" shift. We then shifted the JWST R100 cube to match the F160W astrometry. We also ensured that the R2700 and R100 cubes were consistently aligned by collapsing the [OIII]λ5007 Å emission line visible in both cubes. After verifying the alignment, we shifted the R2700 cube by the same amount as the R100 cube. We assumed that the ALMA astrometry is correct, considering that the ALMA positional accuracy varies from 5% to 20% of the beam depending on the S/N, which corresponds to an accuracy of 0.015"–0.060". The high S/N ensures that the ALMA astrometry is correct within 1 spaxel (i.e., 0.05"). Fig. 1 shows the [OIII] integrated emission along with the contours from the HST F160W image and the [CII] integrated emission after the astrometric correction.

# 3. Integrated analysis

In this section, we describe the analysis of the R100 and R2700 NIRSpec IFS spatially integrated spectra of HZ4. The spatially resolved analysis is reported in Sect. 4.

## 3.1. Spectral fitting of the R2700 cube

In the top panel of Fig. 2 we show the R2700 spectrum integrated over the regions of the cube with S/N > 3 in the wavelength range between 3.27–3.29 µm, which encompasses the [OIII]λ5007 Å emission. The region from which we extracted the spectrum is highlighted in the inset panel as the red contours. The gray shaded area is the error on the cube computed as the quadratic sum of the error extension in the datacube rescaled to match the RMS in the line-free regions of the cube (the error was rescaled by a factor of ~3, see also Übler et al. 2023).

The G395H spectrum spans the wavelengths between 2.87–5.14 µm corresponding to 4388–7859 Å rest-frame. In the integrated spectrum, we detect the brightest rest-frame optical emission lines Hα, Hβ, [OIII]λ5007,4959 Å, [NII]λ6548,6584 Å, [SII]λ6717,6731 Å and HeIλ5875 Å.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** High resolution spectrum of the source. Upper panel: Integrated spectrum of the source extracted from the regions with $S/N > 3$ in the wavelength range 3.27–3.29 µm encompassing the [OIII]λ5007 Å emission. The extraction region is shown as the red contour in the inset panel. Bottom panels: Zoom-in on the integrated spectrum (black) and our best fit (red). The narrow and broad components are shown as the green and blue dashed lines, respectively. The Hβ and [OIII] complex is shown in the left panel, and the Hα, [NII], and [SII]λ6717, 6731 complex is shown in the right panel. The gray shaded areas mark the uncertainties on the data.</figcaption></figure>

We fitted the integrated spectrum by modeling it with the sum of a power-law continuum and two Gaussian components for each line, one narrow and one broad. The [SII]λ6717, 6731 and Hαλ5875 emission lines were fitted only with one narrow line due to the low S/N. We found that one Gaussian profile is enough to reproduce the observed profile of these low S/N lines. We tied the kinematics of each set of line components (narrow and broad) to have the same velocity and line width. For the narrow component, we left the line width free to vary in the range $0 < \sigma < 250$ km s⁻¹, while the broader component was constrained to have a line width larger by at least 20% than the narrow one, and we imposed the amplitude of the broad component to be lower than 50% of the narrow one, following Carniani et al. (2024). The ratios of the intensities of the [OIII] and [NII] doublets have been furthermore constrained to atomic physics motivated ratios (see Osterbrock & Ferland 2006) of [OIII]λ5007/[OIII]λ4959 = 2.98 and [NII]λ6584/[NII]λ6548 = 2.94.

We have also fitted the spectrum by adopting only one Gaussian component for each emission line, but this model has been revealed to be inadequate to reproduce the data. Indeed, the Bayesian information criteria (BIC) test<fnref n="3" /> (Liddle 2007) highly favors the two components model with ΔBIC = BIC₁ − BIC₂ > 10 (Kass & Raftery 1995), where BIC₁ and BIC₂ are the values of the BIC for the fit with one and two components for each emission line, respectively.

In the bottom panels of Fig. 2 we show the best-fitting results and the residuals of the integrated spectrum analysis for the Hβ and [OIII] complex on the left, and the Hα, [NII], and [SII] complex on the right. The flux of each emission line and the FWHM of the narrow and broad components are reported in Table 1. This analysis allows us to compute the integrated properties of

<footnote n="3">³ By assuming a Gaussian noise, BIC = $\chi^2 + k \ln(N)$ where k is the number of free parameters of the fit, and N is the number of data points.</footnote>
<table><caption>**Table 1.** R2700 line fluxes.</caption><thead><tr><th>Line</th><th>Flux × 10<sup>-20</sup> [erg s<sup>-1</sup> cm<sup>-2</sup>]</th></tr></thead><tbody><tr><td>Hα– Narrow</td><td>1617 ± 35</td></tr><tr><td>Hβ– Narrow</td><td>421 ± 18</td></tr><tr><td>[OIII]λ5007 Å – Narrow</td><td>2268 ± 42</td></tr><tr><td>[NII]λ6584 Å – Narrow</td><td>206 ± 17</td></tr><tr><td>Hα– Broad</td><td>986 ± 75</td></tr><tr><td>Hβ– Broad</td><td>313 ± 56</td></tr><tr><td>[OIII]λ5007 Å – Broad</td><td>1297 ± 85</td></tr><tr><td>[NII]λ6584 Å – Broad</td><td>297 ± 41</td></tr><tr><td>HeIλ5875 Å</td><td>97 ± 32</td></tr><tr><td>[SII]λ6717 Å</td><td>122 ± 9</td></tr><tr><td>[SII]λ6731 Å</td><td>107 ± 9</td></tr><tr><td>Component</td><td>FWHM [km s<sup>-1</sup>]</td></tr><tr><td>Narrow</td><td>180 ± 2</td></tr><tr><td>Broad</td><td>430 ± 10</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Prism spectrum integrated from the region marked by red contours in Fig. 2, and best-fit results. In black we show the data, while the gray shaded region represents the $1\sigma$ error. In blue we show the best-fit emission lines, and in red the best-fit model.</figcaption></figure>

the source and compare it with all previous studies that could not resolve the galaxy morphology. From the narrow component, we estimated a redshift of $5.54455 \pm 0.00002$ for the galaxy, which is in agreement with the redshift measured from the [CII] emission (Béthermin et al. 2020; Herrera-Camus et al. 2021). On the other hand, a broader component (FHWM ~ 430 km s⁻¹) is clearly visible and needed to reproduce all the emission line profiles. This component is slightly blueshifted with respect to the narrow component by $v_{\text{broad}} - v_{\text{narrow}} = -24 \pm 3$ km s⁻¹. The FWHM of the broad component is comparable with the [CII] one, observed by Herrera-Camus et al. (2021) in this galaxy, but we measured a blueshift, instead of the redshift measured by Herrera-Camus et al. (2021). We subsequently discuss the possible interpretation of this component in Sect. 7, and we speculate about it being an outflow in Sect. 8, where we compute the outflow properties.

By assuming case B recombination, the theoretical ratio between the Hα and Hβ fluxes is $F_{H\alpha}/F_{H\beta} = 2.86$ (for a temperature $T = 10^4$ K and a density $n_e = 10^2$ cm⁻³, see Osterbrock & Ferland 2006, for details). Assuming the Calzetti et al. (2000) attenuation law ($R_V = 4.05$), we estimated the dust attenuation of the source from the measured Balmer decrement. We obtain a value of $A_V = 1.0 \pm 0.2$ for the narrow component nebular dust attenuation. By using the same method we also computed the obscuration of the broad component, which is $A_V = 0.3^{+0.9}_{-0.3}$.

To compute the star formation rate (SFR) of the galaxy, we used the calibration by Kennicutt & Evans (2012) that exploits the Hα emission line luminosity. We corrected the narrow Hα flux by the dust obscuration reported above, hence we computed an integrated SFR = $77_{-16}^{+19} M_{\odot} yr^{-1}$, which is consistent within the error with the one derived from the SED fitting by Faisst et al. (2020a) by exploiting optical and UV photometry, of $41_{-15}^{+35} M_{\odot} yr^{-1}$.

By exploiting the ratio between the fluxes of [SII]λ6717 Å and [SII]λ6731 Å, we computed the electron density of the source Sanders et al. (2016) and we obtain a value of $n_e = 270_{-200}^{+260}$ cm⁻³. The estimated electron density is, within the large uncertainties, in agreement with the one obtained from [SII]λλ6717, 6731 from a sample of SF galaxies at $2.7 < z < 6.3$ (Reddy et al. 2023), and other high-z galaxies targeted in indi-

vidual NIRSpec-IFS studies (Jones et al. 2024b; Lamperti et al. 2024; Rodríguez Del Pino et al. 2024).

## 3.2. Spectral fitting of the R100 cube

Following the same approach as for the high-resolution cube, we extracted the spectrum from the prism cube from the same aperture as the R2700 one. We show the R100 integrated spectrum in Fig. 3. The prism spans the wavelength range between 0.6 and $5.2 \mu\text{m}$ allowing us to probe both the rest-frame UV and the optical region of the spectrum. In particular, we can now probe the rest-frame wavelengths < 4500 Å, which are inaccessible in the R2700 cube. In the new range of wavelengths, we detect the [OII]λλ3727,3729, [NeIII]λ3870 and Hγ emission lines, as well as a blue UV continuum. We fitted the integrated spectrum as a sum of a continuum emission and a set of emission lines with single Gaussian profiles, since the lower resolution of the spectrum does not allow us to distinguish between the broad and the narrow component, as we did for the R2700 cube. To account for the Balmer break, we modeled the continuum emission as two different power laws with a discontinuity at 3645 Å rest frame, then we convolved the spectrum with the line spread function of the prism (Jakobsen et al. 2022). The velocities of the lines were tied together, but the velocity dispersion was left free to vary in order to allow for the varying spectral resolution of the PRISM with wavelength. It is important to note that the lines are unresolved in the prism spectrum, hence the line width only reflects the instrument's resolution.

In addition, we set the flux ratio of the [OIII] and the [NII] doublets according to atomic physics prescriptions (see Sect. 3.1). We use only one Gaussian component to fit the [SII]λλ6717, 6731 and the [OII]λλ3727,3729 doublets due to the low resolution, which does not allow us to deblend the contribution from each line of the doublet. Therefore, we obtain only the total flux of the doublet. Similarly, the [NII] emission lines are blended with the Hα emission line, hence it is not possible to disentangle their contribution, and the Hα flux is not as reliable as the one estimated from the high-resolution spectrum.

We report the line fluxes measured from the R100 spectrum in Table 2, and note that within the errors, they agree with those reported for the R2700 cube (after summing the contribution from the narrow and broad components). This demonstrates that the flux calibration largely agrees between the two cubes. The only exception is the flux for the Hα and [NII] emission lines, which are blended together, however the summed flux agree within the uncertainties. From the expected line ratio of
$F_{H\gamma}/F_{H\beta} = 0.466$, we estimated a dust extinction of $A_V = 0.7_{-0.7}^{+1.4}$ assuming that the [OIII]\lambda4363 Å emission line flux, which is blended with the Hγ in R100 observations, is zero. We find that this value agrees within the large uncertainties with the one computed from the Hα and Hβ ratio from the grating. From now on we will use as fiducial value of $A_V$ the one inferred from the R2700 grating.

By exploiting the strong-line metallicity diagnostics based on rest-frame optical line ratios, we can calculate the gas phase metallicity of the galaxy. We adopted the calibrations from Curti et al. (2017, 2020), recently revised for high-z, low-metallicity galaxies (Curti et al. 2024; Laseter et al. 2024). In particular, we exploited the diagnostics $R3 = [\text{OIII}]\lambda5007\AA/H\beta$, $R2 = [\text{OIII}]\lambda3727,3729/H\beta$, $R = 0.47 R2 + 0.88 R3$, $O32 = [\text{OIII}]/[\text{OII}]\lambda3727,3729$, $N2 = [\text{NII}]/H\alpha$, and $O3N2 = ([\text{OIII}]/H\beta)/([\text{NII}]/H\alpha)$ to estimate the abundance of oxygen log(O/H).

The $[\text{OII}]\lambda3727,3729$ emission line doublet is only covered by the R100 observations, which do not allow us to disentangle the broad component contribution. Hence, we can estimate a global metallicity by assuming the total line flux for each line (from the PRISM for [OII]\lambda3727,3729, narrow+broad from the grating for the other lines). We corrected each line flux for the dust obscuration computed from the Balmer decrement as the ratio between the total fluxes (narrow+broad) of Hα and Hβ in the grating ($A_V = 0.7 \pm 0.3$). We estimate a value of $12 + \log(O/H) = 8.34 \pm 0.08$ for the entire galaxy. By using only the R3, N2, and O3N2 calibrators (i.e., those provided by the grating spectrum) we are also able to estimate the galaxy versus broad metallicity. We obtain an abundance of $12 + \log(O/H)_{\text{narrow}} = 8.3 \pm 0.1$ for the entire galaxy and $12 + \log(O/H)_{\text{broad}} = 8.4 \pm 0.1$ for the broad component. Both values are consistent with the value we infer by adding the calibrators that exploit the $[\text{OII}]\lambda3727,3729$ lines.

The value that we obtain is in agreement with metallicity obtained for other high-z galaxies with similar stellar mass ($\log(M_*/M_\odot) = 9.75_{-0.10}^{+0.05}$ for HZ4, see Sect. 5) and merging systems (Nakajima et al. 2023; Curti et al. 2024; Sanders et al. 2024; Venturi et al. 2024). Assuming a solar oxygen abundance of $12 + \log(O/H)_{\text{narrow}} = 8.69$ (Asplund et al. 2021), we obtain a metallicity of $Z_{\text{narrow}} = 0.4 \pm 0.1 Z_\odot$ and $Z_{\text{broad}} = 0.5 \pm 0.1 Z_\odot$ for the narrow and broad components, respectively.

# 4. Spatially resolved analysis

## 4.1. Channel maps and continuum emission

The high-resolution NIRSpec R2700 datacube shows that the galaxy has different components at different velocities. This can be seen from the ten channel maps of [OIII] at different velocities shown in Fig. 4. We set the zero velocity at the systemic redshift measured from the narrow component of the integrated R2700 spectrum. In each channel map, we show the contours at 5, 25, and $75\sigma$ of the map at $v = 0$ km s$^{-1}$ to highlight the different morphologies of the [OIII] emission at different velocities.

The $v = 0$ map shows an elongated morphology, with a central peak that coincides with the peak shown in the UV observations (see also Fig. 1) and an extended tail in the southern direction. From the map at $v = 60$ km s$^{-1}$ it is evident that the elongation of the source is due to the presence of three components. In particular, a large southern clump is visible in the redshifted channels ($v > 60$ km s$^{-1}$). We label this component as HZ4-South. The presence of these different components hints that the galaxy may be a galaxy merger or a clumpy galaxy,

<table><caption>Table 2. R100 line fluxes.</caption><thead><tr><th>Line</th><th>Flux × 10<sup>-20</sup> [erg s<sup>-1</sup> cm<sup>-2</sup>]</th></tr></thead><tbody><tr><td>Hα</td><td>3032 ± 110</td></tr><tr><td>Hβ</td><td>761 ± 36</td></tr><tr><td>Hγ</td><td>325 ± 35</td></tr><tr><td>[OIII]λ5007 Å</td><td>3663 ± 71</td></tr><tr><td>[NII]λ6584 Å</td><td>200 ± 91</td></tr><tr><td>HeIλ5875 Å</td><td>158 ± 29</td></tr><tr><td>[NeIII]λ3870 Å</td><td>185 ± 40</td></tr><tr><td>[OII]λλ3727,3729 Å</td><td>1153 ± 90</td></tr><tr><td>[SII]λλ6717,6731 Å</td><td>242 ± 35</td></tr></tbody></table>

as observed for other high-z sources (Carniani et al. 2018; Arribas et al. 2024; Jones et al. 2024a; Venturi et al. 2024). On the other hand, in the blueshifted channels, a component extending in the north direction and peaking in the maps at $v < -182$ km s$^{-1}$ slightly northern than the $75\sigma$ contours peak at $v = 0$ becomes evident. This component hints at the broad component that we will investigate in the next sections.

We also looked at the continuum emission from these sources to identify these three different components. Fig. 5 shows the continuum emission in three different wavelength ranges obtained from the NIRSpec R100 cube. In the left panel, we show the emission of the rest-frame UV continuum (1–2 µm observed wavelength or 1530–3058 Å rest-frame). The UV emission morphology is similar to the one of the [OIII] line (shown as the black contours); the location of the two peaks coincides. Another source appears visible leftward of the target. We identify this source as a foreground source at $z \sim 2.6$ and we discuss in more detail its properties in Appendix A. From the rest-frame optical continuum in the observed range 3.3–4.2 µm (5045–6422 Å rest-frame) and 4.4–5.2 µm (6730–7950 Å rest-frame), shown in the central and right panels of Fig. 5, respectively. Here the southern, redshifted clump, already visible in the channel maps at $v = 121$ km s$^{-1}$, is now evident. This clump is very faint in [OIII] emission (and any other emission line) but appears visible in continuum. The location of the central, small clump already identified in the channel maps at $v = 60$ km s$^{-1}$ is coincident with the location of the peak of the [CII] emission (see Fig. 1) and the peak of the SFR discussed in Sect. 4.

From the maps presented in this section, we can distinguish three different galaxy components. A northern component (HZ4-North) which is very bright in rest-frame UV and optical continuum and [OIII] emission, a central component (HZ4-Central) which is bright in [CII] emission, and a redshifted southern component (HZ4-South), which is faint in [OIII] emission but with an extended tail visible only in the rest-frame optical continuum.

## 4.2. Spaxel-by-spaxel R2700 results

The high S/N of the cube at the spaxel level allows us to perform a spatially resolved analysis of the galaxy. We fitted the cube at a spaxel-by-spaxel level with both a single Gaussian component for each emission line and two Gaussian components, following the same approach discussed in Sect. 3.1. We then selected for each spaxel the best fit by performing the BIC test and selecting the model with two components for each line only when $BIC_1 - BIC_2 > 10$, which represents positive evidence toward the model with two components (Kass & Raftery 1995). In the other cases,
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Channel maps from the R2700 cube targeting the [OIII] λ5007 Å emission line. Black lines are the contour of the 5, 25, and 75 times the RMS level computed in the line-free regions of the cube for the emission at $v = 0$. The velocity marks the central velocity of the channel. The black cross shows the location of the [CII] emission peak (see Fig. 1). The white circles represent the position of the components of HZ4.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Continuum emission from the R100 JWST/NIRSpec datacube. From left to right the average continuum emission between 1–2 µm (1530–3058 Å rest-frame), 3.3–4.2 µm (5045–6422 Å rest-frame), 4.4–5.2 µm (6730–7950 Å rest-frame), respectively. The magenta circles represent the apertures from which we extract the spectrum for the north, central, and southern components for the analysis. The black contours represent the [OIII] emission at 5, 25, and 75σ (see Fig. 1). The black cross shows the location of the [CII] emission peak (see Fig. 1). The *x* and *y* axes represent the displacement in arcsec with respect to the galaxy center (RA = 09h58m28.5s, Dec = +02d03m06.7s, see Sect. 1).</figcaption></figure>

we selected the model with only one component, as the addition of the broader component was not significantly improving the quality of the fit. We excluded all the spaxels with $S/N < 3$ for each emission line. The S/N is computed as the ratio between the flux of the best-fit emission line and the resulting error from the fit.

The results of the fitting (flux, velocity, and velocity dispersion maps) for the narrow and broad components are shown in Fig. 6. The galaxy morphology traced by the rest-frame optical emission lines differs from what is observed in the rest-frame UV continuum and the [CII] FIR emission line, as already seen in the channel maps in Fig. 4. The narrow [OIII] flux map shows a bright clump in the northern region (HZ4-North). This clump is in the same location as the bright emission observed in the F160W filter from HST (see also Fig. 1), which is tracing rest-frame UV emission. Additionally, the [OIII] map shows an extended tail, in the southwest direction. The tail in [OIII] shows the presence of a second faint component (HZ4-Central), which is instead more prominent and visible in the $H\alpha$ map, and is also where the [CII]

peak. This suggests that the second component is likely tracing a region with higher dust obscuration, as we discuss later in this section.

The narrow component velocity field (top-central panel of Fig. 6) shows a different picture from what was seen in ALMA observations. Whereas the [CII] kinematics is consistent with a rotating disk (Herrera-Camus et al. 2021), the ionized gas shows a more complex kinematic pattern, also due to the higher spatial resolution of the JWST data. In the northern region of the galaxy, we observe a relatively shallow velocity gradient ranging from 50 to –20 km s⁻¹, while the southern region is redshifted by ~100–130 km s⁻¹. The clumpy morphology of the galaxy, combined with the observed abrupt velocity gradient in the southern part, suggests that the galaxy is undergoing a merger event rather than being a rotating disk as inferred from ALMA observations (Jones et al. 2021; Herrera-Camus et al. 2022; Parlanti et al. 2023).

In the second row of Fig. 6 we show the broad component maps. In the broad flux map, we also show as green lines the narrow component 5, 25, and 75σ contours. The morphology of
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Results for the spaxel-by-spaxel fitting of the R2700 cube. In the first row, we show from left to right the [OIII] λ5007 Å flux, velocity, and velocity dispersion maps of the narrow component, respectively. In the second row, we show the same but for the [OIII] broad component. In the [OIII] broad flux map, we show the 5, 25, and 75 σ contours of the narrow [OIII] flux in green. In the third row we show from left to right the Hα flux, velocity, and velocity dispersion maps of the narrow component, respectively. In the bottom panels, we show from left to right the flux maps of the Hβ, and [NII] λ6584 Å narrow components, respectively. The plus signs show the location of the three components identified in Figs. 4 and 5.</figcaption></figure>

the broad emission shows an elongated shape in the same direction as the narrow component emission. The brighter region is also in the northern part, but the peaks of the broad and the narrow emission are not in the same location. It is evident that the narrow component peaks ~0.1" to the south of the broad emission, which corresponds to a projected distance of 0.6 kpc. The velocity map shows that in the majority of the spaxels the broad component is blueshifted by up to –130 km s⁻¹ relative to the galaxy-integrated systemic redshift computed from the narrow component. The velocity shift relative to the narrow component increases in the southern part of the broad component. The southern part coincidentally has also the highest velocity dispersion.

In Fig. B.1, we assumed that the broad component is tracing an outflow, and we show the de-projected outflow velocity calculated as $v_{\text{out}} = |\Delta v_{\text{narrow,broad}}| + 2\sigma_{\text{out}}$. The velocity reaches its maximum of approximately 600 km s⁻¹ at the position of the HZ4-Central component. This location coincides with the highest Σ$_{SFR}$, which is indicated by dashed contours in Fig. B.1 and will be discussed at the end of this section. The presence of the higher velocity at the location of the maximum SFR and the fact that the peak flux of the broad component does not coincide with the peak of the narrow one lead us to interpret that this broad component may be an outflow that is launched from the HZ4-Central component. Another possibility could be that we are observing multiple winds driven by the SF regions distributed in
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Top panels: BPT diagram on the left and VO87 diagram on the right. As pink and turquoise crosses we show the location on the BPT diagram of the narrow and the broad component for the integrated spectrum. As red and green points we report the location in the BPT diagram of each spaxel for the narrow and broad component, respectively. The gray shaded region shows the location of SDSS galaxies, while blue diamonds show the location in this diagram of low-metallicity high-redshift ($z \sim 5.5-9.5$) star-forming galaxies from Cameron et al. (2023). Black dots show the location of the MAPPINGS III shock model from Allen et al. (2008) with shock velocities between 200 and 1000 km/s, densities in the range $0.01 \text{ cm}^{-3} \le n \le 1000 \text{ cm}^{-3}$, solar abundances, and magnetic field between $10^{-4}$ and $10 \text{ µG}$. As black dashed and dash-dotted lines, we report the demarcation line for the local star-forming galaxies and AGN samples by Kewley et al. (2001) and Kauffmann et al. (2003). As the dash-dotted black line we show the demarcation line between high-$z$ low-metallicity galaxies and AGNs by Scholtz et al. (2023). AGNs are right-ward of the line, while on the left there is a mixed population of AGNs and SF galaxies. Bottom panels: from left to right the logarithm of N2, R3 in each spaxel for the narrow component, and the logarithm of N2, R3 in each spaxel for the broad component.</figcaption></figure>

the system or, we are witnessing merger flows or shocks. We will discuss the possible scenarios in Sect. 7, while we will explore more in detail the outflow possibility in Sect. 8.

In Fig. 6 we also present the [NII] flux map, which shows a different morphology from the flux maps of the other emission lines. We can highlight the presence of two peaks in the [NII] emission, one associated with the northern component, and one associated with the central one. Differently from the $H\alpha$ and [OIII] the brightest [NII] emission arises from the central region. These differences can be a sign of a different metallicity between the two galaxies or a different ionization mechanism, with the presence of a harder radiation field arising in the center.

To study this latter possibility we exploit the Baldwin-Phillips-Terlevich (BPT; Baldwin et al. 1981) and the VO87 (Veilleux & Osterbrock 1987) diagnostic diagrams, which are sensitive to the ionizing source in a galaxy, between star formation in HII regions and AGN photons. In particular, the first one is based on the line ratios of [NII]/$H\alpha$ and [OIII]/$H\beta$ while the latter uses the [SII]$\lambda6717,6731/H\alpha$ and [OIII]/$H\beta$. In the top-left panel of Fig. 7 we report the spatially integrated measurement of the galaxy in the BPT diagram for the narrow and the broad component, together with the location of each spaxel.

We observe that the integrated emission and the majority of the spaxels, of the narrow component, lie on the Kewley et al. (2001) demarcation line, together with the low-metallicity $z \sim 5.5-9.5$ galaxies from Cameron et al. (2023). On the other hand, the majority of the broad component points lie in the composite region of the BPT between the demarcation lines of Kewley et al. (2001) and Kauffmann et al. (2003) that are based on the SDSS galaxy sample (gray shaded region). Low-redshift observations have shown that in the presence of shocks induced by starburst-driven winds, or by tidal flows due to a merger, the line ratios can shift toward the composite region even in the absence of an AGNs (Monreal-Ibero et al. 2010; Rich et al. 2011, 2015; Medling et al. 2018).

To check for the presence of shock, we over plot the MAPPINGS III shock-ionization models (Allen et al. 2008) with densities of $n = 100$ and $1000 \text{ cm}^{-3}$, velocities ranging from 200 to $1000 \text{ km/s}$ and magnetic fields from $10^{-4}$ to $10 \text{ µG}$. We observe that the shock models cannot reproduce the location of the broad line ratios in the BPT diagram, with the exception of just a few extreme points. We note that the presence of a large-scale galactic outflow can shock-excite the galaxy interstellar medium (Ho et al. 2014), increasing the [NII]/$H\alpha$ line ratio.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Spatially resolved properties of the narrow component. From left to right: gas-phase metallicity (12 + log(O/H)), dust extinction ($A_V$), and dust-corrected star formation rate surface density ($\Sigma_{SFR}$) maps. The plus signs show the location of the three components identified in Figs. 4 and 5. The gray contours represent the 5, 25, and 75σ of the narrow [OIII] flux from Fig. 6.</figcaption></figure>

The demarcation lines used in the BPT diagram for low-$z$ sources have proved to be unreliable at high-$z$ ($z > 3$). Due to the increase of the ionization parameter and the decrease of the metallicity at high redshift, both galaxies and AGNs may end up in the same location in the BPT diagram (Nakajima & Maiolino 2022; Übler et al. 2023; Harikane et al. 2023; Maiolino et al. 2024; Scholtz et al. 2023). We then used the new demarcation line by Scholtz et al. (2023), shown as a dash-dotted line in the plot. According to this new line, there are no points that lie in the AGNs region, but all points lie in the region in which SF galaxies and AGNs may have the same line ratios.

In the bottom-right panel, we show where the narrow component lies in the VO87 diagram, together with the sample of high-$z$ galaxies from Cameron et al. (2023) and the local galaxies from SDSS. The point lies on the demarcation lines by Kewley et al. (2001) and Scholtz et al. (2023), therefore we cannot constrain the ionization source.

Unfortunately, for the wavelength range and spectral resolution available, there are no other diagnostic diagrams that have been proven to be reliable at high redshift (Scholtz et al. 2023; Mazzolari et al. 2024). The absence of HeII and other high-ionization lines rule out the presence of a strong and energetic AGN and favor a SF interpretation, although we cannot completely rule out a contribution to ionization from a weak, sub-dominant AGN.

We display the maps of the log([NiII]/H$\alpha$) and log([OIII]/H$\beta$) ratios for the narrow and broad components in the top row of Fig. 7. For the narrow component, a gradient in the values of log([NiII]/H$\alpha$) is evident, with the ratio increasing from north to south, while the values of log([OIII]/H$\beta$) are almost constant. The broad component, on the other hand, shows a gradient also in the log([OIII]/H$\beta$) values, varying from 0.2 to 0.7, while the log([NiII]/H$\alpha$) values range from -0.8 to -0.2. Since the AGN scenario is disfavored, the observed gradients of the line ratio can be interpreted as a metallicity gradient between the southern and the northern clump with the northern clump being poorer in metals compared to the southern one, given that the [NiII]/H$\alpha$ ratio decreases as the metallicity decreases (Pettini & Pagel 2004; Curti et al. 2017, 2020). Metallicity gradients have been observed at high redshift, especially in clumpy or merging systems (Arribas et al. 2024; Rodríguez Del Pino et al. 2024; Venturi et al. 2024). The study of the metallicity gradients for this galaxy will be presented in Curti et al. (in prep.), but here we show the metallicity maps for the narrow component in the left panel of Fig. 8. We exploit the N2, R3, and O3N2 line ratios and we adopt the calibrations from Curti et al. (2017) to derive the O/H abundances. The estimated abundances vary between 8.2 < 12 + log(O/H) < 8.5. The lowest metallicity values are

near the center of the northern and central components, while for the southern component, the lack of detection of [NiII] emission does not allow us to constrain the metallicity.

Having excluded the presence of an AGN, we can also estimate the SFR surface density ($\Sigma_{SFR}$) at a spaxel level by exploiting the H$\alpha$ dust-corrected luminosity. First, we estimated the dust attenuation at a spaxel level thanks to our detection of H$\alpha$ and H$\beta$ emission lines. By assuming case B recombination ($T_e = 10\,000\,\text{K}$ and $n_e = 100\,\text{cm}^{-3}$), the theoretical ratio between the H$\alpha$ and H$\beta$ flux is $F_{H\alpha}/F_{H\beta} = 2.86$ (Osterbrock & Ferland 2006). We measured a spaxel-by-spaxel ratio ranging from 2.86 to 5.60. Assuming the Calzetti et al. (2000) attenuation law, we estimate a dust attenuation of $0 \le A_V \le 2.5$. In the central panel of Fig. 8 we show the map of dust attenuation across the galaxy.

Interestingly, the northern region, where the line emission and the UV continuum peak, shows little to no dust extinction, similar to what is observed in the southern component. However, there is a zone of high extinction between these two regions where the dust attenuation is reaching values of $A_V \sim 2.5$ which is spatially with the central component.

By correcting the H$\alpha$ flux for the dust extinction and using the calibration from Kennicutt & Evans (2012) to get the SFR from H$\alpha$ luminosity, we computed the resolved SFR surface density map, which is shown in the right panel of Fig. 8. This shows that there is a region of intense and obscured ongoing star formation in the galaxy, which is coincident with the location of the central component. By comparing with the [CII] flux contours shown in Fig. 1, it can be seen that this region is coincident with the peak of the [CII] emission, which indeed is a tracer of star formation unaffected by the presence of dust (De Looze et al. 2011; Pineda et al. 2014; Herrera-Camus et al. 2015).

## 5. SED fitting

The prism data allow us to perform a spectral energy distribution (SED) fitting to study the galaxy's stellar population and infer key galaxy properties such as the stellar mass and the star formation history. We used the SED fitting code *Bagpipes* (Carnall et al. 2018). We adopted the stellar population SED model from Bruzual & Charlot (2003), and we used *Cloudy* (Chatzikos et al. 2023) to obtain the models for the nebular emission with a range in the ionization parameter priors of -3 < log U < 0. We assumed a Calzetti et al. (2000) dust attenuation law, with $A_V$ free to vary between 0 and 4. We assumed the same extinction for the nebular and stellar components. We left the metallicity free to vary within the range 0.1–1 Z⊙. For the star formation history (SFH), we adopt a so-called nonparametric SFH, with continuity priors (Leja et al. 2019) with 10 bins
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** On the left: R100 spectra in black and SED best-fitting result in red extracted from the entire system, the northern clump, the central clump, and the southern clump, from top to bottom, respectively. On the right: best-fit star formation history.</figcaption></figure>

logarithmically spaced from 1 Gyr ($z = 5.5$, the redshift of the source) up to 200 Myr after the Big Bang.

First, we performed the SED fitting on the integrated R100 spectrum of the regions in which the S/N on the [OIII] emission line is larger than 3 (see red region in the inset panel in Fig. 2). We show the spectrum and the best-fit model in the top-left panel of Fig. 9, while we show the SFH in the top-right panel. The properties resulting from the SED fitting are reported in Table 3. We infer a total stellar mass of the system of $log(M_*/M_\odot) = 9.75^{+0.05}_{-0.10}$ which is smaller than the previous estimate from the literature of $log(M_*/M_\odot) = 10.15 \pm 0.15$ of Faisst et al. (2020a), but in agreement with the one from Capak et al. (2015) of $log(M_*/M_\odot) = 9.67 \pm 0.21$.

The SFH for the whole HZ4 system shows an almost constant SFR of a few $M_\odot$ yr$^{-1}$ up to the last 20 Myr, where we see a rapid increase reaching ~100 $M_\odot$ yr$^{-1}$ that could be triggered by the merger (Pearson et al. 2019). Indeed, the integrated spectrum shows the Balmer break at ~2 µm, which is indicative of an older stellar population, while the luminous $H\alpha$ is indicative of SFR in the last 10 Myr.

If we extract the spectra for the northern, central, and southern components (see apertures in Fig. 5) we see that the older population is dominant in the southern component, with a well-defined Balmer break, faint continuum UV emission, and faint emission lines. On the other hand, the spectra extracted in the central and northern regions show a bright UV continuum. This difference in the spectra indicates a different stellar population that reflects different formation histories of the northern, central,

<table><caption>**Table 3.** Galaxy properties inferred from the SED fitting.</caption><thead><tr><th></th><th>log(M<sub>★</sub>/M<sub>☉</sub>)<br>(1)</th><th>A<sub>V</sub><br>(2)</th><th>τ/Myr<br>(3)</th></tr></thead><tbody><tr><td>HZ4 System</td><td>9.75<sup>+0.05</sup><sub>-0.10</sub></td><td>0.39 ± 0.04</td><td>357<sup>+69</sup><sub>-95</sub></td></tr><tr><td>HZ4-North</td><td>8.90 ± 0.05</td><td>0.28 ± 0.03</td><td>27<sup>+43</sup><sub>-17</sub></td></tr><tr><td>HZ4-Center</td><td>8.69 ± 0.06</td><td>0.25 ± 0.04</td><td>50<sup>+70</sup><sub>-32</sub></td></tr><tr><td>HZ4-South</td><td>9.24 ± 0.08</td><td>0.71 ± 0.13</td><td>322<sup>+88</sup><sub>-108</sub></td></tr></tbody></table>

**Notes.** (1) Stellar mass; (2) Dust attenuation; (3) Mass-weighted stellar age.

and southern components (see Fig. 9, right panels). The southern component started to form stars only ~200 Myr after the Big Bang and has continued to grow with a slowly increasing SFR, reaching values of SFR = 8±2 $M_\odot$ yr$^{-1}$ at the time of the observation. From the SED fitting we obtain a mass-weighted age of the southern galaxy of $322^{+88}_{-108}$ Myr. On the other hand, the northern galaxy shows a much quicker rise in the SFR, having assembled most of its mass in the last 30 Myr. The central component has a rising SFH very similar to the one of the northern one. Both the northern and the central components have a very young stellar population, with a mass-weighted age of <100 Myr.

The differences in the SFH and the stellar populations support our interpretation that this galaxy is formed by merging
systems, rather than three star-forming clumps. Cosmological
simulations from Nakazato et al. (2024) show that a merger
event can induce the formation of clumps which lead to a bursty
star formation with timescales of 10–30 Myr, consistent with
what we observe in the northern and central clumps.

6. ALMA analysis

In this section, we describe the analysis of the ALMA obser-
vations of the [CII] emission line and the dust continuum. The
ALMA [CII] cube was analyzed with the same method as the
NIRSpec R2700 cube, to detect and possibly spatially resolve
the broad component seen in the neutral gas phase. Since the
properties of the outflows integrated in one beam are already dis-
cussed in Herrera-Camus et al. (2021), here we focus only on the
spaxel-by-spaxel analysis. To study the outflow component we
exploited the datacube which was reduced with natural weight-
ing, which has a higher S/N at the expense of angular resolu-
tion with respect to the Briggs datacube. We fitted each spaxel
having a S/N of the [CII] emission greater than 3 with both a
single and two Gaussian lines and a constant continuum. The
continuum emission was already removed from the cube during
the data reduction, but we added this component in the fitting to
furthermore remove any residual continuum emission that was
not accounted for. Using a similar method as described for the
JWST datacube, for each spaxel we selected the 2-component
model only when BIC<sub>1</sub> – BIC<sub>2</sub> > 5 and the S/N of both the nar-
row and the broad component were greater than 3, otherwise the
1-component model was adopted. We note that in this case, we
use a slightly lower threshold for the BIC selection due to the
lower S/N of the line, but ΔBIC > 5 still gives positive evidence
in favor of the two Gaussian model. We applied the same pro-
cedure to the Briggs (robust = 0.5) datacube. Unfortunately, the
lower S/N of this latter datacube does not allow us to select the
outflow component in any spaxel due to its S/N falling below
our threshold of 3, even if the ΔBIC was larger than 5 in some
spaxels. But a broad component is still detected integrating into
circular apertures as also shown by Herrera-Camus et al. (2021).

6.1. [CII] spaxel-by-spaxel results

Here we present the results of the spaxel-by-spaxel analysis of
the ALMA data targeting the [CII] emission line. In Fig. 10 we
show the best-fit results for the [CII] emission line in one spaxel.
We present the best-fit results by using one component to fit the
emission line in the right panel, and by using 2 components in the
left panel. The 1-component fit shows an excess around a veloc-
ity of ~200 km s⁻¹. These residuals completely vanish when the
second component is added. The 2-component model is favored
by the BIC test, with a ΔBIC = 14. We perform this test on
each spaxel and we end up selecting the broad component for
115 spaxels.

In Fig. 11 we show the results for the spaxel-by-spaxel fitting
of the [CII] data. In the first two rows we show moment maps
of the narrow and broad components from the natural weight-
ing reduced cube, while in the panels in the bottom row, we
show the moment maps from the datacube reduced with Briggs
weighting, robust = 0.5. The cube reduced with natural weight-
ing has a larger beam, but a higher S/N on a spaxel level with
respect to the Briggs, robust = 0.5 one. The higher S/N allows
us to select a broad component, that we show in the central row.
This component is redshifted with respect to the host galaxy.
This redshifted broad component is what was already observed
in circular apertures by Herrera-Camus et al. (2021), similarly

elongated along the northwest direction. In contrast to the [OIII]
broad component, it is redshifted with respect to the host galaxy
(Δv ~ 35 km s⁻¹), but the projection is co-spatial with the [OIII]
broad (shown as the red contours). The broad [CII] component
appears extended toward the filamentary structure (feature in the
northwest at RA ~ 9h58m28.47s, Dec ~ 2d0.3m07.0s) that we
observe in the narrow component natural flux map and in the
map of the total profile. The moment maps of the narrow com-
ponent show an elongated [CII] emission due to the presence of
two separate spatial components, one peaking where the [OIII]
peaks (the north component), one at the location of the central
component, where the integrated [CII] map from Fig. 1 peaks.
The smeared velocity field resembles that of a rotating galaxy,
as already observed in other studies (Herrera-Camus et al. 2022;
Parlanti et al. 2023).

In the bottom panels of Fig. 11 we show the maps of the nar-
row component obtained from the fitting of the Briggs cube. The
Briggs robust = 0.5 cube unambiguously confirms the presence
of two distinct clumps (northern and central components) in the
amplitude map, through the higher spatial resolution. The peak
[CII] flux is in the location of the central obscured component.
This is in agreement with the [CII] being a tracer of star for-
mation as the maximum luminosity arises from the regions with
higher star formation (see also Sect. 5).

6.2. Dust continuum emission

The ALMA observations also allow us to study the resolved
dust continuum emission at ~160 µm in the galaxy rest frame.
While the optical and UV rest-frame emission allows us to probe
directly the light emitted from the young stars, the FIR contin-
uum emission traces the thermal radiation emitted from the dust
that is heated by the stars. As the UV emission is absorbed by the
dust, we can probe the obscured SF in the system by analyzing
the dust emission.

In Fig. 12, we report the dust emission contours at 3, 6 and
8σ on top of the $A_V$ map. The dust emission morphology fol-
lows the elongation of the source, already seen in [OIII] and
[CII] emission. The dust continuum peaks between the northern
and the central components, where the $A_V$ measured through the
Balmer decrement is 0 instead. The beam size (0.46" × 0.42")
is larger than the separation between the two components. The
difference between the $A_V$ map and the dust continuum emis-
sion can be ascribed to a difference in the column density and
the dust temperature across the field of view. The total flux
of the source at 160 µm is $S_{160\mu m} = 134 \pm 33$ µJy, compara-
ble with the values reported in Herrera-Camus et al. (2021) and
Mitsuhashi et al. (2024) with a similar resolution and Faisst et al.
(2020b) with lower resolution.

The observations of the FIR continuum, together with the
UV continuum and Hα emission line, allow us to compare the
SFR computed with different methods. Before the advent of
JWST, the total SFR of high-z galaxies, due to the lack of opti-
cal rest-frame emission observations, was estimated by sum-
ming the unobscured SFR estimated from UV continuum emis-
sion, and the obscured one estimated from the FIR continuum.
The three, low-resolution (~1''), dust continuum detections in
ALMA bands 6, 7, and 8 of HZ4 allow us to perform a FIR SED
fitting to estimate the obscured SFR. Faisst et al. (2020b) report
SFR<sub>UV</sub> = 29 ± 4 M<sub>☉</sub> yr<sup>-1</sup> exploiting the UV HST observations,
and through SED fitting they estimate a SFR<sub>IR</sub> = 77<sup>+104</sup><sub>-69</sub> M<sub>☉</sub> yr<sup>-1</sup>.
Despite three FIR detections, the SFR<sub>IR</sub> has uncertainties of
~1 dex, principally driven by the uncertainties in determining
the dust temperature ($T_D = 57.3_{-16.6}^{+67.1}$ K; Faisst et al. 2020b) due
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Example of best-fit results for the fitting of the [CII] emission line in one spaxel. On the left, we show the best-fit results by using two components, on the right the best-fit results by using one component. In the bottom panels, we show the residuals of the fitting procedure. The dotted lines show the RMS computed in the line-free regions of the cube.</figcaption></figure>

to the lack of observations near the IR peak, which can strongly constrain the temperature.

Thanks to the JWST observations, we were able to determine the galaxy SFR thanks to the observed $H\alpha$ luminosity, corrected for the dust obscuration according to the Balmer decrement. Assuming that the dust-corrected $H\alpha$ luminosity is able to recover the entire intrinsic SFR ongoing in the system, we obtain that the total SFR of the system is $SFR_{H\alpha} = 77_{-16}^{+19} M_{\odot} yr^{-1}$. Of this, $29 \pm 4 M_{\odot} yr^{-1}$ is unobscured SF, and the remaining is obscured SF of $SFR_{obscured} = 48_{-20}^{+23} M_{\odot} yr^{-1}$. By using the relation between the obscured SFR and $L_{IR}$ by Faisst et al. (2020b; Fudamoto et al. 2023).

Then, we can now perform a SED fitting for the three FIR detections (we adopt the values from Faisst et al. 2020b for the band 6 and 8 fluxes and the value reported above for the band 7 detection), but constraining the $L_{IR}$ from the obscured SFR inferred above. We assumed that the IR and $H\alpha$ trace SF over a similar timescale, a dust emissivity index of $\beta = 2$, and we performed a fit assuming a gray-body curve with a single temperature. With these assumptions, we obtain a dust temperature of $T_D = 70 \pm 12$ K and a dust mass of $log(M_D/M_{\odot}) = 6.62 \pm 0.17$. The dust temperature is in agreement with the values of Faisst et al. (2020b) and Fudamoto et al. (2023), but with smaller uncertainties.

## 7. Discussion

Based on previous observations, HZ4 was considered a typical high-z turbulent rotating disk. Our new JWST/NIRSpec IFS data instead reveal a clumpy structure, a different morphology depending on the tracer UV continuum emission, [OIII]$\lambda$5007 Å and other rest-frame optical emission lines, [CII] or dust continuum, with each clump having different SFHs and different properties,

and an asymmetric velocity field. These new observations hint at this galaxy being a merger between at least two galaxies with different velocities, instead of a rotating disk.

The presence of additional broader components, statistically needed to reproduce the observed line profile, also implies that there are noncircular motions both in the rest-frame optical lines, which trace the ionized gas, and in the re-analyzed [CII] emission line, tracing mainly the neutral gas and confirming the previous results by Herrera-Camus et al. (2021). We discuss below different scenarios that can create the observed broad features (see Fig. 13).

- Different outflows in different directions.

This scenario is supported by the observation of a star forming and obscured clump, which thanks to the SF feedback is able to launch a blue-shifted, ionized gas outflow with the energetics that we observe. On the other hand, the redshifted neutral outflow can be launched by the same region or other star-forming clumps (such as the [OIII] bright northern clump), triggered by the merger event. The ionized and molecular outflows are almost co-spatial, but the [CII] outflow extends in the northwest direction, while the ionized outflow extends in the northeast direction from what we assumed as the launching region. However, the $\Sigma_{SFR}$ is high enough to allow the launch of these outflows from the same regions. We note that the observation of the redshifted outflow only in [CII] does not exclude the possibility of it also hosting ionized gas. If the receding side of the outflow lies behind the disk with respect to our line of sight, that side of the ionized outflow could be obscured by the dust in the system in the optical emission lines observed by NIRSpec, but still be observed in [CII] with ALMA. In that case, we would only be able to detect the blueshifted approaching side. We note that a similar behavior has been observed in Solimano et al. (2024, 2025), where a redshifted [CII] plume is visible in ALMA data, and in the launching region there are signs of blueshifted outflow emission shown by JWST data. Additionally, observation of a local merger and starburst indicate a decoupling in the kinematics
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** [CII] moment maps derived from a spaxel-by-spaxel Gaussian fitting. From left to right we present the maps of the amplitude of the Gaussian, the velocity, and the velocity dispersion, respectively. From top to bottom, we show the narrow and the broad components from the natural cube and the total emission from the 1-component fit from the Briggs robust = 0.5 cube. We note that no broad components were selected from the Briggs cube because they fall below the $3\sigma$ threshold. In the left panels, we show the ALMA beam as a green ellipse and the narrow and broad [OIII] flux from Fig. 6 as lime and red contours, respectively.</figcaption></figure>

of ionized and neutral gas. This suggests that outflows, particularly in starburst-driven mergers, may have a complex geometry that cannot be modeled as simple conical outflows (Shih & Rupke 2010).

- Galactic fountain.

The redshifted [CII] emission can also be explained as a galactic fountain. The ionized wind is not fast enough to escape the potential well of the galaxy as we will discuss in Sect. 8.3, hence eventually the expelled material may cool down and then fall back onto the galaxy itself creating the galactic fountain (Leroy et al. 2015). It is expected that the majority of the gas launched outward by SF-driven winds eventually will fall back as galaxy fountains (Bregman 1980; Oppenheimer et al. 2010; Brook et al. 2012; Übler et al. 2014). Given the low outflow velocity observed in comparison to the escape velocity, eventually this will happen, but it is unclear whether we are witnessing it now in these observations. Unfortunately, the presence of a galactic fountain is hard to prove even in the local Universe and requires a detailed kinematic model of the galaxy and the

outflow to prove or disprove its presence (e.g., Yuan et al. 2023).

- [CII] inflow or merger-induced flows and [OIII] outflow. The accretion of cold gas onto galaxies is one of the main mechanisms of galaxy growth and evolution. Inflow signatures are hard to observe and study at high redshift, and are usually seen thanks to absorption features (Rubin et al. 2012; Bouché et al. 2013; Herrera-Camus et al. 2020), or thanks to hints of radial motion via kinematic analysis (Arribas et al. 2024; Genzel et al. 2023; Übler et al. 2024).

If the most probable scenario is that the [OIII] emission traces outflows, the [CII] emission instead can be associated with the inflow of cold atomic and molecular gas from the CGM and IGM. As shown in Fig. 11, the orientation and location of the [CII] broad component is aligned with the filamentary structure seen in the narrow component flux maps, both in the natural and in the Briggs cubes, which can resemble an accreting filament.

With this scenario the accreting gas must be already enriched in [CII], and thus cannot be pristine gas, hence the gas could
<figure><img src="image_14.png" /><figcaption>Fig. 12. $A_V$ map from narrow component (same as in Fig. 8) and dust continuum emission contours at 160 µm rest-frame. In black we report the dust continuum 3, 6, and 8σ contours. The black plus signs represent the location of the three components identified from the [OIII] channel maps and in optical continuum emission. The purple ellipse represents the beam size of the dust continuum emission.</figcaption></figure>

originate from the galaxy, going back to the fountain scenario.

Although streams of inflowing gas may be present in the system, it is unlikely that they dominate the observed broad line profiles, which imply motions of ~300 km s⁻¹. While these velocities are typical of SF-driven outflows, they are too large for the inflow velocities that are expected to be a fraction of the virial velocity ($v_{vir} = \sqrt{GM_{DM}/r_{vir}} \sim 180$ km s⁻¹; Goerdt & Ceverino 2015).

Another possibility is that the broad wings observed may be explained as due to gravitational interaction (i.e., tidal tails) between these galaxies in the process of merging (e.g., Baron et al. 2024).

- Shock due to the merging.

The merger between the multiple components in the system can also trigger shocks, that increase the velocity dispersion of the gas. In particular, during a major merger, tidal forces can compress and shocks the gas, triggering strong episode of star formation (Pearson et al. 2019). Observationally this can be observed as an increase in line ratios of low-ionization components and an increase in the velocity dispersion usually reaching in the local Universe $\sigma \sim 100-200$ km s⁻¹ (Rich et al. 2011, 2015). In our observations we find velocity dispersion of the broad component compatible with shocked regions in local galaxies (see Fig. 6), however, we do not find line ratios consistent with the presence of shocks (see Fig. 7), at least using the BPT and VO87 diagnostic diagrams. The [OI]λ6300 emission line would be a better tracer of shocks but falls into the detector gap in our high-resolution dataset, and it is undetected in the R100 dataset. However, we can estimate a 3-σ upper limit from the PRISM data of log([OI]λ6300/Hα) < 10⁻¹ for the total profile, which discourages the presence of strong shocks dominating in the ionization of the system (Rich et al. 2011). The current dataset does not allow us to put further constraints on the presence of the shock and their role in driving the observed increase in velocity dispersion of the system.

The present NIRSpec and ALMA observations describe a complex system formed by three merging galaxies, where likely gravitational interactions, inflowing gas, and outflows are present, with the latter explanation being the most probable in explaining the velocities implied by the broad component of the ionized emission line profiles. A better knowledge of its actual geometry would allow us to further constrain the properties of the system and analyze the relative role of these processes.

## 8. Outflow properties

In this section, we assume that the broad emission lines seen in ALMA and JWST data are tracing a neutral and ionized outflow, and we derive its properties. With this assumption, we investigate the origin of the outflows and their launching mechanisms. Moreover, we present and discuss the outflow history.

### 8.1. Ionized outflow properties

Assuming that the measured broad lines are tracing the outflows, we first computed the integrated outflow properties for the ionized gas phase. By using the results from the modeling of the integrated spectrum reported in Table 1, we can compute the ionized gas mass from the flux of the Hα broad emission line component following Genzel et al. (2011), as

$$M_{\text{ion,out}} = 3.2 \times 10^5 \left( \frac{L_{\text{Hα,outflow}}}{10^{40} \text{ erg s}^{-1}} \right) \left( \frac{100 \text{ cm}^{-3}}{n_e} \right) M_{\odot}, \quad (1)$$

and the broad [OIII] emission following Carniani et al. (2024), as

$$M_{\text{ion,out}} = 0.8 \times 10^8 \left( \frac{L_{\text{[OIII],outflow}}}{10^{44} \text{ erg s}^{-1}} \right) \left( \frac{500 \text{ cm}^{-3}}{n_e} \right) \left( \frac{Z_{\odot}}{Z} \right) M_{\odot}, \quad (2)$$

where $L_{Hα,outflow}$ and $L_{[OIII],outflow}$ are the luminosity of the Hα and [OIII] broad components corrected for the dust extinction, $n_e$ the outflow electron density and Z the metallicity of the outflowing gas. We assume $n_e = 200 \text{ cm}^{-3}$, which is consistent with the value derived from the [SII]λ6717, 6731 ratio for the ISM of the galaxy with the integrated fit, and with density measured for outflows in high-z galaxies (e.g., Förster Schreiber et al. 2019), and we adopted an outflow metallicity of Z = 0.5 Z$_{\odot}$ (see Sect. 3.2). By using the value of Hα and [OIII] fluxes reported in Table 1, we corrected for the outflow extinction of $A_V = 0.3^{+0.9}_{-0.3}$ and we estimate an ionized outflow mass of $M_{\text{ion,out}} = (1.8^{+2.0}_{-0.5}) \times 10^7 M_{\odot}$ and $M_{\text{ion,out}} = (2.5^{+4.5}_{-0.8}) \times 10^7 M_{\odot}$ for Hα and [OIII], respectively.

For the outflow velocity $v_{out}$ we adopted the following relation from Rupke & Veilleux (2013) and Fiore et al. (2017):

$$v_{\text{out}} = |\Delta v_{\text{narrow,broad}}| + 2\sigma_{\text{out}} = 389 \pm 12 \text{ km s}^{-1}, \quad (3)$$

where $\Delta v_{\text{narrow,broad}}$ is the velocity shift between the peak of the broad and the narrow components and $\sigma_{\text{out}}$ is the velocity dispersion of the broad component deconvolved for the instrumental broadening. With this assumption, we are accounting for the unknown geometry of the outflow, under the assumption that only the tail of the broad wings trace the outflowing gas directed in (or moving along) our line of sight, and thus its intrinsic velocity.

The mass outflow rate is computed by assuming that the mass outflow rate is constant with time (e.g., Lutz et al. 2020)

$$\dot{M}_{\text{out}} = \frac{v_{\text{out}} M_{\text{out}}}{R_{\text{out}}}, \quad (4)$$
<figure><img src="image_15.png" /><figcaption>Fig. 13. Cartoon illustrating three possible scenarios involving the presence of outflows for the explanation of the observations of the broad lines. In all panels, we draw in red the [CII] noncircular components, and in blue the ionized gas emitting the broad rest-frame optical lines. In the left panel, we represent the multiphase outflow scenario in which the ionized broad line is mainly emitted by the approaching side of the outflow launched by the SF clump. The outflow traced by the [CII] broad component is emitted by outflowing gas on the receding side. In the central panel, we show the galactic fountain scenario, in which the ionized gas is still tracing the outflow, but the [CII] emission is tracing the gas that cannot escape from the potential well of the galaxy, then it cools down and falls back onto the galaxy. The right panel instead depicts the scenario in which the redshifted [CII] emission is tracing the inflow of cold gas from the CGM or IGM onto the galaxy, or is tracing gas stripped by the outer region of the northern and southern galaxies due to gravitational interaction.</figcaption></figure>

where $R_{\text{out}}$ is the extension of the galactic outflow. Based on the spaxel-by-spaxel fit we have inferred a flux-weighted radius of $R_{\text{out}} \sim 2$ kpc (see Fig. 6) that yields $\dot{M}_{\text{out,H}\alpha} = 2.3_{-0.6}^{+2.7} \text{M}_{\odot} \text{yr}^{-1}$ and $\dot{M}_{\text{out,[OIII]}} = 3.3_{-1.1}^{+5.7} \text{M}_{\odot} \text{yr}^{-1}$.

## 8.2. Neutral outflow properties

With the new analysis of the ALMA data, we can also re-estimate the neutral atomic outflow properties under the assumption that the observed broad lines are tracing an outflow. We computed the neutral atomic gas mass in the outflow by using the following equation:

$$M_H = k_{[\text{CII}]}(T, n_H, Z) \times L_{[\text{CII}]} \quad (5)$$

where $k_{[\text{CII}]}$ is a conversion factor that depends on the temperature ($T$), the neutral gas density ($n_H$), and the gas metallicity ($Z$), and $L_{[\text{CII}]}$ is the luminosity of the observed [CII] outflow. Since the outflow temperature and density are unknown we conservatively compute the value of $k_{[\text{CII}]}$ assuming $T = 10^4 \text{ K}$ and $n_H = 10^4 \text{ cm}^{-3}$, which correspond to the maximal excitation case. With this assumption, the mass obtained represents a lower limit on the neutral gas mass traced by the [CII] emission line.

Assuming an oxygen abundance of $12 + \log(\text{O/H}) = 8.34$ ($Z = 0.45 Z_{\odot}$, see Sect. 3.2), the case for maximal excitation gives $k_{[\text{CII}]} = 4.88 M_{\odot}/L_{\odot}$ (Herrera-Camus et al. 2021). The minimum neutral mass in the outflow is then $M_{\text{H,out}} = 2.69 \times 10^8 M_{\odot}$. We compute the flux-weighted radius of the outflow from the map shown in Fig. 11, which gives $R_{\text{H,out}} \sim 1.3$ kpc, and the outflow velocity computed as in Eq. (3) is $v_{\text{H,out}} \sim 260 \text{ km s}^{-1}$. Finally, we compute a lower limit on the mass outflow rate of $\dot{M}_{\text{H,out}} \sim 121 \text{ M}_{\odot} \text{yr}^{-1}$. The value we obtain is higher than the lower limit reported in Herrera-Camus et al. (2021) of $\sim 34 \text{ M}_{\odot} \text{yr}^{-1}$, but this is principally due to the difference in the $k_{[\text{CII}]}$ (Herrera-Camus et al. 2021 assumed $k_{[\text{CII}]} = 1.5 M_{\odot}/L_{\odot}$, which corresponds to maximal excitation assuming solar metallicity).

In HZ4, the neutral atomic mass outflow rate is more than one order of magnitude higher than the ionized mass outflow rate. A similar behavior is observed both at low and high

<figure><img src="image_16.png" /><figcaption>Fig. 14. Comparison between the ionized and neutral mass outflow rate color-coded by redshift. The pink star marks the measurements for HZ4 from this work. The turquoise circles are the local ULIRGs studied in Fluetsch et al. (2021). The triangle, pentagon, diamond, and square display the results obtained for local and high-z AGNs by Perna et al. (2019), Cresci et al. (2023), Belli et al. (2024), D'Eugenio et al. (2024), respectively. The black dashed line represents the 1:1 relation.</figcaption></figure>

redshifts and for AGNs and SF galaxies (Herrera-Camus et al. 2020; Fluetsch et al. 2021; Avery et al. 2022; Cresci et al. 2023; D'Eugenio et al. 2024; Belli et al. 2024). Analyzing a sample of local ULIRGs, Fluetsch et al. (2021) shows that the amount of mass expelled by the galaxy in the ionized phase is usually negligible with respect to the neutral and the molecular phases, the latter of which accounts for the majority of the outflow mass. We compare our results with others from the literature in Fig. 14, where we show that both AGNs and SF galaxies at low and high redshift lie below the 1:1 relation between the ionized mass outflow rate and the neutral one.
## 8.3. Outflow strength and effect on the host galaxy

To understand whether these outflows are able to leave the galaxy, and hence remove gas from the system and enrich the CGM and the IGM, we estimate the escape velocity. We followed the same approach as in Carniani et al. (2024). We modeled the galaxy gravitational potential as the sum of an exponential disk profile potential for the baryonic component with mass $M_\star + M_{\text{gas}}$ (the dust mass, of ~$10^8 M_\odot$, is negligible; Pozzi et al. 2021) and a Navarro-Frenk-White (NFW; Navarro et al. 1996) profile for the potential of the dark matter halo.

We considered as $M_\star$ the stellar mass of the whole system reported in Table 3. For the gas mass, we convert the [CII] luminosity as $M_{\text{gas}} = \alpha_{[CII]} \times L_{[CII]}$. We adopt $\alpha_{[CII]} = 31 M_\odot/L_\odot$ (Zanella et al. 2018). To compute the mass of the dark matter halo, we adopted the stellar-to-halo mass relation from Moster et al. (2013), and we compute the concentration parameter for the NFW profile from the relation between the mass of the dark matter halo and the concentration estimated from Dutton & Macciò (2014).

We find that the velocity necessary for the gas to leave the galaxy and reach infinity is $v_{\text{esc}} \sim 800 \text{ km s}^{-1}$. This velocity is higher than the estimated outflow velocity of the ionized and neutral gas, hence these outflows are leaving the launching site, but eventually, the gas may fall back on the galaxy. We note that these calculations assume that the outflow moves ballistically (i.e., any other slowing mechanisms, such as the drag due to the gas encountered in the way, are not taken into account) and it is launched from the center of these idealized profiles.

To understand the impact that the outflow has on the removal of gas from the galaxy we estimate the mass loading factor $\eta = \dot{M}_{\text{out}}/\text{SFR}$. For the integrated values we obtain $\eta_{\text{ionized}} \sim 3\%$, while $\eta_{\text{neutral}} \sim 150\%$. The ionized value is smaller than what found for other high-z, but less massive galaxies ($\eta_{\text{ionized}} \sim 200\%$ for galaxies with $\log(M_\star/M_\odot) \sim 7.5-8.5$; Carniani et al. 2024). Both cosmological simulations and observations predict that the mass loading factor should decrease with increasing stellar mass for star-formation driven outflows, with mass loading factors reaching less than unity at a mass of around $\log(M_\star/M_\odot) \sim 10$, consistent with what we observe (Muratov et al. 2015; Nelson et al. 2019; Pandya et al. 2021; Carniani et al. 2024).

## 8.4. Outflow driver and launching region

To investigate the outflow driver we first computed the energetics of the outflow. We compute the kinetic and momentum rates as

$$ \dot{E}_{\text{out}} = \frac{1}{2} \dot{M}_{\text{out}} v_{\text{out}}^2 \quad (6) $$

$$ \dot{P}_{\text{out}} = \dot{M}_{\text{out}} v_{\text{out}}. \quad (7) $$

For the ionized outflow, we obtain $\dot{E}_{\text{out,H}\alpha} = 1.16 \times 10^{41} \text{ erg s}^{-1}$, $\dot{E}_{\text{out,[OIII]}} = 1.57 \times 10^{41} \text{ erg s}^{-1}$, and $\dot{P}_{\text{out,H}\alpha} = 5.9 \times 10^{33} \text{ g cm s}^{-2}$, $\dot{P}_{\text{out,[OIII]}} = 8.1 \times 10^{33} \text{ g cm s}^{-2}$ with uncertainties on the order of 0.5 dex. Even though we cannot completely rule out the possible presence of an AGN (see Sect. 3.1), the outflow velocity and the outflow energetics are small in comparison to what we expect from AGN-driven winds (Fiore et al. 2017). On the other hand, the outflow can be driven by the stellar feedback in a particular combination of winds from massive stars, and type 2 supernova (SN) explosions that release energy and momentum into the ISM (Veilleux et al. 2005) as explained below. The SFR threshold for driving an outflow is considered to be around $\Sigma_{\text{SFR}} \sim 1 M_\odot \text{ yr}^{-1} \text{ kpc}^{-2}$ (Newman et al. 2012; Davies et al. 2019). The star-forming obscured clump (HZ4-C)

<figure><img src="image_17.png" /><figcaption>Fig. 15. Ionized mass outflow rate map from the dust-corrected broad H$\alpha$ emission for each radial shell.</figcaption></figure>

has a SFR of ~20 $M_\odot \text{ yr}^{-1}$ in a region of radius ~0.9 kpc, resulting in a $\Sigma_{\text{SFR}} \sim 4 M_\odot \text{ yr}^{-1} \text{ kpc}^{-2}$ which is high enough to launch the ionized outflow.

If the outflow is driven by the stellar feedback from the burst in the central dusty region, we can compare the energy rate of the SN exploding in that region with the kinetic rate of the outflow. We can compute the energy rate driven by SN explosions as $\dot{E}_{\text{SNII}} = \epsilon E_{\text{SNII}} R$, where $\epsilon$ is the supernova efficiency, that is the efficiency of transferring kinetic energy to the ISM, $E_{\text{SNII}}$ is the total energy released in a SN explosion ($10^{51} \text{ erg}$), and $R$ is the supernova rate. Assuming a SFR of ~20 $M_\odot \text{ yr}^{-1}$, for a Kroupa IMF, the rate of SN explosions is related to the SFR as $R \sim 0.01 \text{ yr}^{-1} \times \text{SFR}/(M_\odot \text{ yr}^{-1})$ (Mo et al. 2010). If we assume a SN efficiency of $\epsilon \sim 0.1$, we obtain $\dot{E}_{\text{SNII}} = 7.3 \times 10^{41} \text{ erg s}^{-1}$ which is higher than the outflow kinetic rate ($\dot{E}_{\text{out}} = 1.6 \times 10^{41} \text{ erg s}^{-1}$). Following Heckman et al. (2015), the momentum rate injected by a starburst in the ISM is $\dot{P}_{\text{SB}} = 10^{33.7} \text{ SFR}/(M_\odot \text{ yr}^{-1}) \text{ g cm s}^{-2}$. In the case of the central clump, we obtain $\dot{P}_{\text{SB}} \sim 10^{35} \text{ g cm s}^{-2}$. Both the momentum and energy released by star formation are higher than those inferred for the outflow, thus the SF feedback in the central, obscured clump can drive the ionized outflow that we observe in the rest-frame optical emission lines.

Moreover, we observe that the outflow has higher $\Delta v$ and FWHM near the region of higher SFR (see Figs. 6 and 8), which corresponds to higher outflow velocity corrected for the unknown geometry. This is what we expect if that is the launching region and the outflow moves ballistically and slows down due to gravity, and possibly due to drag by gas encountered.

Since the position of the peak of the outflow flux and the narrow [OIII] flux do not match, we can infer that the outflow peak is not an artifact due to the higher S/N of the emission line, but rather a real higher flux of the outflow in that region. The interpretation of these data is that the dust-obscured clump launches the ionized outflow which moves toward the north in the foreground of the galaxy with respect to our line of sight. A schematic representation of this interpretation is presented in the top-left panel of Fig. 13.

## 8.5. Ionized outflow history

Speculating that the broad component is an outflow, and assuming that the ionized outflow is being launched by the central,
<figure><img src="image_18.png" /><figcaption>Fig. 16. Outflow velocity (top left), mass outflow rate (top right), momentum rate (bottom left), and kinetic rate (bottom right) in each radial shell as a function of radius from the central component (see Fig. 15). Blue and purple points indicate whether the outflow mass was calculated by assuming $A_V = 0$ or the $A_V$ from the outflow Balmer decrement in each ring, respectively. The half-blue, half-purple points have measured $A_V = 0$. The purple and blue shaded regions in top-right panel show the $1\sigma$ region around the mean outflow rate computed with $A_V \neq 0$ and $A_V = 0$, respectively.</figcaption></figure>

obscured, highly star-forming component, we extracted integrated spectra from 6 concentric annuli centered on that clump and spaced in radius by 0.1", which corresponds to 609 pc at $z = 5.54$. For each shell, we analyzed the spectra following the method reported in Sect. 3.1. Then we calculated the outflow properties as a function of radius from the launching region. We computed the mass outflow rate for each shell following Lutz et al. (2020) as

$$ \dot{M}_{\text{out,shell}} = \frac{v_{\text{out,shell}} M_{\text{out,shell}}}{\Delta R} \qquad (8) $$

where $\Delta R = 609$ pc, being the size of each shell, $v_{\text{out,shell}}$ is the velocity of the outflow computed as in Eq. (3), and $M_{\text{out,shell}}$ is the mass in the shell calculated as using Eq. (1). The outflow $H\alpha$ luminosity was corrected for dust obscuration using the Balmer decrement computed in each shell. In this case, the mass outflow rate represents the amount of gas passing through each projected shell. We show the mass outflow rate map in Fig. 15, in which each shell is color-coded by the mass outflow rate.

In Fig. 16 we show the outflow velocity (computed using Eq. (3)), the mass outflow rate, the momentum rate, and the kinetic rate for each radial shell. We report all the quantities,

assuming either no obscuration ($A_V = 0$, blue points) or the inferred obscuration based on the Balmer decrement in each shell ($A_V \neq 0$, purple points). In the first two rings and the last one, the ratio between the broad $H\alpha$ and $H\beta$ components is lower than the theoretical ratio of 2.86, so we assumed $A_V = 0$.

The observed outflow velocity is higher in the two rings closest to the launching region, then it decreases and stays almost constant up to 4 kpc. This can also be seen in the outflow maps from the spaxel-by-spaxel fitting in Figs. 6 and B.1. The mass outflow rate for each shell instead reaches its maximum between 1.4 and 3.5 kpc from the launching region. The momentum and kinetic rate radial profiles show a steady increase with the radius up to 3.5 kpc in both quantities in the case in which we corrected for the dust attenuation of the outflow.

If we assume that each shell has moved at a constant velocity $v_{\text{out}}$ in the plane of the sky since its launch, then the gas in each shell represents the outflow ejected at a time $\tau = v_{\text{out}}/l$, where $l$ is the intrinsic, not projected distance from the launching point to where it is observed now. With this assumption, we deduce that $\dot{M}_{\text{out}}$ was higher in the past since the mass outflow rate is higher at larger radii.
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Mass loading factor of the ionized gas defined as $\dot{M}_{\text{ion,out}}$/SFR as a function of time for a range of assumed inclination angles.</figcaption></figure>

Assuming that the outflow has been launched at a time $\tau$ we can compare the outflow observed at a radius $R$, with the SFR in the central clump at the time $\tau = v_{\text{out}} \times \cos(i)/R$, where $v_{\text{out}}$ is the de-projected outflow velocity, $R$ is the distance along the plane of the sky and $i$ is the assumed inclination angle of the outflow with respect to the plane perpendicular to the line of sight. We performed a SED fitting using Bagpipes with the same parameters as reported in Sect. 5, but with finer SFH bins in the last 20 Myr to properly characterize the variation of the SFR in the last few Myrs when the outflow was launched. In Fig. 17 we plot the mass loading factor $\eta$ as a function of lookback time, defined as $\eta = \dot{M}_{\text{ion,out}}(R)/\text{SFR}(\tau)$, where $\tau$ and $R$ are related by the equation reported above. For $v_{\text{out}}$ we assumed 450 km s$^{-1}$, which is the average outflow velocity across the different shells. Varying the inclination angle $i$ between 30 and 70 degrees we can observe that the outflows must have been stronger in the past, reaching values of $\eta \gtrsim 1$ for $i > 50$. This analysis confirms what already observed by Herrera-Camus et al. (2021) with the [CII] neutral outflow, the mass loading factor of the outflow is low considering the entire galaxy ($\eta_{\text{ion}} \sim 3$%), but becomes significant when considering only the central component.

## 9. Conclusions

In this work, we have presented the JWST NIRSpec/IFS observations of the galaxy HZ4 at $z = 5.544$. This galaxy is the highest redshift galaxy for which we have evidence of a neutral outflow traced by broad [CII] emission lines (Herrera-Camus et al. 2021). We have analyzed both the low- and high-spectral resolution datacubes from JWST NIRSpec IFS to study for the first time the rest-frame UV and optical spectrum of this source. Moreover, we combined our results with a re-analysis of the [CII] ALMA data to constrain the multiphase properties of this galaxy.

In particular, our main results are:

- The JWST high spatial resolution observations reveal that the galaxy is not a regular rotating disk as inferred from ALMA lower resolution observations, but rather it is formed by three different components. In particular, the different components have different properties, such as stellar ages and star formation histories.

- We observe a spatially resolved broad component that we speculate it is tracing an ionized outflow, which is co-spatial with the broad component seen in [CII] in our analysis and also in previous works. If interpreted as an outflow this

would be the highest redshift case in which we are able to probe the SF-driven multiphase outflow.

- Assuming that the broad component is tracing an ionized outflow, this component is extended ~4 kpc from the launching site that we speculate is the highly star-forming region. With this assumption we constructed the ionized outflow history by analyzing the outflow properties in radial shells from the launching site. We observe that the velocity of the outflow decreases with increasing radius, while the mass outflow rate increases. We also computed the mass loading factor history of the outflow and found that the outflow was stronger in the past, reaching values of $\eta \sim 1-2$ at 5–20 Myr, varying with the outflow inclination.

- If both the broad components in [OIII] and [CII] are outflows, we can estimate the multiphase outflow properties, finding that in this $z \sim 5.5$ source the mass of the SF-driven outflow is dominated by the colder gas phases, in agreement with other low-$z$ studies and a handful of studies up to $z \sim 3$.

The synergy between ALMA and JWST observations has allowed us to probe the multiphase ISM properties up to $z \sim 5.5$. These observations highlight also that high-resolution data are crucial for interpreting such complex systems. However, despite the large amount of observations for this source, it remains challenging to fully understand its geometry and pinpoint the origin of the broad components observed in both ALMA and JWST datasets. Further observations and larger samples are needed to better assess the role of mergers and outflows in shaping galaxy properties in the early Universe.

**Acknowledgements.** SC, GV, and SZ acknowledge support from the European Union (ERC, WINGS,101040227). SA and MP acknowledge grant PID2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/ 10.13039/501100011033). RM and FDE acknowledge support by the Science and Technology Facilities Council (STFC), from the ERC Advanced Grant 695671 "QUENCH", RM acknowledges support by funding from a research professorship from the Royal Society. FDE acknowledges support by the Science by the UKRI Frontier Research grant RISEandFALL. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Cavli Foundation through a Newton-Kavli Junior Fellowship. 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" AJB and GCJ 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). IL acknowledges support from grant PRIN-MUR 2020ACSP5K_002 financed by European Union - Next Generation EU. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01605.S, ADS/JAO.ALMA#2012.1.00523.S, ADS/JAO.ALMA#2017.1.00428.L. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

## References

- Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20

- Arribas, S., Perna, M., Rodríguez Del Pino, B., et al. 2024, A&A, 688, A146

- Asplund, M., Amarsi, A. M., & Grevesse, N. 2021, A&A, 653, A141

- Avery, C. R., Wuyts, S., Förster Schreiber, N. M., et al. 2022, MNRAS, 511, 4223

- Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5

- Baron, D., Netzer, H., Lutz, D., Davies, R. I., & Prochaska, J. X. 2024, ApJ, 968, 23

- Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

- Belli, S., Park, M., Davies, R. L., et al. 2024, Nature, 630, 54

- Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2

- Bischetti, M., Piconcelli, E., Feruglio, C., et al. 2019, A&A, 628, A118

- Böker, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82

- Bouché, N., Murphy, M. T., Kacprzak, G. G., et al. 2013, Science, 341, 50
Bregman, J. N. 1980, ApJ, 236, 577

Brook, C. B., Stinson, G., Gibson, B. K., et al. 2012, MNRAS, 419, 771

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Butler, K. M., van der Werf, P. P., Topkaras, T., et al. 2023, ApJ, 944, 134

Calabrò, A., Pentericci, L., Santini, P., et al. 2024, A&A, 690, A290

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, A&A, 677, A115

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Carniani, S., Maiolino, R., Amorin, R., et al. 2018, MNRAS, 478, 1170

Carniani, S., Venturi, G., Parlanti, E., et al. 2024, A&A, 685, A99

Cassata, P., Morselli, L., Faisst, A., et al. 2020, A&A, 643, A6

Chatzikos, M., Bianchi, S., Camilloni, F., et al. 2023, Rev. Mex. Astron. Astrofis., 59, 327

Chisholm, J., Tremonti, C., & Leitherer, C. 2018, MNRAS, 481, 1690

Cicone, C., Brusa, M., Ramos Almeida, C., et al. 2018, Nat. Astron., 2, 176

Coatman, L., Hewett, P. C., Banerji, M., et al. 2019, MNRAS, 486, 5335

Concas, A., Maiolino, R., Curti, M., et al. 2022, MNRAS, 513, 2535

Cresci, G., Tozzi, G., Perna, M., et al. 2023, A&A, 672, A128

Curti, M., Cresci, G., Mannucci, F., et al. 2017, MNRAS, 465, 1384

Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944

Curti, M., Maiolino, R., Curtis-Lake, E., et al. 2024, A&A, 684, A75

Dashyan, G., Silk, J., Mamon, G. A., Dubois, Y., & Hartwig, T. 2018, MNRAS, 473, 5698

Davies, R. L., Förster Schreiber, N. M., Übler, H., et al. 2019, ApJ, 873, 122

De Looze, I., Baes, M., Bendo, G. J., Cortese, L., & Fritz, J. 2011, MNRAS, 416, 2712

Dekel, A., & Silk, J. 1986, ApJ, 303, 39

D'Eugenio, F., Pérez-González, P. G., Maiolino, R., et al. 2024, Nat. Astron., 8, 1443

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32

Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359

Fabian, A. C. 2012, ARA&A, 50, 455

Faisst, A. L., Scherer, D., Lemaux, B. C., et al. 2020a, ApJS, 247, 61

Faisst, A. L., Fudamoto, Y., Oesch, P. A., et al. 2020b, MNRAS, 498, 4192

Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143

Fluetsch, A., Maiolino, R., Carniani, S., et al. 2019, MNRAS, 483, 4586

Fluetsch, A., Maiolino, R., Carniani, S., et al. 2021, MNRAS, 505, 5753

Förster Schreiber, N. M., Übler, H., Davies, R. L., et al. 2019, ApJ, 875, 21

Fudamoto, Y., Inoue, A. K., & Sugahara, Y. 2023, MNRAS, 521, 2962

Gallerani, S., Pallottini, A., Feruglio, C., et al. 2018, MNRAS, 473, 1909

Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101

Genzel, R., Jolly, J. B., Liu, D., et al. 2023, ApJ, 957, 48

Ginolfi, M., Jones, G. C., Béthermin, M., et al. 2020, A&A, 633, A90

Goerdt, T., & Ceverino, D. 2015, MNRAS, 450, 3359

Harikane, Y., Zhang, Y., Nakajima, K., et al. 2023, ApJ, 959, 39

Heckman, T. M., Alexandroff, R. M., Borthakur, S., Overzier, R., & Leitherer, C. 2015, ApJ, 809, 147

Herrera-Camus, R., Bolatto, A. D., Wolfire, M. G., et al. 2015, ApJ, 800, 1

Herrera-Camus, R., Tacconi, L., Genzel, R., et al. 2019, ApJ, 871, 37

Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2020, A&A, 633, L4

Herrera-Camus, R., Förster Schreiber, N., Genzel, R., et al. 2021, A&A, 649, A31

Herrera-Camus, R., Förster Schreiber, N. M., Price, S. H., et al. 2022, A&A, 665, L8

Ho, I.-T., Kewley, L. J., Dopita, M. A., et al. 2014, MNRAS, 444, 3894

Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522

Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80

Janssen, A. W., Christopher, N., Sturm, E., et al. 2016, ApJ, 822, 43

Jones, G. C., Maiolino, R., Caselli, P., & Carniani, S. 2019, A&A, 632, L7

Jones, G. C., Vergani, D., Romano, M., et al. 2021, MNRAS, 507, 3540

Jones, G. C., Übler, H., Perna, M., et al. 2024a, A&A, 682, A122

Jones, G. C., Bunker, A. J., Telikova, K., et al. 2024b, MNRAS, submitted [arXiv:2405.12955]

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121

Koudmani, S., Sijacki, D., Bourne, M. A., & Smith, M. C. 2019, MNRAS, 484, 2047

Koudmani, S., Henden, N. A., & Sijacki, D. 2021, MNRAS, 503, 3568

Lamperti, I., Arribas, S., Perna, M., et al. 2024, A&A, 691, A153

Laseter, I. H., Maseda, M. V., Curti, M., et al. 2024, A&A, 681, A70

Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1

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

Leroy, A. K., Walter, F., Martini, P., et al. 2015, ApJ, 814, 83

Li, M., Bryan, G. L., & Ostriker, J. P. 2017, ApJ, 841, 101

Liddle, A. R. 2007, MNRAS, 377, L74

Lutz, D., Sturm, E., Janssen, A., et al. 2020, A&A, 633, A134

Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2024, A&A, 691, A145

Mallery, R. P., Mobasher, B., Capak, P., et al. 2012, ApJ, 760, 128

Mazzolari, G., Übler, H., Maiolino, R., et al. 2024, A&A, 691, A345

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, Astronomical Society of the Pacific Conference Series, 376, 127

Medling, A. M., Cortese, L., Croom, S. M., et al. 2018, MNRAS, 475, 5194

Mitsuhashi, I., Tadaki, K.-I., Ikeda, R., et al. 2024, A&A, 690, A197

Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge, UK: Cambridge University Press)

Monreal-Ibero, A., Arribas, S., Colina, L., et al. 2010, A&A, 517, A28

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Muratov, A. L., Kereš, D., Faucher-Giguère, C.-A., et al. 2015, MNRAS, 454, 2691

Naab, T., & Ostriker, J. P. 2017, ARA&A, 55, 59

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

Nakajima, K., Ouchi, M., Isobe, Y., et al. 2023, ApJS, 269, 33

Nakazato, Y., Ceverino, D., & Yoshida, N. 2024, ApJ, 975, 238

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Nelson, B. D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Newman, S. F., Genzel, R., Förster-Schreiber, N. M., et al. 2012, ApJ, 761, 43

Oppenheimer, B. D., & Davé, R. 2006, MNRAS, 373, 1265

Oppenheimer, B. D., Davé, R., Kereš, D., et al. 2010, MNRAS, 406, 2325

Osterbrock, D. E., & Ferland, G. J. 2006, in Astrophysics of gaseous nebulae and active galactic nuclei, eds. D. E. Osterbrock, & G. J. Ferland, 2nd edn. (Sausalito, CA: University Science Books)

Pandya, V., Fielding, D. B., Anglés-Alcázar, D., et al. 2021, MNRAS, 508, 2979

Parlanti, E., Carniani, S., Pallottini, A., et al. 2023, A&A, 673, A153

Pearson, W. J., Wang, L., Alpaslan, M., et al. 2019, A&A, 631, A51

Perna, M., Cresci, G., Brusa, M., et al. 2019, A&A, 623, A171

Perna, M., Arribas, S., Marshall, M., et al. 2023, A&A, 679, A89

Pettini, M., & Pagel, B. E. J. 2004, MNRAS, 348, L59

Pineda, J. L., Langer, W. D., & Goldsmith, P. F. 2014, A&A, 570, A121

Planck Collaboration XIII. 2016, A&A, 594, A13

Posses, A., Aravena, M., González-López, J., et al. 2024, A&A, submitted

Pozzi, F., Calura, F., Fudamoto, Y., et al. 2021, A&A, 653, A84

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

Reines, A. E., Condon, J. J., Darling, J., & Greene, J. E. 2020, ApJ, 888, 36

Rich, J. A., Kewley, L. J., & Dopita, M. A. 2011, ApJ, 734, 87

Rich, J. A., Kewley, L. J., & Dopita, M. A. 2015, ApJS, 221, 28

Rodríguez Del Pino, B., Perna, M., Arribas, S., et al. 2024, A&A, 684, A187

Rubin, K. H. R., Prochaska, J. X., Koo, D. C., & Phillips, A. C. 2012, ApJ, 747, L26

Rupke, D. S. N., & Veilleux, S. 2013, ApJ, 768, 75

Sanders, R. L., Shapley, A. E., Kriek, M., et al. 2016, ApJ, 816, 23

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

Sartori, L. F., Schawinski, K., Treister, E., et al. 2015, MNRAS, 454, 3722

Scholtz, J., Maiolino, R., D'Eugenio, F., et al. 2023, A&A, submitted

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Shih, H.-Y., & Rupke, D. S. N. 2010, ApJ, 724, 1430

Silk, J. 2017, ApJ, 839, L13

Solimano, M., González-López, J., Aravena, M., et al. 2024, A&A, 689, A145

Solimano, M., González-López, J., Aravena, M., et al. 2025, A&A, 693, A70

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Spilker, J. S., Aravena, M., Phadke, K. A., et al. 2020, ApJ, 905, 86

Tripodi, R., Lelli, F., Feruglio, C., et al. 2023, A&A, 671, A44

Übler, H., Naab, T., Oser, L., et al. 2014, MNRAS, 443, 2092

Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023, A&A, 677, A145

Übler, H., D'Eugenio, F., Perna, M., et al. 2024, MNRAS, 533, 4287

Veilleux, S., & Osterbrock, D. E. 1987, ApJS, 63, 295

Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769

Venturi, G., Carniani, S., Parlanti, E., et al. 2024, A&A, 691, A19

Villanueva, V., Herrera-Camus, R., González-López, J., et al. 2024, A&A, 691, A133

Villar Martín, M., Perna, M., Humphrey, A., et al. 2020, A&A, 634, A116

Yuan, Y., Krumholz, M. R., & Martin, C. L. 2023, MNRAS, 518, 4084

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976

Zhang, Y., Ouchi, M., Nakajima, K., et al. 2024, ApJ, 970, 19
<figure><img src="image_20.png" /><figcaption>**Fig. A.1.** Top panel: R100 cube integrated between 1.77 and 1.84 µm. Bottom panel: In black we plot the spectrum extracted from the area highlighted as the magenta circle in the panel above. In red we show the best-fitting result from the Bagpipes fitting.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. B.1.** Outflow velocity map defined as $v_{\text{out}} = |\Delta v_{\text{narrow,broad}}| + 2\sigma_{\text{out}}$ (see also Sect. 8). We show as black dashed contours the $\Sigma_{\text{SFR}}$ map contours (see also Fig. 8).</figcaption></figure>

## Appendix A: Other sources in the FOV

While analyzing the R100 cube we detected another source in the field of view, which is distant ~0.8" from the northern bright clump of HZ4. The source is located at RA = 09h58m28.55s, Dec = 02d03m06.69s, and it is undetected in HST or ALMA. The continuum emission of the source is negligible but is visible only when integrating between 1.77 and 1.84 µm as shown in the top panel of Fig. A.1. In the bottom panel of Fig. A.1 we show the spectrum of the source extracted from a circular aperture of radius 0.2" (4 spaxels) highlighted as the pink circle in the map. The spectra show two bright emission lines, that we identify as the [OIII]λλ5007,4959Å, Hβ complex and the Hα and the [NII]λλ6548,6584Å complex, which are unresolved in the prism spectrum. We performed a fit of the emission lines by using the same procedure as reported in 3.1 and implying a single Gaussian component for each emission line. This places the galaxy at a redshift of $z = 2.607 \pm 0.003$.

We also estimate the galaxy properties by fitting the integrated spectrum with Bagpipes.

The best-fit Bagpipes spectrum is shown in the bottom panel of Fig. A.1 as the red line. The galaxy has a stellar mass of $\log (M_{\star}/M_{\odot}) = 7.2 \pm 0.1$.

## Appendix B: Outflow velocity

In this Sect. we show the outflow velocity map computed as Eq. 3 for each spaxel and we compare it with the $\Sigma_{\text{SFR}}$ map (see also Fig. 8).

## Appendix C: Outflow history

In this Sect. we show the individual fit of each spectrum extracted from each ring for the analysis of the outflow history in Fig. C.1. Ring 1 represents the closer one to the launching site, Ring 7 is the farther. We show the [OIII], Hβ complex.
<figure><img src="image_22.png" /><figcaption>**Fig. C.1.** Zoom in on the [OIII] and Hβ complex of each ring. Red and blue dashed lines are the narrow and outflow best-fit models. Please note that the displayed range on the *y* axis is the same for all the panels.</figcaption></figure>