We present restframe optical data of the inlineformulatexmath idTM0001 notationLaTeXzsim 4texmathinlineformula submillimetre galaxy GN20 obtained with the JWST Near Infrared Spectrograph NIRSpec in integral field spectroscopy mode. The Hinlineformulatexmath idTM0002 notationLaTeXalphatexmathinlineformula emission is asymmetric and clumpy and extends over a projected distance of gt15 kpc. To first order the largescale ionized gas kinematics are consistent with a turbulent inlineformulatexmath idTM0003 notationLaTeXsigma sim 90texmathinlineformula km sinlineformulatexmath idTM0004 notationLaTeX1texmathinlineformula rotating disc inlineformulatexmath idTM0005 notationLaTeXvrm rotsim 500texmathinlineformula km sinlineformulatexmath idTM0006 notationLaTeX1texmathinlineformula congruent with previous studies of its molecular and ionized gas kinematics. However we also find clear evidence for noncircular motions in the Hinlineformulatexmath idTM0007 notationLaTeXalphatexmathinlineformula kinematics. We discuss their possible connection with various scenarios such as external perturbations accretion or radial flows. In the centre of GN20 we find broadline emission full width at halfmaximum inlineformulatexmath idTM0008 notationLaTeXsim 10002000texmathinlineformula km sinlineformulatexmath idTM0009 notationLaTeX1texmathinlineformula in the Hinlineformulatexmath idTM0010 notationLaTeXalphatexmathinlineformula  N II complex suggestive of fast active galactic nucleusdriven winds or alternatively of the broadline region of an active black hole. Elevated values of N II inlineformulatexmath idTM0011 notationLaTeXlambda 6583texmathinlineformulaHinlineformulatexmath idTM0012 notationLaTeXalpha  gt 0.4texmathinlineformula and of the H equivalent width EWHinlineformulatexmath idTM0013 notationLaTeXalpha gt 6texmathinlineformula  throughout large parts of GN20 suggest that feedback from the active black hole is able to photoionize the interstellar medium. Our data corroborate that GN20 offers a unique opportunity to observe key processes in the evolution of the most massive presentday galaxies acting in concert over 12 billion years ago.
['Astrophysics - Astrophysics of Galaxies']
GANIFS NIRSpec reveals evidence for noncircular motions and AGN feedback in GN20
{'No Header': ', 1-11 (2024)', 'GA-NIFS: NIRSpec reveals evidence for non-circular motions and AGN feedback in GN20': "Hannah Übler, 1 , 2 ★ Francesco D'Eugenio, 1 , 2 Michele Perna, 3 Santiago Arribas, 3 Gareth C. Jones, 4 Andrew J. Bunker, 4 Stefano Carniani, 5 Stéphane Charlot, 6 Roberto Maiolino, 1 , 2 , 7 Bruno Rodríguez del Pino, 3 Chris J. Willott, 8 Torsten Böker, 9 Giovanni Cresci, 10 Nimisha Kumari, 11 Isabella Lamperti, 3 Eleonora Parlanti, 5 Jan Scholtz, 1 , 2 Giacomo Venturi 5  \n1 Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK  \n- 2 Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HA, UK\n- 3 Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850, Madrid, Spain\n- 4 University of Oxford, Department of Physics, Denys Wilkinson Building, Keble Road, Oxford OX13RH, UK\n- 5 Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy\n- 6 Sorbonne Université, CNRS, UMR 7095, Institut d' Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France\n- 7 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK\n- 8 National Research Council of Canada, Herzberg Astronomy & Astrophysics Research Centre, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada\n- 9 European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore, MD 21218, USA\n- 10 INAF - Osservatorio Astrofisco di Arcetri, largo E. Fermi 5, 50127 Firenze, Italy\n- 11 AURA for the European Space Agency, Space Telescope Science Institute, Baltimore, MD, USA  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'We present rest-frame optical data of the 𝑧 ∼ 4 sub-millimeter galaxy GN20 obtained with JWST /NIRSpec in integral field spectroscopy (IFS) mode. The H 𝛼 emission is asymmetric and clumpy and extends over a projected distance of more than 15 kpc. To first order, the large-scale ionised gas kinematics are consistent with a turbulent ( 𝜎 ∼ 90 km/s), rotating disc ( 𝑣 rot ∼ 500 km/s), congruent with previous studies of its molecular and ionised gas kinematics. However, we also find clear evidence for non-circular motions in the H 𝛼 kinematics. We discuss their possible connection with various scenarios, such as external perturbations, accretion or radial flows. In the centre of GN20, we find broad line emission (FWHM ∼ 1000 -2000 km/s) in the H 𝛼 +[N II] complex, suggestive of fast, AGN-driven winds or, alternatively, of the broad-line region of an active black hole. Elevated values of [N/i.pc/i.pc] 𝜆 6583/H 𝛼 > 0 . 4 and EW(H 𝛼 ) > 6 Å throughout large parts of GN20 suggest that feedback from the active black hole is able to photo-ionise the interstellar medium. Our data corroborates that GN20 offers a unique opportunity to observe key processes in the evolution of the most massive present-day galaxies acting in concert, over 12 billion years ago.  \nKey words: galaxies: kinematics and dynamics - galaxies: active - galaxies: high-redshift', '1 INTRODUCTION': "More than half of the total stellar mass in the present-day Universe is situated in early-type galaxies (ETGs; Fukugita et al. 1998; Hogg et al. 2002; Bell et al. 2003; Renzini 2006). Understanding the formation of these giant systems is one of the key goals of galaxy evolution studies. With estimated ETG formation redshifts of 𝑧 > 3, the discovery of bright sub-millimeter galaxies (SMGs, see review by Blain et al. 2002) at 𝑧 ≳ 4 quickly promoted them to candidate progenitors (e.g. Simpson et al. 2014). These galaxies rapidly assemble high stellar masses through intense star formation (e.g. Hainline et al. 2009; Magnelli et al. 2013; Casey et al. 2014; Swinbank et al. 2014). Besides, some SMGs, including the galaxy studied in this work, are located within cosmic over-densities, hinting at future merger events which would help them to further increase their mass by 𝑧 = 0 (e.g.  \nDaddi et al. 2009; Walter et al. 2012; Casey et al. 2014; Oteo et al. 2018; Pavesi et al. 2018; Álvarez-Márquez et al. 2023; Arribas et al. 2023; Jones et al. 2024). Furthermore, at redshifts 𝑧 > 3 the cores of proto-clusters, which are believed to reside at the nodes of the cosmic web, are expected to be subject to cold gas inflows (e.g. Overzier 2016), a process scarcely studied observationally.  \nAfter the identification of many 𝑧 > 4 SMGs in the past decades (e.g. Smail et al. 1997; Barger et al. 1998; Hughes et al. 1998; Blain et al. 2002; Chapman et al. 2003, 2005; Pope et al. 2005; Simpson et al. 2014; Hodge & da Cunha 2020), the advent of the James Webb Space Telescope ( JWST ) has enabled, for the first time, studies of their rest-frame optical emission-line properties. Within the GA-NIFS survey 1 ('Galaxy Assembly with NIRSpec IFS'; PIs: Santiago Arribas, Roberto Maiolino), as part of the NIRSpec Instrument Science Team Guaranteed Time Observations, we have targeted  \nseveral bright SMGs at 𝑧 ≥ 4 with the integral-field spectroscopic (IFS) mode (Jakobsen et al. 2022; Böker et al. 2022), three of which have already been discussed: ALESS073.1 at 𝑧 = 4 . 76, for which the NIRSpec-IFS data reveal a heavily dust-obscured active galactic nucleus (AGN) (Parlanti et al. 2024); HLFS3 at 𝑧 = 6 . 3, which Jones et al. (2024) identify as a dense galaxy group in the process of merging; and the massive proto-cluster core SPT0311-58 at 𝑧 = 6 . 9, in which Arribas et al. (2023) find evidence for accretion from the cosmic web, inflows and mergers. In this paper, we present highresolution ( 𝑅 ∼ 2700) NIRSpec-IFS data of the rest-frame optical line-emission in GN20, a large SMG at 𝑧 ∼ 4.  \nGN20wasdetected by Pope et al. (2005) as a bright 850 𝜇 msource in the Great Observatories Origins Deep Survey Northern (GOODSN) field (see also Iono et al. 2006; Pope et al. 2006). This dusty star-forming galaxy is located within a proto-cluster environment, and is detected in 1.4 GHz continuum, [C I] and in several CO transitions (1-0, 2-1, 4-3, 5-4, 6-5, 7-6) (Daddi et al. 2009; Carilli et al. 2010; Morrison et al. 2010; Hodge et al. 2013; Cortzen et al. 2020; see also Casey et al. 2009). It has a total infrared luminosity of ∼ 2 -3 × 10 13 𝐿 ⊙ , a star-formation rate of SFR IR ∼ 1800 -3000 𝑀 ⊙ /yr, a stellar mass of 𝑀 ★ ∼ 1 . 1 -2 . 3 × 10 11 𝑀 ⊙ , and a molecular gas mass of 𝑀 mol ∼ 5 -13 × 10 10 𝑀 ⊙ (Daddi et al. 2009; Carilli et al. 2010; Hodge et al. 2012; Tan et al. 2014). A regular velocity gradient spanning about 570 km/s is observed in CO(4-3), consistent with an exceptionally large 𝑟 ∼ 4 kpc and massive rotating disc (Carilli et al. 2010). Hodge et al. (2012) find a clumpy, rotating disc in CO(2-1) with a maximum rotation velocity of 𝑣 rot , max = 575 ± 100 km/s, with a velocity dispersion of 𝜎 = 100 ± 30 km/s. Dynamical mass estimates for GN20 are in the range 2 -6 × 10 11 𝑀 ⊙ (Daddi et al. 2009; Carilli et al. 2010, 2011; Hodge et al. 2012). Intriguingly, restframe ultra-violet (UV) emission tracing young stars as observed with HST/WFC3 F105W imaging is visible only in an extended ( ∼ 9 kpc) stripe to the North-West of the centroids of the CO and 880 𝜇 m emission, beyond the extent of the dust emission (Hodge et al. 2015).  \nRecently, Colina et al. (2023) studied the rest-frame 1.1 𝜇 mimaging of GN20 obtained with the MIRI instrument onboard JWST (see also Crespo Gómez et al. 2024). Their analysis reveals a twocomponent stellar structure composed of an unresolved nucleus offset by 1 kpc from the centre of an extended disc with 𝑅 𝑒 = 3 . 60 kpc. They argue that the offset nucleus may be a result of tidal interactions with other proto-cluster members, or indicate a late-stage merger. Colina et al. (2023) find that the stellar nucleus coincides with the centre of far-infrared continuum emission tracing dust-obscured star formation, while the extended stellar envelope overlaps with the cold molecular gas distribution. Analysing MIRI/MRS-IFU observations, Bik et al. (2024) find clumpy Pa 𝛼 emission out to a radius of 6 kpc. The kinematics of the Pa 𝛼 emission are consistent with a rotating disc, with a maximum rotation velocity of 𝑣 rot = 550 ± 40 km/s and an upper limit on the flux-weighted velocity dispersion ( 𝜎 𝑚 ) of 𝜎 𝑚 = 145 ± 53 km/s. Comparing the unobscured star-formation rate derived from the integrated Pa 𝛼 flux, SFR Pa 𝛼 = 144 ± 9 𝑀 ⊙ /yr, to the infrared star-formation rate SFR IR (Tan et al. 2014), Bik et al. (2024) infer a high average extinction of A 𝑉 = 17 . 2 ± 0 . 4 mag (see also Maseda et al. 2024). Crespo Gómez et al. (2024) infer a lower A 𝑉 ∼ 1 . 5 mag from spectral energy distribution (SED) fitting. They attribute this difference to the presence of either stellar populations older than 10 Myr, or of a buried AGN.  \nWedescribe the JWST /NIRSpec-IFS high-resolution observations of GN20 and our analysis methods in Sections 2 and 3. We present the complex H 𝛼 emission revealed by NIRSpec in Section 4. In Section 5 we discuss the H 𝛼 kinematics and evidence for non-circular motions  \nbased on both our data and dynamical modelling. In Section 6 we show evidence for the presence of an AGN in the centre of GN20. We summarise our findings in Section 7.  \nThroughout this work, we assume a flat Λ CDM cosmology with Ω 𝑚 = 0 . 315 and 𝐻 0 = 67 . 4 km/s/Mpc (Planck Collaboration et al. 2020). With this cosmology, 1 '' corresponds to a transverse distance of 7.07 proper kpc at 𝑧 = 4 . 05.", '2 NIRSPEC-IFU OBSERVATIONS AND DATA PROCESSING': "GN20 was observed in NIRSpec-IFS mode as part of the GA-NIFS survey under programme 1264. The NIRSpec data were taken on February 10, 2023, with a medium cycling pattern of four dither positions and a total integration time of about 2 h with the high-resolution grating/filter pair G395H/F290LP, covering the wavelength range 2 . 87 -5 . 14 𝜇 m(spectral resolution 𝑅 ∼ 2000 -3500; Jakobsen et al. 2022), and about 1 h with PRISM/CLEAR ( 𝜆 = 0 . 6 -5 . 3 𝜇 m, spectral resolution 𝑅 ∼ 30 -300). Within the same programme, MIRI (Rieke et al. 2015; Wright et al. 2015) imaging and spectroscopy of GN20 was obtained. The MIRI imaging data was presented by Colina et al. (2023), Crespo Gómez et al. (2024), and the MRS data by Bik et al. (2024).  \nRaw data files were downloaded from the Barbara A. Mikulski Archive for Space Telescopes (MAST) and subsequently processed with the JWST Science Calibration pipeline 2 version 1.11.1 under the Calibration Reference Data System (CRDS) context jwst\\_1149.pmap. We made several modifications to the default reduction steps to increase data quality, which are described in detail by Perna et al. (2023) and which we briefly summarise here. Countrate frames were corrected for 1 / 𝑓 noise through a polynomial fit. During calibration in Stage 2, we removed regions affected by failed open MSA shutters. We also removed regions with strong cosmic ray residuals in several exposures. Remaining outliers were flagged in individual exposures using an algorithm similar to /l.pc/a.pc/c.pc/o.pc/s.pc/m.pc/i.pc/c.pc (van Dokkum 2001): we calculated the derivative of the count-rate maps along the dispersion direction, normalised it by the local flux (or by three times the rms noise, whichever was highest), and rejected the 95 th percentile of the resulting distribution (see D'Eugenio et al. 2023, for details). The final cube was combined using the 'drizzle' method. The main analysis in this paper is based on the combined cube with a pixel scale of 0 . 05 '' . We used spaxels away from the central source and free of emission features to perform a background subtraction.  \nIn Figure 1, we show a line map of the emission integrated in the range 3 . 313 -3 . 327 𝜇 m, roughly encompassing the H 𝛼 and [N II] 𝜆 6583 lines at 𝑧 ∼ 4 . 055. This will be further discussed in Section 4. In this paper, we focus on the analysis of the high-resolution data. However, to provide some context for our analysis, in Figure 2 we use the prism observations to create a false-colour multiwavelength image of GN20, combining the H 𝛼 +[N II] 𝜆𝜆 6548 , 6583 emission (green) with the emission at 0 . 9 -1 . 2 𝜇 m (blue), roughly corresponding to the HST /WFC3 F105W filter tracing the emission of unobscured young stars, and the emission at 5 . 0 -5 . 3 𝜇 m (red), covering the blue part of the MIRI/F560W filter and tracing stellar mass. The astrometry in Figure 1 is derived by registering our observations to the F560W and F105W images, which are in turn  \nFigure 1. Line map of the emission integrated in the range 3 . 313 -3 . 327 𝜇 m, roughly encompassing the H 𝛼 +[N II] 𝜆 6583 lines at 𝑧 ∼ 4 . 055. GN20 is visible as the extended central source with a complex and clumpy morphology. Weindicate two clumps to the North-West and extended emission in the SouthEast with dashed outlines (see main text for details). We detect a second galaxy to the North-West of GN20, which we call GN20b, at a projected separation of about 12 kpc and a velocity difference of about +750 km/s, respectively a redshift of 𝑧 = 4 . 06532. The white circle indicates the approximate PSF FWHM( GLYPH<31> = 0 . 12 '' ).  \n<!-- image -->  \nregistered to Gaia DR3 (Gaia Collaboration et al. 2023). Based on the positional shifts derived from the MIRI and HST data, and the varying pixel sizes of the observations ranging from 0.05' to 0.1', we estimate an uncertainty of about 0.1' on the astrometry.", '3.1 Emission line fitting and H 𝛼 redshift': 'To analyse the emission-line properties and kinematics of GN20, we first fit a one-component Gaussian model to the H 𝛼 and [N/i.pc/i.pc] 𝜆𝜆 6548 , 6583 emission lines in the grating data, including a constant continuum. Because the signal-to-noise ratio ( 𝑆 / 𝑁 ) varies substantially throughout GN20 due to non-uniform dust obscuration that is particularly high in the central region (Hodge et al. 2015; Colina et al. 2023; Bik et al. 2024), we derive a Voronoi-binned map, using the algorithm /v.pc/o.pc/r.pc/b.pc/i.pc/n.pc (Cappellari & Copin 2003), as implemented in QF/i.pc/t.pc/s.pcV/i.pc/e.pc/w.pc. We first define a mask based on visual inspection of the line emission in the cube, and then require 𝑆 / 𝑁 ≥ 15 in the previously fitted H 𝛼 flux map to derive the Voronoi map. We derive H 𝛼 velocity and velocity dispersion maps by repeating our one-component Gaussian model fitting on the binned data. The resulting maps are shown in the top row of Figure 3.  \nWhile the line emission is relatively narrow in the outer regions of GN20 ( 𝜎 obs ∼ 50 -100 km/s), the H 𝛼 +[N/i.pc/i.pc] emission in the centre is broad. We identify the central regions with broad emission through visual inspection (see aperture indicated in the bottom panels of Figure 5, diameter of ∼ 4 kpc). We fit the spectra extracted from the central regions based on the Voronoi tessellation with a set of narrow components plus a set of broad components for H 𝛼 and [N/i.pc/i.pc] 𝜆𝜆 6548 , 6583. Here, we require FWHM broad > 400 km / s > FWHM narrow . We stitch the narrow component results to our initial one-component maps. As can be seen from the example spectra in Figure 5, the narrow components in the central region are distinct from the underlying broad emission. It is therefore unlikely  \nFigure 2. Comparisonofrest-frame optical, UV, and near-infrared emission in GN20, obtained through collapsing the NIRSpec prism cube in three different wavelength regions (not PSF-matched). The H 𝛼 +[N II] emission is shown in green, blue colors trace emission at 0 . 9 -1 . 2 𝜇 m(rest-frame 0 . 18 -0 . 24 𝜇 m; roughly corresponding to the HST /WFC3 F105W filter), and red colors show emission at 5 . 0 -5 . 3 𝜇 m(rest-frame 0 . 99 -1 . 05 𝜇 m; covering the blue part of the MIRI/F560W filter). The nuclear H 𝛼 +[N II] emission overlaps with the nucleus detected by MIRI (Colina et al. 2023; Crespo Gómez et al. 2024), while the diffuse rest-frame near-infrared emission extends until the outer H 𝛼 +[N II] loop. We note that the Pa 𝛼 emission detected by MIRI has a similar extent (Bik et al. 2024). The rest-frame UV emission overlaps with the Western region detected in H 𝛼 +[NII], including the location of the NorthWestern clumps.  \n<!-- image -->  \nthat we over-subtract emission associated with the narrow component through this procedure. The resulting maps are shown in the bottom row of Figure 3. In Section 6 we discuss alternative approaches to model the broad nuclear emission.  \nTo measure a redshift for GN20 from the H 𝛼 data, we rely on the kinematics, since the one-dimensional, integrated emission line profile is skewed by the bright, off-centre region in the North-West. After initially deriving the H 𝛼 kinematics as described above assuming 𝑧 = 4 . 055 (Daddi et al. 2009), we shift the systemic velocity such that we find comparable absolute maximum and minimum velocities, and the (one-component) dispersion peak coincides with 𝑣 obs ∼ 0 km/s along the kinematic major axis. We adopt a velocity shift of +75 km/s. Given the median uncertainties on the fitted velocities of 11 km/s (including in the regions of minimum and maximum velocities), this corresponds to a redshift of 𝑧 H 𝛼 = 4 . 05374 ± 0 . 00075. 3', '3.2 Dynamical modelling': "To construct dynamical models, we use DysmalPy (Davies et al. 2004a,b, 2011; Cresci et al. 2009; Wuyts et al. 2016; Genzel et al. 2017; Lang et al. 2017; Übler et al. 2018; Price et al. 2021; Genzel et al. 2023; L. L. Lee et al. subm.), a 3D forward-modelling code that takes into account the instrumental effects of beam-smearing, line broadening, and finite spatial resolution. We build a mass model  \n3 Asdescribed in the Section 3.2, in our dynamical modelling we allow for an additional velocity shift of ± 100 km/s. Our fiducial model has a velocity shift of + 28 ± 2 km/s. This would correspond to a redshift of 𝑧 H 𝛼, model = 4 . 05325 ± 0 . 00075. However, we caution that the velocity shift can be degenerate with the adopted centre (fixed in our model). In addition, as discussed in Section 5, the model does not capture all kinematic features of the data.  \nFigure 3. Top: Voronoi-binned maps (target 𝑆 / 𝑁 ≥ 15) of the observed H 𝛼 flux (left), velocity (middle) and velocity dispersion (right) from a one-component fit. In the bottom row, we show corresponding maps where we have stitched the central region including only the narrow component of a two-component Gaussian fit (see Section 3). The large-scale kinematics are consistent with a rotating disc. The high dispersion values in the central region in the top right panel trace a nuclear broad flux component. Removing this component, we find high average values of 𝜎 obs ∼ 50 -150 km/s throughout GN20 (bottom right panel). We also find evidence for non-circular motions in the velocity field through patterns that deviate from the classical spider diagram expected for rotating discs: to guide the eye, we mark 𝑣 obs ∼ 0 km/s as a dashed line in the bottom middle panel (repeated in the bottom right panel for reference).  \n<!-- image -->  \nfor GN20 informed through existing multi-wavelength constraints, in particular the recent MIRI observations at mid-infrared wavelengths that trace the stellar light distribution (Colina et al. 2023). The model includes a baryonic disc and bulge component together with a spherical dark matter halo. For the baryonic mass distribution, we account for a finite flattening following Noordermeer (2008).  \nWe fix the following structural parameters to the best-fit values derived by Colina et al. (2023): the disc effective radius 𝑅 𝑒, disc = 3 . 6 kpc, the disc Sérsic index 𝑛 S , disc = 0 . 42, the minor-to-major axis ratio 𝑏 / 𝑎 = 0 . 80. We adopt the convention that in the face-on case, the galaxy rotates in counter-clockwise direction for 𝑖 = 0 · , and in clockwise direction for 𝑖 = 180 · . Assuming a thick disc with ratio of scale height to scale length of 𝑞 0 = 0 . 2 (e.g. Wuyts et al. 2016; Genzelet al. 2017), 𝑏 / 𝑎 = 0 . 80corresponds to an inclination of either 38 · (projected counter-clockwise rotation; see also Crespo Gómez et al. 2024) or 142 · (projected clockwise rotation). Assuming that the regions in the North-West which are brightest in H 𝛼 and UV emission are closer to us (the 'near side'), this would suggest that GN20 rotates in a clockwise direction (see also stellar light distribution; Colina et al. 2023). We therefore fix 𝑖 = 142 · . We include a Gaussian prior on the total baryonic mass centred on log ( 𝑀 bar , tot / 𝑀 ⊙ ) = 11 . 4, with a standard deviation of 0.3 and bounds of [10.4; 12.4]. This is motivated by estimates of the stellar mass ( 𝑀 ★ = 1 . 1 × 10 11 𝑀 ⊙ ; Tan et al. 2014) and of the molecular gas mass ( 𝑀 H2 = 1 . 3 × 10 11 𝑀 ⊙ ; Hodge et al. 2012, here assuming 𝛼 CO = 0 . 8). Colina et al. (2023) identify an unresolved nuclear component (upper size limit 0.8 kpc). The authors conclude that this component corresponds to an obscured nuclear starburst, but they also discuss the alternative options of a massive stellar bulge or an AGN. To account in our model for this central mass component, we use their stellar mass estimate derived under the hypothesis of a bulge, 2 . 5 × 10 10 𝑀 ⊙ . This is about 10 per cent of our total baryonic mass estimate. We measure the effective radius of  \nthe nuclear component through a two-dimensional Gaussian fit to the H 𝛼 +[N II] linemap, and find 0.76-0.80 kpc, consistent with the upper limit derived by Colina et al. (2023) from the MIRI F560W imaging. We model a round bulge with 𝑅 𝑒, bulge = 0 . 8 kpc, 𝑛 S , bulge = 4 and a bulge-to-total baryonic mass ratio of 𝐵 / 𝑇 = 0 . 1. We include a Navarro et al. (1996) dark matter halo with a concentration parameter of 𝑐 = 3 . 1 (the expected value for a halo of log ( 𝑀 halo / 𝑀 ⊙ ) = 13 at this redshift, following Dutton & Macciò 2014). We fit for the dark matter fraction within the effective radius, 𝑓 DM ( 𝑟 < 𝑅 𝑒 ) , from which the total halo mass is inferred (see Price et al. 2021). We use a flat prior 𝑓 DM ( < 𝑅 𝑒 ) = [ 0; 1 ] . Further to the priors described above, we include a flat prior on the intrinsic velocity dispersion, 𝜎 0 = [ 30; 150 ] km/s, which is assumed to be constant and isotropic throughout the disc. In addition to the total baryonic mass 𝑀 bar , the dark matter fraction at the effective radius 𝑓 DM ( < 𝑅 𝑒 ) , and the intrinsic velocity dispersion 𝜎 0 , we allow the position angle to vary in the range PA=[0 · ; 50 · ], and we allow the model to adjust the systemic velocity by ± 100 km/s.  \nWe model the point-spread function (PSF) as a Gaussian with FWHM = 0 . 12 '' , corresponding to the approximate PSF (along slicers) at the wavelength of H 𝛼 (see D'Eugenio et al. 2023). For the instrumental dispersion we adopt 𝜎 inst = 57 km/s, which corresponds to the nominal spectral resolution at the wavelength of H 𝛼 (Jakobsen et al. 2022).  \nAs input for our fiducial DysmalPy runs, we use the Voronoi maps of the narrow H 𝛼 velocity and velocity dispersion. We fix the centre of the model such that it falls onto the Voronoi bin with 𝑣 obs ∼ 0 km/s along the axis connecting the observed velocity minimum and maximum. This coincides roughly with the centre of the disc component identified by Colina et al. (2023). To account for the Voronoi binning during the minimization, we scale the uncertainties by a factor equal to the square root of the number of spaxels per bin.  \nWefind the best fit through MCMC sampling, using 200 walkers and 500 steps after a burn-in of 200 steps (Foreman-Mackey et al. 2013, emcee) . Our final chains are longer than 10 times the autocorrelation time for the individual fit parameters. For visual comparison of the best-fit model to our data in the figures below, we show the median values of the model per Voronoi bin.  \nAs motivated in Section 6, we further construct a model with the same parameters as described above, but which additionally includes a uniform planar, radial inflow. The inflow velocity is not fitted for, but added as a fixed parameter (see Price et al. 2021). The preferred inflow velocity is initially identified through a grid search by injecting different values of 𝑣 𝑟 in steps of 10 km/s from zero to 200 km/s, using least-squares minimization (Markwardt 2009). We repeat the fit with the so identified preferred inflow velocity ( 𝑣 𝑟 = 130 km/s) using MCMC. The results of this model are presented in Appendix A.", '4 COMPLEX H 𝛼 EMISSION': 'In Figure 1 we show a line map of the emission integrated in the range 3 . 313 -3 . 327 𝜇 m in our NIRSpec observations. This range covers the H 𝛼 +[N II] 𝜆 6584 flux in GN20 (extended central object). Within the IFS field of view, we detect another galaxy North-West of GN20, which we call GN20b. It is at a projected separation of about 12 kpc from the centre of GN20 (also seen in the HST data; see Figure 2) and a velocity difference of about +750 km/s, respectively a redshift of 𝑧 H 𝛼, GN20b = 4 . 06532 ± 0 . 00003. GN20b is also detected in continuum in the prism observations.  \nThe map for GN20 reveals a complex, clumpy, loop-like structure. The emission is brightest in the North-Western region, close to the location of the rest-frame UV emission tracing young stars, previously detected with HST (see Figure 2). Further to the North-West of this region, we detect two fainter clumps (dashed circles in Figure 1). In the remaining regions of GN20, the H 𝛼 emission is much fainter. A central flux concentration largely overlaps with the nuclear component identified in the MIRI imaging (Colina et al. 2023, the black-edged white circle in their Figure 2), which is also prominent in the NIRSpec data above ∼ 5 𝜇 m(see Figure 2). Immediately North of this region H 𝛼 is barely detected. In the South-East, we detect clumpy H 𝛼 emission extending to the South (dashed ellipse in Figure 1).  \nThe peculiar surface brightness distribution of the H 𝛼 emission in GN20 could spark doubts regarding its disc-like nature. However, we know from previous observations in the rest-frame far-infrared continuum that the central regions of GN20 are very dusty (Hodge et al. 2015). This has also been shown through the recent analyses by Bik et al. (2024); Crespo Gómez et al. (2024), finding high attenuation in the central parts of the galaxy. In addition, the analysis of the stellar light distribution from MIRI imaging is consistent with a fairly smooth stellar disc plus a compact nuclear component (Colina et al. 2023; Crespo Gómez et al. 2024). Clumpy H 𝛼 emission originating from smooth, disc-like mass distributions has also been observed at lower redshift (Wuyts et al. 2012). Notably, most of these galaxies show smooth H 𝛼 velocity fields consistent with disc rotation. Indeed, this appears to be case as well for GN20: the H 𝛼 kinematics, much like previous studies of GN20 utilising CO and Pa 𝛼 as kinematic tracers, support the interpretation of a large rotating disc. We discuss the large-scale rotation and deviations from circular motions in Section 5.  \nThe two H 𝛼 clumps in the North-West appear also visible in CO(21) and in the UV (Hodge et al. 2015), and there is also faint emission seen in the MIRI F560W image, on the level of a few per cent of the peak flux (see Colina et al. 2023, and Figure 2). Some F560W flux  \nextends from the main disc along the same direction of the SouthEast extension (see Colina et al. 2023). Clumpy UV emission is also seen here, however this likely stems from a 𝑧 ∼ 1 . 74 galaxy identified through rest-frame optical line emission in our prism data. The H 𝛼 clumps and extension could resemble real substructure embedded within GN20, e.g. young star-forming regions in the outer disc or a protruding spiral arm. Alternatively, they could correspond to fresh material in the process of accretion onto GN20. We note faint H 𝛼 flux between GN20 and GN20b, which could further indicate a past or ongoing interaction.', '5 LARGE-SCALE ROTATION AND NON-CIRCULAR MOTIONS': 'In the middle and right panels of Figure 3 we show H 𝛼 velocity and velocity dispersion maps based on our fits to the Voronoi-binned cube (top: one-component; bottom: stitched narrow component, see Section 3 for details). We find a large velocity gradient across GN20 of about Δ 𝑣 obs ∼ 610 km/s, and observed velocity dispersion values of 𝜎 obs ∼ 50 -150 km/s in the Voronoi bins in the outer disc. The large-scale H 𝛼 velocity field is consistent with disc rotation, and the kinematic major axis, observed velocities and velocity dispersions are in general agreement with earlier results based on molecular gas kinematics (Carilli et al. 2010; Hodge et al. 2012). However, the data reveal clear deviations from circular motions in the velocity field. In the bottom middle panel of Figure 3, we indicate 𝑣 obs ∼ 0 km/s by dashed lines, illustrating a twist of the apparent kinematic minor axis. Along the North-Western 𝑣 obs ∼ 0 km/s twist we observe slightly elevated velocity dispersions, as evident from the bottom right panel of Figure 3.  \nIn the following, we further investigate the deviations from circular motions through residual analysis of a pure circular motion dynamical model. We also compare our modelling results to previous studies of the kinematics in GN20 based on other kinematic tracers.', '5.1 Dynamical modelling results and comparison to other kinematic tracers': 'We fit a dynamical model including a rotating disc, bulge, and dark matter halo, as described in Section 3, to our data. The best-fit model and residuals are shown in Figure 4. The model provides a reasonable fit to the data, with median offsets in velocity and velocity dispersion of Δ 𝑣 med = 5 . 5 km/s and Δ 𝜎 med = 14 . 9 km/s. From this model, we measure a rotation velocity 𝑣 rot ( 𝑅 𝑒 = 3 . 6kpc ) = 469 km/s (circular velocity 𝑣 𝑐 ( 𝑅 𝑒 ) = 496 km/s), reaching a maximum value of 𝑣 max = 531 km/s ( 𝑣 𝑐, max = 574 km/s) at 𝑟 = 6 . 2 kpc ( 𝑟 = 6 . 8 kpc). The circular velocity accounts for the effects of pressure support from turbulent motions, and is defined as 𝑣 2 𝑐 ( 𝑟 ) = 𝑣 2 rot ( 𝑟 ) + 2 𝜎 2 0 𝑟 / 𝑅 𝑑 , with 𝑅 𝑑 being the disc scale length (see Burkert et al. 2010, 2016). We constrain an intrinsic velocity dispersion of 𝜎 0 = 89 km/s. For a galaxy at 𝑧 ∼ 4 . 055, from the relation by Übler et al. (2019), calibrated based on ground-based observations up to 𝑧 ∼ 3 . 5, we would expect an intrinsic ionised gas velocity dispersion of 𝜎 0 ∼ 63 ± 19 km/s, somewhat lower than our best-fit result. However, we note that some regions in the outer disc of GN20 show comparable dispersion values in 𝜎 obs . Evaluating the ratio of maximum rotation velocity to intrinsic velocity dispersion, we find 𝑣 rot ( 𝑅 𝑒 )/ 𝜎 0 = 5 . 3 ( 𝑣 rot , max / 𝜎 0 = 6 . 0). The dynamical mass enclosed within one 𝑅 𝑒 is log ( 𝑀 dyn ( < 𝑅 𝑒 )/ 𝑀 ⊙ ) = 11 . 3, and within two 𝑅 𝑒 is log ( 𝑀 dyn ( < 2 𝑅 𝑒 )/ 𝑀 ⊙ ) = 11 . 7. The inferred values from our best-fit model are reported in Table 1.  \nFigure 4. Our best-fit dynamical model including a disc, bulge and halo. Top (bottom) panels show the observed velocity (velocity dispersion) field (left), the best-fit model (middle left), residuals (data-model; middle right), and goodness-of-fit (data-model)/uncertainties (right). The cross and line indicate the centre and best-fit PA. For comparison, the dotted grey line in the top left panel indicates the axis connecting the observed velocity minimum and maximum. The model provides a reasonable fit to the data, with a median velocity offset of Δ 𝑣 med = 5 . 5 km/s and a rms velocity difference of Δ 𝑣 rms = 49 . 3 km/s. Corresponding values for the velocity dispersion are Δ 𝜎 med = 14 . 9 km/s and Δ 𝜎 rms = 34 . 7 km/s. Yet, we observe strong residuals particularly in the North-Western region, indicating a disturbance of the kinematics for instance through accretion or external perturbations.  \n<!-- image -->  \nTo compare our modelling results to previous results on the kinematics of GN20 obtained in the literature, we repeat our dynamical modelling with the inclination fixed to the values adopted in the literature studies. Assuming 𝑖 = 45 · , Carilli et al. (2010) find a rotation velocity of 𝑣 rot ( 𝑟 ∼ 4kpc ) = 570 km/s based on CO(65) molecular gas kinematics, with an enclosed dynamical mass of 𝑀 dyn ( 𝑟 < 4kpc ) = 3 × 10 11 𝑀 ⊙ . Repeating our modelling with the inclination fixed to 𝑖 = 135 · (i.e. corresponding to their 𝑖 = 45 · ), we find a lower value of 𝑣 rot ( 𝑟 ∼ 4kpc ) = 442 km/s, but a comparable enclosed total mass of 𝑀 dyn ( 𝑟 < 4kpc ) = 2 × 10 11 𝑀 ⊙ .  \nHodge et al. (2012) and Bik et al. (2024) adopt an inclination of 𝑖 = 30 · . From the CO(2-1) kinematics, Hodge et al. (2012) inferred 𝑣 max = 575 ± 100 km/s and 𝜎 0 = 100 ± 30 km/s. Bik et al. (2024) found 𝑣 max ( 𝑟 < 4 kpc ) = 550 ± 40 km/s and 𝜎 𝑚 = 145 ± 53 km/s from Pa 𝛼 kinematics, where 𝜎 𝑚 represents an upper limit on the intrinsic flux-weighted velocity dispersion. Repeating our modelling with the inclination fixed to 𝑖 = 150 · (i.e. corresponding to their 𝑖 = 30 · ), we find a rotation velocity of 𝑣 rot ( 𝑅 𝑒 = 3 . 6kpc ) = 551 km/s ( 𝑣 max = 646 km/s), and 𝜎 0 = 89 km/s. Our modelling results obtained by fixing 𝑖 = 150 · are in good agreement with what has been inferred from these studies analysing molecular gas kinematics and Pa 𝛼 kinematics on similar radial scales. Comparing the amount of rotational support, Hodge et al. (2012) find 𝑣 max / 𝜎 0 | CO ( 2 -1 ) = 5 . 8, and Bik et al. (2024) find a somewhat lower value of 𝑣 max / 𝜎 𝑚 | Pa 𝛼 = 3 . 8 ± 1 . 4, however based on the upper limit on the flux-weighted velocity dispersion. Using 𝑖 = 150 · , we find even higher values of 𝑣 max / 𝜎 0 = 7 . 3, yet at the effective radius (fixed to 𝑅 𝑒 = 3 . 6 kpc), we find 𝑣 rot ( 𝑅 𝑒 )/ 𝜎 0 = 6 . 2. This is in general agreement with the literature results. Hodge et al. (2012) infer a dynamical mass of 𝑀 dyn = 5 . 4 ± 2 . 4 × 10 11 𝑀 ⊙ . Fixing 𝑖 = 150 · , we find 𝑀 dyn ( 𝑟 < 𝑅 𝑒 ) = 2 . 6 × 10 11 𝑀 ⊙ and 𝑀 dyn ( 𝑟 < 2 𝑅 𝑒 ) = 7 . 0 × 10 11 𝑀 ⊙ , in broad agreement with the results by Hodge et al. (2012).  \nMulti-phase gas kinematic measurements at high redshift are still rare. Recently, Parlanti et al. (2024) measured 𝑣 rot ( 𝑅 𝑒 = 3kpc )/ 𝜎 0 | H 𝛼 ∼ 9 . 2 in an obscured AGN at 𝑧 = 4 . 76, with a higher value of 𝑣 rot ( 𝑅 𝑒 )/ 𝜎 0 | [ CII ] ] ∼ 15 . 8 measured from the [C II]  \nTable 1. Results from our best-fit dynamical model including a thick disc, bulge and dark matter halo. For 𝜎 0 , we adopt minimum uncertainties of 10 km/s. The effective radius is fixed to 𝑅 𝑒 = 3 . 6 kpc.  \nline (Lelli et al. 2021). In this case, the authors speculate that the AGN might have deposited energy in the ionised gas phase, leading to higher dispersion values and therefore lower rotational support. In SPT0311-58, which is at the core of a massive protocluster, Arribas et al. (2023) find 𝑣 rot , max / 𝜎 0 | [ OIII ] ∼ 1 based on [O III] 𝜆 5007 kinematics. Their median intrinsic velocity dispersion is 𝜎 0 = 113 ± 19 km/s, and in some regions 2-3 times higher than the velocity dispersion measured in [C II]. This likely indicates that SPT0311-58 at 𝑧 = 6 . 9 is in an earlier phase of disc formation compared to GN20.', '5.2 Deviations from circular motions': "Despite the relatively good fit of our circular motion model to the H 𝛼 kinematics of GN20, and the agreement with literature results, we find strong residuals in particular in the North-Western region of the galaxy. This can be seen from the middle-right and right panels in Figure 4. These residuals are close to the H 𝛼 clumps identified in Figure 1. Furthermore, the best-fit kinematic major axis appears  \ntilted with respect to the velocity minimum and maximum in the data (see top left panel in Figure 4). Deviations from circular motions in disc galaxies can be associated with a variety of phenomena: perturbations through past or ongoing interaction with neighbours, smooth accretion, streaming motions due to non-axisymmetric sub-structure like a bar or spiral arms, or outflows (e.g. van der Kruit & Allen 1978; Roberts et al. 1979; van Albada & Roberts 1981; Shlosman et al. 1989; Athanassoula 1992; Wada & Habe 1992; Bournaud & Combes 2002; Binney & Tremaine 2008; Tsukui et al. 2024).  \nGN20 is known to be located within a proto-cluster environment together with the AGN GN20.2a and GN20.2b, located at projected distances of about ∼ 150 kpc and ∼ 200 kpc (Daddi et al. 2009). Based on the projected separations and stellar masses of GN20.2a and GN20.2b, we do not expect a significant effect on the kinematics of GN20: we find a tidal strength parameter (Dahari 1984; Verley et al. 2007) of 𝑄 = -3 . 9, while sizeable effects are generally expected for 𝑄 ≥ -2. However, as mentioned in Section 4, we detect another galaxy (GN20b) to the North-West of GN20 which could be involved in an ongoing or past interaction. It is therefore possible that deviations from circular motions in GN20 are triggered by past or ongoing interactions with neighbours. This could also potentially explain the fairly large velocity dispersions found in some regions in GN20. We note that offset nuclear emission detected in the MIRI imaging by Colina et al. (2023) and Crespo Gómez et al. (2024) may also suggest a past gravitational interaction.  \nAs mentioned in Section 3 and further discussed in Section 6, we detect broad emission in the central region of GN20 which may be associated with high-velocity outflows. However, the impact of these non-circular motions is already removed in the kinematic maps used for the dynamical modelling (see Figure 3). Still, it is possible that past feedback activity in the centre of GN20, either by an AGN or an intense starburst, could have deposited energy in the interstellar medium, thus increasing the intrinsic velocity dispersion in the central parts of the galaxy (see e.g. Harrison et al. 2016; Übler et al. 2019; Marasco et al. 2023; Parlanti et al. 2024). We do not expect any further impact on the narrow-component kinematics from outflows.  \nFinally, we consider non-circular motions in the form of radial flows. Those could be triggered through smooth accretion from the circumgalactic medium or internal streaming motions caused by disc instabilities or substructures such as spiral arms or a bar. The stellar light distribution derived through the MIRI F560W data by Colina et al. (2023) does not provide evidence for a bar, though some substructure is revealed through their residual analysis which may indicate spiral-arm or ring-like features in the outer parts of the galaxy. Based on our assumption of clockwise rotation (i.e. the bright NW side being the 'near side', and the fainter SE side being the 'far side'), the signature of an axisymmetric, planar radial inflow would be a mirrored S-shaped twist in the iso-velocity contours (see van der Kruit & Allen 1978, and the recent examples at 𝑧 ∼ 2 by Genzel et al. 2023 and Price et al. 2021). We observe a similar effect in the velocity field of GN20 (see bottom middle panel of Figure 3). This motivates us to explore a second dynamical model including a uniform, planar radial inflow. This is clearly a simplified assumption for GN20, but it may still provide us with informative clues about the nature of the observed deviations from circular motions. We show the results from this second model in Figure A1. This model, which includes an inflow with 𝑣 𝑟 = 130 km/s, indeed provides a somewhat better fit to the observed kinematics ( Δ 𝜒 2 red = 0 . 33), suggesting that gas may be flowing from the outskirts to the centre of GN20. However, the strong residuals in the NW largely remain, indicating that additional processes for instance recent or ongoing interaction disturb the kinematics in GN20.", '6 BROAD NUCLEAR EMISSION AND AGN SIGNATURES': "As discussed in Sections 3, we find broad emission in the H 𝛼 +[N II] complex in the central region of GN20. To derive the narrow emission line maps we have used for the dynamical modelling, we have fitted this broad emission as a set of broad Gaussians for H 𝛼 and [N II] 𝜆𝜆 6548 , 6583 in addition to a set of narrow components, i.e. interpreting this broad emission as an outflow component. However, the data are very noisy due to the high obscuration in the centre. In fact, if we alternatively fit the central regions with only one broad component for H 𝛼 in addition to the narrow line components for H 𝛼 +[N II] 𝜆𝜆 6548 , 6583, we get spectral fits of similar quality. We, therefore, cannot exclude, given the quality of the data, that the broad emission in the centre traces the broad-line region (BLR) of an accreting black hole. Furthermore, more complex scenarios including both a potential BLR and an outflow component may be possible. In the following we briefly discuss the two 'extreme' scenarios of a pure outflow and a pure BLR, however we stress that we cannot robustly distinguish between these two scenarios based on our data.", '6.1 The outflow scenario': 'Starting with the outflow interpretation, we show examples of our two-component fits to the Voronoi bin in the central region showing the brightest emission (top), and to the integrated spectrum extracted over a larger aperture where broad emission is visible (bottom) in the left panels of Figure 5. In both regions (and indeed in all individual regions encompassed by the larger aperture) the broad components are blue-shifted. Their FWHM of 820-1340 km/s are typical of AGNdriven winds (e.g. Veilleux et al. 2005, 2020; Fabian 2012; Heckman &Best 2014; Genzel et al. 2014; Harrison et al. 2016; Carniani et al. 2015; Rupke et al. 2017; Förster Schreiber et al. 2019, and references therein). From the large integrated aperture, we calculate a maximum outflow velocity of 𝑣 out = ⟨ 𝑣 broad ⟩ + 2 𝜎 broad = 970 ± 220 km/s, where 𝜎 broad is corrected for instrumental resolution (e.g. Genzel et al. 2011; Davies et al. 2019). The large outflow velocities indicate that the outflow is driven by an AGN (see also discussion by Maiolino et al. 2024).  \nAssuming a photo-ionised, constant-velocity spherical outflow (Genzel et al. 2011; Newman et al. 2012; Förster Schreiber et al. 2019; Davies et al. 2019, 2020; Cresci et al. 2023) of the extent 𝑅 out equal to the aperture radius ( ∼ 2 kpc), with an electron density of 𝑛 𝑒, out = 1000/cm 3 (e.g. Perna et al. 2017; Kakkad et al. 2018; Förster Schreiber et al. 2019), we would find a low mass-outflow rate of / 𝑀 out , ion = 0 . 7 ± 0 . 3 𝑀 ⊙ /yr. We note that this estimate is uncertain due to the unknown outflow geometry and electron density in the outflow (the data in the central region is too noisy for a direct measurement from the [S II] doublet). Certainly, this value would correspond to a lower limit, considering the high obscuration in the centre of GN20, and because we are only tracing the warmionisedgasphase(e.g. Rupke & Veilleux 2013; Herrera-Camus et al. 2019; Roberts-Borsani 2020; Fluetsch et al. 2021; Avery et al. 2022; Baron et al. 2022; Cresci et al. 2023; Belli et al. 2024; Davies et al. 2024). For comparison, in a recent study of another obscured AGN at 𝑧 = 4 . 76, Parlanti et al. (2024) find a mass outflow rate of / 𝑀 out , ion = 11 + 57 -5 . 5 𝑀 ⊙ /yr.', '6.2 The BLR scenario': 'In the right panels of Figure 5, we show the corresponding fits including only one broad component fixed at the position of the narrow H 𝛼 emission in the central Voronoi bin (see upper panel). The  \nFigure 5. Two-component fits to integrated spectra in the central region of GN20. The top panels show fits to the central Voronoi bin with the brightest broad emission flux, while the bottom panel shows a fit to the larger aperture encompassing regions with broad emission based on visual inspection (see black contours in insets). Left: fits including a set of broad components for H 𝛼 +[N II] 𝜆𝜆 , i.e. interpreting the broad emission as an outflow. The broad emission is blue-shifted with respect to the narrow emission, with FWHM typical of AGN-driven winds, particularly in the central Voronoi bin (FWHM broad > 1000 km/s). Right: fits including a broad component only for H 𝛼 , i.e. interpreting the broad emission as the broad-line region of an accreting black hole. The two fitting setups have a comparable goodness-of-fit. This is also apparent from the very similar residuals shown in the bottom part of the panels, res.=(data-model)/uncertainties. The integrated spectra are not corrected for the velocity field.  \n<!-- image -->  \nfit to the larger aperture is shown in the bottom, where we keep the position and FWHM of the broad component fixed to the fitting results from the central Voronoi bin. The fits are of comparable quality to the fits interpreting the broad emission as an outflow ( Δ BIC outflow -BLR , central = 14, Δ BIC outflow -BLR , large aperture = 9). From the properties of the (unobscured) putative BLR component, we would derive a black hole mass of log ( 𝑀 · / 𝑀 ⊙ ) = 7 . 3 ± 0 . 4 following Reines & Volonteri (2015) and using the fit to the large aperture. For a galaxy of the stellar mass of GN20 ( 𝑀 ★ ∼ 1 . 1 × 10 11 𝑀 ⊙ ), a value of log ( 𝑀 · / 𝑀 ⊙ ) ∼ 7 . 3 would fall well within the scatter of local BLR AGN. However, this estimate does not account for extinction towards the BLR, and may therefore correspond to a lower limit, log ( 𝑀 · / 𝑀 ⊙ ) ≳ 7 . 3.  \nAssuming that GN20 hosts an obscured AGN, Riechers et al. (2014) estimated an Eddington limit for its black hole mass of log ( 𝑀 Edd · / 𝑀 ⊙ ) = 8 . 1 -8 . 5 based on its 6 𝜇 m continuum luminosity and an upper limit on the 2 -10 keV 𝑋 -ray luminosity. This value is consistent with our lower limit.', '6.3 Emission line diagnostics': 'In addition to the broad linewidths in the central region of GN20, suggesting either large outflow velocities or the presence of an actively accreting black hole, we find high narrow-component ratios of [N II] 𝜆 6583/H 𝛼 > 0 . 6 throughout large parts of GN20, and especially in the Northern nuclear region and toward the South-West. This  \nFigure 6. Left: map of the narrow-component [N II] 𝜆 6583/H 𝛼 ratio derived from fits to the H 𝛼 Voronoi-binned map (see Section 3). High ratios of [N II] 𝜆 6583/H 𝛼 > 0 . 6 throughout large parts of GN20, and especially in the Northern central regions and towards the South-West, indicate likely photo-ionisation by an AGN. Right: GN20 as classified through the WHAN diagnostic diagram (see Section 6.3). GN20 is largely consistent with the Seyfert regime, while some regions in the outskirts are consistent with being photo-ionised by star formation (SF). The black cross indicates the galaxy centre as adopted for the dynamical modelling.  \n<!-- image -->  \nis illustrated in the left panel of Figure 6, where we show a map of the narrow-component [N II] 𝜆 6583/H 𝛼 ratio based on the Voronoibinned map (see Section 3). Values of [N II] 𝜆 6583/H 𝛼 > 0 . 6 may be indicative of contributions through shocks, 𝛼 -enhanced evolved stellar populations, or photo-ionisation by an AGN (Baldwin et al. 1981; Kewley et al. 2001; Kauffmann et al. 2003; Byler et al. 2019).  \nThe  \nBPT  \ndiagram  \n(Baldwin  \net  \nal.  \n1981)  \ncomparing  \n[N II] 𝜆 6583/H 𝛼 to [O III] 𝜆 5007/H 𝛽 can provide insights into the dominant ionisation mechanism in galaxies. The [O III] 𝜆 5007 and H 𝛽 lines are not covered in our grating observations. They are covered in the prism observations, yet due to the high obscuration they are undetected in most regions (see also Maseda et al. 2024). As an alternative, we use the WHAN diagnostic diagram, which utilises the equivalent width of H 𝛼 in combination with the [N II] 𝜆 6583/H 𝛼 ratio (Cid Fernandes et al. 2010, 2011). In this diagram, for EW(H 𝛼 ) > 6 Å and log([N II] 𝜆 6583/H 𝛼 ) > -0 . 4, galaxies fall into the Seyfert regime. Indeed, we find that most regions in GN20 are consistent with the Seyfert regime, while some regions in the outskirts are classified as star-forming (SF; EW(H 𝛼 ) > 3 Å and log([N II] 𝜆 6583/H 𝛼 ) < -0 . 4). We show the classification based on the WHAN diagram throughout GN20 in the right panel of Figure 6.', '7 SUMMARY AND DISCUSSION': "Wehave presented an analysis of the rest-frame optical kinematics in the large 𝑧 ∼ 4 . 055 sub-mm galaxy GN20 using JWST /NIRSpec-IFU data. In the following, we discuss and summarise our main findings.  \nH 𝛼 morphology. We find extended, but clumpy and asymmetric H 𝛼 emission, that is largely arranged in a loop around a central emission peak. H 𝛼 is brightest along the loop in the North-Western region (see Figure 1). The central H 𝛼 emission peak broadly coincides with the compact rest-frame 1 . 1 𝜇 m nuclear (stellar) emission found through MIRI/F560W imaging by Colina et al. (2023) (see also Crespo Gómez et al. 2024). The edge of the extended envelope identified in their analysis overlaps with the H 𝛼 loop. The brightest H 𝛼 emission region partly overlaps with the rest-frame UV emission obtained with HST in F105W, but is generally located further to the North-East. We compare the rest-frame optical line emission obtained through our NIRSpec-IFS observations with rest-frame UV and near-infrared wavebands in Figure 2.  \nLarge-scale H 𝛼 kinematics. We find smooth, large-scale rotation in H 𝛼 ( Δ 𝑣 obs ∼ 610 km/s; 𝜎 obs , disc ∼ 50 -150 km/s; Figure 3). We construct a three-dimensional dynamical model including a thick disc, bulge, and dark matter halo, and simultaneously model the observed velocity and velocity dispersion maps by minimization on the projected two-dimensional kinematics. From our best-fit model (Figure 4), we find 𝑣 rot ( 𝑅 𝑒 ) = 469 km/s, 𝜎 0 = 89 km/s, 𝑣 rot ( 𝑅 𝑒 )/ 𝜎 0 = 5 . 3, and log ( 𝑀 dyn ( < 2 𝑅 𝑒 )/ 𝑀 ⊙ ) ∼ 11 . 7. These results are largely consistent with previous analyses of the molecular and ionised gas kinematics constrained by VLA observations of CO(2-1), PdBI observations of CO(6-5), and MIRI/MRS observations of Pa 𝛼 (Carilli et al. 2010; Hodge et al. 2012; Bik et al. 2024).  \nNon-circular motions. The H 𝛼 velocity field shows deviations from circular motions. These deviations are particularly apparent in the North-Western region of GN20, as also confirmed through a residual analysis based on our best-fit dynamical model (Figure 4). Non-circular motions may arise from a variety of phenomena, and several may be acting simultaneously in GN20, potentially including tidal interactions, accretion and radial motions. To test one scenario, we construct a second dynamical model including a planar radial inflow of 𝑣 𝑟 = 130 km/s, which indeed provides a slightly better fit to the data ( Δ 𝜒 2 red = 0 . 33). However, significant residuals in the North-West remain also with this model, indicating that additional causes for perturbation of simple rotation exist in GN20 (Figure A1).  \nAGN signatures. We find broad (FWHM ∼ 1000 -2000 km/s) emission in the H 𝛼 +[N II] complex in the central region of GN20. This broad emission can be modelled as a blue-shifted, high-velocity  \noutflow in H 𝛼 and [N II] ( 𝑣 out = 970 ± 220 km/s), or as a broadline region in H 𝛼 emission, corresponding to a black hole of mass log ( 𝑀 · / 𝑀 ⊙ ) ≳ 7 . 3 (Figure 5). We cannot robustly distinguish between these scenarios (or a combination of both) based on the data quality, however, both scenarios provide evidence for the presence of an active black hole in GN20. Further evidence is provided through high values of narrow-component [N II] 𝜆 6583/H 𝛼 > 0 . 4 together with EW(H 𝛼 ) > 6 Å throughout large regions of the galaxy, consistent with photo-ionisation by an AGN (Figure 6).  \nEvidence for black hole feeding and feedback. Our kinematic and dynamical analysis provides tentative evidence that gas may be channeled into the nuclear region of GN20. These inflows could potentially fuel central star formation and/or accretion onto the central black hole. Indeed we find evidence for an AGN in the centre of GN20 in the form of broad nuclear emission, signaling either AGNdriven outflows or an actively accreting black hole, and Seyfert-like ionisation throughout the galaxy. In the local Universe, evidence for AGN feeding through molecular gas inflows has been observed, for instance through gravity torques induced by spiral or bar structure (e.g. Maiolino et al. 2000; Hunt et al. 2008; Casasola et al. 2011; Combes et al. 2014, 2019; Speights & Rooke 2016; Venturi et al. 2018; Audibert et al. 2019; Izumi et al. 2023). At higher redshift, 1 < 𝑧 < 2 . 5, signatures of both ionised and molecular gas inflows have been revealed in deep observations of massive main-sequence galaxies (Price et al. 2021; Genzel et al. 2023).  \nIn summary, our analysis suggests that GN20, an intriguingly large and massive SMG 1.5 billion years after the Big Bang, shows evidence for both non-circular motions (potentially fuelling AGN activity), and high ionisation and possibly high-velocity outflows (i.e. feedback by the AGN), or the broad-line region of an actively accreting black hole. This work together with the recent studies by, for instance, Arribas et al. (2023); Jones et al. (2024); Parlanti et al. (2024) highlights the potential of JWST /NIRSpec-IFS to unveil detailed, spatially-resolved properties of 𝑧 > 4 dusty star-forming galaxies, such as the source(s) powering nebular emission and high-accuracy kinematic measurements, even in cases of extremely high dust attenuation. While SMGs were previously mainly characterised through their bright sub-mm emission, JWST opened up new avenues to study their formation and evolution. As argued by several authors, the rapid evolution, high ongoing star formation, and dense environment make GN20 a prime progenitor candidate for present-day massive earlytype galaxies, those systems that contain the majority of the stellar mass existing today. Deeper NIRSpec observations of GN20 could help to further uncover the source of its broad nuclear emission, and to investigate in greater detail the deviations from circular motions in the H 𝛼 kinematics. This might shed further light on the assembly and evolution history of today's most massive galaxies.", 'ACKNOWLEDGEMENTS': "We are grateful to the anonymous referee for a constructive report that helped to improve the quality of this manuscript. We thank Chris Carilli and Luis Colina, Arjan Bik and the MIRI GTO team for sharing data. We thank Arjan Bik and Takafumi Tsukui for comments on an earlier version of this manuscript. HÜ thanks the Infrared/Submillimeter Group at MPE for making the latest version of DysmalPy available prior to the code becoming public. Special thanks go to Stavros Pastras and Juan Manuel Espejo Salcedo for assistance with the code installation. AJB and GCJ acknowledge funding from the 'FirstGalaxies' Advanced Grant from the European Research Council (ERC) under the European Union's Horizon 2020  \nresearch and innovation program (Grant agreement No. 789056). BRdP, MP and SA acknowledge grant PID2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research (MICIN/AEI/ 10.13039/501100011033). FDE, JS and RM acknowledge support by the Science and Technology Facilities Council (STFC), by the ERC through Advanced Grant 695671 'QUENCH', and by the UKRI Frontier Research grant RISEandFALL. GV and SCa acknowledge support from the European Union (ERC, WINGS,101040227). HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton-Kavli Junior Fellowship. IL acknowledges support from PID2022-140483NB-C22 funded by AEI 10.13039/501100011033 and BDC 20221289 funded by MCIN by the Recovery, Transformation and Resilience Plan from the Spanish State, and by NextGenerationEU from the European Union through the Recovery and Resilience Facility.", 'DATA AVAILABILITY': 'The NIRSpec data used in this research has been obtained within the NIRSpec-IFU GTO programme GA-NIFS and is publicly available at MAST. Data presented in this paper will be shared upon reasonable request to the corresponding author.', 'REFERENCES': 'Data Analysis Software and Systems XVIII. p. 251 ( arXiv:0902.2850 ) Maseda M. V., et al., 2024, arXiv e-prints, p. arXiv:2403.05506  \nMorrison G. E., Owen F. N., Dickinson M., Ivison R. J., Ibar E., 2010, ApJS, 188, 178  \nNavarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563  \nNewman S. F., et al., 2012, ApJ, 752, 111  \nNoordermeer E., 2008, MNRAS, 385, 1359  \nOteo I., et al., 2018, ApJ, 856, 72  \nOverzier R. A., 2016, A&ARv, 24, 14  \nParlanti E., et al., 2024, A&A, 684, A24  \nPavesi R., et al., 2018, ApJ, 861, 43  \nPerna M., Lanzuisi G., Brusa M., Cresci G., Mignoli M., 2017, A&A, 606,  \nA96  \nPerna M., et al., 2023, A&A, 679, A89  \nPlanck Collaboration et al., 2020, A&A, 641, A6  \nPope A., Borys C., Scott D., Conselice C., Dickinson M., Mobasher B., 2005,  \nMNRAS, 358, 149  \nPope A., et al., 2006, MNRAS, 370, 1185  \nPrice S. H., et al., 2021, ApJ, 922, 143  \nReines A. E., Volonteri M., 2015, ApJ, 813, 82  \nRenzini A., 2006, ARA&A, 44, 141  \nRiechers D. A., et al., 2014, ApJ, 786, 31  \nRieke G. H., et al., 2015, PASP, 127, 665  \nRoberts-Borsani G. W., 2020, MNRAS, 494, 4266  \nRoberts W. W. J., Huntley J. M., van Albada G. D., 1979, ApJ, 233, 67  \nRupke D. S. N., Veilleux S., 2013, ApJ, 768, 75  \nRupke D. S. N., Gültekin K., Veilleux S., 2017, ApJ, 850, 40  \nShlosman I., Frank J., Begelman M. C., 1989, Nature, 338, 45  \nSimpson J. M., et al., 2014, ApJ, 788, 125  \nSmail I., Ivison R. J., Blain A. W., 1997, ApJ, 490, L5  \nSpeights J. C., Rooke P. C., 2016, ApJ, 826, 2  \nSwinbank A. M., et al., 2014, MNRAS, 438, 1267  \nTan Q., et al., 2014, A&A, 569, A98  \nToomre A., 1964, ApJ, 139, 1217  \nTsukui T., Wisnioski E., Bland-Hawthorn J., Mai Y., Iguchi S., Baba J.,  \nFreeman K., 2024, MNRAS, 527, 8941  \nÜbler H., et al., 2018, ApJ, 854, L24  \nÜbler H., et al., 2019, ApJ, 880, 48  \nVeilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 43, 769  \nVeilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, A&ARv, 28, 2  \nVenturi G., et al., 2018, A&A, 619, A74  \nVerley S., et al., 2007, A&A, 472, 121  \nWada K., Habe A., 1992, MNRAS, 258, 82  \nWalter F., et al., 2012, Nature, 486, 233  \nWright G. S., et al., 2015, PASP, 127, 595  \nWuyts S., et al., 2012, ApJ, 753, 114  \nWuyts S., et al., 2016, ApJ, 831, 149  \nvan Albada G. D., Roberts W. W. J., 1981, ApJ, 246, 740  \nvan Dokkum P. G., 2001, PASP, 113, 1420  \nvan der Kruit P. C., Allen R. J., 1978, Annual Review of Astronomy and Astrophysics, 16, 103', 'APPENDIX A: DYNAMICAL MODEL INCLUDING A FIXED RADIAL INFLOW': 'In Figure A1 we show the results of our dynamical modelling including a uniform, planar radial inflow of fixed velocity 𝑣 𝑟 = 130 km/s. I.e., we assume that gas is inflowing through the disc in an axisymmetric fashion and with radially constant velocity. The inflow velocity is not fitted for, but added as a fixed parameter (see Price et al. 2021). The preferred inflow velocity is initially identified through a grid search by injecting different values of 𝑣 𝑟 in steps of 10 km/s from zero to 200 km/s, using least-squares minimization. We then obtain the best-fit through MCMC techniques as described in Section 3, with the same free parameters as our fiducial model plus the radial  \nTable A1. Results from our best-fit dynamical model including a thick disc, bulge and dark matter halo, and a planar radial inflow of velocity 𝑣 𝑟 = 130 km/s. For 𝜎 0 , we adopt minimum uncertainties of 10 km/s, and for the total baryonic mass we adopt minimum uncertainties of 0.1 dex. For 𝑓 DM we give the 3 𝜎 upper limit.  \ninflow fixed at 𝑣 𝑟 = 130 km/s. Our modelling results are listed in Table A1.  \nWhile we still find strong residuals with this second model specifically in the North-Western region of GN20, the overall magnitude of the velocity residuals decreases. The reduced chi-squared statistics are improved compared to our fiducial model presented in Section 5 by Δ 𝜒 2 red = 0 . 33. This shows that the inclusion of a radial inflow can help to explain some of the non-circular motions observed in GN20.  \nWe point out the change in the kinematic position angle of the best-fit model including inflow with respect to our fiducial model. The best-fit rotation velocities are reduced. In this model, we find little contribution of dark matter to the dynamics on galactic scales.  \nWe compare the inflow velocity identified through our grid search to analytical estimates, as recently presented by Genzel et al. (2023). We follow their Equation 11 and parameter choices to obtain an analytical expectation for the inflow velocity. For this, we evaluate the Toomre -𝑄 parameter (Toomre 1964) following Binney & Tremaine (2008); Escala & Larson (2008); Dekel et al. (2009) at the radius where our best-fit circular velocity curve reaches its maximum (see Übler et al. 2019), and find 𝑄 gas ∼ 0 . 49. Following Genzel et al. (2023), we then obtain 𝑣 𝑟 = 102 km/s, which is somewhat lower than our model-derived value. If we approximate the inflow velocity through 𝑣 𝑟 ∼ 𝑓 2 gas 𝑣 𝑐 (Genzel et al. 2023), with 𝑣 𝑐 measured at the effective radius, we obtain 𝑣 𝑟 = 138 km/s, comparable to our modelderived value. In general, as pointed out by Genzel et al. (2023), analytical estimates of inflow velocities can vary by factors of a few (e.g. Dekel et al. 2009, 2013; Krumholz & Burkert 2010; Krumholz et al. 2018). Our comparison to analytical estimates shows that the model-derived inflow velocity is not unexpected for a galaxy with the mass and kinematics of GN20.  \nEvidence for inflows in massive, high-redshift ( 𝑧 ∼ 1 -2 . 5) starforming disc galaxies has been found in deep ground-based observations (Price et al. 2021; Genzel et al. 2023). The radial velocities in those cases are measured to be 𝑣 𝑟 ∼ 30 -120 km/s, comparable to the putative inflow velocity of 𝑣 𝑟 = 130 km/s in GN20.  \nIn Figure A2, we show the difference between our fiducial velocity model as presented in the top middle-left panel of Figure 3, and the velocity model including the radial inflow (top middle-left panel of Figure A1). The differences between the two models are most pronounced in the centre and in the outer regions.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure A1. Same as Figure 4, but including a planar radial inflow of 𝑣 𝑟 = 130 km/s in addition to our fiducial disc, bulge and halo model. Top (bottom) panels show the observed velocity (velocity dispersion) field (left), the best-fit model (middle left), residuals (data-model; middle right), and goodness-of-fit (data-model)/uncertainties (right). The cross and line indicate the centre and best-fit PA. The median velocity offset is Δ 𝑣 med = 8 . 8 km/s with a rms velocity difference of Δ 𝑣 rms = 40 . 0 km/s. Corresponding values for the velocity dispersion are Δ 𝜎 med = 15 . 4 km/s and Δ 𝜎 rms = 34 . 6 km/s. While we observe strong residuals particularly in the North-Western region, as in our fiducial model, the magnitude of the residuals in the velocity field has decreased in this model, indicating that part of the deviations from circular motions could be explained by an inflow component.  \n<!-- image -->  \nFigure A2. Difference in the model velocities of our fiducial model compared to the model including the radial inflow of 𝑣 𝑟 = 130 km/s. The differences are most pronounced in the centre and in the outer regions.  \n<!-- image -->'}