© The Authors 2024

1. Introduction

Observations have revealed that the mode of galaxy evolution in the first 2 Gyr of the Universe (z ≳ 3) was drastically different from that of following epochs. This can be seen in the similar evolution of the global star formation rate (SFR) density and molecular gas mass, both of which increased with cosmic time during this epoch and decreased for z < 1 (e.g. Bouwens et al. 2015; Decarli et al. 2019). At higher redshifts, the major merger rate is higher (e.g. Duncan et al. 2019), as are the mean SFR and black hole accretion rate for a given stellar mass (e.g. Speagle et al. 2014; Yang et al. 2018; Popesso et al. 2023). Together, these results suggest that galaxies in the first 2 Gyr formed rapidly, accreting gas via filaments and mergers, which resulted in the buildup of stellar and black hole masses on shorter timescales than in the following epochs.

The unique environment of the early Universe allowed galaxies to start forming stars quickly, resulting in high-redshift galaxies that exhausted their gas supply and stopped forming stars by z ≳ 3 (e.g. see the recent JWST results of Looser et al. 2023; Strait et al. 2023; Carnall et al. 2023). While this quenched appearance could be caused by a temporary minimum in a stochastic star formation history (e.g. Arata et al. 2020; Dome et al. 2024), there have also been a number of z > 4 hyper-luminous infrared galaxies (HyLIRGs) detected with large infrared luminosities (L_IR > 10¹³ L_⊙) and SFRs of ≳10³ M_⊙ yr⁻¹ (e.g. Wagg et al. 2014; Venemans et al. 2019; Carniani et al. 2019; Riechers et al. 2020; Cheng et al. 2020) that could deplete gas reservoirs rapidly. The elevated SFR of these sources was believed to be the result of ongoing hierarchical merging (e.g. Hopkins et al. 2006), although the discovery of ordered rotation in some HyLIRGs (e.g. Tsukui & Iguchi 2021) suggests that they could be fuelled by secular accretion.

One of the most extreme high-redshift HyLIRGs is 1HERMES S350 J170647.8+584623, or HFLS3. This source lies in the Spitzer First Look Survey (FLS; Fadda et al. 2004) extragalactic field (a four square degree field centred on 17h18m00s +59°30′00.0″) and was included in the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012). Based on very red Herschel Spectral and Photometric Imaging Receiver (SPIRE) colours (S​_500 μm > S​_350 μm > S​_250 μm and S​_500 μm/S​_350 μm > 1.3), Riechers et al. (2013) identified this source as a high-redshift dusty galaxy with a nearby low-redshift companion to the north (G1B, z ∼ 2.092 based on the detection of CIVλ1549 with Keck/LRIS). A comprehensive suite of observations – including the Plateau de Bure Interferometer (PdBI), the Very Large Array (VLA), Keck, and the William Herschel Telescope – enabled the creation of a near-infrared/radio spectral energy distribution (SED) and the measurement of a number of emission lines for HFLS3 (z_spec = 6.3369 ± 0.0009). Based on their analysis of the line and continuum emission, Riechers et al. (2013) report a high mass of molecular gas (M_H2 = (1.04 ± 0.09)×10¹¹/α_CO M_⊙), dust ($M_{\mathrm{dust}}=1.{31}_{-0.30}^{+0.32}\times {10}^9{M}_{\odot }$), and stars (M_* ∼ 3.7 × 10¹⁰ M_⊙), as well as a high far-infrared (FIR) luminosity ($L_{\mathrm{FIR}}=2.{86}_{-0.31}^{+0.32}\times {10}^{13}$ L_⊙), a high SFR_FIR of ∼2900 M_⊙ yr⁻¹, and evidence for a velocity gradient from PdBI [CII]158 μm¹. Some of the spectral lines were asymmetric, which Riechers et al. (2013) interpret as evidence for a possible close separation merger. However, some of the data were taken at low resolution (> 2″) and thus did not allow for detailed source differentiation and characterisation.

Follow-up Keck/NIRC2 and Hubble Space Telescope (HST) Wide-Field Camera 3 (WFC3) imaging by Cooray et al. (2014) revealed the presence of three close-by companions: two to the north (G1 and G2) and one to the south (R1). The northern sources were assumed to be at the same redshift (z ∼ 2.0), while a photometric redshift of z_phot ∼ 6 was found for R1. Gravitational lens modelling suggests that HFLS3 may be composed of two sources that are slightly magnified (μ ∼ 2.2 ± 0.3) by the foreground sources G1 and G2.

The idea that the starburst nature of HFLS3 could be caused by a history of major mergers in a dense environment inspired multiple searches for an overdensity of galaxies. Robson et al. (2014) used the Submillimetre Common-User Bolometer Array 2 (SCUBA-2) on the James Clerk Maxwell Telescope (JCMT) to perform a low-resolution (∼14″), wide-area (67 arcmin²) search for significant sub-millimetre emission at two wavelengths (450 μm and 850 μm), but found no overdensity. This was followed by a search at optical wavelengths with HST and the Gran Telescopio Canarias (GTC; Laporte et al. 2015); no large-scale overdensity was found, but a possible small-scale (∼36 kpc) overdensity of faint objects was.

These observations suggest that HFLS3 is a gas-rich, starbursting, dusty galaxy that had a velocity gradient (suggesting merging activity) when the Universe was only ∼850 Myr old. Because of these exceptional properties, it was chosen as a target for the JWST Near Infrared Spectrograph (NIRSpec) guaranteed time observation (GTO) programme Galaxy Assembly with NIRSpec IFS (GA-NIFS; PI: S. Arribas & R. Maiolino). This programme aims to observe (in cycles 1 and 3) a sample of 55 galaxies at z ∼ 3 − 11 spanning a variety of types (e.g. quasi-stellar objects, active galactic nuclei, star-forming galaxies, strongly lensed galaxies, quenched sources, major mergers), to show the power of the integral field unit (IFU; Böker et al. 2022) on NIRSpec (Jakobsen et al. 2022) for exploring resolved kinematics and gas properties. The survey is ongoing, and detailed studies of z = 3 − 7 active galactic nuclei (AGN) and quasi-stellar objects have been published (Marshall et al. 2023; Übler et al. 2023; Perna et al. 2023).

Here we present the JWST/NIRSpec IFU observations of HFLS3. The high resolution and sensitivity of these observations reveal a more complex system than implied by previous data, with six strongly detected sources at z ≃ 6.3 − 6.5 and two low-redshift interlopers (z ∼ 2.0 − 3.5) within a 3″ × 3″ field. We present the details of our dataset in Sect. 2 and characterise the field in Sect. 3. Section 4 contains further analysis (i.e. morpho-kinematic maps, gravitational lens modelling, line ratio-based excitation conditions, SFR derivation, and galaxy merger discussion), and we conclude in Sect. 5.

We use a standard concordance cosmology (h_o, Ω_m, Ω_Λ = 0.7, 0.3, 0.7) throughout, where 1″ corresponds to ∼5.5 kpc at z ∼ 6.34, ∼7.3 kpc at z ∼ 3.48, and ∼8.4 kpc at z ∼ 2.0. To match the notation of other works (e.g. Curti et al. 2020), emission lines are named based on their air wavelength, while we use their vacuum wavelengths for analysis (e.g. λ_[OIII]λ5007 = 5008.24 Å). We adopt a Salpeter stellar initial mass function (IMF; Salpeter 1955).

2. JWST/NIRSpec IFU data description

HFLS3 was observed with the JWST/NIRSpec IFU using two disperser/filter combinations (G395H/290LP, hereafter R2700; and PRISM/CLEAR, hereafter R100) on 1 September 2022 as part of Program ID 1264 (PI: L. Colina; see Table 1). Both settings used the improved reference sampling and subtraction pattern IRS² of Rauscher et al. (2017) and adopted a 4-point dither with a medium “cycling” pattern and a starting point of 1.

We used v1.8.2 of the GTO pipeline under the Calibration Reference Data System (CRDS) context jwst_1105.pmap to create R2700 and R100 cubes with drizzle weighting (Fruchter & Hook 2002). This context accounts for flux leakage from 68 failed open shutters in the multi-shutter array (MSA; ∼0.03% of the full MSA) into the IFU aperture. A patch was included to correct some important bugs affecting this specific pipeline version (see details in Perna et al. 2023). We corrected count-rate images for 1/f noise via a polynomial fit. During stage 2, we removed all data in regions of known failed open MSA shutters. We also masked pixels at the edge of the slices (one pixel wide) to conservatively exclude pixels with unreliable flat-field reference files for the spectrograph optics (i.e. SFLAT files), and implemented the outlier rejection of D’Eugenio et al. (2023).

