# GA-NIFS: NIRSpec reveals evidence for non-circular motions and AGN feedback in GN20

Hannah Übler, Francesco d'Eugenio, Michele Perna, Santiago Arribas, Gareth

C. Jones, Andrew J. Bunker, Stefano Carniani, Stéphane Charlot, Roberto

Maiolino, Bruno Rodríguez del Pino, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Hannah Übler, Francesco d'Eugenio, Michele Perna, Santiago Arribas, Gareth C. Jones, et al., GA-NIFS: NIRSpec reveals evidence for non-circular motions and AGN feedback in GN20. Monthly Notices of the Royal Astronomical Society, 2024, 533, pp.4287-4299. 10.1093/mnras/stae1993. insu-04851506

HAL Id: insu-04851506

https://insu.hal.science/insu-04851506v1

Submitted on 21 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# GA-NIFS: NIRSpec reveals evidence for non-circular motions and AGN feedback in GN20

Hannah Übler<sup>1,2*</sup> Francesco D'Eugenio,<sup>1,2</sup> Michele Perna,<sup>3</sup> Santiago Arribas,<sup>3</sup> Gareth C. Jones,<sup>4</sup>

Andrew J. Bunker,<sup>4</sup> Stefano Carniani<sup>5</sup>,<sup>5</sup> Stéphane Charlot<sup>6</sup>,<sup>6</sup> Roberto Maiolino,<sup>1,2,7</sup>

Bruno Rodríguez del Pino<sup>1,3</sup>, Chris J. Willott,<sup>8</sup> Torsten Böker<sup>1,9</sup>,<sup>9</sup> Giovanni Cresci,<sup>10</sup>

Nimisha Kumari<sup>1,11</sup>, Isabella Lamperti<sup>1,3</sup>,<sup>3</sup> Eleonora Parlanti,<sup>5</sup> Jan Scholtz<sup>1,2</sup> and Giacomo Venturi<sup>1,5</sup>

Accepted 2024 August 15. Received 2024 July 21; in original form 2024 March 6

## ABSTRACT

We present rest-frame optical data of the $z \sim 4$ submillimetre galaxy GN20 obtained with the JWST Near Infrared Spectrograph (NIRSpec) in integral field spectroscopy mode. The $H\alpha$ emission is asymmetric and clumpy and extends over a projected distance of $>15$ kpc. To first order, the large-scale ionized gas kinematics are consistent with a turbulent ($\sigma \sim 90$ km s$^{-1}$), rotating disc ($v_{rot} \sim 500$ km s$^{-1}$), congruent with previous studies of its molecular and ionized gas kinematics. However, we also find clear evidence for non-circular motions in the $H\alpha$ 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 (full width at half-maximum $\sim 1000-2000$ km s$^{-1}$) in the $H\alpha + [N\text{II}]$ complex, suggestive of fast, active galactic nucleus-driven winds or, alternatively, of the broad-line region of an active black hole. Elevated values of $[N\text{II}] \lambda 6583/H\alpha > 0.4$ and of the $H\alpha$ equivalent width $EW(H\alpha) > 6$ Å 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 present-day galaxies acting in concert, over 12 billion years ago.

**Key words:** galaxies: active – galaxies: high-redshift – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

