11institutetext:  Department of Physics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, UK 22institutetext: Kavli Institute for
Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
33institutetext: Cavendish Laboratory, University of Cambridge, 19 JJ Thomson
Avenue, Cambridge CB3 0HE, UK 44institutetext: Centro de Astrobiología (CAB),
CSIC–INTA, Cra. de Ajalvir Km. 4, 28850- Torrejón de Ardoz, Madrid, Spain
55institutetext: Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126
Pisa, Italy 66institutetext: Sorbonne Université, CNRS, UMR 7095, Institut
d’Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France
77institutetext: Department of Physics and Astronomy, University College
London, Gower Street, London WC1E 6BT, UK 88institutetext: NRC Herzberg, 5071
West Saanich Rd, Victoria, BC V9E 2E7, Canada 99institutetext: Jodrell Bank
Centre for Astrophysics, Department of Physics and Astronomy, School of
Natural Sciences, The University of Manchester, Manchester, M13 9PL, UK
1010institutetext: European Space Agency, c/o STScI, 3700 San Martin Drive,
Baltimore, MD 21218, USA 1111institutetext: INAF - Osservatorio Astrofisco di
Arcetri, largo E. Fermi 5, 50127 Firenze, Italy 1212institutetext: European
Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany
1313institutetext: AURA for European Space Agency, Space Telescope Science
Institute, 3700 San Martin Drive. Baltimore, MD, 21210, USA

# GA-NIFS: JWST/NIRSpec IFU observations of HFLS3 reveal a dense galaxy group
at $z\sim 6.3$

Gareth C. Jones E-mail<EMAIL_ADDRESS>Hannah Übler 2233
Michele Perna 44 Santiago Arribas 44 Andrew J. Bunker 11 Stefano Carniani 55
Stephane Charlot 66 Roberto Maiolino 223377 Bruno Rodríguez Del Pino 44 Chris
Willott 88 Rebecca A. A. Bowler 99 Torsten Böker 1010 Alex J. Cameron 11
Jacopo Chevallard 11 Giovanni Cresci 1111 Mirko Curti 1212 Francesco D’Eugenio
2233 Nimisha Kumari 1313 Aayush Saxena 1177 Jan Scholtz 2233 Giacomo Venturi
55 Joris Witstok 2233

(Received X / Accepted Y)

Massive, starbursting galaxies in the early Universe represent some of the
most extreme objects in the study of galaxy evolution. One such source is
HFLS3 ($z\sim 6.34$), which was originally identified as an extreme starburst
galaxy with mild gravitational magnification ($\mu\sim 2.2$). Here, we present
new observations of HFLS3 with the JWST/NIRSpec IFU in both low (PRISM/CLEAR;
$R\sim 100$) and high spectral resolution (G395H/290LP; $R\sim 2700$), with
high spatial resolution ($\sim 0.1^{\prime\prime}$) and sensitivity. Thanks to
the combination of the NIRSpec data and a new lensing model with accurate
spectroscopic redshifts, we find that the $3^{\prime\prime}\times
3^{\prime\prime}$ field is crowded, with a lensed arc (C, $z=6.3425\pm
0.0002$), two galaxies to the south (S1 and S2, $z=6.3592\pm 0.0001$), two
galaxies to the west (W1, $z=6.3550\pm 0.0001$; W2, $z=6.3628\pm 0.0001$), and
two low-redshift interlopers (G1, $z=3.4806\pm 0.0001$; G2, $z=2.00\pm 0.01$).
We present spectral fits and morpho-kinematic maps for each bright emission
line (e.g., [OIII]$\lambda 5007$, H$\alpha$, [NII]$\lambda 6584$) from the
R2700 data for all sources except G2 (whose spectral lines fall outside the
observed wavelengths of the R2700 data). From a line ratio analysis, the
galaxies in component C are likely powered by star formation, while we cannot
rule out or confirm the presence of AGN in the other high-redshift sources. We
perform gravitational lens modelling, finding evidence for a two-source
composition of the lensed central object and a comparable magnification factor
($\mu=2.1-2.4$) to previous work. The projected distances and velocity offsets
of each galaxy suggest that they will merge within the next $\sim 1$ Gyr.
Finally, we examine the dust extinction-corrected SFRHα of each $z>6$ source,
finding that the total star formation ($510\pm 140$ M⊙ yr-1, magnification-
corrected) is distributed across the six $z\sim 6.34-6.36$ objects over a
region of diameter $\sim 11$ kpc. Altogether, this suggests that HFLS3 is not
a single starburst galaxy, but instead is a merging system of star-forming
galaxies in the Epoch of Reionisation.

###### Key Words.:

galaxies: high-redshift, galaxies: star formation, galaxies: kinematics and
dynamics, gravitational lensing: strong

##  1 Introduction

Observations have revealed that the mode of galaxy evolution in the first $2$
Gyr of the Universe ($z\gtrsim 3$) was drastically different from that of
following epochs. This can be seen in the similar evolution of the global star
formation rate density (SFRD) 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. 2019a), as well as the mean star formation rate
(SFR) and black hole accretion rate (BHAR) 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, resulting in the buildup of stellar and black hole
mass on shorter timescales than the following epochs.

The unique environment of the early Universe allowed galaxies to start forming
stars quickly, resulting in high-redshift galaxies that have exhausted their
gas supply and have already stopped forming stars by $z\gtrsim 3$ (e.g., see
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^{13}$ L⊙) and
SFRs of $\gtrsim 10^{3}$ M⊙ yr-1 (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${}^{\circ}30^{\prime}00.0^{\prime\prime}$), and was
included in the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et
al. 2012). Based on very red Herschel/SPIRE colours ($S_{500\mu m}>S_{350\mu
m}>S_{250\mu m}$ and $S_{500\mu m}/S_{350\mu 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\sim 2.092$ based on detection of
CIV$\lambda 1549$ with Keck/LRIS). A comprehensive suite of observations
(e.g.; Plateau de Bure Interferometer, PdBI; the Very Large Array, VLA; Keck;
William Herschel Telescope) enabled the creation of a NIR-radio spectral
energy distribution (SED) and the measurement of a number of emission lines
for HFLS3 ($z_{spec}=6.3369\pm 0.0009$). Based on their analysis of the line
and continuum emission, Riechers et al. (2013) report a high mass of molecular
gas ($\rm M_{H_{2}}=(1.04\pm 0.09)\times 10^{11}/\alpha_{CO}$ M⊙), dust ($\rm
M_{dust}=1.31^{+0.32}_{-0.30}\times 10^{9}$ M⊙), and stars ($\rm M_{*}\sim
3.7\times 10^{10}$ M⊙), as well as a high FIR luminosity ($\rm
L_{FIR}$=$2.86^{+0.32}_{-0.31}\times 10^{13}$ L⊙) and SFR${}_{\rm FIR}\sim
2900$ M⊙ yr-1 and evidence for a velocity gradient from PdBI [CII]158
$\mu$m111The values reported by Riechers et al. (2013) assume no gravitational
lensing.. 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^{\prime\prime}$) and thus
did not allow for detailed source differentiation and characterisation.

Follow-up Keck/NIRC2 and Hubble Space Telescope (HST) WFC3 imaging by Cooray
et al. (2014) revealed the presence of three close-by companions: two to the
north (G1 & G2) and one to the south (R1). The northern sources were assumed
to be at the same redshift ($z\sim 2.1$), while a photometric redshift of
$z_{phot}\sim 6$ was found for R1. Gravitational lens modelling confirmed that
HFLS3 may be composed of two sources, which are slightly magnified ($\mu\sim
2.2\pm 0.3$) by the foreground sources G1 & 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 JCMT/SCUBA-2 to perform a
low-resolution ($\sim 14^{\prime\prime}$), wide-area (67 arcmin2) search for
significant submm emission at two wavelengths ($450\,\mu$m and $850\,\mu$m),
but found no overdensity. This was followed by a search at optical wavelengths
with HST and GTC (Laporte et al. 2015), finding no large-scale overdensity but
a possible small-scale ($\sim 36$ kpc) overdensity of faint objects.

