© The Authors 2024

1. Introduction

According to observations, the star formation rate density (SFRD) of the Universe peaked at around 3.5 − 6 Gyr after the Big Bang (1 < z < 2) and about 25% of the present-day stellar mass was already formed by then (Madau & Dickinson 2014). Consequently, galaxies observed at earlier epochs (< 3.5 Gyr) are expected to be actively star-forming and to be undergoing processes such as galaxy mergers and/or the accretion of large amounts of gas and stars, which could lead to significant mass growth and a global increase in the SFRD. The evolutionary paths that these early galaxies follow can be constrained through the study of the imprints that different drivers of galaxy evolution leave on galaxies. In particular, the ongoing star formation activity, properties of the interstellar medium (ISM) such as gas metallicity, dust content, and ionisation conditions, as well as the global kinematics of galaxies, can be probed through the rest-frame optical emission. However, until recently, detailed, spatially resolved studies of galaxies at high redshifts have only been possible with ground-based near-infrared integral field spectroscopic observations accessing the Hα line up to z ∼ 2.6 (Förster Schreiber et al. 2009, 2019; Stott et al. 2016; Wisnioski et al. 2019; Curti et al. 2020) and the [O III] λ5007 line up to z ∼ 3.5 (Troncoso et al. 2014; Turner et al. 2017). These previous works have revealed that galaxies at these epochs present a variety of metallicity gradients (Cresci et al. 2010; Troncoso et al. 2014; Wuyts et al. 2016; Curti et al. 2020), with a high incidence of ionised outflows in massive galaxies (Genzel et al. 2011; Förster Schreiber et al. 2019; Gupta et al. 2023; Llerena et al. 2023) but inconspicuous outflows in lower-mass ones (e.g. Concas et al. 2022), and an increase in the intrinsic velocity dispersion (higher pressure support) with redshift (Turner et al. 2017; Übler et al. 2019). Unfortunately, comparable studies of galaxies beyond z ∼ 3 are impossible because the bright Hα line is redshifted outside the near-infrared window observable from the ground.

Thankfully, after a successful launch and commissioning, the James Webb Space Telescope (JWST; Gardner et al. 2023) has opened a new era in the exploration of the early epochs of galaxy evolution through its access to near- and mid-infrared wavelengths. Since the start of scientific operations, JWST has provided numerous results on the properties of galaxies beyond z ∼ 3, which, compared to lower-redshift ones, are found to be more metal-poor (Schaerer et al. 2022; Curti et al. 2024; Nakajima et al. 2023a; Trump et al. 2023), display harder ionising spectra (Cameron et al. 2023; Sanders et al. 2023), and exhibit higher electron densities (Isobe et al. 2023; Reddy et al. 2023a).

Despite the wealth of information already obtained from spectroscopic observations of large samples of galaxies, most studies have been based on their integrated emission (e.g. Cameron et al. 2023; Carniani et al. 2024; Sanders et al. 2023), and a spatially resolved analysis of their properties is still missing. In this regard, the advent of integral field spectroscopy (IFS) instrumentation aboard JWST provides, for the first time, the possibility of spatially resolving the full rest-frame optical suite of emission lines, from [O II] to Hα + [N II] + [S II], in galaxies beyond z > 3, to study the physical properties of galaxies at early epochs. In particular, the NIRSpec instrument (Jakobsen et al. 2022) provides IFS capabilities (Böker et al. 2022) to study, in the near-infrared regime, the rest-frame optical emission from galaxies up to z ∼ 9. However, until now, only a few spatially resolved studies have been performed on galaxies within this redshift range (e.g. Perna et al. 2023; Übler et al. 2023; Wylezalek et al. 2022).

In this paper, we shed more light on the internal properties of high-redshift galaxies through the spatially resolved study, using NIRSpec IFS observations, of a galaxy group at z ∼ 3.7 that includes the massive, star-forming galaxy GS_4891. Previously observed in CANDELS (Grogin et al. 2011; Koekemoer et al. 2011) as part of the GOODS-South campaign, GS_4891 lies at a favourable redshift to cover a wide range of spectral features from [O II] λλ3726, 29 to [S II] λλ6716, 31 in a single NIRSpec spectroscopic band at a high spectral resolution.

This paper is organised as follows. In Sect. 2 we describe the NIRSpec observations, the reduction, and the analysis of the data. In Sect. 3 we present the results obtained on the ISM properties and kinematics in the GS_4891 system, together with a discussion of their implications. Finally, the summary of the main results and concluding remarks are presented in Sect. 4.

Throughout this paper, we adopt the initial mass function (IMF; 0.1 − 100 M_⊙) of Chabrier (2003) and a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7, and Ω_m = 0.3. Emission lines are referred to using their rest-frame air wavelengths, although for the analysis we use their vacuum wavelengths.

2. Observations and data reduction

2.1. JWST/NIRSpec IFS observations

The data presented in this paper are part of the NIRSpec IFS GTO program “Galaxy Assembly with NIRSpec IFS” (GA-NIFS; PIs: S. Arribas, R. Maiolino) and are included in proposal #1216 (PI: N. Luetzgendorf). Observations of GS_4891 (RA 3^h32^m18.29^s, Dec −27°51′59.17″, J2000) were executed on August 13 2022, with the IFS mode of NIRSpec covering a contiguous area of 3.1″ × 3.2″ with a native spaxel size of 0.1″ (Böker et al. 2022; Rigby et al. 2023). IFS observations of GS_4891 were obtained at high- (R2700) and low-resolution (R100) configurations, applying a medium (0.5″) cycling pattern of eight dithers. The R2700 observations were obtained with a total integration time of 15872.7 seconds (4.4 h) using the grating-filter pair G235H/F170LP, which provides a spectral resolution, R ∼ 1900 − 3500 between 1.7 μm and 3.15 μm (Jakobsen et al. 2022). A total of 3968.2 s (1.1 h) were devoted to the R100 observations (R ∼ 30 − 330 between 0.6 μm and 5.3 μm).

Since GS_4891 resides in the GOODS-South field, it has been previously observed with multi-wavelength observations, including WFC3/HST and Spitzer/IRAC. Based on these previous observations, GS_4891 has a spectroscopic redshift of 3.7027, a stellar mass of approximately 1.3 × 10¹⁰ M_⊙, and a star formation rate (SFR) of 49.5 M_⊙ yr⁻¹ (Guo et al. 2013; Pérez-González et al. 2005).

In Fig. 1 we show a false-colour image of the group of galaxies highlighting the field of view of the NIRSpec IFS observations presented here. The main galaxy, GS_4891, has a very close companion to the north, which we refer to as GS_4891_n, and another galaxy to the east at a projected distance of ∼15 kpc identified in CANDELS as GS_28356 (photo-z = 3.646). The green and red channels correspond to NIRCam F182M + F210M and F410M images, respectively, from the JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023), whereas the blue channel is a combination of ACS/HST archival data (F606W, F775W, F814W, F850LP; Giavalisco et al. 2004), re-processed following Illingworth et al. (2016) and Whitaker et al. (2019).

Fig. 1. False-colour image of the field around GS_4891 highlighting the position of the NIRSpec IFS field of view. Blue colours are a combination of ACS/HST images (Giavalisco et al. 2004), whereas green and red correspond to NIRCam images (Williams et al. 2023).

2.2. Data reduction

The raw data were reduced with the JWST calibration pipeline version 1.8.2 under CRDS context jwst_1068.pmap. Several modifications to the default reduction were introduced to improve the data quality; these are described in detail in Perna et al. (2023). The final cube was combined using the “drizzle” method, for which we used an official patch to correct a known bug¹. The main analysis in this paper is based on the combined R2700 cube with a pixel scale of 0.05″. During the analysis we noticed that the noise provided in the datacube (“ERR” extension) is underestimated compared to the actual noise in the data. For this reason, the noise vector in each spaxel was re-scaled based on the standard deviation in the continuum in regions free of emission lines, increasing the noise by a factor that varies with wavelength and ranges between 1.5 and 2. No background subtraction was performed on the R2700 data because in this paper we focus on the analysis of emission lines and any background contribution is taken into account in the spectral modelling. The same steps were applied to reduce the R100 data, with the addition of the background subtraction using spaxels away from the sources.

2.3. Analysis

2.3.1. R2700 data – emission line analysis

