{'fraction_non_alphanumeric': 0.039019006504795575, 'fraction_numerical': 0.030179907409216104, 'mean_word_length': 4.241744739671326, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A High-Resolution Investigation of the Multi-Phase ISM in a Galaxy during the First Two Billion Years

We have carried out the first spatially-resolved investigation of the multi-phase interstellar medium (ISM) at high redshift, using the z=4.24 strongly-lensed sub-millimetre galaxy H-ATLASJ142413.9+022303 (ID141). We present high-resolution (down to ~350 pc) ALMA observations in dust continuum emission and in the CO(7-6), H_2O (2_{1,1} - 2_{0,2}), CI(1-0) and CI(2-1) lines, the latter two allowing us to spatially resolve the cool phase of the ISM for the first time. Our modelling of the kinematics reveals that the system appears to be dominated by a rotationally-supported gas disk with evidence of a nearby perturber. We find that the CI(1-0) line has a very different distribution to the other lines, showing the existence of a reservoir of cool gas that might have been missed in studies of other galaxies. We have estimated the mass of the ISM using four different tracers, always obtaining an estimate in the range (3.2-3.8) x 10^{11} M_sol, significantly higher than our dynamical mass estimate of (0.8-1.3) x 10^{11} M_sol. We suggest that this conflict and other similar conflicts reported in the literature is because the gas-to-tracer ratios are ~4 times lower than the Galactic values used to calibrate the ISM in high-redshift galaxies. We demonstrate that this could result from a top-heavy initial mass function and strong chemical evolution. Using a variety of quantitative indicators, we show that, extreme though it is at z=4.24, ID141 will likely join the population of quiescent galaxies that appears in the Universe at z~3.

lengths in the early universe. SMGs are typically highly opticallyobscured by dust that originates from low/intermediate mass stars and supernovae (e.g., Rowlands et al. 2014) in regions undergoing extreme rates of star formation. These stellar nurseries are most abundant at redshifts around = 2−4 when the Universe was most rapidly forming the galaxies and clusters we see around us today. SMGs have therefore become a valuable probe of the physics of the growth of these structures, particularly at the more massive end of the galaxy mass function (e.g., Casey et al. 2014;Hodge & da Cunha 2020).
The most extreme SMGs are those discovered in the wide-area surveys carried out with Herschel (Negrello et al. 2017;Bakx et al. 2018), the South Pole Telescope (Vieira et al. 2013;Reuter et al. 2020) and Planck (Cañameras et al. 2015), including ones that have star-formation rates many thousands of solar masses per year, which have been found out to a redshift of ∼ 6 (e.g., Fudamoto et al. 2017). An obvious question is, what do these extreme objects evolve into at more recent epochs? In particular, do these galaxies evolve into the population of quiescent galaxies that appear in the Universe at ∼ 3? (e.g., Glazebrook et al. 2017;Girelli et al. 2019;Valentino et al. 2020;Forrest et al. 2020b,a;Merlin et al. 2019).
These huge star-formation rates must be fueled by large amounts of gas and many studies have found that SMGs house substantial molecular gas reservoirs of masses generally in the range 10 10 − 10 11 M (e.g., Ivison et al. 2011;Bothwell et al. 2013;Huynh et al. 2017;Dye et al. 2018;Kaasinen et al. 2019). The bulk constituent of this gas, molecular hydrogen, is not rotationally excited at the typical temperatures in the molecular phase in galaxies. As such, investigations have been forced to use tracers, either dust or a variety of molecular and atomic species, thus requiring knowledge of the gas-to-tracer ratio. There is a complex hierarchy of calibrations but most of the values assumed for these ratios are ultimately based on observations in the Milky Way, the only place where it is possible to measure both the luminosity in the tracer and the mass of the underlying gas (Dunne et al. 2021).
As for galaxies in the nearby Universe, the most frequently used tracer of the molecular phase of the ISM, because it has so many bright spectral lines in atmospheric windows, is the CO molecule (Tacconi et al. 2018;Birkin et al. 2021). As a tracer of the underlying gas, CO has some well known disadvantages (Bolatto et al. 2013), in particular the fact that the J=1-0 line, the best for tracing the coldest gas, is optically thick. In addition, observations with Fermi (Abdo et al. 2010), Planck (Planck Collaboration et al. 2011) and Herschel (Pineda et al. 2013) all suggest that approximately one third of the molecular gas in the Galaxy does not emit CO, probably because of photo-dissociation of the CO molecule (Planck Collaboration et al. 2011). For high-resolution studies of the ISM in high-redshift galaxies there is the additional problem that the Atacama Large Millimetre Array (ALMA) is only able to observe the higher-excitation CO lines, making it possible that these studies have missed much of the gas (although lower-frequency ALMA bands will become available in the near future).
The other bright molecular lines in high-redshift galaxies are those of the water molecule, which have been used to observe both SMGs Spilker et al. 2014;Yang et al. 2016;Jarugula et al. 2019) and AGN (e.g. van der Werf et al. 2011). The water lines are excited by far-IR photons (compared to collisional excitations of CO by H 2 ), explaining the correlation between line luminosity and far-infrared luminosity over three orders of magnitude (Yang et al. 2013). The water lines are therefore probably more useful as a tracer of star formation than the ISM (Jarugula et al. 2019).
Arguably the best tracers of the molecular phase of the ISM in highredshift galaxies are the fine-structure lines of atomic carbon, despite being generally fainter. The [ C ] ( 3 1 − 3 0 ) and [ C ] ( 3 2 − 3 1 ) lines have similar excitation temperatures and critical densities to the low-J CO lines, but they are optically thin and they are emitted in all molecular regions including CO-dark regions Papadopoulos & Greve 2004;Dunne et al. 2021).
The high-resolution configurations of ALMA, plus the magnification (often greater than a factor of 10) provided by gravitational lensing, have huge potential for investigating the physical conditions in the ISM in galaxies in the early Universe. This potential has been demonstrated quite dramatically using observations of the gravitationally-lensed SMG SDP81 that were carried out in the Science Verification Phase of the ALMA long baselines (ALMA Partnership et al. 2015). After being corrected for gravitational lens distortion, the reconstructed spectral-line and continuum images had a resolution of ∼50 pc (Dye et al. 2015;Rybak et al. 2015a,b). The images revealed a smoothly-rotating but clumpy disk, with a Toomre Q-parameter suggesting that the disk is in a state of collapse. The resolution, approximately the size of a giant molecular cloud, made it possible to compare the Larson relations for the clumps with those in the Milky Way (Swinbank et al. 2015).
The drawback of the observations of SDP81, however, is that they were in high-J CO lines, making it uncertain whether a cool molecular phase had been missed; it is possible that these high-J CO clumps simply trace regions in the disk where the gas is particularly warm rather than genuine agglomerations of gas. In addition, at the redshift of SDP81, several spectral lines do not fall within atmosphericallytransparent windows. To circumvent these shortfalls, we have therefore embarked on a second high-resolution study of the ISM using a galaxy at a redshift such that all of the key spectral lines are visible to ALMA.
One of the brightest strongly-lensed SMGs, H-ATLAS J142413.9 + 022303, hereafter referred to as 'ID141' following its label assigned by Cox et al. (2011) was discovered in the Herschel Astrophysical Terrahertz Large-Area Survey (Eales et al. 2010;Cox et al. 2011;Negrello et al. 2017). ID141 is not only bright because it is strongly-lensed but also because it is intrinsically very luminous. After correction for the lens magnification factor of ∼ 6, the starformation rate in ID141, estimated from the bolometric luminosity of the dust (and assuming the initial mass function of Kroupa & Weidner 2003), is 2400 M yr −1 , five times higher than SDP81. With a redshift of = 4.24 (Cox et al. 2011), ID141 is also one of the most distant strongly-lensed SMGs. At this redshift, both [ C ] lines, the bright [ C ] line at 158 m, multiple CO lines (including the CO(1-0) line in the near future) and water lines, and many other lines are visible to ALMA.
As a result of its interest, ID141 has been observed many times in the submillimetre/millimetre wavebands (Cox et al. 2011;Bussmann et al. 2012Bussmann et al. , 2013Omont et al. 2013;Dye et al. 2018;Cheng et al. 2020). Lens modelling by Dye et al. (2018, D18 hereafter) and Enia et al. (2018) revealed evidence of a disturbed morphology in the reconstructed continuum emission of the source (see also Bussmann et al. 2012 who showed that significant image residuals remain after subtracting the best fit lensed image of a single smooth source). A relatively high level of ionisation in the source was inferred from a high far-IR luminosity to gas mass ratio measured by Cox et al. (2011). This, along with very large line widths of ∼ 800 km s −1 observed in CO emission led the authors to suggest the SMG is a merger, like most other SMGs (e.g., Engel et al. 2010).
In this paper, we present high-resolution (∼ 0.1 arcsec) ALMA observations of ID141 in bands 3 and 4, complementing the ALMA band 7 observations presented in D18. The new observations cover the emission lines CO(7-6), H 2 O(2 1,1 − 2 0,2 ), [ C ] ( 3 2 − 3 1 ) and [ C ] ( 3 1 − 3 0 ), which have for the first time made it possible to carry out a high-resolution investigation of the ISM in a high-redshift galaxy that includes both the warm and the cool gas. The spectral-line observations have also allowed us to make a detailed kinematic analysis of this system, in particular addressing the question of whether ID141 is a merger.
One of the key findings of this paper is that the dynamical mass of ID141 is approximately four times lower than the gas mass inferred from Galactic gas-to-tracer ratios. SMGs with inferred gas masses that are equal to or larger than their dynamical mass have been reported in the literature for over a decade (e.g. Tacconi et al. 2008;Hodge et al. 2012;Bothwell et al. 2013). More recently, there has been an acceleration in the prevalence of such conflicting systems (e.g. Dye et al. 2015;Yang et al. 2017;Calistro Rivera et al. 2018;Neeleman et al. 2020;Mizukoshi et al. 2021). A common approach to reconcile these anomalies is to reduce the conversion factor (based on local Universe observations) that scales CO line intensity to total gas mass. In this work, we find that four different tracers yield the same excessive gas mass. It is unlikely that all four tracer conversion factors would need to be scaled by the same amount and so we have sought a more natural explanation in the form of a chemical evolution model which solves the discrepancy between the gas mass and dynamical mass in ID141.
The layout of this paper is as follows. In Section 2 we describe the acquisition and the reduction of the ALMA data. Section 3 presents our lens modelling and source reconstruction. Section 4 describes our investigation of the basic properties of the source, including the kinematic modelling and estimates of the mass of the ISM in the galaxy made five different ways. We discuss our findings in Section 5 and introduce a chemical evolution model of ID141. Section 6 is a summary of both our basic observational results and inferences from these results. Throughout this work, we have adopted the flat cosmological model with parameters ℎ = 0.678, Ω m = 0.308 and Ω Λ = 0.692 (Planck Collaboration et al. 2016).