These observations suggest that HFLS3 is a gas-rich, starbursting, dusty
galaxy with a velocity gradient (suggesting merging activity) when the
Universe was only $\sim 850$ Myr old. Because of these exceptional properties,
it was chosen as a target for the JWST/NIRSpec 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\sim 3-11$
spanning a variety of types (e.g., QSOs, AGN, SFGs, 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$ AGN and QSOs are now 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\simeq
6.3-6.5$ and two low-redshift interlopers ($z\sim 2.0-3.5$) within a
$3^{\prime\prime}\times 3^{\prime\prime}$ field. We present the details of our
dataset in Section 2, and characterise the field in Section 3. Section 4
contains further analysis (i.e., morpho-kinematic maps, gravitational lens
modelling, line ratio-based excitation conditions, star formation rate
derivation, and galaxy merger discussion), and we conclude in Section 5.

We use a standard concordance cosmology ($h_{o}$, $\Omega_{m}$,
$\Omega_{\Lambda}$ = 0.7, 0.3, 0.7) throughout, where $1^{\prime\prime}$
corresponds to $\sim$5.5 kpc at $z\sim 6.34$, $\sim$7.3 kpc at $z\sim 3.48$,
and $\sim$8.3 kpc at $z\sim 2.1$. 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.,
$\lambda_{[OIII]\lambda 5007}=5008.24\textup{\AA}$). 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 PID 1264 (PI: L. Colina; see Table 1). Both
settings used the improved reference sampling and subtraction pattern IRS2 of
Rauscher et al. (2017) and adopted a 4-point dither with a medium ‘cycling’
pattern and a starting point of 1.

Table 1: JWST NIRSpec/IFU observation properties. Grating/Filter | Readout | Groups/Int | Ints/Exp | Exposures | Total Time
---|---|---|---|---|---
| Pattern |  |  |  | [s]
G395H/290LP | NRSIRS2 | 25 | 1 | 4 | 7352.801
PRISM/CLEAR | NRSIRS2RAPID | 60 | 1 | 4 | 3559.689

We use v1.8.2 of the GTO pipeline with CRDS context 1105 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; $\sim 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 have corrected count-rate
images for 1/f noise through a polynomial fit. During stage 2, we have removed
all data in regions of known failed open MSA shutters. We have 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 implement 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 FWHM is wavelength dependent
(FWHM${}_{\rm major}\sim 0.12-0.20^{\prime\prime}$ between $\lambda_{\rm
obs}=3.0-5.0$ $\mu$m). The combination of a dither and drizzle weighting
allowed us to sub-sample the detector pixels, resulting in cube spaxels of
$0.05^{\prime\prime}$. The resolving power of the R2700 data increases from
$\rm R_{2700}\sim 2000-3500$ from $\lambda_{\rm obs}=3.0-5.0$ $\mu$m, while
the resolving power of the R100 data decreases from $\rm R_{100}\sim 90-30$
over $\lambda_{\rm obs}=0.5-1.2$ $\mu$m and increases from $\rm R_{100}\sim
30-290$ over $\lambda_{\rm obs}=1.2-5.0$ $\mu$m.

###  2.1 Astrometric correction