The main emission lines ([O II] λλ3726, 29, Hβ, [O III] λ5007, Hα, [N II] λλ6548, 6584, and [S II] λλ6716, 31) were modelled, in the spectrum of each spaxel and in the integrated spectrum of the regions defined in Sect. 3.1, with a set of individual Gaussian functions tied together to have the same kinematics. In this fit, the ratios [O II] λ3729/λ3726 and [S II] λ6716/λ6731 are constrained within [0.3839, 1.4558] and [0.4375, 1.4484], respectively, corresponding to the theoretical limits of low- (1 cm⁻³) and high- (10⁵ cm⁻³) density regimes estimated by Sanders et al. (2016). We also set a minimum value for the Hα/Hβ flux ratio of 2.86, which corresponds to an electron temperature, T_e = 10⁴ K, and an electron density, N_e = 100 cm⁻³, for a case B recombination (Osterbrock & Ferland 2006).

We explored the presence of additional kinematic components in the spectra by including a second set of Gaussian functions, with a larger velocity dispersion, in a separate fit to the brighter [O III] λ5007 and Hα lines, also tying their kinematics together. The two-component model was preferred when the Bayesian information criterion statistics (Schwarz 1978) was significantly better (> 60) than for the one-component model and the broad component presented a maximum flux that was at least two times larger than the scatter in the neighbouring continuum. In these cases, the rest of the emission lines ([O II] λλ3726, 29, Hβ and [N II] λλ6548, 6584) were modelled again, forcing them to have the same two kinematic components as in the fit to [O III] λ5007 and Hα, and only allowing the flux to vary in each component. The narrow kinematic component of the fit was associated with the galaxy and used to study its ISM properties and gas kinematics; the broad component was used to explore the presence and properties of a galactic outflow (see Sect. 3.3.2).

The signal-to-noise (S/N) in each line was computed as TotalFlux/RMS × $\sqrt{2\times FWHM}$, where TotalFlux is the summed flux within 2 × FWHM (full width at half maximum) of the line and the RMS is measured in windows of 20 pixels in the continuum around the emission lines. The adopted FWHM was taken from an initial fit performed on the brightest lines (i.e. [O III] λ5007 and Hα) on a spaxel-by-spaxel basis. In order to extract physical information from regions with low S/N, in some of the emission lines the spectra from different spaxels were combined following the Voronoi binning method from Cappellari & Copin (2003). It should be noted that the binning changes depending on the target emission line; in particular, analyses that require the [N II] λ6584 line (e.g. N2 metallicity indicator, excitation mechanisms) demand larger bins than those that use Hβ (e.g. nebular extinction, R2 and R3 metallicity indicators). Errors in the physical parameters derived in this paper were computed with Markov chain Monte Carlo techniques using the EMCEE software by Foreman-Mackey et al. (2013), assuming uniform priors. The reported 1-σ uncertainties on each quantity were calculated as half the difference between the 16th and 84th percentiles of the posterior probability density distribution. The spectral modelling took into account the wavelength-dependent spectral resolution of NIRSpec (Jakobsen et al. 2022).

2.3.2. R100 data – stellar mass estimation

Previous stellar mass estimates of GS_4891 were based on aperture photometry and also included the galaxy to the north, GS_4891_n (see Fig. 1), since they appeared to be two connected structures in HST imaging. In this work we study them separately as we find evidence that they might be different structures in an early stage of interaction. Thus, we needed to estimate their individual stellar masses.

The SED fitting to estimate stellar masses was done using the R100 NIRSpec IFS data (0.6 μm to 5.3 μm). For GS_4891, the modelling was performed at the spaxel level in the same region used to model the kinematics of the system (see central panel in the middle row of Fig. A.1), whereas for GS_4891_n we used the spaxels within the region demarcated in Fig. 2. The method is described in Pérez-González et al. (2023) and D’Eugenio et al. (2023). Briefly, the clear-prism NIRSpec/IFS observations were compared to stellar population models from the Bruzual & Charlot (2003) library, assuming a star formation history described by a delayed exponential characterised by a timescale, τ (taking values from 1 Myr to 1 Gyr in 0.1 dex steps), and age, t₀ (ranging from 1 Myr to the age of the Universe at the redshift of the galaxy). The stellar metallicity, Z, was left as a free parameter, enabling us to take all the discrete values provided by the Bruzual & Charlot (2003) library from 2% to 2.5 times solar. Nebular (continuum and line) emission was taken into account, asin Pérez-González et al. (2003, 2008). The attenuation of the stellar and nebular emission was modelled with a Calzetti et al. (2000) law, with A_V values ranging from 0 to 3 mag in 0.1 mag steps. The stellar mass was obtained by scaling the mass-normalised stellar model to the spectrum.

Fig. 2. Maps of integrated fluxes in the emission lines [O II] λλ3726, 29, [O III] λ5007, and Hα, and the continuum around 1.6 μm for the GS_4891 system. The three main galaxies identified in the system are highlighted with red rectangles in the top-right panel. The black cross indicates the centroid of continuum emission (see text for details). North is up and east is to the left.

Integrating the light in the regions around GS_4891_ext and GS_4891_n shown in the top-right panel of Fig. 2, we obtain stellar masses of 7.1 ± 0.3 × 10⁹ M_⊙ and 1.3 ± 0.1 × 10⁹ M_⊙, respectively. We note that the sum of these stellar masses is slightly lower than the previous estimations that included the two systems (1.3 × 10¹⁰ M_⊙); however, we ascribe this disagreement to small differences in the area encompassed by our individual apertures and the one used in the HST data.

3. Results and discussion

3.1. The complex neighbourhood around GS_4891

In Fig. 2 we show the maps of the distribution of the total [O II] λλ3726, 29, [O III] λ5007, and Hα fluxes, together with the continuum flux in the NIRSpec IFS R100 data in a wavelength window equivalent to the HST F160W/WFC3 filter, tracing the rest-frame UV emission. The black cross in this figure marks the centroid of continuum emission computed with a 2D Gaussian fit. The three galaxies display spatially resolved emission in the two lines over several kiloparsecs. The red rectangles shown in the top-right panel of Fig. 2 were chosen to encompass most of the emission from each system and to study their individual properties. The integrated NIRSpec G235H/F170LP spectra from these apertures are shown in Fig. 4, which demonstrates that the emission lines are robustly detected with high S/N in a wide spectral range. The similar observed wavelength of the emission lines in their spectra confirms that they reside at the same redshift and are part of a galaxy group. The three galaxies exhibit clear detections of strong emission lines such as the [O II] λλ3726, 29 doublet, [O III] λλ4959, 5007, and the hydrogen Balmer recombination lines, Hβ and Hα. Additionally, GS_4891 also displays other emission lines such as higher-order Balmer transitions lines (Hγ, Hδ, Hϵ) as well as [Ne III] λ3869, He Iλ5876, [O I] λ6300, [N II] λλ6548, 6584, and [S II] λλ6716, 31.

Apart from the three main structures, the map of the [O III] λ5007 flux displays evidence of substructure to the southeast of GS_4891, which is less clearly seen in Hα. The presence of this and further substructure is investigated in more detail through the flux maps shown in Fig. 3 that correspond to windows of two spectral pixels (∼50 km s⁻¹ wide) centred at 25 km s⁻¹ (left panel) and −170 km s⁻¹ (right panel) with respect to the position of [O III] λ5007 at the systemic redshift (z = 3.7027). The plot on the left shows a clump of ionised gas (C1) to the southeast of GS_4891 at the same velocity, whereas an additional clump (C2) is seen in the southwest moving at approximately −150 km s⁻¹ with respect to the main galaxy. These two clumps can also be identified in the composite HST+NIRCam image shown in Fig. 1 through their blue colours, indicative of strong rest-frame UV radiation likely coming from young stars. This composite image also shows a compact blue clump in the north of the galaxy, closer to its centre than the other two clumps, which does not stand out in the NIRSpec data. The observed internal substructure suggests either in situ star formation or ongoing minor mergers.

Fig. 3. Flux maps in windows of two spectral pixels (∼50 km s−1 wide) centred at 25 km s−1 and −170 km s−1 with respect to the position of [O III] λ5007 at the systemic redshift (z = 3.7027) where the C1 and C2 separate clumps (magenta circles) are identified. The origin of the coordinates corresponds to RA 3h32m18.29s, Dec −27°51′59.17″, J2000. North is up and east is to the left.

3.2. Properties of the ISM

3.2.1. Electron densities

The detection of the [O II] λλ3726, 29 and [S II] λλ6716, 31 doublets in the integrated spectra of GS_4891 (see Fig. 4) offers the possibility of computing the electron density of the ionised gas (Osterbrock & Ferland 2006). Although each individual set of emission lines can be used for this computation, the detection of both of them can be leveraged to provide a better constraint on the electron density. Assuming a typical electron temperature for HII regions, T_e = 10⁴ K, and following the expressions derived by Sanders et al. (2016), the [O II] λλ3726, 29 and [S II] λλ6716, 31 in the integrated spectra of GS_4891 were simultaneously modelled, forcing them to have the same kinematics and to agree on the value of the electron density. In this fit, we also constrained the line ratios within the same theoretical limits from Sanders et al. (2016) stated in Sect. 2.3. This modelling, shown in Fig. 5, yields an electron density, n_e = 776 ± 307 cm⁻³, whereas the independent modelling of [O II] λλ3726, 29 and [S II] λλ6716, 31 gives electron densities of 533 ± 357 cm⁻³ and 357 ± 309 cm⁻³, respectively. The three estimates are consistent within the (large) uncertainties and provide a range of electron densities compatible with our data.