DATA ACQUISITION AND PROCESSING
The ALMA data analysed in this paper were acquired under two separate ALMA programmes. The first, 2016.1.00450.S, obtained band 4 spectral imaging on 20th August 2017 in four spectral windows (SPWs) centred on the frequencies 141.2, 143.1, 153.3 and 154.9 GHz, each with bandwidths of approximately 2 GHz. The total on-source integration time was 2360s using 44 12m antennas giving a spatial resolution of 0.15 arcsec, a maximum recoverable scale of 2.1 arcsec and a sensitivity of 0.34 mJy beam −1 assuming 10 km s −1 channels.
The second programme, 2017.1.00029.S, obtained band 3 spectral imaging on the 14th and 16th of November 2017 when 46 and 44 12m antennas were used respectively. Combining the measurement sets from both dates, the resulting images have a spatial resolution of 0.11 arcsec, a maximum recoverable scale of 2.0 arcsec and a sensitivity of 0.56 mJy beam −1 assuming 10 km s −1 channels. These observations were again carried out in four SPWs with central frequencies of 93.9, 95.8, 105.8 and 107.7 GHz with bandwidths of ∼2 GHz and with a total on-source integration time of 3560s. This configuration of SPWs was chosen specifically to target the emission lines of CO, H 2 O and [ C ] (see below). Table 1 gives more specific details on the spectral set up.
Calibration of both the band 3 and band 4 data was carried out using the scripts provided by the ALMA science archive using version 5.1.1 for band 3 and version 4.7.2 for band 4 of the Common Astronomy Software Applications (CASA) software package (Mc-Mullin et al. 2007). Bandpass and phase calibration was provided via the calibrators J1337-1257 and J1410+0203 respectively. For band 3, the two measurement sets (MSs) were independently calibrated and then combined into one MS using CASA. We manually verified that no erroneous visibilities remained after calibration using CASA's plotms tool. All observed images of ID141 displayed in this paper were produced by cleaning these MSs using CASA with natural weighting and a pixel scale of 0.024 arcsec. To isolate emission lines (see Section 2.1), we used the uvcontsub tool within CASA for each SPW to fit a linear spectrum to channels not containing line emission and then subtract this fit from the full visibility data. A separate set of visibilities was created for analysis of the continuum emission by excluding channels containing line emission. In addition to the band 3 and 4 ALMA data detailed above, we have also used the ALMA band 7 data acquired under the ALMA programme 2013.1.00358.S and analysed by D18 (see Fig. 1). The observations in this programme amount to a total on-source integration time of 151s within four spectral windows spanning 336 to 352 GHz. The spectral windows contain no detectable emission lines and so we use the resulting 0.12 arcsec resolution continuum image solely for the purposes of morphological comparison of the reconstructed lensed SMG. We refer the reader to D18 for further details regarding these observations. Upon re-cleaning this dataset, we identified an error made in D18 such that the flux was over-estimated by approximately 30 per cent. The updated flux measured from the lensed image is 91 ± 5 mJy.

Emission lines
We identified four emission lines: H 2 O(2 1,1 − 2 0,2 ) ( rest = 752.033 GHz), CO(7-6) ( rest = 806.652 GHz), [ C ] ( 3 2 − 3 1 ) ( rest = 809.342 GHz, referred to as [ C ] (2-1) hereafter) and [ C ] ( 3 1 − 3 0 ) ( rest = 492.161 GHz, referred to as [ C ] (1-0) hereafter). Observed images of these lines are shown in Fig. 2 and the corresponding observed lensed line profiles (measured within a mask determined by selecting significant pixels) are shown in Fig. 3. Table 2 lists the four line fluxes measured in the observed (i.e. lensed) images. The CO(7-6), [ C ] (1-0) and [ C ] (2-1) lines were also detected by Cox et al. (2011). The line flux measured for the [ C ] (1-0) line measured by Cox et al. (2011) is consistent with our measurement, but the Cox et al. fluxes of the other two lines are nearly a factor of two lower than those measured by ALMA. Although the discrepancies are larger than expected, the errors quoted in Cox et al. (2011) do not include any uncertainty due to continuum subtraction which are large enough to eliminate any significant disagreement (R. Neri, private communication). Fig. 3 also shows the emission lines determined from the reconstructed source plane emission (see Section 3.2). In this case, we extracted the line emission from the source plane region that, according to the lens model, maps to the image plane mask used to determine the image plane emission. It is clear from the three lines with the highest signal-to-noise ratio (SNR), i.e. H 2 O, CO(7-6) and [ C ] (2-1), that emission lines extracted from the image have slightly stronger flux at negative velocities whereas the reverse is true in the lines extracted from the source plane. Our lensing model confirms that this is due to differential magnification; the region of the source responsible for the emission of the negative side of the lines lies closer to the lens caustic where the magnification is higher.
Taking the median observed frequency of each of the H 2 O, CO(7-6) and [ C ] (2-1) lines measured from the reconstructed source and calculating the mean of these weighted by their SNR results in an inferred source redshift of = 4.24417 ± 0.00003 (see Table 2 -note that we did not include the [ C ] (1-0) line in this measurement due to its significantly lower SNR). This measurement is consistent with the redshift of = 4.243 ± 0.001 determined by Cox et al. (2011). Table 2. Emission lines detected in the band 3 and band 4 ALMA data. The rest-frame frequencies, rest , are measured from the reconstructed source cubes and therefore eliminate differential magnification of the emission lines. The SNR-weighted mean of the redshifts implied by the H 2 O, CO(7-6) and [ C ] (2-1) lines is = 4.24417 ± 0.00003. The final column gives the line flux in the observed (i.e. lensed) image.

LENS MODELLING
To recover an undistorted image of the lensed SMG, we applied the modelling process of D18. Using a given model for the distribution of mass in the foreground lensing galaxy, this process reconstructs a pixelised map of the source surface brightness (see Warren & Dye 2003). In this paper, we modelled solely the interferometric visibilities. The process works by computing a set of model visibilities from the lensed image of each source pixel. We computed these visibilities using uvmodel contained within the MIRIAD software suite (Sault et al. 2011). The linear combination of the sets of model visibilities that best fits the observed visibilities then recovers the source surface brightness distribution for a given lens model. The lens model is optimised iteratively, repeating the linear inversion with every iteration of the lens model. In this work, lens model optimisation was performed with NeuralNest (Moss 2020) which uses a neural network to generate highly efficient proposals for Markov Chain Monte Carlo sampling.

Lens Model
Following the lens modelling in D18, we used an elliptical powerlaw density profile for the foreground lens, using the form given by Rest-frame frequency (GHz) CI(1-0) Figure 3. Emission lines, H 2 O(2 1,1 − 2 0,2 ), CO(7-6), [ C ] (2-1) and [ C ] (1-0), measured within a masked area determined by selecting significant pixels in their cleaned observed images (grey) and the corresponding line profile in their reconstructed source plane cubes (black). Image plane line profiles have been corrected for total lens magnification. Differential magnification causes the lensed image profiles to contain more flux at negative velocities even though intrinsically the source emits more flux at positive velocities. Kassiola & Kovner (1993) which has a surface mass density where 0 is the normalisation surface mass density and is the powerlaw index of the volume mass density profile. The elliptical radial co-ordinate,˜, is defined in terms of the Euclidean co-ordinates, ( . ) relative to the ellipse centre, and the elongation, (i.e., the ratio of semi-major to semi-minor axis length), by˜2 = 2 + 2 / 2 . The orientation of the semi-major axis, is measured east of north and the co-ordinates of the centre of the lens in the image plane are ( , ). Initial lens modeling tests with a model including external shear were attempted but it was found that negligible shear was required in the best fit model. This agrees with the findings of D18 and hence a parameterisation of external shear was not included in the final lens model which comprised six parameters in total. We also tried a lens model comprising two separate mass components as used by Bussmann et al. (2012) and as in D18, we found no significant improvement in the fit compared to a single ellipsoidal power-law. Although the two models give rise to a difference in the total source plane magnification of approximately 5 per cent, this is negligible in comparison to the errors on the intrinsic source properties we derive (see Section 5.2 for further discussion). Additionally, the two models yield almost identical source morphology and kinematical properties.
To determine our optimal lens model, we simultaneously fitted the visibilities from all four SPWs in the band 4 continuum data from this paper and the band 7 visibilities from D18. Combining bands and including the continuum in this way maximises the signal-tonoise ratio (SNR) to give the best possible lens model. We used the band 7 visibilities as presented in D18 but we time-averaged the more numerous band 4 visibilities into bins of 20 s to improve computational efficiency. We found that the size of this temporal bin can be changed by a factor of several without making any noticeable difference to the reconstructed lens model or source. We obtained a model that improves upon the model presented in D18 when applied solely to the band 7 visibilities. We attribute this to a combination of the new optimisation strategy (i.e., NeuralNest) used in this work and the complicated parameter space that results from randomising the source plane pixelisation (see Section 3.2).
The parameters of our best-fit lens model are: 0 = (0.60 ± 0.01) × 10 10 M kpc −2 , = 1.14±0.02, = (96±2) • , = 2.01±0.03. Comparing against the model in D18, there are two significant changes; the elongation has increased from 1.07 to 1.14 and the lens orientation has increased from 85 • to 96 • east of north. The source structure qualitatively looks very similar (see Fig. 1) but the total magnifica-tion of the source continuum has fallen from 6.4 to 6.1 in band 7. The Einstein radius of 0.97 ± 0.04 arcsec remains the same.

Source Reconstruction
We reconstructed the lensed source on a 0.8×0.8 arcsec source plane consisting of 400 Voronoi pixels formed using the method described in D18. The Voronoi tessellation is formed within the source plane region traced by mapping the region enclosed within an annular mask applied to the ring images as can be seen in Figs. 1 and 2. When optimising the lens model, we computed seven randomisations of the source plane pixelisation for each trial set of lens model parameters and returned the median likelihood of these seven randomisations to NeuralNest. This alleviates the effect reported in Nightingale et al. (2018) whereby lens model optimisation can occasionally stagnate when exceptionally good pixelisations occur by chance.
When reconstructing source cubes, we extracted the set of visibilities for each channel and applied the best-fit lens model determined in Section 3.1. We reconstructed the visibilities corresponding to each slice of each cube 50 times, randomising over source plane pixelisations every time. All 50 random Voronoi pixelisations were then individually binned onto a fine grid of 100 × 100 square pixels across the source plane and their mean computed to arrive at the final source reconstruction for a given cube slice.
The reconstruction scheme applied to cube slices results in an effective point spread function (PSF) that varies with lens magnification and therefore with position in the source plane. To quantify the resolution of our reconstructions, we computed the spatial variation of this PSF in the following way. First, we generated a set of model visibilities of a lensed image of a point source created using the best fit lens model. These visibilities were obtained using MIRIAD's uvmodel task which evaluates them at the same uv co-ordinates as the observed visibilities. We then applied our reconstruction method, including averaging the 50 random Voronoi source planes as described above to arrive at a reconstructed point source. This procedure was repeated for a grid of 3 × 3 point source positions in the source plane (arranged to cover the actual observed source light) and the resulting superposition of these recovered point sources is shown in Fig. 4. The dimensions of the PSFs obtained in this way for band 3 and 4 are very similar and these are 0.17 × 0.10 arcsec, 0.08 × 0.05 arcsec and 0.12 × 0.08 arcsec, in the south-eastern, north-western and middle of the source respectively where the dimensions quoted are full widths at half maximum (FWHM) along the major and minor axes. At the source redshift (where 1 arcsec corresponds to a physical scale of 6.90 kpc), the smallest physical scale we are able to probe is therefore approximately 350 pc. The median resolution across the source, For each point source, we made a model set of visibilities and reconstructed the set with the same source pixelisation scheme used for reconstructing the lensed SMG (see main text for more details). The beam at the centre of the grid has a FWHM of 0.12 and 0.08 arcsec along its major and minor axis respectively. For the band 3 data, the PSF shows the same qualitative variation across the source plane but has a linear scale that is approximately 15 per larger than the band 4 PSF.
computed as the geometric mean of the PSF FWHMs, is approximately 550 pc.
To determine the unlensed continuum flux of ID141, for each band, we applied the source reconstruction method to the emission line-subtracted visibilities obtained as outlined in section 2. The reconstructed source continuum in bands 3, 4 and 7 is shown in Fig. 1. Note that the same morphological features seen in the band 7 continuum source shown in this figure are also seen in D18, despite the different lens model and the different source pixelisation scheme. The quantitative difference is a 5 per cent reduction in magnification compared to D18 which results in a visually imperceptible increase in the linear scale of the source shown in Fig. 1 of 2-3 per cent. To determine the unlensed fluxes of the four identified emission lines, we reconstructed the source using the visibilities containing line emission, for each line. These are shown in Fig. 2.
With the source reconstructed in this way, the flux was then measured in the source plane aperture that maps to the image plane annulus for both the continua and emission lines. Line fluxes are given in Table 4 and the continuum fluxes measured in this way are 0.26±0.04 mJy, 1.16±0.14 mJy and 14.9±0.8 mJy in band 3, 4 and 7 respectively. The continuum and line magnifications were computed in more restricted apertures. We formed an image-plane mask by selecting all pixels with a significance of at least 3 and then mapped this to the source plane. Magnifications were then determined as the ratio of the image flux to source flux in these apertures. Table 3 lists the magnifications obtained in this way for the emission lines and band 3, 4 and 7 continua.

