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

A Bik, J Álvarez-Márquez, L Colina, A. Crespo Gómez, F Peissker, F Walter, L.A Boogaard, G Östlin, T.R Greve, G Wright, et al.

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

► To cite this version:

A Bik, J Álvarez-Márquez, L Colina, A. Crespo Gómez, F Peissker, et al. Clumpy star formation and an obscured nuclear starburst in the luminous dusty z=4 galaxy GN20 seen by MIRI/JWST. Astronomy & Astrophysics - A&A, 2024, 686, pp.A3. 10.1051/0004-6361/202348845. hal-04356069

HAL Id: hal-04356069

https://hal.science/hal-04356069v1

Submitted on 20 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Clumpy star formation and an obscured nuclear starburst in the luminous dusty z = 4 galaxy GN20 seen by MIRI/JWST

A. Bik¹, J. Álvarez-Márquez², L. Colina², A. Crespo Gómez², F. Peissker³, F. Walter⁴, L. A. Boogaard⁴, G. Östlin¹, T. R. Greve⁵,⁶,⁷, G. Wright⁸, A. Alonso-Herrero⁹, K. I. Caputi¹⁰,⁶, L. Costantin², A. Eckart³, S. Gillman⁵,⁶, J. Hjorth¹¹, E. Iani¹⁰, I. Jermann⁵,⁶, A. Labiano¹²,², D. Langeroodi¹¹, J. Melinder¹, P. G. Pérez-González², J. P. Pye¹³, P. Rinaldi¹⁰, T. Tikkanen¹³, P. van der Werf¹⁴, M. Güdel¹⁵,⁴,¹⁶, Th. Henning⁴, P. O. Lagage¹⁷, T. Ray¹⁸, and E. F. van Dishoeck¹⁴

*(Affiliations can be found after the references)*

Received 5 December 2023 / Accepted 2 March 2024

## ABSTRACT