Fig. 4. Integrated spectra of GS_4891 and the two other systems, GS_4891_n and GS_28356, identified in the NIRSpec/IFS observations. The integrated spectra were obtained by adding all the spaxels within the apertures highlighted in Fig. 2.

Fig. 5. Simultaneous spectral modelling (orange lines) of the observed [O II] λλ3726, 29 and [S II] λλ6716, 31 (black lines) to constrain the electron density of the ionised gas in the integrated spectrum of GS_4891.

Our high n_e measurements in GS_4891 are in agreement with an overall increase in electron density with redshift observed in galaxies at z < 3 (Steidel et al. 2014; Sanders et al. 2016; Kashino et al. 2017; Davies et al. 2021). However, our values are also higher than the median electron densities, n_e ∼ 300 cm⁻³, measured using [O II] λλ3726, 29 in a sample of ten galaxies with masses around 10⁹ M_⊙ at 4 ≲ z ≲ 6 by Isobe et al. (2023) and much higher than in the sample of 24 massive (mean log(M_⋆/M_⊙) = 9.86) galaxies at 2.7 < z < 6.3 studied in Reddy et al. (2023a), where they find a mean n_e < 100 cm⁻³ using [S II] λλ6716, 31. Nevertheless, our estimates of the electron density are still within the large scatter [40 − 600 cm⁻³] observed in these previous works.

In their work on a total of 48 galaxies at 2.7 < z < 6.3, Reddy et al. (2023a) find that the electron density increases for higher values of the O32 ([O III] λ5007/([O II] λλ3726, 29)), with the highest densities, n_e > 500 cm⁻³, being associated with O32 > 2. Such results are interpreted as a consequence of the positive correlation between the SFR surface density (Σ_SFR, whose increment leads to higher electron densities) and the ionisation parameter (log U, traced by O32), both found to increase with redshift. At high masses they find that galaxies have low electron densities, n_e < 100 cm⁻³, and O32 < 6, due to the lower Σ_SFR of more massive systems. In the case of GS_4891, the extinction-corrected estimate of O32, 1.36 is lower than the O32 values [2 − 8] measured for electron densities > 500 cm⁻³ in Reddy et al. (2023a), whereas Σ_SFR = 6.6 ± 0.1 M_⊙ yr⁻¹ kpc⁻², estimated within the aperture highlighted in Fig. 2, is in agreement with the range they obtain [2 − 20 M_⊙ yr⁻¹ kpc⁻²] in the high-density regime.

Overall, a comparison with previous works indicates that GS_4891 presents an electron density that is higher than expected for a massive galaxy at z = 3.7, although the uncertainties are large. A possible explanation could be the presence of an ionised outflow in the central regions of GS_4891 (see Sect. 3.3.2), since galactic outflows are generally found to have higher-density gas (Arribas et al. 2014; Villar Martín et al. 2014; Mingozzi et al. 2019; Rodríguez Del Pino et al. 2019). The electron density measurements in GS_4891 provide an estimate of the physical conditions of the ISM, reinforcing the assumption of a lower limit of Hα/Hβ = 2.86 (a value that changes to 2.85 only at N_e = 10⁴ cm⁻³, Groves et al. 2012) that is used in the spectral modelling of our data (Sect. 2.3) and the computation of the nebular extinction (Sect. 3.2.2).

3.2.2. Extinction

The availability of several hydrogen Balmer emission lines in the integrated spectrum of GS_4891 enables us to evaluate the consistency of their flux ratios with standard attenuation curves. Using the Hδ, Hγ, Hβ and Hα fluxes we derive the nebular extinction, A_V, assuming the attenuation laws from Cardelli et al. (1989, R_V = 3.1) and Calzetti et al. (2000, R_V = 4.05). The A_V values derived from the flux ratios between Hα and the higher-order Balmer lines for each attenuation law are shown in Fig. 6. All estimations of A_V present good agreement within the errors and are consistent with A_V ∼ 0.9 − 1 mag, indicating that in GS_4891 both attenuation laws are compatible. Such consistency with local attenuation laws was also observed in stacks of star-forming galaxies at 1.4 < z < 2.6 by Reddy et al. (2020) and, together with our results, suggests that at high redshifts (at least up to z ∼ 3.7) the effects of dust along the line of sight are similar to those in the local Universe. Thus, throughout this work we computed the nebular extinction from the Balmer decrement (Hα/Hβ) following the method described in Domínguez et al. (2013), which assumes Calzetti’s attenuation law and a minimum Hα/Hβ = 2.86.

Fig. 6. AV values for the integrated spectrum of GS_4891 obtained using different pairs of hydrogen Balmer lines. The purple and orange points were obtained using the attenuation curves from Cardelli et al. (1989) and Calzetti et al. (2000), respectively.

A map of the nebular extinction for the galaxies in our system is shown in panel a of Fig. 7, after binning the data to have S/N ≥ 3 in Hβ. The integrated A_V of each system within the apertures defined in Fig. 2 are included in Table 1. Throughout the whole system, the nebular extinction varies by more than two magnitudes, including regions with almost no extinction. GS_4891 presents A_V of roughly one magnitude at the centre, consistent with its integrated value. In the case of GS_4891_n, the integrated nebular extinction is quite high (A_V ∼ 2 mag), but displays a large internal scatter. Interestingly, the region between GS_4891 and GS_4891_n presents high extinction, A_V ≥ 1.5, despite it not hosting significant star formation activity (black contours corresponding to the Hα emission). On the contrary, clumps C1 and C2 present very low extinction, A_V < 0.5. The spatially resolved and integrated A_V estimates were used to correct for nebular attenuation the measured line fluxes in the combined spectra. However, for the map of the SFR shown in Fig. 7, we kept a spaxel-by-spaxel approach, applying the same extinction correction to all the spaxels within a given bin.

Fig. 7. Maps of the main ISM parameters derived for the GS_4891 galaxy group from emission line analysis. (a) Nebular extinction, AV, derived from the Balmer decrement (Hα/Hβ). (b) SFRs (M⊙ yr−1) estimated from the Hα luminosity. (c) Gas oxygen abundance obtained with the R2, R3, R23, and O32 metallicity calibrations from Curti et al. (2017). (d) Ionisation parameter, U, computed with the gas metallicity and the [O II] λλ3726, 29 lines following Eq. (5) in Díaz et al. (2000). The data were spatially binned to have a S/N ≥ 3 in the emission lines involved. The dashed black contours correspond to the Hα emission, whereas the black cross indicates the centroid of continuum emission.

3.2.3. Excitation mechanisms

The NIRSpec IFS observations of GS_4891 allow us, for the first time, to spatially resolve the sources of ionisation in a galaxy at z ∼ 3.7 using rest-frame optical emission lines. This is an intermediate redshift between the well-studied range z < 2.5 via ground-based observations (e.g. Strom et al. 2017) and the recently analysed cases at z > 5 with JWST (Cameron et al. 2023; Sanders et al. 2023). Here, we explore the sources of ionisation throughout the galaxies in the GS_4891 system by using the standard diagnostic diagram (Baldwin et al. 1981) that includes the rest-frame optical emission line ratios [O III] λ5007/Hβ and [N II] λ6584/Hα. We also consider an alternative diagram that uses [S II] λλ6716, 31 instead of [N II] λ6584 (Veilleux & Osterbrock 1987) for the integrated spectrum of GS_4891 ([O I] λ6300 could also be used instead but is barely detected with S/N < 2). The [N II]- and [S II]-based Baldwin, Phillips, and Terlevich (BPT) diagrams are shown in Fig. 8 together with a map of the distribution of the different sources of ionisation in our system following the [N II]-based classification. We note that clump C2 is not shown in the BPT diagram because the [N II] λ6584 line has S/N < 3. The diagrams include reference lines from studies of local galaxies that demarcate the regimes that could be explained considering only ionisation from star formation (below the black, dashed line from Kauffmann et al. 2003) and the one where theoretical models require the presence of active galactic nucleus (AGN) activity to explain the ionisation state of the gas (black, dotted lines from Kewley et al. 2001). The “intermediate” regime (between both lines) corresponds to ionisation that cannot be explained either by star formation or AGN activity alone (Cid Fernandes et al. 2010). The solid black line is aimed at separating AGNs from LIERs (low-ionisation emission-line regions; Monreal-Ibero et al. 2006; Belfiore et al. 2016), as was suggested by Cid Fernandes et al. (2010). In addition, we also show as a cyan line the best-fit curves tracing the location of the star-forming population at z ∼ 2.3 from Strom et al. (2017).

