Clumpy star formation and an obscured nuclear starburst in the luminous dusty z=4 galaxy GN20 seen by MIRI/JWST

Dusty star-forming galaxies emit most of their light at far-IR to mm wavelengths as their star formation is highly obscured. Far-IR and mm observations have revealed their dust, neutral and molecular gas properties. The sensitivity of JWST at rest-frame optical and near-infrared wavelengths now allows the study of the stellar and ionized gas content. We investigate the spatially resolved distribution and kinematics of the ionized gas in GN20, a dusty star forming galaxy at z = 4.0548. We present deep MIRI / MRS integral field spectroscopy of the near-infrared rest-frame emission of GN20. We detect spatially resolved Pa α , out to a radius of 6 kpc, distributed in a clumpy morphology. The star formation rate derived from Pa α (144 ± 9 M ⊙ yr − 1 ) is only 7.7 ± 0 . 5% of the infrared star formation rate (1860 ± 90 M ⊙ yr − 1 ). We attribute this to very high extinction (A V = 17.2 ± 0.4 mag, or A V , mixed = 44 ± 3 mag), especially in the nucleus of GN20, where only faint Pa α is detected, suggesting a deeply buried starburst. We identify four, spatially unresolved, clumps in the Pa α emission. Based on the double peaked Pa α profile we find that each clump consist of at least two sub-clumps. We find mass upper limits consistent with them being formed in a gravitationally unstable gaseous disk. The UV bright region of GN20 does not have any detected Pa α emission, suggesting an age of more than 10 Myrs for this region of the galaxy. From the rotation profile of Pa α we conclude that the gas kinematics are rotationally dominated and the v rot /σ m = 3 . 8 ± 1 . 4 is similar to low-redshift LIRGs. We speculate that the clumps seen in GN20 could contribute to building up the inner disk and bulge of GN20.


Introduction
Dusty Star Forming Galaxies (DSFGs) are extreme infrared luminous galaxies with short intense starburst episodes, forming stars at rates from 500 to several 1000 M ⊙ yr −1 (see Casey et al. 2014, for a review).They are the most luminous starbursts in the Universe and are considered the progenitors of the massive quiescent galaxies between z ∼ 2 − 3 (Valentino et al. 2020).These galaxies emit most of their light at (far) infrared to radio wavelengths and are optically obscured due to a significant dust content which attenuates the light.They are identified as submillimeter galaxies (SMGs) in deep sub-millimeter surveys (e.g.Hughes et al. 1998;Borys et al. 2003;Pope et al. 2006;Vieira et al. 2013) and represent an important phase of cosmic star formation missed by the deep rest-frame optical and ultraviolet surveys (Bouwens et al. 2021(Bouwens et al. , 2022;;Finkelstein et al. 2023;Pérez-González et al. 2023b).More recent surveys with ALMA have concluded that the infrared bright galaxies dominate the cosmic star formation rate below z ∼ 4 (Pérez-González et al. 2005;Gruppioni et al. 2013;Zavala et al. 2021), while at higher redshift their contribution goes down to 10-30% at z ∼ 6 − 7 (Gruppioni et al. 2020;Algera et al. 2022;Barrufet et al. 2023a).
There are several scenarios presented in the literature for the triggering mechanism of the extreme starburst in DSFGs.Many DSFGs are thought to be starburst dominated major mergers.ALMA observations reveal the disturbed morphology and kinematics in several of those galaxies (e.g.Gómez-Guijarro et al. 2018;Riechers et al. 2020;Ginolfi et al. 2020;Spilker et al. 2022;Álvarez Márquez et al. 2023b).Several of these galaxies are located in over densities or are members of proto clusters, making interactions a likely cause for their extreme starburst (e.g.Daddi et al. 2009;Drake et al. 2020).Their merger origin would make the DSFGs high redshift analogues of ultra luminous infrared galaxies (ULIRGs) in the local universe.
The other scenario is driven by the observations of starforming galaxies at high redshift showing gas rich disks (e.g.Genzel et al. 2008;Elmegreen et al. 2009b).In this scenario the star formation is fuelled by cold mode accretion (e.g.Dekel et al. 2009a,b), where the high star formation in the gas rich disk is maintained by smooth infall and accretion of gas.Due to the high gas densities, the disk becomes unstable due to gravitational instabilities (Krumholz et al. 2018), resulting in the formation of giant star forming clumps in the disk (Dekel et al. 2009b;Romeo et al. 2010;Ceverino et al. 2010;Mandelker et al. 2014;Romeo & Agertz 2014) with kpc sizes and masses of ∼10 9 M ⊙ (Genzel et al. 2008;Elmegreen et al. 2009b;Ceverino et al. 2010).These clumps could then migrate further into the center of the galaxy due to dynamical friction over timescales of several hundred Myr (e.g.Mandelker et al. 2014) and contribute to the formation of the thick disk and the bulge growth (Elmegreen et al. 2008(Elmegreen et al. , 2009a;;Ceverino et al. 2010).Observations of individual galaxies support this scenario by finding trends between galactocentric distance and age of the clumps (Adamo et al. 2013;Cava et al. 2018).
Due to the fact that the DSFGs have high extinction, almost all we know about these galaxies is derived from their molec-ular gas and dust content (e.g.Hodge et al. 2012Hodge et al. , 2015) ) tracing mostly the neutral interstellar medium.[CII] observations of DS-FGs trace both the PDR and ionized gas (e.g.Spilker et al. 2022).However, their stellar distribution and ionized gas (without contamination of the PDR emission) properties remained hidden.Many of these galaxies are dark at HST wavelengths, or only reveal very faint emission.Only now with the arrival of JWST we are able to reveal the restframe optical and near-infrared wavelengths via imaging (Colina et al. 2023;Gillman et al. 2023;Barrufet et al. 2023b;Gómez-Guijarro et al. 2023;Cochrane et al. 2023;Pérez-González et al. 2023a;Zavala et al. 2023;Akins et al. 2023) and spatially resolved spectroscopy, where we have for the first time access to emission lines originating in the ionized gas, tracing the most recent star formation (Álvarez Márquez et al. 2023b;Arribas et al. 2023;Jones et al. 2023a;Parlanti et al. 2023).
In this paper we analyze the ionized gas distribution and kinematics of the DSFG GN20.GN20 was identified as a bright 850 µm source in GOODS-North based on deep SCUBA imaging (Pope et al. 2006).GN20 is a DSFG at redshift 4.0548 (Carilli et al. 2011), located in a proto-cluster or galaxy overdensity (Daddi et al. 2009) and has an infrared luminosity of 1.86 × 10 13 L ⊙ .The star formation rate derived from fitting the observed spectral energy distribution (SED), assuming a Chabrier initial mass function (IMF), adds up to 1860 M ⊙ yr −1 (Tan et al. 2014).With a total stellar mass of 1.1 × 10 11 M ⊙ (Tan et al. 2014), GN20 is located well above the z=4 star formation main sequence (Speagle et al. 2014;Caputi et al. 2017).The specific star formation rate of GN20 computes to sSFR = 10 −7.7 yr −1 , just below the starburst line for z = 4 -5 derived by Caputi et al. (2017).
Deep Expanded Very Large Array (EVLA) observations of GN20 in the CO (2-1) line (Carilli et al. 2011;Hodge et al. 2012) revealed a large, rotating molecular gas disk, visible under an inclination of 30 • (M dyn = 5.4 ± 2.4 × 10 11 M ⊙ ), extending as far at 14 kpc in diameter with a very clumpy spatial distribution.Hodge et al. (2012) conclude that the presence of this disk would point to another process than a major merger as a driver for the starburst.This accretion could be enhanced by the interaction with the other galaxies in the overdensity.
Recently Colina et al. (2023) presented deep MIRI F560W imaging observations, probing the rest-frame near-infrared view of GN20.This study revealed for the first time the presence of a large stellar disk off-centred from the bright nucleus.This would, however, suggest that GN20 is involved in an interaction or a merger event, as supported by the presence of a secondary nucleus in the F560W image.
In this paper we present MIRI/MRS observations of GN20, aiming at detecting and characterizing the emission from the ionized gas, tracing the most recent star formation.The MRS observations are focused on the ionized hydrogen emission lines Paα (λ rest = 1.8756 µm) and Paβ (λ rest = 1.2822 µm).Even though intrinsically fainter than Hα, these lines suffer significantly less from extinction and should therefore be easier to detect in galaxies with significant dust content such as DSFGs.In Sect. 2 we present the observations, data reduction and post processing as well as the ancillary data.In Sect. 3 we outline our analysis methods.In Sect. 4 we present the line emission map and analyse the spectral and spatial distribution.Our discussion in the context of the evolution of DSFGs and our conclusions are presented in Section 5 and 6 respectively.
In this paper we adopt the cosmological parameters from the Planck Collaboration et al. ( 2016) and a redshift of GN20 of z=4.0548.With these parameters, 1 ′′ is equal to 7.08 kpc.We use vacuum emission line wavelengths.

MRS observations
GN20 was observed by JWST on November 24th, 2022 using the Mid-Infrared Instrument (MIRI, Rieke et al. 2015;Wright et al. 2015Wright et al. , 2023) ) as part of the European Consortium MIRI Guaranteed Time (proposal ID 1264, PI: L. Colina).This paper presents the data obtained with the Medium Resolution Spectrograph (MRS, Wells et al. 2015;Argyriou et al. 2023).Colina et al. (2023) present the F560W imaging of GN20 obtained while MRS was performing background observations.
The data were obtained in the SLOWR1 readout mode using a 4-point extended source dither pattern, repeated twice, resulting in eight different exposures.Each exposure contained one integration.Each integration contained 40 groups, resulting in 955.6 seconds per integration, summing up to a total of 7645 sec.Additionally, a background observation on an empty sky region, comprising of a 2-point dither with the same observing setup was taken to aid in the background subtraction.
Stage one of the MRS pipeline (Morrison et al. 2023) is executed with mostly standard settings, producing the rate files.As our deep observations contain many groups and are taken in the SLOWR1 readout mode, special care has to be taken to remove the cosmic ray (CR) showers present in the data.The developers version we used contained the latest algorithms to minimise the effect of the CR showers.We turned on the find_showers keyword and, following Bosman et al. (2023) set the time_masked_after_shower and max_extended_radius to 60 seconds and 400 pixels respectively.This procedure successfully removes a number of faint CR showers, however the brightest showers remain present.
From the rate files, we remove the warm pixels by calculating the median of all source and background frames.This median frame has all the cosmic ray events removed and only events persistent over many frames remain present.The source emission is very faint and not detected in each single frame.To remove the warm pixels we applied sigma clipping and updated the data quality (DQ) frame of each individual rate file.The GN20 observations have a lower dark current than the reference dark present in the CRDS used by the pipeline.This results in an offset of the values of the rate files.We measure the median of the pixels in the inter-channel area of the detector and correct for the found offset (∼ 0.16 DN/s).
The second stage of the pipeline (Argyriou et al. 2023;Patapis et al. 2023;Gasman et al. 2023) is run with default parameters, apart from skipping the background subtraction step.After stage two we performed a pixel-by-pixel background subtraction.We calculate a master background frame by median averaging both object and background frames.Due to the different dither positions median average does not contain any source flux.
As a test we created a background frame from the two background observations and compared the rms in the datacube to that created from the background with all the frames.The lowest rms was reached by using the median of all the frames as the background.We also extracted a spectrum of the Paα emission in both cubes and found the same, showing that no source emission is subtracted.The background frame was subtracted from each object frame.Finally, stage three of the pipeline was run, skipping the master background subtraction step.In the outlier rejection step we use a kernel_size of '11, 1' pixels, and a threshold_percent of 99.8 %.By combining all the exposures, 4 3D datacubes are produced (one for each band) using the drizzle algorithm (Law et al. 2023).This produces datacubes with a spatial resolution of 0.13 ′′ per pixel (channel 1) and 0.17 ′′ per pixel (channel 2).The FWHM of the MRS observations is calculated as 0.31× λ (µm)/ 8µm, resulting in 0.37 ′′ (2.6 kpc) at the wavelength of Paα (Argyriou et al. 2023).The spectral resolution of the observation is R= 3774 for Paβ and R=3390 for Paα (Labiano et al. 2021;Jones et al. 2023b).
As already stated, our data on GN20 are taken with the SLOWR1 readout and consist of long integrations, cosmic ray showers affect a significant part of the frames.Even though the updated pipeline significantly improves the removal of the CR showers, not all are removed and due to the faint emission of the target they considerably affect the extracted spectra to be analysed.The CR shower affect rather large parts of the detector and have elliptical shapes, Due to their spatial extent on the rate files, cosmic rays manifest themselves as stripes in the x-direction in the reduced data cubes.As a final step of the data reduction process we employed the algorithm developed by Spilker et al. (2023) to remove the striping due to the CR showers which are not removed by the pipeline.