THE PROPERTIES OF THE LENSED SOURCE
In this section we give consideration to the reconstructed morphology of ID141, its spectral energy distribution (SED) including the inferred quantities of dust mass, dust temperature and star formation rate and the gas mass of the ISM. Section 4.6 presents our kinematical analysis. We discuss the implications of these results in Section 5.

Source Morphology
The reconstructed emission line images of ID141 in Fig. 2 show clear morphological differences to the reconstructed continuum emission in Fig. 1. Qualitatively, the line emission is more irregular than the continuum emission. This is partly an effect of the lower SNR of the reconstructed line flux but not entirely; the lensed emission line images in the leftmost column of Fig. 2 clearly have higher noise than the lensed continuum images in Fig. 1 but many of the clumpy features seen in the reconstructed line images in the rightmost column of Fig. 2 are multiply-imaged in the observed ring and can be seen at different wavelengths. One particularly striking feature is an offset in the peak emission seen in the [ C ] (1-0) line relative to the other lines and relative to the continuum emission. This is apparent in both the reconstructed source and the lensed image. The bottom-left panel of Fig. 2 illustrates this offset by over-plotting contours of the band 3 continuum flux on the line emission map. Since [ C ] (1-0) has a lower excitation temperature than the other lines, this suggests the presence of cooler gas in the southern parts of the galaxy. We return to this point in Section 4.3.

Far-IR Spectral Energy Distribution
Fig. 5 plots measurements of the continuum flux of ID141 at different rest-frame wavelengths. The filled points are taken from D18 but reduced by 5 per cent to account for the updated lens magnification and from Cheng et al. (2020). We also updated the 850 m flux originally taken from Bakx et al. (2018) with the flux of 90 mJy following Bakx et al. (2020). The new band 3 and 4 continuum fluxes are shown in this figure with empty circles.
Following D18, we fitted the rest-frame photometry with a modified black-body SED (referred to hereafter as the 'optically thick SED' ) comprising a single dust temperature component as well as an optically thin SED with two dust temperature components. This is an over-simplistic treatment of the expected continuous range of dust temperatures which we incorporate in later SED fitting (see Section 4.5). Nevertheless, these two simple SED models provide a good estimate of the range of possible dust masses. We measured the dust mass predicted by each SED using the method described in Dunne et al. (2011). For this, we took the observed ALMA 880 m flux density of 91 ± 5 mJy from the band 7 data and a dust mass absorption coefficient at the corresponding rest-frame wavelength, 168 m, of 1.7 m 2 kg −1 from James et al. (2002). Both SEDs were fitted with a fixed emissivity index of 2.0 (Smith et al. 2013). The open circles are the band 3 and band 4 continuum fluxes from the present work. The solid grey and black lines show an optically thin SED with two temperature components and an optically thick single-temperature SED fit respectively (see text for more details). The dashed lines are the two separate components that comprise the optically thin SED.
We allowed the temperature and normalisation of both components to vary in the fit for the optically thin SED giving four degrees of freedom. For the optically thick SED, temperature, normalisation and the 100 m-opacity, 100 , were varied in the fit, giving three degrees of freedom. The best-fit parameters are warm = 60 ± 6 K and cold = 34 ± 3 K with a ratio of cold-to-warm dust mass of 74 ± 13 for the optically thin SED. For the optically thick SED, we found = 57 ± 1 K and 100 = 4.8. The predicted dust masses are (1.4 ± 0.3) × 10 9 M and (7.4 ± 0.8) × 10 8 M for the optically thin and optically thick SEDs respectively.
Integrating these SEDs from 8-1000 m gives a far-IR luminosity, FIR , of 1.6 × 10 13 L and 1.5 × 10 13 L for the optically thin and thick SEDs respectively. These values are consistent with Cox et al. (2011) once corrected for lens magnification and places ID141 in the category of hyper-luminous infra-red galaxies (Rowan-Robinson 2000). Taking the mean of these luminosities and converting it to a star formation rate (SFR) according to the prescription given by Kennicutt & Evans (2012), which uses the initial mass function given by Kroupa & Weidner (2003), gives a value of 2400 ± 500 M yr −1 . This is consistent with the SFR calculated in D18 but is 5 per cent higher due to the 5 per cent reduction in magnification of the new lens model. As discussed in D18, given its dust mass, this places ID141 significantly above the range of SFRs for < 0.5 H-ATLAS galaxies and in the upper envelop of high redshift SMGs (Rowlands et al. 2014).

The Temperature Distribution in the ISM
There is a significant offset of the peak emission seen in the [ C ] (1-0) line compared to the peak of emission seen in the other three lines. Initially concerned that this may simply be the result of an offset in the world co-ordinate system between the band 3 data (where the [ C ] (1-0) line is located) and band 4 data, we confirmed that the continuum images between both bands align perfectly. Indeed, the [ C ] (1-0) emission is offset significantly from the band 3 continuum (see Fig. 2) and this is clearly apparent in both the reconstructed The advantages of the [ C ] lines for tracing the cool ISM over the traditional CO method has become increasingly recognised in the last few years Dunne et al. 2021, Dunne et al. in preparation ). The [ C ] lines have similar excitation temperatures and critical densities to the low-J CO lines, but their advantages over the CO lines are that the [ C ] lines are mostly optically thin and atomic carbon is present everywhere in the molecular phase . In contrast, a significant and variable fraction of the molecular phase in galaxies appears to contain no CO, presumably because of photo-dissociation (Abdo et al. 2010;Pineda et al. 2013;Planck Collaboration et al. 2011;Bigiel et al. 2020;Dunne et al. 2021).
Since the [ C ] lines are optically thin, it is reasonably straightforward to use the ratio of the [ C ] (2-1) and [ C ] (1-0) lines to investigate the variation of the temperature in the ISM. Therefore, we divided the [ C ] (2-1) and [ C ] (1-0) source images in Fig. 2 into a stack of five slices running east-west, each 0.2 arcsec long and 0.08 arcsec wide, with the stack positioned to cover the main [ C ] emission. We calculated the ratio of the luminosity (measured in K km s −1 pc 2 , Solomon & Vanden Bout 2005) of the two lines, = (2 − 1)/ (1 − 0), in each of these slices. Assuming local thermodynamic equilibrium, we then converted the line ratio into an excitation temperature using T ex = 38.8K/ln(2.11/ CI ) (Stutzki et al. 1997). Fig. 6 shows the variation of temperature in the stacked slices moving northwards across ID141. A 2.5 increase in temperature is seen rising from 40 K in the three southern bins to 100 K in the two northern bins. This trend is in accord with the other images in Fig. 2 because the water line and the CO(7-6) line also peak in the north, and both lines are likely to be from high-temperature gas (Papadopoulos et al. 2012a;Liu et al. 2017). The excitation temperature of the galaxy as a whole, estimated from the global fluxes for the [ C ] lines (Table 4) and the equation above, is 33 K. To put this into context, in the sample of 49 galaxies of Dunne et al. (in preparation), ID141 would rank as the fifth highest [ C ] excitation temperature.
The mass of the molecular phase of the ISM is proportional to (1 − 0))/ 10 , where 10 is the partition function for the [ C ] (1-0) line (Solomon & Vanden Bout 2005), which depends on the excitation temperature. Fig. 7 shows the fraction of the molecular mass in each of our five slices calculated from the [ C ] (1-0) flux and our estimate of the excitation temperature in each slice, based on the assumption of local thermodynamic equilibrium. This demonstrates that although the gas is much hotter in the north, most of the mass of the molecular phase is in the colder gas to the south.

First Estimates of the Mass of the ISM
Any estimate of the mass of the molecular ISM in a galaxy other than our own requires the use of an ISM tracer such as CO molecules or dust. This in turn requires an estimate of the constant of proportionality (the calibration factor) between the luminosity of the tracer and the mass of the ISM. In this section, we estimate the mass of the ISM in ID141 using four different tracers and five different methods. The uncertainties quoted on the ISM mass estimates include the flux uncertainty and the uncertainty of the conversion factors used.

The gas mass from CO
We have two measurements of CO lines, a CO(7-6) line luminosity which has been corrected by the lensing model derived in this paper, and a CO(1-0) line luminosity which has not been corrected for lensing because the observations did not have high enough resolution (Dunne et al. in preparation). Since most of the work on calibrating the CO method has used the CO(1-0) line and since the ratio of high-J to low-J line flux has a large dispersion (Cañameras et al. 2018), we have estimated the gas mass from the CO(1-0) measurement.
The total CO(1-0) line flux is 0.61 ± 0.09 Jy km s −1 (Dunne et al. in preparation). We have corrected this using the magnification factor of 5.5 derived for band 3. Dunne et al. (2021) have derived the CO(1-0) calibration factor from a sample of H-ATLAS galaxies with ∼ 0.35 with ALMA measurements of CO(1-0), [ C ] (1-0) and the continuum dust emission. By combining the measurements, they derive calibration factors for three methods of estimating the mass of molecular hydrogen. All three of these methods scale linearly with the mean gas-to-dust ratio (GDR) in the Universe, which they assume is 135, similar to the estimates of the value in the Milky Way. They derive a calibration factor CO = 3.0 ± 0.5 M (K km s −1 pc 2 ) −1 . Using this value and including a factor of 1.36 to allow for the elements other than hydrogen gives an estimate of the mass of the ISM of (3.3 ± 0.7) × 10 11 M (Table 5). This mass is also consistent with Cox et al. (2011) after correcting for lens magnification and using the value of CO adopted above.