Fig. 8. Left: map of excitation mechanisms throughout GS_4891 and GS_4891_n according to the location in the [N II]-based BPT diagnostic diagram. The dashed black contours correspond to the Hα emission, whereas the black cross indicates the centroid of continuum emission. Centre and right: [N II]- and [S II]-based BPT diagrams containing the spatially binned regions (red dots), the integrated regions of the system where the relevant lines are detected with S/N ≥ 3 (red square and diamond, and pink circle), and the broad kinematic component associated with the ionised outflow identified in GS_4891 (blue circle). Clump C2 is not shown in any of the diagrams due to the low S/N in Hβ and [N II]. As a reference, we have added the stack of galaxies at z ∼ 6 from Cameron et al. (2023, purple plus), the stacks in different high-redshift bins from Sanders et al. (2023, yellow symbols), the obscured AGN at z ∼ 3.3, “Jil2”, from Perna et al. (2023, lightgreen star), and the broad-line AGN, GS_3073, at z ∼ 5.5 studied in Übler et al. (2023, darkgreen star). The shaded grey contours demarcate the regions in the diagram encompassed by 70, 80, and 90% of the local SDSS sample (Abazajian et al. 2009). Dashed (Kauffmann et al. 2003) and dotted black lines (Kewley et al. 2001) demarcate the separation between ionisation coming from star formation and AGN activity based on studies of local galaxies (see Sect. 3.2.3 for details). The cyan lines correspond to the best-fit curves to the star-forming population at z ∼ 2.3 from Strom et al. (2017).

According to these demarcation lines, there is no evidence for the presence of an AGN in either GS_4891 or GS_4891_n, since all the points (with the exception of one region away from the nucleus of GS_4891 that is probably spurious) lie below the Kewley et al. (2001) line. However, we have to take into account that standard BPT diagrams are designed for local galaxies and have been proven to be less reliable in constraining excitation mechanisms in higher-redshift galaxies due to their lower metallicities. In fact, photoionisation models (e.g. Nakajima et al. 2023b) predict that low-metallicity AGNs have lower [O III]/Hβ and [N II]/Hα ratios, locating them in regions of the BPT diagram that correspond to “intermediate” and star-forming regimes. Such an effect can be observed in high-redshift AGNs selected through other tracers such as the detection of a broad-line region in the Balmer emission lines (Harikane et al. 2023; Kocevski et al. 2023; Maiolino et al. 2023; Oesch et al. 2023) or alternative diagnostic diagrams involving other emission lines that require harder ionising radiation such as HeIIλ4686 (Übler et al. 2023). As an example, we show in Fig. 8 (dark green star) the low-metallicity AGN (GS_3073; Z_gas/Z_⊙ ∼ 0.21) at z ∼ 5.5 studied in Übler et al. (2023) that, despite being classified as star-forming following the standard BPT diagram, was confirmed to be an AGN based on the presence of a broad-line region and the diagnostic diagram involving the HeIIλ4686 line. We also show the obscured AGN at z ∼ 3.3, “Jil2” (light green star), that lies very close to the demarcation line from Kewley et al. (2001) but for which Perna et al. (2023) inferred a HeIIλ4686 consistent with AGN ionisation. Although these previous works demonstrate that the BPT diagram must be used with caution in high-redshift studies, our NIRSpec IFS spectra do not present evidence of a broad-line region in either Hβ or Hα, and the higher ionisation emission line HeIIλ4686 is not detected (Fig. 4). Moreover, the gas metallicities we measure are Z_gas/Z_⊙ ∼ 0.35 − 0.6 (Sect. 3.2.5), a regime where excitation diagnostics using the BPT diagram should still be reliable. Finally, the [S II]-based BPT classification (right panel in Fig. 8) for the integrated spectrum of GS_4891 is also consistent with ionisation from star formation.

Figure 8 shows that the spatially resolved and integrated regions in our system are clearly offset from the parameter space covered by local star-forming galaxies (SDSS; Abazajian et al. 2009). If our system were at z ∼ 0, such deviations could be interpreted as signs of harder ionisation mechanisms such as those originating in shocks induced by galactic outflows (e.g. Ho et al. 2014) or by mergers (e.g. Perna et al. 2020). In fact, as is described in Sect. 3.3.2, we find evidence of an ionised outflow in the central regions of GS_4891 oriented towards the northwest, a direction in which there is an increase in the global velocity dispersion (right panel of Fig. 11) and in the [N II] λ6584/Hα ratio (panel d in Fig. 9), a correlation expected in regions ionised by shocks (e.g. Ho et al. 2014). In the case of GS_4891_n, the [N II] λ6584/Hα is even higher and it presents velocity dispersions (σ ∼ 60 km s⁻¹) that are lower but still consistent with the turbulence expected from shocks (Johnston et al. 2023). Moreover, the substructure observed in GS_4891 and the close distance to GS_4891_n are compatible with ongoing or past mergers that could induce shocks. However, despite shocks possibly contributing to the observed line ratios, it is already well established that star-forming galaxies at high redshifts tend to display higher [O III] λ5007/Hβ at fixed [N II]/Hα ratios than local galaxies due to changes in their physical conditions (Brinchmann et al. 2008; Kewley et al. 2013; Shapley et al. 2015; Sanders et al. 2016; Steidel et al. 2016; Strom et al. 2017). This tendency is clearly illustrated by the location in Fig. 8 of the star-forming population at z ∼ 2.3 from Strom et al. (2017) and the higher-redshift stacks from Cameron et al. (2023) and Sanders et al. (2023). The GS_4891 data points cluster together slightly above the star-forming population at z ∼ 2.3, although they are still consistent within the intrinsic scatter (Fig. 5 in Strom et al. 2017). Interestingly, clump C1 is offset from the rest of the GS_4891 system due to its lower metallicity and higher ionisation, being the only region close to the z ∼ 3.3 stack (yellow triangle). However, we note that this data point from Sanders et al. (2023) corresponds to a wide redshift range, 2.7 < z < 4.0, with a larger scatter in the BPT diagram (see Fig. 3 in their paper) that is consistent with the line ratios observed in the GS_4891 system.

Fig. 9. Maps of emission line ratios that trace gas metallicity throughout the GS_4891 system. The data have been spatially binned to have a S/N ≥ 3 in the emission lines involved. The black cross indicates the centroid of continuum emission, whereas dashed contours in the right panel correspond to the Hα emission. Black boxes and circles highlight the individual systems identified in Sect. 3.1.

In summary, the GS_4891 system occupies a region in the BPT diagram in between z ∼ 2 and z > 4 galaxies, consistent with an evolution in the ionisation state with redshift. Moreover, the possible merging activity of GS_4891 and the presence of an ionised outflow indicate that shocks could also contribute to the observed line ratios.

3.2.4. Star formation rates

The results from the previous section demonstrate that the ionisation in the GS_4891 is dominated by star formation. Therefore, the extinction-corrected Hα fluxes can be used to estimate the ongoing SFRs following the calibrations listed in Kennicutt (1998), consistent with a Salpeter IMF (Salpeter 1955). A map of the SFR is shown in panel b of Fig. 7, whereas the integrated values of each individual region are included in Table 1. As expected, the distribution of the SFR is very similar to that of the Hα emission, peaking at the central regions of GS_4891 and GS_4891_n, with a gap of low SFR between both systems. The total integrated SFR (Hα) of these two systems is approximately 52 M_⊙ yr⁻¹, which is very similar to the one previously estimated from SED modelling of HST+IRAC/Spitzer photometry, 49.5 M_⊙ yr⁻¹ (Pérez-González et al. 2005), also using a Salpeter IMF. However, as can be noted, using the individual apertures highlighted in Fig. 2 we are not including all the Hα flux of the system. If we sum all the Hα flux in the GS_4891 and GS_4891_n structures, we obtain a total SFR of ∼125 M_⊙ yr⁻¹, doubling previous SED estimations. This difference in the SFRs estimated using Hα and SED fitting might be due to the fact that Hα emission traces star formation on shorter time scales than the UV and IR emission traced in the SED fitting, implying that GS_4891 is experiencing a burst of star formation at the time of observations. Taking into account their stellar masses and SFRs, GS_4891 and GS_4891_n would lie approximately 0.34 dex and 0.43 dex, respectively, above the main-sequence at z = 3.7 (Speagle et al. 2014), indicating that these are starburst systems. Regarding GS_28356, the Hα-based SFR, 2.2 ± 0.5 M_⊙ yr⁻¹, is the same as the one previously estimated from SED fitting, 2.2 M_⊙ yr⁻¹.