A comparison of the JWST/NIRSpec IFU data with archival HST data222HST images
from Programme 13045 (PI: A. R. Cooray) were retrieved from the MAST archive.
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 DR3 frame,
we find an offset of $0.17\pm 0.07^{\prime\prime}$ (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 extract a mean spectrum from the R100 and R2700 data
cubes using a 25 spaxel ($12.5^{\prime\prime}$) diameter aperture (MPDAF task
aperture; Bacon et al. 2016) from a signal-free region to the southeast.

We assume that the background emission is uniform across the field of view
(FoV) and create a correction by subtracting this background spectrum from the
spectrum of each spaxel. Note that for the R2700 spectrum, we only include
channels that are not affected by chip gap/edge issues ($\lambda_{obs}\approx
2.87057-3.9778$, $4.2078-5.26643\,\mu$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\sim 6.34$), two low-redshift sources to the north (G1 and G2,
reported to be $z\sim 2.0$), and a fourth source only detected in rest-frame
UV emission to the southeast (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 proceed 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$\alpha$ redshifted
to $z=6.34$ (i.e., $\lambda_{obs}=4.79954\,-\,4.84467\,\mu$m, or
$z_{H\alpha}\sim 6.31-6.38$; left 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.

If we instead integrate the R100 cube over a wavelength range similar to the
HST/WFC3 F160W transmission (i.e., $\lambda_{obs}=1.4\,-\,1.6\,\mu$m), the
distribution is drastically changed (Fig. 1, right panel). Note that this
distribution is identical to that seen in previous HST/WFC3 F160W observations
(see Fig. 13), 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 define 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 use 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 field of view, which we explore
in Appendix B.

Figure 1: Integrated emission of the HFLS3 field, using two illustrative
wavelength ranges: redshifted H$\alpha$ for $z\sim 6.34$ for the R2700 cube
($\lambda_{obs}=4.79954\,-\,4.84467\,\mu$m, left panel) and the approximate
wavelength range of HST/WFC3 F160W for the R100 cube
($\lambda_{obs}=1.4\,-\,1.6\,\mu$m; right panel). The adopted source masks are
shown with white boundaries. North is up and east is to the left. The right
panel is nearly identical to the observed HST/WFC3 F160W image, but at higher
spatial resolution (see Fig. 13).

###  3.2 R2700 integrated spectra

To investigate the distribution of flux in this field, we extract spectra from
the R2700 cube using each of the spatial masks (see top panel of Fig. 1)
except ‘G2’, which will be 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\,\mu$m value:

$\rm F(\lambda_{obs})=F_{4\mu m}\left(\frac{\lambda_{obs}}{4\mu
m}\right)^{\alpha}$ (1)

The models are convolved with the line spread function (LSF)333As recorded in
the JWST documentation; https://jwst-docs.stsci.edu/jwst-near-infrared-
spectrograph/nirspec-instrumentation/nirspec-dispersers-and-filters before
comparison to the data. The velocity widths of line pairs (i.e.,
[SII]$\lambda\lambda 6716,6731$, [NII]$\lambda\lambda 6548,6584$,
[OIII]$\lambda\lambda 4959,5007$, and each Balmer line) are fixed to be
identical. We adopt the standard assumptions of [NII]$\lambda
6584$/[NII]$\lambda 6548$ $=2.94$ (e.g., Dojčinović et al. 2023) and
[OIII]$\lambda 5007$/[OIII]$\lambda 4959$ $=2.98$ (e.g., Dimitrijević et al.
2007). Following other NIRSpec IFU investigations (Übler et al. 2023; Perna et
al. 2023), we rescale 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 are estimated using the standard
error output from lmfit. Only lines with $>1\sigma$ emission at the expected
centroid wavelength are fit.

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

With a redshift of $z_{G1}=3.4805\pm 0.0001$, the northern galaxy ‘G1’ is
well-detected in H$\alpha$, [NII]$\lambda\lambda 6548,6584$,
[SII]$\lambda\lambda 6716,6731$, [SIII]$\lambda 9532$, and HeI$\lambda 10829$.
Because HeI$\lambda 10829$ may be resonantly scattered (e.g., Rudy et al.
1989; Belli et al. 2023), we allow the centroid of this line to vary
(resulting in an offset of $137\pm 19$ km s-1). In addition, [SIII]$\lambda
9532$ may be blended with Pa-$\epsilon$ (e.g., Kehrig et al. 2006), so we
allow the model of this line to feature a velocity offset ($54\pm 14$ km s-1).
The best-fitting model shows minor ($\sim 1\sigma$) residuals for the
H$\alpha$ complex, possibly suggesting the presence of an additional galaxy or
outflow.

The central component ‘C’ ($z_{C}=6.3425\pm 0.0002$) features well-detected
H$\alpha$, H$\beta$, [OIII]$\lambda 5007$, and [NII]$\lambda 6584$ emission
with broad profiles (FWHM$\sim 500-700$ km s-1). We find evidence for weak
[OIII]$\lambda 4959$, [NII]$\lambda 6548$, [SII]$\lambda\lambda 6716,6731$,
with no H$\gamma$. The [SII]$\lambda\lambda 6716,6731$ emission is broad and
low-level, so fits with a free FWHM returned unphysical values ($\sim 1500$ km
s-1). To better constrain these lines, we fixed the width of each to be
identical to that of H$\alpha$. The $4\,\mu$m continuum is strong but features
no significant slope in $\rm F_{\lambda}$. 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 Section
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$\gamma$, H$\beta$,
[OIII]$\lambda\lambda 4959,5007$, and H$\alpha$. But while the S emission is
well-fit by a single redshift ($z_{S}=6.3592\pm 0.0001$), the double-peaked
appearance of the W component suggests two galaxies ($z_{W1}=6.3549\pm 0.0001$
and $z_{W2}=6.3627\pm 0.0001$, with a velocity difference of $368\pm 7$ km
s-1) that may each be characterised directly from the integrated spectrum. No
significant [NII]$\lambda\lambda 6548,6584$ or [SII]$\lambda\lambda 6716,6731$
emission is detected in the integrated W spectrum.

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

Figure 2: Integrated spectra of the R2700 cube using the masks of Fig. 1
(magenta; excluding G2), with $1\sigma$ errors from the associated error
spectrum shown as 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.

Figure 3: (Cont.) Figure 4: (Cont.) Table 2: Best-fit redshift and continuum
flux levels for each of the HFLS3 field components.

| Redshift |  $4\,\mu$m Continuum Flux | Spectral
---|---|---|---
|  | [$10^{-20}$ erg s-1 cm-2 Å -1] | Slope
G1 | $3.4805\pm 0.0001$ | $1.26\pm 0.02$ | $-0.7\pm 0.1$
G2 | $2.00\pm 0.01$ | - | -
C | $6.3425\pm 0.0002$ | $2.54\pm 0.06$ | $-0.3\pm 0.1$
S1 | $6.3592\pm 0.0001$ | $0.25\pm 0.03$ | $-0.9\pm 0.8$
S2 | $6.3593\pm 0.0001$ | $0.25\pm 0.02$ | $3.8\pm 0.5$
W1 | $6.3549\pm 0.0001$ | $0.69\pm 0.05$ | $1.1\pm 0.4$
W2 | $6.3627\pm 0.0001$ | — | —
444Fits derived through a multi-line and continuum simultaneous fit of R2700
spectra (except for G2, which used the R100 spectrum). Spectral slope given
for $F_{\lambda}(\lambda_{obs})$ (see equation 1). Since the results of $\rm
W1$ and $\rm W2$ come from the same fit, the best-fit continuum properties
listed for $\rm W1$ describe the combined continuum.

Table 3: Best fit FWHM and integrated fluxes of observed spectral lines.

|  | H$\gamma$ | H$\beta$ | [OIII]$\lambda$4959 | [OIII]$\lambda$5007 | [NII]$\lambda$6548 | H$\alpha$ | [NII]$\lambda$6584 | [SII]$\lambda$6716 | [SII]$\lambda$6731 | [SIII]$\lambda$9532 | HeI$\lambda$10829
---|---|---|---|---|---|---|---|---|---|---|---|---
G1 | FWHM | $-$ | $-$ | $-$ | $-$ | $282\pm 85$ | $333\pm 53$ | $282\pm 85$ | $194\pm 70$ | $194\pm 70$ | $352\pm 18$ | $432\pm 39$
| I | $-$ | $-$ | $-$ | $-$ | $63\pm 13$ | $1035\pm 65$ | $184\pm 39$ | $64\pm 14$ | $71\pm 19$ | $286\pm 11$ | $402\pm 34$
C | FWHM | (500) | $679\pm 43$ | $476\pm 53$ | $476\pm 53$ | $571\pm 101$ | $679\pm 43$ | $571\pm 101$ | $679\pm 43$ | $679\pm 43$ | $-$ | $-$
| I | $<172$ | $310\pm 37$ | $165\pm 14$ | $492\pm 42$ | $144\pm 28$ | $1747\pm 190$ | $423\pm 82$ | $290\pm 70$ | $265\pm 68$ | $-$ | $-$
S1 | FWHM | $116\pm 15$ | $116\pm 15$ | $119\pm 5$ | $119\pm 5$ | $98\pm 58$ | $116\pm 15$ | $98\pm 58$ | (150) | (150) | $-$ | $-$
| I | $95\pm 14$ | $216\pm 10$ | $498\pm 9$ | $1484\pm 26$ | $22\pm 6$ | $802\pm 29$ | $64\pm 16$ | $<61$ | $<64$ | $-$ | $-$
S2 | FWHM | $125\pm 14$ | $125\pm 14$ | $65\pm 5$ | $65\pm 5$ | (150) | $125\pm 14$ | (150) | (150) | (150) | $-$ | $-$
| I | $75\pm 18$ | $130\pm 10$ | $235\pm 7$ | $701\pm 20$ | $<53$ | $458\pm 30$ | $<44$ | $<55$ | $<64$ | $-$ | $-$
W1 | FWHM | $332\pm 27$ | $332\pm 27$ | $324\pm 14$ | $324\pm 14$ | (300) | $332\pm 27$ | (300) | (300) | (300) | $-$ | $-$
| I | $92\pm 29$ | $310\pm 40$ | $452\pm 24$ | $1348\pm 73$ | $<205$ | $1223\pm 158$ | $<227$ | $<236$ | $<199$ | $-$ | $-$
W2 | FWHM | $172\pm 16$ | $172\pm 16$ | $135\pm 8$ | $135\pm 8$ | (200) | $172\pm 16$ | (200) | (200) | (200) | $-$ | $-$
| I | $59\pm 19$ | $169\pm 26$ | $555\pm 21$ | $1655\pm 64$ | $<136$ | $891\pm 120$ | $<151$ | $<157$ | $<133$ | $-$ | $-$
X1 | FWHM | (150) | $270\pm 69$ | $131\pm 28$ | $131\pm 28$ | (150) | $270\pm 69$ | (150) | (150) | (150) | $-$ | $-$
| I | $<27$ | $39\pm 12$ | $30\pm 3$ | $88\pm 9$ | 54 | $324\pm 53$ | $<51$ | $<59$ | $<55$ | $-$ | $-$

555FWHM and integrated fluxes given in units of [km s-1[] and [$10^{-20}$erg
s-1 cm-2], respectively. Values found using simultaneous fit of a power law
continuum and 1-D Gaussian for each line. For undetected lines, we list the
$3\sigma$ upper limit on the integrated flux based on the error spectrum and
an assumed FWHM (see value in parentheses). No correction for dust reddening
or gravitational magnification has been applied. We include the best-fit line
fluxes for a weaker candidate galaxy (X1, see Appendix B).

###  3.3 R100 integrated spectra

Due to the coarse LSF of the R100 data ($\sim 10^{4}$ km s-1; 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 line spread
function (e.g., Heintz et al. 2023; Jones et al. 2023; Umeda et al. 2023).
While the R100 data is therefore not useful for kinematics or precise redshift
derivations, we may examine the spectrum of ‘G2’ to find a general redshift.

‘G2’ features two clear emission lines in a wavelength regime ($\lambda_{\rm
obs}=1.4-2.0\,\mu$m) that is not covered by the R2700 data. To fit this, we
first create a higher-resolution model spectrum ($\delta\lambda=0.001\,\mu$m)
that is populated with the strongest expected lines (i.e., H$\beta$,
[OIII]$\lambda\lambda 4959,5007$, [NII]$\lambda\lambda 6548,6584$, and
H$\alpha$)666We find no evidence for significant [OII]$\lambda\lambda
3726,3729$ emission in G1.. Since the lines are unresolved, we adopt Gaussian
profiles instead of the more detailed pseudo-Voigt profiles used in the R2700
analysis. We find evidence for a spectral break at $\lambda_{\rm obs}\sim
1.1\,\mu$m, which we tentatively interpret as a Balmer break. The continuum is
thus modelled as a power law red-wards of this break ($\lambda_{\rm
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 rebinned to the R100 spectral bins. We use lmfit
with a least-square minimisation to find the best-fit continuum and line
parameters and redshift, resulting in the fit of Fig. 5. The error spectrum
was rescaled to match the standard deviation of the residual spectrum.

Figure 5: Integrated spectrum extracted from the R100 cube using the ‘G2’ mask
of Fig. 1 with illustrative $1\sigma$ 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$\beta$, [OIII]$\lambda\lambda 4959,5007$,
[NII]$\lambda 6548$, H$\alpha$, and [NII]$\lambda 6584$. The best-fit redshift
is $z_{G2}=2.00\pm 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\pm 0.01$,
which is comparable to the CIV-based redshift of Riechers et al. (2013)
($z\sim 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 $\lambda_{obs}\sim 1.72$ $\mu$m. While the line amplitude is
quite weak ($\sim 1\sigma$), the centroid wavelength of this line is in
agreement with the faint auroral line [NII]$\lambda 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 investigate 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$\alpha$ and [NII]$\lambda\lambda
6548,6584$) we instead use a parametric model. The standard assumption of a
Gaussian model resulted in some low-level residuals, so we instead assume 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 Fig. 6, 7, 8, and 9. Some weaker lines (e.g.,
[SII]$\lambda 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.

The maps for component G1 are shown in Fig. 6. For the two strongest lines
(H$\alpha$ and [NII]$\lambda 6584$), we see that the emission is circular with
a slight concentration in the northeast-southwest 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]$\lambda 9532$ and
HeI$\lambda 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. 7), 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 Figure 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 will investigate further in Section 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. 8) with continuum and H$\alpha$ intensity minima between the
galaxies. The brighter galaxy (to the southeast; S1) was originally identified
as ‘R1’ by Cooray et al. (2014).

The W component features a strong velocity gradient (Figure 9). 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\sim 3.481$; G2,
$z\sim 2.01$), a lensed source with complex kinematics suggesting two
components (C1 and C2, $z\sim 6.342$), two galaxies to the south (S1 and S2,
$z\sim 6.359$), and two galaxies to the west (W1, $z\sim 6.355$; W2, $z\sim
6.363$).

Figure 6: 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]$\lambda\lambda 6716,6731$ are not presented due to
low S/N. Due to the construction of our model, [NII]$\lambda 6548$ has the
same integrated velocity field, velocity dispersion map, and continuum map as
[NII]$\lambda 6584$, but an integrated intensity map that is a factor 2.98
lower. For [SIII]$\lambda 9532$ and HeI$\lambda 10829$, $v_{50}$ values are
given with respect to their best-fit redshifts ($z=3.4810$ and 3.4821,
respectively).

Figure 7: 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$\beta$, [OIII]$\lambda\lambda 4959,5007$, and
[SII]$\lambda\lambda 6716,6731$ are not presented due to low S/N. Due to the
construction of our model, [NII]$\lambda 6548$ has the same integrated
velocity field, velocity dispersion map, and continuum map as [NII]$\lambda
6584$, but an integrated intensity map that is a factor 2.98 lower.

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

Figure 9: Morpho-kinematic maps of the strong emission lines detected in the W
components ($z_{W1}=6.3550$, $z_{W2}=6.3628$). Integrated intensities are
given per spaxel. The corresponding maps of H$\gamma$ are not presented due to
low S/N. Due to the construction of our model, [OIII]$\lambda 4959$ has the
same velocity field, velocity dispersion map, and continuum map as
[OIII]$\lambda 5007$, but an integrated intensity map that is a factor 2.94
lower.

###  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\sim 2.1$) and used the marginally resolved PdBI [CII]158$\mu$m image of
HFLS3 to derive a low magnification factor ($\mu\sim 2.2$). With the precise
spectroscopic redshifts of G1 and G2 and a resolved map of HFLS3, we refine
this model using the public lens modelling software pyautolens
777https://github.com/Jammy2211/PyAutoLens (Nightingale & Dye 2015;
Nightingale et al. 2018, 2021).

####  4.2.1 Methods

Using pyautolens, we may derive the intrinsic (source-plane) mass and light
profiles for each source in the field of view. Each light profile is given by
a Sérsic profile (Sérsic 1963), while we assume elliptical isothermal mass
profiles. The mass and light profiles of each object are assumed to have
shared centres, axis ratios, and position angles. We model the PSF as a
circular Gaussian with FWHM$\sim 0.1^{\prime\prime}$. The fitting process
results in the best-fit centroid, intrinsic axis ratio ($q\equiv a/b$), and
position angle ($\phi$) of the mass and light profiles, the effective radius
($\rm r_{eff}$)and Sérsic index ($n$) of the light profile, and the Einstein
radius ($\rm r_{Ein}$) of the mass profile. By dividing the total fluxes of
the best-fit image-plane and source-plane models for each component, we
calculate a total magnification factor.

Due to the complexity of the field, we begin by examining the two low-redshift
sources (G1 and G2). These are isolated by collapsing the R100 data cube over
$\lambda_{\rm obs}=0.8-1.1\,\mu$m (which is not covered by the R2700 cube).
This collapsed image (top left panel of Fig. 10) has little contribution from
the $z>6$ components. Because of this, we may examine how G2 ($z\sim 2.0$, the
northernmost source in the top row of Fig. 10) lenses G1 ($z\sim 3.5$, the
elongated central source in the top row of Fig. 10). 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. 10, with parameters listed in Table 8 and source-
plane models of G1 and G2 given in Appendix D.

Figure 10: Results of pyautolens lensing analysis of collapsed images of the
R100 cube ($\lambda_{\rm obs}=0.8-1.1\,\mu$m; top row) and the R2700 cube
($\lambda_{\rm obs}=4.80-4.85\,\mu$m) assuming that the component C is
composed of a single source (middle row) or two sources (C1 to the west and C2
to the east; 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 turn to a wavelength range in the R2700 cube containing H$\alpha$ at
$z\sim 6.34$ ($\lambda_{obs}=4.80-4.85\,\mu$m) to examine whether component C
is better modelled by a single or double source. This wavelength range
contains strong H$\alpha$ emission from component C, moderate continuum from
G1, and no strong emission from G2. We model 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. 10, with parameters listed in Table 8.

Table 4: Best-fit parameters from pyautolens analysis.

|  | ($\delta$x,$\delta$y) | $q$ | $\phi$ | $\rm r_{eff}$ | $n$ | $\mu$ | $\rm r_{Ein}$
---|---|---|---|---|---|---|---|---
|  | [$"$,$"$] |  | [∘] | [$"$] |  |  | [$"$]
R100 | G1 | ($0.14\pm 0.004,0.301\pm 0.009$) | $1.43\pm 0.01$ | $57.02\pm 1.51$ | $1.028\pm 0.018$ | $5.06\pm 0.16$ | 1.2 | -
| G2 | ($0.194\pm 0.002,0.004\pm 0.002$) | $2.39\pm 0.05$ | $-30.52\pm 0.95$ | $0.251\pm 0.007$ | $1.22\pm 0.06$ | - | $0.12\pm 0.011$
R2700 (Single C) | G1 | ($0.181\pm 0.002,0.268\pm 0.004$) | $1.94\pm 0.08$ | $36.13\pm 0.85$ | $1.072\pm 0.041$ | $3.19\pm 0.1$ | 1.12 | $0.145\pm 0.008$
| G2 | - | - | - | - | - | - | $0.11\pm 0.01$
| C | ($0.13\pm 0.005,-0.017\pm 0.007$) | $2.26\pm 0.04$ | $69.65\pm 0.77$ | $0.298\pm 0.002$ | $1.19\pm 0.04$ | 1.89 | -
R2700 (Double C) | G1 | ($0.176\pm 0.003,0.298\pm 0.009$) | $1.77\pm 0.09$ | $29.2\pm 0.42$ | $2.673\pm 0.26$ | $4.73\pm 0.23$ | 1.1 | $0.169\pm 0.006$
| G2 | - | - | - | - | - | - | $0.12\pm 0.01$
| C1 | ($0.102\pm 0.005,0.012\pm 0.004$) | $2.00\pm 0.05$ | $57.1\pm 1.2$ | $0.293\pm 0.006$ | $1.34\pm 0.07$ | 2.09 | -
| C2 | ($0.293\pm 0.009,0.004\pm 0.002$) | $1.49\pm 0.04$ | $46.22\pm 4.7$ | $0.211\pm 0.019$ | $3.17\pm 0.26$ | 2.43 | -

888Parameters include the offset from central position (17h06m47.799s,
$+58^{\circ}46^{\prime}23.729^{\prime\prime}$) in arcseconds, axis ratio
($q\equiv a/b$), position angle ($\phi$), effective radius ($\rm r_{eff}$),
Sérsic index ($n$), lensing magnification ($\mu$), and Einstein radius ($\rm
r_{Ein}$). We present the lensing analysis results for a collapsed region of
the R100 cube that contains mainly emission from G1 and G2 ($\lambda_{\rm
obs}=0.8-1.1\,\mu$m, top section) and for a collapsed region of the R2700 cube
that primarily contains H$\alpha$ from C and continuum emission from G1
($\lambda_{obs}=4.80-4.85\,\mu$m). The middle section shows the results
assuming C is composed of one source, while the lower section assumes two
components. For each R2700 result, we assume that the mass profile of G2 is
identical to the best-fit light profile of the R100 image, but with a variable
Einstein radius.

####  4.2.2 Results

From this analysis, we find that the lowest-redshift source (G2) is elongated
($q\sim 2.4$), compact ($r_{\rm eff}\sim 0.25^{\prime\prime}$, or $\sim 2.1$
kpc at $z=2.00$), and not centrally peaked ($n\sim 1.2$). On the other hand,
the source-plane morphology of G1 is characterised as elongated ($q\sim
1.4-1.9$), extended ($r_{\rm eff}\sim 1.0-2.7^{\prime\prime}$, or $\sim
7.3-19.8$ kpc at $z=3.48$), centrally peaked ($n\sim 3-5$), and only
moderately magnified by G2 ($\mu\sim 1.1-1.2$).

Component C in the collapsed R2700 image appears as a curved arc with two
bright regions. We first examine whether this may be explained as a single
lensed source, yielding a best-fit source-plane morphology that is elongated
($q\sim 2.3$), compact ($r_{\rm eff}\sim 0.3^{\prime\prime}$, or $\sim 1.6$
kpc at $z=6.34$), with an exponential profile ($n\sim 1.2$), and a low
magnification factor ($\mu=1.89$). As seen in the middle row of Fig. 10, 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 assume that component C is composed of two spatially separate Sérsic
profiles (C1 to the west and C2 to the east). The best-fit model returns one
source (C1) that is similar to the similar to the best-fit single-component
model, with $q\sim 2.0$, $r_{\rm eff}\sim 0.3^{\prime\prime}$, $n\sim 1.2$,
and a higher magnification factor ($\mu=2.09$). However, this model contains
an additional component to the west that is more circular ($q\sim 1.5$),
smaller ($r_{\rm eff}\sim 0.2^{\prime\prime}$), more centrally peaked ($n\sim
3$), and is more magnified ($\mu=2.43$). The resulting image-plane morphology
of this best-fit model is double-peaked, resulting in a better representation
(i.e., a lower $\chi^{2}$). The two best-fit components have a small source-
plane projected distance ($0.2^{\prime\prime}$, or $\sim 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^{\prime\prime}$ (see Section 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$\alpha$ emission in the central component
suggests that it is composed of two spatially separate ($\sim 1$ kpc projected
offset) galaxies. By examining the best-fit morphologies (Fig. 18), 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. 7), we see that each line shows a
weak ($\delta v\sim 100$ km s-1) 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. 10), 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 which peaked in the southeast. 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\pm 0.0009$). From our best-fit redshift ($z_{C}=6.3425\pm
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 bottom-centre and bottom-right panels
of Figure 18), 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 the component C of HFLS3 is composed of two
components. The best-fit magnification of our model is also similar ($\mu\sim
2.1-2.4$, compared to $\mu=2.2\pm 0.3$).

###  4.3 Excitation conditions

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

$N2=log_{10}([NII]\lambda 6584/H\alpha)$ (2) $S2=log_{10}([SII]\lambda\lambda
6716,6731/H\alpha)$ (3) $R3=log_{10}([OIII]\lambda 5007/H\beta)$ (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. 11) for each
high-redshift ($z>6$) source in the field. We examine 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 best-fit source-
plane images in Figure 18), 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 do not apply corrections for dust reddening or gravitational
lensing.

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\sim 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]$\lambda 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\sim 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 C and S
components are composed of two sub-components each (see Section 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.

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

###  4.4 Distribution of star formation

HFLS3 was originally found to be an extreme starburst, with an SFR${}_{\rm
FIR}=2900$ M⊙ yr-1 (Riechers et al. 2013). Further studies found that the
central galaxy was lensed, resulting in a lower intrinsic $\rm SFR\sim 1300$
M⊙ yr-1 or a 100-Myr averaged $\rm SFR\sim 660$ M⊙ yr-1 (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^{\prime\prime}$,
Spitzer/IRAC with $\sim 2^{\prime\prime}$) 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 $\sim 11$ kpc? This may be explored using the H$\alpha$ luminosity
(which mainly traces star formation on timescales of $\sim 10-20$ Myr; e.g.,
Kennicutt 1998; Glazebrook et al. 1999), corrected for dust attenuation.

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

$E(B-V)_{BD}=\frac{2.5}{k(\lambda_{H\beta})-k(\lambda_{H\alpha})}log_{10}\left(\frac{F_{H\alpha,obs}/F_{H\beta,obs}}{2.86}\right)$
(5)

where $k(\lambda)$ is the assumed dust attenuation curve (Calzetti et al.
2000), $F_{H\alpha,obs}/F_{H\beta,obs}$ is the observed line flux ratio, and
we have assumed an intrinsic line ratio of 2.86 (Osterbrock 1989). We then use
this value to estimate the intrinsic line flux of H$\alpha$, assuming case B
recombination and $T_{e}=10^{4}$ K:

$F_{H\alpha,int}=F_{H\alpha,obs}10^{k(\lambda_{H\alpha})E(B-V)_{BD}/2.5}$ (6)

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

$SFR_{H\alpha}=\frac{4\pi D_{L}^{2}F_{H\alpha,int}}{1.26\times
10^{41}\mathrm{erg\,s^{-1}}}$ (7)

where $D_{L}$ is the luminosity distance. The resulting values are presented
in Table 9. 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
($\lesssim 1\sigma$ difference) values of SFR, while the high-mass, high-
redshift analogue law of Salim et al. (2018) returns slightly lower ($\lesssim
2\sigma$ difference) values of SFR.

From this analysis, we see that the central component indeed features a high
SFR. But as found in Section 4.2, this emission is lensed, resulting in a
dust- and magnification-corrected SFR of $160\pm 80$ M⊙ yr-1 (assuming an
average $\mu\sim 2.26\pm 0.17$). The combined SFR of the $z\sim 6.3$ galaxies
is $510\pm 140$ M⊙ yr-1, which is comparable to the previous 100 Myr-averaged
SFR of Cooray et al. (2014), or $654_{-90}^{104}$ M⊙ yr-1. But if S1, S2, W1,
and/or W2 contain AGN (which is not ruled out by their line ratios, see
Section 4.3), then their SFRHα 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$\alpha$) is not concentrated in
a single source but is distributed in multiple objects: a central lensed arc
that has two components ($\sim$32$\%$ of the SFR), two galaxies to the south
($\sim$17$\%$ of the SFR), and two galaxies to the west ($\sim$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 submm
observations to be confirmed), the current data imply that HFLS3 is not a
single starbursting galaxy, as previously reported.

