© The Authors 2025

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_⊙) 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 feedback, 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.

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⁴ K), which is usually detected as blueshifted broad wings in the brightest optical emission lines, such as [O III]λ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 [C II]λ158 μm¹ 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₂0). 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 ≤ z ≤ 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, [C II] (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 [C II] emission of samples of galaxies at z ∼ 4 − 5 (Gallerani et al. 2018; Ginolfi et al. 2020). So far, broad [C II] 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 [C II] emission line (Herrera-Camus et al. 2021). HZ4 was identified as a Lyman Break Galaxy (LBG) at z ∼ 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 [C II] 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 [C II] 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 [C II] 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₀ = 67.7 km s⁻¹ Mpc⁻¹, Ω_m = 0.307, and Ω_Λ = 0.691, giving 1″ = 6.09 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² 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 ∼ 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 ∼ 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 jwst1097.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_detector1” that accounts for detector level corrections, we also corrected 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 [C II] 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 [O III]λ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 [O III] integrated emission along with the contours from the HST F160W image and the [C II] integrated emission after the astrometric correction.

Fig. 1. Integrated [O III]λ5007 Å flux map of HZ4 from the JWST NIRSpec/IFS R2700 cube. We show the 5, 25, and 75σ of [O III] as green contours. Red contours illustrate the [C II] 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.

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 [O III]λ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).

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 [O III]λ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 [O III] complex is shown in the left panel, and the Hα, [N II], and [S II]λλ 6717,6731 complex is shown in the right panel. The gray shaded areas mark the uncertainties on the data.

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β, [O III]λλ5007,4959 Å, [N II]λλ6548,6584 Å, [S II]λλ 6717,6731 Å and He Iλ 5875 Å.

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 [S II]λλ 6717,6731 and He Iλ 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 < σ < 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 [O III] and [N II] doublets have been furthermore constrained to atomic physics motivated ratios (see Osterbrock & Ferland 2006) of [O III]λ5007/[O III]λ4959 = 2.98 and [N II]λ6584/[N II]λ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³ (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 [O III] complex on the left, and the Hα, [N II], and [S II] 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 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 ± 0.00002 for the galaxy, which is in agreement with the redshift measured from the [C II] 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_broad − v_narrow = −24 ± 3 km s⁻¹. The FWHM of the broad component is comparable with the [C II] 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α/F_Hβ = 2.86 (for a temperature T = 10⁴ K and a density n_e = 10² 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 ± 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.3}^{+0.9}$.

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 $\mathrm{SFR}={77}_{-16}^{+19}$ M_⊙ yr⁻¹, 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_⊙ yr⁻¹.

By exploiting the ratio between the fluxes of [S II]λ6717 Å and [S II]λ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 [S II]λλ 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 individual 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 μ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 [O II]λλ3727,3729, [Ne III]λ 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.

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σ error. In blue we show the best-fit emission lines, and in red the best-fit model.

In addition, we set the flux ratio of the [O III] and the [N II] doublets according to atomic physics prescriptions (see Sect. 3.1). We use only one Gaussian component to fit the [S II]λλ 6717,6731 and the [O II]λλ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 [N II] 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 [N II] emission lines, which are blended together, however the summed flux agree within the uncertainties. From the expected line ratio of F_Hγ/F_Hβ = 0.466, we estimated a dust extinction of $A_V=0.7_{-0.7}^{+1.4}$ assuming that the [O III]λ4363 Å 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 = [O III]λ5007 Å/Hβ, R2 = [O II]λλ3727,3729/Hβ, $\widehat{\mathrm{R}}=0.47\mathrm{R}2+0.88\mathrm{R}3$, O32 = [O III]/[O II]λλ3727,3729, N2 = [N II]/Hα, and O3N2 = ([O III]/Hβ)/([N II]/Hα) to estimate the abundance of oxygen log(O/H).

The [O II]λλ3727,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 [O II]λλ3727,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 ± 0.3). We estimate a value of 12 + log(O/H) = 8.34 ± 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)_narrow = 8.3 ± 0.1 for the entire galaxy and 12 + log(O/H)_broad = 8.4 ± 0.1 for the broad component. Both values are consistent with the value we infer by adding the calibrators that exploit the [O II]λλ3727,3729 lines.