Archival datasets
In this paper we compare the MRS observations with a multiwavelength dataset to obtain a comprehensive picture of GN20.We use the CO(2-1) Very Large Array (VLA) observations presented in Carilli et al. (2011) and Hodge et al. (2012), with a synthesized beam size of 0.21 ′′ , to compare the distribution and kinematics of the ionized gas to that of the molecular gas.We also use the 880 µm continuum emission obtained with the Plateau de Bure Interferometer (PdBI, Hodge et al. 2015).These data have a synthesized beam size of 0.2 ′′ x 0.3 ′′ .
To relate the ionized gas properties to the stellar distribution in GN20 we use the JWST/MIRI F560W imaging (1.1 µm rest-frame, Colina et al. 2023).Additionally we use HST imaging (PI: Faber, ID: 12442) with the ACS camera (F606W, F775W and F814W) and the WFC3 camera (F105W, F125W and F160W), tracing the UV-optical rest-frame emission.The reduced HST images were downloaded from the Mikulski Archive for Space Telescopes (MAST), astrometically aligned to GAIA DR3 (Gaia Collaboration et al. 2023).

Astrometry
To improve the absolute astrometry of the observations we followed the recommended procedure and obtained simultaneous imaging in the F770W filter.This image was reduced using the JWST pipeline v1.8.2 and context jwst_1030.pmapfollowing the procedure presented in Álvarez Márquez et al. (2023b).The reduced F770W image is astrometrically calibrated using GAIA DR3 (Gaia Collaboration et al. 2023), resulting in an astrometric accuracy of ∼ 70 mas.The offset from the initial astrometry provided by the telescope was calculated and applied to the header of the MRS data cube.
Both the F814W image from HST and the F560W MIRI image are astrometrically calibrated with GAIA DR3 with an uncertainty less than ∼70 mas.The 880mm and CO(2-1) observations were phase referenced to quasars and their absolute astrometry is estimated to be a tenth of their synthesized beam, i.e. 20 mas.This makes the total accuracy of all the data (70 mas) well within the size of a pixel in the MRS data cubes.
As a final check we stacked the cube of channel 2 in the wavelength direction to get a pseudo broadband image, detecting the bright nucleus seen in the F560W imaging Colina et al. (2023) in the continuum.We found that the continuum emission in the MRS cube overlaps with the bright nucleus in the F560W image.