Dusty star-forming galaxies emit most of their light at far-infrared to millimeter wavelengths as their star formation is highly obscured. Far-infrared and millimeter 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_\odot$ yr⁻¹) is only 7.7 ± 0.5% of the infrared star formation rate (1860 ± 90 $M_\odot$ yr⁻¹). We attribute this to very high extinction ($A_V$ = 17.2 ± 0.4 mag, or $A_{V,\text{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 consists of at least two sub-clumps. We find mass upper limits consistent with them being formed in a gravitationally unstable gaseous disk. The ultraviolet bright region of GN20 does not have any detected Paα emission, suggesting an age of more than 10 Myr for this region of the galaxy. From the rotation profile of Paα, we conclude that the gas kinematics are rotationally dominated and the $v_{\text{rot}}/\sigma_m$ = 3.8 ± 1.4 is similar to low-redshift luminous infrared galaxies. From the Paα kinematics, we cannot distinguish between a rotational profile of a large disk and a late stage merger mimicking a disk. We speculate that GN20 is in the late stage of a major merger, where the clumps in a large gas-rich disk are created by the major merger, while the central starburst is driven by the merger event.

**Key words.** galaxies: high-redshift – galaxies: ISM – galaxies: individual – galaxies: kinematics and dynamics – galaxies: starburst

## 1. 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_\odot$ yr⁻¹ (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 that attenuates the light. They are identified as submillimeter galaxies (SMGs) in deep submillimeter 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, 2022; Finkelstein et al. 2023; Pérez-González et al. 2023a). More recent surveys with the Atacama Large Millimeter/submillimeter Array (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. 2023). Several of these galaxies are located in overdensities or are members of protoclusters, 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 analogs of ultra-luminous infrared galaxies (ULIRGs) in the local Universe.

The other scenario is driven by the observations of star-forming galaxies at high redshift showing gas-rich disks (e.g., Genzel et al. 2008; Elmegreen et al. 2009a). 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 the 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 kilo-parsec sizes and masses of ~$10^9 M_\odot$ (Genzel et al. 2008; Elmegreen et al. 2009a; 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 million years (e.g., Mandelker et al. 2014) and contribute
to the formation of the thick disk and the bulge growth (Elmegreen et al. 2008, 2009b; Ceverino et al. 2010). Observations of individual galaxies support this scenario by finding trends between the galactocentric distance and age of the clumps (Adamo et al. 2013; Cava et al. 2018).

Due to the fact that DSFGs have high extinction, almost all we know about these galaxies is derived from their molecular gas and dust content (e.g., Hodge et al. 2012, 2015) tracing mostly the neutral interstellar medium. [CII] observations of DSFGs trace both the photodissociation regions (PDRs) 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 wavelengths probed by the Hubble Space Telescope (HST), or only reveal very faint emission. Only now with the arrival of JWST we are able to reveal the rest-frame 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. 2024; Pérez-González et al. 2023b; 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. 2023; Arribas et al. 2023; Jones et al. 2024; Parlanti et al. 2024).

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 Submillimetre Common-User Bolometer Array (SCUBA) imaging (Pope et al. 2006). GN20 is a DSFG at redshift 4.0548 (Carilli et al. 2011), located in a protocluster or galaxy overdensity (Daddi et al. 2009) and has an infrared luminosity of $1.86 \times 10^{13}$ L<sub>☉</sub>. 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<sub>☉</sub> yr<sup>-1</sup> (Tan et al. 2014). With a total stellar mass of $1.1 \times 10^{11}$ M<sub>☉</sub> (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<sup>-1</sup>, 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 \pm 2.4 \times 10^{11}$ M<sub>☉</sub>), extending as far as 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-centered 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α ($\lambda_{rest} = 1.8756$ µm) and Paβ ($\lambda_{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 Sects. 5 and 6 respectively.

In this paper we adopt the cosmological parameters from the Planck Collaboration XIII (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.

# 2. Observations and data reduction

## 2.1. MRS observations

GN20 was observed by JWST on November 24, 2022, using the Mid-Infrared Instrument (MIRI, Rieke et al. 2015; Wright et al. 2015, 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 MRS observations were obtained in the MEDIUM configuration. The MEDIUM configuration provides observations in four different wavelength intervals: 5.66–6.63 µm, 8.67–10.13 µm, 13.34–15.57 µm, and 20.69–24.48 µm (Labiano et al. 2021; Argyriou et al. 2023). This paper describes the analysis of MRS channels 1 and 2 in the MEDIUM configuration, our observations do not have enough sensitivity to detect any emission lines in channels 3 and 4, these bands are not discussed further. The observations mainly cover the Paα emission line ($\lambda_{rest} = 1.87561$ µm), but channels 1 and 2 also cover additional ionized gas tracers as; Paβ ($\lambda_{rest} = 1.282$ µm) and [SiVI] ($\lambda_{rest} = 1.963$ µm) as well as two H₂ lines tracing molecular gas, namely H₂ (1–0) S(7) ($\lambda_{rest} = 1.747$ µm) and H₂ (1–0) S(3) ($\lambda_{rest} = 1.957$ µm).

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 s per integration, summing up to a total of 7645 s. 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.

## 2.2. Calibration

The MRS data were calibrated using the developers' version of the JWST calibration pipeline, version 1.11.1.dev3+g04332ea2 (Bushouse et al. 2023) with context 1097 of the Calibration Reference Data System (CRDS). We follow the same data reduction strategy as described in Álvarez Márquez et al. (2023, 2024), Bosman et al. (2023).

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 minimize the effect of the CR showers. We turned on the find_shower keyword and, following Bosman et al. (2023) set the time_masked_after_shower and max_extend radius to 60 s 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. 2024; 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 noise in the datacube to that created from the background with all the frames. The lowest rms noise 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. 2023).

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 that are not removed by the pipeline.

## 2.3. Archival datasets

In this paper we compare the MRS observations with a multi-wavelength 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" × 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), astrometrically aligned to Gaia DR3 (Gaia Collaboration 2023).

## 2.4. 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.pmap following the procedure presented in Álvarez Márquez et al. (2023). The reduced F770W image is astrometrically calibrated using Gaia DR3 (Gaia Collaboration 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 880 mm and CO(2–1) observations were phase referenced to quasars and their absolute astrometry is estimated to be a tenth of their synthesized beam, that is 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.

# 3. 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.

## 3.1. 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
<figure><img src="image_3.png" /></figure>

**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).

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

**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 show the corresponding background map from which we derive the error on the emission line map.

spectra look fairly similar, with the exception that the Paα emission shows blue shifted emission to −500 km s⁻¹, while the CO emission is only detected to −400 km s⁻¹.

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<sub>Paα</sub> = 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. (2024). 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.

**3.2. Construction of emission line maps**

The integrated Paα spectrum (Fig. 1) shows emission at velocities between −450 and +400 km s⁻¹. 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α emis-

sion 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⁻¹) 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₂ and [SVI] lines are also located in channel 2. In order to extract emission line maps, we shift the map with the central wavelength derived for Paα to the central wavelength of the line of interest. For the Paβ aperture in channel 1, we scale the extraction aperture taking into account the different pixel size and spectral resolution between channel 1 and 2.

**4. Results**

**4.1. Emission line maps**

The extracted Paα emission line map is shown in Fig. 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 four clumps detected in the outer regions of the disk and very little emission toward 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⁻¹⁹ erg s⁻¹ cm⁻² 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⁻¹⁸ erg s⁻¹ cm⁻². Figure 3 shows the 2, 3 and 4σ contours of Paα, derived using this standard deviation. For Paβ, we derive
<figure><img src="image_5.png" /><figcaption>**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 gray 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 four clumps discussed in Sect. 4.4 are annotated.</figcaption></figure>

<table><caption>**Table 1.** Integrated flux measurements.</caption><thead><tr><th>Line</th><th>Flux<br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>SFR<sub>Av=0</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Paα</td><td>1.9 ± 0.11 × 10<sup>-17</sup></td><td>144 ± 9</td></tr><tr><td>Paβ</td><td>&lt;2.6 × 10<sup>-18</sup> (a)</td><td>&lt;43 (a)</td></tr></tbody></table>

Notes. (a) 1σ upper limit.

a standard deviation of $1.86 \times 10^{-19}$ erg s$^{-1}$ cm$^{-2}$ per pixel, translating to a $3\sigma$ detection limit in a aperture equal to the FWHM (0.3") of the PSF of $2.3 \times 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 \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$ for [SiVI] and $1.4 \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$ for both the H$_2$ lines.

## 4.2. 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 \pm 0.11 \times 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 \pm 9$ M$_{☉}$ yr$^{-1}$, assuming zero extinction.

To derive an upper limit 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\sigma$ upper limit for the integrated flux of Paβ of $3 \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$ (equivalent to an SFR$_{Paβ} < 43$ M$_{☉}$ yr$^{-1}$ in the case of no extinction).

The $3\sigma$ detection limit of Paβ is $7.89 \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$, comparing this to the observed Paα flux allows us to put a con-

straint 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\sigma$ 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\sigma$ 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.

## 4.3. Multiwavelength 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. Toward the bright central point sources we detect very little Paα emission. Northwest of the brightest clump 4 the $3\sigma$ contour shows an extension toward 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 rest-frame 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 inward (Hodge et al. 2015).

As discussed in Sect. 5.1, this multiwavelength 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).

## 4.4. 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.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Multiwavelength comparison of GN20 with the Paα emission line map (cyan contours). (a): F560W (rest-frame 1.1 µm) image from Colina et al. (2023), (b): HST/ACS F814W (rest-frame 0.15 µm) image, (c): 880 µm (rest-frame 160 µm) image from Hodge et al. (2015), (d): CO(2–1) emission from Hodge et al. (2012). The red circles show the position of the 5 CO clumps identified by Hodge et al. (2012).</figcaption></figure>

We calculate the SFR<sub>Paα</sub> 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<sub>☉</sub> yr<sup>-1</sup> (Table 2). The total SFR<sub>Paα</sub> by summing up the values of the clumps is 225 ± 42 M<sub>☉</sub> yr<sup>-1</sup>, 55% higher than the integrated SFR<sub>Paα</sub> 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–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 ionized gas by quadratically subtracting the width of the instrumental line spread function (Labiano et al. 2021; Jones et al. 2023). The FWHM of the instrumental broadening at 9.47 µm is 88 km s<sup>-1</sup> (σ<sub>inst</sub> = 38 km s<sup>-1</sup>). From the derived velocity dispersion (Table 2), we can derive the virial mass using the isotropic virial estimator (e.g., Hodge et al. 2012):

$$M_{\text{virial}} = \frac{C\sigma^2 R_g}{G} (M_{\odot}). \qquad (1)$$

The observed velocity dispersion is represented by σ in km s<sup>-1</sup>, the gravitational radius by R<sub>g</sub> in pc and the gravitational G = 1/232 pc (km s<sup>-1</sup>)<sup>2</sup> M<sub>☉</sub><sup>-1</sup> is the gravitational constant. The factor C depends on the mass distribution and varies from C = 7.4 for
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Paα spectra and channel maps of the clumps in GN20. Left column: integrated spectra of the clumps identified in the Paα line map. The gray histograms represent the $1\sigma$ errors. The vertical dotted lines represent the borders in which the flux and integrated star formation rate is calculated (Table 2). The green histograms show the best Gaussian fits to the emission profile, with the dotted line the individual Gaussians for clumps 1–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\sigma$) peak at +125 km s$^{-1}$ in the spectrum of clump 1. 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.</figcaption></figure>

an exponential profile to $C = 4$ for a Vaucouleurs 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.37'', 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 and -60 km s$^{-1}$). In the figures we highlight the clump for which the channel maps
<table><caption>**Table 2.** Clump properties.</caption><thead><tr><th>Clump</th><th>RA</th><th>Dec</th><th>Flux<br>(10<sup>-18</sup> erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>SFR<sub>Av=0</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Velocity<br>(km s<sup>-1</sup>)</th><th>σ<br>(km s<sup>-1</sup>)</th><th>M<sub>virial</sub><br>(10<sup>9</sup> M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>1</td><td>12h37m11.91s</td><td>62°22′12.8"</td><td>6.6 ± 2.5</td><td>50 ± 19</td><td>—</td><td>—</td><td>—</td></tr><tr><td>1a</td><td>12h37m11.89s</td><td>62°22′13.0"</td><td>4.6 ± 3.1</td><td>35 ± 24</td><td>-348 ± 16</td><td>77 ± 16</td><td>&lt;6.8 ± 2.8</td></tr><tr><td>1b</td><td>12h37m11.87s</td><td>62°22′12.8"</td><td>1.9 ± 2.4</td><td>14 ± 18</td><td>-134 ± 16</td><td>37 ± 15</td><td>&lt;0.05 ± 0.04</td></tr><tr><td>2</td><td>12h37m11.98s</td><td>62°22′12.4"</td><td>6.3 ± 2.4</td><td>48 ± 18</td><td>—</td><td>—</td><td>—</td></tr><tr><td>2a</td><td>12h37m11.90s</td><td>62°22′12.8"</td><td>2.3 ± 2.5</td><td>18 ± 19</td><td>-371 ± 20</td><td>47 ± 15</td><td>&lt;1.2 ± 0.8</td></tr><tr><td>2b</td><td>12h37m11.87s</td><td>62°22′12.9"</td><td>4.4 ± 3.9</td><td>33 ± 29</td><td>-170 ± 21</td><td>77 ± 24</td><td>&lt;7.8 ± 4.8</td></tr><tr><td>3</td><td>12h37m11.94s</td><td>62°22′11.7"</td><td>7.5 ± 3.3</td><td>57 ± 25</td><td>—</td><td>—</td><td>—</td></tr><tr><td>3a</td><td>—</td><td>—</td><td>1.5 ± 2.0</td><td>12 ± 15</td><td>-74 ± 18</td><td>41 ± 15</td><td>&lt;0.4 ± 0.3</td></tr><tr><td>3b</td><td>12h37m11.91s</td><td>62°22′11.9"</td><td>6.4 ± 3.7</td><td>49 ± 28</td><td>175 ± 20</td><td>119 ± 23</td><td>&lt;1.9 ± 8</td></tr><tr><td>4</td><td>12h37m11.86s</td><td>62°22′11.9"</td><td>9.3 ± 3.1</td><td>70 ± 23</td><td>94 ± 18</td><td>138 ± 18</td><td>&lt;2.6 ± 7</td></tr></tbody></table>

are made by a blue circle. The same is shown for clump 2 and 3. 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 northeast from the red-shifted emission (clump 1a). For clump 2, the red-shifted emission (clump 2a) seems to be related to the redshifted emission in clump 1, while the blue shifted emission (clump 2b) peaks at the southeast 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 that 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. 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.

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⁻¹. The largest value (138 ± 18 km s⁻¹) 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⁻¹, 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.

### 4.5. 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 signal-to-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 PV-diagram (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 modeling 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⁻¹ and -400 km s⁻¹. The velocity dispersion varies between 85 and 190 km s⁻¹, with a typical error of 40 km s⁻¹. 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 for example 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
<figure><img src="image_8.png" /><figcaption>**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.</figcaption></figure>

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⁻¹. 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} = \frac{1}{2}(v_{max} - v_{min}) = 550 \pm 40$ km s⁻¹. This value agrees with the dynamical modeling of the CO(2–1) emission ($575 \pm 100$ km s⁻¹ 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): $\sigma_m = 145 \pm 53$ km s⁻¹. 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}/\sigma_m = 3.8 \pm 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) that have $v_{rot}/\sigma_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}/\sigma_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}/\sigma_m$ of GN20 is consistent with that of low-z LIRGs. Bellocchi et al. (2013) derive an average $v_{rot}/\sigma_m = 3.4$ for LIRGs and $v_{rot}/\sigma_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}/\sigma_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 $\sigma_m$ with theoretical predictions for $\sigma_m$ in the context of a galactic disk model. Krumholz et al. (2018) shows that $\sigma_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 $\sigma_m$ vs SFR diagram (Fig. 4 in Krumholz et al. 2018) show that the high $\sigma_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 $\sigma_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 $\sigma_m$ due to the higher mass of the clumps.