A recent in-depth investigation of the NIRSpec IFU point spread function (PSF) by D’Eugenio et al. (2023) found that it may be non-circular. Using a standard star, they determine that the major axis full width at half maximum (FWHM) is wavelength dependent (FWHM_major ∼ 0.12 − 0.20″ between λ_obs = 3.0 − 5.0 μm). The combination of a dither and drizzle weighting allowed us to sub-sample the detector pixels, resulting in cube spaxels of 0.05″. The resolving power of the R2700 data increases from R₂₇₀₀ ∼ 2000 − 3500 from λ_obs = 3.0 − 5.0 μm, while the resolving power of the R100 data decreases from R₁₀₀ ∼ 90 − 30 over λ_obs = 0.5 − 1.2 μm and increases from R₁₀₀ ∼ 30 − 290 over λ_obs = 1.2 − 5.0 μm.

2.1. Astrometric correction

A comparison of the JWST/NIRSpec IFU data with archival HST data² revealed a positional offset between images at comparable wavelengths. Astrometric errors like this are currently common in NIRSpec IFU data at the time of analysis (e.g. Wylezalek et al. 2022; Perna et al. 2023). Since we wish to estimate spatially resolved quantities, it is crucial to have each image aligned to a common reference frame. By shifting to the Gaia Data Release 3 (DR3) frame, we find an offset of 0.17 ± 0.07″ (see Appendix A.1 for details).

2.2. Background subtraction

Since these observations did not include a dedicated background exposure, no background subtraction was performed by the pipeline. To estimate the background emission, we extracted a mean spectrum from the R100 and R2700 data cubes using a 25 spaxel (12.5″) diameter aperture (MPDAF task aperture; Bacon et al. 2016) from a signal-free region to the south-east.

We assumed that the background emission is uniform across the field of view (FoV) and created a correction by subtracting this background spectrum from the spectrum of each spaxel. We note that for the R2700 spectrum, we only included channels that are not affected by chip gap or edge issues (λ_obs ≈ 2.87057 − 3.9778, 4.2078 − 5.26643 μm).

3. Source distribution

Previous analyses of the HFLS3 field found evidence for multiple sources (Riechers et al. 2013; Cooray et al. 2014): the primary lensed starburst (HFLS3, z ∼ 6.34), two low-redshift sources to the north (G1 and G2, reported to be z ∼ 2.0), and a fourth source only detected in rest-frame UV emission to the south-east (R1). Initial exploration of the IFU data cubes shows evidence for a complex distribution of flux with multiple spatial components. It is not clear a priori if these are the same as previously detected, so we proceeded with an uninformed search for emission.

3.1. Field characterisation

As an example of the complex source distribution, we present an integrated map of the R2700 cube over the wavelength range that contains Hα redshifted to z = 6.34 (i.e. λ_obs = 4.79954 − 4.84467 μm, or z_Hα ∼ 6.31 − 6.38; upper panel in Fig. 1). There are four separate areas of emission: a bright source to the north, an east–west arc near the centre, a southern source, and a double-lobed source to the west.

Fig. 1. Integrated emission of the HFLS3 field, using two illustrative wavelength ranges: redshifted Hα for z ∼ 6.34 for the R2700 cube (λobs = 4.79954 − 4.84467 μm, upper panel) and the approximate wavelength range of HST/WFC3 F160W for the R100 cube (λobs = 1.4 − 1.6 μm; lower panel). The adopted source masks are shown with white boundaries. North is up and east is to the left. The lower panel is nearly identical to the observed HST/WFC3 F160W image, but at a higher spatial resolution (see Fig. A.2).

If we instead integrate the R100 cube over a wavelength range similar to the HST/WFC3 F160W transmission (i.e. λ_obs = 1.4 − 1.6 μm), the distribution is drastically changed (Fig. 1, lower panel). We note that this distribution is identical to that seen in previous HST/WFC3 F160W observations (see Fig. A.2), but at higher spatial resolution. The northern source is still present, and there is an extended object with a bright core between this source and the central arc (as seen in the HST images; Cooray et al. 2014). The other three sets of sources are much weaker.

From these maps, we defined five regions (see Fig. 1): the bright northern component (G1), the core of the extended object in the R100 data (G2), the arc (C), the southern sources (S), and the western sources (W). The regions G1 and G2 are identical to those of Cooray et al. (2014), while the C region contains the HFLS3 source. The previously found source R1 is the brightest region within the S mask, but we used a larger area to encompass nearby emission. Laporte et al. (2015) identified the S and W sources as possible faint sources (ID2 and ID3, respectively) in HST maps. We also note that weaker line emission is detected in the FoV, which we explore in Appendix B.

3.2. R2700 integrated spectra

To investigate the distribution of flux in this field, we extracted spectra from the R2700 cube using each of the spatial masks (see the top panel of Fig. 1) except G2, which is examined in the next subsection. These spectra are fit with multiple pseudo-Voigt profiles (one pseudo-Voigt profile per line). This profile is an adjustable combination of a Gaussian and Lorentzian, which represent line broadening due to pressure and Doppler (or thermal) motion, respectively (Armstrong 1967; see Appendix C.1 for profile details). All lines are fit simultaneously using LMFIT with a least-squared minimisation, assuming a single redshift for most lines and a power law continuum normalised to the 4 μm value:

$$\begin{array}{c}\hfill F({\lambda }_{\mathrm{obs}})=F_{4\mathrm{\mu }\mathrm{m}}{\left(\frac{{\lambda }_{\mathrm{obs}}}{4\mathrm{\mu }\mathrm{m}}\right)}^{\alpha }·\end{array}$$(1)

The models are convolved with the line spread function (LSF)³ before comparison to the data. The velocity widths of line pairs (i.e. [SII]λλ6716, 6731, [NII]λλ6548, 6584, [OIII]λλ4959, 5007, and each Balmer line) are fixed to be identical. We adopted the standard assumptions of [NII]λ6584/[NII]λ6548 = 2.94 (e.g. Dojčinović et al. 2023) and [OIII]λ5007/[OIII]λ4959 = 2.98 (e.g. Dimitrijević et al. 2007). Following other NIRSpec IFU investigations (Übler et al. 2023; Perna et al. 2023), we rescaled the error spectrum generated by the calibration pipeline (“ERR”) using the standard deviation of the line-free regions in the observed spectrum. Uncertainties on each parameter were estimated using the standard error output from LMFIT. Only lines with > 1σ emission at the expected centroid wavelength are fit.

The best-fit spectra are shown in Fig. 2. The redshifts and continuum fluxes for each source are presented in Table 2, and the line properties are listed in Table 3.

Fig. 2. Integrated spectra of the R2700 cube using the masks of Fig. 1 (magenta; excluding G2), with 1σ errors from the associated error spectrum shown as the shaded region. Best-fit models (line emission and continuum) are shown by green lines. The centroids of each line are depicted by dashed lines, with red lines indicating that the spectral line was not fit. The lower panel shows the residual.

Fig. 2. continued.

Fig. 2. continued.

With a redshift of z_G1 = 3.4805 ± 0.0001, the northern galaxy G1 is well detected in Hα, [NII]λλ6548, 6584, [SII]λλ6716, 6731, [SIII]λ9532, and HeIλ10829. Because HeIλ10829 may be resonantly scattered (e.g. Rudy et al. 1989; Belli et al. 2023), we allowed the centroid of this line to vary (resulting in an offset of 137 ± 19 km s⁻¹). In addition, [SIII]λ9532 may be blended with Pa-ϵ (e.g. Kehrig et al. 2006), so we allowed the model of this line to feature a velocity offset (54 ± 14 km s⁻¹). The best-fitting model shows minor (∼1σ) residuals for the Hα complex, possibly suggesting the presence of an additional galaxy or outflow.

The central component C (z_C = 6.3425 ± 0.0002) features well-detected Hα, Hβ, [OIII]λ5007, and [NII]λ6584 emission with broad profiles (FWHM ∼ 500 − 700 km s⁻¹). We find evidence for weak [OIII]λ4959, [NII]λ6548, [SII]λλ6716, 6731, with no Hγ. The [SII]λλ6716, 6731 emission is broad and low-level, so fits with a free FWHM returned unphysical values (∼1500 km s⁻¹). To better constrain these lines, we fixed the width of each to be identical to that of Hα. The 4 μm continuum is strong but features no significant slope in F_λ. While the best-fit model spectrum does not show significant residuals (Fig. 2), we note that this component is found to feature two discrete sources with a velocity offset (see Sect. 4.2.3) with a velocity separation smaller than our LSF. With this in mind, the large linewidth is due to the blending of the two galaxies.

The integrated spectra of the S and W components are quite similar to each other, with very strong detections in Hγ, Hβ, [OIII]λλ4959, 5007, and Hα. But while the S emission is fit well by a single redshift (z_S = 6.3592 ± 0.0001), the double-peaked appearance of the W component suggests two galaxies (z_W1 = 6.3549 ± 0.0001 and z_W2 = 6.3627 ± 0.0001, with a velocity difference of 368 ± 7 km s⁻¹) that may each be characterised directly from the integrated spectrum. No significant [NII]λλ6548, 6584 or [SII]λλ6716, 6731 emission is detected in the integrated W spectrum.