More than half of the total stellar mass in the present-day Universe is situated in early-type galaxies (ETGs; Fukugita, Hogan & Peebles 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 $z > 3$, the discovery of bright submillimetre galaxies (SMGs; see review by Blain et al. 2002) at $z \gtrsim 4$ quickly promoted them to candidate progenitors (e.g. Simpson et al. 2014). These galaxies rapidly assemble high stellar masses through intense star formation (SF; e.g. Hainline et al. 2009; Magnelli et al. 2013; Casey, Narayanan & Cooray 2014; Swinbank et al. 2014). Besides, some

SMGs, including the galaxy studied in this work, are located within cosmic overdensities, hinting at future merger events which would help them to further increase their mass by $z=0$ (e.g. Daddi 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. 2024; Jones et al. 2024). Furthermore, at redshifts $z > 3$, the cores of protoclusters, 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.

After the identification of many $z > 4$ SMGs in the past decades (e.g. Smail, Ivison & Blain 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 JWST has enabled, for the first time, studies of their rest-frame optical emission-line properties. Within the Galaxy Assembly with NIRSpec Integral Field Spectroscopy (GA-NIFS)

<footnote n="1"><sup>1</sup>Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="2"><sup>2</sup>Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HA, UK</footnote>

<footnote n="3"><sup>3</sup>Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, Torrejón de Ardoz, E-28850 Madrid, Spain</footnote>

<footnote n="4"><sup>4</sup>Department of Physics, Denys Wilkinson Building, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="5"><sup>5</sup>Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy</footnote>

<footnote n="6"><sup>6</sup>Institut d' Astrophysique de Paris, UMR 7095, Sorbonne Université, CNRS, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="8"><sup>8</sup>National Research Council of Canada, Herzberg Astronomy & Astrophysics Research Centre, 5071 West Saanich Road, Victoria BC V9E 2E7, Canada</footnote>

<footnote n="9"><sup>9</sup>European Space Agency, c/o STScI, 3700 San Martin Drive, Baltimore MD 21218, USA</footnote>

<footnote n="10"><sup>10</sup>INAF - Osservatorio AstrofSCO di Arcetri, largo E. Fermi 5, I-50127 Firenze, Italy</footnote>

<footnote n="11"><sup>11</sup>AURA for the European Space Agency, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote>*E-mail: hu215@cam.ac.uk</footnote>
survey<fnref n="1" /> (PIs: Santiago Arribas and Roberto Maiolino), as part of the NIRSpec Instrument Science Team Guaranteed Time Observations, we have targeted several bright SMGs at $z \ge 4$ with the integral field spectroscopy (IFS) mode (Böker et al. 2022; Jakobsen et al. 2022), three of which have already been discussed: ALESS073.1 at $z = 4.76$, for which the NIRSpec-IFS data reveal a heavily dust-obscured active galactic nucleus (AGN; Parlanti et al. 2024); HLFS3 at $z = 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 $z = 6.9$, in which Arribas et al. (2024) find evidence for accretion from the cosmic web, inflows, and mergers. In this paper, we present high-resolution ($R \sim 2700$) NIRSpec-IFS data of the rest-frame optical line-emission in GN20, a large SMG at $z \sim 4$.

GN20 was detected by Pope et al. (2005) as a bright 850 µm source in the Great Observatories Origins Deep Survey Northern field (see also Iono et al. 2006; Pope et al. 2006). This dusty star-forming (SF) galaxy is located within a proto-cluster environment, and is detected in 1.4 GHz continuum, [C]1 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 $\sim 2–3 \times 10^{13}$ L$_{\odot}$, a star-formation rate (SFR) of $SFR_{IR} \sim 1800–3000$ M$_{\odot}$ yr$^{-1}$, a stellar mass of $M_{\star} \sim 1.1–2.3 \times 10^{11}$ M$_{\odot}$, and a molecular gas mass of $M_{mol} \sim 5–13 \times 10^{10}$ M$_{\odot}$ (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$^{-1}$ is observed in CO(4–3), consistent with an exceptionally large $r \sim 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 $v_{rot,max} = 575 \pm 100$ km s$^{-1}$, with a velocity dispersion of $\sigma = 100 \pm 30$ km s$^{-1}$. Dynamical mass estimates for GN20 are in the range of $2–6 \times 10^{11}$ M$_{\odot}$ (Daddi et al. 2009; Carilli et al. 2010, 2011; Hodge et al. 2012). Intriguingly, rest-frame ultraviolet (UV) emission tracing young stars as observed with Hubble Space Telescope (HST)/WFC3 F105W imaging is visible only in an extended ($\sim 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).

Recently, Colina et al. (2023) studied the rest-frame 1.1 µm imaging of GN20 obtained with the MIRI instrument onboard JWST (see also Crespo Gómez et al. 2024). Their analysis reveals a two-component stellar structure composed of an unresolved nucleus offset by 1 kpc from the centre of an extended disc with $R_e = 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 SF, while the extended stellar envelope overlaps with the cold molecular gas distribution. Analysing integral-field unit (IFU) observations taken with the medium-resolution spectrometer (MRS) of the Mid-Infrared Instrument (MIRI) on board JWST, Bik et al. (2024) find clumpy Pa$\alpha$ emission out to a radius of 6 kpc. The kinematics of the Pa$\alpha$ emission are consistent with a rotating disc, with a maximum rotation velocity of $v_{rot} = 550 \pm 40$ km s$^{-1}$ and an upper limit on the flux-weighted velocity dispersion ($\sigma_m$) of $\sigma_m = 145 \pm 53$ km s$^{-1}$. Comparing the unobscured SFR derived from the integrated Pa$\alpha$ flux, $SFR_{Pa\alpha} = 144 \pm 9$ M$_{\odot}$ yr$^{-1}$, with the infrared SFR i.e. $SFR_{IR}$ (Tan et al. 2014), Bik et al. (2024) infer a high average extinction of $A_V = 17.2 \pm 0.4$ mag (see also Maseda et al. 2024). Crespo

Gómez et al. (2024) infer a lower $A_V \sim 1.5$ mag from spectral energy distribution fitting. They attribute this difference to the presence of either stellar populations older than 10 Myr, or of a buried AGN.

We describe the JWST/NIRSpec-IFS high-resolution observations of GN20 and our analysis methods in Sections 2 and 3. We present the complex H$\alpha$ emission revealed by NIRSpec in Section 4. In Section 5, we discuss the H$\alpha$ kinematics and evidence for non-circular motions based 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 summarize our findings in Section 7.

Throughout this work, we assume a flat Lambda cold dark matter cosmology with $\Omega_m = 0.315$ and $H_0 = 67.4$ km s$^{-1}$ Mpc$^{-1}$ (Planck Collaboration VI 2020). With this cosmology, 1 arcsec corresponds to a transverse distance of 7.07 proper kpc at $z = 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 2023 February 10, 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 of 2.87–5.14 µm (spectral resolution $R \sim 2000–3500$; Jakobsen et al. 2022), and about 1 h with PRISM/CLEAR ($\lambda = 0.6–5.3$ µm, spectral resolution $R \sim 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 were presented by Colina et al. (2023) and Crespo Gómez et al. (2024), and the MRS data were presented by Bik et al. (2024).

Raw data files were downloaded from the Barbara A. Mikulski Archive for Space Telescopes and subsequently processed with the JWST Science Calibration pipeline<fnref n="2" /> version 1.11.1 under the Calibration Reference Data System 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 summarize here. Count-rate frames were corrected for $1/f$ noise through a polynomial fit. During calibration in Stage 2, we removed regions affected by failed open MSA (micro-shutter assembly) 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 LACOSMIC (van Dokkum 2001): We calculated the derivative of the count-rate maps along the dispersion direction, normalized it by the local flux (or by three times the rms noise, whichever was highest), and rejected the 95th 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 arcsec. We used spaxels away from the central source and free of emission features to perform a background subtraction.

In Fig. 1, we show a line map of the emission integrated in the range 3.313–3.327 µm, roughly encompassing the H$\alpha$ and [N II] λ6583 lines at $z \sim 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 Fig. 2 we use the prism observations to create a false-colour multiwavelength image of GN20, combining the H$\alpha$ + [N II] λλ6548, 6583 emission (green) with the emission at 0.9–1.2 µm (blue), roughly corresponding to the

<footnote n="1">¹https://ga-nifs.github.io</footnote>

<footnote n="2">²https://jwst-pipeline.readthedocs.io/en/stable/jwst/introduction.html</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Line map of the emission integrated in the range of 3.313–3.327 µm, roughly encompassing the $H\alpha + [N\text{II}]\lambda6583$ lines at $z \sim 4.055$. GN20 is visible as the extended central source with a complex and clumpy morphology. We indicate two clumps to the north-west and extended emission in the south-east 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$^{-1}$, respectively, a redshift of $z = 4.06\,532$. The white circle indicates the approximate PSF full width at half-maximum (FWHM; $\sigma = 0.12$ arcsec).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Comparison of rest-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\alpha + [N\text{II}]$ emission is shown in green, blue colours 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 colours 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\alpha + [N\text{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\alpha + [N\text{II}]$ loop. We note that the $Pa\alpha$ 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\alpha + [N\text{II}]$, including the location of the north-western clumps.</figcaption></figure>

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 Fig. 1 is derived by registering our observations to the F560W and F105W images, which are in turn registered to Gaia DR3 (Gaia Collaboration 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 arcsec, we estimate an uncertainty of about 0.1 arcsec on the astrometry.

# 3 FITTING AND DYNAMICAL MODELLING

## 3.1 Emission line fitting and $H\alpha$ redshift

To analyse the emission-line properties and kinematics of GN20, we first fit a one-component Gaussian model to the $H\alpha$ and [NII] $\lambda\lambda6548, 6583$ emission lines in the grating data, including a constant continuum. Because the signal-to-noise ratio (S/N) 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 VORBIN (Cappellari & Copin 2003), as implemented in QFITSVIEW. We first define a mask based on visual inspection of the line emission in the cube, and then require $S/N \ge 15$ in the previously fitted $H\alpha$ flux map to derive the Voronoi map. We derive $H\alpha$ 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 Fig. 3.

While the line emission is relatively narrow in the outer regions of GN20 ($\sigma_{\text{obs}} \sim 50–100$ km s$^{-1}$), the $H\alpha + [N\text{II}]$ 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 Fig. 5, diameter of $\sim 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\alpha$ and [NII] $\lambda\lambda6548, 6583$. Here, we require $\text{FWHM}_{\text{broad}} > 400$ km s$^{-1}$ $> \text{FWHM}_{\text{narrow}}$. We stitch the narrow component results to our initial one-component maps. As can be seen from the example spectra in Fig. 5, the narrow components in the central region are distinct from the underlying broad emission. It is therefore unlikely that we oversubtract emission associated with the narrow component through this procedure. The resulting maps are shown in the bottom row of Fig. 3. In Section 6, we discuss alternative approaches to model the broad nuclear emission.

To measure a redshift for GN20 from the $H\alpha$ data, we rely on the kinematics, since the one-dimensional (1D), integrated emission line profile is skewed by the bright, off-centre region in the north-west. After initially deriving the $H\alpha$ kinematics as described earlier assuming $z = 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 $v_{\text{obs}} \sim 0$ km s$^{-1}$ along the kinematic major axis. We adopt a velocity shift of +75 km s$^{-1}$. Given the median uncertainties on the fitted velocities of 11 km s$^{-1}$ (including in the regions of minimum and maximum velocities), this corresponds to a redshift of $z_{H\alpha} = 4.05\,374 \pm 0.00\,075.$$^3$

<footnote>$^3$As described in the Section 3.2, in our dynamical modelling we allow for an additional velocity shift of $\pm 100$ km s$^{-1}$. Our fiducial model has a velocity shift of $+28 \pm 2$ km s$^{-1}$. This would correspond to a redshift of $z_{H\alpha, \text{model}} = 4.05\,325 \pm 0.00\,075$. 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.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top panel: Voronoi-binned maps (target S/N ≥ 15) of the observed Hα flux (left panel), velocity (middle panel), and velocity dispersion (right panel) 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 σ<sub>obs</sub> ~ 50–150 km s<sup>-1</sup> 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 v<sub>obs</sub> ~ 0 km s<sup>-1</sup> as a dashed line in the bottom-middle panel (repeated in the bottom-right panel for reference).</figcaption></figure>

## 3.2 Dynamical modelling

To construct dynamical models, we use DYSMALPY (Davies, Tacconi & Genzel 2004a, b; Cresci et al. 2009; Davies et al. 2011; Wuyts et al. 2016; Genzel et al. 2017, 2023; Lang et al. 2017; Übler et al. 2018; Price et al. 2021; Lee et al., submitted), 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 for GN20 informed through existing multiwavelength 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).

We fix the following structural parameters to the best-fitting values derived by Colina et al. (2023): the disc effective radius $R_{e,disc} = 3.6$ kpc, the disc Sérsic index $n_{S,disc} = 0.42$, and the minor-to-major axis ratio $b/a = 0.80$. We adopt the convention that in the face-on case, the galaxy rotates in counter-clockwise direction for $i = 0°$, and in clockwise direction for $i = 180°$. Assuming a thick disc with ratio of scale height to scale length of $q_0 = 0.2$ (e.g. Wuyts et al. 2016; Genzel et al. 2017), $b/a = 0.80$ corresponds 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 $i = 142°$. We include a Gaussian prior

on the total baryonic mass centred on $log(M_{bar,tot}/M_{⊙}) = 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 ($M_{*} = 1.1 \times 10^{11} M_{⊙}$; Tan et al. 2014) and of the molecular gas mass ($M_{H2} = 1.3 \times 10^{11} M_{⊙}$; Hodge et al. 2012, here assuming $\alpha_{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 \times 10^{10} M_{⊙}$. This is about 10 per cent of our total baryonic mass estimate. We measure the effective radius of the nuclear component through a 2D Gaussian fit to the Hα + [N II] line map, 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 $R_{e,bulge} = 0.8$ kpc, $n_{S,bulge} = 4$, and a bulge-to-total baryonic mass ratio of $B/T = 0.1$. We include a Navarro, Frenk & White (1996) dark matter halo with a concentration parameter of $c = 3.1$ (the expected value for a halo of $log(M_{halo}/M_{⊙}) = 13$ at this redshift, following Dutton & Macciò 2014). We fit for the dark matter fraction within the effective radius, $f_{DM}(r < R_{e})$, from which the total halo mass is inferred (see Price et al. 2021). We use a flat prior $f_{DM}(< R_{e}) = [0; 1]$. Further to the aforementioned priors, we include a flat prior on the intrinsic velocity dispersion, $\sigma_0 = [30; 150] \text{ km s}^{-1}$, which is assumed to be constant and isotropic throughout the disc. In addition to the total baryonic mass $M_{bar}$, the dark matter fraction at the effective radius $f_{DM}(< R_{e})$, and the intrinsic velocity dispersion $\sigma_0$,
we allow the position angle (PA) to vary in the range PA = [0°; 50°], and we allow the model to adjust the systemic velocity by ±100 km s⁻¹.

We model the point-spread function (PSF) as a Gaussian with FWHM = 0.12 arcsec, corresponding to the approximate PSF (along slicers) at the wavelength of Hα (see D'Eugenio et al. 2023). For the instrumental dispersion, we adopt σ<sub>inst</sub> = 57 km s⁻¹, which corresponds to the nominal spectral resolution at the wavelength of Hα (Jakobsen et al. 2022).

As 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 on to the Voronoi bin with $v_{\text{obs}} \sim 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. We find the best fit through Markov chain Monte Carlo (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-fitting model to our data in the figures below, we show the median values of the model per Voronoi bin.

As motivated in Section 6, we further construct a model with the same parameters as described earlier, 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 $v_r$ in steps of 10 km s⁻¹ from 0 to 200 km s⁻¹, using least-squares minimization (Markwardt 2009). We repeat the fit with the so identified preferred inflow velocity ($v_r = 130$ km s⁻¹) using MCMC. The results of this model are presented in Appendix A.

## 4 COMPLEX Hα EMISSION

In Fig. 1, we show a line map of the emission integrated in the range of 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 Fig. 2) and a velocity difference of about +750 km s⁻¹, respectively, a redshift of $z_{\text{Hα,GN20b}} = 4.06\,532 \pm 0.00\,003$. GN20b is also detected in continuum in the prism observations.

The 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 Fig. 2). Further to the north-west of this region, we detect two fainter clumps (dashed circles in Fig. 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 fig. 2), which is also prominent in the NIRSpec data above ~5 µm (see Fig. 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 Fig. 1).

The 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) and 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 utilizing 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.

The two Hα clumps in the north-west appear also visible in CO(2–1) 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 Fig. 2). Some F560W flux extends from the main disc along the same direction of the south-east extension (see Colina et al. 2023). Clumpy UV emission is also seen here; however, this likely stems from a $z \sim 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 SF regions in the outer disc or a protruding spiral arm. Alternatively, they could correspond to fresh material in the process of accretion on to 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 Fig. 3, we show Hα velocity and velocity dispersion maps based on our fits to the Voronoi-binned cube (top panel: one component; bottom panel: stitched narrow component, see Section 3 for details). We find a large velocity gradient across GN20 of about $\Delta v_{\text{obs}} \sim 610$ km s⁻¹, and observed velocity dispersion values of $\sigma_{\text{obs}} \sim 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 Fig. 3, we indicate $v_{\text{obs}} \sim 0$ km s⁻¹ by dashed lines, illustrating a twist of the apparent kinematic minor axis. Along the north-western $v_{\text{obs}} \sim 0$ km s⁻¹ twist, we observe slightly elevated velocity dispersions, as evident from the bottom-right panel of Fig. 3.

In 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-fitting
model and residuals are shown in Fig. 4. The model provides a reasonable fit to the data, with median offsets in velocity and velocity dispersion of $\Delta v_{\text{med}} = 5.5 \text{ km s}^{-1}$ and $\Delta \sigma_{\text{med}} = 14.9 \text{ km s}^{-1}$. From this model, we measure a rotation velocity $v_{\text{rot}}(R_e = 3.6 \text{ kpc}) = 469 \text{ km s}^{-1}$ (circular velocity $v_c(R_e) = 496 \text{ km s}^{-1}$), reaching a maximum value of $v_{\text{max}} = 531 \text{ km s}^{-1}$ ($v_{c,\text{max}} = 574 \text{ km s}^{-1}$) at $r = 6.2 \text{ kpc}$ ($r = 6.8 \text{ kpc}$). The circular velocity accounts for the effects of pressure support from turbulent motions, and is defined as $v_c^2(r) = v_{\text{rot}}^2(r) + 2\sigma_0^2 r/R_d$, with $R_d$ being the disc scale length (see Burkert et al. 2010, 2016). We constrain an intrinsic velocity dispersion of $\sigma_0 = 89 \text{ km s}^{-1}$. For a galaxy at $z \sim 4.055$, from the relation by Übler et al. (2019), calibrated based on ground-based observations up to $z \sim 3.5$, we would expect an intrinsic ionized gas velocity dispersion of $\sigma_0 \sim 63 \pm 19 \text{ km s}^{-1}$, somewhat lower than our best-fitting result. However, we note that some regions in the outer disc of GN20 show comparable dispersion values in $\sigma_{\text{obs}}$. Evaluating the ratio of maximum rotation velocity to intrinsic velocity dispersion, we find $v_{\text{rot}}(R_e)/\sigma_0 = 5.3$ ($v_{\text{rot,\text{max}}}/\sigma_0 = 6.0$). The dynamical mass enclosed within one $R_e$ is $\log(M_{\text{dyn}}(< R_e)/M_{\odot}) = 11.3$, and within two $R_e$ is $\log(M_{\text{dyn}}(< 2R_e)/M_{\odot}) = 11.7$. The inferred values from our best-fitting model are reported in Table 1.

To compare our modelling results with 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 $i = 45^{\circ}$, Carilli et al. (2010) find a rotation velocity of $v_{\text{rot}}(r \sim 4 \text{ kpc}) = 570 \text{ km s}^{-1}$ based on CO(6–5) molecular gas kinematics, with an enclosed dynamical mass of $M_{\text{dyn}}(r < 4 \text{ kpc}) = 3 \times 10^{11} M_{\odot}$. Repeating our modelling with the inclination fixed to $i = 135^{\circ}$ (i.e. corresponding to their $i = 45^{\circ}$), we find a lower value of $v_{\text{rot}}(r \sim 4 \text{ kpc}) = 442 \text{ km s}^{-1}$, but a comparable enclosed total mass of $M_{\text{dyn}}(r < 4 \text{ kpc}) = 2 \times 10^{11} M_{\odot}$.

Hodge et al. (2012) and Bik et al. (2024) adopt an inclination of $i = 30^{\circ}$. From the CO(2–1) kinematics, Hodge et al. (2012) inferred $v_{\text{max}} = 575 \pm 100 \text{ km s}^{-1}$ and $\sigma_0 = 100 \pm 30 \text{ km s}^{-1}$. Bik et al. (2024) found $v_{\text{max}}(r < 4 \text{ kpc}) = 550 \pm 40 \text{ km s}^{-1}$ and $\sigma_m = 145 \pm 53 \text{ km s}^{-1}$ from Paα kinematics, where $\sigma_m$ represents an upper limit on the intrinsic flux-weighted velocity dispersion. Repeating our modelling with the inclination fixed to $i = 150^{\circ}$ (i.e. corresponding to their $i = 30^{\circ}$), we find a rotation velocity of $v_{\text{rot}}(R_e = 3.6 \text{ kpc}) = 551 \text{ km s}^{-1}$ ($v_{\text{max}} = 646 \text{ km s}^{-1}$), and $\sigma_0 = 89 \text{ km s}^{-1}$. Our modelling results obtained by fixing $i = 150^{\circ}$ 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 $v_{\text{max}}/\sigma_0|_{\text{CO(2-1)}} = 5.8$, and Bik et al. (2024) find a somewhat lower value of $v_{\text{max}}/\sigma_m|_{\text{Paα}} = 3.8 \pm 1.4$, however based on the upper limit on the flux-weighted velocity dispersion. Using $i = 150^{\circ}$, we find even higher values of $v_{\text{max}}/\sigma_0 = 7.3$, yet at the effective radius (fixed to $R_e = 3.6 \text{ kpc}$), we find $v_{\text{rot}}(R_e)/\sigma_0 = 6.2$. This is in general agreement with the literature results. Hodge et al. (2012) infer a dynamical mass of $M_{\text{dyn}} = 5.4 \pm 2.4 \times 10^{11} M_{\odot}$. Fixing $i = 150^{\circ}$, we find $M_{\text{dyn}}(r < R_e) = 2.6 \times 10^{11} M_{\odot}$ and $M_{\text{dyn}}(r < 2R_e) = 7.0 \times 10^{11} M_{\odot}$, in broad agreement with the results by Hodge et al. (2012).

Multiphase gas kinematic measurements at high redshift are still rare. Recently, Parlanti et al. (2024) measured $v_{\text{rot}}(R_e = 3 \text{ kpc})/\sigma_0|_{\text{H\alpha}} \sim 9.2$ in an obscured AGN at $z = 4.76$, with a higher value of $v_{\text{rot}}(R_e)/\sigma_0|_{\text{CII}} \sim 15.8$ measured from the [C II] line (Lelli et al. 2021). In this case, the authors speculate that the AGN might have deposited energy in the ionized gas phase, leading to higher dispersion values and therefore lower rotational

support. In SPT0311–58, which is at the core of a massive proto-cluster, Arribas et al. (2024) find $v_{\text{rot,\text{max}}}/\sigma_0|_{\text{CII}} \sim 1$ based on [O III] $\lambda5007$ kinematics. Their median intrinsic velocity dispersion is $\sigma_0 = 113 \pm 19 \text{ km s}^{-1}$, and in some regions 2–3 times higher than the velocity dispersion measured in [C II]. This likely indicates that SPT0311–58 at $z = 6.9$ is in an earlier phase of disc formation compared with 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 Fig. 4. These residuals are close to the Hα clumps identified in Fig. 1. Furthermore, the best-fitting kinematic major axis appears tilted with respect to the velocity minimum and maximum in the data (see top-left panel in Fig. 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 substructure like a bar or spiral arms, or outflows (e.g. van der Kruit & Allen 1978; Roberts, Huntley & van Albada 1979; van Albada & Roberts 1981; Shlosman, Frank & Begelman 1989; Athanassoula 1992; Wada & Habe 1992; Bournaud & Combes 2002; Binney & Tremaine 2008; Tsukui et al. 2024).

GN20 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 $\sim 150$ and $\sim 200 \text{ 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 $Q = -3.9$, while sizeable effects are generally expected for $Q \ge -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 MIR imaging by Colina et al. (2023) and Crespo Gómez et al. (2024) may also suggest a past gravitational interaction.

As 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 Fig. 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.

Finally, 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 MIR F560W data by Colina et al. (2023) does not provide evidence for a bar, though some substructure is revealed through their residual analysis which
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Our best-fitting dynamical model including a disc, bulge, and halo. Top (bottom) panels show the observed velocity (velocity dispersion) field (left panel), the best-fitting model (middle-left panel), residuals (data-model; middle-right panel), and goodness-of-fit (data-model)/uncertainties (right panel). The cross and line indicate the centre and best-fitting 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 $\Delta v_{\text{med}} = 5.5$ km s$^{-1}$ and an rms velocity difference of $\Delta v_{\text{rms}} = 49.3$ km s$^{-1}$. Corresponding values for the velocity dispersion are $\Delta\sigma_{\text{med}} = 14.9$ km s$^{-1}$ and $\Delta\sigma_{\text{rms}} = 34.7$ km s$^{-1}$. Yet, we observe strong residuals particularly in the north-western region, indicating a disturbance of the kinematics for instance through accretion or external perturbations.</figcaption></figure>

<table><caption>**Table 1.** Results from our best-fitting dynamical model including a thick disc, bulge, and dark matter halo. For $\sigma_0$, we adopt minimum uncertainties of 10 km s$^{-1}$. The effective radius is fixed to $R_e = 3.6$ kpc.</caption><thead><tr><th></th><th>Fitted parameters</th></tr></thead><tbody><tr><td>log(M<sub>bar</sub>/M<sub>☉</sub>)</td><td>11.42<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>σ<sub>0</sub> (km s<sup>-1</sup>)</td><td>89 ± 10</td></tr><tr><td>f<sub>DM</sub>(&lt; R<sub>e</sub>)</td><td>0.30 ± 0.07</td></tr><tr><td>PA<sub>kin</sub></td><td>34.7<sup>+0.6</sup><sub>-0.3</sub></td></tr><tr><td>Velocity shift (km s<sup>-1</sup>)</td><td>29 ± 2</td></tr><tr><td>χ<sup>2</sup><sub>red</sub></td><td>4.9728</td></tr><tr><td></td><td>Derived values</td></tr><tr><td>v<sub>rot</sub>(R<sub>e</sub>) (km s<sup>-1</sup>)</td><td>469</td></tr><tr><td>v<sub>circ</sub>(R<sub>e</sub>) (km s<sup>-1</sup>)</td><td>496</td></tr><tr><td>v<sub>rot,max</sub> (km s<sup>-1</sup>)</td><td>531</td></tr><tr><td>v<sub>rot</sub>(R<sub>e</sub>)/σ<sub>0</sub></td><td>5.3</td></tr><tr><td>log(M<sub>dyn</sub>(&lt; 2R<sub>e</sub>)/M<sub>☉</sub>)</td><td>11.68</td></tr></tbody></table>

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 $z \sim 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 Fig. 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 Fig. A1. This model, which includes an inflow with $v_r = 130$ km s$^{-1}$, indeed provides a somewhat better fit to the observed kinematics ($\Delta\chi^2_{\text{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\alpha + [N\text{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\alpha$ and $[N\text{II}]\lambda\lambda6548, 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\alpha$ in addition to the narrow line components for $H\alpha + [N\text{II}]\lambda\lambda6548, 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 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 panel), and to the integrated spectrum extracted over a larger aperture where broad emission is visible (bottom panel) in the left panels of Fig. 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$^{-1}$ are typical of AGN-driven winds (e.g. Veilleux, Cecil & Bland-
<figure><img src="image_7.png" /><figcaption>**Figure 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 panel: Fits including a set of broad components for $H\alpha$ + [N II] $\lambda\lambda6548, 6583$, 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$^{-1}$). Right panel: Fits including a broad component only for $H\alpha$, i.e. interpreting the broad emission as the broad-line region (BLR) of an accreting black hole. The two fitting set-ups 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.</figcaption></figure>

Hawthorn 2005; Fabian 2012; Genzel et al. 2014; Heckman & Best 2014; Carniani et al. 2015; Harrison et al. 2016; Rupke, Gültekin & Veilleux 2017; Förster Schreiber et al. 2019; Veilleux et al. 2020, and references therein). From the large integrated aperture, we calculate a maximum outflow velocity of $v_{out} = (v_{broad}) + 2\sigma_{broad} = 970 \pm 220$ km s$^{-1}$, where $\sigma_{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).

Assuming a photoionized, constant-velocity spherical outflow (Genzel et al. 2011; Newman et al. 2012; Davies et al. 2019, 2020; Förster Schreiber et al. 2019; Cresci et al. 2023) of the extent $R_{out}$ equal to the aperture radius (~2 kpc), with an electron density of $n_{e,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 $\dot{M}_{out,ion} = 0.7 \pm 0.3$ M$_\odot$ yr$^{-1}$. 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 are 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 warm ionized gas phase (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 $z = 4.76$, Parlanti et al. (2024) find a mass outflow rate of $\dot{M}_{out,ion} = 11_{-5.5}^{+57}$ M$_\odot$ yr$^{-1}$.

## 6.2 The BLR scenario

In the right panels of Fig. 5, we show the corresponding fits including only one broad component fixed at the position of the narrow $H\alpha$ emission in the central Voronoi bin (see upper panel). The fit 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 with the fits interpreting the broad emission as an outflow ($\Delta$BIC$_{outflow-BLR,central} = 14$, $\Delta$BIC$_{outflow-BLR,large aperture} = 9$). From the properties of the (unobscured) putative BLR component, we would derive a black hole mass of $\log(M_{\bullet}/M_{\odot}) = 7.3 \pm 0.4$ following Reines & Volonteri (2015) and using the fit to the large aperture. For a galaxy of the stellar mass of GN20 ($M_{\bullet} \sim 1.1 \times 10^{11}$ M$_\odot$), a value of $\log(M_{\bullet}/M_{\odot}) \sim 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(M_{\bullet}/M_{\odot}) \gtrsim 7.3$.

Assuming that GN20 hosts an obscured AGN, Riechers et al. (2014) estimated an Eddington limit for its black hole mass of $\log(M_{\bullet}^{Edd}/M_{\odot}) = 8.1-8.5$ based on its 6 µm continuum luminosity
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left panel: 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 photoionization by an AGN. Right panel: 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 photoionized by SF. The black cross indicates the galaxy centre as adopted for the dynamical modelling.</figcaption></figure>

and an upper limit on the 2–10 keV X-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 towards the south-west.
This is illustrated in the left panel of Fig. 6, where we show a map of
the narrow-component [N II] λ6583/Hα ratio based on the Voronoi-
binned map (see Section 3). Values of [N II] λ6583/Hα > 0.6
may be indicative of contributions through shocks, α-enhanced
evolved stellar populations, or photoionization by an AGN (Baldwin,
Phillips & Terlevich 1981; Kewley et al. 2001; Kauffmann et al. 2003;
Byler et al. 2019).

The BPT diagram (Baldwin, Phillips & Terlevich; Baldwin et al. 1981) comparing [N II] λ6583/Hα with [O III] λ5007/Hβ can provide insights into the dominant ionization 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 diagnostic diagram, we use the equivalent width of Hα versus the [N II] λ6583/Hα ratio (the 'WHAN' diagram; 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 consistent with being photoionized by 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 Fig. 6.

**7 SUMMARY AND DISCUSSION**

We have presented an analysis of the rest-frame optical kinematics
in the large z ~ 4.055 SMG GN20 using JWST/NIRSpec-IFU data.
In the following, we discuss and summarize our main findings.

**7.1 Hα 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 Fig. 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 Fig. 2.

**7.2 Large-scale Hα kinematics**

We find smooth, large-scale rotation in Hα (Δ$v_{obs} \sim 610$ km s⁻¹ and
$\sigma_{obs,disc} \sim 50–150$ km s⁻¹; Fig. 3). We construct a 3D 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 2D kinematics. From our
best-fitting model (Fig. 4), we find $v_{rot}(R_e) = 469$ km s⁻¹, $\sigma_0 =
89$ km s⁻¹, $v_{rot}(R_e)/\sigma_0 = 5.3$, and log(Mdyn(< 2R_e)/M☉) ~ 11.7.
These results are largely consistent with previous analyses of
the molecular and ionized gas kinematics constrained by Very
Large Array (VLA) observations of CO(2–1), Plateau de Bure
Interferometer (PdBI) observations of CO(6–5), and MIRI/MRS
observations of Paα (Carilli et al. 2010; Hodge et al. 2012; Bik
et al. 2024).

**7.3 Non-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-fitting dynamical model (Fig. 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
$v_r = 130$ km s⁻¹, which indeed provides a slightly better fit to the
data (Δχ²=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 (Fig. A1).

**7.4 AGN 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 outflow in Hα and [N II] (v<sub>out</sub> = 970 ± 220 km s⁻¹), or as a BLR in Hα emission, corresponding to a black hole of mass log(M<sub>•</sub>/M☉) ≳ 7.3 (Fig. 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 photoionization by an AGN (Fig. 6).
## 7.5 Evidence for black hole feeding and feedback

Our kinematic and dynamical analysis provides tentative evidence that gas may be channelled into the nuclear region of GN20. These inflows could potentially fuel central SF and/or accretion on to the central black hole. Indeed we find evidence for an AGN in the centre of GN20 in the form of broad nuclear emission, signalling either AGN-driven outflows or an actively accreting black hole, and Seyfert-like ionization 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 < z < 2.5$, signatures of both ionized and molecular gas inflows have been revealed in deep observations of massive main-sequence galaxies (Price et al. 2021; Genzel et al. 2023).

In 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 ionization and possibly high-velocity outflows (i.e. feedback by the AGN), or the BLR of an actively accreting black hole. This work together with the recent studies by, e.g. Arribas et al. (2024), Jones et al. (2024), and Parlanti et al. (2024), highlight the potential of JWST/NIRSpec-IFS to unveil detailed, spatially resolved properties of $z > 4$ dusty SF 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 characterized through their bright submillimetre emission, JWST opened up new avenues to study their formation and evolution. As argued by several authors, the rapid evolution, high ongoing SF, and dense environment make GN20 a prime progenitor candidate for present-day massive ETGs, 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, 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 research 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 RISEand-FALL. 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 have 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

* Álvarez-Márquez J. et al., 2023, A&A, 671, A105

* Arribas S. et al., 2024, A&A, 688, A146

* Athanassoula E., 1992, MNRAS, 259, 345

* Audibert A. et al., 2019, A&A, 632, A33

* Avery C. R. et al., 2022, MNRAS, 511, 4223

* Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5

* Barger A. J., Cowie L. L., Sanders D. B., Fulton E., Taniguchi Y., Sato Y., Kawara K., Okuda H., 1998, Nature, 394, 248

* Baron D., Netzer H., Lutz D., Prochaska J. X., Davies R. I., 2022, MNRAS, 509, 4457

* Bell E. F., McIntosh D. H., Katz N., Weinberg M. D., 2003, ApJS, 149, 289

* Belli S. et al., 2024, Nature, 630, 54

* Bik A. et al., 2024, A&A, 686, A3

* Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, Princeton

* Blain A. W., Smail I., Ivison R. J., Kneib J. P., Frayer D. T., 2002, Phys. Rep., 369, 111

* Böker T. et al., 2022, A&A, 661, A82

* Bournaud F., Combes F., 2002, A&A, 392, 83

* Burkert A. et al., 2010, ApJ, 725, 2324

* Burkert A. et al., 2016, ApJ, 826, 214

* Byler N., Dalcanton J. J., Conroy C., Johnson B. D., Choi J., Dotter A., Rosenfield P., 2019, AJ, 158, 2

* Cappellari M., Copin Y., 2003, MNRAS, 342, 345

* Carilli C. L. et al., 2010, ApJ, 714, 1407

* Carilli C. L., Hodge J., Walter F., Riechers D., Daddi E., Dannerbauer H., Morrison G. E., 2011, ApJ, 739, L33

* Carniani S. et al., 2015, A&A, 580, A102

* Casasola V., Hunt L. K., Combes F., García-Burillo S., Neri R., 2011, A&A, 527, A92

* Casey C. M. et al., 2009, MNRAS, 400, 670

* Casey C. M., Narayanan D., Cooray A., 2014, Phys. Rep., 541, 45

* Chapman S. C., Blain A. W., Ivison R. J., Smail I. R., 2003, Nature, 422, 695

* Chapman S. C., Blain A. W., Smail I., Ivison R. J., 2005, ApJ, 622, 772

* Cid Fernandes R., Stasińska G., Schlickmann M. S., Mateus A., Vale Asari N., Schoenell W., Sodré L., 2010, MNRAS, 403, 1036

* Cid Fernandes R., Stasińska G., Mateus A., Vale Asari N., 2011, MNRAS, 413, 1687

* Colina L. et al., 2023, A&A, 673, L6

* Combes F. et al., 2014, A&A, 565, A97

* Combes F. et al., 2019, A&A, 623, A79

* Cortzen I. et al., 2020, A&A, 634, L14

* Cresci G. et al., 2009, ApJ, 697, 115

* Cresci G. et al., 2023, A&A, 672, A128

* Crespo Gómez A. et al., 2024, preprint (arXiv:2402.18672)

* D'Eugenio F. et al., 2023, preprint (arXiv:2308.06317)
Daddi E. et al., 2009, ApJ, 694, 1517

Dahari O., 1984, AJ, 89, 966

Davies R. I., Tacconi L. J., Genzel R., 2004a, ApJ, 602, 148

Davies R. I., Tacconi L. J., Genzel R., 2004b, ApJ, 613, 781

Davies R. et al., 2011, ApJ, 741, 69

Davies R. L. et al., 2019, ApJ, 873, 122

Davies R. L. et al., 2020, ApJ, 894, 28

Davies R. L. et al., 2024, MNRAS, 528, 4976

Dekel A., Sari R., Ceverino D., 2009, ApJ, 703, 785

Dekel A., Zolotov A., Tweed D., Cacciato M., Ceverino D., Primack J. R., 2013, MNRAS, 435, 999

Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359

Escala A., Larson R. B., 2008, ApJ, 685, L31

Fabian A. C., 2012, ARA&A, 50, 455

Fluetsch A. et al., 2021, MNRAS, 505, 5753

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Förster Schreiber N. M. et al., 2019, ApJ, 875, 21

Fukugita M., Hogan C. J., Peebles P. J. E., 1998, ApJ, 503, 518

Gaia Collaboration, 2023, A&A, 674, A1

Genzel R. et al., 2011, ApJ, 733, 101

Genzel R. et al., 2014, ApJ, 785, 75

Genzel R. et al., 2017, Nature, 543, 397

Genzel R. et al., 2023, ApJ, 957, 48

Hainline L. J., Blain A. W., Smail I., Frayer D. T., Chapman S. C., Ivison R. J., Alexander D. M., 2009, ApJ, 699, 1610

Harrison C. M. et al., 2016, MNRAS, 456, 1195

Heckman T. M., Best P. N., 2014, ARA&A, 52, 589

Herrera-Camus R. et al., 2019, ApJ, 871, 37

Hodge J. A., da Cunha E., 2020, R. Soc. Open Sci., 7, 200556

Hodge J. A., Carilli C. L., Walter F., de Blok W. J. G., Riechers D., Daddi E., Lentati L., 2012, ApJ, 760, 11

Hodge J. A., Carilli C. L., Walter F., Daddi E., Riechers D., 2013, ApJ, 776, 22

Hodge J. A., Riechers D., Decarli R., Walter F., Carilli C. L., Daddi E., Dannerbauer H., 2015, ApJ, 798, L18

Hogg D. W. et al., 2002, AJ, 124, 646

Hughes D. H. et al., 1998, Nature, 394, 241

Hunt L. K. et al., 2008, A&A, 482, 133

Iono D. et al., 2006, ApJ, 640, L1

Izumi T. et al., 2023, Science, 382, 554

Jakobsen P. et al., 2022, A&A, 661, A80

Jones G. C. et al., 2024, A&A, 682, A122

Kakkad D. et al., 2018, A&A, 618, A6

Kauffmann G. et al., 2003, MNRAS, 346, 1055

Kewley L. J., Dopita M. A., Sutherland R. S., Heisler C. A., Trevena J., 2001, ApJ, 556, 121

Krumholz M., Burkert A., 2010, ApJ, 724, 895

Krumholz M. R., Burkhart B., Forbes J. C., Crocker R. M., 2018, MNRAS, 477, 2716

Lang P. et al., 2017, ApJ, 840, 92

Lelli F., Di Teodoro E. M., Fraternali F., Man A. W. S., Zhang Z.-Y., De Breuck C., Davis T. A., Maiolino R., 2021, Science, 371, 713

Magnelli B. et al., 2013, A&A, 553, A132

Maiolino R., Alonso-Herrero A., Anders S., Quillen A., Rieke M. J., Rieke G. H., Tacconi-Garman L. E., 2000, ApJ, 531, 219

Maiolino R. et al., 2024, Nature, 627, 59

Marasco A. et al., 2023, A&A, 670, A92

Markwardt C. B., 2009, in Bohlender D. A., Durand D., Dowler P., eds, ASP Conf. Ser. Vol. 411, Astronomical Data Analysis Software and Systems XVIII. Astron. Soc. Pac., San Francisco, p. 251

Maseda M. V. et al., 2024, preprint (arXiv:2403.05506)

Morrison G. E., Owen F. N., Dickinson M., Ivison R. J., Ibar E., 2010, ApJS, 188, 178

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Newman S. F. et al., 2012, ApJ, 752, 111

Noordermeer E., 2008, MNRAS, 385, 1359

Oteo I. et al., 2018, ApJ, 856, 72

Overzier R. A., 2016, A&AR, 24, 14

Parlanti E. et al., 2024, A&A, 684, A24

Pavesi R. et al., 2018, ApJ, 861, 43

Perna M., Lanzuisi G., Brusa M., Cresci G., Mignoli M., 2017, A&A, 606, A96

Perna M. et al., 2023, A&A, 679, A89

Planck Collaboration VI, 2020, A&A, 641, A6

Pope A., Borys C., Scott D., Conselice C., Dickinson M., Mobasher B., 2005, MNRAS, 358, 149

Pope A. et al., 2006, MNRAS, 370, 1185

Price S. H. et al., 2021, ApJ, 922, 143

Reines A. E., Volonteri M., 2015, ApJ, 813, 82

Renzini A., 2006, ARA&A, 44, 141

Riechers D. A. et al., 2014, ApJ, 786, 31

Rieke G. H. et al., 2015, PASP, 127, 665

Roberts-Borsani G. W., 2020, MNRAS, 494, 4266

Roberts W. W., Jr, Huntley J. M., van Albada G. D., 1979, ApJ, 233, 67

Rupke D. S. N., Veilleux S., 2013, ApJ, 768, 75

Rupke D. S. N., Gültekin K., Veilleux S., 2017, ApJ, 850, 40

Shlosman I., Frank J., Begelman M. C., 1989, Nature, 338, 45

Simpson J. M. et al., 2014, ApJ, 788, 125

Smail I., Ivison R. J., Blain A. W., 1997, ApJ, 490, L5

Speights J. C., Rooke P. C., 2016, ApJ, 826, 2

Swinbank A. M. et al., 2014, MNRAS, 438, 1267

Tan Q. et al., 2014, A&A, 569, A98

Toomre A., 1964, ApJ, 139, 1217

Tsukui T., Wisnioski E., Bland-Hawthorn J., Mai Y., Iguchi S., Baba J., Freeman K., 2024, MNRAS, 527, 8941

Übler H. et al., 2018, ApJ, 854, L24

Übler H. et al., 2019, ApJ, 880, 48

van Albada G. D., Roberts W. W. Jr, 1981, ApJ, 246, 740

van Dokkum P. G., 2001, PASP, 113, 1420

van der Kruit P. C., Allen R. J., 1978, ARA&A, 16, 103

Veilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 43, 769

Veilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, A&AR, 28, 2

Venturi G. et al., 2018, A&A, 619, A74

Verley S. et al., 2007, A&A, 472, 121

Wada K., Habe A., 1992, MNRAS, 258, 82

Walter F. et al., 2012, Nature, 486, 233

Wright G. S. et al., 2015, PASP, 127, 595

Wuyts S. et al., 2012, ApJ, 753, 114

Wuyts S. et al., 2016, ApJ, 831, 149

## APPENDIX A: DYNAMICAL MODEL INCLUDING A FIXED RADIAL INFLOW

In Fig. A1, we show the results of our dynamical modelling including a uniform, planar radial inflow of fixed velocity $v_r = 130 \text{ km s}^{-1}$, 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 $v_r$ in steps of $10 \text{ km s}^{-1}$ from 0 to 200 km s$^{-1}$, 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 inflow fixed at $v_r = 130 \text{ km s}^{-1}$. Our modelling results are listed in Table A1.

While 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 with our fiducial model presented in Section 5 by $\Delta\chi^2_{\text{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.

We point out the change in the kinematic position angle of the best-fitting model including inflow with respect to our fiducial model. The best-fitting rotation velocities are reduced. In this model, we find little contribution of dark matter to the dynamics on galactic scales.

We 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-Q parameter (Toomre 1964) following Binney & Tremaine (2008), Escala & Larson (2008), and Dekel, Sari & Ceverino (2009), at the radius where our best-fitting circular velocity curve reaches its maximum (see Übler et al. 2019), and find $Q_{gas} \sim 0.49$. Following Genzel et al. (2023), we then obtain $v_r = 102$ km s$^{-1}$, which is somewhat lower than our model-derived value. If we approximate the inflow velocity through $v_r \sim f_{gas}^2 v_c$ (Genzel et al. 2023), with $v_c$ measured at the effective radius, we obtain $v_r = 138$ km s$^{-1}$, comparable with our model-derived 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.

Evidence for inflows in massive, high-redshift ($z \sim 1-2.5$) SF 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 $v_r \sim 30-120$ km s$^{-1}$, comparable with the putative inflow velocity of $v_r = 130$ km s$^{-1}$ in GN20.

In Fig. A2, we show the difference between our fiducial velocity model as presented in the top-middle-left panel of Fig. 3, and the velocity model including the radial inflow (top-middle-left panel of Fig. A1). The differences between the two models are most pronounced in the centre and in the outer regions.

<table><caption>**Table A1.** Results from our best-fitting dynamical model including a thick disc, bulge, and dark matter halo, and a planar radial inflow of velocity $v_r = 130$ km s$^{-1}$. For $\sigma_0$, we adopt minimum uncertainties of 10 km s$^{-1}$, and for the total baryonic mass we adopt minimum uncertainties of 0.1 dex. For $f_{DM}$, we give the 3$\sigma$ upper limit.</caption><thead><tr><th></th><th>Fitted parameters</th></tr></thead><tbody><tr><td>log(M<sub>bar</sub>/M<sub>☉</sub>)</td><td>11.5 ± 0.1</td></tr><tr><td>σ<sub>0</sub> (km s<sup>-1</sup>)</td><td>89 ± 10</td></tr><tr><td>f<sub>DM</sub>(&lt; R<sub>e</sub>)</td><td>&lt; 0.07</td></tr><tr><td>PA<sub>kin</sub></td><td>17.1<sup>+0.4</sup><sub>-0.5</sub></td></tr><tr><td>Velocity shift (km s<sup>-1</sup>)</td><td>31 ± 2</td></tr><tr><td>χ<sup>2</sup><sub>red</sub></td><td>4.6414</td></tr><tr><td></td><td>Derived values</td></tr><tr><td>v<sub>rot</sub>(R<sub>e</sub>) (km s<sup>-1</sup>)</td><td>443</td></tr><tr><td>v<sub>circ</sub>(R<sub>e</sub>) (km s<sup>-1</sup>)</td><td>472</td></tr><tr><td>v<sub>rot,max</sub> (km s<sup>-1</sup>)</td><td>483</td></tr><tr><td>v<sub>rot</sub>(R<sub>e</sub>)/σ<sub>0</sub></td><td>5.0 ± 0.6</td></tr><tr><td>log(M<sub>dyn</sub>(&lt; 2R<sub>e</sub>)/M<sub>☉</sub>)</td><td>11.51</td></tr></tbody></table>
<figure><img src="image_9.png" /></figure>

**Figure A1.** Same as Fig. 4, but including a planar radial inflow of $v_r = 130 \text{ km s}^{-1}$ in addition to our fiducial disc, bulge, and halo model. Top (bottom) panels show the observed velocity (velocity dispersion) field (left panel), the best-fitting model (middle-left panel), residuals (data-model; middle-right panel), and goodness-of-fit (data-model)/uncertainties (right panel). The cross and line indicate the centre and best-fitting PA. The median velocity offset is $\Delta v_{\text{med}} = 8.8 \text{ km s}^{-1}$ with an rms velocity difference of $\Delta v_{\text{rms}} = 40.0 \text{ km s}^{-1}$. Corresponding values for the velocity dispersion are $\Delta\sigma_{\text{med}} = 15.4 \text{ km s}^{-1}$ and $\Delta\sigma_{\text{rms}} = 34.6 \text{ km s}^{-1}$. 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.

<figure><img src="image_10.png" /></figure>

**Figure A2.** Difference in the model velocities of our fiducial model compared with the model including the radial inflow of $v_r = 130 \text{ km s}^{-1}$. The differences are most pronounced in the centre and in the outer regions.

This paper has been typeset from a TEX/LATEX file prepared by the author.