# 5. 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.

## 5.1. 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 \pm 9 M_{\odot} \text{ yr}^{-1}$ (Table 1). As shown in Sect. 4.2, the observed $3\sigma$ upperlimit of Paβ, does not provide a strict constraint on the extinction. Tan et al. (2014) derived the infrared star formation rate (SFRIR) based on 1.2–3.3 mm continuum observations, assuming a Chabrier (2003) IMF, and find a value for GN20 of SFRIR = $1860 \pm 90 M_{\odot} \text{ yr}^{-1}$. This means that the observed Paα emission only reveals $7.7 \pm 0.5\%$ of the total star formation rate as detected from the mm observations.

By comparing the SFR<sub>Paα</sub> to the SFR<sub>IR</sub> 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<sub>Paα</sub> corrected for extinction is typically lower than the SFR<sub>IR</sub>. Giménez-Arteaga et al. (2022) find again a better correlation between the SFR derived from Paβ with the SFR<sub>IR</sub>.

The disagreement between SFR<sub>Paα</sub> and SFR<sub>IR</sub> 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<sub>IR</sub>, but not to SFR<sub>Paα</sub> (Piqueras-López et al. 2016), (iii) SFR<sub>IR</sub> 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. 2023).

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). There-
fore 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 emis-
sion 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<sub>Pαα</sub> traces the most recent star formation over the last 10 Myr, while SFR<sub>IR</sub> 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<sub>Pαα</sub> and SFR<sub>IR</sub> to both these effects. The difference between SFR<sub>Pαα</sub> and SFR<sub>IR</sub> 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<sub>Pαα</sub> and SFR<sub>IR</sub>. Using the Cardelli et al. (1989) extinction law with a R<sub>V</sub> = 4.05, we find a total extinction of A<sub>V</sub> = 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 Eq. (1) in Calabrò et al. (2018), our ratio of observed SFRs gives a total visual attenuation of A<sub>V,mixed</sub> = 44 ± 3 mag.