While component S is fit well by a single velocity component, we find that it is composed of two spatially separated sources at the same redshift (see Sect. 4.1). To investigate each component (S1 and S2), we created two sub-masks (see the dashed line in Fig. 1) and fitted the integrated spectra of each. S1 is brighter in line emission (with significant [NII]λλ6548, 6584 emission), but the continuum value of each component is approximately equal.

3.3. R100 integrated spectra

Due to the coarse LSF of the R100 data (∼10⁴ km s⁻¹; Böker et al. 2023), lines are frequently blended. In order to fit models to these data, one must create a high-resolution model and convolve it with the LSF (e.g. Heintz et al. 2023; Jones et al. 2024; Umeda et al. 2023). While the R100 data are therefore not useful for kinematics or precise redshift derivations, we can examine the spectrum of G2 to find a general redshift.

G2 features two clear emission lines in a wavelength regime (λ_obs = 1.4 − 2.0 μm) that is not covered by the R2700 data. To fit this, we first created a higher-resolution model spectrum (δλ = 0.001 μm) that is populated with the strongest expected lines (i.e. Hβ, [OIII]λλ4959, 5007, [NII]λλ6548, 6584, and Hα)⁴. Since the lines are unresolved, we adopted Gaussian profiles instead of the more detailed pseudo-Voigt profiles used in the R2700 analysis. We find evidence for a spectral break at λ_obs ∼ 1.1 μm, which we tentatively interpret as a Balmer break. The continuum is thus modelled as a power law redwards of this break (λ_rest = 364.5 nm) and a separate power law bluewards of this value. The combined line + continuum model is then convolved with the resolving power to account for the LSF, and is re-binned to the R100 spectral bins. We used LMFIT with a least-square minimisation to find the best-fit continuum and line parameters and redshift, resulting in the fit of Fig. 3. The error spectrum was rescaled to match the standard deviation of the residual spectrum.

Fig. 3. Integrated spectrum extracted from the R100 cube using the G2 mask of Fig. 1, with illustrative 1σ errors shown as shaded region. A basic fit to the data is shown by the green line. The lower panel shows the residual, with the centroids of each line depicted by dashed lines. From left to right, the lines are Hβ, [OIII]λλ4959, 5007, [NII]λ6548, Hα, and [NII]λ6584. The best-fit redshift is zG2 = 2.00 ± 0.01.

Due to the heavy blending of the lines, we do not present best-fit widths or integrated fluxes. However, we find a best-fit redshift of z = 2.00 ± 0.01, which is comparable to the CIV-based redshift of Riechers et al. (2013; z ∼ 2.092). The redshift difference may be ascribed to a different extraction aperture or the simplicity of our current model. Indeed, the best-fit model features residuals around the putative Balmer break and around a possible line at λ_obs ∼ 1.72 μm. While the line amplitude is quite weak (∼1σ), the centroid wavelength of this line is in agreement with the faint auroral line [NII]λ5755.

We note that the rich R100 spectrum of each source (which contain more emission lines) will be modelled with advanced codes in future works, and our current model is solely used to confirm the approximate redshift of the G2 object.

4. Analysis

4.1. Morpho-kinematic maps

Next, we investigated the morpho-kinematics of these components using the R2700 data. This is commonly done using non-parametric measures (e.g. cumulative velocity distributions or moment maps), which are well suited to sources with isolated, asymmetric lines. Since the emission lines of our sources feature overlapping lines in some spectra (e.g. Hα and [NII]λλ6548, 6584) we instead used a parametric model. The standard assumption of a Gaussian model resulted in some low-level residuals, so we instead assumed that each line may be modelled by a pseudo-Voigt profile (see details in Appendix C.2).

The resulting integrated intensity, velocity offset, velocity dispersion, and continuum maps are shown in Figs. 4–7. Some weaker lines (e.g. [SII]λ6731 in G1) are detected in the integrated spectrum, but lack the S/N to be significantly detected in the spectrum of multiple individual spaxels.

Fig. 4. Morpho-kinematic maps of the strong emission lines detected in the component G1 (z = 3.4806). Integrated intensities are given per spaxel. The corresponding maps of [SII]λλ6716, 6731 are not presented due to low S/N. Due to the construction of our model, [NII]λ6548 has the same integrated velocity field, velocity dispersion map, and continuum map as [NII]λ6584, but an integrated intensity map that is a factor of 2.98 lower. For [SIII]λ9532 and HeIλ10829, v50 values are given with respect to their best-fit redshifts (z = 3.4810 and 3.4821, respectively).

Fig. 5. Morpho-kinematic maps of the strong emission lines detected in the C component (z = 6.3425). Integrated intensities are given per spaxel. The corresponding maps of Hβ, [OIII]λλ4959, 5007, and [SII]λλ6716, 6731 are not presented due to low S/N. Due to the construction of our model, [NII]λ6548 has the same integrated velocity field, velocity dispersion map, and continuum map as [NII]λ6584, but an integrated intensity map that is a factor of 2.98 lower.

Fig. 6. Morpho-kinematic maps of the strong emission lines detected in the S components (z ∼ 6.3592). Integrated intensities are given per spaxel. The corresponding maps of Hγ are not presented due to low S/N. Due to the construction of our model, [OIII]λ4959 has the same velocity field, velocity dispersion map, and continuum map as [OIII]λ5007, but an integrated intensity map that is a factor of 2.94 lower.

Fig. 7. Morpho-kinematic maps of the strong emission lines detected in the W components (zW1 = 6.3550, zW2 = 6.3628). Integrated intensities are given per spaxel. The corresponding maps of Hγ are not presented due to low S/N. Due to the construction of our model, [OIII]λ4959 has the same velocity field, velocity dispersion map, and continuum map as [OIII]λ5007, but an integrated intensity map that is a factor of 2.94 lower.

The maps for component G1 are shown in Fig. 4. For the two strongest lines (Hα and [NII]λ6584), we see that the emission is circular with a slight concentration in the north-east–south-west diagonal. There is no strong velocity gradient, and both the velocity dispersion and continuum emission feature a central peak. The two lower-S/N lines ([SIII]λ9532 and HeIλ10829) are detected in the core, with a similar distribution to the core of the bright lines. This elongation may be influenced by a non-circular PSF (e.g. D’Eugenio et al. 2023). These combined morpho-kinematics may be interpreted either as a face-on disk or a dispersion-dominated galaxy.

From the maps of component C (Fig. 5), an east–west velocity gradient is apparent. This has previously been seen in PdBI [CII] observations (Riechers et al. 2013), and the higher spatial resolution of our data allows for a more in-depth investigation. We find that while the continuum emission peaks in the centre (possibly with a contribution from G2; see Fig. 1), the line intensity peaks in the east and west sides. In addition, the velocity dispersion does not feature a central maximum. This argues for the presence of two separate galaxies, which we investigate further in Sect. 4.2.

The S component is clearly composed of two spatially separated galaxies at a similar redshift (i.e. low velocity dispersion, no apparent velocity gradient; Fig. 6) with continuum and Hα intensity minima between the galaxies. The brighter galaxy (to the south-east; S1) was originally identified as R1 by Cooray et al. (2014).

The W component features a strong velocity gradient (Fig. 7). However, the spatially offset integrated intensity peak (to the north) and higher velocity dispersion in the south argues against rotation. Instead, this appears to be two galaxies that are spatially and spectrally separated. This is supported by the asymmetric double-peak nature seen in the integrated spectrum of this component (Fig. 2) and the two spatial peaks seen in Fig. 1.

Altogether, our analysis suggests that the HFLS3 field is composed of several components: two low-redshift sources to the north (G1, z ∼ 3.481; G2, z ∼ 2.00), a lensed source with complex kinematics suggesting two components (C1 and C2, z ∼ 6.342), two galaxies to the south (S1 and S2, z ∼ 6.359), and two galaxies to the west (W1, z ∼ 6.355; W2, z ∼ 6.363).

4.2. Gravitational lens modelling

The gravitational lens model of the HFLS3 field was previously analysed by Cooray et al. (2014), who assumed that G1 and G2 were at the same redshift (z ∼ 2.0) and used the marginally resolved PdBI [CII]158 μm image of HFLS3 to derive a low magnification factor (μ ∼ 2.2). With the precise spectroscopic redshifts of G1 and G2 and a resolved map of HFLS3, we refined this model using the public lens modelling software PYAUTOLENS⁵ (Nightingale & Dye 2015; Nightingale et al. 2018, 2021).

4.2.1. Methods