Table 5: Observed B-V colour excess derived from Balmer-decrement, as well as
integrated H$\alpha$ flux and resulting SFR corrected for dust reddening and
gravitational magnification.

Component | E(B-V)BD | $F_{H\alpha,intrinsic}$ | SFRHα
---|---|---|---
C | $0.58\pm 0.14$ | $4.6\pm 2.3$ | $160\pm 80$
S1 | $0.22\pm 0.05$ | $1.6\pm 0.3$ | $58\pm 10$
S2 | $0.18\pm 0.09$ | $0.8\pm 0.2$ | $29\pm 9$
W1 | $0.27\pm 0.16$ | $2.8\pm 1.6$ | $100\pm 60$
W2 | $0.52\pm 0.17$ | $4.4\pm 2.8$ | $160\pm 100$
999Flux and SFR values are in units of [$10^{-17}$ erg s-1 cm-2] and [M⊙
yr-1], respectively. We assume a Salpeter (1955) IMF and Calzetti et al.
(2000) extinction law.

###  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\approx 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 $\gtrsim 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 ($\sim 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$\alpha$
emitter space density at $z>6$, which is not yet well predicted (e.g.,
Pozzetti et al. 2016). Alternatively, we may compare the observed galaxy
distribution to a cosmological simulation that includes H$\alpha$ emission and
dust extinction (e.g., Hashimoto et al. 2023). This will be deferred to a
future work.