The multiwavelength 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<sub>Pαα</sub> of the clumps will increase this fraction, for example increasing this fraction to 20% would require an A<sub>V</sub> = 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 behavior where large extinction
variations, with areas with A<sub>V</sub> as high as 20–30 mag, are found
(see also Bohn et al. 2023). These (U)LIRGs typically have an
obscured center. However, in most local (U)LIRGs the extinc-
tion corrected SFR traced by the Paschen lines follows pretty
well the SFR<sub>IR</sub> (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 mid-infrared diag-
nostics, Haas et al. (2001) derive values of A<sub>V,mixed</sub> = 110 mag for

the central starburst in Arp 220 and A<sub>V,mixed</sub> = 50 mag in the cen-
ter of UGC 5105 (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 toward 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¹⁰ 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 multiwavelength 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. 2024; 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 emis-
sion longward of 170 µm is optically thick, consistent with a very
compact starburst in the nucleus of GN20. Similar properties are
observed for other high-redshift starbursts (Jin et al. 2022).

5.2. 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 empiri-
cal PSF derived by Colina et al. (2023) that is based on close-by
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

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 10 000 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). 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.

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 northeast 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 \ge 1$ (e.g., Elmegreen et al. 2009a; Förster-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 \pm 2.4 \times 10^{11} M_{\odot}$. The upper limits we derive on the clump masses ($\sim 10^9 M_{\odot}$, 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. (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.

### 5.3. 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 SED of GN20, Tan et al. (2014) find an SFR<sub>IR</sub>/SFR<sub>UV</sub> of ~13, this would lead to an SFR<sub>UV</sub> of $140 M_{\odot} \text{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<sup>2</sup> in size, matching the size and morphology of the UV-bright emission ($1.6'' \times 0.8''$, PA: 0°).

Under the assumption of no extinction we calculate the SFR<sub>UV</sub> from the measured F814W flux using the relation detailed by Murphy et al. (2011) and find SFR<sub>UV</sub> = $15 \pm 2 M_{\odot} \text{yr}^{-1}$. As a next step we perform a fit using CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019) to the integrated UV-optical 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 and 500 Myr. For the young stellar population (relevant for the UV emission) we find a SFR = $135 \pm 70 M_{\odot} \text{yr}^{-1}$ and an extinction of $A_V = 0.8$ mag.

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_\odot$ yr$^{-1}$ for the young population, with a similar extinction as the previous fit.

Based on the standard deviation derived for the Paα image (Sect. 4.2) we derive a $3\sigma$ detection limit of $2.5 M_\odot$ yr$^{-1}$ per pixel. Assuming a constant spatial distribution and an MRS pixel size of 0.029 arcsec$^2$ gives a $3\sigma$ upper limit of $15 M_\odot$ 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$ mag, 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_\odot$ 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 toward 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.

## 5.4. The origin of the starburst in GN20

Two main scenarios have 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 & Hernquist (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.

## 6. 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 \pm 9 M_\odot$ yr$^{-1}$, assuming no extinction. The SFR$_{Paα}$ is only a small fraction of the SFR$_{IR}$ ($7.7 \pm 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 \pm 0.4$ mag. Such high values are also measured in the nearby ULIRGs Arp 220 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 sub-clumps 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 toward 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.

**Acknowledgements.** We thank the anonymous referee for helpful and constructive comments which improved this paper. We thank Jacqueline Hodge for providing the CO(2–1) data cube and moment maps. A.B. thanks Angela Adamo and Sean Linden for discussions on clumps and ULIRGs. A.B., J.M., and G.Ö. acknowledge support from the Swedish National Space Administration (SNSA). L.C., J.A.-M., and A.C.G. acknowledge support by grant PIB2021-127718NB-100 from the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe". F.P. and A.E. acknowledge support through the German Space Agency DLR 500S1501 and DLR 500S2001 from 2015 to 2023. K.I.C and E.I. acknowledge funding from the Netherlands Research
School for Astronomy (NOVA). K.I.C. acknowledges funding from the Dutch Research Council (NWO) through the award of the Vici Grant VI.C.212.036. A.A.H. acknowledges support from grant PID2021-124665NB-100 funded by MCIN/AEI/10.13039/501100011033 and by ERDF A way of making Europe. S.G. acknowledges financial support from the Villum Young Investigator grant 37440 and 13160 and the Cosmic Dawn Center (DAWN). The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. T.R.G. is grateful for support from the Carlsberg Foundation via grant No. CF20-0534. J.P.P. and T.T. acknowledge financial support from the UK Science and Technology Facilities Council, and the UK Space Agency. J.H and D.L. were supported by a VILLUM FONDEN Investigator grant (project number 16599) P.G.P.-G. acknowledges support from grant PID2022-139567NB-100 funded by Spanish Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033, FEDER, UE. P.-O.L. acknowledges financial support from CNES. T.R. acknowledges support from ERC grant no. 743029 (EASY). The work presented is the effort of the entire MIRI team and the enthusiasm within the MIRI partnership is a significant factor in its success. MIRI draws on the scientific and technical expertise of the following organisations: Ames Research Center, USA; Airbus Defence and Space, UK; CEA-Irfu, Saclay, France; Centre Spatial de Liège, Belgium; Consejo Superior de Investigaciones Científicas, Spain; Carl Zeiss Optronics, Germany; Chalmers University of Technology, Sweden; Cosmic Dawn Center (DAWN), DTU Space, Technical University of Denmark, Denmark; Dublin Institute for Advanced Studies, Ireland; European Space Agency, The Netherlands; ETCA, Belgium; ETH Zurich, Switzerland; Goddard Space Flight Center, USA; Institute d'Astrophysique Spatiale, France; Instituto Nacional de Técnica Aeroespacial, Spain; Institute for Astronomy, Edinburgh, UK; Jet Propulsion Laboratory, USA; Laboratoire d'Astrophysique de Marseille (LAM), France; Leiden University, The Netherlands; Lockheed Advanced Technology Center (USA); NOVA Opt-IR group at Dwingeloo, The Netherlands; Northrop Grumman, USA; Max-Planck Institut für Astronomie (MPIA), Heidelberg, Germany; Laboratoire d'Études Spatiales et d'Instrumentation en Astrophysique (LESIA), France; Paul Scherrer Institut, Switzerland; Raytheon Vision Systems, USA; RUAG Aerospace, Switzerland; Rutherford Appleton Laboratory (RAL Space), UK; Space Telescope Science Institute, USA; Toegepast- Natuurwetenschappelijk Onderzoek (TNO-TPD), The Netherlands; UK Astronomy Technology Centre, UK; University College London, UK; University of Amsterdam, The Netherlands; University of Arizona, USA; University of Cardiff, UK; University of Cologne, Germany; University of Ghent; University of Groningen, The Netherlands; University of Leicester, UK; University of Leuven, Belgium; University of Stockholm, Sweden; Utah State University, USA. A portion of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We would like to thank the following National and International Funding Agencies for their support of the MIRI development: NASA; ESA; Belgian Science Policy Office; Centre Nationale d'Études Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luft und Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economía y Comunicidad; Netherlands Research School for Astronomy (NOVA); Netherlands Organisation for Scientific Research (NWO); Science and Technology Facilities Council; Swiss Space Office; Swedish National Space Board; UK Space Agency. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. Some data were obtained from the *Mikulski* Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST, and from the European JWST archive (eJWST, https://jwst.esac.esa.int/archive/) operated by the ESDC.

## References

Aalto, S., Martín, S., Costagliola, F., et al. 2015, A&A, 584, A42

Adamo, A., Östlin, G., Bastian, N., et al. 2013, ApJ, 766, 105

Adamo, A., Hollyhead, K., Messa, M., et al. 2020, MNRAS, 499, 3267

Akins, H. B., Casey, C. M., Allen, N., et al. 2023, ApJ, 956, 61

Algera, H. S. B., Inami, H., Oesch, P. A., et al. 2022, MNRAS, 518, 6142

Alonso-Herrero, A., Rieke, G. H., Rieke, M. J., et al. 2006, ApJ, 650, 835

Álvarez Márquez, J., Gómez, A. C., Colina, L., et al. 2023, A&A, 671, A105

Álvarez Márquez, J., Colina, L., Gómez, A. C., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202347946

Argyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111

Arribas, S., Colina, L., Bellocchi, E., Maiolino, R., & Villar-Martín, M. 2014, A&A, 568, A14

Arribas, S., Perna, M., Pino, B. R. D., et al. 2023, A&A, submitted [arXiv:2312.00899]

Barrufet, L., Oesch, P. A., Bouwens, R., et al. 2023a, MNRAS, 522, 3926

Barrufet, L., Oesch, P. A., Weibel, A., et al. 2023b, MNRAS, 522, 449

Bellocchi, E., Arribas, S., Colina, L., & Miralles-Caballero, D. 2013, A&A, 557, A59

Bellocchi, E., Pereira-Santaella, M., Colina, L., et al. 2022, A&A, 664, A60

Bik, A., Östlin, G., Hayes, M., Melinder, J., & Menacho, V. 2022, A&A, 666, A161

Birrer, S., & Amara, A. 2018, Phys. Dark Univ., 22, 189

Blecha, L., Snyder, G. F., Satyapal, S., & Ellison, S. L. 2018, MNRAS, 478, 3056

Bohn, T., Inami, H., Diaz-Santos, T., et al. 2023, ApJ, 942, L36

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Borys, C., Chapman, S., Halpern, M., & Scott, D. 2003, MNRAS, 344, 385

Bosman, S. E. I., Álvarez Márquez, J., Colina, L., et al. 2023, ArXiv e-prints [arXiv:2307.14414]

Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, AJ, 162, 47

Bouwens, R. J., Illingworth, G., Ellis, R. S., Oesch, P., & Stefanon, M. 2022, ApJ, 940, 55

Bradley, L. 2023, https://doi.org/10.5281/zenodo.7946442

Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, https://doi.org/10.5281/zenodo.7577320

Calabrö, A., Daddi, E., Cassata, P., et al. 2018, ApJ, 862, L22

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Caputi, K. I., Deshmukh, S., Ashby, M. L. N., et al. 2017, ApJ, 849, 45

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Carilli, C. L., Hodge, J., Walter, F., et al. 2011, ApJ, 739, L33

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

Cava, A., Schaerer, D., Richard, J., et al. 2018, Nat. Astron., 2, 76

Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151

Chabrier, G. 2003, PASP, 115, 763

Chomiuk, L., & Povich, M. S. 2011, AJ, 142, 197

Cochrane, R. K., Anglés-Alcázar, D., Cullen, F., & Hayward, C. C. 2024, ApJ, 961, 37

Colina, L., Gómez, A. C., Álvarez Márquez, J., et al. 2023, A&A, 673, L6

Cortzen, I., Magdis, G. E., Valentino, F., et al. 2020, A&A, 634, L14

Costantin, L., Pérez-González, P. G., Vega-Ferrero, J., et al. 2023, ApJ, 946, 71

Daddi, E., Dannerbauer, H., Stern, D., et al. 2009, ApJ, 694, 1517

Dekel, A., Birnboim, Y., Engel, G., et al. 2009a, Nature, 457, 451

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

Dessauges-Zavadsky, M., & Adamo, A. 2018, MNRAS, 479, L118

Donnan, F. R., Rigopoulou, D., García-Bernete, I., et al. 2023, ApJ, 669, A87

Drake, A. B., Walter, F., Novak, M., et al. 2020, ApJ, 902, 37

Elmegreen, B. G., Bournaud, F., & Elmegreen, D. M. 2008, ApJ, 688, 67

Elmegreen, D. M., Elmegreen, B. G., Marcus, M. T., et al. 2009a, ApJ, 701, 306

Elmegreen, B. G., Elmegreen, D. M., Fernandez, M. X., & Lemonias, J. J. 2009b, ApJ, 692, 12

Engel, H., Davies, R. I., Genzel, R., et al. 2011, ApJ, 729, 58

Finkelstein, S. L., Leung, G. C. K., Bagley, M. B., et al. 2023, ApJ, submitted [arXiv:2311.04279]

Fujimoto, S., Ouchi, M., Shibuya, T., & Nagai, H. 2017, ApJ, 850, 83

Förster-Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661

Förster-Schreiber, N. M., Genzel, R., Bouché, N., et al. 2020, ApJ, 706, 1364

Förster-Schreiber, N. M., Shapley, A. E., Genzel, R., et al. 2011, ApJ, 739, 45

Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

Gasman, D., Argyriou, I., Sloan, G. C., et al. 2023, A&A, 673, A102

Genzel, R., Burkert, A., Bouché, N., et al. 2008, ApJ, 687, 59

Gillman, S., Gullberg, B., Brammer, G., et al. 2023, A&A, 676, A26

Giménez-Arteaga, C., Brammer, G. B., Marchesini, D., et al. 2022, ApJS, 263, 17

Ginolfi, M., Jones, G. C., Béthermin, M., et al. 2020, A&A, 643, A7

Girard, M., Dessauges-Zavadsky, M., Schaerer, D., et al. 2018, A&A, 613, A72

Glazebrook, K. 2013, PASA, 30, 7

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Gómez-Guijarro, C., Toft, S., Karim, A., et al. 2018, ApJ, 856, 121

Gómez-Guijarro, C., Magnelli, B., Elbaz, D., et al. 2023, A&A, 677, A34

Haas, M., Klaas, U., Müller, S. A. H., Chini, R., & Coulson, I. 2001, A&A, 367, L9

Hodge, J. A., Carilli, C. L., Walter, F., et al. 2012, ApJ, 760, 11

Hodge, J. A., Riechers, D., Decarli, R., et al. 2015, ApJ, 798, L18

Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103

Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Nature, 394, 241

Hung, C.-L., Rich, J. A., Yuan, T., et al. 2015, ApJ, 803, 62

Inoue, S., Dekel, A., Mandelker, N., et al. 2015, MNRAS, 456, 2052

Jin, S., Daddi, E., Magdis, G. E., et al. 2022, A&A, 665, A3

Jones, O. C., Álvarez-Márquez, J., Sloan, G. C., et al. 2023, MNRAS, 523, 2519

Jones, G. C., Übler, H., Perna, M., et al. 2024, A&A, 682, A122

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kroupa, P. 2001, MNRAS, 322, 231

Krumholz, M. R., Burkhart, B., Forbes, J. C., & Crocker, R. M. 2018, MNRAS, 477, 2716
Labiano, A., Argyriou, I., Álvarez Márquez, J., et al. 2021, A&A, 656, A57

Larson, K. L., Díaz-Santos, T., Armus, L., et al. 2020, ApJ, 888, 92

Law, D. R., Steidel, C. C., Erb, D. K., et al. 2009, ApJ, 697, 2057

Law, D. R., Morrison, J. E., Argyriou, I., et al. 2023, AJ, 166, 45

Livermore, R. C., Jones, T., Richard, J., et al. 2012, MNRAS, 427, 688

Lucy, L. B. 1974, AJ, 79, 745

Mandelker, N., Dekel, A., Ceverino, D., et al. 2014, MNRAS, 443, 3675

Messa, M., Dessauges-Zavadsky, M., Richard, J., et al. 2022, MNRAS, 516, 2420

Meštrić, U., Vanzella, E., Zanella, A., et al. 2022, MNRAS, 516, 3532

Morrison, J. E., Dicken, D., Argyriou, I., et al. 2023, PASP, 135, 075004

Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67

Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793

Ormerod, K., Conselice, C. J., Adams, N. J., et al. 2023, MNRAS, 527, 6110

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Sausalito: University Science Books)

Parlanti, E., Carniani, S., Übler, H., et al. 2024, A&A, 684, A24

Patapis, P., Argyriou, I., Law, D. R., et al. 2024, A&A, 682, A53

Peißker, F., Eckart, A., Zajaček, M., & Britzen, S. 2022, ApJ, 933, 49

Pereira-Santaella, M., Colina, L., García-Burillo, S., et al. 2016, A&A, 587, A44

Piqueras-López, J., Colina, L., Arribas, S., & Alonso-Herrero, A. 2013, A&A, 553, A85

Piqueras-López, J., Colina, L., Arribas, S., Pereira-Santaella, M., & Alonso-Herrero, A. 2016, A&A, 590, A67

Planck Collaboration XIII. 2016, A&A, 594, A13

Pope, A., Scott, D., Dickinson, M., et al. 2006, MNRAS, 370, 1185

Popping, G., Pillepich, A., Calistro Rivera, G., et al. 2022, MNRAS, 510, 3321

Pérez-González, P. G., Rieke, G. H., Egami, E., et al. 2005, ApJ, 630, 82

Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023a, ApJ, 951, L1

Pérez-González, P. G., Barro, G., Annunziatella, M., et al. 2023b, ApJ, 946, L16

Ricci, C., Bauer, F. E., Treister, E., et al. 2017, MNRAS, 468, 1273

Riechers, D. A., Pope, A., Daddi, E., et al. 2014, ApJ, 786, 31

Riechers, D. A., Hodge, J. A., Pavesi, R., et al. 2020, ApJ, 895, 81

Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584

Rizzo, F., Vegetti, S., Fraternali, F., Stacey, H. R., & Powell, D. 2021, MNRAS, 507, 3952

Romeo, A. B., & Agertz, O. 2014, MNRAS, 442, 1230

Romeo, A. B., Burkert, A., & Agertz, O. 2010, MNRAS, 407, 1223

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Spilker, J. S., Hayward, C. C., Marrone, D. P., et al. 2022, ApJ, 929, L3

Spilker, J. S., Phadke, K. A., Aravena, M., et al. 2023, Nature, 618, 708

Springel, V., & Hernquist, L. 2005, ApJ, 622, L9

Sturm, E., Lutz, D., Genzel, R., et al. 1996, A&A, 315, L133

Sánchez-García, M., Pereira-Santaella, M., García-Burillo, S., et al. 2022, A&A, 659, A102

Tan, Q., Daddi, E., Magdis, G., et al. 2014, A&A, 569, A98

Tateuchi, K., Konishi, M., Motohara, K., et al. 2015, ApJS, 217, 1

Ueda, J., Iono, D., Yun, M. S., et al. 2014, ApJS, 214, 1

Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93

Vieira, J. D., Marrone, D. P., Chapman, S. C., et al. 2013, Nature, 495, 344

Ward, E., de la Vega, A., Mobasher, B., et al. 2024, ApJ, 962, 176

Wells, M., Pel, J.-W., Glasse, A., et al. 2015, PASP, 127, 646

Wright, G. S., Wright, D., Goodson, G. B., et al. 2015, PASP, 127, 595

Wright, G. S., Rieke, G. H., Glasse, A., et al. 2023, PASP, 135, 048003

Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, ApJ, 909, 165

Zavala, J. A., Buat, V., Casey, C. M., et al. 2023, ApJ, 943, L9

<footnote n="1">¹ Department of Astronomy, Stockholm University, Oscar Klein Centre, AlbaNova University Centre, 106 91 Stockholm, Sweden
e-mail: arjan.bik@astro.su.se</footnote>

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

<footnote n="3">³ I. Physikalisches Institut der Universität zu Köln, Zülpicher Str. 77, 50937 Köln, Germany</footnote>

<footnote n="4">⁴ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="5">⁵ DTU Space, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="6">⁶ Cosmic Dawn Centre (DAWN), Copenhagen, Denmark</footnote>

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

<footnote n="8">⁸ UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="9">⁹ Centro de Astrobiología (CAB), CSIC-INTA, Camino Viejo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain</footnote>

<footnote n="10">¹⁰ Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands</footnote>

<footnote n="11">¹¹ DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark</footnote>

<footnote n="12">¹² Telespazio UK for the European Space Agency, ESAC, Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Spain</footnote>

<footnote n="13">¹³ School of Physics & Astronomy, Space Research Centre, Space Park Leicester, University of Leicester, 92 Corporation Road, Leicester LE4 5SP, UK</footnote>

<footnote n="14">¹⁴ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="15">¹⁵ University of Vienna, Department of Astrophysics, Türkenschanzstrasse 17, 1180 Vienna, Austria</footnote>

<footnote n="16">¹⁶ Institute of Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-Str 27, 8093 Zurich, Switzerland</footnote>

<footnote n="17">¹⁷ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="18">¹⁸ Dublin Institute for Advanced Studies, Astronomy & Astrophysics Section, 31 Fitzwilliam Place, Dublin D02 XF86, Ireland</footnote>