Using PYAUTOLENS, we could derive the intrinsic (source-plane) mass and light profiles for each source in the FoV. Each light profile is given by a Sérsic profile (Sérsic 1963), while we assumed elliptical isothermal mass profiles. The mass and light profiles of each object are assumed to have shared centres, axis ratios, and position angles. We modelled the PSF as a circular Gaussian with FWHM ∼ 0.1″. The fitting process results in the best-fit centroid, intrinsic axis ratio (q ≡ a/b), and position angle (ϕ) of the mass and light profiles, the effective radius (r_eff) and Sérsic index (n) of the light profile, and the Einstein radius (r_Ein) of the mass profile. By taking the ratio of the image-plane and source-plane fluxes of the best-fit model for each component, we calculated a total magnification factor.

Due to the complexity of the field, we began by examining the two low-redshift sources (G1 and G2). These are isolated by collapsing the R100 data cube over λ_obs = 0.8 − 1.1 μm (which is not covered by the R2700 cube). This collapsed image (upper-left panel of Fig. 8) has little contribution from the z > 6 components. Because of this, we could examine how G2 (z ∼ 2.0, the northernmost source in the top row of Fig. 8) lenses G1 (z ∼ 3.5, the elongated central source in the top row of Fig. 8). Our model contains the mass and light profile of G2 and the light profile of G1. The resulting fit is shown in the top row of Fig. 8, with parameters listed in Table 4 and source-plane models of G1 and G2 given in Appendix D.

Fig. 8. Results of PYAUTOLENS lensing analysis of collapsed images of the R100 cube (λobs = 0.8 − 1.1 μm; top row) and the R2700 cube (λobs = 4.80 − 4.85 μm) assuming that component C is composed of a single source (middle row) or two sources (C1 to the east and C2 to the west; lower row). The observed data, best-fit model, and residual (all in the image plane) are presented from left to right. Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.

We next turned to a wavelength range in the R2700 cube containing Hα at z ∼ 6.34 (λ_obs = 4.80 − 4.85 μm) to examine whether component C is better modelled by a single or double source. This wavelength range contains strong Hα emission from component C, moderate continuum from G1, and no strong emission from G2. We modelled the light and mass profile of G1, assuming the best-fit mass model for G2 from the R100 collapsed image but with a variable Einstein radius. Component C is fit using either a single Sérsic profile or two spatially offset Sérsic profiles. The resulting fits for a single and double component for C are shown in the middle and lower row of Fig. 8, with parameters listed in Table 4.

4.2.2. Results

From this analysis, we find that the lowest-redshift source (G2) is elongated (q ∼ 2.4), compact (r_eff ∼ 0.25″, or ∼2.1 kpc at z = 2.00), and not centrally peaked (n ∼ 1.2). On the other hand, the source-plane morphology of G1 is characterised as elongated (q ∼ 1.4 − 1.9), extended (r_eff ∼ 1.0 − 2.7″, or ∼7.3 − 19.8 kpc at z = 3.48), centrally peaked (n ∼ 3 − 5), and only moderately magnified by G2 (μ ∼ 1.1 − 1.2).

Component C in the collapsed R2700 image appears as a curved arc with two bright regions. We first examined whether this may be explained as a single lensed source, yielding a best-fit source-plane morphology that is elongated (q ∼ 2.3), compact (r_eff ∼ 0.3″, or ∼1.6 kpc at z = 6.34), with an exponential profile (n ∼ 1.2), and a low magnification factor (μ = 1.89). As seen in the middle row of Fig. 8, the image-plane morphology of this source is indeed arced, but has a central peak rather than two spatially separate clumps. This results in a significant residual at the location of the western spot in the observed arc.

Next, we assumed that component C is composed of two spatially separate Sérsic profiles (C1 to the east and C2 to the west). The best-fit model returns one source (C1) that is similar to the similar to the best-fit single-component model, with q ∼ 2.0, r_eff ∼ 0.3″, n ∼ 1.2, and a higher magnification factor (μ = 2.09). However, this model contains an additional component to the west that is more circular (q ∼ 1.5), smaller (r_eff ∼ 0.2″), more centrally peaked (n ∼ 3), and is more magnified (μ = 2.43). The resulting image-plane morphology of this best-fit model is double-peaked, resulting in a better representation (i.e. a lower χ²). The two best-fit components have a small source-plane projected distance (0.2″, or ∼1.1 kpc at z = 6.34) that is comparable to their effective radius. This is discussed in the next subsection.

We note that this analysis suggests that none of the emission peaks in the observed fields represent multiply imaged objects. The z < 6 sources (G1 and G2) feature significantly different redshifts and are quite compact, implying that they are individual galaxies. The C emission is better fit by two separate components, rather than a single galaxy imaged twice. The effects of the gravitational lens are quite small at larger spatial separations from G2, so it is unlikely that the S and W components are magnified or multiply imaged (as supported by the different morpho-kinematics of each sub-component). If the PSF is truly larger than 0.1″ (see Sect. 2), then our general conclusions are unchanged, although C1 and C2 feature smaller best-fit source-plane effective radii with higher magnification.

4.2.3. Nature of component C

This lensing analysis of the Hα emission in the central component suggests that it is composed of two spatially separate (∼1 kpc projected offset) galaxies. By examining the best-fit morphologies (Fig. D.3), it is clear that their proximity to each other results in overlapping emission in the image-plane. So while the far east and west sides of the central arc in the image plane are dominated by the emission from each source, the space between the two apparent galaxies contains emission from both sources.

Returning to the morpho-kinematic maps (Fig. 5), we see that each line shows a weak (δv ∼ 100 km s⁻¹) east–west velocity gradient (as seen in previous [CII] observations; Riechers et al. 2013). Since none of the integrated intensity or velocity dispersion maps show a central peak, as expected for a single rotating disk (see the best-fit single component model in Fig. 8), these maps instead imply multiple galaxies separated in velocity.

This is supported by the previous morpho-kinematic investigation of FIR lines by Riechers et al. (2013). An east–west velocity gradient was found in [CII], but with a high velocity dispersion that peaked in the south-east. In addition, the integrated spectrum exhibited two velocity components in multiple FIR lines (z = 6.3335 and z = 6.3427, resulting in a weighted average of z = 6.3369 ± 0.0009). From our best-fit redshift (z_C = 6.3425 ± 0.0003), this suggests that the rest-frame optical emission is dominated by the higher redshift component. Due to the overlapping image-plane distributions of the two components (see the bottom-centre and bottom-right panels of Fig. D.3), they are not easily separable without lens modelling.

To summarise, the lensing and morpho-kinematic analyses of our new high-resolution data agree with previous results (Riechers et al. 2013; Cooray et al. 2014), who claimed that component C of HFLS3 is composed of two sub-components. The best-fit magnification of our model is also similar (μ ∼ 2.1 − 2.4, compared to μ = 2.2 ± 0.3).

4.3. Excitation conditions

Using the well-detected lines for our sample, we could place constraints on the excitation conditions of each source using line ratio diagnostic diagrams. To do this, we first calculated three line ratios:

$$\begin{array}{cc}& \mathrm{N}2={log}_{10}([\mathrm{NII}]\lambda 6584/\mathrm{H}\alpha )\hfill \end{array}$$(2)

$$\begin{array}{cc}& \mathrm{S}2={log}_{10}([\mathrm{SII}]\lambda \lambda 6716,6731/\mathrm{H}\alpha )\hfill \end{array}$$(3)

$$\begin{array}{cc}& \mathrm{R}3={log}_{10}([\mathrm{OIII}]\lambda 5007/\mathrm{H}\beta ).\hfill \end{array}$$(4)

These are used to create the [NII]-BPT (R3 vs. N2; Baldwin et al. 1981) and [SII]-VO87 (R3 vs. S2; Veilleux & Osterbrock 1987) plots (see Fig. 9) for each high-redshift (z > 6) source in the field. We examined the line ratios of both galaxies in component S (S1 and S2) and W (W1 and W2) separately. On the other hand, the galaxies in C are blended in the image plane (see the best-fit source-plane images in Fig. D.3), so we present the combined line ratios of each galaxy in this close pair. Since we are examining ratios of nearby emission line fluxes, we did not apply corrections for dust reddening or gravitational lensing.

Fig. 9. [SII]-VO87 (upper) and [NII]-BPT (lower) plots created using best-fit line fluxes for each source (see the values in Table 3). The distributions of low-redshift galaxies from SDSS (MPA-JHU DR8 catalogue; Kauffmann et al. 2003b; Brinchmann et al. 2004) are shown as background points. Left-facing arrows represent 3σ upper limits on N2 and S2. We include the demarcation lines for the z ∼ 0 galaxies of Kewley et al. (2001, solid lines) and Kauffmann et al. (2003a, dashed line). We compare our results with values from z ∼ 5.5 − 7.0 galaxies observed with JWST/NIRSpec MSA as part of the JADES survey (Cameron et al. 2023; black symbols).