The value that we obtain is in agreement with metallicity obtained for other high-z galaxies with similar stellar mass ($log(M_{\star }/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)_narrow = 8.69 (Asplund et al. 2021), we obtain a metallicity of Z_narrow = 0.4 ± 0.1 Z_⊙ and Z_broad = 0.5 ± 0.1 Z_⊙ 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 [O III] 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σ of the map at v = 0 km s⁻¹ to highlight the different morphologies of the [O III] emission at different velocities.

Fig. 4. Channel maps from the R2700 cube targeting the [O III]λ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 [C II] emission peak (see Fig. 1). The white circles represent the position of the components of HZ4.

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⁻¹ 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⁻¹). 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, 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⁻¹ slightly northern than the 75σ 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 [O III] 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 ∼ 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⁻¹, is now evident. This clump is very faint in [O III] 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⁻¹ is coincident with the location of the peak of the [C II] emission (see Fig. 1) and the peak of the SFR discussed in Sect. 4.

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 [O III] emission at 5, 25, and 75σ (see Fig. 1). The black cross shows the location of the [C II] 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).

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 [O III] emission, a central component (HZ4-Central) which is bright in [C II] emission, and a redshifted southern component (HZ4-South), which is faint in [O III] 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₁ − BIC₂ > 10, which represents positive evidence toward the model with two components (Kass & Raftery 1995). In the other cases, 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 [C II] FIR emission line, as already seen in the channel maps in Fig. 4. The narrow [O III] 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 [O III] map shows an extended tail, in the southwest direction. The tail in [O III] shows the presence of a second faint component (HZ4-Central), which is instead more prominent and visible in the Hα map, and is also where the [C II] peak. This suggests that the second component is likely tracing a region with higher dust obscuration, as we discuss later in this section.

Fig. 6. Results for the spaxel-by-spaxel fitting of the R2700 cube. In the first row, we show from left to right the [O III]λ5007 Å flux, velocity, and velocity dispersion maps of the narrow component, respectively. In the second row, we show the same but for the [O III] broad component. In the [O III] broad flux map, we show the 5, 25, and 75σ contours of the narrow [O III] 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 [N II]λ6584 Å narrow components, respectively. The plus signs show the location of the three components identified in Figs. 4 and 5.

The narrow component velocity field (top-central panel of Fig. 6) shows a different picture from what was seen in ALMA observations. Whereas the [C II] 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 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_out = |Δv_{narrow, broad}|+2σ_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 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 [N II] 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 [N II] emission, one associated with the northern component, and one associated with the central one. Differently from the Hα and [O III] the brightest [N II] 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 [N II]/Hα and [O III]/Hβ while the latter uses the [S II]λλ 6717,6731/Hα and [O III]/Hβ. 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 ∼ 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).

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 ∼ 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 cm−3 ≤ n ≤ 1000 cm−3, solar abundances, and magnetic field between 10−4 and 10 μ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.

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 cm⁻³, velocities ranging from 200 to 1000 km/s and magnetic fields from 10⁻⁴ to 10 μ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 [N II]/Hα line ratio.

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 He II 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, subdominant AGN.

We display the maps of the log([N II]/Hα) and log([O III]/Hβ) 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([N II]/Hα) is evident, with the ratio increasing from north to south, while the values of log([O III]/Hβ) are almost constant. The broad component, on the other hand, shows a gradient also in the log([O III]/Hβ) values, varying from 0.2 to 0.7, while the log([N II]/Hα) 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 [N II]/Hα 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 [N II] emission does not allow us to constrain the metallicity.

Having excluded the presence of an AGN, we can also estimate the SFR surface density (Σ_SFR) at a spaxel level by exploiting the Hα dust-corrected luminosity. First, we estimated the dust attenuation at a spaxel level thanks to our detection of Hα and Hβ emission lines. By assuming case B recombination (T_e = 10 000 K and n_e = 100 cm³), the theoretical ratio between the Hα and Hβ flux is F_Hα/F_Hβ = 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 ≤ A_V ≤ 2.5. In the central panel of Fig. 8 we show the map of dust attenuation across the galaxy.

Fig. 8. Spatially resolved properties of the narrow component. From left to right: gas-phase metallicity (12 + log(O/H)), dust extinction (AV), and dust-corrected star formation rate surface density (Σ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 [O III] flux from Fig. 6.

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 ∼ 2.5 which is cospatial with the central component.

By correcting the Hα flux for the dust extinction and using the calibration from Kennicutt & Evans (2012) to get the SFR from Hα 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 [C II] flux contours shown in Fig. 1, it can be seen that this region is coincident with the peak of the [C II] 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 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 [O III] 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_{\star }/M_{\odot })=9.{75}_{-0.10}^{+0.05}$ which is smaller than the previous estimate from the literature of log(M_⋆/M_⊙) = 10.15 ± 0.15 of Faisst et al. (2020a), but in agreement with the one from Capak et al. (2015) of log(M_⋆/M_⊙) = 9.67 ± 0.21.

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.