Analysis
In this section we describe the procedure we have developed to extract emission line spectra of both extended and point source emission.Additionally we describe the procedure to create the emission line maps we use to study the spatial distribution of the ionized gas.

Extraction of the spectra
In order to get the integrated emission line profiles we extract the spectrum by integrating over a circular aperture using photutils (Bradley 2023) in each wavelength bin.For extraction of the integrated Paα emission line we use an aperture with a radius of 6 pixels (1.02 ′′ , 7 kpc), covering the observed extent of GN20 in CO(2-1) (Hodge et al. 2012) and the F560W continuum (Colina et al. 2023).No aperture correction is applied as the aperture is relatively large compared to the size of the PSF (0.37 ′′ ) and the emission is spatially extended.The resulting integrated Paα emission line profile is shown in Fig. 1.We calculate the error by measuring the standard deviation (plotted in green) in a sky aperture with the same size, multiplied by the square-root of the number of pixels.We show in red the CO(2-1) spectrum of GN20 from Hodge et al. (2012).The two spectra look fairly similar, with the exception that the Paα emission shows blue shifted emission to −500 km s −1 , while the CO emission is only detected to −400 km s −1 .
To measure the redshift of GN20 from the Paα emission, we fit a single Gaussian profile to the integrated spectrum (Fig. 1).We find z Paα = 4.0553 ± 0.0006, within 1σ from the value derived from the CO(2-1) observations by Carilli et al. (2011).For consistency with the literature we use z=4.0548 throughout the paper.
For the extraction of unresolved sources (clumps) we extract Paα spectra with a radius equal to the instrumental PSF (0.37 ′′ ).We measure the aperture correction for this aperture on the MRS over-sampled and drizzled PSF models from Patapis et al. (2023).We find a correction factor of 1.42.For error determination of the point source extractions we extracted sky spectra at 9 places in the data cube where no galaxy emission is present and use standard deviation of these spectra as an error estimate.

Construction of emission line maps
The integrated Paα spectrum (Fig. 1) shows emission at velocities between −450 to +400 km s −1 .The CO(2-1) emission of GN20 shows a similar range of velocities and shows a velocity gradient due to the global kinematics of the galaxy (Hodge et al. 2012).
In order to reduce the background noise in the emission line map, and keep the extraction window as narrow as possible, we first create a position velocity (PV) diagram using the Paα emission as guide.This line is the brightest and can be detected in the data cube.We sum up the emission over a pseudo-long slit placed over the major axis of the galaxy with a width of 8 pixels (1.4 ′′ ).Hodge et al. (2012) find a position angle (PA) of +25 • .We find that the same PA gives the brightest Paα emission in the PV diagram (Fig. 2) and find no evidence for a difference in PA between the CO and Paα emission.
Based on the observed shape of the Paα emission in the PV diagram, we construct an aperture where the velocity scales as the arctan of the radius, commonly used to describe galaxy velocity profiles (Glazebrook 2013).We adjust the parameters such that all the Paα emission in the PV diagram is included.This function provides us with a central velocity of the Paα emission as a function of position along the major axis.From this we construct a map of central wavelengths of the Paα emission for each spatial pixel in the datacube.We construct the line map summing the flux in an aperture of 9 pixels in the wavelength direction (370 km s −1 ) centered on the derived central wavelength.Parallel to the aperture centered on the Paα line we construct a "background" aperture by shifting the same function with 10 pixels along the radius axis from which we derive the noise statistics.
As Paα is the brightest line, we use the Paα aperture shape to extract line maps for other expected emission lines.The H 2 and [SVI] lines are also located in channel 2. In order to extract emission line maps, we shift the map with the central wavelength

Emission line maps
The extracted Paα emission line map is shown in figure 3. We detect Paα emission in GN20 over a large extent of the galaxy, up to a radius of 6 kpc.The Paα emission shows a very clumpy structure with 4 clumps detected in the outer regions of the disk and very little emission towards the central regions of GN20.We do not detect any significant Paβ emission or any of the other emission lines covered in the two channels.
We derive the error on the detected Paα flux by measuring the standard deviation in a large aperture in the "background" image as described in Sect.3.2.For the Paα line map we find the standard deviation of the background image to be 1.06 × 10 −19 erg s −1 cm −2 per pixel, integrated over 9 wavelength elements.This would translate to a 3 σ detection limit in an extraction aperture equal to the FWHM of the PSF (0.37 ′′ ) of 1.2 × 10 −18 erg s −1 cm −2 .Fig. 3 shows the 2, 3 and 4-σ contours of Paα, derived using this standard deviation.For Paβ, we derive a standard deviation of 1.86 × 10 −19 erg s −1 cm −2 per pixel, translating to a 3 σ detection limit in a aperture equal to the FWHM (0.3 ′′ ) of the PSF of 2.3 × 10 −18 erg s −1 cm −2 .For the other potential emission lines we derive the standard deviation in the same way and find 3 sigma detection limits of 1.5 × 10 −18 erg s −1 cm −2 for [SiVI] and 1.4 × 10 −18 erg s −1 cm −2 for both the H 2 lines.

Integrated star formation rate
We calculate the integrated Paα flux of GN20 by integrating the emission line map over a circular aperture with a radius of 6 pixels (1.02 ′′ ).This results in a total flux of 1.9 ± 0.11 × 10 −17 erg s −1 cm −2 (Table 1).We use the star formation rate calibration from Kennicutt & Evans (2012).Kennicutt & Evans (2012) use a Kroupa (2001) IMF in their calibration, while in the remainder of the paper we use the IMF of Chabrier (2003), yielding identical results (Chomiuk & Povich 2011).We assume an intrinsic Paα over Hα ratio of 0.116 (Osterbrock & Ferland 2006) and we calculate a corresponding unobscured star formation rate (SFR Paα ) of 144 ± 9 M ⊙ yr −1 , assuming zero extinction.
Notes. (a) 1 σ upper limit To derive an upperlimit of the Paβ flux we integrate the flux in the Paβ map in a similar circular aperture (8 pixels, 1.04 ′′ ).This results in a 1-σ upper limit for the integrated flux of Paβ of 3 × 10 −18 erg s −1 cm −2 (equivalent to an SFR Paβ < 43 M ⊙ yr −1 in the case of no extinction).
The 3 σ detection limit of Paβ is 7.89 × 10 −18 erg s −1 cm −2 , comparing this to the observed Paα flux allows us to put a constraint on the foreground extinction.We assume an intrinsic Paα over Paβ ratio of 2.075 (Osterbrock & Ferland 2006), for T = 10 4 K and an electron density (N e ) of 10 2 cm −3 (for N e = 10 4 , this value would be 2.051).The observed (3 σ lower limit) Paα over Paβ ratio is 2.37, a factor 1.08 higher than the theoretical ratio, this would result in a lower limit of A V = 4.3 mag, assuming a Cardelli et al. (1989) extinction law with an R V = 4.05.Taking the 1 σ error (Sect.4.2), this lower limit in A V would increase to A V = 39.9 mag.In Sect.5.1 we compare the derived SFR to values derived via other tracers.

