Detailed study of a rare hyperluminous rotating disk in an Einstein ring 10 billion years ago

Abstract

Hyperluminous infrared galaxies (HyLIRGs) are the rarest and most extreme starbursts and found only in the distant Universe (z ≳ 1). They have intrinsic infrared (IR) luminosities L_IR ≥ 10¹³ L_⊙ and are commonly found to be major mergers. Recently, the Planck All-Sky Survey to Analyze Gravitationally-lensed Extreme Starbursts project (PASSAGES) searched ~10⁴ deg² of the sky and found ~20 HyLIRGs. We describe a detailed study of PJ0116-24, the brightest (μL_IR ≈ 2.6 × 10¹⁴ L_⊙, magnified with μ ≈ 17) Einstein-ring HyLIRG in the southern sky, at z = 2.125, with observations from the near-IR integral-field spectrograph VLT/ERIS and the submillimetre interferometer ALMA. We detected Hα, Hβ, [N ii] and [S ii] lines and obtained an extreme Balmer decrement (Hα/Hβ ≈ 8.73 ± 1.14). We modelled the molecular-gas and ionized-gas kinematics with CO(3–2) and Hα data at ~100–300 pc and (sub)kiloparsec delensed scales, respectively, finding consistent regular rotation. We found PJ0116-24 to be highly rotationally supported (v_rot/σ_{0, mol. gas} ≈ 9.4) with a richer gaseous substructure than other known HyLIRGs. Our results imply that PJ0116-24 is an intrinsically massive (M_baryon ≈ 10^11.3 M_⊙) and rare starbursty disk (star-formation rate, SFR = 1,490 M_⊙ yr⁻¹) probably undergoing secular evolution. This indicates that the maximal SFR (≳1,000 M_⊙ yr⁻¹) predicted by simulations could occur during a galaxy’s secular evolution, away from major mergers.

Main

All-sky surveys conducted with infrared (IR) satellites, such as the Infrared Astronomical Satellite (IRAS), the Wide-field Infrared Survey Explorer (WISE) and the Planck satellite, have led to great breakthroughs in finding rare hyperluminous infrared galaxies (HyLIRGs)^1,2,3,4,5,6. The brightest HyLIRGs have IR luminosities μL_IR > 10¹⁴ L_⊙ and are often found to be gravitationally lensed. High-resolution imaging with the Hubble Space Telescope (HST) in the optical and near-IR and with (sub)millimetre and radio interferometers like the Very Large Array (VLA), the Northern Extended Millimeter Array (NOEMA) and the Atacama Large Millimeter/Submillimeter Array (ALMA) have revealed that HyLIRGs are almost exclusively gas-rich mergers^{7,8,9,10,11,12}. In most cases, they also host an active galactic nucleus (AGN).

A widely accepted scenario is that HyLIRGs are the distant higher-luminosity tail of the local ultra-luminous IR galaxies (ULIRGs; L_IR > 10¹² L_⊙) with extreme starburst activities triggered by major mergers. Thus, they fit into the well-known merger-driven ULIRG and quasar paradigm¹³. Theoretically, however, a second possibility has been raised that HyLIRGs are very young, ‘primaeval’ galaxies that are going through their maximal star formation periods¹⁴. Hydrodynamical simulations^15,16 of massive, turbulent, gas-rich galaxies have provided evidence that extreme star-formation rates (SFR ≈ 1,000 M_⊙ yr⁻¹) can be achieved in the brightest submillimetre galaxies at redshift 2–3. Star formation is fuelled by cold gas rapidly inspiralling towards the centre during their secular evolution, that is, away from any major merger. The theories are still not robustly constrained by observations due to the rarity of HyLIRGs with high-resolution, high-quality, multi-tracer observations.

A unique HyLIRG under the microscope of strong lensing

In recent years, hundreds of apparent HyLIRGs have been discovered from ≥10⁴ deg² of Planck, Herschel and South Pole Telescope (SPT) surveys, including PASSAGES, which searched ~11,500 deg² of sky^4,5,6. Among these, there are only a few tens of HyLIRGs with μL_IR > 10¹⁴ L_⊙, and those that are lensed into Einstein rings (thus, with large magnifications) are even rarer. Here we select the brightest Einstein-ring HyLIRG in the southern sky from PASSAGES, PJ011646.8-243702 (hereafter PJ0116-24), at z = 2.125 and with μL_IR ≈ 2.5 × 10¹⁴ L_⊙ (see physical properties in Table 1), for a pioneering and comprehensive study of its heavy dust attenuation and kinematics, about which there is presently little information for the HyLIRG population. We obtained high-resolution CO(3–2) observations with ALMA to trace the cold-gas distribution and kinematics at ~0.16″ angular resolution and high-quality near-IR integral field unit (IFU) observations with a new versatile instrument, the enhanced resolution imager and spectrograph (ERIS)¹⁷, installed on the Very Large Telescope (VLT)’s UT4, to trace the ionized-gas distribution and kinematics at ~0.68″ angular resolution. We also measured the global Balmer decrement and metallicity using the observed rest-frame optical spectral lines (Table 2).

Table 1 Physical properties of PJ0116-24

Table 2 Line fluxes, Balmer decrement, dust attenuation and SFRs derived from ERIS IFU data

As shown in Fig. 1a, PJ0116-24 is well resolved into a nearly complete Einstein ring in the HST 1.6 µm image (ref. ⁶ and J. Lowenthal et al., manuscript in preparation), which provides rich structural information: two bright emission peaks A1 and A2; several stretched arcs; and four clump images B1 to B4. The distribution of cold gas, as traced by the CO(3–2) line-integrated emission (Methods), exhibits extended structures that surround the HST emission peaks.

Fig. 1: Stellar, cold-gas and ionized-gas emission distributions of the rare HyLIRG Einstein ring, PJ0116-24.

a, Lensed distribution of the stellar light as traced by HST 1.6 µm emission to a Sauron colour scale with blue-white-green-yellow. The foreground lens has been subtracted, and the PSF ≈ 0.2″ is indicated at bottom left. Also shown in red is the lensed distribution of the cold gas as traced by the ALMA CO(3–2) line intensity. The beam, 0.19″ × 0.16″, is indicated at bottom right. b, CO(3–2) line-centre velocity map (−330 to +330 km s⁻¹ from dark blue to light pink). c, Comparison of the distributions of the cold gas (red) and the ionized gas as traced by the ERIS/VLT Hα line intensity (white). d, Source-plane distributions (delensed from the full Einstein ring) of emission from stars, cold gas and ionized gas. The AGN at the galaxy centre is lensed into the positions A1 and A2 as labelled in a and b. The clump at the north-east ~4 kpc from the centre is lensed to the positions B1 to B4 as labelled in a and c. (The further north Hα blob is an artefact of delensing due to the noise peaks, for example, around B2, scattered across the caustic curves.) Caustic curves from our lens model are indicated by the yellow lines in d, and white lines mark the CO arm-like substructures. The CO(3–2) ‘delensed beams’ at the galaxy centre reconstructed from positions A1 and A2 are indicated by the red ellipses at bottom left (FWHM ≈ 0.09″ × 0.06″ and ~0.15″ × 0.04″, respectively; the delensed beam is even smaller closer to the caustic lines).