For the [NII]-BPT, sources that lie above the solid Kewley et al. (2001) line are believed to be dominated by AGN excitation, while those that lie beneath the dashed Kauffmann et al. (2003a) line are mainly star-forming, and those that lie between the lines are a combination of the two excitation sources (i.e. “composite”). Similarly, sources below (above) the solid Kauffmann et al. (2003a) line in the [SII]-VO87 are thought to be star-forming (AGN-dominated).

We find that while the sources in S and W have high R3 ratios (comparable to the R3 ratios of z ∼ 5.5 − 7.0 galaxies observed with JWST/NIRSpec MSA as part of the JADES survey; Cameron et al. 2023), the high upper limits on N2 and S2 do not allow us to rule out or confirm the presence of AGN for most components. However, the detection of [NII]λ6584 in S1 shows that this point lies on the demarcation line for AGN activity. Component C lies in the star-forming regime, but with large errors that may place it in the composite region.

We note that these demarcation lines were derived for z ∼ 0 galaxies with approximately solar metallicity. Recent results suggest that these lines may not separate SF- and AGN-driven ionisation for high-redshift, low-metallicity sources (e.g. Kewley et al. 2013; Nakajima & Maiolino 2022; Übler et al. 2023). This interpretation is further complicated by the fact that the component C is composed of two sub-components (see Sects. 4.1 and 4.2), which may exhibit different ionisation sources and/or metallicities. These line ratios may also be affected by shocks within galaxies (e.g. Allen et al. 2008), which we explore further in Appendix E.

These results suggest that the line ratios of components S and W do not allow us to robustly claim the presence or absence of AGN, while component C is likely composed of star-forming galaxies but may include an AGN.

4.4. Distribution of star formation

HFLS3 was originally found to be an extreme starburst, with an SFR_FIR = 2900 M_⊙ yr⁻¹ (Riechers et al. 2013). Further studies found that the central galaxy was lensed, resulting in a lower intrinsic SFR ∼ 1300 M_⊙ yr⁻¹ or a 100-Myr averaged SFR ∼ 660 M_⊙ yr⁻¹ (Cooray et al. 2014). However, both of these SFRs are based on fits to SEDs that contain low-resolution observations (e.g. Herschel/SPIRE with > 20″, Spitzer/IRAC with ∼2″), which would capture emission from multiple components of the HFLS3 system (i.e. the lensed emission in C, both galaxies in S, and both galaxies in W). This raises the question of how the star formation activity is distributed in the HFLS3 system: a single dominant starburst or multiple star-forming galaxies within an area of projected diameter ∼11 kpc. This may be explored using the Hα luminosity (which mainly traces star formation on timescales of ∼10 − 20 Myr; e.g. Kennicutt 1998; Glazebrook et al. 1999), corrected for dust attenuation.

To investigate this, we first derived an estimate of the B − V colour excess from the Balmer decrement (e.g. Domínguez et al. 2013):

$$\begin{array}{c}\hfill E{(B-V)}_{\mathrm{BD}}=\frac{2.5}{k({\lambda }_{\mathrm{H}\beta })-k({\lambda }_{\mathrm{H}\alpha })}{log}_{10}\left(\frac{F_{\mathrm{H}\alpha ,\mathrm{obs}}/F_{\mathrm{H}\beta ,\mathrm{obs}}}{2.86}\right),\end{array}$$(5)

where k(λ) is the assumed dust attenuation curve (Calzetti et al. 2000), F_Hα, obs/F_Hβ, obs is the observed line flux ratio, and we have assumed an intrinsic line ratio of 2.86 (Osterbrock 1989). We then used this value to estimate the intrinsic line flux of Hα, assuming case B recombination and T_e = 10⁴ K:

$$\begin{array}{c}\hfill F_{\mathrm{H}\alpha ,\mathrm{int}}=F_{\mathrm{H}\alpha ,\mathrm{obs}}{10}^{k({\lambda }_{\mathrm{H}\alpha })E{(B-V)}_{\mathrm{BD}}/2.5}.\end{array}$$(6)

Finally, this dust-corrected flux is used to estimate the SFR assuming a Salpeter IMF (Kennicutt 1998):

$$\begin{array}{c}\hfill {\mathrm{SFR}}_{\mathrm{H}\alpha }=\frac{4\pi D_{\mathrm{L}}^2{F}_{\mathrm{H}\alpha ,\mathrm{int}}}{1.26\times {10}^{41}\mathrm{erg}{\mathrm{s}}^{-1}},\end{array}$$(7)

where D_L is the luminosity distance. The resulting values are presented in Table 5. The Calzetti et al. (2000) dust attenuation curve was chosen to ease comparison to previous results, but other curves could be adopted. For example, the quadratic law of Reddy et al. (2020) results in slightly higher (≲1σ difference) values of SFR, while the high-mass, high-redshift analogue law of Salim et al. (2018) returns slightly lower (≲2σ difference) values of SFR.

From this analysis, we see that the central component indeed features a high SFR. But as found in Sect. 4.2, this emission is lensed, resulting in a dust- and magnification-corrected SFR of 160 ± 80 M_⊙ yr⁻¹ (assuming an average μ ∼ 2.26 ± 0.17). The combined SFR of the z ∼ 6.3 galaxies is 510 ± 140 M_⊙ yr⁻¹, which is comparable to the previous 100 Myr-averaged SFR of Cooray et al. (2014), or ${654}_{-90}^{104}$ M_⊙ yr⁻¹. But if S1, S2, W1, and/or W2 contain AGN (which is not ruled out by their line ratios; see Sect. 4.3), then their SFR_Hα values may be contaminated (e.g. Garn et al. 2010). On the other hand, if these sources contain star formation in optically thick regions, then the true SFR may be higher.

In summary, the star formation (as traced by Hα) is not concentrated in a single source but is distributed in multiple objects: a central lensed arc that has two components (∼32% of the SFR), two galaxies to the south (∼17% of the SFR), and two galaxies to the west (∼51% of the SFR). While the absolute SFR of this system is dependent on the presence of dust and AGN (which require further high-resolution radio and sub-millimetre observations to be confirmed), the current data imply that HFLS3 is not a single starbursting galaxy, as previously reported.

4.5. HFLS3: A galaxy group

In the previous sections, we have found that the HFLS3 field contains at least six strongly detected sources at z ≈ 6.35. Here, we discuss the implications of this apparent high density of galaxies.

Previous large-scale optical searches around HFLS3 found little evidence for an overdensity on the scale of ≳100 kpc (Robson et al. 2014; Laporte et al. 2015). However, Laporte et al. (2015) detected three faint optical companions to HFLS3 on small scales (∼36 kpc), based on HST imaging with comparable spatial resolution to the current work. Two of these sources are detected in our data (S and W), while the third source falls outside the IFU FoV. Since we have found that each of these sources and the central source are composed of two galaxies, the true density is likely higher. Whether this represents a true galaxy overdensity requires knowledge of the Hα emitter space density at z > 6, which is not yet well predicted (e.g. Pozzetti et al. 2016). Alternatively, we could compare the observed galaxy distribution to a cosmological simulation that includes Hα emission and dust extinction (e.g. Hashimoto et al. 2023). This will be deferred to a future work.

At this time, we could use the spatial and velocity offsets of each galaxy to see if they are likely gravitationally interacting. Since these values are projected, they represent lower limits on the true three-dimensional offset and relative velocity of galaxy pairs, but may be used as a first test of close association. Studies of galaxy pairs usually adopt criteria of Δr ≲ 20 kpc and Δv ≲ 500 km s⁻¹ for “close” pairs (e.g. Duncan et al. 2019; Ventou et al. 2019; Romano et al. 2021). The six z > 6 galaxies detected here are all located within a circle of diameter ∼2″ (∼11 kpc at z = 6.34) with a maximum redshift difference of δz ∼ 0.02 (i.e. a line-of-sight velocity difference of ∼800 km s⁻¹). So while each galaxy in this field firmly meets the distance criterion, the total velocity difference is high.

However, each galaxy is a member of a pair with smaller spatial offsets and redshift offsets: C1 and C2 (δz ∼ 0.0092, using the redshifts of Riechers et al. 2013), S1 and S2 (δz ∼ 0.0001), and W1 and W2 (δz ∼ 0.0078). Since the distance and velocity criteria are met for each pair, they are likely interacting. In addition, the galaxy pairs in the S and W regions are also likely interacting, while C features a slightly higher velocity offset (Δv = 591 ± 11 km s⁻¹ with respect to z_W1). Based on these small velocity and spatial offsets, it is likely that these galaxies represent an interacting system that will merge within ∼1 Gyr (e.g. Conselice 2006).

5. Conclusion