Multi-wavelength comparison
The Paα emission is observed over an area with a radius of ∼ 6 kpc (Fig. 3).The spatial distribution is very clumpy with most of the emission coming from the outer regions of the galaxy.In this section we compare the Paα spatial distribution with images tracing the stellar content, the cold dust, and the cold molecular gas in GN20.In Fig. 4 we compare the morphology of the Paα line emission, represented by the cyan contours to the tracers discussed above.
Comparison with the stellar distribution traced by the F560W image (panel a, Colina et al. 2023) shows that the Paα emission has a similar spatial extent as the F560W emission.The Paα clumps we detect are all originating in the outer areas of the stellar disk.Towards the bright central point sources we detect very little Paα emission.North-west of the brightest clump 4 the 3 σ contour shows an extension towards the nucleus (Fig. 3).
Panel b shows the comparison with the HST F814W image, tracing the rest-frame 0.15 µm emission.We find that there is virtually no overlap between the Paα emission and the restframe UV emission.The F814W reveals the UV-bright part of the galaxy where extinction is expected to be very low.Surprisingly, the UV emission is not associated with Paα nebular emission (see Sect. 5.3 for further discussion about the nature of the UV emission).
Comparison with the 880 µm emission shows very little spatial correlation, most of the 880 µm emission comes from the central area of the galaxy, while Paα emission is found further out.Comparison with the CO (2-1) emission line map (panel d) shows again that the CO emission is mostly coming from the central area of the galaxy where Paα is not detected.The Paα emission overlaps with the fainter outer areas of the CO emission.Similar to the Paα spatial distribution, also the CO emission shows a clumpy structure, although the clumps are located further inwards (Hodge et al. 2015).
As discussed in Sect.5.1, this multi-wavelength picture of GN20 suggests very high and strongly spatially varying extinction in GN20, where Paα and UV emission only emerge from the outskirts of the disk.Large extinction variations on (sub)kpc scales are commonly observed in local LIRGs and ULIRGs (e.g.Alonso-Herrero et al. 2006;Piqueras-López et al. 2013).

Paα clumps
The Paα emission line map shows that the Paα emission is concentrated in four clumps, annotated in Fig. 3.We extract the spectra of the clumps with an circular aperture with a radius of 2.15 pixels (FWHM of the PSF at the observed wavelength of Paα).The spectra are displayed in the left column of Fig. 5.All clumps show a spectrally resolved Paα emission profile, with clumps 1, 2 and 3 showing a double peaked profile.Clump 4 shows one single broad profile.
We calculate the SFR Paα of each of the individual clumps by integrating the spectrum over the velocity range where the clumps are detected (vertical dashed-dotted lines in Fig. 5).We find values between 48 and 70 M ⊙ yr −1 (Table 2).The total SFR Paα by summing up the values of the clumps is 225 ± 42 M ⊙ yr −1 , 55% higher than the integrated SFR Paα from the emission line map.The discrepancy between these two values could be caused by the fact that we applied an aperture correction to the clump spectra, but not to the integrated galaxy spectrum as that is not a point source and is extracted over a much larger aperture.Additionally, the clumps are located close to each other with an average separation of ∼0.5 -1.0 ′′ , making the wings of the PSF of the clumps overlap each other.What we can conclude is that the clumps dominate the total emission in the galaxy and no significant diffuse Paα emission is detected.
To derive the kinematics of the clumps we fit Gaussian profiles to the extracted spectra using the astropy modeling package.Clumps 1, 2 and 3 are fit with two gaussians simultaneously.Clump 4 shows a hint of a second component on the red side of the Paα profile, but the components are too close to allow a double Gaussian fit.We calculate the velocity dispersion of the ion-ized gas by quadratically subtracting the width of the instrumental line spread function (Labiano et al. 2021;Jones et al. 2023b).The FWHM of the instrumental broadening at 9.47 µm is 88 km s −1 (σ inst = 38 km s −1 ).From the derived velocity dispersion (Table .2), we can derive the virial mass using the isotropic virial estimator (e.g.Hodge et al. 2012): The observed velocity dispersion is represented by σ in km s −1 , the gravitational radius by R g in pc and the gravitational G = 1/232 pc (km s −1 ) 2 M −1 ⊙ is the gravitational constant.The factor C depends on the mass distribution and varies from C = 7.4 for an exponential profile to C = 4 for a Vau- couleurs profile (e.g.Bellocchi et al. 2013), introducing a factor two of uncertainty.In order to compare our results directly with the properties of the CO clumps, we follow Hodge et al. (2012) and choose C = 5 (a uniform sphere).For the gravitational radius we use the FWHM of the MRS PSF at 9.48 micron (0.367 ′′ , 2.6 kpc).Given that the Paα clumps are not resolved, the effective radii derived in this way will only be upper limits, and the derived masses are thus also upper limits.The resulting values are presented in Table 2.
The velocity profile of the clumps suggests that they are not single massive star forming clumps, but they consist of at least two sub-clumps shifted in velocity.Even clump 4 shows a hint of a second component.To study the nature of the clumps in more detail we construct channel maps of the Paα emission for each of the three double peaked clumps (middle and right columns, Fig. 5).We make channel maps for each of the clumps separating the emission of each of the components.In the top row of Fig. 5, we show first the extracted spectrum of clump 1 with the double peaked emission, the middle figure shows the channel map of the red peak between −500 and −220 km s −1 , while the right figure shows the channel map of the blue peak (−220 -−60 km s −1 ).In the figures we highlight the clump for which the channel maps In the case of clump 4 we cannot separate the two peaks and we do not show the channel map.For all three clumps the emission in the red and blue channels are spatially slightly shifted from each other.We fit a two dimensional 2D gaussian profile in Qfitsview to determine the center of the clumps.In the case of clump 1 the blue-shifted emission (clump 1b, Table 2) is shifted by 0.45 ′′ (3.2 kpc) to the north-east from the red-shifted emission (clump 1a).For clump 2, the redshifted emission (clump 2a) seems to be related to the redshifted emission in clump 1, while the blue shifted emission (clump 2b) peaks at the south-east of the extraction aperture.Clump 3a is very faint and we could not measure the position of this clump in the channel map.The bottom image in the middle column of Fig. 5 shows the Paα integrated line map (Fig. 3) with superimposed the location of the individual sub-clumps.This analysis shows that the Paα clumps are not single, very massive clumps, but more complex giant star forming regions blended by the spatial resolution of the MRS spectrograph, using both spatial and velocity information we can start separating the clumps into their different sub-clumps.
We compare the properties of these clumps to the clumps found in CO(2-1).Hodge et al. (2012) found five clumps in the CO emission, located in the central regions of GN20.The Paα clumps are located at further distances from the center (Fig. 4), apart from Paα clump 4, which spatially coincides with CO clump nr 3 (Hodge et al. 2012).
Due to the difference in spatial resolution between the CO and Paα maps it is hard to associate more clumps to each other.The kinematic properties of the Paα clumps are similar to those of the CO clumps.The derived velocities reflect the large scale rotation curve of GN20 (Sect.4.5).The measured velocity dispersion range between ∼ 40 and 138 km s −1 .The largest value (138 ± 18 km s −1 ) is measured for clump 4. Inspecting the emission line profile suggests that clump 4 also has two components, but they are too close in velocity to each other to fit them with a double gaussian.This artificially broadens the single Gaussian fit.The velocity dispersion found for the CO clumps is ∼65 km s −1 , in the same range as our Paα clumps.The CO map has a slightly higher spatial resolution than our Paα map, but also the CO clumps are spatially unresolved, resulting in upper limits for the dynamical mass.Within the limitations of both data sets the properties of both sets of clumps are similar.