The SFH for the whole HZ4 system shows an almost constant SFR of a few M_⊙ yr⁻¹ up to the last 20 Myr, where we see a rapid increase reaching ∼100 M_⊙ yr⁻¹ 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α 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, 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_⊙ yr⁻¹ at the time of the observation. From the SED fitting we obtain a mass-weighted age of the southern galaxy of ${322}_{-108}^{+88}$ 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 observations of the [C II] emission line and the dust continuum. The ALMA [C II] 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 discussed 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 weighting, which has a higher S/N at the expense of angular resolution with respect to the Briggs datacube. We fitted each spaxel having a S/N of the [C II] 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₁ − BIC₂ > 5 and the S/N of both the narrow 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 procedure 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. [C II] spaxel-by-spaxel results

Here we present the results of the spaxel-by-spaxel analysis of the ALMA data targeting the [C II] emission line. In Fig. 10 we show the best-fit results for the [C II] 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 velocity 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.

Fig. 10. Example of best-fit results for the fitting of the [C II] 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.

In Fig. 11 we show the results for the spaxel-by-spaxel fitting of the [C II] data. In the first two rows we show moment maps of the narrow and broad components from the natural weighting 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 weighting 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 [O III] broad component, it is redshifted with respect to the host galaxy (Δv ∼ 35 km s⁻¹), but the projection is co-spatial with the [O III] broad (shown as the red contours). The broad [C II] 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 component show an elongated [C II] emission due to the presence of two separate spatial components, one peaking where the [O III] peaks (the north component), one at the location of the central component, where the integrated [C II] 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).

Fig. 11. [C II] 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σ threshold. In the left panels, we show the ALMA beam as a green ellipse and the narrow and broad [O III] flux from Fig. 6 as lime and red contours, respectively.

In the bottom panels of Fig. 11 we show the maps of the narrow 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 [C II] flux is in the location of the central obscured component. This is in agreement with the [C II] being a tracer of star formation 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 continuum 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 follows the elongation of the source, already seen in [O III] and [C II] 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 emission 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 μm = 134 ± 33 μJy, comparable 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.

Fig. 12. AV 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 [O III] channel maps and in optical continuum emission. The purple ellipse represents the beam size of the dust continuum emission.

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 optical rest-frame emission observations, was estimated by summing the unobscured SFR estimated from UV continuum emission, 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_UV = 29 ± 4 M_⊙ yr⁻¹ exploiting the UV HST observations, and through SED fitting they estimate a ${\mathrm{SFR}}_{\mathrm{IR}}={77}_{-69}^{+104}$ M_⊙ yr⁻¹. Despite three FIR detections, the SFR_IR 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 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α luminosity, corrected for the dust obscuration according to the Balmer decrement. Assuming that the dust-corrected Hα luminosity is able to recover the entire intrinsic SFR ongoing in the system, we obtain that the total SFR of the system is ${\mathrm{SFR}}_{\mathrm{H}\alpha }={77}_{-16}^{+19}$ M_⊙ yr⁻¹. Of this, 29 ± 4 M_⊙ yr⁻¹ is unobscured SF, and the remaining is obscured SF of ${\mathrm{SFR}}_{\mathrm{obscured}}={48}_{-20}^{+23}$ M_⊙ yr⁻¹. By using the relation between the obscured SFR and L_IR by Kennicutt & Evans (2012) log(SFR/[M_⊙ yr⁻¹]) = log(L_IR/[erg s⁻¹]) − 43.41, we infer an IR luminosity of log(L_IR/L_⊙) = 11.5 ± 0.2. The IR luminosity obtained is comparable within the errors to the value obtained in other works from the SED fitting (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α trace SF over a similar timescale, a dust emissivity index of β = 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 ± 12 K and a dust mass of log(M_D/M_⊙) = 6.62 ± 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, [O III]λ5007 Å and other rest-frame optical emission lines, [C II] 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 [C II] 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 [O III] bright northern clump), triggered by the merger event. The ionized and molecular outflows are almost co-spatial, but the [C II] 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 Σ_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 [C II] 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 [C II] 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 [C II] 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 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 [C II] 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).