In this work we present JWST/NIRSpec IFU observations of a field containing the z = 6.34 source HFLS3, as part of the GA-NIFS programme. By exploring both the low (R ∼ 100) and high spectral resolution (R ∼ 2700) data, we find a crowded field, with two low-redshift sources (G1, z ∼ 3.4806; G2, z ∼ 2.00), a central gravitationally lensed arc that is composed of two sources (C, z ∼ 6.3425), two close galaxies at the same redshift to the south (S1 and S2, z ∼ 6.3592), and two close galaxies with a velocity offset to the west (W1, z ∼ 6.3550; W2, z ∼ 6.3628). All of the z > 6 galaxies are located within an area of ∼2″, or ∼11 kpc at z ∼ 6.34.

The spectral fits and morpho-kinematic map analysis of these data reveal a variety of kinematic features. G1 has no strong velocity gradient and is likely a dispersion-dominated source. However, C features a strong velocity gradient across the length of its arc, hinting at possible merging activity. The galaxies in component S are distinct, but are at similar redshifts with no velocity gradient. This is different from the galaxies in component W, which are closely associated but have a strong velocity gradient, as seen in the double-peak profile of each line. The red component (W2) is brighter but has a lower velocity dispersion. Because of the asymmetry in the integrated intensity and velocity dispersion maps, this likely represents a merger between two galaxies, rather than a single rotating disk.

Next, we used our high-quality IFU data and updated source redshifts to examine the gravitational lensing of the HFLS3 field. Our best-fit models show that G1 is moderately magnified by G2 (μ ∼ 1.1 − 1.2), while the Hα emission of C is more strongly magnified (μ ∼ 2.1 − 2.4, comparable to the value of ∼2.2 from Cooray et al. 2014) and is composed of two closely separated components.

The integrated line fluxes are plotted on [NII]-BPT and [SII]-VO87 plots, showing that component C is likely powered by star formation. We cannot rule out or confirm the presence of AGN in the S and W components.

The observed Balmer decrements were then used to derive extinction-corrected SFR_Hα values for each source. This shows that the star formation is distributed across the field (SFR_Hα = 510 ± 140 M_⊙ yr⁻¹, corrected for lensing and extinction), with the largest contribution (∼51%) from the W galaxies. However, the presence of AGN in some sources may inflate these SFRs.

We used the projected spatial offsets and relative line-of-sight velocities of each z > 6 galaxy to investigate whether they likely represent a closed, merging group. Using the standard “close pair” criteria of δr < 20 kpc and δv < 500 km s⁻¹, we find that the field contains multiple likely mergers: C1/C2, S1/S2, and W1/W2. In addition, the combined W and S components meet these criteria and will likely merge. While the C components feature a slightly higher velocity offset (∼600 km s⁻¹), the HFLS3 group contains multiple close pairs.

Taken together, our results require a drastic reinterpretation of the HFLS3 field. It is composed of at least six distinct sources within z ∼ 6.34 − 6.36 that are lensed by two foreground galaxies, at z ∼ 2.0 and z ∼ 3.5. All of the z > 6 galaxies feature strong velocity gradients and/or are closely associated with another galaxy, implying ongoing interaction. This behaviour has been seen at high redshift (e.g. Ginolfi et al. 2022), including four galaxies at z ∼ 7.9 detected with the NIRSpec IFU (Hashimoto et al. 2023). Thus, HFLS3 is likely not an extreme starburst, but instead represents one of the densest groups of interacting star-forming galaxies within the first 1 Gyr of the Universe. Recent and ongoing high-resolution observations with JWST/MIRI, NOEMA, and the JVLA (as well as a future in-depth study of the R100 data cube briefly explored here) will help further characterise this unique field.

Acknowledgments

G.C.J., A.J.B., A.J.C., J.C., and A.S. 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. M.P., S.A., and B.Rd.P. 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 Programa Atracción de Talento de la Comunidad de Madrid via grant 2018-T2/TIC-11715. S.C. and G.V. acknowledge support from the European Union (ERC, WINGS,101040227). R.M., F.D.E., J.S., and J.W. acknowledge support by the Science and Technology Facilities Council (STFC), from the ERC Advanced Grant 695671 “QUENCH”, and by the UKRI Frontier Research grant RISEandFALL. R.M. also acknowledges funding from a research professorship from the Royal Society. R.B. acknowledges support from an STFC Ernest Rutherford Fellowship [grant number ST/T003596/1]. G.C. 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”. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

Appendix A: Astrometry correction

A.1. HST astrometry verification

In order to compare the JWST/NIRSpec IFU data with archival data, a common astrometric reference frame is required (i.e. Gaia DR3; Gaia Collaboration 2016, 2021). We retrieved HST images from the MAST archive, but it is not clear if all of them have been aligned to the Gaia DR3 frame (i.e. some are lacking this comment in their headers). To verify that they have been correctly aligned, we retrieved the locations of the five closest objects in the Gaia archive⁶ (see Table A.1) and examined the locations in each HST image.

As seen in Fig. A.1, the HST images show significant emission at all five Gaia locations with small centroid offsets (< 0.2″). From this, we conclude that the images have been properly aligned.

Fig. A.1. Verification that HST images have been aligned to the Gaia DR3 reference frame. Each panel contains a 1″ × 1″ view of an HST image (see the labels for each row) that is focused on a given location from the Gaia archive (see the labels for each column). The Gaia position is shown as a 0.1″ diameter red circle. The locations of each Gaia source are listed in Table A.1.

A.2. Application of spatial shift

To shift the JWST/NIRSpec data cubes to the Gaia frame, we first created narrowband JWST images by convolving the R100 data cube with the transmission profiles of the three HST/WFC3 filters used by Cooray et al. (2014) (see the colour maps in Fig. A.2)⁷. These maps were then compared with the corresponding HST images (see the contours in Fig. A.2). An offset is estimated by finding the brightest spaxel in the JWST and HST maps for each filter and determining the distance between the centres of these spaxel.

Fig. A.2. JWST/NIRSpec IFU R100 data cubes integrated over the corresponding HST filter bandpasses: F105W (left), F125W (centre), and F160W (right). In each panel, the collapsed JWST emission is depicted as the background colours, while the HST data are shown as red contours. The JWST emission in the top row is shown without the astrometric correction, while the lower row includes the alignment to Gaia DR3 data.

In all three cases, this offset is 0.17″. Because the spaxel sizes are 0.05″ (JWST) and 0.12825″ (HST/WFC3-IR), we estimated an uncertainty on this offset based on half the spaxel sizes, added in quadrature (±0.07″). Thus, the offset is consistent with the pointing accuracy of ∼0.1″ (Rigby et al. 2023). We verified that the R100 and R2700 cubes are aligned to the same initial reference frame, so throughout this work we applied this offset correction to the R100 and R2700 IFU data cubes.

We briefly note that Cooray et al. (2014) aligned each image to the Sloan Digital Sky Survey (SDSS) frame rather than the Gaia frame. Comparing objects near the HFLS3 field in SDSS DR9 and Gaia DR3, we find a negligible offset between these frames. While the previous HST data featured a small astrometric uncertainty (∼0.05″), the Keck/NIRC2 and PdBI data have larger uncertainties (∼0.1″ and 0.1″ − 0.3″, respectively). When combined with the ∼0.1″ pointing uncertainty of our NIRSpec/IFU data, we find no disagreement between previous positions and our work.

Appendix B: Additional candidate galaxies

The HFLS3 field contains multiple galaxies that are strongly detected in line and/or continuum emission. In addition, there are several areas of weak emission that we do not include in our analysis. We present additional details of each below.

B.1. North-east extension (X1)

The first candidate emission area is visible as an elongated stretch of emission in the upper panel of Fig. B.1 with a slightly brighter clump at the western edge. This morphology (a bar aligned with the axis of the IFU slices used to create the image cube) is similar to artefacts introduced by bright sources that saturate a detector slice (e.g. Böker et al. 2022), but our image is unsaturated. Instead, the high-contrast image of Fig. B.1 shows that multiple slices show a speckled noise pattern, possibly reflecting a low-level calibration issue in this wavelength range.

Fig. B.1. Integrated emission of the HFLS3 field, using the same illustrative wavelength ranges as in Fig. 1: redshifted Hα for z ∼ 6.34 for the R2700 cube (λobs = 4.79954  −  4.84467 μm, upper panel) and the approximate wavelength range of HST/WFC3 F160W for the R100 cube (λobs = 1.4  −  1.6 μm; lower panel). The colour scale has been adjusted to highlight weak emission. Low-level candidate galaxies are shown with cyan boundaries. North is up and east is to the left.