Ionized gas kinematics
We use the spatially resolved Paα emission to derive the global kinematics of the ionized gas in GN20 and compare that to the kinematics of the neutral gas as traced by CO(2-1).The signalto-noise ratio of the Paα emission is too low to perform a Gaussian fit to each pixel in the data cube.We therefore use the PVdiagram (Fig. 2), to derive the rotational velocity as a function of position along the semi-major axis of GN20.
We construct the rotation curve along the semi-major axis with a position angle of 25 • by fitting a Gaussian profile to each spatial pixel of the position velocity diagram (Fig. 2).We correct the derived velocities for the inclination of GN20 of 30±15 • (using the thin-disk approximation), as derived by kinematic modelling of the CO(2-1) emission (Hodge et al. 2012), assuming the inclination of the CO is the same as Paα.The center of the rotation curve (radius = 0) is determined by collapsing the PV-diagram along the wavelength direction.This results in the integrated flux along the semi-major axis of GN20 and reveals the bright nucleus as seen in the F560W image.We fit a Gaussian profile to the emission peak to determine the center of the rotation curve.
The resulting rotation curve is shown in Fig. 6.We measure a velocity profile, increasing with radius to +300 km s −1 and -400 km s −1 .The velocity dispersion varies between 85 and 190 km s −1 , with a typical error of 40 km s −1 .The spatial resolution of the Paα observation is 2.6 kpc.This means that beam smearing strongly affects the shape of the observed velocity profile as well as the velocity dispersion.The gradual change of the rotational velocity in the central two -three kpc remains unresolved and cannot be distinguished from e.g. a step function.The velocity dispersion is affected by overlapping clumps, as can be seen in the PV-diagram (Fig. 2), where at ∼ -0.7 ′′ two clumps (clump 3 and 4) create an artificially broad Paα profile.
As comparison we constructed the rotation velocity profile from the CO (2-1) data.We use the derived velocity map (Fig. 3 in Hodge et al. 2012) to measure the rotational velocity.First we correct for the difference in spatial resolution between the CO(2-1) maps (0.19 ′′ ) and the Paα map (0.365 ′′ ).We convolve the CO(2-1) velocity map with a Gaussian kernel with σ = 0.32 ′′ (the quadratic difference between the two FWHM values).Using the astrometric information, we regrid the convolved velocity map to the same pixel scale as the Paα linemap to compare them on the same spatial resolution.We collapse the velocity map along the semi-major axis and use the maximum observed velocity at each distance as the radial velocity measurement, the standard deviation of the velocity values is used as estimate of the uncertainty.The resulting CO(2-1) velocity curve is plotted in red in Fig. 6.
Both the CO(2-1) and Paα rotation curves reach the same minimum and maximum velocity.However, they also show differences.The Paα rotation curve is much steeper than the CO(2-1) rotation curve, and both at -4 kpc as well as 1 kpc, the two rotation curves differ by ∼ 200 km s −1 .This could be related to the different spatial distribution of the Paα emission with respect to the CO(2-1) emission (Fig. 4.) The Paα clumps are situated around the nucleus and the Paα rotation curve is dominated by a rotating "ring" of clumps.The CO distribution is much smoother and stronger in the center and emitting more light at lower rotational velocities.
From the rotation curve of Paα we find a v rot = 1 2 (v max − v min ) = 550 ± 40 km s −1 .This value agrees with the dynamical modeling of the CO(2-1) emission (575 ± 100 km s −1 Hodge et al. 2012).From the measured velocity dispersion we calculate the flux weighted velocity dispersion, a measure of the random gas motions in the galaxy (e.g.Glazebrook 2013): σ m = 145 ± 53 km s −1 .This value should be interpreted as an upper limit to the real flux weighted velocity dispersion as we integrate the emission over the minor axis and overlapping clumps will create an artificial broadening of the emission line profile.From these quantities we calculate v rot /σ m = 3.8 ± 1.4.This shows that the overall gas kinematics of GN20 are rotationally dominated in contrast to many (high-z) starburst galaxies (e.g.Förster-Schreiber et al. 2009;Law et al. 2009;Bik et al. 2022), which have v rot /σ m values below the dividing line between dispersion and rotationally dominated of 1.83 (Förster-Schreiber & Wuyts 2020).Rizzo et al. (2021) present a kinematic study of 5 strongly lensed DSFGs around z=4.5 with ALMA and found high v rot /σ m values for all galaxies, suggesting that they are rotationally dominated disk galaxies.Comparing GN20 values with those of (U)LIRGs shows that the v rot /σ m of GN20 is consistent with that of low-z LIRGs.Bellocchi et al. (2013) derive an average v rot /σ m = 3.4 for LIRGs and v rot /σ m = 2.0 for ULIRGs.Bellocchi et al. (2013) divide their (U)LIRG sample in different morphological types corresponding to the different phases along the merging process.The measured v rot /σ m makes GN20 consistent with interacting systems with a perturbed disk.
Considering that the Paα and CO(2-1) kinematics as well as the morphology suggest that GN20 is a large disk galaxy, we compare the derived value for σ m with theoretical predictions for σ m in the context of a galactic disk model.Krumholz et al. (2018) shows that σ m in disk galaxies is determined by both star formation feedback processes and turbulence due to gravitational instabilities.Comparing the location of GN20 in the σ m vs SFR diagram (Fig. 4 in Krumholz et al. 2018) show that the high σ m values for GN20 derive from gravity driven turbulence and shares the location with ULIRGs in the more nearby universe (see also Arribas et al. 2014).Interestingly, these ULIRGs are not disk galaxies, but galaxies undergoing a major merger event.However, our derived σ m should be seen as an upper limit due to the overlapping clumps.Additionally, GN20 is significantly more massive and likely has a much more gas rich disk than ULIRG galaxies, resulting in higher σ m due to the higher mass of the clumps.

Discussion
Using deep MIRI/MRS spectroscopy, we have detected spatially resolved Paα emission in GN20.In this section we discuss the implications of our findings and compare them to auxiliary data.