The CO(3–2) line-centre velocity map, as shown in Fig. 1b, clearly exhibits an ordered rotation with line-of-sight velocities ranging from about −330 to +330 km s⁻¹ relative to the systematic redshift velocity of PJ0116-24. We performed detailed kinematic modelling and quantified the rotation-to-dispersion ratio as v_rot/σ_{0, mol. gas} = 9.4 ± 0.9 for the molecular gas and v_rot/σ_{0, ion. gas} = 5.7 ± 0.6 for the ionized gas at the baryonic-mass effective radius of R_e ≈ 3.6 kpc (corrected for beam-smearing, line spread function and inclination; Methods). This rapid and ordered rotation is rarely found in HyLIRGs and is best constrained among all known HyLIRGs. The HST emission peaks, A1 and A2, are at the near-zero velocity positions and map to the same source-plane (intrinsic) position, confirming that they are both the photometric and kinematic centres of the galaxy.

The ionized-gas distribution, as traced by the Hα line-integrated emission from the ERIS IFU data (Methods) in Fig. 1c, further shows an interesting configuration. The gas extends across the disk, tends to be more pronounced on the blueshifted side of the galaxy and has consistent detections at all B1–B4 locations. Given that Hα emission is sensitive to dust attenuation, the blue/red-side asymmetry could be the combined effect of spatially varying dust attenuation, the underlying distribution of star formation across the galaxy and lensing.

Through detailed new lens modelling using high-resolution HST and CO channel maps and mesh-grid-based delensing (Methods), we reconstructed the source-plane CO(3–2) and Hα line-integrated emission distributions as shown in Fig. 1d. The magnification of the total flux of the resolved HST and CO(3–2) emission is estimated as ~16 and ~17, respectively, using the observed and delensed maps. Due to the lensing effect, the spatial resolution (beam or point spread function (PSF)), which is uniform in the observed maps, varies with location in the source-plane delensed maps. For instance, when delensing a Gaussian beam at the A2 image of the galaxy centre, the delensed beam has an approximately Gaussian shape with a minor-axis full width at half-maximum (FWHM) of ~300 pc and an axis ratio of ~0.3, elongated along north-east to south-west. This delensed beam varies across the source plane. A maximum resolution of ≤100 pc was achieved along the Einstein ring when crossing the critical curve (for example, between A1–B4 and A2–B2; Methods).

Under the microscope of strong lensing, we found that PJ0116-24 has the intrinsic characteristics of massive star-forming disk galaxies at z ≈ 2 that have been well studied with adaptive optics at kiloparsec scales¹⁸. That is, it has a central stellar concentration, an extended, massive, cold-gas-rich disk, and giant (kiloparsec scale) star-forming clump(s) that formed due to the strong Toomre instability of a gas-rich disk¹⁸. The CO-traced gas structures are within about one effective radius of ~3.6 kpc and have a Toomre’s Q parameter Q_gas ≈ 0.3 (based on equation (3) of ref. ¹⁹, Q_gas = a / f_gas × (σ₀ / v_rot), where a ≈ 1.4 is a factor for a general flat rotation curve, f_gas = 0.56 is the gas fraction and v_rot/σ_{0, mol. gas} = 9.4 is from our kinematic modelling). This Q_gas is well below the typical stability criterion for a thick disk Q_crit ≈ 0.67 (ref. ¹⁹). Thus, the gas naturally fragments into structures at (sub) Toomre scales or (sub)kiloparsec sizes (λ_T ≈ f_gasR_disk; ref. ²⁰). In our work, we resolved very rich substructures down to ~100–300 pc within the most rotationally supported HyLIRG.

We performed rest-frame ultraviolet (UV) to millimetre spectral energy modelling (Methods) and found that PJ0116-24 is forming stars at a rate of SFR_{UV + IR} ≈ 1,490 ± 400 M_⊙ yr⁻¹, rivalling the small number of known highest-SFR systems (~1,000–3,000 M_⊙ yr⁻¹) at high redshift^21,22,23,24. The non-major-merger scenario of PJ0116-24 is further supported by the cold-gas distribution, which exhibits a deficit within the central ~500 pc but is enriched at 1–2 kpc and extends out to ~3–4 kpc with stream- or arm-like features (as outlined by the white lines in Fig. 1d). These features can naturally form in turbulent, gas-rich, secularly evolving disks when gas is accreted and transported into the centre²⁵.

A giant clump at ~4 kpc in projection north-east of the galaxy centre in the source plane is lensed into the B1–B4 emission peaks. It is detected in all stellar, Hα and CO tracers and is probably a kiloparsec-scale, young, star-forming clump with either stellar or gas mass accounting for only ≤3% of that of PJ0116-24 (Methods). With FWHM ≈ 1.2 kpc for a two-dimensional (2D) Gaussian fit of the delensed HST image, it is ×10 to ×100 larger than the star-forming complexes in the Milky Way but typical of giant clumps in massive star-forming galaxies at z ≈ 2–3. These clumps also naturally form due to the galaxies’ high gas fractions (f_gas ≡ M_gas/(M_gas + M_*) ≳ 0.5)²⁰ and increased Toomre instability (refs. ^{19,20,26,27,28,29}; see also simulations in refs. ^25,30,31 and Methods).

Figure 2 provides a clearer view of the delensed distributions of cold and ionized gas. The two gas tracers exhibit remarkably consistent velocity, velocity dispersion and position–velocity (P–V) maps, albeit apparently affected by corresponding beam-smearing, which thus needs detailed dynamical modelling. We performed state-of-the-art dynamical modelling based on direct fitting of the high-resolution CO data in the image plane (Methods), which reproduced well the P–V diagrams of both CO and Hα under different beam-smearing conditions (Fig. 2g,h). The centrally peaked apparent velocity dispersion is also strong evidence for an ordered disk rather than a major merger¹⁸. Additionally, the clump appears to follow the systematic rotation and has an indistinguishable cold-gas velocity dispersion as the disk, again in line with a non-merger scenario.

Fig. 2: CO and Hα line intensity and kinematic maps of PJ0116-24 after correcting for lensing.

a,b Line intensity maps for CO (a) and Hα (b). c,d, Velocity maps for CO (c) and Hα (d). e,f, Velocity dispersion maps (vel. disp. maps) for CO (e) and Hα (f). g,h, P–V diagrams for CO (g) and Hα (h). HST stellar light distribution is indicated by the contours in a–f. The angular resolution for reconstructions from images A1 and A2 is indicated by the ellipses at the lower left corners of a and b. Colour bar next to c and d displays the velocities in units of km s⁻¹. Colour bar next to e and f displays the velocity dispersion in units of km s⁻¹. In g and h, contours represent our best-fitting kinematic model’s P–V diagrams with the same beam-smearing effect as the data. Dotted lines indicate our best-fitting intrinsic rotation curves without beam-smearing but with the inclination effect. Thus, they match the contours where the rotation curve flattens and the beam-smearing is minimized. Dashed lines indicate the intrinsic rotation curve after beam-smearing and inclination corrections. In h, the [N ii] λ6584, Hα and [N ii] λ6548 lines (from top to bottom) are shown in the same map due to their proximity in wavelength and velocity. The vertical axis of h indicates the velocity with zero centred at the Hα line.

Balmer decrement, metallicity and ionization source

Our ERIS Science Verification Program (principal investigator (PI) D. Liu) observed PJ0116-24 in two near-IR bands: the H-band covering the redshifted Hβ and [O iii] λλ4959,5007 lines and the K-band covering the Hα, [N ii] λλ6548,6584 and [S ii] λλ6717,6731 lines. All lines were detected with signal-to-noise ratio (S/N) ≈ 2–20. This dataset provides a full suite of strong-line diagnostics for studying the galaxy’s dust attenuation, metallicity and ionization source.