3.2.5. Gas metallicities

The study of the spatially resolved metal content of the ISM is critical to identifying regions in a galaxy that have been more metal-enriched due to subsequent episodes of star formation and those where diluted, pristine gas from the circumgalactic medium is being accreted. Leveraging the detection of several emission lines throughout our spatially resolved spectra, we can apply well-established calibrations (e.g. Curti et al. 2017) to constrain the oxygen abundance in the GS_4891 group of galaxies. After correcting for dust attenuation, we use flux ratios between specific emission lines that correspond to standard metallicity indicators: R2 ([O II] λλ3726, 29/Hβ), R3 ([O III] λ5007/Hβ), O32 ([O III] λ5007/[O II] λλ3726, 29), R23 (([O II] λλ3726, 29 + [O III] λλ4959, 5007)/Hβ), and N2 ([N II] λ6584/Hα). The spatial distributions of four of these indicators, obtained after spatially binning the data to have S/N ≥ 3 in the emission lines involved (Hβ for the first three panels and [N II] λ6584 for the last one), are shown in Fig. 9. The metallicity calibrations from Curti et al. (2017) for the R2, R3, R23, and O32 parameters were fitted together to generate the oxygen abundance map presented in panel c of Fig. 7 (see Curti et al. 2024 for details). The N2 indicator is not included in this estimation as it requires coarser spatial binning. The integrated oxygen abundances in each system were also estimated using the R2, R3, R23, and O32 indicators, with the addition of N2 for GS_4891 and C1 where [N II] λ6584 is detected. Their values are presented in Table 1. We note that, since GS_4891 is at z = 3.7, the calibrations provided by Sanders et al. (2024) based on galaxies at 2 < z < 9 should, in principle, be more appropriate to computing the metallicities than the ones from Curti et al. (2017) that are based on local galaxies. However, the calibrations from Sanders et al. (2024) do not provide good constraints at the regime of relatively high oxygen abundances (12 + log(O/H) > 8.3) estimated in the GS_4891 system.

The three galaxies in the group display a wide range of oxygen abundances, from 12 + log(O/H) ∼ 8.1 to ∼8.5. In particular, GS_4891 presents a clear gradient in metallicity from the south (including clump C1) to the northwest. The value at the centre, 12 + log(O/H) ∼ 8.3, is 0.2 dex higher than in the south and ∼0.15 dex lower than in the northwest. Regarding the drop in metallicity at the south, the low velocity dispersion (∼50 km s⁻¹) estimated from the kinematic analysis performed in Sect. 3.3.1 (see also Fig. 11) suggests that this is a region where lower-metallicity gas is likely being accreted through a minor merger. The positive gradient observed from the centre of GS_4891 to the northwest resembles those previously reported in galaxies at z ≥ 3 by Curti et al. (2020) and Cresci et al. (2010), and recently at z ∼ 7 by Arribas et al. (2023), which have generally been ascribed to the dilution of metallicity in the central regions as a consequence of new episodes of star formation fuelled by cold flows of pristine gas (Dekel et al. 2009). Overall, the gradients observed in the oxygen abundance map are also visible in the individual line ratios presented in Fig. 9. In particular, the N2 indicator, which was not included in the estimation of the oxygen abundance and which correlates with metallicity (Pettini & Pagel 2004; Curti et al. 2017), also displays a clear increase from the southeast to the northwest of GS_4891.

In addition to the observed metallicity gradient, probably the most striking result is the difference of approximately 0.2 dex between the metallicities at the centre of GS_4891 and GS_4891_n, given that the former is four times more massive and that it consequently should be more efficient at retaining metals (e.g., Tremonti et al. 2004). In fact, while GS_4891 has an integrated metallicity (12 + log(O/H) = 8.4) in agreement with the mass–metallicity relation at z ∼ 3.3 (Sanders et al. 2021) given its mass (5.5 ± 0.6 × 10⁹ M_⊙), the metallicity of GS_4891_n (12 + log(O/H) = 8.5) is ≥0.3 dex higher than the mean value in galaxies of similar mass (1.3 ± 0.1 × 10⁹ M_⊙). A possible explanation for this result is that the most massive galaxy, GS_4891, is the main beneficiary of inflows of lower-metallicity gas following the baryon cycle that maintains the mass–metallicity relation. Such a scenario is supported by the higher SFR density (the same as the SFR in Fig. 2) in GS_4891 that is expected in the case of the accretion of pristine gas (Ceverino et al. 2016). Furthermore, the higher metallicity in GS_4891_n could be a consequence of gas stripping and “starvation” (or “strangulation” Larson et al. 1980; Peng et al. 2015) by the surrounding environment. Considering GS_4891_n as a satellite galaxy and GS_4891 as the central one, the difference in their metallicities resembles the findings at lower redshifts (z < 0.1) where satellite galaxies are invariably more metal-rich than centrals of the same mass (e.g. Pasquali et al. 2012; Peng & Maiolino 2014; Lian et al. 2019; Schaefer et al. 2019). Two scenarios have been proposed to explain this trend: either satellites accrete relatively more metal-rich gas (pre-processed by their more massive centrals, e.g. Peng & Maiolino 2014; Gupta et al. 2018), or, alternatively, satellites are more metal-rich due to a combination of starvation (which truncates gas accretion), ram-pressure stripping (which preferentially removes low-metallicity gas situated at large radii), and a lower impact of metal-loaded outflows (e.g. Pasquali et al. 2012; Rodríguez Del Pino et al. 2023). In the case of GS_4891_n, we see enrichment of a galaxy relative to its own central, which could favour the second scenario. The metallicity difference between GS_4891 and its satellite could also explain the positive radial metallicity gradient observed from the centre to the northwest of GS_4891: if the satellite already passed the pericentre, some of the high-metallicity gas to the northwest of GS_4891 could be stripped of material originally belonging to GS_4891_n. However, it would be difficult to explain the lower velocity dispersions found in the region between the two galaxies. We cannot draw general conclusions from a single system that may even represent the exception instead of the rule but, if our interpretation is correct, this system may be one of the earliest cases of environmental effects known. In any case, our findings demonstrate that JWST may help long-standing problems in the local Universe to be solved.

Finally, the NIRSpec IFS observations of GS_4891 provide evidence of gas metallicities varying significantly within single systems at high redshifts. This fact is very relevant for studies based on multi-object spectroscopic observations that only probe regions of galaxies falling within the slits and not the whole systems, indicating that variations in metallicity within single systems might also contribute to the observed scatter between galaxies.

3.2.6. Ionisation parameter

Our dataset also allows the exploration of the ionisation parameter, U, which is the ratio of the ionising photon density to the particle density and which provides a measure of the degree of ionisation of the gas. This parameter can be estimated directly through O32, following Eq. (7) in Díaz et al. (2000); however, since this indicator is used in the computation of the gas metallicity, the use of the previous equation leads by default to an anti-correlation between log U and the gas metallicity. Instead, we employed Eq. (5) from Díaz et al. (2000), which combines the gas metallicity, which has been well constrained above using several indicators, with the O2 ratio. A map of log U is shown in panel d of Fig. 7, whereas the integrated values for the three galaxies in the group and the two stellar clumps are included in Table 1.

GS_4891 presents the highest log U values >  − 2.6 of the system in the south, in the region around the C1 clump, the centre, and in the northwest. In GS_4891_n the values are generally lower, log U < −2.6, whereas GS_28356 displays an intermediate value of log U = −2.8. Although log U is found to anti-correlate with gas metallicity in local galaxies (e.g. Pérez-Montero 2014; Sánchez et al. 2015), such behaviour is less clear at higher redshifts (Strom et al. 2018; Reddy et al. 2023b). In GS_4891 we observe anti-correlations between log U and metallicity in clump C1, where log U peaks and the metallicity is the lowest of the system. However, there are regions of high metallicity, such as the northwest of the galaxy, that also have a relatively high ionisation parameter. These results suggest that the physical conditions to hold an anti-correlation between log U and metallicity might be present in individual regions of high-redshift galaxies, but not in the rest of the system.