Extinction in GN20
We find a large discrepancy between the SFR derived from the Paα emission and the SFR derived from the infrared luminosity.Based on the Paα emission, under the assumption of no extinction, we derive SFR Paα = 144 ± 9 M ⊙ yr −1 (Table 1).As shown in Sect.4.2, the observed 3σ upperlimit of Paβ, does not provide a strict constraint on the extiction.Tan et al. (2014) derived the infrared star formation rate (SFR IR ) based on 1.2-3.3mm continuum observations, assuming a Chabrier (2003) IMF, and find a value for GN20 of SFR IR = 1860 ± 90 M ⊙ yr −1 .This means that the observed Paα emission only reveals 7.7 ± 0.5 % of the total star formation rate as detected from the mm observations.By comparing the SFR Paα to the SFR IR for a sample of local LIRG galaxies, Tateuchi et al. (2015) find that, after extinction correction, both SFR measurements agree within a scatter of 0.27 dex.Piqueras-López et al. ( 2016) study the same relation, but with a sample of ULIRGs included.They find that in the case of the ULIRGs the SFR Paα corrected for extinction is typically lower than the SFR IR .Giménez-Arteaga et al. ( 2022) find again a better correlation between the SFR derived from Paβ with the SFR IR .
The disagreement between SFR Paα and SFR IR can be caused by (a combination of) reasons: (i) The mm emission is not (fully) caused by star formation, but has a contribution from an active galactic nucleus (AGN), (ii) Ionizing photons can be directly absorbed by dust, contributing to SFR IR but not to SFR Paα (Piqueras-López et al. 2016), (iii) SFR IR traces the star formation history over a longer time (up to 100 Myr) than Paα (10 Myr) (Kennicutt & Evans 2012), and (iv) The extinctions are so large that even Paα does not trace the full amount of ionized gas (Álvarez Márquez et al. 2023b).
A significant fraction of SMGs are major mergers (e.g.Gillman et al. 2023), as their local ULIRG counterparts, which makes the presence of an AGN in the center a realistic possibility.GN20 is located in an overdensity and shows evidence for a double nucleus, suggesting that GN20 is interacting or in an advanced stage of a major merger (Colina et al. 2023).Therefore the presence of an AGN can be expected (Ricci et al. 2017;Blecha et al. 2018).Based on the detection of the PAH 6.2 µm emission, Riechers et al. (2014) conclude that the infrared emission from the nucleus could contain a significant contribution of a buried AGN.The total observed infrared SED is, however, dominated by the obscured starburst.Also Colina et al. (2023) attribute the bright nuclear emission seen in the F560W imaging to a strong nuclear starburst.
The center of GN20 contains a large amount of dust and gas, therefore a fraction of the ionizing photons will be absorbed by dust directly, before ionizing the gas.Additionally, SFR Paα traces the most recent star formation over the last 10 Myr, while SFR IR is sensitive to star formation over a much longer time period as intermediate mass stars are also capable of heating the dust, but do not produce ionizing photons in sufficiently large quantities.Piqueras-López et al. (2016) attribute the observed difference between SFR Paα and SFR IR to both these effects.The difference between SFR Paα and SFR IR in GN20 is, however, much larger than observed by Piqueras-López et al. (2016).Hodge et al. (2015) shows that the obscured, strong, starburst in the center of GN20 has a gas depletion time of ∼ 100 Myr.If no Paα would be emitted from such a starburst, that would mean that the extreme starburst in the center suddenly stopped forming stars >10 Myr ago, while there is still a huge reservoir of molecular material present.
Therefore, we attribute the discrepancy to the fact that we do not trace all the star formation with Paα due to the high extinction.As we do not detect any other hydrogen recombination line, we cannot derive the extinction for Paα.We can, however, derive an estimate from the comparison between SFR Paα and SFR IR .Using the Cardelli et al. (1989) extinction law with a R V = 4.05, we find a total extinction of A V = 17.2 ± 0.4 mag.As a more realistic approach we use an extinction model assuming the gas and stars are mixed.We use the relation derived by Calabrò et al. (2018) for a sample of dusty starburst galaxies at z=0.5-0.9 in combination with the extinction coefficients from Cardelli et al. (1989) extinction law to get an estimate of the attenuation.Using equation 1 in Calabrò et al. (2018), our ratio of observed SFRs gives a total visual attenuation of A V,mixed = 44 ± 3 magnitudes.
The multi-wavelength comparison shown in Fig. 4 shows that the mm continuum emission is concentrated in the center of GN20, while Paα is originating from the clumps at a few kpc distance from the center.Without extinction correction these clumps account for 7.7% of the total star formation.An extinction corrected SFR Paα of the clumps will increase this fraction, for example increasing this fraction to 20% would require an A V = 7 mag for the Paα clumps.
We attribute the remaining fraction of the star formation to be in the center, which would suggest an even higher extinction in the center of GN20.Piqueras-López et al. (2013) show that most (U)LIRGs display a similar behaviour where large extinction variations, with areas with A V as high as 20-30 magnitudes, are found (see also Bohn et al. 2023).These (U)LIRGs typically have an obscured center.However, in most local (U)LIRGs the extinction corrected SFR traced by the Paschen lines follows pretty well the SFR IR (Giménez-Arteaga et al. 2022), suggesting that extinctions in most (U)LIRGs are not as extreme as in GN20.Two ULIRGs show very high extinction in their nuclei.Based on near-infrared extinction measurements, Engel et al. (2011) derive very high extinction in the two nuclei of Arp 220.They show that some of the gas is so obscured that it is not detected in the near-infrared emission lines.Based on MIR diagnostics, Haas et al. (2001) derive values of A V,mixed = 110 mag for the central starburst in Arp 220 and A V,mixed = 50 mag in the center of UGC5105 (see also Sturm et al. 1996).This shows that such high extinction as seen in GN20 is also detected in the more extreme ULIRG galaxies, where a very dusty extreme starburst is driven by the merger process.
The high extinction measured towards GN20 will have strong implication on stellar mass determinations from optical and near-infrared rest frame images.Colina et al. (2023) derive a disk mass of 5.2×10 10 M ⊙ (assuming zero extinction), a factor of 2 lower than the stellar mass derived from SED fitting (Tan et al. 2014).To quantify the effect of the spatially variable extinction on the stellar mass determination of GN20 is beyond the scope of this paper.The extinction most likely affects the youngest populations the strongest, while the older population is carrying the mass of the galaxy.A detailed multi-wavelength analysis utilizing HST, future NIRCAM and MIRI imaging covering the UV to near-IR rest-frame will address this issue.
As shown in Fig. 4 the different tracers in GN20 show a very different spatial morphology.The high extinction makes the galaxy almost invisible and results in a significant suppression of the Paα emission from the central starburst.This would bias the measured half-light radii to larger values, affecting the mass-size relation of starburst galaxies (e.g. Ward et al. 2023;Costantin et al. 2023;Ormerod et al. 2023).Observations of local LIRGs show that the galaxies are much more compact in CO than in ionized gas or stars (Bellocchi et al. 2022).Additionally, Bellocchi et al. (2022) compare the effective radii from Hα and Paα and find that the Hα radii are significantly larger, suggesting that Hα is more affected by extinction in the center of the LIRGS.Fujimoto et al. (2017) find similar trends based on HST and ALMA observations of z = 1 -5 starburst galaxies.Popping et al. (2022) demonstrate that galaxies in the TNG50 cosmological simulation have larger half light radii in the optical and near-infrared light compared to the dust emission and attribute this difference to extinction.
High spatial resolution observations of U(LIRGs) have revealed extremely compact (<100 pc) obscured nuclei (CONs) in many ULIRGs and LIRGs (e.g. Aalto et al. 2015;Donnan et al. 2023).Deeply buried in these CONs is either a rapidly growing super massive black hole or a compact nuclear starburst, hidden behind very high column densities of dust and gas.In GN20 we do not have the spatial resolution to find such a compact region, however, the highly obscured and extreme nuclear starburst we find in GN20 might resemble what we observe as CONs in nearby (U)LIRGs.This is consistent with the findings of Cortzen et al. (2020) who show that the emission longward of 170 µm is optically thick, consistent with a very compact starburst in the nucleus of GN20.Similar properties are observed for other highredshift starbursts (Jin et al. 2022).