The western clump in X1 is brighter than this noise pattern, but is substantially weaker than the primary components (e.g. C, W). An integrated spectrum over this rectangular region returns strong detections in Hβ, [OIII]λλ4959, 5007, and Hα, yielding a redshift of z_X1 = 6.3438 ± 0.0001 (see Table 3 and Fig. B.2). The unique redshift and line properties suggest that this represents a true galaxy rather than a noise artefact or a lensed image of another galaxy. However, the weak detection of Hβ and the non-detection of [NII]λλ6548, 6584 and [SII]λλ6716, 6731 limits our ability to characterise the source. Its detection suggests that the HFLS3 field contains other low-level galaxies that may contribute to the star formation activity but are undetected by our analysis, as suggested by cosmological simulations (e.g. Pallottini et al. 2019).

Fig. B.2. Integrated spectra of the R2700 cube using the masks of Fig. B.1, with 1σ errors from the associated error spectrum shown as the shaded region. Best-fit models (line emission and continuum) are shown as green lines. The centroids of each line are depicted by dashed lines, with red lines indicating that the spectral line was not fit. The lower panel shows the residual.

B.2. South-west clumps (X2)

Figure B.1 shows two low-level clumps of emission to the south-west of component S in the collapsed R2700 image. These points lie along the same IFU slice, and a spectrum extracted over a region containing both reveals moderate emission that may correspond to Hα emission at z ∼ 6.36. However, the amplitude of this emission is quite low (< 2σ) with the exception of two narrow peaks that do not align with the centres of lines. Because of the low amplitude and artefact-like morphology, we did not investigate this region further.

B.3. G2 extension (X3)

When the R100 cube is collapsed over a wavelength range that is dominated by continuum emission from the z < 6 components (λ_obs = 0.8 − 1.1 μm; see the upper-left panel of Fig. 8), there is a circular source of emission to the south-east of G2 (hereafter X3). This emission is not captured by the Sérsic mode that we fitted to the data (see the residual in the upper-right panel of Fig. 8). There is a significant emission line in this region in the R100 data at λ_obs = 4.822 μm, but a moment zero image shows that this is simply low-level Hα from the outskirts of component C. Our lens modelling analysis does not show an image-plane component to the south, so X3 is likely not a lensed image of C. This extension of G2 also appears in archival HST images (red contours of Fig. A.2), suggesting that it is not an artefact of our JWST data. Further high spatial resolution observations are required for a spectral redshift of this source, in order to determine its nature.

Appendix C: Pseudo-Voigt profile fitting

C.1. Profile details

A standard Voigt profile is a convolution of a Gaussian and Lorentz profile. To simplify the computation of this profile, we adopted the pseudo-Voigt profile as implemented in LMFIT. This begins with a normalised Gaussian:

$$\begin{array}{c}\hfill G(\lambda )=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-{(\lambda -{\lambda }_C)}^2/(2{\sigma }^2)},\end{array}$$(C.1)

which is centred on λ_C and features a FWHM=$2\sqrt{2ln(2)}\sigma$. Next, we considered a normalised Lorentzian with the same centroid and FWHM:

$$\begin{array}{c}\hfill L(\lambda )=\frac{1}{\pi }\frac{\sqrt{2ln(2)}\sigma }{{(\lambda -{\lambda }_C)}^2+2ln(2){\sigma }^2}.\end{array}$$(C.2)

A pseudo-Voigt profile with integrated area A is then given by

$$\begin{array}{c}\hfill V(\lambda )=(1-\alpha )AG(\lambda )+\alpha AL(\lambda ),\end{array}$$(C.3)

where α is the fraction of emission from the Lorentzian component. We note that this form is flexible, with the ability to encapsulate Gaussian profiles (α = 0), Lorentzian profiles (α = 1), and all combinations of these profiles such that the integrated intensity, centroid, and FWHM are conserved (0 < α < 1).

C.2. Morpho-kinematic map creation

In order to examine the morpho-kinematics of each source and line, we can produce maps of integrated intensity, velocity offset, and velocity dispersion. These are usually derived by extracting spectra from each spaxel (f(v)) and measuring a normalised cumulative velocity distribution ($\mathrm{F}(\mathit{v})={\int }_{-\infty }^{\mathit{v}}\mathrm{f}(\mathit{v}\prime )\mathrm{d}\mathit{v}\prime /{\int }_{-\infty }^{\infty }\mathrm{f}(\mathit{v}\prime )\mathrm{d}\mathit{v}\prime$; e.g. Zakamska & Greene 2014), which is used to calculate the velocity at which N% of the flux is captured (v_N; e.g. v₅₀) as well as associated widths (e.g. w₈₀ ≡ v₉₀ − v₁₀). Since there are overlapping lines in some spectra, we adopted a pseudo-Voigt-based approach.

For each spaxel, we extracted a spectrum and fit a model containing a flat continuum and one or more pseudo-Voigt lines using LMFIT with a least-squares minimiser. Lines that are closely associated or related through a flux ratio are fit concurrently (i.e. [OIII]λ4959/[OIII]λ5007, [NII]λ6548/Hα/[NII]λ6584, and [SII]λ6716/[SII]λ6731) with identical continuum values. Line pairs (i.e. [OIII]λ4959/[OIII]λ5007, [NII]λ6548/[NII]λ6584, and [SII]λ6716/[SII]λ6731) are assumed to have identical kinematics. We adopted the standard assumptions of [NII]λ6584/[NII]λ6548 = 2.94 (e.g. Dojčinović et al. 2023) and [OIII]λ5007/[OIII]λ4959 = 2.98 (e.g. Dimitrijević et al. 2007) for each fit. The LSF is accounted for when calculating the linewidths. Continuum maps are generated using the best-fit constant continuum value for each fit and spaxel. Only fits with r² (i.e. the coefficient of determination) values > 0.5 are presented.

Using the best-fit parameters of each model pseudo-Voigt profile (Integrated intensity A [erg s cm⁻²], centroid wavelength λ_C [μm], FWHM [μm], and Lorentzian fraction α), we could generate morpho-kinematic maps. Since the Gaussian and Lorentzian components of each pseudo-Voigt profile are normalised, the integrated intensity map is simply the best-fit A:

$$\begin{array}{c}\hfill \frac{I(x,y)}{{\mathrm{erg}\mathrm{s}\mathrm{cm}}^{-2}}=A(x,y).\end{array}$$(C.4)

The velocity field (v₅₀) represents the dominant line-of-sight velocity for a given spaxel. Since the Gaussian and Lorentzian components of the pseudo-Voigt profile are symmetric around the same centroid, v₅₀ may be expressed for a redshifted rest-frame wavelength of λ_o as

$$\begin{array}{c}\hfill \frac{v_{50}(x,y)}{{\mathrm{km}\mathrm{s}}^{-1}}=c(\frac{{\lambda }_C(x,y)}{{\lambda }_o}-1),\end{array}$$(C.5)

where c[km s⁻¹] is the speed of light.

The velocity dispersion is more complex. While the pseudo-Voigt profile is constructed such that its FWHM is constant with respect to the Lorentzian fraction α, the non-parametric width w₈₀ is dependent on α. To calculate this value, we generated a cumulative velocity distribution of each best-fit model spectrum and extract v₁₀ and v₉₀. Using these,

$$\begin{array}{c}\hfill \frac{w_{80}(x,y)}{{\mathrm{km}\mathrm{s}}^{-1}}=v_{90}-v_{10}.\end{array}$$(C.6)

In this process, we assumed that each spectral line may be fit with a single pseudo-Voigt profile. While this is motivated by the lack of an obvious broad component (as seen in e.g. Marshall et al. 2023), it is possible that a spaxel may contain contributions from narrow and broad emission. Since there are no strong AGN (from previous observations) included in the HFLS3 field, this is expected to be a small effect.

Appendix D: Best-fit lens models of components

In Sect. 4.2 we use PYAUTOLENS to examine the gravitational lensing effect created by the z < 6 sources (G1 and G2) on the C source. The observed maps were compared with the best-fit image-plane model. Here, we show source- and image-plane maps of each individual component: G1 (Fig. D.1), G2 (Fig. D.2), and C (Fig. D.3).

Fig. D.1. Best-fit source-plane (left column) and image-plane PYAUTOLENS models (right column) of G1 using the R100 cube (λ = 0.8 − 1.1 μm; top row) and the R2700 cube (λ = 4.80 − 4.85) assuming that C is composed of a single component (middle row) or two components (lower row). Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.

Fig. D.2. Best-fit PYAUTOLENS model of the unlensed source G2 using the R100 cube (λobs = 0.8 − 1.1 μm). Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.

Fig. D.3. Best-fit source-plane (top row) and image-plane PYAUTOLENS models (bottom row) of component C using the R2700 cube (λobs = 4.80 − 4.85) assuming that it is composed of a single component (left column) and the individual models of the two-component model (centre and right columns). Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.

Appendix E: Effects of shocks on line ratio diagnostics

In Sect. 4.3 we briefly discuss the applicability of low-redshift demarcation lines for each high-redshift source. Here, we explore the effects of shocks, as done for the local galaxy group Stephan’s Quintet (Duarte Puertas et al. 2021). This previous work applied the MAPPINGS III (Allen et al. 2008) shock models, which contain expected line ratios for a variety of shock velocities, initial density, and magnetic field strength. The resulting model grids with our data are shown in Figure E.1, adopting the same assumptions of solar metallicity, low density (n ∼ 0.1 cm⁻³), and no shock precursor.