In Fig. 10 we compare R23 and O32 in the different regions of our galaxies with those found in local galaxies (Abazajian et al. 2009) at z ∼ 2 (Sanders et al. 2016) and up to z > 9.5 (Cameron et al. 2023; Sanders et al. 2023). First of all, it is worth noticing the wide area of the diagram covered by the different regions in our galaxies, encompassing more than 0.2 and 0.5 dex in R23 and O32, respectively. Despite the large errors, GS_4891_n (red diamond) stands out from the rest of the sample due to its low O32 but high R23. In comparison with galaxies at 2 < z < 4, regions in GS_4891 are spread similarly in O32 but are more concentrated towards higher values of R23. Such high values of R23 are also shared by the high-redshift population, but they present ≥0.2 dex higher O32. Considering that the ionisation parameter correlates with O32 and the gas metallicity anti-correlates with R23, the results shown in Fig. 10 indicate that the GS_4891 galaxy group at z ∼ 3.7 is at an intermediate stage sharing roughly similar metallicities with the higher-redshift sample but lower ionisation, whereas galaxies at z ≤ 4 continue evolving towards lower ionisation but their gas metallicities increase.

Fig. 10. R23−O32 diagram for the GS_4891 system in comparison with other works. Spatially binned regions are shown as red dots, whereas integrated values are shown as red polygons and empty and filled magenta circles. As a reference, we have added the stacks of galaxies at z ∼ 2.5, 3.5 from Kriek et al. (2015, light and dark blue points, respectively), the stacks in different redshift bins from Sanders et al. (2023, yellow symbols), and the individual galaxies at z ∼ 6, 8 from Cameron et al. (2023, purple and brown pluses, respectively). Shaded contours demarcate the regions in the diagram encompassed by 70, 80, and 90% of the local SDSS sample (Abazajian et al. 2009).

3.3. Ionised gas kinematics

3.3.1. Global kinematics

We extracted the kinematic information of the ionised gas from the spectral modelling of the [O III] λ5007 and Hα lines performed in Sect. 2.3 tying together their velocities. These are the strongest lines in the spectra and allow us to constrain the extended kinematics throughout our system. Figure 11 contains the velocity maps at the 2th, 50th, and 98th percentiles (V₀₂, V₅₀, V₉₈) and the velocity dispersion traced by the non-parametric estimate W80 (the line width that contains 80% of the flux, obtained as the difference between the velocities at the 90th and 10th percentiles). First of all, the velocity maps confirm that the three galaxies lie at the same redshift, indicating that they belong to a common group. In the mid panel of Fig. 11, the V₅₀ map shows velocity gradients in GS_4891 and GS_4891_n suggestive of rotating gas in both systems. As was already identified in Fig. 3, clump C2 is moving at blue-shifted velocities of more than 100 km s⁻¹ with respect to GS_4891, whereas clump C1 displays similar kinematics to the rest of the galaxy. The maximum V₅₀ differences are 175 km s⁻¹ in GS_4891 and 98 km s⁻¹ in GS_4891_n. The W80 map shows a clear increase in the northwest region of GS_4891, reaching values higher than 300 km s⁻¹, indicative of higher turbulence, compared to the lower values observed in other regions of GS_4891 and GS_4891_n. This increase in W80 suggests that the ionised outflow detected in the nuclear parts of the galaxy potentially extends towards the northwest, a scenario that is explored in detail in the next section. GS_28356 displays a small variation in velocities (< 50 km s⁻¹) and low W80 (< 150 km s⁻¹).

Fig. 11. Velocity maps of the ionised gas in the GS_4891 system obtained from the fitting of the [O III] λ5007 and Hα lines. The three left panels show the velocities at the 2nd, 50th, and 98th percentiles of the velocity distributions. The right panel shows the W80 map. The black cross indicates the centroid of continuum emission, whereas the dashed black contours correspond to the Hα emission.

Motivated by the velocity gradients observed in GS_4891, we attempted to model the gas kinematics with 3D-BAROLO (Di Teodoro & Fraternali 2015) following the method described in Perna et al. (2022). From this analysis (presented in Appendix A) we obtained a rotational velocity, v_rot, of ∼120 km s⁻¹ and a velocity dispersion of the gas, σ_gas, of ∼90 km s⁻¹ at a radius of 2.2 kpc. Using Eq. (3) in Perna et al. (2022), we first estimated the circular velocity, v_c ∼ 180 km s⁻¹, assuming an asymmetric drift term η = 3; then, using Eq. (4) and assuming that the source of the gravitational potential is spherically distributed, we estimated that the dynamical mass within a radius of 2.2 kpc is M_dyn ∼ 1.7 × 10¹⁰ M_⊙. We estimated a stellar mass within the area used for the kinematic modelling of 5.5 ± 0.3 × 10⁹ M_⊙, ∼2 times smaller than the dynamical mass, indicating that this galaxy is dominated by dark matter.

3.3.2. Detection of an ionised outflow

Studies of galaxies at z ≥ 1 have revealed that galactic outflows are a relatively common phenomenon in massive (≥10¹⁰ M_⊙) galaxies (Perna et al. 2018; Förster Schreiber et al. 2019; Swinbank et al. 2019; Concas et al. 2022), although at these masses outflows are more generally found in AGN systems (e.g. Concas et al. 2022). The presence of ionised outflows in the GS_4891 system is explored through the modelling of the spectra using two kinematic components, as was described in Sect. 2.3. This analysis reveals the presence of a second, broader kinematic component in several individual spaxels around the nuclear parts of GS_4891. The maps of the mean velocity (v_broad) and velocity dispersion (σ_broad) of the broad component are shown in Fig. 12, superimposed on the Hα flux distribution in GS_4891. Motivated by the observed increase in W80 towards the northwest of GS_4891, we performed an additional kinematic modelling of the integrated spectrum extracted from the region with W80 > 200 km s⁻¹, highlighted with filled, dashed contours in the maps. In this integrated spectrum we also find evidence for a secondary, broad component in the [O III] λ5007 and Hα lines according to the criteria defined in Sect. 2.3.

Fig. 12. Velocity and velocity dispersion maps of the broad kinematic component superimposed on the Hα flux distribution in GS_4891 (grey background). The black cross indicates the centroid of continuum emission.

The region where the outflow was detected in individual spaxels spans a total of ∼4.4 kpc², from which we estimated a circularised radius of the outflow of 1.2 kpc. However, given that a broad component was detected in the integrated region to the northwest, we considered this radius to be a lower limit and adopted an alternative radius of the outflow, r_out = 1.5 kpc, which corresponds to the distance from the centre of the galaxy to half the width (approximately two spaxels) of the integrated region across the northwest direction. The velocity of the broad gas component, v_broad, shows a gradient from the southeast (∼50 km s⁻¹) to the northwest (∼ − 90 km s⁻¹) of the galaxy centre, whereas the velocity dispersion, σ_broad, is fairly constant throughout the extension of the outflow and only increases in the integrated region to the northwest. The extension and kinematic properties indicate that the broad emission is likely associated with a blue-shifted, galactic outflow originating at the centre of GS_4891. This outflow would have a maximum projected velocity, v_out, of ∼400 km s⁻¹, estimated as the highest value of ⟨v_broad⟩+2σ_broad (e.g. Fiore et al. 2017; Carniani et al. 2024; Übler et al. 2023), considering only the regions where a broad component is detected.

We explored the properties of the outflowing gas by integrating all the individual spaxels where a broad component is detected. The region to the northwest was not included because the lower S/N in emission lines different to [O III] λ5007 and Hα reduces the quality of the fit. These line fluxes were used to estimate the nebular attenuation (A_V, out) and the gas metallicity (12 + log(O/H)_out) in the outflow, whereas the associated SFR and kinematic properties were estimated including the region to the northwest also. We performed the two kinematic component modelling of this integrated spectrum of the outflow in the [O II] λλ3726, 29, Hβ, [O III] λ5007, Hα and [N II] λλ6548, 6584 emission lines, which is shown in Fig. 13. This spectral modelling allowed us to constrain the broad component flux in the different lines that, together with the main parameters derived for the outflow, are reported in Table 2. As can be inferred from the line fluxes, the contributions of the broad component to the total line flux in [O III] λ5007 and Hα are 23% and 33%, respectively, in agreement with the values found in ionised outflows of massive galaxies (> 10¹¹ M_⊙) at 0.6 ≤ z ≤ 2.7 (Förster Schreiber et al. 2019; Concas et al. 2022) and lower-mass (< 10¹⁰ M_⊙) systems at z = 3 − 4 (Gupta et al. 2023; Llerena et al. 2023).

Fig. 13. Spectral modelling of the main emission lines in the integrated spectrum of the region hosting the ionised outflow. The blue and red lines correspond to the narrow and broad components, respectively, which summed together (green line) provide the best fit to the data (black line).

In the BPT diagram shown in Fig. 8, the fluxes in the broad component place the ionisation of the outflowing gas (blue circle) in a region of higher ionisation than local, star-forming galaxies but compatible with ionisation by star formation at the redshift of GS_4891 (see Sect. 3.2.3). Under the assumption of star formation-driven ionisation of the outflowing gas, we can estimate its metallicity using the same prescriptions from Curti et al. (2017) as for the host galaxy (Sect. 3.2.5), in this case with the N2, R2, R3, R23, and O32 indicators. We obtain 12 + log(O/H) ∼ $8.5_{-0.1}^{+0.1}$, which is 0.2 dex higher than the one measured in the integrated spectrum of GS_4891 (Sect. 1).