The ratio of Hα to Hβ line fluxes, known as the Balmer decrement, is a reliable indicator of the dust attenuation of the emitting source. In case b recombination with an electron temperature T_e = 10⁴ K and density n_e ≈ 10²–10⁴ cm⁻³, appropriate for star-forming regions, the intrinsic, unattenuated Hα/Hβ flux ratio is 2.86. In comparison, we observed an attenuated Hα/Hβ flux ratio of 8.73 ± 1.14 (corrected for stellar Balmer absorption). Assuming a Calzetti et al.³² attenuation curve with R_V = 3.1 for the nebular line, we obtained a nebular reddening of A_{V, neb} = 2.95 ± 0.13 mag. This is consistent with the dust attenuation determined through our fitting of the spectral energy distribution (SED; Methods). Converting the Hα line flux to an intrinsic, attenuation-corrected SFR³³, we obtained SFR_{Hα, corr} ≈ 470 ± 60 M_⊙ yr⁻¹, ~31% ± 9% of the total SFR_{UV + IR}. The apparent, uncorrected SFR_{Hα, not corr} ≈ 58 ± 2 M_⊙ yr⁻¹ is as low as ~3.9% ± 1.1%. The small fractions indicate that a substantial amount of Hα emission is too obscured by dust to be detected. The high magnification in this Einstein ring helps to unveil the less obscured parts of the galaxy, but the measured Balmer decrement places only a lower limit on the total dust obscuration for the whole galaxy.

Figure 3 shows the ERIS H-band Hβ, [O iii] doublet, and K-band Hα and [N ii] doublet spectra, as well as a Balmer decrement versus stellar mass diagram. PJ0116-24 occupies a unique parameter space distinct from the vast majority of local galaxies in the Sloan Digital Sky Survey (SDSS) and z ≈ 2 main-sequence star-forming galaxies in the MOSDEF survey (see Methods for detailed comparisons). The lack of galaxies in this unique parameter space of high M_* and high Hα/Hβ is for several reasons: the decline of stellar mass function at the massive end, the decline of the IR luminosity function at the bright end and the challenge of performing rest-frame optical spectroscopy for heavily obscured galaxies like PJ0116-24. The PASSAGES sample will be ideal to fill this parameter space thanks to the strong lensing that makes rest-frame optical IFU imaging spectroscopy possible. Figure 4 further presents several sets of metallicity diagnostics using various line ratios as labelled in each panel. Using the Curti et al.³⁴ calibrations, which are shown as the solid blue lines in Fig. 4, we interpolated the metallicity 12 + log(O/H) from the observed line ratios. These diagnostics indicate solar to supersolar metallicity. This is much higher than in non-starburst galaxies at the same redshifts, which is typically ~1/10 solar³⁵. The N2 (Hα/[N ii]) and RS32 ([S ii]/Hβ + [O iii]/Hα) diagnostics indicate an ~0.2 dex supersolar metallicity. Whether the supersolar metallicity from the only two calibrations is reliable or not is still an open question awaiting more HyLIRGs studied in the same way.

Fig. 3: ERIS spectra of PJ0116-24.

a–c, Rest-frame optical lines integrated over the entire Einstein ring. a, Hβ. b, [O iii] λλ4959,5007. c, Hα and [N ii] λλ6548,6584. We fitted the Hα and [N ii] λλ6548,6584 lines each with a narrow (FWHM ≈ 240 km s⁻¹) and a broad component (FWHM ≈ 500 km s⁻¹). d, Balmer decrement (Hα/Hβ line flux ratio) versus stellar mass of PJ0116-24 compared with SDSS z ≈ 0 galaxies (intensity image with colour indicating the number density), MOSDEF z ≈ 1–3 galaxies, and another HyLIRG G165 Arc 1a and its companion galaxy G165 NS 46 at z ≈ 2 (see references and discussion in Methods). We show PJ0116-24’s Hα/Hβ uncorrected for the stellar Balmer absorption (st.uncorr.) as an open pentagon for comparison. Error bars indicate observational uncertainties for individual data points (mean values ± standard deviation (s.d.)).

Fig. 4: Strong-line diagnostics of PJ0116-24’s gas-phase metallicity.

Each panel indicates a strong-line diagnostic method using the Hα, Hβ, [N ii], [S ii] and [O iii] line-integrated flux measurements from the ERIS H- and K-band IFU data (Methods). Lines indicate empirical relations for star-forming galaxies from ref. ³⁴ (Methods), with shading indicating the scatter (s.d.) from ref. ³⁴. Solar metallicity is indicated by the vertical line. Error bars indicate observational uncertainties for individual data points (mean values ± s.d.).

Our analysis using the BPT diagnostic diagram (as detailed in Methods) further reveals that PJ0116-24 is dominated by star formation, with a low contribution from the AGN, in line with recent studies using X-ray spectroscopy³⁶ and high-resolution radio observations⁶. For the highest-S/N Hα and [N ii] lines, their line profiles are spectroscopically well resolved, as shown in Fig. 3, and our analysis indicates that the complex line profile composed of a narrow component (FWHM ≈ 240 km s⁻¹) and a broad component (FWHM ≈ 510 km s⁻¹) is mainly due to the asymmetric lensing of the rotating disk (with the blue side being more magnified; Methods).

Central exhaustion of cold gas and the growth of the bulge

The centre of PJ0116-24 exhibits characteristics of negative ‘feedback’. That is, star formation is suppressed due to the lack of cold gas within the central half kiloparsec, which could be attributed to past activities of an AGN^37,38 or a circumnuclear starburst^37,39,40. The central stellar mass surface density of PJ0116-24 is ${\varSigma }_{\ast , < 0.5\,{\rm{kpc}}} \approx3\times {10}^{10}\,M_{\odot }\,{{\rm{kpc}}}^{-2}$ within ~0.5 kpc, or ${\varSigma }_{\ast , < 1\,{\rm{kpc}}} \approx1.5\times {10}^{10}\,M_{\odot }\,{{\rm{kpc}}}^{-2}$ within ~1 kpc, estimated from its total stellar mass and assuming that the delensed HST 1.6 µm emission traces the Σ_* distribution. This is about 1/3 of Σ_* in the densest stellar systems of massive compact elliptical cores in the local Universe⁴¹. The cold molecular gas surrounding the centre at galactocentric radii ~1–2 kpc is abundant and dense, with ${\varSigma }_{\text{cold gas},1-2\,\text{kpc}} \approx$ $2\times {10}^{10}\,M_{\odot }\,{\text{kpc}}^{-2}$, and the SFR surface density there is ${\varSigma }_{\text{SFR},1-2\,\text{kpc}} \approx$ $100\text{--}200\,M_{\odot }\,{\text{yr}}^{-1}\,{\text{kpc}}^{-2}$, based on the delensed distribution of dust emission at ~2 mm (with a comparable resolution as CO). These Σ_{cold gas, 1-2kpc} and Σ_{SFR, 1-2 kpc} are more than ten times higher than those of most local spiral galaxies⁴² and are comparable to the maximal local gas and SFR surface densities predicted in secularly evolving galaxy simulations¹⁵. This picture is probably consistent with the formation path of compact, quenched bulges. Simulations have predicted a ‘wet disc compaction’ mechanism^43,44,45 in which disks fed by cold streams and undergoing violent disc instability will drive gas to inspiral towards the centre due to dissipative processes. Soon after a high central mass surface density is accumulated (${\varSigma }_{* , < 1\,\text{kpc}} \approx{10}^{9\text{--}10}\,M_{\odot }\,{\text{kpc}}^{-2}$)⁴⁴, internal quenching occurs, which consumes the cold gas inside-out and forms a massive, compact central bulge^25,31. PJ0116-24 appears to have been caught in the middle stage of this wet compaction, namely the inside-out quenching. This scenario is further supported by the supersolar metallicity of this galaxy, given that about half of the stars in the Galactic bulge have indeed supersolar metallicity⁴⁶.