The gas mass from [CI]
A major problem with using the CO molecule as a tracer of the ISM is the mounting evidence that a significant fraction of the molecular phase of the ISM does not contain any CO molecules, probably as the result of photo-dissociation. For example, there is evidence from Fermi (Abdo et al. 2010), Planck (Planck Collaboration et al. 2011) and Herschel (Pineda et al. 2013) that roughly one third of the molecular gas in the Milky Way does not contain any CO. There is also increasing evidence that the fraction of CO-dark molecular gas varies from galaxy to galaxy (Bigiel et al. 2020;Dunne et al. 2021). The [ C ] (1-0) line is an attractive alternative to CO because it does not suffer from the problem of photo-dissociation, is optically thin, and has a similar critical density and excitation temperature to CO Geach & Papadopoulos 2012;Tomassetti et al. 2014;Bourne et al. 2019). We have estimated the gas mass from the line flux of the [ C ] ( 3 P 1 − 3 P 0 ) transition, corrected with our lensing model (Table 4). Using their calibration sample (see above), Dunne et al. (2021) have derived a value for the [ C ] calibration factor of [CI] = 18.8 M (K km s −1 pc 2 ) −1 . The excitation function for the line ( 10 ) varies only weakly with temperature (Dunne et al. 2021) and so we have used this calibration factor rather than trying to make any correction to it using our estimate of the temperature from the ratio of the two [CI] lines (Section 4.3). We estimate that the mass of the ISM in ID141, including all elements, is (3.8 ± 0.8) × 10 11 M (Table 5).

The gas mass from [CII]
There is a strong correlation between the luminosity of the finestructure line of [ C ] at 158 m and the mass of molecular gas in the ISM (Zanella et al. 2018). Zanella et al. derive a [ C ] to 2 conversion ratio of CII 30 M /L . We do not have a highresolution map of the [ C ] line and therefore no specifically-derived lensing model for it. Instead, we simply used the global line flux given by Cox et al. (2011), corrected it for the magnification factor of 5.9 we measured for the dust continuum emission in band 4, the closest in wavelength to the [ C ] line, and then applied the calibration factor of Zanella et al. We obtained a total gas mass of (3.2±0.6) ×10 11 M (Table 5).