The properties of the outflow and its impact on the host galaxy can be further characterised with the estimation of the outflowing mass rate, Ṁ_out (the rate at which the gas is being ejected by the outflow), and the energetics associated. The total mass of ionised gas associated with the outflow is estimated following Eq. (1) from Cresci et al. (2023), which uses the Hα luminosity of the broad component and the electron density of the outflow, N_e. In our data, the only set of emission lines that can be used to constrain the electron density in the outflow is [O II] λλ3726, 29, since the [S II] λλ6716, 31 doublet presents low S/N to be modelled with two kinematic components. We note here that the NIRSpec spectral resolution barely allows us to resolve the [O II] λλ3726, 29 and that constraining the fluxes of each kinematic component is very challenging. Therefore, instead of estimating the electron density from the broad component, we assume that the outflowing gas has the same electron density as all the gas in GS_4891 (n_e = 776 ± 307 cm⁻³) estimated in Sect. 3.2.1.

Thus, taking into account all the light within the region hosting the ionised outflow (individual spaxels and the integrated region to the northwest), we derive an outflowing gas mass, M_g = ∼1.3 × 10⁷ M_⊙. Considering the maximum velocity of the outflowing gas (v_out) and the size of the outflow (r_out), we estimate an outflowing mass rate, Ṁ_out ∼ 4 M_⊙ yr⁻¹. Based on the Hα flux and nebular extinction (1.0 ± 0.3 mag, derived as in Sect. 3.2.2) of the total (narrow + broad) line fluxes, we estimate an associated SFR, following Kennicutt (1998), of 57.6 ± 7.4 M_⊙ yr⁻¹. This SFR is higher than the one reported for GS_4891 in Table 1 (37.5 ± 0.6 M_⊙ yr⁻¹), because the area spanned by the outflow is larger than the one defined to extract the integrated spectrum from GS_4891 (Fig. 2) and now we are also including the Hα flux in the broad component. The ratio between the outflowing mass rate and the total SFR in the host region, known as the mass-loading factor (η_out = Ṁ_out/SFR), is equal to ∼0.04. This low η value indicates that the ionised outflow is not having a significant impact on the star formation activity in the host galaxy, although we have to bear in mind that we are only tracing the ionised phase of the gas and more outflow mass could be in colder and hotter gas phases (e.g. Belli et al. 2023). Our measured η is in agreement with those reported in massive (> 10¹¹ M_⊙) star-forming galaxies at 0.6 ≤ z ≤ 2.7 (η < 0.2; Förster Schreiber et al. 2019; Concas et al. 2022) but lower than in less massive (< 10¹⁰ M_⊙) systems at z = 3 − 4 (0.05 < η < 3.26; Gupta et al. 2023; Llerena et al. 2023) and at z = 3 − 9 (0.5 < η < 3.6; Carniani et al. 2024). Finally, the corresponding kinetic power of the ionised outflow is Ė_kin = 1.8 × 10⁴¹ erg s⁻¹ and the momentum rate, Ṗ_kin = 9.2 × 10³³ dyne. As was mentioned in Sect. 3.2.1, the ionised outflow might also contribute to an increase in the electron density of the gas in the galaxy, which could explain the high electron densities measured in GS_4891.

Finally, we explored the possibility that the outflowing gas might escape from the host galaxy by estimating the escape velocity of the system using the dynamical mass estimated in Sect. 3.3.1. Following Arribas et al. (2014), we computed the escape velocity at r_out = 1.5 kpc for an isothermal sphere truncated at 2.2 kpc (the radius within which the dynamical mass is computed), obtaining v_esc ∼ 240 km s⁻¹. Given that v_out can reach values up to ∼400 km s⁻¹, at least part of the outflowing gas could escape the galaxy or be re-distributed throughout the galaxy. This result would imply that outflows might also be contributing to the regulation of metals at z = 3.7. Moreover, given that the outflow is oriented towards the northwest of GS_4891, in the case in which the gas again falls on the galaxy (in the form of a galactic fountain; Fraternali & Binney 2008), it could contribute to explaining the observed positive metallicity gradient from the nucleus of GS_4891 to the northwest region (see panel c in Fig. 7). Thus, the inside-out distribution of metals via outflows would provide an alternative scenario to explain positive metallicity gradients in high-redshift galaxies, in addition to inflows of pristine gas directly into the central regions (Cresci et al. 2010; Curti et al. 2020). Nevertheless, as was discussed in Sect. 3.2.5, in the case of GS_4891 the higher metallicities at the northwest could also be a consequence of the interaction with GS_4891_n, which leads to the mixing of high-metallicity gas.

4. Summary and conclusions

In this work we present the JWST/NIRSpec IFS data of a highly star-forming galaxy group around the massive GS_4891 galaxy at z ∼ 3.7, observed as part of the GA-NIFS program. We used the low-resolution data (R100) to estimate stellar masses, whereas we used the high-resolution (R2700) observations that include a wealth of emission lines, from [O II] λλ3726, 29 to [S II] λλ6716, 31, to study with an unprecedented spatial and spectral resolution the ISM properties and ionised gas kinematics of GS_4891 and its close environment within a region of 20 × 20 kpc². The wide range of emission lines in our spectra have allowed us to study the spatially resolved distribution of the nebular extinction through the Balmer decrement, the excitation mechanisms through standard BPT diagnostic diagrams, the gas metallicity using a total of five indicators, and the ionisation parameter. We have also studied the kinematics of the galaxy group, from which we have estimated the dynamical mass of GS_4891 and detected the presence of an ionised outflow in this galaxy. Our main findings are:

• GS_4891 is part of a small group of galaxies together with one galaxy to the north, GS_4891_n, which could be interacting with GS_4891, and another one to the east, GS_28356, at a projected distance of less than 15 kpc. Additionally, GS_4891 displays evidence of internal substructure in the form of two star-forming clumps in the southeast and west that could be indicative of minor mergers and/or gas accretion.

• The nebular extinction estimates in GS_4891, computed using the Hδ, Hγ, Hβ, and Hα hydrogen Balmer emission lines, are in agreement with local attenuation curves (Cardelli et al. 1989; Calzetti et al. 2000), suggesting that the effects of dust along the line of sight at z ∼ 3.7 are similar to those in the local Universe.

• The excitation mechanisms throughout GS_4891 and GS_4891_n are consistent with star formation, with no signs of AGN activity. These systems occupy regions in the BPT and R23−O32 diagrams in between z ∼ 2 and z > 4 galaxies, consistent with an increase in the ionisation parameter and a decrease in gas metallicity with redshift.

• The gas metallicity in GS_4891_n is 0.2 dex higher than at the centre of GS_4891 despite being four times less massive. There is a drop in metallicity to the south of GS_4891, close to clump C1, which we associate with a minor merger with a lower metallicity companion or with the accretion of relatively metal-poor gas. GS_4891 presents a positive metallicity gradient (0.15 dex) towards the northwest, resembling those observed at slightly lower redshifts (Curti et al. 2020; Cresci et al. 2010) and recently also at higher redshifts (Arribas et al. 2023), which are likely associated with inflows of cold gas in the central regions. The 0.2 dex difference in metallicity between GS_4891 and GS_4891_n could be a consequence of inflows of pristine gas being more favoured in the more massive galaxy (“starvation” Larson et al. 1980; Peng et al. 2015), potentially leading to a differential build-up of the mass–metallicity relation between satellites and centrals. Assuming star formation is truncated earlier in satellites, these observations could help to explain the relative enrichment of satellite galaxies with respect to central galaxies of the same mass at low redshifts (Pasquali et al. 2012; Peng & Maiolino 2014). Preferential accretion of pristine gas by the central galaxy would also explain the positive metallicity gradient observed in GS_4891. Alternative explanations to this positive gradient include the mixing of metal-rich gas from GS_4891_n and the distribution of metal-rich gas by the ionised outflow detected in GS_4891.

• GS_4891 hosts a metal-rich, ionised outflow around the central regions that extends out to r_out = 1.5 kpc. The ionised outflow reaches maximum velocities r_out ∼ 400 km s⁻¹, yielding an outflowing mass rate, Ṁ_out ∼ 4, which is significantly smaller than the SFR in the host region, 57.6 ± 7.4 M_⊙ yr⁻¹. The corresponding low mass-loading factor, η ∼ 0.04, implies that the outflow has little effect on the SFRs of the host.