At this time, we may 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 $\Delta
r\lesssim 20$ kpc and $\Delta$v$\lesssim 500$ km s-1 for ‘close’ pairs (e.g.,
Duncan et al. 2019b; Ventou et al. 2019; Romano et al. 2021). The six $z>6$
galaxies detected here are all located within a circle of diameter $\sim
2^{\prime\prime}$ ($\sim 11$ kpc at $z=6.34$) with a maximum redshift
difference of $\delta z\sim 0.02$ (i.e., a line of sight velocity difference
of $\sim 800$ km s-1). 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 ($\delta z\sim 0.0092$, using the redshifts of
Riechers et al. 2013), S1 and S2 ($\delta z\sim 0.0001$), and W1 and W2
($\delta z\sim 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 ($\Delta$v$=591\pm 11$ km s-1 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 $\sim 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\sim 100$) and high spectral resolution ($R\sim 2700$) data, we
find a crowded field, with two low-redshift sources (G1, $z\sim 3.4806$; G2,
$z\sim 2.00$), a central gravitationally lensed arc that is composed of two
sources (C, $z\sim 6.3425$), two close galaxies at the same redshift to the
south (S1 & S2, $z\sim 6.3592$), and two close galaxies with a velocity offset
to the west (W1, $z\sim 6.3550$; W2, $z\sim 6.3628$). All of the $z>6$
galaxies are located within an area of $\sim 2^{\prime\prime}$, or $\sim
11\,$kpc at $z\sim 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 use 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 ($\mu\sim 1.1-1.2$), while the H$\alpha$
emission of C is more strongly magnified ($\mu\sim 2.1-2.4$, comparable to the
value of $\sim 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, while we may not
rule out or confirm the presence of AGN in the S and W components.

The observed Balmer decrements are then used to derive extinction-corrected
SFRHα values for each source. This shows that the star formation is
distributed across the field ($\rm SFR_{\rm H\alpha}=510\pm 140$ M⊙ yr-1,
corrected for lensing and extinction), with the largest contribution ($\sim
51\%$) from the W galaxies. However, the presence of AGN in some sources may
inflate these SFRs.

We use 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 $\delta r<$20 kpc
and $\delta$v$<500$ km s-1, 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 ($\sim 600$ km s-1), 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\sim
6.34-6.36$ that are lensed by two foreground galaxies at $z\sim 2.1$ and
$z\sim 3.5$. All of the $z>6$ galaxies either 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\sim 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
to further characterise this unique field.

## Acknowledgements

GCJ, AJB, AJC, JC, AS acknowledges 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. MP, SA, and BRdP
acknowledge grant PID2021-127718NB-I00 funded by the Spanish Ministry of
Science and Innovation/State Agency of Research (MICIN/AEI/
10.13039/501100011033). MP also acknowledges Programa Atracción de Talento de
la Comunidad de Madrid via grant 2018-T2/TIC-11715. SC and GV acknowledge
support from the European Union (ERC, WINGS,101040227). RM, FDE, JS, and JW
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. RM also acknowledges funding from a research professorship
from the Royal Society. RB acknowledges support from an STFC Ernest Rutherford
Fellowship [grant number ST/T003596/1]. GC acknowledges the support of the
INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle
up to Cosmic Dawn with the latest generation IFU facilities”. 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

  * Allen et al. (2008) Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20
  * Arata et al. (2020) Arata, S., Yajima, H., Nagamine, K., Abe, M., & Khochfar, S. 2020, MNRAS, 498, 5541
  * Armstrong (1967) Armstrong, B. 1967, J. Quant. Spec. Radiat. Transf., 7, 61
  * Bacon et al. (2016) Bacon, R., Piqueras, L., Conseil, S., Richard, J., & Shepherd, M. 2016, MPDAF: MUSE Python Data Analysis Framework, Astrophysics Source Code Library, record ascl:1611.003
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Belli et al. (2023) Belli, S., Park, M., Davies, R. L., et al. 2023, arXiv e-prints, arXiv:2308.05795
  * Böker et al. (2022) Böker, T., Arribas, S., Lützgendorf, N., et al. 2022, A&A, 661, A82
  * Böker et al. (2023) Böker, T., Beck, T. L., Birkmann, S. M., et al. 2023, PASP, 135, 038001
  * Bouwens et al. (2015) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Cameron et al. (2023) Cameron, A. J., Saxena, A., Bunker, A. J., et al. 2023, A&A, 677, A115
  * Carnall et al. (2023) Carnall, A. C., McLeod, D. J., McLure, R. J., et al. 2023, MNRAS, 520, 3974
  * Carniani et al. (2019) Carniani, S., Gallerani, S., Vallini, L., et al. 2019, MNRAS, 489, 3939
  * Cheng et al. (2020) Cheng, C., Cao, X., Lu, N., et al. 2020, ApJ, 898, 33
  * Conselice (2006) Conselice, C. J. 2006, ApJ, 638, 686
  * Cooray et al. (2014) Cooray, A., Calanog, J., Wardlow, J. L., et al. 2014, ApJ, 790, 40
  * Curti et al. (2020) Curti, M., Mannucci, F., Cresci, G., & Maiolino, R. 2020, MNRAS, 491, 944
  * Decarli et al. (2019) Decarli, R., Walter, F., Gónzalez-López, J., et al. 2019, ApJ, 882, 138
  * D’Eugenio et al. (2023) D’Eugenio, F., Perez-Gonzalez, P., Maiolino, R., et al. 2023, arXiv e-prints, arXiv:2308.06317
  * Dimitrijević et al. (2007) Dimitrijević, M. S., Popović, L. Č., Kovačević, J., Dačić, M., & Ilić, D. 2007, MNRAS, 374, 1181
  * Dojčinović et al. (2023) Dojčinović, I., Kovačević-Dojčinović, J., & Popović, L. Č. 2023, Advances in Space Research, 71, 1219
  * Dome et al. (2024) Dome, T., Tacchella, S., Fialkov, A., et al. 2024, MNRAS, 527, 2139
  * Domínguez et al. (2013) Domínguez, A., Siana, B., Henry, A. L., et al. 2013, ApJ, 763, 145
  * Duarte Puertas et al. (2021) Duarte Puertas, S., Vilchez, J. M., Iglesias-Páramo, J., et al. 2021, A&A, 645, A57
  * Duncan et al. (2019a) Duncan, K., Conselice, C. J., Mundy, C., et al. 2019a, ApJ, 876, 110
  * Duncan et al. (2019b) Duncan, K., Conselice, C. J., Mundy, C., et al. 2019b, ApJ, 876, 110
  * Fadda et al. (2004) Fadda, D., Jannuzi, B. T., Ford, A., & Storrie-Lombardi, L. J. 2004, AJ, 128, 1
  * Fruchter & Hook (2002) Fruchter, A. S. & Hook, R. N. 2002, PASP, 114, 144
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Garn et al. (2010) Garn, T., Sobral, D., Best, P. N., et al. 2010, MNRAS, 402, 2017
  * Ginolfi et al. (2022) Ginolfi, M., Piconcelli, E., Zappacosta, L., et al. 2022, Nature Communications, 13, 4574
  * Glazebrook et al. (1999) Glazebrook, K., Blake, C., Economou, F., Lilly, S., & Colless, M. 1999, MNRAS, 306, 843
  * Hashimoto et al. (2023) Hashimoto, T., Álvarez-Márquez, J., Fudamoto, Y., et al. 2023, ApJ, 955, L2
  * Heintz et al. (2023) Heintz, K. E., Watson, D., Brammer, G., et al. 2023, arXiv e-prints, arXiv:2306.00647
  * Hopkins et al. (2006) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80
  * Jones et al. (2023) Jones, G. C., Bunker, A. J., Saxena, A., et al. 2023, arXiv e-prints, arXiv:2306.02471
  * Kauffmann et al. (2003a) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003a, MNRAS, 346, 1055
  * Kauffmann et al. (2003b) Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003b, MNRAS, 341, 33
  * Kehrig et al. (2006) Kehrig, C., Vílchez, J. M., Telles, E., Cuisinier, F., & Pérez-Montero, E. 2006, A&A, 457, 477
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541
  * Kewley et al. (2013) Kewley, L. J., Dopita, M. A., Leitherer, C., et al. 2013, ApJ, 774, 100
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
  * Laporte et al. (2015) Laporte, N., Pérez-Fournon, I., Calanog, J. A., et al. 2015, ApJ, 810, 130
  * Looser et al. (2023) Looser, T. J., D’Eugenio, F., Maiolino, R., et al. 2023, arXiv e-prints, arXiv:2302.14155
  * Marshall et al. (2023) Marshall, M. A., Perna, M., Willott, C. J., et al. 2023, A&A, 678, A191
  * Nakajima & Maiolino (2022) Nakajima, K. & Maiolino, R. 2022, MNRAS, 513, 5134
  * Nightingale & Dye (2015) Nightingale, J. W. & Dye, S. 2015, MNRAS, 452, 2940
  * Nightingale et al. (2018) Nightingale, J. W., Dye, S., & Massey, R. J. 2018, MNRAS, 478, 4738
  * Nightingale et al. (2021) Nightingale, J. W., Hayes, R. G., Kelly, A., et al. 2021, J. Open Source Softw., 6, 2825
  * Oliver et al. (2012) Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pallottini et al. (2019) Pallottini, A., Ferrara, A., Decataldo, D., et al. 2019, MNRAS, 487, 1689
  * Perna et al. (2023) Perna, M., Arribas, S., Marshall, M., et al. 2023, A&A, 679, A89
  * Popesso et al. (2023) Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526
  * Pozzetti et al. (2016) Pozzetti, L., Hirata, C. M., Geach, J. E., et al. 2016, A&A, 590, A3
  * Rauscher et al. (2017) Rauscher, B. J., Arendt, R. G., Fixsen, D. J., et al. 2017, PASP, 129, 105003
  * Reddy et al. (2020) Reddy, N. A., Shapley, A. E., Kriek, M., et al. 2020, ApJ, 902, 123
  * Riechers et al. (2013) Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Nature, 496, 329
  * Riechers et al. (2020) Riechers, D. A., Hodge, J. A., Pavesi, R., et al. 2020, ApJ, 895, 81
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001
  * Robson et al. (2014) Robson, E. I., Ivison, R. J., Smail, I., et al. 2014, ApJ, 793, 11
  * Rodrigo & Solano (2020) Rodrigo, C. & Solano, E. 2020, in XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society, 182
  * Rodrigo et al. (2012) Rodrigo, C., Solano, E., & Bayo, A. 2012, SVO Filter Profile Service Version 1.0, IVOA Working Draft 15 October 2012
  * Romano et al. (2021) Romano, M., Cassata, P., Morselli, L., et al. 2021, A&A, 653, A111
  * Rudy et al. (1989) Rudy, R. J., Cohen, R. D., Rossano, G. S., Puetter, R. C., & Chapman, S. C. 1989, ApJ, 341, 120
  * Salim et al. (2018) Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Speagle et al. (2014) Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15
  * Strait et al. (2023) Strait, V., Brammer, G., Muzzin, A., et al. 2023, ApJ, 949, L23
  * Tsukui & Iguchi (2021) Tsukui, T. & Iguchi, S. 2021, Science, 372, 1201
  * Übler et al. (2023) Übler, H., Maiolino, R., Curtis-Lake, E., et al. 2023, A&A, 677, A145
  * Umeda et al. (2023) Umeda, H., Ouchi, M., Nakajima, K., et al. 2023, arXiv e-prints, arXiv:2306.00487
  * Veilleux & Osterbrock (1987) Veilleux, S. & Osterbrock, D. E. 1987, ApJS, 63, 295
  * Venemans et al. (2019) Venemans, B. P., Neeleman, M., Walter, F., et al. 2019, ApJ, 874, L30
  * Ventou et al. (2019) Ventou, E., Contini, T., Bouché, N., et al. 2019, A&A, 631, A87
  * Wagg et al. (2014) Wagg, J., Carilli, C. L., Aravena, M., et al. 2014, ApJ, 783, 71
  * Wylezalek et al. (2022) Wylezalek, D., Vayner, A., Rupke, D. S. N., et al. 2022, ApJ, 940, L7
  * Yang et al. (2018) Yang, G., Brandt, W. N., Vito, F., et al. 2018, MNRAS, 475, 1887
  * Zakamska & Greene (2014) Zakamska, N. L. & Greene, J. E. 2014, MNRAS, 442, 784

##  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 et
al. 2016, 2021). We retrieve HST images from the STScI 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 retrieve the locations of the five closest objects in
the Gaia archive101010https://gea.esac.esa.int/archive/ (see Table 6) and
examine the locations in each HST image.

As seen in Fig. 12, the HST images show significant emission at all five Gaia
locations with small centroid offsets ($<0.2^{\prime\prime}$). From this, we
conclude that the images have been properly aligned.

Figure 12: Verification that HST images have been aligned to Gaia DR3 reference frame. Each panel contains a $1^{\prime\prime}\times 1^{\prime\prime}$ view of an HST image (see labels for each row) that is focused on a given location from the Gaia archive (see labels for each column). The Gaia position is shown as a $0.1^{\prime\prime}$ diameter red circle. The locations of each Gaia source is listed in Table 6. Table 6: Positions of five objects in the Gaia archive that are closest to HFLS3. Number | RA | Dec
---|---|---
1 | 256.7108238422430 | 58.767424792645000
2 | 256.6918696535650 | 58.7610752035986
3 | 256.68468210287400 | 58.77599807468070
4 | 256.6942094613340 | 58.78483710562300
5 | 256.7115031731110 | 58.785864102389900

###  A.2 Application of spatial shift

To shift the JWST/NIRSpec data cubes to the Gaia frame, we first create
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
colour maps in Fig. 13)111111Filter profiles were obtained from the SVO Filter
Profile Service (Rodrigo et al. 2012; Rodrigo & Solano 2020). We only use the
WFC3 filters due to a higher S/N in the IFU data over the corresponding
wavelengths.. These maps were then compared with the corresponding HST images
(see contours in Fig. 13). 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.

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