The gas mass from dust emission
Many studies have made the case that dust grains are a better way of tracing the ISM than the CO molecule (e.g. Eales et al. 2012;Scoville et al. 2014Scoville et al. , 2016Scoville et al. , 2017. Some of the advantages are that the dust grains do not suffer photo-dissociation like CO molecules, the dust emission is optically thin, and the relationship between the calibration factor and metallicity appears to be simpler (Eales et al. 2012(Eales et al. , 2018. We used the calibration factor between dust luminosity at 850 m and ISM mass derived by (Scoville et al. 2016, their Eq. A11): We used the global flux density measured by Cox et al. (2011) at 3290 m of 1.2 ± 0.1 mJy, which makes it the closest flux measurement in rest-frame wavelength (628 m) to the wavelength used in the calibration factor above, minimising the uncertainty in extrapolating over wavelength. We estimated the luminosity at rest-frame 850 m, using a dust temperature of 25K and an emissivity index of ( = 1.8), the same values used by Scoville et al. (2016). We then calculated the ISM mass using the magnification factor of 5.5 derived from our band 3 lensing model and the calibration factor above, obtaining a total ISM mass of (3.6 ± 0.7) × 10 11 M (Table  5). Dunne et al. (2021) have used several samples of calibration galaxies to derive an optimal technique for estimating the calibration factors ( CO , [ ] , dust ) for an individual galaxy and for producing a best estimate of the ISM mass in the galaxy that uses all three. After correcting for lens magnification, using this technique on ID141 gives an estimate of the ISM mass of (3.6 ± 0.7) × 10 11 M (Table 5). Table 5 lists the ISM mass estimates from the five different methods. We note a few important points. One caveat is that we could have used the lower value of CO that has been claimed for ultra-luminous infrared galaxies. This would lower the ISM mass estimated from CO by a factor of a few, although there are now some compelling arguments that the claim that ultra-luminous infrared galaxies have a different CO than other galaxies is incorrect (Scoville et al. 2016, Dunne et al. in preparation). A second caveat is that the only estimate in which the observations of the tracer have been corrected for gravitational lensing from the observations themselves is the estimate from the [ C ] (1-0) line. For every other estimate, we have corrected the tracer using a total magnification value derived at some other wavelength. There is therefore the possibility of differential magnification (see, for example Serjeant 2012;Omont et al. 2013). However, the similarity between the line profiles measured in the image plane and those in the source plane, together with the similarity of the seven magnification values presented in Table 3 (from four line images and three continuum images), indicate that this is only a minor effect. We have quantified a typical size of the effect by computing the magnification factor of a source that has twice the diameter of the band 3 source shown in Fig. 1, finding this to be only 3 per cent lower than the original magnification. Another important point is that all ISM estimates are only for the molecular phase of the ISM because of the way the methods are calibrated. Finally, da Cunha et al. (2013) and Zhang et al. (2016) have shown that the effect of the cosmic microwave background radiation will make it more difficult to detect the full extent of the ISM in high-redshift galaxies, with this effect depending on the frequency of the observations.

Results
The mass estimates lie between 3.2 × 10 11 M and 3.8 × 10 11 M , which are very large values and suggest we are unlikely to be missing much of the ISM. The agreement between the values is also quite remarkable considering that four of the methods use different tracers of the ISM and we are observing the galaxy only 1.5 billion years after the Big Bang. However, none of the methods are truly independent of each other because they are all calibrated from observations of the Milky Way. We discuss this further in Section 5.2.

The stellar mass
An additional means of establishing the evolutionary stage of a galaxy is by measuring its mass in stars. We return to the theme of stellar mass in Section 5.3 where we consider a chemical evolution model for ID141, but in this section, we outline our estimation of stellar mass of ID141 from photometry following the procedure outlined in Hopwood et al. (2011). The key ingredient is the mid-IR source flux density which we have determined using Spitzer Space Telescope observations (proposal ID. 80156) acquired with the Infra Red Array Camera (IRAC) as described below.
To begin, we measured the surface brightness profiles of the lensing mass from the s -band data presented in Bussmann et al. (2012) using GALFIT (Peng et al. 2002). As Bussmann et al. (2012) found, subtracting the GALFIT model light profiles from the observed sband image showed negligible residual lens light and no evidence of the lensed source image. We then used these profiles to construct model IRAC 3.6 m and 4.5 m surface brightness profiles and subtracted these from the IRAC images to leave residuals resembling the ALMA image (see Fig. A1). We repeated this removal procedure for three other galaxies selected to have similar morphology and flux in the IRAC data and found only negligible residuals. This demonstrates that an accurate IRAC PSF was obtained and that the near-IR traces the mid-IR well in these three systems, although we are unable to rule out the possibility of there being a mis-match in the specific case of the lensing mass in ID141. A final measure of the lensed source IRAC flux was then extracted from an annulus surrounding the observed residual flux (see Fig. A1) resulting in fluxes of 8.7 ± 2.2 Jy and 8.1 ± 1.8 Jy at 3.6 m and 4.5 m respectively. The errors here do not include any uncertainty arising from lens light removal.
In the final part of the process, following Hopwood et al. (2011), we used version 2 of the SED fitting code MAGPHYS (Battisti et al. 2020;da Cunha et al. 2008) and supplied this with the two IRAC fluxes as described above, a limiting s -band flux, the far-IR fluxes as described in Section 4.2 and Wide-field Infrared Survey Explorer fluxes at 12 m and 22 m (WISE; Lang et al. 2016). We refer the reader to Section A for further details.
To de-magnify the quantities output by MAGPHYS, we assumed a magnification of 5.8, equal to the mean of the continuum magnifications given in Table 3. In addition, since MAGPHYS assumes a Chabrier initial mass function (Chabrier 2003), we scaled the stellar mass, dust mass and SFR output by MAGPHYS by a factor of 1.41 for comparison with the values determined in Section 4.2. This factor of 1.41 was determined from the scalings derived in Madau & Dickinson (2014). The resulting stellar mass is M ★ = 4.8 +2.5 −2.3 × 10 11 M . In addition to stellar mass, MAGPHYS also gives a dust mass of 1.8 ± 0.3 × 10 9 M and a SFR of 2100 ± 500 M yr −1 . The 1 errors quoted are taken from the distributions output by MAGPHYS (see Appendix A for more details). These results are consistent with those obtained in Section 4.2 based on our own SED derived from solely the far-IR photometry. Table 4. Emission line fluxes used in the estimation of the gas mass of ID141. The CO(7-6), H 2 O and two [ C ] lines are from the present work and have been determined in the source plane. Errors quoted for these lines include uncertainty due to the lens model and continuum subtraction. The CO(1-0) and [ C ] lines were taken from Dunne et al. in preparation and Cox et al. (2011) respectively and have been corrected for lens magnification using the band 3 and band 4 magnification factors respectively.

The kinematics of the source
The double-peaked nature of our highest SNR emission line profile, the H 2 O line, shown in Fig. 3 suggests that either the source is a rotating disk-like system with a prominent outer ring (see, for example Stewart et al. 2014), or that it is composed of two or more individual merging components.
To better understand the dynamical characteristics of the source, we investigated the spatial extent of the kinematics using the kinematic modelling code 3D Barolo (Di Teodoro & Fraternali 2015). Instead of fitting a kinematical model to the 2D velocity field, this code fits directly to the full datacube. This avoids the effect of artificial flattening of velocity fields due to inadequate spatial resolution and allows for the fitting of additional galaxy properties (i.e., velocity dispersion, disk scale height, brightness profile).

The kinematical model
First, we identified the signal in the data cube by using the 'SEARCH' routine of 3D Barolo, which is based on the DUCHAMP source finder (Whiting 2012). In short, this routine automatically identifies the noise level of the data cube, searches each channel for spaxels that are above a user-provided SNR (here 3.5 ), creates a mask by extrapolating outwards to a second user-provided S/N threshold (here 2.0 ), and combines overlapping masks in different channels into discrete 'sources'. All sources that are too spectrally thin (i.e., possessing only one channel) or spatially compact (i.e., smaller than the synthesized beam of the data cube) are rejected. This final threedimensional mask contains all significant emission, while excluding noise peaks that would be included using a flat S/N threshold.
With the signal identified, we next performed a preliminary analysis of the masked data cube in order to produce physically motivated estimates for further dynamical fitting. The CASA toolkit task image.moments was applied to the masked cubes to create maps of integrated intensity (zeroth-moment), line-of-sight velocity (first-moment), and velocity dispersion (second-moment). The zeroth-moment map was fitted with a two-dimensional Gaussian using the CASA toolkit task image.fitcomponents, resulting in a central position, morphological position angle, and FWHM of the minor and major axes. The initial estimate of the rotational velocity at this central position was then extracted from the first-moment map.
Before performing the fitting, it is crucial to determine the size of each ring, as well as the maximum radius of the model. In order to avoid over-resolving the data, we adopted a minimum ring width of a factor of 0.4 times the FWHM of the minor axis of the restoring beam. The factor of 0.4 here is comparable to factors used in other similar works (e.g., Talia et al. 2018;Fan et al. 2019). Similarly, we adopted a maximum model radius of a factor of 0.8 times the FWHM of the major axis of the 2D Gaussian fitted to the zeroth-moment map. The factor of 0.8 is the largest value that does not return non-physical results at large radii.
With these initial parameters and estimates in hand, we engaged the 3DFIT routine of 3D Barolo to fit the masked data cube with a tilted ring model comprising a series of concentric rings which each rotate with an independent velocity. The number of rings and ring width were determined as detailed in the previous paragraph, while the model center was fixed to the central position based on the 2D Gaussian fitted to the zeroth-moment map. We assumed a thin disk (i.e., = 0.01 arcsec). This routine begins by populating a model physical space (i.e., X, Y, Z) with gas clumps, so that these initial estimates of galaxy properties are replicated. This physical space is then translated into an observed-space data cube (i.e., RA, Dec, velocity), and the variable galaxy properties (i.e., rotational velocity, velocity dispersion, inclination, position angle, and systemic velocity) are perturbed until a minimum residual is met. This process is repeated for each ring, after which most variables are averaged, and the fitting process is repeated with a smaller number of variables (i.e., rotational velocity and velocity dispersion). This process results in best-fit values of each parameter and a model data cube that may be compared to the original data cube on a pixel-by-pixel basis.
The normalisation of the flux density within each ring can be either 'local', meaning that the integral along the spectral axis of a given model pixel in the ring must match that of the corresponding datacube pixel or, 'azimuthal', whereby the model is normalised to the azimuthally-averaged flux in each ring. Local normalisation allows for non-axially symmetric distributions of emission across the disk that can bias fitting of azimuthally-averaged fits. In this work, we have used both normalisation types for comparison.
We applied 3D Barolo to the two emission lines with the highest SNR, namely the CO(7-6) and H 2 O lines.

The results of the kinematic model
By applying this tilted ring fitting code to the CO(7-6) and H 2 O data cubes, we found good fits (as seen in the comparisons of Fig.  8). This implies that ID141 is a strongly rotation-dominated disk, a fact that is also apparent from the rotation curves shown in Fig.  9. The resulting dynamical masses and galaxy properties are listed in Table 6 and show close agreement with each other, regardless of the line or normalization employed (local or azimuthal). Comparing these dynamical mass estimates with our ISM mass estimates in the preceding section, we find a very significant discrepancy to the extent that the ISM masses are a factor of approximately four times those of the dynamical mass estimates. We discuss possible remedies for this disparity in Section 5.
One of the strengths of three-dimensional dynamical characterization is the ability to compare the data and best-fit model on a pixel-by-pixel basis. Inspecting Fig. 8, it is clear that both CO and H 2 O are poorly fit by models using azimuthal averaging (note the strong residuals in zeroth-moment maps and spectra). However, the velocity fields, velocity dispersion maps, and general shape of each position-velocity diagram are well-reproduced. This raises an interesting issue: the kinematics of the galaxy are well-fit by a rotating disk, but the distribution of emission (i.e., morphology) is not. This feature will be discussed in Section 5.1 where we consider the possibility of ID141 being a merger of two or more galaxies.

One galaxy or several?
While the strong emission lines in this paper (i.e., CO and H 2 O) are well-fit by 3D Barolo (see Section 4.6.2) and return similar dynamical masses, the source morphologies are complex, resulting in imperfect fits. One key piece of information in interpreting this issue is the lack of symmetry between the velocity field structure to the north and south. While the velocity field to the north is a continuation of the central gradient, it extends to the east, and lacks a strong southern counterpart. When combined with the strong line emission and the H 2 O dispersion peak in the northern region, this raises the possibility that the northern emission is not an extension of ID141, but a separate galaxy.
To investigate this possibility, we masked the northern section of the CO and H 2 O data cubes, and fitted this "main" galaxy using 3D Barolo. In each case, we assumed azimuthal normalization and the same model central position as in the unmasked data cubes. The fit to these masked data significantly improve upon the azimuthal fit to the unmasked data with a disk-like model, as shown in Fig 10. Next, we took the best-fit models of these masked CO and H 2 O data cubes and subtracted them from the unmasked data cube (CASA task immath) to isolate the signal to the north (hereafter referred to as the 'secondary'). As seen in the lower row of Fig. 10, the secondary signal shows a continuous velocity gradient. The corresponding dynamical mass of this secondary is (3.9 ± 0.9) × 10 10 M and when added to the dynamical mass of the main galaxy, increases the total dynamical mass by approximately 50 per cent compared to the single-galaxy model. The secondary exhibits an offset velocity dispersion peak and two separate components in each zeroth-moment map. Because of these properties, 3D Barolo returns a poor fit using azimuthal normalization, and a better but still unsatisfactory fit using local normalization (which we also show in Fig. 10). While 3D Barolo is able to reproduce the brightness distribution and overall velocity gradient of each data cube, the significant residuals in each moment map and position-velocity diagram suggest that this secondary source is not a single, well-ordered rotating galaxy.
Taking these findings into account, there are three obvious possibilities for the nature of ID141: (i) The system may simply be a single rotating disk. In this case, the asymmetry to the north would be a minor disturbance.
(ii) The source to the north may be a separate, rotating galaxy that is closely associated to ID141. This is supported by the continuous velocity gradients shown in the unmasked data cube (see Fig. 8) and reasonable fits to the H 2 O emission from the 'main' and 'secondary' sources (see right column of Fig. 10).
(iii) Because the 'secondary' source shows two kinematicallyseparated clumps in both CO and H 2 O emission, it is conceivable that the ID141 system is composed of three significant clumps: the main galaxy, a source to the north that is strong in CO emission, and a third source to the northeast that is strong in H 2 O emission.
Due to the diversity in observed continuum (Fig, 1) and line (Fig.  2) morphologies, we are unable to ascertain which of these scenarios, if any 1 , is the most likely without further observations although most bright SMGs are known to comprise mergers (e.g. Engel et al. 2010).

Why do the ISM mass estimates agree so well and why are they greater than the total mass?
The five gas masses derived in Section 4.4 (see Table 5) have a remarkably small range of (3.2 − 3.8) × 10 11 M , especially since four of the methods use different tracers. They are also remarkably larger than the total mass estimated from the kinematics of the galaxy; both application of 3D Barolo and a simple estimate from 2 / with velocity = 430 km s −1 and radius = 2 kpc give a dynamical mass of ∼ 0.8 × 10 11 M . We first consider possible observational explanations of the mass discrepancy and then suggest an astrophysical explanation of both points. One possible explanation for the mass discrepancy is that the lens model and therefore the lens magnification is incorrect. In this case, the derived intrinsic source properties will be systematically biased. We can estimate the extent to which such an error will propagate into the ratio of gas mass to dynamical mass. Regarding the gas mass, since ∝ int ∝ obs / , where int and obs are the intrinsic and observed (i.e., lensed) line luminosities, the gas mass scales inversely proportionally with magnification. Under the assumption of circular orbits used by 3D Barolo, the dynamical mass scales as dyn ∝ 2 int , where is the circular velocity and int is some intrinsic physical length scale. Since the observed is invariant under lensing and the intrinsic physical length scales as int ∝ obs / √ , the dynamical mass scales with magnification as dyn ∝ −1/2 and so / dyn ∝ −1/2 . Therefore, the magnification we have obtained from our lens modelling would have to be underestimated by a factor of 16 to obtain a gas mass that is equal to the dynamical mass.
Although the statistical uncertainties we quote are orders of magnitude beneath this, larger errors can arise if the lens model we have adopted is not a faithful representation of the true lensing mass. As previously stated, changing our lens model by adding an additional mass profile to match a possible second lensing galaxy only modifies the magnification by approximately 5 per cent. An alternative is to use a mass profile that does not follow a power law. Schneider & Figure 8. Comparison of observed emission lines and best-fit 3D Barolo models for local normalization (top row) and azimuthal normalization (bottom row), as well as both CO (left column) and H 2 O (right column). In each set of panels, we show the zeroth-moment (i.e., integrated intensity) maps (top row), first-moment (i.e., line-of-sight velocity) maps (second row; solid and dashed line indicates major and minor kinematical axis respectively), and second-moment (i.e., velocity dispersion) maps (third row) for the data (left column), model (middle column), and the residual between the two (right row). We also show the position-velocity diagrams taken along the major axis (bottom row, left column) and minor axis (bottom row, middle column) of the data cube (colour) and model cube (contours). Finally, we compare the spectra of the data (blue), model (red), and residual (green) in the lower right panel, showing the 1 noise level by the shaded grey region.  Sluse (2013) present an example of a composite lens model comprising two significantly different density profiles with a resulting magnification that differs from the power-law model by 40 per cent but still broadly reproduces the observed lensed images. Even this rather extreme example is over an order of magnitude beneath the size of shift required to reconcile our dynamical and gas mass estimates.
Another possible explanation for the discrepancy is that the total (dynamical) mass has been estimated from the CO (7-6) and water lines while the ISM estimates have been obtained from observations of the continuum and different spectral lines. The [ C ] (1-0) line, which is probably the best tracer of the cool ISM but which is too weak for us to perform a full kinematic analysis, does have a different spatial distribution from the lines used to investigate the kinematics (Fig. 2). However, the velocity extent of all three lines is very similar (Fig. 3) and since the kinematic analysis extends out to a similar radius to the [ C ] (1-0) emission (Fig. 3), it is unlikely that our estimate of the total mass is too low. A related issue that has been argued in the literature (e.g., da Cunha et al. 2013; Zhang et al. 2016) is that the mass estimates of high-redshift galaxies may be too low because the increasing energy density in the cosmic microwave background with redshift makes it increasingly hard to detect extended gas. This is unlikely to remedy such a large discrepancy because, to first order, this effect simultaneously reduces our measurements of the total mass (by causing an under-estimation of dynamical radius) and the ISM mass.
An additional concern to consider is in the kinematic analysis from Section 5.1. We were ultimately unable to confirm that ID141 is a single rotating disk rather than several galaxies. The stronglylensed SMG SMM J02399-0136, initially identified as a large disk galaxy (Genzel et al. 2003) and subsequently resolved into a complex system of several components (Ivison et al. 2010) is an example of why caution must be taken. Nevertheless, in the case of ID141, this is unlikely to explain the large discrepancy. If we consider option (ii) in Section 5.1, that ID141 consists of two rotating disks, this gives the largest of our dynamical mass estimates of 1.3 × 10 11 M which is still significantly short of explaining the discrepancy. Figure 10. Comparison of observed emission lines and best-fit 3D Barolo models for masked CO (left column) and H 2 O (right column) data cubes. In each set of panels, we present the zeroth-moment (i.e., integrated intensity) maps (top row), first-moment (i.e., line-of-sight velocity) maps (second row), and second-moment (i.e., velocity dispersion) maps (third row) for the data (left column), model (middle column), and the residual between the two (right row). We also present the position-velocity diagrams taken along the major axis (bottom row, left column) and minor axis (bottom row, middle column) of the data cube (colour) and model cube (contours). Finally, we compare the spectra of the data (blue), model (red), and residual (green) in the lower right panel, showing the 1 noise level by the shaded grey region.
Finally, we note that the discrepancy becomes even worse if we include stellar mass. Our estimate is that the stellar mass in this system is M ★ = 4.8 +2.5 −2.3 × 10 11 M , which in itself is greater than the dynamical mass. However, the uncertainty on this measurement is high and furthermore, does not include any uncertainty arising from removal of the lens light.
We now consider a possible astrophysical solution, which might also partly explain the surprising agreement between the mass estimates made with four different tracers, a solution which arises because none of the methods is really independent. The dust method of Scoville et al. (2016) is calibrated using CO(1-0) measurements of galaxies and is based on the CO(1-0) calibration factor derived for giant molecular clouds in the Milky Way; the calibration factor for [ C ] (Zanella et al. 2018) is based on dust measurements, which are again ultimately based on the CO(1-0) calibration factor from galactic giant molecular clouds; and all of the calibration factors discussed in Dunne et al. (2021) are based on the assumption that the GDR is the same as in the Milky Way. All the methods are ultimately benchmarked to the same astrophysical source, for the obvious reason that it is only in the Milky Way (and to a much lesser extent in galaxies in the Local Group) where it is possible to make the direct measurements of the ISM masses that are needed to calibrate the tracer methods.
A possible solution is if the calibration factors for ID141 are very different to the Galactic ones. To explain the agreement observed between the different mass estimates, all of the gas-to-tracer ratios would have to be 4 times lower than in the Milky Way. Since all of the tracers are composed of heavy elements, the most direct way to bring some insight to the problem is with a chemical evolution model.

Chemical Evolution in the Epoch of Galaxy Formation
In this section, we describe the results of a chemical evolution model we have constructed to investigate the likely evolution of the interstellar dust, gas and stellar mass of ID141. The details of this model can be found in Appendix B. We have used this model to investigate the evolution of the GDR. Although we have not modelled the evolution of the other tracers (gas-to-CO etc.), it is likely that the results will be very similar since all of the tracers must depend on the abundances of the metals in the galaxy. The similarity of the ISM estimates for four different tracers provides strong circumstantial evidence in support of this.
In the model, which is an extension of the model in De Vis et al. (2017), the galaxy starts as a cloud of gas with no heavy elements that is gradually converted into stars. The model incorporates gas flowing in to the galaxy via accretion from the cosmic web and out from it via stellar feedback and active galactic nuclei. Dust is produced by stars, both asymptotic giant branch stars and supernovae, and grain growth and it is destroyed by astration, supernova shocks and outflows.
We considered two star formation histories, one with continuous star formation starting at some formation redshift, f , and a second where we are currently observing 40 per cent of the way through a burst lasting 200 Myr 2 . Hereafter, we refer to these as the 'conSF' and the 'SB' model.
To account for possible variations in the IMF (Appendix B), we used three functions that span most of the range of IMFs considered in the literature, namely a Salpeter IMF (Salpeter 1955), the Chabrier IMF (Chabrier 2003) and the more top-heavy IMF proposed by Zhang et al. (2018), required to explain isotope ratios at high redshifts 3 . In Section 4.2 we estimated a SFR for ID141 from its far-infrared luminosity of 2400 ± 500 M yr −1 , based on the assumption of an IMF from Kroupa & Weidner (2003). If we instead use the IMFs considered here, we obtain an SFR of 2800, 1800 and 1300 M yr −1 for the Salpeter, Chabrier and top-heavy IMF. In our model we set f = 4.8, 5.5 and 7 respectively (the higher observed current-day SFRs for the Salpeter and Chabrier IMFs leads to a more rapid conversion of gas into stars, and so their star formation needs to begin later). The parameters of each model have been chosen so that at the redshift of ID141 ( = 4.24) the SFR predicted by the model is similar to the SFR estimated for ID141 from its far-infrared luminosity.
For the conSF star formation model, we assume a SFR that follows the Kennicutt-Schmidt law (Kennicutt 1998), i.e., ( ) = ( ) where is the star formation efficiency, is the total gas mass and is set to unity. is estimated using = / dyn where dyn is the dynamical (free-fall) timescale dyn = √︁ 3 /2 (Walter et al. 2009;Swinbank et al. 2010) and is the observed SFR (see also Appendix B). Given an effective galaxy radius of 1.4 kpc (measured from the reconstructed source continuum) and an observed ISM mass of 3.5 × 10 11 M , the maximum SFE is 1 per cent for the Salpeter IMF. If we instead assume that the true ISM mass must be less than this (because of the limit from the dynamical mass) and set the ISM mass equal to the dynamical mass of the one-disk model ( dyn 0.8 × 10 11 M ), the range of maximum SFE obtained is 5 and 2 per cent for the Salpeter and top-heavy IMFs respectively. Although ID141 is one of the most extreme objects in the universe in terms of its SFR, both estimates of SFE suggest that it is not forming stars at a significantly more efficient rate than normal systems. The SFE is similar to other high redshift SMGs Riechers et al. 2013) but an order of magnitude lower than the = 6.4 hyper-starburst from Walter et al. (2009). Using the current observed level of SFR and assuming that ≈ dyn when the galaxy started forming stars, the gas reservoir would be exhausted within 40 Myr for a Salpeter IMF if there is no additional replenishment of gas via mergers or inflows.
In Table 7, we compare the model outputs predicted at the redshift of ID141 of = 4.24. The Salpeter IMF models produce lower dust masses (since there are fewer high mass stars producing metals and dust and there is more gas being locked up in stellar remnants) and higher stellar masses than the Chabrier and top heavy IMF models. The resulting gas fractions are also lower, indicative of a much later evolutionary stage. Davé et al. (2011) predict gas fractions of ∼ 0.4 at = 4 for a current day halo mass of ℎ ( = 0) = 10 14 . This is larger than the Salpeter IMF results but more in line with the results from the starburst, top heavy IMF model. The total baryonic mass ( + * ) predicted by the conSF models is ∼ 5 times larger than the observed dynamical mass of ID141, though the starburst models predict a more comparable mass (∼ 1 − 3×).  Table 7. The values predicted by our chemical evolution models for the properties of our ID141-like galaxy at the redshift at which we observe ID141: = 4.24. We give results for two different star-formation histories -starburst ('SB') and continuous ('conSF') -and for the Salpeter, Chabrier and top heavy initial mass functions. is the gas fraction. GDR denotes the gas-to-dust ratio and M tot is the total (gas mass + stellar mass) returned by the model. The metallicity, , is defined as the ratio of the metal mass and gas mass. By design, the SFR predicted by each model is similar to the SFR estimated for ID141 from its far-infrared luminosity.

Evolution of the Gas to Dust Ratio
In Fig. 11 we show the evolution of the model GDR as a ratio of the Milky Way GDR, ( / ) model /( / ) MW , for each model plotted against = /( + ★ ), the fraction of baryonic mass that is gas. For clarity, we only show and discuss the results for the Salpeter and top heavy IMFs as the results assuming a Chabrier IMF sits mid-way between the two. The shaded area around each line shows the effect of changing the parameter governing the rate of dust grain growth in the dense clouds (Appendix B) from 700 to 15,000, where the former value is appropriate for the Milky Way today and the latter indicates extremely rapid grain growth (Mattsson & Andersen 2012;Asano et al. 2013;De Vis et al. 2017). The GDR is high at high gas fractions because the abundance of metals (and therefore dust) is still low. The GDR decreases as the gas is converted into stars and as dust is formed in stellar winds, supernovae and, after reaching some critical metallicity, via grain growth.
We have assumed a value for the GDR ratio of the Milky Way that is appropriate for the Scoville et al. (2016) calibration used in Section 4.4.4 to convert dust mass and submm continuum luminosity to gas mass. We have estimated this by assuming = 1.8, d = 25 K, 850 = 0.077 m 2 kg −1 (James et al. 2002), and ISM = 4.8 × 10 −4 m 2 kg −1 (Scoville et al. 2014). Since the models account for the total gas, whereas the calibration of dust-to-gas ( 850 ) provided by Scoville et al. is for the molecular component only, we follow Scoville et al. (2014) and assume that there are roughly equal amounts of H and H 2 . We also account for Helium. This gives a GDR for the Milky Way of 160 given the assumed 850 (a factor of ∼1.2 higher than that assumed in Dunne et al. 2021). Table 7 shows that, for all models, the GDR reached by redshift 4.24 is < 160 (i.e., the ratio of the model GDR to the Milky Way GDR is less than unity), suggesting it is not appropriate to assume a Milky Way calibration when estimating gas masses from the dust emission in SMGs like ID141. Such low GDRs have been measured in other high redshift SMGs (see, for example, Berta et al. 2021;Yang et al. 2017, and references therein). Fig. 11 also shows that the GDR can vary by over an order of magnitude during the SMG phase of a galaxy. Therefore it is also dangerous to assume a constant value for the GDR when estimating the ISM mass from the submm continuum emission. However, we can use the model to our advantage: the results in Fig. 11 provide a way to use the model itself to estimate the true GDR, and thus the mass of the ISM. We can also use it to probe how far through ID141 is in its conversion of gas into stars i.e., its evolutionary phase. Defining (obs) to be the ISM mass we have estimated on the assumption of the Galactic GDR, the true ISM mass, including the contribution of the atomic gas, is given by: The total mass of the model galaxy is then given by: We have made the assumption that there is no contribution from dark matter (see Genzel et al. 2017), but if there is a contribution from dark matter it will strengthen the conclusions below. We show the evolution of true,tot in Fig. 11 as dashed lines. The shaded horizontal region indicates where the total baryonic mass (gas plus stars) is consistent with the dynamical mass ( + ★ ≤ dyn ) within the uncertainty, defined as 3 .
We can immediately draw the following conclusions: (i) it is only possible to make the masses consistent if the IMF is top-heavy, adding to the evidence from isotopic abundances (Zhang et al. 2018) that the IMF is different in high-redshift star-forming galaxies.
(ii) It is possible to make the masses consistent when the galaxy is roughly halfway through its evolution ( ∼ 0.5). Although it is often assumed that galaxies like ID141 must have very high gas fractions, the discrepancy between the mass estimates gets rapidly worse above ∼ 0.6.
We also notice that the true total mass calculated using Eq. 4 overestimates dyn in comparison to the model tot when assuming a Salpeter IMF. This is due to a combination of the faster consumption of gas due to the higher SFR and the lower dust masses formed per stellar population with a Salpeter IMF. For example, to reach a total mass 1 × 10 11 M using Eq. 4, the conSF Salpeter model would require a GDR/ g 50, or alternatively, for a galaxy with g < 0.5, a GDR of 25 is needed. Such low GDRs are not reached in the Salpeter IMF models.
Although the conSF top heavy IMF appears to result in Eqn. 4 producing a total mass within the shaded region, this is misleading. The variation of the GDR ratio with from this model is low enough that scaling the observed gas mass produces a true total mass from Eqn. 4 that lies below the dynamical mass limit (in the shaded region). However, with the continuous SF model, such low GDRs can only be Figure 11. The gas-to-dust ratio (solid lines), as a fraction of the Milky Way value, predicted by standard dust-evolution models plotted against , the fraction of baryon mass in the form of gas for top: the continuous star forming (conSF) model and bottom: the starburst (SB) model. The pink and blue colours denote the top-heavy and Salpeter IMFs respectively and the circles (filled and solid) show the values at = 4.24. The dashed lines show the 'true' total mass ( true,tot ) derived from the GDRs from the models and the observed gas mass assuming a MW-like GDR (Eqn. 4) as the galaxy evolves. The black horizontal line indicates a GDR equal to the Milky Way. The shaded grey region indicates + ★ dyn , where we use the maximal value ( dyn = 1.3 × 10 11 M ) with a 3 uncertainty. produced in a model that has already consumed a lot of its gas and built up a large reservoir of stars, hence violating the mass criteria we are trying to satisfy ( tot (conSF, top heavy) > 4 × 10 11 M in Table 7).
Finally, we note that if we increase the amount of gas accreted in the model, although this would increase the gas fraction, it would also increase both the GDR and the stellar mass, hence would not solve this issue.

Metallicity Evolution
The low GDRs we predict for ID141 imply these galaxies will reach supersolar metallicities by = 4.24 (where = metals / = 1.5 − 4.4 4 , Table 7) depending on the assumed IMF and SFH.
To check if such high metallicities are expected in these galaxies, we next compare the evolution of the gas-to-dust ratios ( Fig. 11) with metallicity in Fig. 12. We convert the model oxygen abundance assuming a top-heavy IMF is 12 + log(O/H) = 9.33 and 9.16 for the constant SF and starburst respectively. Local GDR-scaling relations from Herschel surveys, e.g. Rémy-Ruyer et al. (2014); De Vis et al. (2019), are shown in Fig. 12. These relationships imply values of 12+log(O/H) that are consistent with the predicted model GDRs of 25-50 for ID141. At the same metallicity as the SB model with a top heavy IMF, the Rémy-Ruyer et al. (2014) relationships predict GDRs ranging from 28-55. The properties of local galaxies from the DustPedia sample of De Vis et al. (2019) also reach low gas-to-dust ratios, but at lower metallicities than those predicted for ID141 (Fig. 12). This sample is made up mostly of late type galaxies, with stellar masses > 10 8 making up only 4 per cent.
We also show the relationship from Santini et al. (2014) derived by combining observed dust masses of galaxies out to = 2.5 with metallicities derived from the stellar mass-SFR-metallicity relationship of Mannucci et al. (2010). Although a number of massmetallicity relations exist in the literature (eg Bothwell et al. 2013;Peng et al. 2015, Maiolino & Mannucci 2019 and references therein) we caution here that there are very few studies carried out at high enough stellar mass, at redshifts > 3, and with SFRs > 100 M yr −1 , making it difficult to apply these relationships to ID141.
What we can say, is that if we assume that ID141 will shortly turn into a massive quiescent early-type galaxy with very little future evolution, then our models predict high mass galaxies with super solar metallicities ∼ 1.3 − 3 at present day. This is in agreement with observations of massive ellipticals in the local universe (Gallazzi et al. 2005;Peng et al. 2015, see also Maiolino & Mannucci 2019, their Fig. 12).

Broader Implications of Low Gas-to-Dust Ratios
The conflict between dynamical mass and ISM mass when estimating the latter using CO has been previously reported in the literature (e.g., Tacconi et al. 2008;Hodge et al. 2012;Bothwell et al. 2013;Calistro Rivera et al. 2018). A common conclusion is to reduce the CO conversion factor, CO 1 (but, in contrast, see Boogaard et al. 2021). However, even with this low value, the tension between the dynamical mass and the ISM mass is sometimes still seen, requiring an additional explanation (e.g. in changing the structure of the galaxy Bothwell et al. 2013). Low values of CO seemed plausible given the evidence for a low value in low-redshift Ultraluminous Infrared Galalaxies (ULIRGs) (Downes & Solomon 1998). More recently, this evidence has weakened (Papadopoulos et al. 2012c,b, also see Dunne et al. in preparation), which makes the low CO explanation of the mass tension given in these previous studies less likely. The advantage of this work is that we have used four different tracers of the ISM: [ C ] , [ C ] , dust and CO. Therefore, although we might resolve the tension in our ISM estimate from the CO tracer by reducing CO by a factor of ∼ 4, we would also have to reduce the three other calibration factors by the same amount.
Given the quite startling agreement between the ISM mass estimates for ID141, which has been noticed for other galaxies (e.g., Dye et al. 2015), it seems more natural to assume that the true explanation of the mass discrepancy is due to chemical evolution rather than the assumption that the physical conditions in the ISM are producing a low value of CO . The former naturally explains why the calibration factors for all four gas tracers are lower in ID141.
This may seem counter-intuitive at first since galaxies in the first few billion years might be expected to have higher GDRs than in lowredshift galaxies because there has been less time for heavy elements to form (e.g., see Péroux & Howk 2020). There is one strong argument why, in this situation, common sense is probably wrong. All the highredshift galaxies observed have been found because they contain a lot of dust or because they already contain a large mass of stars (Scoville et al. 2014(Scoville et al. , 2017Tacconi et al. 2018;Millard et al. 2020), and so there is a strong selection effect that will mean that samples of high-redshift galaxies are likely to contain galaxies in which the stellar mass or dust mass is already high. The simplest closed-box models for the formation of dust (Eales & Edmunds 1996, Dunne et al. 2003, for example, predict that the dust mass in a galaxy is at a peak when ∼ 0.5, and thus the galaxies in a sample compiled from a submm continuum survey are likely to have high dust masses and to have already turned much of their gas into stars. If our conclusion that the gas-to-tracer ratio is 4 times lower than in the Milky Way is true for high redshift galaxies in general, it would have a dramatic effect on our assessment of their evolutionary states, reducing all ISM mass estimates by a factor of 4. However, the variation in the GDR seen in our chemical-evolution model (Fig. 11) shows that calibrating all the tracers using ID141 would be just as risky as using the Galactic calibrations. The indication from Fig. 11 is that there is no perfect calibration, and it is better to think of a gas-totracer ratio as being a new useful tool for assessing the evolutionary state of a galaxy.

How quickly will ID141 become a quiescent galaxy?
In this section, we consider the evolution of ID141 itself. In particular, we examine whether the evolution might be fast enough for it to join the population of quiescent galaxies that is seen at ∼ 3 (Glazebrook et al. 2017;Girelli et al. 2019;Valentino et al. 2020;Forrest et al. 2020b,a;Merlin et al. 2019). For reference, the time interval between the two redshifts is 7.2 × 10 8 yr.
To begin addressing this question, we turn to our kinematic analysis (Section 4.6) which shows evidence for a large rotating disk, with possibly one or two additional merging systems. Since high-redshift galaxies are often separated into two classes -disks or mergers -the kinematics of ID141 suggest that both are often going on at the same time. The escape velocity of the system is given by: Given a dynamical mass of 1.3 × 10 11 M and a radius of 1.4 kpc, the escape velocity is ∼ 910 km s −1 . This is similar to the velocity range of the system (Fig. 3), which suggests that most of the material we see is destined to end up in a single galaxy. There are two other observational results that are relevant to the dynamical evolution of the system. First, the [ C ] (1-0) image (Fig. 1), which likely gives us the most unbiased view of the gas, shows that the gas is divided into a number of clumps. The existence of clumps in cleaned interferometric data needs to be treated with care (e.g., see Ivison et al. 2020) since there are observational effects that can lead to them being falsely identified. However, when clumps are multiplyimaged in a strong lensing system such as ID141, these observational effects are mitigated. A particularly striking example of this is in the system SDP81 which unequivocally shows multiply-imaged clumps in high resolution, high signal-to-noise ALMA data (Dye et al. 2015).
Clumps are expected to migrate into the centre of the galaxy through dynamical friction on a timescale given by (Genzel et al. 2011 where 0 is the velocity dispersion and is the circular velocity. The ratio of and 0 is 10, giving a migration timescale of inspiral 100 dyn 2.3 × 10 8 yr. Alternatively, we can consider the dynamical evolution of the system by calculating the Toomre Q-parameter (Toomre 1964), which quantifies the stability of rotating, thin symmetric disks against gravitational collapse (Toomre 1964). For a flat rotation curve, the Qparameter is given by (Genzel et al. 2011): Using the results of our chemical evolution model, we assume gas / tot 0.5. With this value of the mass ratio, the results from the kinematical modelling of the CO line ( Fig. 9 -left panels) imply that 0.6 in the central few hundred parsecs of the galaxy, falling to a lower value beyond a radius of 1 kpc. The results from the kinematical modelling of the H 2 O line ( Fig. 9 -right panels) imply much lower values of gas . Our results therefore show that the gas disk is unstable. We reached the same conclusion for the disk of the other gravitationally-lensed Herschel source that has been observed at high resolution with ALMA, SDP81.
All of these results show that the dynamic evolution of ID141 is very fast and suggest that it could easily turn into a single fairly undisturbed galaxy by ∼ 3.
Turning now to the evolution of ID141's gas, the only chemical evolution models that produce galaxies consistent with the tight mass limits are those with a top-heavy IMF and a starburst ( Table 7). The conSF model has a gas mass and star-formation rate at = 4.24 of 1.16 × 10 11 M and 1336 M yr −1 respectively, but the model predicts a stellar mass greater than the dynamical mass limit. However, if we assume that the galaxy continues to form stars at the same rate, and that the galaxy accretes no fresh gas, the gas will be consumed in 9 × 10 7 yr. If we perform the same calculation with the SB model, we obtain a timescale of 6 × 10 7 yr.
In conclusion, both the dynamics and the other properties of ID141 are consistent with the idea that, although it is an extreme system at = 4.24, by a redshift of ∼ 3 it will have evolved into a single, fairly undisturbed galaxy, containing little gas.

CONCLUSIONS
ID141 is a gravitationally-lensed hyper-luminous SMG at a redshift of 4.24 with a star-formation rate of 2400 M yr −1 , making it an extreme example of an already extreme class. The system was chosen as a target for a comprehensive high-resolution investigation of the phases of the ISM in high-redshift galaxies because at its redshift, the key ISM spectral lines fall in the ALMA bands. In this paper we have presented high-resolution observations in the [ C ] (1-0) line, the [ C ] (2-1) line, the CO(7-6) line and the H 2 O(2 1,1 − 2 0,2 ) lines. After the lensing reconstruction, the images in the source plane reach down to a resolution of 350 pc. In this section, we outline the observational results and our inferences from these results.
Our observational results are: (i) The gas is distributed over a region of 3 kpc. The [ C ] (1-0) emission is offset by 1 kpc from the other lines, including the [ C ] (2-1) line. This shows that the galaxy contains a large reservoir of cool gas that is missed by the other spectral lines that we have observed with ALMA, although the estimates of the ISM mass from the [ C ] (1-0) and CO(1-0) lines are remarkably similar (see below).
(ii) We have used the ratio of the [ C ] (2-1) and [ C ] (1-0) lines to measure the variation in excitation temperature over the galaxy. We find that the temperature increases from ex ∼ 40K where the [ C ] (1-0) emission peaks to ex 100K where the CO and water lines are at a peak.
(iii) The gas disk appears to be broken into clumps.
(iv) We have used the CO and water lines to carry out a kinematic analysis of the galaxy. We find strong evidence of a rotating system but cannot distinguish between there being a single rotating galaxy or there also being one or two additional galaxies. Our estimate of the dynamical mass depends on which of these scenarios we choose, but our estimates lie in the range 0.8 − 1.3 × 10 11 M .
(vi) After subtracting the lens from a Spitzer image, we estimate that the stellar mass of ID141 is M ★ = 4.8 +2.5 −2.3 × 10 11 M .
From these basic observational results we make the following inferences: (i) All five of our ISM estimates are very similar and are substantially larger than the dynamical mass of the system (without even considering the stellar mass and any dark matter). Inconsistencies like this have been observed before and have always been attributed to a low value of CO , which is then attributed to physical conditions in the ISM of an SMG probably being different from the Milky Way. This can not be the explanation in ID141 because we are using four different tracers, each of which will be affected by different physical conditions in the ISM in different ways. The similarity of the ISM mass estimates using the four different tracers suggests that chemical evolution could explain the mass discrepancy because the abundances of all the tracers will depend on the abundances of heavy elements in the galaxy.
(ii) We have used a chemical evolution model to show that the GDR does change substantially during the evolution of an SMG. It can reach values much lower than the Milky Way value as long as the IMF is top-heavy in SMGs for which there is now good evidence from spectroscopy of molecules containing rare isotopes (Zhang et al. 2018).
(iii) Although there is a complex web of calibrations, all estimates of the ISM masses in high-redshift galaxies are ultimately calibrated in the Milky Way. Our conclusion that the gas-to-tracer ratio for four different tracers is four times lower than in the Milky Way puts all these estimates into doubt. On a positive note, we suggest that the gas-to-tracer ratio may be a useful new tool for assessing the evolutionary state of a galaxy.
(iv) We have used ID141 to test whether extreme SMGs like this might evolve into the population of quiescent galaxies that is seen at ∼ 3. Several indicators (the Toomre Q-parameter, the inspiral time for the clumps, and the escape velocity system of the system) suggest that ID141 should be a single, fairly undisturbed galaxy by ∼ 3. The depletion time for the gas also implies that, on the assumption of no gas infall, the galaxy will be devoid of gas by this redshift. Therefore, extreme though it is at = 4.2, it seems plausible that ID141 will evolve into an unexceptional galaxy by ∼ 3. I00/DOI:10.13039/501100011033).This paper makes use of the following ALMA data: ADS/JAO.ALMA #2013.1.00358.S, ADS/JAO.ALMA #2016.1.00450.S and ADS/JAO.ALMA #2017.1.00029.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. In addition, publications from NA authors must include the standard NRAO acknowledgement: The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

DATA AVAILABILITY
All data presented in this paper can be obtained from the ALMA science archive under the project IDs 2013.1.00358.S, 2016.1.00450.S and 2017.1.00029.S.

APPENDIX A: STELLAR MASS ESTIMATION
Our estimate of the stellar mass of ID141 follows the procedure used by Hopwood et al. (2011) who used MAGPHYS to fit an SED to observed optical, mid-IR and far-IR photometry. In this paper, we have used version 2 of MAGPHYS (see Battisti et al. 2020).
The dominant constraints on the SED in the mid-IR region are provided by Spitzer Space Telescope data which were obtained under proposal ID 80156 using the Infra Red Array Camera (IRAC) operating at 3.6 m and 4.5 m. To obtain a measure of the source emission, it was first necessary to remove the lens emission since the two are blended due to the relatively large IRAC PSF. This was done by first determining the surface brightness profiles of the two foreground galaxies using the s -band data presented in Bussmann et al. (2012) using GALFIT (Peng et al. 2002). Subtracting the model profiles of these two well-resolved foreground galaxies which lie well within the lensed source light observed by ALMA left only negligible residuals and no evidence of any source emission. We then used these profiles to construct model surface brightness profiles of the lens galaxies as they would be observed by IRAC at 3.6 m and 4.5 m. To achieve this, we extracted the PSF at both wavelengths using five stars in the field of view of the IRAC observations and convolved these with the s -band profiles. Subtracting these model light profiles from the observed IRAC profiles left residuals which resemble the lensed source image as seen by ALMA (see Fig. A1).
To test whether these residuals were the result of a poorly-estimated PSF or because the mid-IR does not closely follow the near-IR light, we repeated the procedure for three galaxies of similar flux and morphology selected at random within the field of view common to both the s -band and IRAC observations. Removing the model IRAC profiles left no significant residuals in any of the three galaxies we tested. We therefore assumed that the residuals left from removing the lens profiles were indeed the source image although we can not rule out the possibility that the light profile of the lens in the near-IR is different to that in the mid-IR in the ID141 system. To obtain a final measure of the lensed source IRAC flux, we summed the residuals within an annulus surrounding the observed residual flux as shown in Fig. A1. This resulted in fluxes of 8.7±2.2 Jy and 8.1±1.8 Jy at 3.6 m and 4.5 m respectively. The errors quoted here are estimated from the image noise and do not include any uncertainty arising from lens light removal.
To constrain the far-IR part of the SED fitted by MAGPHYS, we used the photometry as discussed in Section 4.2. The optical/near-IR part of the SED was constrained using the Keck s -band data. Since the lensed source was not detected in the lens-subtracted s -band image, we set the flux in this band to be three times the standard deviation of the background flux and assumed a generous 1 error equal to this value. In this way, the s flux acts like an upper limit in the SED fitting. Finally, we also included Wide-field Infrared Survey Explorer fluxes at 12 m and 22 m (WISE; Lang et al. 2016). Since the WISE fluxes are an unknown combination of the lens and source light, we attributed 50 per cent of the flux in each passband to the source and took a 1 error equal to this value so that again, the WISE data act like upper limits. To convert the quantities derived by MAGPHYS from observed (i.e., lensed) to intrinsic quantities of the unlensed source, we assumed a magnification of 5.8, equal to the mean of the continuum magnifications given in Table 3. The resulting stellar mass we obtained is 3.4 +1.8 −1.7 × 10 11 M where the quoted 1 errors are taken from the distribution of masses output by MAGPHYS. Fig. A2 shows the corresponding MAGPHYS SED.

APPENDIX B: CHEMICAL EVOLUTION MODELLING
In this Section, we describe a one-zone chemical evolution model that we developed to compare the evolution of the interstellar dust, gas and stellar mass of a ID141-like system. The model is based on De Vis et al. (2017): the galaxy starts as a cloud of gas, which is converted into stars via the initial mass function and ongoing star formation. It includes a realistic description of the amount of gas flowing in (accretion from cosmic web) and out (supernovae and AGN) from the galaxy. The dust prescription accounts for dust formation by stars (both low-intermediate mass stars in their asymptotic giant branch phase and massive-star supernovae). It also includes dust grain growth in the diffuse and dense environments of the galaxy, thought to be a dominant production route of dust in galaxies above a metallicity threshold (Asano et al. 2013;Dunne et al. 2011;Zhukovska 2014;De Vis et al. 2017). Dust is removed from the ISM via astration (locked up in stellar remnants), dust destruction by SN shocks and outflows.
As detailed in the main text, we used two star formation histories, one with a fixed continuous star formation rate and one with a starburst. For the latter, we assumed that we are currently observing a starburst at 40 per cent of the way through a burst lasting 200 Myr. This is motivated by the fact that the duration of high-starbursts are thought to range from 50-250 Myr (e.g. Danielson et al. 2017); typical duty cycles of high-redshift SMGs range from 40-100 Myr (Coppin et al. 2008;Toft et al. 2014;Narayanan et al. 2015).
The three chosen IMFs adhere to the form / log ∝ − where is stellar mass. The Salpeter IMF is described by = 1.7. Both the Chabrier and top-heavy IMF of Zhang et al. (2018) are broken power laws with exponents Chab = [0.4, 1.35] and TH = [0.3, 1.1] (see Table B1). Both IMFs are normalised to unity between stellar masses of 0.1 M and 100 M .
Following Zhukovska (2014), we set the gas inflow rate to ( ) ∝ − / where is the infall timescale. We normalise this such that the galaxy accretes a total mass of gas acc over its lifetime (see below). To check that this is a realistic level of inflow, we compared  Table B1. The IMFs used in the model to show extremes of the gas and dust evolution. The IMF slopes are in terms of mass, and 0 and 2 are the upper and lower limits of stellar masses for which the IMF is normalised to unity. 1 indicates the mass at which the slope may change. The top-heavy IMF used here is from Zhang et al. (2018), see also Romano et al. (2017); Cai et al. (2020). This is a less extreme form than the Baugh et al. (2005) top-heavy IMF, where = 0 across the mass range.
with the halo baryonic mass accretion rate expected for a massive star-forming galaxy such as ID141 in the following way. Extreme star-forming galaxies with SFRs > 1000 M yr −1 at ∼ 4 exist in dark matter halo masses of a few ×10 12 M (Behroozi et al. 2013). is the baryonic fraction (0.17) and in denotes an efficiency. As an example, a galaxy with ℎ ( = 0) = 1 × 10 14 M (Narayanan et al. 2015) would have a mass accretion rate of ext ∼ 700 M −1 at = 4.2, and ∼ 400 M −1 at = 7. A current day halo mass of 10 13 M would have a mass accretion rate of roughly 8 times lower at these redshifts 5 . For the top heavy IMF with conSF SFR prescription, accretion of a total mass of 5 × 10 11 M of gas over the lifetime of the galaxy in our model, produces very similar inflow rates to those predicted from the halo mass growth described above. For the Salpeter IMF, The unattenuated spectrum is shown as the blue line. Red points and black circles show observed and predicted fluxes respectively. The bottom plot shows the residuals calculated as observed minus predicted flux. Note that all photometry in this plot is as observed, i.e., not lens de-magnified. the mass accretion rates need to be approximately twice as high to sustain the higher SFRs.
For the star formation histories, we also account for the different redshifts at which different mass halos undergo their peak SFE. Generally, this is at higher redshifts for higher mass halos. We scale the star formation efficiency by redshift such that the peak efficiency occurs at = 3 − 5 but decreases at higher and lower redshifts. For the halo described above, Behroozi et al. (2013) (their Fig. 8) shows that this varies from (0.4 − 1) × ( = 4). We fit a polynomial to their function.
We set the dimensionless parameter controlling the rate of dust grain growth in the dense interstellar medium to 7000. Recent work by De Vis et al. (2017) and De Vis et al. in review suggests an appropriate range of values for the rate of grain growth is 4000-7000. As discussed in the main text, this was allowed to vary from 700-15,000 in our analysis of ID141. In addition, we allow 60 per cent of the metals in the cloud to be 'available' for dust formation in order to produce a higher dust mass (lower gas to dust ratio). This factor is predicted to be ∼0.5 in nearby spirals (De Vis et al. in press). The grain growth parameter depends on the galaxy metal mass fraction, the current dust-to-metal ratio and the SFE, as well as the fraction of the ISM in dense clouds (we set this to 0.6, Inoue 2003; Yoshida et al. 2002). Finally, we let the galaxy evolve until (i) the stellar mass reaches 5 × 10 11 M where we assume an outflow clears the galaxy of the ISM (Romano et al. 2017), (ii) in the case of the starburst model, we stop the model when the burst has ended or (iii) when the gas fraction falls to less than 3 per cent.