Nature of the star forming clumps
We identify four individual clumps based on the Paα emission line map, using their kinematics we find that they are not single clumps, but consist of at least two sub-clumps in each of them.We search for continuum emission from these clumps in the F560W image (Colina et al. 2023).The F560W image has a spatial resolution of 0.24 ′′ (1.7 kpc), slightly higher than that of the Paα linemap.
Following Colina et al. (2023), we apply the Lucy-Richardson deconvolution algorithm (Lucy 1974) to the image data.For the deconvolution of the input data, we use the empirical PSF derived by Colina et al. (2023) which is based on close-by stars in the MIRI imaging field of view.Due to the confusion-free detection of GN20 in the MIRI data, we do not apply any background subtraction.In total, we use 10000 iteration steps to minimize the chance of artifacts.In the final step, we convolve the resulting delta map with a 3-pixel Gaussian Kernel filter adapting the procedure described in Peißker et al. (2022).We find that some of the Paα clumps can be attributed to clumpy emission in the F560W image, clumps 2b, 3a and 3b show a matching F560W clump.Based on their location and kinematics we suggest in Sect.4.4 that clump 1a and 2a might be tracing the same clump.In the F560W image, there is a clumpy structure located north-east of clump 2a, making this a potential counterpart.For clumps 1a and 4 on the other hand, no clear counterpart in the F560W can be found.Future higher spatial resolution NIRCAM imaging of GN20 will reveal more details on the clumpy structure of the GN20 disk.
Clumps are commonly observed in star forming galaxies at z ≥ 1 (e.g.Elmegreen et al. 2009b;Schreiber et al. 2011;Livermore et al. 2012;Adamo et al. 2013;Dessauges-Zavadsky & Adamo 2018).The clumpiness of high-z galaxies is attributed to the higher gas fraction in galaxy disks (e.g Genzel et al. 2008;Girard et al. 2018).The disks are fed by infall and accretion of gas-rich material, where clumps form via violent disk instabilities (Dekel et al. 2009b;Ceverino et al. 2010;Mandelker et al. 2014;Inoue et al. 2015).In this scenario, the giant clumps that form have masses of a few percent of the mass of the disk (Ceverino et al. 2010).Hodge et al. (2012) derive a dynamical mass for GN20 of 5.4 ± 2.4 × 10 11 M ⊙ .The upper limits we derive on the clump masses (∼ 10 9 M ⊙ , Table 2) are consistent with the predictions of the violent disk instabilities scenario.The masses derived for the CO clumps by Hodge et al. (2012) are consistent with this as well.Similarly, Spilker et al. (2022) measure high velocity dispersion in the clumps detected in [CII] in the merging DSFG SPT0311-58.This also suggests very massive clumps under the assumption that the measured velocity dispersion is driven by gravitational motions.Higher spatial resolution observations might well break the observed clumps down in several sub-clumps (e.g.Messa et al. 2022;Meštrić et al. 2022), making the measured velocity dispersion only upper limits.
Additionally, merging galaxies trigger strong starburst resulting in the formation of clumps.

Nature of the UV-bright region
The HST images of GN20 (Fig. 4) reveal a UV-bright elongated structure west of the center of the galaxy.It is detected in HST images, from F606W (rest-frame 0.12µm) to F160W (rest-frame 0.32µm).Also in the MIRI F560W image the arm is visible (Colina et al. 2023).The convolved F560W image reveals 4-5 clumps associated with the UV-bright emission.This arm is not detected in CO, 880µm and Paα.The fact that this region is detected in the UV would suggest low extinction.However, UV-bright would also imply a young stellar population, and therefore we would expect Paα from this low extinction, young region in GN20.
Based on fitting of the observed UV to the radio spectral energy distribution (SED) of GN20, Tan et al. (2014) find an SFR IR /SFR UV of ∼13, this would lead to an SFR UV of 140 M ⊙ yr −1 .However, the UV emission only traces a small piece of the galaxy, therefore we reanalyze the HST data, by measuring the integrated flux with an elliptical aperture of 1 arcsec 2 in size, matching the size and morphology of the UV-bright emission (1.6 ′′ ×0.8 ′′ , PA: 0 • ).
Under the assumption of no extinction we calculate the SFR UV from the measured F814W flux using the relation detailed by Murphy et al. (2011) and find SFR UV = 15 ± 2 M ⊙ yr −1 .As a next step we perform a fit using CIGALE (Burgarella et al. 2005;Noll et al. 2009;Boquien et al. 2019)  Based on the standard deviation derived for the Paα image (Sect.4.2) we derive a 3 σ detection limit of 2.5 M ⊙ yr −1 per pixel.Assuming a constant spatial distribution and an MRS pixel size of 0.029 arcsec 2 gives a 3σ upper limit of 15 M ⊙ yr −1 for a 1 arcsec 2 aperture.This Paα upper limit is consistent with the directly measured SFR UV , but only under the assumption of no extinction.The CIGALE fits show a significant extinction of A V = 0.8 magnitudes, increasing the SFR UV significantly.Under the assumption of the A V = 0.8, the UV derived star formation rate would increase to 43 M ⊙ yr −1 .
In summary, Paα traces the star formation over the last 10 Myr in GN20, while the star formation traced by the UV emission is sensitive to the last 100 Myr (Murphy et al. 2011;Kennicutt & Evans 2012).The derived SFR values are consistent with a stellar population older than 10 Myr and younger than 100 Myr, such that it results in a bright UV region without strong ionized gas emission.This is comparable with what is observed in the local LIRG galaxy NGC 1640 (Adamo et al. 2020), where very little Hα emission is detected towards a UV-bright arm.Adamo et al. (2020) find a cluster population in this arm of ∼ 20-30 Myr.This population still contributes strongly to the UV flux, but is too old for significant production of ionizing photons.

What triggered the starburst in GN20?
Two main scenarios been proposed in the literature for the presence of starbursts in DSFG.The first scenario is that of a major merger, similar to local ULIRGS, while the second scenario is violent disk instabilities in gas rich galaxy disks.In this section we summarize the properties of GN20 derived in this and previous papers and speculate about the origin of the starburst in GN20.
The offset between the nucleus and the center of the outer isophotes of the disk as well as the secondary nucleus found in GN20 by Colina et al. (2023) would suggest that GN20 is at the final stage of a major merger.The extremely obscured nucleus (Sect.5.1) and the very efficient star formation (Hodge et al. 2015) make GN20 display similar properties to the most extreme ULIRG galaxies.In this scenario, the star forming clumps would then be a byproduct of the merger process as seen in local merging (U)LIRGs (e.g.Pereira-Santaella et al. 2016;Piqueras-López et al. 2016;Adamo et al. 2020).Additionally, GN20 is known to be located in an overdensity (Daddi et al. 2009), making it likely that the star formation is at least affected by gravitational interaction.
The kinematics of both the CO and Paα are consistent with a rotationally supported, gas-rich disk.Additionally, our Paα observations and the CO(2-1) observations of Hodge et al. (2012) show that the star formation happens over a much larger spatial extent than in ULIRG galaxies (Hodge et al. 2016;Bellocchi et al. 2022), and covers the entire observed disk of GN20.These observations would hint that GN20 is a disk galaxy whose star formation is fuelled by accretion of neutral gas from intergalactic matter, where the clumps would form in the gas rich disk by violent disk instabilities (Dekel et al. 2009b;Ceverino et al. 2010;Inoue et al. 2015).Hung et al. (2015) study how the kinematics of a sample of (U)LIRGS would look like when redshifted to z=1.5 with a spatial resolution of 900 pc.They show that the kinematics of galaxies in a late stage of the merger process are indistinguishable from a disk-like rotational profile.On the other hand Ueda et al. (2014) studied 24 local (U)LIRGS and find that extended molecular disks are common in the late stages of mergers.This is consistent with major merger simulations of Springel & Hern-quist (2005), predicting the formation of extended gas rich disks at the end of the merging process.
Beam smearing, due to the spatial resolution of the MRS, prevents us from deriving a well sampled rotation curve from the Paα kinematics, and makes it impossible to determine whether the kinematics are determined by a real gas rich galaxy disk, possibly rebuilt after the major merger, or a late-stage merger mimicking a disk.In the case of a rebuilt gas disk during the merger processes (Ueda et al. 2014), the Paα clumps could possibly be formed in the disk due to violent disk instabilities, while the central starburst is the result of the actual merger.