• [C II] inflow or merger-induced flows and [O III] 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 [O III] emission traces outflows, the [C II] 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 [C II] 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 [C II], and thus cannot be pristine gas, hence the gas could 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_{\mathrm{vir}}=\sqrt{G{M}_{\mathrm{DM}}/r_{\mathrm{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 σ ∼ 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 [O I]λ 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([O I]λ 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.

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 [C II] 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 [C II] 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 [C II] 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 [C II] 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.

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

$$\begin{array}{c}\hfill M_{\mathrm{ion},\mathrm{out}}=3.2\times {10}^5\left(\frac{{\mathrm{L}}_{\mathrm{H}\alpha ,\mathrm{outflow}}}{{10}^{40}{\mathrm{erg}\mathrm{s}}^{-1}}\right)\left(\frac{100{\mathrm{cm}}^{-3}}{n_e}\right)M_{\odot },\end{array}$$(1)

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

$$\begin{array}{c}\hfill M_{\mathrm{ion},\mathrm{out}}=0.8\times {10}^8\left(\frac{{\mathrm{L}}_{[\mathrm{OIII}],\mathrm{outflow}}}{{10}^{44}{\mathrm{erg}\mathrm{s}}^{-1}}\right)\left(\frac{500{\mathrm{cm}}^{-3}}{n_e}\right)\left(\frac{Z_{\odot }}{Z}\right)M_{\odot },\end{array}$$(2)

where L_{Hα, outflow} and L_{[OIII],outflow} are the luminosity of the Hα and [O III] 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 cm⁻³, which is consistent with the value derived from the [S II]λλ 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_⊙ (see Sect. 3.2). By using the value of Hα and [O III] fluxes reported in Table 1, we corrected for the outflow extinction of $A_V=0.3_{-0.3}^{+0.9}$ and we estimate an ionized outflow mass of $M_{\mathrm{ion},\mathrm{out}}=(1.8_{-0.5}^{+2.0})\times {10}^7{M}_{\odot }$ and $M_{\mathrm{ion},\mathrm{out}}=(2.5_{-0.8}^{+4.5})\times {10}^7{M}_{\odot }$ for Hα and [O III], respectively.

For the outflow velocity v_out we adopted the following relation from Rupke & Veilleux (2013) and Fiore et al. (2017):

$$\begin{array}{c}\hfill v_{\mathrm{out}}=|\mathrm{\Delta }{v}_{\mathrm{narrow},\mathrm{broad}}|+2{\sigma }_{\mathrm{out}}=389\pm 12{\mathrm{km}\mathrm{s}}^{-1},\end{array}$$(3)

where Δv_{narrow, broad} is the velocity shift between the peak of the broad and the narrow components and σ_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)

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}=\frac{v_{\mathrm{out}}{M}_{\mathrm{out}}}{R_{\mathrm{out}}},\end{array}$$(4)

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

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:

$$\begin{array}{c}\hfill M_H=k_{[\mathrm{CII}]}(T,n_H,Z)\times L_{[\mathrm{CII}]}\end{array}$$(5)

where k_[CII] is a conversion factor that depends on the temperature (T), the neutral gas density (n_H), and the gas metallicity (Z), and L_[CII] is the luminosity of the observed [C II] outflow. Since the outflow temperature and density are unknown we conservatively compute the value of k_[CII] assuming T = 10⁴ K and n_H = 10⁴ cm⁻³, 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 [C II] emission line.

Assuming an oxygen abundance of 12 + log(O/H) = 8.34 (Z = 0.45 Z_⊙, see Sect. 3.2), the case for maximal excitation gives k_[CII] = 4.88 M_⊙/L_⊙ (Herrera-Camus et al. 2021). The minimum neutral mass in the outflow is then M_H, out = 2.69 × 10⁸ M_⊙. We compute the flux-weighted radius of the outflow from the map shown in Fig. 11, which gives R_H, out ∼ 1.3 kpc, and the outflow velocity computed as in Eq. (3) is v_H, out ∼ 260 km s⁻¹. Finally, we compute a lower limit on the mass outflow rate of Ṁ_H,out ∼ 121 M_⊙ yr⁻¹. The value we obtain is higher than the lower limit reported in Herrera-Camus et al. (2021) of ∼34 M_⊙ yr⁻¹, but this is principally due to the difference in the k_[CII] (Herrera-Camus et al. 2021 assumed k_[CII] = 1.5 M_⊙/L_⊙, 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 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.

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.

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_⋆ + M_gas (the dust mass, of ∼10⁸ M_⊙, 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_⋆ the stellar mass of the whole system reported in Table 3. For the gas mass, we convert the [C II] luminosity as M_gas = α_[CII] × L_[CII]. We adopt α_[CII] = 31 M_⊙/L_⊙ (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_esc ∼ 800 km s⁻¹. 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 η = Ṁ_out/SFR. For the integrated values we obtain η_ionized ∼ 3%, while η_neutral ∼ 150%. The ionized value is smaller than what found for other high-z, but less massive galaxies (η_ionized ∼ 200% for galaxies with log(M_⋆/M_⊙)∼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_⋆/M_⊙)∼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

$$\begin{array}{cc}\hfill {\dot{E}}_{\mathrm{out}}& =\frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\dot{P}}_{\mathrm{out}}& ={\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}.\hfill \end{array}$$(7)

For the ionized outflow, we obtain Ė_out,Hα = 1.16 × 10⁴¹ erg s⁻¹, Ė_out[OIII] = 1.57 × 10⁴¹ erg s⁻¹, and Ṗ_out,Hα = 5.9 × 10³³ g cm s⁻², Ṗ_out,[OIII] = 8.1 × 10³³ g cm s⁻² 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 Σ_SFR ∼ 1 M_⊙ yr⁻¹ kpc⁻² (Newman et al. 2012; Davies et al. 2019). The star-forming obscured clump (HZ4-C) has a SFR of ∼20 M_⊙ yr⁻¹ in a region of radius ∼0.9 kpc, resulting in a Σ_SFR ∼ 4 M_⊙ yr⁻¹ kpc⁻² 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 Ė_SNII = ϵE_SNIIR, where ϵ is the supernova efficiency, that is the efficiency of transferring kinetic energy to the ISM, E_SNII is the total energy released in a SN explosion (10⁵¹ erg), and R is the supernova rate. Assuming a SFR of ∼20 M_⊙ yr⁻¹, for a Kroupa IMF, the rate of SN explosions is related to the SFR as R ∼ 0.01 yr⁻¹ × SFR/(M_⊙ yr⁻¹) (Mo et al. 2010). If we assume a SN efficiency of ϵ ∼ 0.1, we obtain Ė_SNII = 7.3 × 10⁴¹ erg s⁻¹ which is higher than the outflow kinetic rate (Ė_out = 1.6 × 10⁴¹ erg s⁻¹). Following Heckman et al. (2015), the momentum rate injected by a starburst in the ISM is Ṗ_SB = 10^33.7 SFR/(M_⊙ yr⁻¹) g cm s⁻². In the case of the central clump, we obtain Ṗ_SB ∼ 10³⁵ g cm s⁻². 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 Δ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 [O III] 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, 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

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out},\mathrm{shell}}=\frac{v_{\mathrm{out},\mathrm{shell}}{M}_{\mathrm{out},\mathrm{shell}}}{\mathrm{\Delta }R}\end{array}$$(8)

where ΔR = 609 pc, being the size of each shell, v_out, shell is the velocity of the outflow computed as in Eq. (3), and M_out, shell is the mass in the shell calculated as using Eq. (1). The outflow Hα 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 ≠ 0, purple points). In the first two rings and the last one, the ratio between the broad Hα and Hβ components is lower than the theoretical ratio of 2.86, so we assumed A_V = 0.

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 AV = 0 or the AV from the outflow Balmer decrement in each ring, respectively. The half-blue, half-purple points have measured AV = 0. The purple and blue shaded regions in top-right panel show the 1σ region around the mean outflow rate computed with AV ≠ 0 and AV = 0, respectively.

Fig. 15. Ionized mass outflow rate map from the dust-corrected broad Hα emission for each radial shell.

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.

Fig. 17. Mass loading factor of the ionized gas defined as Ṁion,out/SFR as a function of time for a range of assumed inclination angles.

If we assume that each shell has moved at a constant velocity v_out in the plane of the sky since its launch, then the gas in each shell represents the outflow ejected at a time τ = v_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 Ṁ_out was higher in the past since the mass outflow rate is higher at larger radii.

Assuming that the outflow has been launched at a time τ we can compare the outflow observed at a radius R, with the SFR in the central clump at the time τ = v_out × cos(i)/R, where v_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 η as a function of lookback time, defined as η = Ṁ_ion,out(R)/SFR(τ), where τ and R are related by the equation reported above. For v_out we assumed 450 km s⁻¹, 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 η ≳ 1 for i > 50. This analysis confirms what already observed by Herrera-Camus et al. (2021) with the [C II] neutral outflow, the mass loading factor of the outflow is low considering the entire galaxy (η_ion ∼ 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 [C II] 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 [C II] 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 [C II] 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 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 η ∼ 1 − 2 at 5–20 Myr, varying with the outflow inclination.

• If both the broad components in [O III] and [C II] are outflows, we can estimate the multiphase outflow properties, finding that in this z ∼ 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 ∼ 3.

The synergy between ALMA and JWST observations has allowed us to probe the multiphase ISM properties up to z ∼ 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.

Acknowledgments

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 Kavli 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

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 [O III]λλ5007,4959 Å, Hβ complex and the Hα and the [N II]λλ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 ± 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_⋆/M_⊙) = 7.2 ± 0.1.

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.

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 Σ_SFR map (see also Fig. 8).

Fig. B.1. Outflow velocity map defined as vout = |Δvnarrow, broad|+2σout (see also Sect. 8). We show as black dashed contours the ΣSFR map contours (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 [O III], Hβ complex.

Fig. C.1. Zoom in on the [O III] 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.

All Tables

All Figures

Fig. 1. Integrated [O III]λ5007 Å flux map of HZ4 from the JWST NIRSpec/IFS R2700 cube. We show the 5, 25, and 75σ of [O III] as green contours. Red contours illustrate the [C II] 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.

In the text

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 [O III]λ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 [O III] complex is shown in the left panel, and the Hα, [N II], and [S II]λλ 6717,6731 complex is shown in the right panel. The gray shaded areas mark the uncertainties on the data.

In the text

	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σ error. In blue we show the best-fit emission lines, and in red the best-fit model.
In the text

Fig. 4. Channel maps from the R2700 cube targeting the [O III]λ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 [C II] emission peak (see Fig. 1). The white circles represent the position of the components of HZ4.

In the text

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 [O III] emission at 5, 25, and 75σ (see Fig. 1). The black cross shows the location of the [C II] 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).

In the text

Fig. 6. Results for the spaxel-by-spaxel fitting of the R2700 cube. In the first row, we show from left to right the [O III]λ5007 Å flux, velocity, and velocity dispersion maps of the narrow component, respectively. In the second row, we show the same but for the [O III] broad component. In the [O III] broad flux map, we show the 5, 25, and 75σ contours of the narrow [O III] 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 [N II]λ6584 Å narrow components, respectively. The plus signs show the location of the three components identified in Figs. 4 and 5.

In the text

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 ∼ 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 cm−3 ≤ n ≤ 1000 cm−3, solar abundances, and magnetic field between 10−4 and 10 μ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.

In the text

	Fig. 8. Spatially resolved properties of the narrow component. From left to right: gas-phase metallicity (12 + log(O/H)), dust extinction (AV), and dust-corrected star formation rate surface density (Σ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 [O III] flux from Fig. 6.
In the text

	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.
In the text

	Fig. 10. Example of best-fit results for the fitting of the [C II] 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.
In the text

Fig. 11. [C II] 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σ threshold. In the left panels, we show the ALMA beam as a green ellipse and the narrow and broad [O III] flux from Fig. 6 as lime and red contours, respectively.

In the text

	Fig. 12. AV 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 [O III] channel maps and in optical continuum emission. The purple ellipse represents the beam size of the dust continuum emission.
In the text

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 [C II] 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 [C II] 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 [C II] 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 [C II] 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.

In the text

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.

In the text

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 AV = 0 or the AV from the outflow Balmer decrement in each ring, respectively. The half-blue, half-purple points have measured AV = 0. The purple and blue shaded regions in top-right panel show the 1σ region around the mean outflow rate computed with AV ≠ 0 and AV = 0, respectively.

In the text

	Fig. 15. Ionized mass outflow rate map from the dust-corrected broad Hα emission for each radial shell.
In the text

	Fig. 17. Mass loading factor of the ionized gas defined as Ṁion,out/SFR as a function of time for a range of assumed inclination angles.
In the text

	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.
In the text

	Fig. B.1. Outflow velocity map defined as vout = |Δvnarrow, broad|+2σout (see also Sect. 8). We show as black dashed contours the ΣSFR map contours (see also Fig. 8).
In the text

	Fig. C.1. Zoom in on the [O III] 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.
In the text