A secular evolution pathway for massive HyLIRGs

PJ0116-24 is a unique hyperluminous rotating disk appearing as the brightest HyLIRG Einstein ring in the southern sky, and it is one of the rare cases studied in detail for several gas phases. Unlike almost all other extreme HyLIRGs, which are major mergers^{7,8,9,10,11,12,21,47,48}, PJ0116-24 does not obviously have massive companions or disturbed kinematics as evidence for major mergers. Recently, only one other lensed HyLIRG Einstein ring, named 9io9 or PJ020941.3 + 001559 (at z = 2.55, with μL_IR ≈ 1.3 × 10¹⁴ L_⊙ and μ ≈ 10, also in the PASSAGES selection)^4,49,50 has shown evidence of circular rotation, with ${v}_{\max }/{\sigma }_{0,{\mathrm{mol}}.\,{\mathrm{gas}}} \approx4.9\pm 0.7$ from the cold-gas tracer CO(4–3) observations at ~360 pc delensed scales⁵⁰. To date, these are the only two known HyLIRGs with evidence of circular rotation over ~10⁴ deg² in the sky. Compared to 9io9 (ref. ⁵⁰), PJ0116-24’s CO data have a slightly higher spatial resolution and reveal much richer substructures, and PJ0116-24’s near-IR IFU observations probe the heavily attenuated rest-frame optical lines that have rarely been imaged in other HyLIRGs.

The robust confirmation of PJ0116-24 as the most rotationally supported HyLIRG (${v}_{\text{rot}}/{\sigma }_{0,{\mathrm{mol}}.\,{\mathrm{gas}}}=9.4\pm 0.9$ and ${v}_{\text{rot}}/{\sigma }_{0,{\mathrm{ion}}.\,{\mathrm{gas}}}=$ $5.7\pm 0.6$) from this work is key evidence suggesting that secular evolution, that is, without recent major mergers, can be responsible for maximal star formation in the Universe (intrinsic $\text{SFR} >\text{1,000}\,{M}_{\odot }\,{\text{yr}}^{-1}$, in our case PJ0116-24 has $\text{SFR} \approx\text{1,500}\,{M}_{\odot }\,{\text{yr}}^{-1}$). This scenario has been well predicted in simulations before¹⁵ and is a decisive complement to the established paradigm of the evolution of massive (${M}_{* } > {10}^{10.5}\,{M}_{\odot }$), gas-rich, turbulent, Toomre-unstable, disk galaxies at z ≈ 1–3 (refs. ^{19,26,27,28,29}) studied at kiloparsec scales.

In our state-of-the-art three-dimensional (3D) dynamical modelling directly in the image plane (Methods), we corrected for the lensing, spatial beam-smearing, line instrumental broadening and 3D projection effects and robustly constrained that PJ0116-24 is genuinely a massive, turbulent, rotating disk with baryonic mass ${M}_{\text{disk},\text{dyn}} \approx{10}^{11.3\pm 0.3}\,{M}_{\odot }$ and intrinsic cold-molecular and warm-ionized turbulence dispersion ${\sigma }_{0,{\mathrm{mol}}. {\mathrm{gas}}} \approx43\pm 5\,\text{km}\,{\text{s}}^{-1}$ and ${\sigma }_{0,{\mathrm{ion}}.\,{\mathrm{gas}}} \approx69\pm 8\,\text{km}\,{\text{s}}^{-1}$. Moreover, it is one of the rarest HyLIRG Einstein rings in the Universe. This study presents a unique case that not only exemplifies the secular evolution pathway for HyLIRGs but also serves as the most detailed high-resolution, multi-tracer study of the HyLIRG population, which is nowadays routinely found in existing and future wide-area surveys.

Methods

VLT ERIS IFU Science Verification Program