Figure 13: 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 is shown by 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.

We briefly note that Cooray et al. (2014) aligned each image to the 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 ($\sim
0.05^{\prime\prime}$), the Keck/NIRC2 and PdBI data have larger uncertainties
($\sim 0.1^{\prime\prime}$ and $0.1^{\prime\prime}-0.3^{\prime\prime}$,
respectively). When combined with the $\sim 0.1^{\prime\prime}$ 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.

Figure 14: Integrated emission of the HFLS3 field, using the same illustrative
wavelength ranges as in Fig. 1: redshifted H$\alpha$ for $z\sim 6.34$ for the
R2700 cube ($\lambda_{obs}=4.79954\,-\,4.84467\,\mu$m, left panel) and the
approximate wavelength range of HST/WFC3 F160W for the R100 cube
($\lambda_{obs}=1.4\,-\,1.6\,\mu$m; right 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.

###  B.1 Northeast extension (X1)

The first candidate emission area is visible as an elongated stretch of
emission in the left panel of Fig. 14 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. 14 shows that
multiple slices show a speckled noise pattern, possibly reflecting a low-level
calibration issue in this wavelength range.

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$\beta$,
[OIII]$\lambda\lambda 4959,5007$, and H$\alpha$, yielding a redshift of
$z_{X1}=6.3438\pm 0.0001$ (see Table 5 and Fig. 15). 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$\beta$ and the non-detection of [NII]$\lambda\lambda 6548,6584$ and
[SII]$\lambda\lambda 6716,6731$ limits our ability to characterise the source.
Its detection suggests that the HFLS3 field contains other low-level galaxy
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).