The results obtained in this work indicate that the system around GS_4891 might be in a pre-merger stage, based on the large differences in metallicities and the lack of a high gas disturbance between GS_4891 and GS_4891_n that would have been indicative of an ongoing merger. This work showcases the potential of NIRSpec IFS observations to study the internal structure and close environment of high-redshift systems, which is the main goal of the GA-NIFS survey.

Acknowledgments

We would like to acknowledge the whole JWST mission and the instrument science teams for a successful launch and commissioning of the observatory. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1216, as part of the NIRSpec Galaxy Assembly IFS GTO program. B.R.P., M.P. and S.A. acknowledge grant PID2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/10.13039/501100011033). M.P. also acknowledges the Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. I.L. acknowledges support from PID2022-140483NB-C22 funded by AEI 10.13039/501100011033 and BDC 20221289 funded by MCIN by the Recovery, Transformation and Resilience Plan from the Spanish State, and by NextGenerationEU from the European Union through the Recovery and Resilience Facility. PGP-G acknowledges support from grants PGC2018-093499-B-I00 and PID2022-139567NB-I00 funded by Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033, FEDER, UE. G.C., G.V. and M.P. acknowledge 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”. G.V. and S.C. acknowledge support from the European Union (ERC, WINGS, 101040227). A.J.B., G.C.J. and J.C. 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). H.Ü. gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. F.D.E., J.S. and R.M. acknowledge support by the Science and Technology Facilities Council (STFC), from the ERC Advanced Grant 695671 “QUENCH”. F.D.E. and J.S. also acknowledge support by the UKRI Frontier Research grant RISEandFALL. R.M. also acknowledges funding from a research professorship from the Royal Society. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2018). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work has made use of the Rainbow Cosmological Surveys Database, which is operated by the Centro de Astrobiología (CAB), CSIC-INTA, partnered with the University of California Observatories at Santa Cruz (UCO/Lick, UCSC).

References

Appendix A: 3D-Barolo fit

The best-fit model of the ionised gas kinematics in GS_4891 performed with 3D-BAROLO is shown in Figure A.1. The velocity maps can be well reproduced with a rotating disc model.

Fig. A.1. Disc kinematic best-fit of the total flux, velocity, and velocity dispersion of GS_4891 (first to third rows). The analysis was performed with 3D-BAROLO (Di Teodoro & Fraternali 2015) on the kinematics maps obtained from a spaxel-by-spaxel emission-line fitting (see Section 3.3 and Figure 11). The black and green lines identify the major axis and the zero-velocity curve, respectively.

All Tables

All Figures

	Fig. 1. False-colour image of the field around GS_4891 highlighting the position of the NIRSpec IFS field of view. Blue colours are a combination of ACS/HST images (Giavalisco et al. 2004), whereas green and red correspond to NIRCam images (Williams et al. 2023).
In the text

	Fig. 2. Maps of integrated fluxes in the emission lines [O II] λλ3726, 29, [O III] λ5007, and Hα, and the continuum around 1.6 μm for the GS_4891 system. The three main galaxies identified in the system are highlighted with red rectangles in the top-right panel. The black cross indicates the centroid of continuum emission (see text for details). North is up and east is to the left.
In the text

	Fig. 3. Flux maps in windows of two spectral pixels (∼50 km s−1 wide) centred at 25 km s−1 and −170 km s−1 with respect to the position of [O III] λ5007 at the systemic redshift (z = 3.7027) where the C1 and C2 separate clumps (magenta circles) are identified. The origin of the coordinates corresponds to RA 3h32m18.29s, Dec −27°51′59.17″, J2000. North is up and east is to the left.
In the text

	Fig. 4. Integrated spectra of GS_4891 and the two other systems, GS_4891_n and GS_28356, identified in the NIRSpec/IFS observations. The integrated spectra were obtained by adding all the spaxels within the apertures highlighted in Fig. 2.
In the text

	Fig. 5. Simultaneous spectral modelling (orange lines) of the observed [O II] λλ3726, 29 and [S II] λλ6716, 31 (black lines) to constrain the electron density of the ionised gas in the integrated spectrum of GS_4891.
In the text

	Fig. 6. AV values for the integrated spectrum of GS_4891 obtained using different pairs of hydrogen Balmer lines. The purple and orange points were obtained using the attenuation curves from Cardelli et al. (1989) and Calzetti et al. (2000), respectively.
In the text

Fig. 7. Maps of the main ISM parameters derived for the GS_4891 galaxy group from emission line analysis. (a) Nebular extinction, AV, derived from the Balmer decrement (Hα/Hβ). (b) SFRs (M⊙ yr−1) estimated from the Hα luminosity. (c) Gas oxygen abundance obtained with the R2, R3, R23, and O32 metallicity calibrations from Curti et al. (2017). (d) Ionisation parameter, U, computed with the gas metallicity and the [O II] λλ3726, 29 lines following Eq. (5) in Díaz et al. (2000). The data were spatially binned to have a S/N ≥ 3 in the emission lines involved. The dashed black contours correspond to the Hα emission, whereas the black cross indicates the centroid of continuum emission.

In the text

Fig. 8. Left: map of excitation mechanisms throughout GS_4891 and GS_4891_n according to the location in the [N II]-based BPT diagnostic diagram. The dashed black contours correspond to the Hα emission, whereas the black cross indicates the centroid of continuum emission. Centre and right: [N II]- and [S II]-based BPT diagrams containing the spatially binned regions (red dots), the integrated regions of the system where the relevant lines are detected with S/N ≥ 3 (red square and diamond, and pink circle), and the broad kinematic component associated with the ionised outflow identified in GS_4891 (blue circle). Clump C2 is not shown in any of the diagrams due to the low S/N in Hβ and [N II]. As a reference, we have added the stack of galaxies at z ∼ 6 from Cameron et al. (2023, purple plus), the stacks in different high-redshift bins from Sanders et al. (2023, yellow symbols), the obscured AGN at z ∼ 3.3, “Jil2”, from Perna et al. (2023, lightgreen star), and the broad-line AGN, GS_3073, at z ∼ 5.5 studied in Übler et al. (2023, darkgreen star). The shaded grey contours demarcate the regions in the diagram encompassed by 70, 80, and 90% of the local SDSS sample (Abazajian et al. 2009). Dashed (Kauffmann et al. 2003) and dotted black lines (Kewley et al. 2001) demarcate the separation between ionisation coming from star formation and AGN activity based on studies of local galaxies (see Sect. 3.2.3 for details). The cyan lines correspond to the best-fit curves to the star-forming population at z ∼ 2.3 from Strom et al. (2017).

In the text

	Fig. 9. Maps of emission line ratios that trace gas metallicity throughout the GS_4891 system. The data have been spatially binned to have a S/N ≥ 3 in the emission lines involved. The black cross indicates the centroid of continuum emission, whereas dashed contours in the right panel correspond to the Hα emission. Black boxes and circles highlight the individual systems identified in Sect. 3.1.
In the text

Fig. 10. R23−O32 diagram for the GS_4891 system in comparison with other works. Spatially binned regions are shown as red dots, whereas integrated values are shown as red polygons and empty and filled magenta circles. As a reference, we have added the stacks of galaxies at z ∼ 2.5, 3.5 from Kriek et al. (2015, light and dark blue points, respectively), the stacks in different redshift bins from Sanders et al. (2023, yellow symbols), and the individual galaxies at z ∼ 6, 8 from Cameron et al. (2023, purple and brown pluses, respectively). Shaded contours demarcate the regions in the diagram encompassed by 70, 80, and 90% of the local SDSS sample (Abazajian et al. 2009).

In the text

	Fig. 11. Velocity maps of the ionised gas in the GS_4891 system obtained from the fitting of the [O III] λ5007 and Hα lines. The three left panels show the velocities at the 2nd, 50th, and 98th percentiles of the velocity distributions. The right panel shows the W80 map. The black cross indicates the centroid of continuum emission, whereas the dashed black contours correspond to the Hα emission.
In the text

	Fig. 12. Velocity and velocity dispersion maps of the broad kinematic component superimposed on the Hα flux distribution in GS_4891 (grey background). The black cross indicates the centroid of continuum emission.
In the text

	Fig. 13. Spectral modelling of the main emission lines in the integrated spectrum of the region hosting the ionised outflow. The blue and red lines correspond to the narrow and broad components, respectively, which summed together (green line) provide the best fit to the data (black line).
In the text

	Fig. A.1. Disc kinematic best-fit of the total flux, velocity, and velocity dispersion of GS_4891 (first to third rows). The analysis was performed with 3D-BAROLO (Di Teodoro & Fraternali 2015) on the kinematics maps obtained from a spaxel-by-spaxel emission-line fitting (see Section 3.3 and Figure 11). The black and green lines identify the major axis and the zero-velocity curve, respectively.
In the text