Through the ERIS Science Verification Program (ID 110.258 S, PI D. Liu), we observed PJ0116-24 in the H- and K-bands in seeing-limited conditions with the high-spectral-resolution H_short and K_short gratings and the 250 mas pixel scale. The two gratings have a spectral resolution of R ≈ 11,200 and cover λ ≈ 1.46–1.67 µm (1.93–2.22 µm) in the H-band (K-band). The observations were taken on 2 December 2022 during the end of ERIS commissioning (https://www.eso.org/sci/activities/vltsv/erissv.html), as for the other ERIS Science Verification Programs, which aimed to verify the capabilities of this new instrument.

The observation blocks included a first acquisition to a nearby bright ‘cal_PSF’ star (Gaia DR3 5040854651181560064) and then a sequence of object exposures (J2000 right ascension 01 h 16 min 46.7792 s and declination −24° 37' 02.518″) and sky exposures (J2000 right ascension 01 h 16 min 49.1979 s and declination −24° 37′ 24.918″) with small dithering offsets (±0.5″) to allow for a better combined spatial sampling and to avoid the effects of bad or hot detector pixels. For each H- and K-band observation, we obtained eight object exposures and four sky exposures with a sequence of object–sky–object and object–sky–object dithering. The dither pattern of the object exposures was [(−0.5″, –0.5″), (+0.0″, −0.5″), (+0.5″, −0.5″), (+0.5″, +0.0″), (+0.5″, +0.5″), (+0.0″, +0.5″), (−0.5″, +0.5″), (−0.5″, +0.0″)]. The sky exposures pointed towards the same blank sky. Each exposure had 300 s of Detector Integration Time (DIT).

The latest ERIS pipeline was used to perform flat fielding, distortion correction, wavelength axis calibration, atmospheric OH line fitting and subtraction, sky exposure subtraction from object exposure data, and correction for differential atmospheric refraction. We calibrated the astronomical data unit to physical flux scale conversion based on ‘telluric’ stars observed with the same filter, spectral resolution and pixel scale as our observations. We collected all available telluric star observations from April 2022 to July 2023 and determined a mean physical flux scale conversion spectrum by normalizing to a telluric star’s theoretical black body (with temperature from stellar type) or template SED (see the European Southern Observatory (ESO) website; for example, https://www.eso.org/sci/facilities/paranal/decommissioned/isaac/tools/lib.html). The variation of the flux scale conversion spectrum was ~20%–50%, which could be due to the varying atmospheric conditions of observations not considered for in the pipeline (the conditions considered were air mass and DIT but not the seeing or water vapour).

After obtaining individual calibrated 300 s exposure data cubes, we combined each cube with a custom script that performs sigma clipping for bad pixels, image reprojection and weighted-averaging channel by channel. The weighting was based on the spectral pixel (spaxel) errors determined from sigma-clipped statistics. The final combined image cube was smoothed by a 2 pixel boxcar kernel to increase the S/N. The PSF was determined by fitting a 2D Gaussian to the ‘cal_PSF’ star, which was also smoothed by a 2 pixel kernel. The PSF FWHM was 0.7″ in the K-band data and 1.2″ in the H-band data.

Emission line spectra were extracted within a polygon region that corresponds to the full Einstein-ring area seen in the HST stellar emission when smeared to the low resolution (Supplementary Fig. 1). The uncertainties for the spectra were extracted correspondingly using the aforementioned error cube, as shown in red in Supplementary Fig. 1 (right). Spikes indicate channels contaminated by the atmospheric OH line (for example, ref. ⁵¹). Then the one-dimensional flux-calibrated spectra were fitted simultaneously with lines and a polynomial continuum, for each H- and K-band. The line flux errors were determined based on the fitting method, for which we used the astropy.modeling software package^52,53,54 with a fitting type of LevMarLSQFitter.

The Hα and [N ii] λλ6548,6584 lines were also fitted with a broad-line and a narrow-line component for each of the lines, as shown in Fig. 3b, due to their high S/N. We assumed that the broad-line components have the same linewidth and line-centre relative offsets, whereas the narrow-line components have a different linewidth and set of line-centre offsets. We found redshifted, broad-line components for all Hα and [N ii] λλ6548,6584 data points. The broad-line FWHM was ~510 km s⁻¹, about twice of that of the narrow-line FWHM, which was ~240 km s⁻¹. Instead of considering it as evidence of an outflow, we verified that it was mostly due to the asymmetric lensing magnification of the blueshifted and redshifted parts of the galaxy, as described in the main text. In our later rotating-disk forward modelling, we were able to recover a very similar asymmetric line profile in the image plane from a modelled pure Sérsic-profile rotating disk in the source plane. Therefore, the global asymmetric spectrum does not provide convincing information of an inflow or outflow. A future high-resolution near-IR IFU follow-up will shed light on any inflow or outflow signal for the ionized gas across the galaxy and at the AGN.

Then, we corrected the stellar Balmer absorption for the Hα and Hβ line fluxes using our SED models based on BC03, Bruzual and Charlot’s stellar population synthesis code (‘Fitting the SED’)⁵⁵. We used the high-resolution stellar library to generate a model SED with stellar Balmer absorption and measured the Hα and Hβ stellar absorption equivalent widths as 4.1 and 5.4 Å, respectively. These equivalent widths were multiplied with an emission filling fraction of 0.36 (0.23) for Hα (Hβ)⁵⁶, multiplied with the corresponding stellar continuum flux densities and then added to the directly measured line fluxes.

ALMA observations

We obtained several ALMA observing programmes (ID 2017.1.01214.S, PI M. Yun; ID 2019.1.01197.S, PI P. Kamieneski; and ID 2021.1.00353.S, PI K. Harrington) to observe the molecular gas CO line emission and dust continuum in the PASSAGES strongly lensed HyLIRG sample. The CO(3–2) data analysed in this study was combined with data from programme 2017.1.01214.S with 7.08 min on-source integration at 0.31″ angular resolution and with data from programme 2019.1.01197.S with 4.56 min on-source integration at 0.35″ resolution and 19.08 min on-source integration at 0.10″ resolution.

In this work, we calibrated the raw ALMA measurement sets with the standard calibration pipeline with the Common Astronomy Software Applications (CASA)⁵⁷. We then regridded the visibility data to a common spectral axis using the CASA mstransform task with a channel width of 30 km s⁻¹ and subtracted the continuum with the CASA uvcontsub task. The CO(3–2) image cube was produced by running the CASA tclean task, first with multiscale models with a threshold down to 3σ and then by adding hogbom (single-scale) models down to 1σ, where σ is the root mean square of the cube determined after a few major cycle iterations. A clean mask was applied based on the smoothed version of a first-run CO(3–2) image cube. This tclean approach has been found very useful for imaging galaxies with extended emission structures⁵⁸. We used the natural weighting, which gives the best S/Ns at a sufficiently small, synthesized beam (FWHM ≈ 0.19″ × 0.16″) for our analysis. The maximum recoverable scale was 7.5″, which was sufficient for recovering the extended emission from PJ0116-24. We also produced Briggs-weighting clean cubes, which have higher angular resolutions and better suppressed sidelobes but are noisier. Thus, we used them mainly for visually inspecting the emission peaks in our lensing modelling.

As our CO dataset includes several array configurations, the synthesized clean beam deviates from a 2D Gaussian shape. There is a factor of 2–3 difference between the effective area of the clean beam and the dirty beam, leading to an incorrect flux scale in the clean residual (that is, emission that was not contributed from any clean model component) and affecting the total flux measurements by a few tens of percent depending on the clean depth and the clean-to-dirty beam area ratio ε. Therefore, a correction was needed to scale down the flux within the clean residual by the factor ε, which is called the JvM correction⁵⁹. We applied this correction to all our imaged cubes by computing ${I}_{\text{clean image}}-\left(1-\varepsilon \right){I}_{\text{clean residual}}$, where ε is measured from the ratio between the clean beam area and the area of the dirty beam image integrated from the centre to the first null radius.

The ALMA dust continuum images at various bands were produced similarly and used to constrain the SED of PJ0116-24 (‘Fitting the SED’). We compared the dust continuum images and fluxes made by us and made by previous studies^5,6 and found good agreement.

VLT MUSE observations confirmed the lens redshift

Next we report VLT MUSE IFU observations around PJ0116-24, which spectroscopically confirmed the foreground lens to be a massive early-type galaxy (J2000 right ascension 01 h 16 min 46.7963 s and declination −24° 37′ 02.234″) at z = 0.554. The MUSE observations were obtained in a filler programme (programme code 110.23, PI F. Bian) under poor seeing weather conditions during 5 to 6 October 2022, with a seeing of 0.86″–1.80″. The MUSE data cover a wavelength range of 4,800–9,300 Å with a spectral resolution of R ≈ 1,770–3,590. The wide-field mode offers a field of view of 1′ × 1′ with a pixel scale of 0.2″. The data reduction used the standard MUSE data reduction pipeline provided by ESO (https://www.eso.org/sci/software/pipelines/; see B. Alcalde Pampliega et al., manuscript in preparation).

Supplementary Fig. 2 shows the MUSE spectrum extracted at this foreground lens within an aperture of 1.3″ and a fitted spectrum using the pyplatefit Python package’s fit_spec function⁶⁰ and with identified line names labelled. Strong absorption lines from old stellar populations are clearly visible in the spectrum, especially the prominent Ca doublet. Thus, we determined the redshift of the massive lens galaxy to be z = 0.554. A few more sources have a MUSE spectrum indicating that they are at z ≈ 0.554. Thus, they could be in the same galaxy group or cluster as the massive galaxy acting as the primary lens of PJ0116-24 (B. Alcalde Pampliega et al., manuscript in preparation). However, the closest second member of this potential galaxy group or cluster is at 8.2″ away from the massive lens galaxy, so that they have little impact on the lens modelling for PJ0116-24. We considered a dark matter halo plus an external shear in addition to the massive lens galaxy itself (see the next section). Thus, we took the effect of this potential galaxy group or cluster into account.

HST and Gemini optical and near-IR data

The HST WFC3 F160W data are from the observation programme GO-14653 led by PI J. Lowenthal (J. Lowenthal et al., manuscript in preparation) and ref. ⁶ has details of these observations. The image was dithered with five positions and combined to achieve a pixel size of 0.065″. The PSF is ~0.15″. The point-source 5σ depth is m_AB ≈ 27.4. There are no other HST or James Webb Space Telescope (JWST) observations for PJ0116-24 as at the time of writing.

For PJ0116-24, there is also Gemini multi-object spectrograph r′- and z′-band imaging from the Gemini south programme GS-2018A-Q-216 (PI J. Lowenthal). The Gemini imaging for the whole PASSAGES sample will be presented in O. Cooper et al. (manuscript in preparation). A summary has also been presented in ref. ⁶. The images have a pixel size of 0.16″ and the seeing-limited PSF is ~0.7″–0.8″.

Lens modelling

PJ0116-24 is primarily lensed by a foreground elliptical galaxy and its dark matter halo. We used the following steps to construct our lens model.

• (1) We corrected the HST astrometry using the Gaia DR3 star catalogue as a reference and then performed source fitting of the HST images. We manually added prior sources at HST emission peaks and used GALFIT for the prior source fitting with Sérsic models. The prior sources include both the foreground elliptical galaxy and the lensed PJ0116-24. As the foreground elliptical galaxy has strong and extended HST emission, we used two Sérsic models to model it best. Once we had obtained the best-fitting parameters, we ran GALFIT again but only to model the foreground elliptical galaxy without fitting to obtain the foreground-subtracted HST image only with PJ0116-24’s emission. The foreground-subtracted image is shown in Fig. 1 and Supplementary Fig. 3 (see also refs. ^5,6, and J. Lowenthal et al., manuscript in preparation where the original HST image is shown).

• (2) We identified emission peaks (‘knots’) in HST and ALMA channel CO(3–2) maps and used their sky coordinates as the input for our later lens model optimization. As demonstrated in ref. ⁶¹, which analysed a ×10 less luminous, massive, strongly lensed galaxy J0901-1814 in a similar way, it is very important to have both HST and ALMA knots as they provide highly complementary spatial coverage. The ALMA channel maps are critical because they provide unique velocity information for the pairing of knots. Supplementary Fig. 3 shows the knots identified for HST and CO emission peaks as well as the CO and dust intensity and CO velocity maps.

• (3) With the knot positions as the input, we used the glafic software^62,63 to optimize the lens model parameters. Our model set included a dark matter halo with a Navarro–Frenk–White profile⁶⁴, the foreground elliptical galaxy as a singular isothermal ellipsoid model and external shear to allow for better optimization (see also ref. ⁶¹). We ran the glafic default minimum chi-square algorithm to determine the best-fitting lens model parameters. The glafic input files and lens model parameters are publicly available with this work.

• (4) After obtaining the best-fitting lens model parameters with the minimum chi-square algorithm, we used a Markov chain Monte Carlo (MCMC) algorithm to explore the parameter space and to verify that the parameters were tightly constrained, following ref. ⁶¹. The MCMC fitting gives the probability distribution of each lens model parameter and verifies that the minimum chi-square lens model is within the 1σ confidence regions. The MCMC-based fitting results are also publicly available with this work. As discussed in ref. ⁶¹, MCMC fitting is useful for obtaining the probabilities and uncertainties, but the minimum chi-square lens model provides the best fits to the knots. Therefore, we used the minimum chi-square lens model as our best model for the further analyses. We visually verified the reconstructed source-plane maps to ensure the reliability of our best model (as can be seen in Figs. 1 and 2 for the different tracers). In addition, we also compared our delensed maps with those from Kamieneski et al.⁶ using the lenstool software and based on the HST and ALMA 0.4″ resolution 1.1 mm dust continuum maps. We found very similar delensed maps. Kamieneski et al.⁶ reported a magnification of μ_1.1mm ≈ 7 for the 0.4″ resolution dust continuum map, and their model also implies μ_HST ≈ 13.8 for the same 0.15″ resolution HST map as used in this work (private communication). Their μ_HST is, thus, within 13% of the μ_HST ≈ 15.9 in this work, which was based on the 0.15″–0.19″ resolution HST and ALMA CO channel maps.

After obtaining the best lens model, we converted the image-plane data into the source plane, namely delensing. We used the mesh grid from the best-fitting lens model, which describes how rectangle footprints (‘cells’) in the image plane are mapped into polygons in the source plane. We looped over mesh grid cells and co-added the emission from the image-plane pixels (in each cell’s rectangle footprint) into pixels in the source plane (in each cell’s polygon footprint) with linear interpolation inside each cell. For the CO and Hα cubes, we always first reconstructed the source-plane data cubes then extracted the line amplitude, velocity and velocity dispersion properties with one-dimensional Gaussian profile fitting for each spaxel in the cube, with either the scipy.optimize.curve_fit or pymc3 Python package. Note that the source-plane images were mostly for visualization and were not directly used for our kinematic modelling. Instead, the lens model’s mesh grid was used by our kinematic modelling software DysmalPy+Lensing to allow direct fitting of the image-plane data (‘Dynamical modelling with DysmalPy+Lensing’).

As previously mentioned, we measured the ‘delensed beam’, that is, the angular resolution in the source plane, by generating 2D Gaussian beams in the image plane at various positions then delensing each of them into the source plane and analysing their shapes. For image-plane 2D Gaussian beams (FWHM ≈ 0.16″) at positions A1 and A2, we measured a delensed beam FWHM of 0.088″ × 0.055″ and 0.149″ × 0.039″, respectively. Their minor-axis FWHM, thus, correspond to a physical scale of 450 and 320 pc, respectively, at z ≈ 2.125 (8.18 kpc arcsec⁻¹). Close to the critical curves, the best delensed angular resolution can reach ≤100 pc along the northern CO-bright arc, for instance, along the C2–C1–B4 knots in Supplementary Fig. 3. This happens to be the location of the CO arm-like substructures, which allowed us to spatially resolve these disk substructures even better than the centre, as shown in Fig. 1d.

In Fig. 1d, we also see some other Hα blobs having no HST or CO counterparts. The brightest one is north of clump B, as mentioned in the caption. We made tests to verify that this Hα blob is probably an artefact of delensing. The relatively poor angular resolution of our ERIS data created beam-smeared emission and noise fluctuations. Those that are scattered right across the critical curves, for example, around B2, were mapped to the exact position of the artificial Hα blob. Higher angular resolution ERIS observations with adaptive optics will probably solve this issue, as beam-smearing across critical curves will be largely reduced.

Comparison with Balmer decrements in the literature

Figure 3 compare the Balmer decrement of PJ0116-24 with the following galaxy samples:

• local SDSS galaxies from the SDSS DR8 MPA-JHU catalogue^65,66,67, with type ‘GALAXY’ and with S/N ≥ 10 in both Hα and Hβ line fluxes

• MOSDEF high-z galaxies from the MOSDEF survey⁶⁸, with S/N ≥ 4 in both Hα and Hβ line fluxes

• G165 Arc 1a, a HyLIRG that is also in the PASSAGES sample, and its companion galaxy NS 46; they are both in the G165 field⁶⁹

There are also Balmer decrement measurements for several hundred galaxies of different types, for example, local ULIRGs⁷⁰ and star-forming galaxies at z ≈ 0.75–7 (refs. ^{71,72,73,74,75,76}). Local ULIRGs are rare, merger-driven starbursts. A few known ULIRGs (Arp220, UGC05101 and IRAS12112 + 0305) have an extreme Hα/Hβ ≈ 9–20 (ref. ⁷⁰). Nevertheless, the majority of high-z galaxies with rest-frame optical line measurements do not have a global Hα/Hβ as high as that of PJ0116-24, and their stellar masses and SFRs are, in general, much lower. An exceptional case was reported recently by Frye et al.⁶⁹, G165 Arc 1a, a HyLIRG. This HyLIRG is also in the PASSAGES sample and has JWST NIRSpec multi-object spectroscopy indicating that Hα/Hβ ≈ 11.7 ± 3.1 (corrected for stellar Balmer absorption). It is probably less massive (intrinsic stellar mass ${M}_{* } \approx{10}^{10.0}\,{M}_{\odot }$) and much less magnified (μ ≈ 2.7)⁶⁹ than PJ0116-24. It also appears to be interacting with a companion galaxy named NS 46, which is 1,400 km s⁻¹ offset from the systematic velocity of G165 Arc 1a and is an even less-massive galaxy (intrinsic stellar mass of ~10^9.1 M_⊙; ref. ⁶⁹). Both G165 Arc 1a and PJ0116-24 are HyLIRGs with extremely high Hα/Hβ, but their masses and merger history are very different. Thus, both are unique cases for understanding extreme dust attenuation.

Ionization source and gas-phase metallicity diagnostics

With the detected nebular emission lines in PJ0116-24, we analysed the ionization source using a BPT diagram⁷⁷. Supplementary Fig. 4 shows the [N ii]-BPT⁷⁷ and [S ii]-BPT⁷⁸ diagnostics. For the [N ii]-BPT, we adopted the starburst-AGN classification curve from ref. ⁷⁹, the star-forming and composite star-forming division curve from ref. ⁸⁰ and the separation curve between Seyfert and LINER-type AGNs from ref. ⁸¹ (see also ref. ⁸²). For the [S ii]-BPT, we used the ref. ⁸² classifications.

Our galaxy falls into the star-forming/composite regime. That is, the emission lines are mostly powered by star-forming H ii regions at global scales. This diagnostic is conservative, as refs. ^83,84 and many studies have found that the classification lines for high-redshift galaxies need to move towards higher [N ii]/Hα and [O iii]/Hβ, which puts our galaxy further inside the star-forming regime.

Interstellar medium mass from CO and [C i] line modelling

The mass of the interstellar medium (ISM), which is considered to be molecular gas in this study, was measured for PJ0116-24 as $\log (\mu {M}_{\text{ISM}}/{M}_{\odot })={12.49}_{-0.02}^{+0.32}$ from an updated MCMC run of the radiative transfer modelling in Harrington et al.⁸⁵ (to be presented in K. Harrington et al., manuscript in preparation). We re-imaged all available ALMA CO and [C i] line data using the cleaning and JvM correction method as mentioned above in ‘ALMA observations’. This included the highest-resolution CO(3–2) data used in this study and lower-resolution CO(4–3), CO(7–6), CO(9–8), [C i](1–0) and [C i](2–1) data observed by ALMA programmes 2017.1.01214.S (PI M. Yun), 2019.1.01197.S (PI P. Kamieneski), 2021.1.00353.S (PI K. Harrington) and 2021.2.00062.S (PI D. Riechers). The single-dish CO(1–0), CO(3–2), CO(6–5), CO(7–6), CO(8–7), CO(9–8) and [C i](2–1) line fluxes reported in Harrington et al.⁸⁵ were also used, and these were consistent with our interferometric measurements when both exist.

The Harrington et al.⁸⁵ radiative transfer modelling assumes a turbulent gas density probability distribution function to describe the ISM. CO and [C i] molecular line emission was modelled by integrating the contributions from different gas densities. The gas kinetic temperature (T_kin) was modelled as a function of the H₂ density (${n}_{{\text{H}}_{2}}$). The parameters of the temperature–density function, CO and [C i] abundances, ${n}_{{\text{H}}_{2}}$, cloud radius properties and turbulence velocity dispersion ∆V_turb were all free parameters. Harrington et al.⁸⁵ reported $\log (\mu {M}_{\text{ISM}}/{M}_{\odot })=12.67\pm 1.15$ for PJ0116-24, with a relatively large uncertainty mostly from the limited CO and [C i] lines and the degeneracy among parameters in that study. Our updated modelling included more data so that the resulting M_ISM has a much tighter constraint.

There are still a few types of uncertainty that our modelling did not account for. For instance, the assumed monotonic form of the temperature–density function is probably too simple to account for the real ISM conditions, as there could be a range of gas temperatures for a given gas density. Indeed, by running a simpler two-component radiative transfer model following Liu et al.⁸⁶, we found a slightly smaller M_ISM ($\log (\mu {M}_{\text{ISM}}/{M}_{\odot }) \approx12.38\pm 0.11$, which is at the lower 1σ boundary of the updated Harrington et al. radiative transfer modelling) due to a warm and very dense gas component in addition to a cold and dense component. However, the two-component modelling achieves poorer fitting to the observed data, especially at mid-J CO lines. Thus, we used the updated Harrington et al.⁸⁵ radiative transfer modelling result.

The scenario of a giant star-forming clump

Clump B, which is labelled in Fig. 1, was detected in the HST and Gemini images, as well as in the ERIS Hα map and CO channel map but not in the dust continuum images. It is quadruply lensed into the B1–B4 counterimages. As discussed in Kamieneski et al.⁶, this lensing configuration indicates that the clump should be at the same redshift as the Einstein ring. The Hα line moment-0 map (Fig. 1d) clearly confirms this clump as all four counterimages were detected, and it indicates that it is forming stars. The CO(3–2) emission was also detected at a significance of ~2–4σ in the spectra extracted from the clump at the B1–B4 locations, as shown in Supplementary Fig. 6. The integrated CO(3–2) fluxes of the clump at the B1–B4 locations were 0.19 ± 0.10, 0.29 ± 0.10, 0.40 ± 0.16 and 0.23 ± 0.06 Jy km s⁻¹, respectively. They are in line with the magnification factors estimated for the B1–B4 locations (based on our best-fitting lens model magnification map smoothed to the CO(3–2) angular resolution), which are ~9–10 at B1 and B4 and ~14–15 at B2 and B3. Given that the total integrated CO(3–2) flux of PJ0116-24 was 66.4 ± 8.5 Jy km s⁻¹ from our new data reduction (which is consistent within ~1σ to the measurement 51.10 ± 10.22 Jy km s⁻¹ reported in Harrington et al.⁸⁵), the clump’s total CO emission is only ~1.6% of the total galaxy’s CO emission. If the CO emitting gas in the clump and across the PJ0116-24 disk has similar excitation conditions and, thus, a CO-to-H₂ conversion factor, then the CO line flux ratio is also a molecular gas mass ratio, indicating that the cold ISM mass is negligible in the clump compared to the entire galaxy.

The HST 1.6 μm emission of the clump is about 13% (1/7) of the whole PJ0116-24, but this is a rest-frame 5,000 Å light ratio instead of a mass ratio. Compared to the central kiloparsec of PJ0116-24, the clump is ~1 mag bluer in the z − H colour (roughly based on the Gemini z-band and HST H-band images, although their angular resolutions are very different). Based on the BC03 (ref. ⁵⁵) stellar population synthesis models, ~1 mag bluer z − H colour corresponds to a H-band mass-to-light ratio lower by a factor of 5 under solar metallicity (for example, based on the SED templates in Liu et al.⁸⁷; see also the appendix in ref. ⁸⁸). Thus, the clump probably has a stellar mass ratio of less than 1–3%, which also is negligible compared to the entire galaxy.

The [N ii] line map does not show detections at the four clump counterimages as clear as in the Hα line map, probably due to the relative line S/N. The [N ii]/Hα ratio could be useful for probing the metallicity of the clump, but realizing an accurate ratio is subject to a future high-resolution IFU follow-up study.

Given the aforementioned environment of PJ0116-24 and the detection of stellar light and cold and ionized gas in the clump, we believe that the origin of the clump is very probably in situ. It may have formed inside the massive turbulent gas disk due to gravitational instability and was fed by inflowing gas from the dark matter halo^{25,30,31,89,90,91,92,93,94,95,96,97,98,99,100,101,102}. Such Hα-bright, giant, star-forming clumps have been found in many massive unlensed galaxies at z ≈ 1–3, as studied with adaptive-optics-assisted near-IR IFU observations^{19,26,27,28,29,103,104} (see also the review by Förster Schreiber and Wuyts¹⁸). An ex situ scenario is less probable, as it requires the clump to be a low-mass (~10⁹ M_⊙) blue starburst with a high SFR but little dust while having a relative motion that is co-rotating with the disk of PJ0116-24.

Dynamical modelling with DysmalPy+Lensing

We performed 3D forward kinematic modelling that incorporates the lensing effect using the Dysmal software (refs. ^{19,26,27,28,105,106,107,108,109,110,111}) and the recent DysmalPy version¹¹² (see also refs. ^{61,104,113,114,115,116,117,118}) with the Lensing module developed by ref. ⁶¹. The advantage of this forward-modelling approach is that the beam-smearing and lensing are simultaneously taken into account. This avoids the complication of modelling a system to fit the delensed data that have a spatially varying PSF or inaccuracies resulting from approximating the variable resolution with a uniform value. We briefly describe the DysmalPy+Lensing forward modelling below (see ref. ⁶¹ for more details).

The source-plane galaxy model is a combination of two geometrically thick Sérsic-profile components, namely a disk and a bulge, and a dark matter halo with a Navarro–Frenk–White profile. The disk and bulge components represent the baryonic matter in the galaxy and have a systematic rotation and an intrinsic velocity dispersion. These two Sérsic-profile components may not physically represent a clearly distinct disk and bulge as in local galaxies, but this allows for a wide flexibility of the mass distribution model and, hence, the shape of the rotation curve. The 3D geometry of the disk + bulge is defined by an inclination angle (i) and a position angle (PA), which we determined from the delensed HST image and CO velocity field. The disk and bulge half-mass radii R_e,disk and R_e,bulge were initialized with values obtained from fitting the delensed HST + CO mass radial distribution with two Sérsic profiles (see, for example, ref. ⁶¹). Here, R_e,disk ≈ 3 kpc and R_e,bulge ≈ 0.6 kpc. The disk size is at the lower boundary of the size–mass relation of main-sequence star-forming galaxies¹¹⁹ and is in the regime for the transition to compact quiescent galaxies of the same mass. The disk has a scale height that is fixed to 1/4 of its scale length, which is a typical value for a thick disk at z ≈ 2 (refs. ^112,117). The bulge-profile mass component was found to contribute little to the kinematically fitted baryonic mass as the fitted bulge-to-disk ratio (B/T) was less than 0.1.

We performed the DysmalPy+Lensing fitting for a spatially and spectrally smoothed (0.4″) ALMA CO(3–2) data cube to increase the sensitivity to the extended, outer rotation curve. Supplementary Fig. 7 shows the 2D projected velocity maps of the data and the best-fitting model, which agree well within the uncertainty. The velocity dispersion in the data was generally fitted well in the disk but less well near the galaxy centre, where our idealized model does not account for any noncircular motion, which could increase the dispersion and also slightly twist the velocity field¹⁰⁴.

Our fitting method used MCMC to sample the high-dimensional parameter space. This approach allowed the disk + bulge mass M_{disk, dyn}, R_eff, B/T, intrinsic velocity dispersion σ₀, the dark matter halo virial mass M_{DM, dyn}, i and PA to vary (as well as allowing the geometric centre in the source plane to shift slightly). We applied a Gaussian prior on log(M_{disk, dyn}) (with a mean of ~11.5 and s.d. of ~0.35) guided by the photometrically measured ISM and stellar masses. We found that using a flat prior resulted in a slightly lower log(M_{disk, dyn}) within 1σ error. We also applied Gaussian priors on sin i and PA based on photometrically determined values.

The DysmalPy+Lensing fitting algorithm generated a 3D galaxy model and circular velocity per voxel in the source plane, which was projected onto the sky and deflected to the image plane. Then, it generated the effects of the PSF and the line spread function. It extracted the 2D velocity and dispersion fields, as for real data, and finally, it calculated the differences between the synthetic model data and the real data (velocity and dispersion fields). We ran the MCMC fitting with 250 walkers, each with 300 steps after about 30 burn-in steps. Supplementary Fig. 8 shows the MCMC posterior probability distributions of our fitting, with a few key parameters highlighted. All parameters were reasonably well constrained, and the residuals could not be further improved given our idealized mass model assumptions (introducing noncircular motion could be a future direction to explain the residuals seen in Supplementary Fig. 7; see ref. ¹⁰⁴).

Our obtained disk baryonic mass is ${M}_{\text{disk},\text{dyn}} \approx{10}^{11.29\pm 0.3}\,{M}_{\odot }$, which is not far from the sum of the ISM and stellar masses $({M}_{\text{ISM}}+{M}_{* }) \approx{10}^{11.51}\,{M}_{\odot }$, where ${M}_{\text{ISM}} \approx{10}^{11.26\pm 0.3}\,{M}_{\odot }$ and ${M}_{* } \approx{10}^{11.16\pm 0.3}\,{M}_{\odot }$. The excess of (M_ISM + M_*) compared to M_{disk, dyn} is at ≲ 1σ significance. This could be due to several assumptions or uncertainties unaccounted for, such as the star-formation history in our spectral energy fitting, the temperature–density coupling function in our ISM mass modelling, the initial mass function (IMF), which may change in high-z starbursts¹²⁰, and the contribution of the AGN to the continuum and line emission. All these uncertainties may add up to ~0.5 dex for the masses, so that the marginal difference between M_{disk, dyn} and (M_ISM + M_*) is not a significant issue in this study.

Fitting the SED

We ran multi-wavelength SED fitting including optical, IR, (sub)millimetre and radio bands to constrain the stellar and dust properties of the galaxy using the MiChi² code⁸⁷. The photometry data were mostly from refs. ^4,5,6, except that we reduced the ALMA 180, 260 and 395 GHz dust continuum consistently (Supplementary Fig. 3).

The SED model is a composition of a stellar component generated using the BC03 code⁵⁵ with a constant star-formation history and solar metallicity, a mid-IR AGN using empirical templates from ref. ¹²¹, a Draine and Li dust model heated by a physically driven interstellar radiation field¹²² (updated in ref. ¹²³) and a simple power-law synchrotron radio component. We assumed a Chabrier IMF¹²⁴. A bottom-heavier Salpeter IMF¹²⁵ would result in a stellar mass and SFR a factor of 1.7 larger. A top-heavier IMF would result in a smaller stellar mass and SFR. The MiChi² code⁸⁷ fitted the four components simultaneously and explored the multi-dimensional parameter space with Monte Carlo realizations to sample the chi-square-based probability distribution function of each parameter.

Supplementary Fig. 10 shows the SED of PJ0116-24 and the probability distributions of the apparent (magnified) stellar mass, stellar continuum dust attenuation E(B − V)_cont., IR luminosity and dust mass (Table 1). Adopting a usual stellar continuum to nebular dust attenuation ratio η_neb ≈ 0.5 (refs. ^{126,127,128,129}), the SED-based ${E(B-V)}_{\text{neb},\text{SED}}=1.00\pm 0.26$, which is fully consistent with the Balmer-decrement-based measurement reported in the main text.