Figure 15: Integrated spectra of the R2700 cube using the masks of Fig. 14,
with $1\sigma$ errors from the associated error spectrum shown as 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.

###  B.2 Southwest clumps (X2)

Fig. 14 shows two low-level clumps of emission to the southwest 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$\alpha$ emission at $z\sim 6.36$. However, the
amplitude of this emission is quite low ($<2\sigma$) 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 do 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 ($\lambda_{obs}=0.8-1.1\,\mu$m;
see upper left panel of Fig. 10), there is a circular source of emission to
the southeast of G2 (hereafter X3). This emission is not captured by the
Sérsic mode that we fit to the data (see residual in upper right panel of Fig.
10). There is a significant emission line in this region in the R100 data at
$\lambda_{obs}=4.822\,\mu$m, but a moment zero image shows that this is simply
low-level H$\alpha$ 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. 13), 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 adopt the pseudo-Voigt profile
as implemented in lmfit. This begins with a normalised Gaussian:

$G(\lambda)=\frac{1}{\sigma\sqrt{2\pi}}e^{-(\lambda-\lambda_{C})^{2}/(2\sigma^{2})}$
(8)

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

$L(\lambda)=\frac{1}{\pi}\frac{\sqrt{2\ln(2)}\sigma}{(\lambda-\lambda_{C})^{2}+2\ln(2)\sigma^{2}}$
(9)

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