Conclusions
We present deep MIRI/MRS spectroscopy of Paα and Paβ of the dusty star forming galaxy GN20 at z=4.0548.We detect for the first time ionized gas Paα emission out to a radius of 6 kpc, distributed in a clumpy morphology.The Paβ line is not detected.From the integrated Paα flux we derive a SFR Paα = 144 ± 9 M ⊙ yr −1 , assuming no extinction.The SFR Paα is only a small fraction of the SFR IR (7.7 ± 0.5 %, Tan et al. 2014).We attribute this difference as due to the high extinction in GN20, especially in the central starburst and find an average extinction of A V = 17.2 ± 0.4 mag.Such high values are also measured in the nearby ULIRGs Arp220 and UGC 5105.
We detect four spatially unresolved clumps in the Paα emission line map.By studying their Paα emission profile we find that each clump consists of at least two sub-clumps.These subclumps are also identified in the convolved F560W image.Their velocity dispersions are in the same range as measured for the CO clumps (Hodge et al. 2012).The CO clumps are located more in the central region of GN20.
We construct a Paα rotational profile and find it broadly consistent with the CO kinematics (Hodge et al. 2012).Within the spatial resolution of the Paα observations (2.5 kpc) the kinematics is consistent with that of a rotationally dominated disk galaxy.This suggests that star formation is driven by gravitational instabilities in the gas rich galactic disk.The observed differences between the rotation curves could be related to the differences in spatial distribution between the CO and Paα emission.
We do not detect any Paα emission towards the UV-bright emission of GN20 as seen in HST imaging.Based on the SED fitting and the Paα upper limit, we conclude that this region is older than 10 Myr, where the stellar population is still bright in the UV but does not produce significant LyC radiation anymore.
Finally, we speculate on the nature of the starburst in GN20.Evidence suggests that GN20 is in the late stage of a major merger.Based on the Paα kinematics we cannot distinguish between a late stage merger mimicking the disk-like rotation profile and a large extended disk, possibly re-created toward the end of the merger process.

Fig. 1 .
Fig. 1.Integrated spectrum of GN20 centered on Paα.Plotted in green is the standard deviation derived from a background aperture with the same size as the galaxy spectrum.The spectrum plotted in red color shows the CO(2-1) emission extracted from the data presented in Hodge et al. (2012), overplotted in blue is the single Gaussian fit to derive the redshift (see text).The systemic velocity corresponds to a redshift of z = 4.0548(Carilli et al. 2011).

Fig. 2 .
Fig. 2. Position-velocity diagram of Paα along the major axis of the galaxy (PA = 25 • ).The offset is relative to the position of the bright nucleus of GN20 (Colina et al. 2023).The red solid lines show the extraction apertures used for creating the Paα emission line map, while the red dashed lines shows the corresponding background map from which we derive the error on the emission line map.

Fig. 3 .
Fig. 3. Paα emission line map after extraction based on the PV diagram (Fig.2, see text).The 2,3 and 4-σ contours are shown in cyan.The black contours show the bright nucleus of GN20 detected in the F560W image (Colina et al. 2023).The MRS point spread function at 9.47µm is displayed as a grey circle.The PSF is slightly asymmetric and the size shows the FWHM in the α (x)-and β (y)-direction, rotated according to the observed position angle.The 4 clumps discussed in Sect.4.4 are annotated.

Fig. 5 .
Fig. 5. Left column: Integrated spectra of the clumps identified in the Paα line map.The grey histograms represent the 1σ errors.The vertical dotted lines represent the borders in which the flux and integrated star formation rate is calculated (Tab.2).The green histograms show the best gaussian fits to the emission profile, with the dotted line the individual gaussians for clumps 1,2 and 3. Middle and right columns: Paα channel map for each clump annotated by the blue circle.The two channel maps for each clump show the location of the two different peaks in the Paα spectra.Clump 4 is fitted by 1 Gaussian only and no channel maps are shown.In the emission channel map related to the red peak of clump 3 (row 3, right channel map), emission is visible at the location of clump 4 due to the similar velocity range.The peak at the location of clump 1 is caused by a (2-σ) peak at +125km s −1 in the spectrum of clump1.The bottom image in the middle column displays the integrated Paα linemap (Fig. 3) with the location of the sub-clumps overlayed as blue asterisks.
ALMA and HST observations of the nearby LIRG IC 4687 shows a similar pattern, where a large fraction of the clumps is detected in only one of the two tracers (Paα or CO(2-1), Pereira-Santaella et al. 2016).Studying a larger sample of LIRGS, Sánchez-García et al. (2022) reached similar conclusions.

Fig. 6 .
Fig. 6.Rotational velocity of GN20 derived from a single gaussian fit to the Paα PV diagram (Fig. 2) along the semi-major axis compared to the rotation profile from CO(2-1) derived from the velocity map in Hodge et al. (2012).The blue shaded area shows the range affected by beam smearing.

Fig. 7 .
Fig. 7. Convolved image of GN20 in F560W.The white contours stand for the contours of the Paα emission line map (Fig. 3).The blue asterisks show the position of the clumps identified in Sect.4.4.The red circles show the location of the CO clumps identified by Hodge et al. (2012).The beam size of the CO observations is shown in the bottom left and the PSF of the MIRI F560W image on the bottom right.
Figure 7 shows the resulting image with the Paα emission map contours in white.The blue asterisks show the different Paα clumps identified in Sect.4.4 and the red circles show the position of the CO clumps identified in Hodge et al. (2012).The final convolved F560W image shows a very clumpy structure, similar to the residual image (Fig. 1 in Colina et al. 2023), after removing the best fitted two component model, representing the smooth light distribution in GN20.Considering the clumpy nature of both the CO and the Paα, both tracing current star formation, the diffuse component of the F560W image could be related to more mature stellar populations having a smoother spatial distribution.
(U)LIRGs in the local universe show clumpy star formation.Piqueras-López et al. (2016) finds clumps in Paα in local (U)LIRGs of sizes between 300 pc and 1500 pc in the ULIRGs and a bit smaller in the LIRGs, while Larson et al. (2020) find Paα clumps in local LIRGs and find sizes between 90 and 900 pc.
to the integrated UVoptical SED as seen by HST.We use a Chabrier (2003) IMF with a Calzetti et al. (2000) extinction law and model the SED with two stellar populations; a young population (< 40 Myr) with constant star formation history and a older single stellar population between 100 -500 Myr.For the young stellar population (relevant for the UV emission) we find a SFR = 135 ± 70 M ⊙ yr −1 and an extinction of A V = 0.8 magnitudes.Next, we measure the flux of the UV-bright region in the MIRI F560W (Colina et al. 2023) and F770W (Crespo Gómez et al, in prep) with the same aperture as for the HST images.We measure the flux in the images with the 2 component Lenstronomy (Birrer & Amara 2018) fit removed (see Colina et al. 2023, for details).Repeating the CIGALE fit with the two extra MIRI points gives a similar SFR = 160 ± 100 M ⊙ yr −1 for the young population, with a similar extinction as the previous fit.

Table 1 .
Integrated flux measurements