Fig. E.1. [SII]-VO87 (top) and [NII]-BPT (lower) plots created using best-fit line fluxes for each source (see the values in Table 3), as seen in Figure 9. We now include the model grids from MAPPINGS III (Allen et al. 2008). Red-shaded lines display models with a constant magnetic parameter ($\mathrm{MP}=B/\sqrt{n}$ in units of μG cm3/2), while blue-shaded lines represent models with a constant shock velocity. We assume solar metallicity, low density (n ∼ 0.1 cm−3), and no shock precursor.

It is clear that the high-R3 sources (S, W) are not captured by these models, while component C falls within the low-velocity (v ∼ 200 km s⁻¹), low-magnetic field (MP ∼ 0.001 μG cm^3/2) portion of the grid. Many emission line systems in Stephan’s Quintet also fall into this region (Duarte Puertas et al. 2021), which is interpreted as evidence for the presence of shocks. While this may also be true for component C of HFLS3, we note that these models were derived for relatively high metallicities (Z′ = 1), and lower metallicity models may be required for the high-redshift galaxies.

All Tables

All Figures

Fig. 1. Integrated emission of the HFLS3 field, using two illustrative wavelength ranges: redshifted Hα for z ∼ 6.34 for the R2700 cube (λobs = 4.79954 − 4.84467 μm, upper panel) and the approximate wavelength range of HST/WFC3 F160W for the R100 cube (λobs = 1.4 − 1.6 μm; lower panel). The adopted source masks are shown with white boundaries. North is up and east is to the left. The lower panel is nearly identical to the observed HST/WFC3 F160W image, but at a higher spatial resolution (see Fig. A.2).

In the text

Fig. 2. Integrated spectra of the R2700 cube using the masks of Fig. 1 (magenta; excluding G2), with 1σ errors from the associated error spectrum shown as the shaded region. Best-fit models (line emission and continuum) are shown by green lines. The centroids of each line are depicted by dashed lines, with red lines indicating that the spectral line was not fit. The lower panel shows the residual.

In the text

	Fig. 2. continued.
In the text

	Fig. 2. continued.
In the text

Fig. 3. Integrated spectrum extracted from the R100 cube using the G2 mask of Fig. 1, with illustrative 1σ errors shown as shaded region. A basic fit to the data is shown by the green line. The lower panel shows the residual, with the centroids of each line depicted by dashed lines. From left to right, the lines are Hβ, [OIII]λλ4959, 5007, [NII]λ6548, Hα, and [NII]λ6584. The best-fit redshift is zG2 = 2.00 ± 0.01.

In the text

Fig. 4. Morpho-kinematic maps of the strong emission lines detected in the component G1 (z = 3.4806). Integrated intensities are given per spaxel. The corresponding maps of [SII]λλ6716, 6731 are not presented due to low S/N. Due to the construction of our model, [NII]λ6548 has the same integrated velocity field, velocity dispersion map, and continuum map as [NII]λ6584, but an integrated intensity map that is a factor of 2.98 lower. For [SIII]λ9532 and HeIλ10829, v50 values are given with respect to their best-fit redshifts (z = 3.4810 and 3.4821, respectively).

In the text

Fig. 5. Morpho-kinematic maps of the strong emission lines detected in the C component (z = 6.3425). Integrated intensities are given per spaxel. The corresponding maps of Hβ, [OIII]λλ4959, 5007, and [SII]λλ6716, 6731 are not presented due to low S/N. Due to the construction of our model, [NII]λ6548 has the same integrated velocity field, velocity dispersion map, and continuum map as [NII]λ6584, but an integrated intensity map that is a factor of 2.98 lower.

In the text

Fig. 6. Morpho-kinematic maps of the strong emission lines detected in the S components (z ∼ 6.3592). Integrated intensities are given per spaxel. The corresponding maps of Hγ are not presented due to low S/N. Due to the construction of our model, [OIII]λ4959 has the same velocity field, velocity dispersion map, and continuum map as [OIII]λ5007, but an integrated intensity map that is a factor of 2.94 lower.

In the text

Fig. 7. Morpho-kinematic maps of the strong emission lines detected in the W components (zW1 = 6.3550, zW2 = 6.3628). Integrated intensities are given per spaxel. The corresponding maps of Hγ are not presented due to low S/N. Due to the construction of our model, [OIII]λ4959 has the same velocity field, velocity dispersion map, and continuum map as [OIII]λ5007, but an integrated intensity map that is a factor of 2.94 lower.

In the text

Fig. 8. Results of PYAUTOLENS lensing analysis of collapsed images of the R100 cube (λobs = 0.8 − 1.1 μm; top row) and the R2700 cube (λobs = 4.80 − 4.85 μm) assuming that component C is composed of a single source (middle row) or two sources (C1 to the east and C2 to the west; lower row). The observed data, best-fit model, and residual (all in the image plane) are presented from left to right. Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.

In the text

Fig. 9. [SII]-VO87 (upper) and [NII]-BPT (lower) plots created using best-fit line fluxes for each source (see the values in Table 3). The distributions of low-redshift galaxies from SDSS (MPA-JHU DR8 catalogue; Kauffmann et al. 2003b; Brinchmann et al. 2004) are shown as background points. Left-facing arrows represent 3σ upper limits on N2 and S2. We include the demarcation lines for the z ∼ 0 galaxies of Kewley et al. (2001, solid lines) and Kauffmann et al. (2003a, dashed line). We compare our results with values from z ∼ 5.5 − 7.0 galaxies observed with JWST/NIRSpec MSA as part of the JADES survey (Cameron et al. 2023; black symbols).

In the text

	Fig. A.1. Verification that HST images have been aligned to the Gaia DR3 reference frame. Each panel contains a 1″ × 1″ view of an HST image (see the labels for each row) that is focused on a given location from the Gaia archive (see the labels for each column). The Gaia position is shown as a 0.1″ diameter red circle. The locations of each Gaia source are listed in Table A.1.
In the text

Fig. A.2. JWST/NIRSpec IFU R100 data cubes integrated over the corresponding HST filter bandpasses: F105W (left), F125W (centre), and F160W (right). In each panel, the collapsed JWST emission is depicted as the background colours, while the HST data are shown as red contours. The JWST emission in the top row is shown without the astrometric correction, while the lower row includes the alignment to Gaia DR3 data.

In the text

Fig. B.1. Integrated emission of the HFLS3 field, using the same illustrative wavelength ranges as in Fig. 1: redshifted Hα for z ∼ 6.34 for the R2700 cube (λobs = 4.79954  −  4.84467 μm, upper panel) and the approximate wavelength range of HST/WFC3 F160W for the R100 cube (λobs = 1.4  −  1.6 μm; lower panel). The colour scale has been adjusted to highlight weak emission. Low-level candidate galaxies are shown with cyan boundaries. North is up and east is to the left.

In the text

	Fig. B.2. Integrated spectra of the R2700 cube using the masks of Fig. B.1, with 1σ errors from the associated error spectrum shown as the shaded region. Best-fit models (line emission and continuum) are shown as green lines. The centroids of each line are depicted by dashed lines, with red lines indicating that the spectral line was not fit. The lower panel shows the residual.
In the text

Fig. D.1. Best-fit source-plane (left column) and image-plane PYAUTOLENS models (right column) of G1 using the R100 cube (λ = 0.8 − 1.1 μm; top row) and the R2700 cube (λ = 4.80 − 4.85) assuming that C is composed of a single component (middle row) or two components (lower row). Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.

In the text

	Fig. D.2. Best-fit PYAUTOLENS model of the unlensed source G2 using the R100 cube (λobs = 0.8 − 1.1 μm). Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.
In the text

Fig. D.3. Best-fit source-plane (top row) and image-plane PYAUTOLENS models (bottom row) of component C using the R2700 cube (λobs = 4.80 − 4.85) assuming that it is composed of a single component (left column) and the individual models of the two-component model (centre and right columns). Mass and light profiles are assumed to be Sérsic and isothermal ellipse profiles, respectively. The outline of the spatial mask is shown by red markers.

In the text

Fig. E.1. [SII]-VO87 (top) and [NII]-BPT (lower) plots created using best-fit line fluxes for each source (see the values in Table 3), as seen in Figure 9. We now include the model grids from MAPPINGS III (Allen et al. 2008). Red-shaded lines display models with a constant magnetic parameter ($\mathrm{MP}=B/\sqrt{n}$ in units of μG cm3/2), while blue-shaded lines represent models with a constant shock velocity. We assume solar metallicity, low density (n ∼ 0.1 cm−3), and no shock precursor.

In the text