$V(\lambda)=(1-\alpha)AG(\lambda)+\alpha AL(\lambda)$ (10)

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

###  C.2 Morpho-kinematic map creation

In order to examine the morpho-kinematics of each source and line, we may
produce maps of integrated intensity, velocity offset, and velocity
dispersion. These are usually derived by extracting spectra from each spaxel
($\rm f(v)$) and measuring a normalised cumulative velocity distribution ($\rm
F(v)=\int_{-\infty}^{v}f(v^{\prime})dv^{\prime}/\int_{-\infty}^{\infty}f(v^{\prime})dv^{\prime}$;
e.g., Zakamska & Greene 2014), which is used to calculate the velocity at
which N$\%$ of the flux is captured ($\rm v_{N}$; e.g., $\rm v_{50}$) as well
as associated widths (e.g., $\rm w_{80}\equiv v_{90}-v_{10}$). Since there are
overlapping lines in some spectra, we will adopt a pseudo-Voigt-based
approach.

For each spaxel, we extract 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]$\lambda 4959$/[OIII]$\lambda 5007$,
[NII]$\lambda 6548$/H$\alpha$/[NII]$\lambda 6584$, and [SII]$\lambda
6716$/[SII]$\lambda 6731$) with identical continuum values. Line pairs (i.e.,
[OIII]$\lambda 4959$/[OIII]$\lambda 5007$, [NII]$\lambda 6548$/[NII]$\lambda
6584$, and [SII]$\lambda 6716$/[SII]$\lambda 6731$) are assumed to have
identical kinematics. We adopt the standard assumptions of [NII]$\lambda
6584$/[NII]$\lambda 6548$ $=2.94$ (e.g., Dojčinović et al. 2023) and
[OIII]$\lambda 5007$/[OIII]$\lambda 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^{2}$ (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-2], centroid wavelength $\lambda_{C}$ $[\mu$m], FWHM
$[\mu$m], and Lorentzian fraction $\alpha$), we may 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$:

$\frac{I(x,y)}{\rm erg\,s\,cm^{-2}}=A(x,y)$ (11)

The velocity field ($\rm v_{50}$) 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, $\rm v_{50}$
may be expressed for a redshifted rest-frame wavelength of $\lambda_{o}$ as:

$\frac{v_{50}(x,y)}{\rm
km\,s^{-1}}=c\left(\frac{\lambda_{C}(x,y)}{\lambda_{o}}-1\right)$ (12)

where $c[$km s${}^{-1}]$ 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 $\alpha$, the non-parametric width $w_{80}$ is dependent on $\alpha$.
To calculate this value, we generate a cumulative velocity distribution of
each best-fit model spectrum and extract $v_{10}$ and $v_{90}$. Using these:

$\frac{w_{80}(x,y)}{\rm km\,s^{-1}}=v_{90}-v_{10}$ (13)

In this process, we assume 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 Section 4.2, we used 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. 16), G2 (Fig. 17), and
C (Fig. 18).

Figure 16: Best-fit source-plane (left column) and image-plane pyautolens
models (right column) of G1 using the R100 cube ($\lambda=0.8-1.1\,\mu$m; top
row) and the R2700 cube ($\lambda=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. Figure
17: Best-fit pyautolens model of the unlensed source G2 using the R100 cube
($\lambda_{obs}=0.8-1.1\,\mu$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.

Figure 18: Best-fit source-plane (top row) and image-plane pyautolens models
(bottom row) of component C using the R2700 cube ($\lambda_{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 Section 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 19, adopting the same
assumptions of solar metallicity, low density ($\rm n\sim 0.1\,cm^{-3}$), and
no shock precursor.

Figure 19: [SII]-VO87 (top) and [NII]-BPT (lower) plots created using best-fit
line fluxes for each source (see values in Table 5), as seen in Figure 11. We
now include the model grids from MAPPINGS III (Allen et al. 2008). Red-shaded
lines display models with constant magnetic parameter (MP$=B/\sqrt{n}$ in
units of $\mu$G cm3/2), while blue-shaded lines represent models with constant
shock velocity. We assume solar metallicity, low density ($\rm n\sim
0.1\,cm^{-3}$), and no 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\sim 200$ km s-1), low-
magnetic field ($\rm MP\sim 0.001\,\mu$G cm3/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 metallicity ($Z^{\prime}=1$),
and lower metallicity models may be required for the high-redshift galaxies.