High resolution spectral imaging of CO(7-6), [CI](2-1) and continuum of three high-z lensed dusty star-forming galaxies using ALMA

Our sample of three galaxies was selected from an initial program (PIs: Béthermin and Strandet, project number: 097.A-0973) targeting CO(7-6) and [CI](2-1) in 8 DSFGs from the SPT-SMG sample with the APEX/SEPIA instrument (on-source observation time are 492 mins for SPT0103-45, 559 mins for SPT2147-50 and 412 mins for SPT2357-51), a precursor of the ALMA band-5 receivers (Belitsky et al., 2018).

We reduced the data using the CLASS package of the GILDAS software¹¹1https://www.iram.fr/IRAMFR/GILDAS/. The data quality was estimated using an automatic procedure (Zhang et al. in prep.) and the bad data were removed from the analysis. We produced spectra combining all the other data after masking the region of the line and subtracting the baselines.

We use an antenna gain of 38 Jy K^-1 to obtain the APEX fluxes. We first estimate the line fluxes from the APEX spectra by summing the fluxes for all the channels in a manually-defined narrow range of velocities around the line peak. The noise is estimated by computing the standard deviation of the signal ($S_{\nu}\Delta v$ in Jy km/s) in line-free channels multiplied by $\sqrt{\rm N}$, where N is the number of channels used to estimate the line flux. [CI](2-1) is tentatively detected (3.24 $\sigma$) in SPT0103-45, but not CO(7-6). No line is detected in SPT2147-50. For SPT2357, we tentatively detect CO(7-6) at 4.1 $\sigma$, but not [CI](2-1). The results are summarised in Table 4 and are compared to ALMA fluxes in Sect. 4.2.

In the APEX sample, 5 out of 8 DSFGs have at least one line tentatively detected. From these, we selected the sources with a reliable lens model to aid the interpretation (Spilker et al., 2016). Three objects match our criteria: SPT0103-45, SPT0125-50 and SPT2357-51. These sources exhibit very diverse properties. SPT0103-45 has a bright [CI](2-1) line, SPT0125-50 is brighter in CO(7-6), whereas for SPT2357-51, the CO(7-6) line is unusually broad in the APEX/SEPIA spectra.²²2Note that the [CI](2-1) line was also tentatively detected and broad in the initial data reduction of the APEX spectra used to write the ALMA proposal. The APEX spectra are shown in Fig. 3 (grey dotted line). We also required the sources to be observable with ALMA receivers. Unfortunately, the lines of SPT0125-50 were between the ALMA Band 4 and Band 5 frequencies and thus could not be observed. We thus observed two sources from this sample, SPT0103-45 and SPT2357-51 at high-resolution with ALMA, with the Band-5 receivers (PI: Béthermin, 2017.1.01018.S). The target sensitivity on these observations is derived from the flux estimated on the APEX spectra (see, Sect. 2.1).

Our third source, SPT2147-50, was not detected with APEX despite the second data reduction, showing a weak signal for CO(7-6) (Fig. 3). However, SPT2147-50 is a part of the Targeting Extremely Magnified Panchromatic Lensed Arcs and Their Extended Star formation (TEMPLATES) James Webb Space Telescope (JWST) Early Release Science program³³3https://www.stsci.edu/jwst/observing-programs/approved-ers-programs/program-1355. The source was thus observed with ALMA in CO(7-6), [CI](2-1) and dust continuum to aid in the interpretation of the upcoming JWST data (PI: Vieira, 2018.1.01060.S).

Figure 1 shows the line luminosity of [CI](2-1) and CO(7-6) against the IR luminosity. We compare our sources with the sample compilation presented in Valentino et al. (2020). Our sources are at the bright end of the SMG cloud. Our final sample thus consists of three sources: SPT0103-45, SPT2147-50 and SPT2357-51. Our sources are presented in Table 1.

All our targets were observed with the ALMA 12-metre array using the Band-5 receivers and in the C43-5 configuration (shortest baseline: 15 m; longest baseline: 1400 m). The requested angular resolution was 0.3” for SPT0103-45 and SPT2357-51 and 0.4” for SPT2147-50. This was the best possible angular resolution allowed in Band 5 during Cycle 5 and allows us to resolve the sources over several synthesised beams.

To observe CO(7-6) and [CI](2-1), two spectral windows with 1.875 GHz bandwidth and 7.813 MHz resolution (corresponding to 11–13  km s^-1 velocity resolution) were used⁴⁴4 SPT2147-50 has only one spectral window centered on the lines with 7.813 MHz resolution. The second spectral window in the sideband has a coarser resolution (31.25 MHz).. For each source, we also placed two other spectral windows in the other sideband to measure the continuum. The same spectral resolution was used for SPT0103-45 and SPT2357-51, while 31.25 MHz was used for SPT2147-50.

To characterise the dynamics of these systems, we need to securely constrain the position of the peak of at least one line in each synthesised beam. For SPT0103-45 and SPT2357-51, we chose the target sensitivity to obtain a 10$\,\sigma$ detection when integrated over the full line width. We estimated the signal and the line width from the APEX spectra (Sect. 2.1) and divided the flux into the number of synthesised beams covered by the source from the previous 870 $\mu$m continuum imaging and lens model (Spilker et al., 2016). We derived a target sensitivity of 80 $\mu$Jy beam^-1 in a 320 km s^-1 channel for SPT0103-45 and 67$\mu$Jy beam^-1 in a 510 km s^-1 channel for SPT2357-51. For SPT2147-50, the sensitivity was estimated by rescaling the observed continuum profile to the estimated CO and [CI] flux, aiming for a 6$\,\sigma$ detection over 1/3 of the line width at the peak of the detection, which is equivalent to 10.4$\,\sigma$ over the full line width. This corresponds to 100 $\mu$Jy beam^-1 in a bandwidth of 100 km s^-1 (1/3 of the line width) or 58 $\mu$Jy beam^-1 in a bandwidth of 300 km s^-1. The target sensitivities for the three sources are thus rather similar in the end. The details of the observations are presented in Table 2.

Initial calibrations were done by the observatory using the standard ALMA pipeline based on the Common Astronomy Software Applications (CASA) (McMullin et al., 2007). Since our sources are bright, we attempted to self-calibrate the phase solutions. However, many long baselines were poorly fit and were thus flagged. The beam size was consequently degraded without a significant gain in sensitivity, so we discarded the self-calibration.

To image the data, we used the CLEAN routine in CASA. The data were imaged using a pixel size of 0.04^′′ in order to finely sample the synthesised beam ($\sim$10 pixels per beam). A Briggs weighting (Briggs, 1995) of the visibilities, with a robust parameter of 0.5 was chosen to give the best compromise between S/N and angular resolution.

To measure the continuum, we combine the two spectral windows in the line-free sideband and the channels without line contamination in the sideband where CO(7-6) and [CI](2-1) are located. These channels are identified using a first preliminary line imaging (see Sect. 3.3). The continuum images are then generated using a multi-frequency synthesis (MFS) method (Conway et al., 1990). A first imaging is made to estimate the noise level (hereafter $\sigma_{\rm noise}$). We then perform a second imaging using 3 $\sigma_{\rm noise}$ from the previous imaging as the Clean threshold. The first row of Fig. 2 shows the resulting continuum maps of our sources.

Using similar parameters as for the continuum, we first image the spectral window(s) covering CO(7-6) and [CI](2-1). We then subtract the continuum directly in the uv plane using the uvcontsub CASA task. The channels that are line-free in our initial image cube are used to fit and subtract the continuum. We assume a constant continuum level (as a function of the frequency) for SPT2357-51 and a first-order polynomial for the two other sources as the first-order polynomial gave a better fit of the continuum. We then re-imaged the data to obtain a continuum-subtracted line cube. No rebinning from the native $\sim 12$ km/s spectral resolution was used in the process, since the S/N per channel is sufficiently high.

For the continuum, the synthesised beam major axis size is between 0.43^′′ and 0.48^′′, while the minor axis size is between 0.35^′′ and 0.41^′′. The axis ratio varies from 1.17 to 1.21.

Concerning the line cubes, for SPT0103-45 and SPT2147-50 the synthesised beam is slightly larger than for the continuum, with the major axis ranging from 0.47^′′ to 0.48^′′ and the minor axis ranging from 0.39^′′ to 0.42^′′. This is because the spectral lines are located in the lower sideband of the frequency coverage. Conversely for SPT2357-51, the continuum has a coarser resolution because the lines are in the upper sideband, so the continuum-only data is at a lower average frequency. The beam axis ratio varies from 1.12–1.22. Overall, all three sources have similar, nearly circular beams. We estimate the noise level at the phase center by computing the standard deviation in the non-primary beam corrected map after masking the source. The average continuum sensitivity varies from 13–25 $\mu$Jy beam^-1. The same procedure is used channel by channel for the line cubes. The mean sensitivity per channel and the range of sensitivities per channel for every source are summarised in Table 3.

The dust continuum images are shown in the upper panels of Fig. 2.

SPT0103-45 is ring-shaped, which is characteristic of strongly lensed sources (Spilker et al., 2016). The brightest region is located in the northwest (henceforth, the gem) and another bright region can be seen in the southeast. They correspond to two counter-images of the same region of the source (see Sect. 6).

SPT2147-50 exhibits an extended arc in the northeast and a compact counter-image in the southwest (see Sect. 6). The northeast arc is brighter in its most eastern part and there is a small gap with no detected continuum at the most northern part of the arc.

SPT2357-51 is also gravitationally lensed and has a very compact continuum emission compared to the other sources. The compact emission is due to the small Einstein radius of the foreground lens galaxy ($\theta_{\rm E}\sim 0.2^{\prime\prime}$, Spilker et al. 2016). The resolved structure of SPT2357-51 is difficult to see visually in our default continuum image using all the line-free channels and a Briggs weighting with a robustness parameter of 0.5. In order to demonstrate the lensing geometry of this source more clearly, we re-imaged the continuum selecting only visibilities on baselines $>$200 k$\lambda$, shown in Fig.  18. We note that this image is for demonstration purposes only, as our lens modelling procedure (Sect.  6.2) fits directly to the interferometric visibilities and is not influenced by choices in imaging parameters.

The total continuum flux density is estimated from these maps by integrating the flux in an elliptical region encompassing the entire source. The continuum fluxes are tabulated in Table 4. To estimate the uncertainties of these fluxes, we select a polygonal region near the source and measure the noise RMS in this region. We then re-scale this noise by the factor of $\sqrt{\rm N}$, where N is the number of synthesised beams in the elliptical region where the source flux is estimated. We combine this statistical uncertainty with the estimated 10$\%$ absolute flux calibration typical of ALMA data.

We extract the integrated spectra from the line cubes after manually selecting the emission region. This process was repeated multiple times to ensure the robustness of the method. Figure 3 shows the integrated spectra of all our sources. There is a clear detection, S/N ¿ 14, for both the CO(7-6) and [CI](2-1) lines in all sources.

The line profiles of SPT0103-45 show a well-defined peak along with a redshifted tail in both CO and [CI] (Fig. 3, top panel). These tails extend to about $+$1000  km s^-1. The uni-directionality of the high-velocity tail is interesting and will be further explored in Sect. 5.1. The [CI](2-1) tail is slightly blended with the CO(7-6) peak emission and cannot be well isolated. We assessed the impact of this line blending by fitting the spectrum of SPT0103-45 with the sum of three Gaussian profiles each, for CO and [CI], fixing the center velocity relative to each line’s systemic velocity and the width of each pair of components (i.e. with a total of 12 free parameters to fit both lines simultaneously: three center velocities, three Gaussian line widths and six amplitudes). We find that accounting for the line blending increases the [CI]/CO line ratio by $\sim$5% compared to the simple division illustrated in Fig .3. We increase the uncertainty on the measured line ratio by this amount, added in quadrature to the statistical uncertainty from the noise in the data. The deblending procedure and results are further demonstrated in Appendix. B.

In SPT2147-50, both the lines are broad (FWHM $\sim$ 550 km s^-1) and asymmetric, with the redshifted side brighter than the blueshifted side of the spectra (see Fig. 3, middle panel).

For SPT2357-51, the line profiles show double peaks in both lines (see Fig. 3, bottom panel). The blueshifted peak is slightly brighter ($\sim$2 mJy) than the redshifted peak. The nature of this double peak feature will be further explored in the Sect. 5.3.

In Fig. 3, we compare the ALMA spectral profile with the APEX/SEPIA observations. We also re-estimate the APEX line fluxes using a broader integration window corresponding to the entire width of the line seen by ALMA (blue and yellow shaded regions of Fig. 3). The results are tabulated in Table 4. While the shape of the APEX and ALMA spectra generally agree, there are differences up to 2$\,\sigma$ in the integrated line fluxes. It is unlikely to be an interferometric effect (e.g., filtered-out large-scale emission) since the largest difference is found for SPT2357-51, the most compact source of our sample. It is more likely that the APEX spectrum is simply of lower quality than the ALMA spectra.

We construct integrated intensity (i.e. moment-0) maps of the [CI](2-1) and CO(7-6) lines by summing the flux of every channel covering the line using the CASA immoments task. We define a frequency range for each line as the integration window to compute the integrated intensity from the continuum-subtracted datacubes. The integration windows are shown in Fig. 3, with CO(7-6) in yellow and [CI](2-1) in blue. For SPT0103-45, the most redshifted CO(7-6) emission is blended with the most blueshifted [CI](2-1) emission, as noted above. We thus use a narrower range to integrate the [CI](2-1) emission in order to avoid the contamination of CO(7-6) emission, but note that the line blending introduces some unavoidable uncertainty into the moment maps.

The moment-0 maps for the [CI](2-1) and CO(7-6) lines are shown in the second and third rows of Fig. 2 respectively. The line emission is more extended than the continuum emission for all sources. The compactness of the dust continuum with respect to spectral lines that trace cold gas is very common in DSFGs (e.g., Spilker et al., 2015; Calistro Rivera et al., 2018; Dong et al., 2019; Apostolovski et al., 2019).

We estimate the line intensities of the sources by selecting the region encompassing the source manually and obtaining the intensities in this region. The uncertainties are estimated similarly to the continuum emission (Sect. 3.2). The intensities are tabulated in Table 4.

The two most intriguing features in this source are the high-velocity redshifted tail visible in both lines (see Fig. 3, top panel) and the presence of a particularly bright continuum feature, which we call the ‘gem’.

Such high-velocity tails are often seen as signatures of molecular outflows driven by AGNs (e.g., Maiolino et al., 2012; Cicone et al., 2014; Florez et al., 2021; Pantoni et al., 2021) and/or supernovae (Ginolfi et al., 2020; Spilker et al., 2020b, a). However, we see high-velocity tails only in one direction which makes it unlikely to be a signature of an outflow. We further explore other explanations such as magnification effects that could result in distorting the line profiles (see Sect. 6.7).

To try to understand the nature of SPT0103-45, we extract the spectra at the pixel where the continuum flux density is maximal, i.e. at the position of the gem. Figure 4 shows the spectrum at the position of the gem, clearly showing very broad and blended lines. The high-velocity tail is almost as bright as the peak of the emission near the systemic velocity. We also reconstruct the integrated intensity map of the high-velocity emission of CO(7-6) using the integration range shown in yellow in the center panel of Fig. 4. The 3, 5 and 10 $\sigma$ contours of the continuum are also plotted on the high-velocity intensity map shown in the bottom-right panel of Fig. 4. This reveals a small spatial offset between the high-velocity regions and the brightest continuum region.

We then build the velocity map of the source. To avoid artifacts caused by magnification or high-velocity components, we mapped the peak velocity of each line for each line of sight instead of a moment-1 map. To do that, we extract the spectrum at each pixel and fit a linear spline to reduce the impact of the noise. We then identify the peak of emission and the velocity corresponding to it. We search for the lines in the 197.5–198.1 GHz range (727.1 to -257.5 km s^-1) for [CI](2-1) and 196.5–197.38 GHz (866.2 to -45.6 km s^-1) for CO(7-6). These intervals were chosen to avoid cross-contamination between the two lines in the region, where low-velocity CO(7-6) and high-velocity [CI](2-1) could overlap. To ensure that this peak is real, we require that the maximal flux determined by the spline fit should be above 3 $\sigma$.

In Fig. 5, we show maps of the mean velocity of each line (i.e. moment-1 maps). To allow a better visualisation of the geometry of the source, the 3, 5 and 10 $\sigma$ contours of the continuum are plotted on the velocity maps. For both lines, we observe a very similar velocity gradient from east to west. However, the high-velocity regions corresponding to the redshifted tail are all located west of the gem.

To further analyse the kinematics of the source, we construct position-velocity (PV) diagrams along the bright arc passing through the gem. The slice is indicated as the white line in the left panel of Fig. 4. The PV diagram is shown in Fig. 6. The zero velocity corresponds to the observed frequency of the CO(7-6) line. The PV diagrams of both lines are similar in profile. From a position 0 to 1.5 arcsec, the velocity is almost flat followed by a quick increase between 1.5 and 2 arcsec and finally a second plateau around $\sim$800 km s^-1. The second plateau is much fainter, which could be due to differential magnification (see Sect. 6.4).

The [CI](2-1) and CO(7-6) lines of SPT2147-50 are broad and asymmetric with more flux on the redshifted side (see Fig. 3). To understand the origin of the asymmetric and broad lines, we construct a pixel-wise velocity map.

Since the lines are not always symmetric even in a given pixel, we locate the line peak velocity using a similar approach as for SPT0103-45 (see Sect. 5.1). As only one single peak per line can be found along any given line of sight, we define a single velocity range from -400 to +400km s^-1 to search for the peak.

These peak velocities are then used to construct the velocity maps seen in Fig. 7. We also plotted the 3, 5 and 10 $\sigma$ contours of the continuum emission to better visualize the geometry of the source. The velocity maps measured independently for each line are strikingly similar. The peak velocity varies from $-$300 to 150 km s^-1. We observe a smooth velocity gradient from east to west in the northern arc and the opposite in the south arc, as expected for a gravitationally-lensed arc and counterimage (see further discussions in Sect. 6.5).

The line profile of both lines has a doubly-peaked profile (see Fig. 3). A similar profile was observed for the [CII] line of SPT0346-52, which has been identified as a major merger (Litke et al., 2019). To assess whether SPT2357-51 may also be a merging system, we first examine the position-velocity diagram extracted along the right ascension axis (Fig. 8), chosen to maximize the spatial separation between the two components. It shows two clearly-separated components: a western component at blueshifted velocities (henceforth, the blue component) and an eastern component at redshifted velocities (henceforth, the red component). Their peaks of emission are separated spatially by about $\sim$1^′′ and spectrally by $\sim$400km s^-1.

To characterise both components seen in the position-velocity diagram, we performed a pixel-wise decomposition of the spectra into two Gaussian components for each line.

We extract the continuum-subtracted spectrum for each pixel and fit it with four Gaussian profiles (two for each line). To ensure a physical and meaningful result, we applied several constraints. The amplitude is forced to be positive, since we do not expect to see CO(7-6) or [CI](2-1) in absorption in such a system. The line width is allowed to vary from 0 to 250km s^-1 to avoid unphysically broad components. Finally, we define a velocity range associated with each component in order to avoid overfitting one component with two Gaussians if one component is very faint. The bounds for the central velocity of the peaks were defined from the velocity separation seen in the position-velocity diagram. The red component is fitted in a velocity range between 0 and +250km s^-1 and the blue component between $-$150 and $-$300km s^-1 for both [CI](2-1) and CO(7-6) lines.

To derive uncertainties on the amplitude, width and position of our Gaussians, we measure the noise in each channel and take it into account with the fitting tool (curvefit from the scipy package). In the rest of the analysis, we only consider results for which the significance of the amplitude is higher than 3$\,\sigma$. The residuals are below 2.5$\,\sigma$ for every fit.

Figure 9 shows the velocity map for each of the components of the [CI](2-1) and CO(7-6) lines along with the 3, 5 and 10$\,\sigma$ contours from the continuum map. The velocity gradients obtained for CO(7-6) and [CI](2-1) are remarkably similar. The blue component shows a smooth velocity gradient suggestive of rotation. The velocity structure of the red component is less clear, showing essentially no coherent structure.

The continuum is centered on the red component while the blue component is an offset from the peak of the dust continuum. Since these two components have similar apparent line luminosities, this could indicate that the red component has a larger dust emission and/or stronger obscured star formation. This could originate from different gas metallicities or temperatures between the two components. This difference of ISM properties between the two velocity components is a clue that this object could be a merger between two galaxies. It is unlikely that differential magnification could cause the low dust emission of the blue component. Lensing is a purely geometric effect, so any given individual region must have the same magnification for the lines and the continuum. Further discussion on the geometry of this source based on lens modelling will be done in Sect. 13.

We use the lens modelling code visilens from Spilker et al. (2016). It directly models Fourier-plane interferometric visibilities instead of the inverted image plane to limit the biases caused by correlated noise generated by the imaging algorithm. This also helps us better account for the antenna delay and the mismatched absolute flux calibration that could arise due to combining multiple observations taken on different days, which is the case for SPT2147-50. The angular resolution of our data does not allow us to make a reliable pixelated reconstruction of the source-plane.

Our model describes the foreground lens with one singular isothermal ellipsoid (SIE) profile. The profile is parametrised by the following free parameters: the position of the lens relative to the phase center ($x_{L}$,$y_{L}$), the angular Einstein radius $\theta_{E,L}$, the ellipticity $e_{L}$ and the position angle of the major axis $\phi_{L}$. The background sources are modelled with Sérsic (1968) and/or Gaussian profiles, for both the continuum emission and each spectral channel. The position of the source ($x_{S}$,$y_{S}$) is defined relative to the lens position. The other parameters are total flux density ($S$), effective half-light major axis ($a_{S}$), axis ratio ($b_{S}$/$a_{S}$), Sérsic index ($n_{S}$) and position angle ($\phi_{S}$). We note that models with substantially more freedom in the source morphology, such as a pixellated source plane (e.g., Dye et al., 2015; Hezaveh et al., 2016) are not warranted for our moderate-resolution data. Visilens uses a Monte Carlo Markov chain (MCMC) fitting procedure to sample the parameter space (Spilker et al., 2016). The best-fit parameters are then determined after the chains converge. We finally produce model dirty images, high-resolution models and source plane models from the best-fit model (see Fig. 10).

To perform the lens modelling of the continuum, we select only the data from the line-free spectral windows. We then average the two orthogonal polarisations. We choose a first estimate of the lens and source parameters using the previous modelling performed on low-resolution data by Spilker et al. (2016) and define loose reasonable priors (e.g., positive fluxes, $0<n_{S}<4$). We use a single component for the source profile of SPT2147-50 and SPT2357-51 but two components for SPT0103-45, as the model gave us lower residuals. We use 50 chains of 10 000 usable steps and 1000 burn-in steps to run the MCMC chains and verify the chains converged. We also check that the marginalized parameter distributions do not return a high probability for values close from the prior limits. The best-fit lens parameters for all our sources are given in Table 7. These parameters are in agreement with the best-fit lens parameters from Spilker et al. (2016).

In order to model the line emission of [CI](2-1) and CO(7-6), we use the continuum-free visibilities described in Sect. 3.3. We divide the lines into frequency bins with widths that ensure that the line emission is detected at a peak S/N$\gtrsim$8 in each channel. For SPT2147-50 and SPT2357-51, we use 50 MHz bins ($\sim$80 km s^-1). For SPT0103-45, we use 100 MHz bins ($\sim$150 km s^-1) for the redshifted tail where the emission is faint and narrower bins of 40 MHz width ($\sim$50 km s^-1) at the line peak where the S/N is high.

Since the S/N of the continuum is much higher than any of the lines, we simply fix the parameters of the lensing potential for all velocity channels of each source to those found in the continuum modelling. The source parameter starting values are also derived from the continuum analysis, but these parameters are allowed to vary. We use 40 chains of 10 000 usable steps and 1000 burn-in steps for each of the runs. We verify that the chains converged for every run. The best-fit parameters for the line modelling of SPT0103-45, SPT2147-50 and SPT2357-51 are listed in Table 8, Table 9 and Table 10 respectively.

For the continuum of SPT0103-45, a single source component was insufficient to properly fit the visibilities, thus we use two Sérsic components to model the source. The best-fit model converged to a solution with an extended component ($n_{S}=0.994\pm 0.011$, see Table 8) and a second very compact component ($r_{eff}=0.055\pm 0.002$ arcsec). The maximum continuum residuals were $\sim$8.3 $\sigma$, approximately 4% of the peak S/N in the data. This illustrates how difficult it is to model our very high S/N continuum data even with two components and suggests that a future approach with more freedom in the source plane (e.g. pixellated modelling) could be more successful. The dirty image of the data and best-fit model are shown in Fig. 10 (top-right panel) together with a high-resolution model (neglecting the Fourier-plane sampling) and the source-plane emission. The best-fit parameters of both the components are given in Table 8. The intrinsic flux ratio of these two components is about 17, with the extended component being brighter. The nature of this second component is unclear and we cannot conclude based on our data if it is a clump in the disk, a compact neighbour, or even an artifact of the model to mimic an asymmetry of the disk in the source plane.

To model the CO(7-6) and [CI](2-1) line data in each velocity bin, a single Sérsic source was sufficient. We start our MCMC chains from the best-fit parameters of the extended continuum component. The highest residuals are ¡ 5$\,\sigma$, except for one velocity bin (centered at $\sim 136$ km s^-1) for which the maximum residual was $\sim 7\,\sigma$. In Fig. 21 and Fig. 22, we compare the best-fit model and the observed data for the [CI](2-1) line and the CO(7-6) line, respectively.

In Fig. 11, we plot the best-fit source position and size for every velocity bin for each of the lines. We can see a smooth velocity gradient across the source, consistent for both the lines. Emission at blueshifted velocities is close to the inner diamond caustic, while the most redshifted emission is near the outer caustic. The results thus favor a scenario in which SPT0103-45 is a rotating disk.

In the case of SPT2147-50, the source was observed on three different days (see, Table 2). At each step of the MCMC code, we separately compute the likelihood of the data corresponding to each day and then combined them. To take into consideration the different atmospheric conditions and the small calibration offsets that could occur in these observations, we introduce two multiplicative amplitude re-scaling as free parameters for the latter two observations. In simple terms, these parameters account for the fact that SPT2147-50 is detected with S/N$\gg$20 while the absolute flux scale of ALMA data is quoted by the observatory to be accurate at the $\sim 5-10$% level.

For both the continuum and the line models, we use a single component with a Sérsic profile. This relatively simple model is sufficient to fit our source. The continuum residuals are lower than 1.5$\,\sigma$ and the residuals for the line modelling are always below 4$\,\sigma$ in any velocity channel. In Fig. 10 (middle row), we show the continuum dirty image and the modelled dirty images together with the high-resolution model and the source plane model. The comparison between the model and the dirty image in each velocity channel is presented in Appendix 23 and Appendix 24.

The best-fit value of the Sérsic index for the continuum is 1.36$\pm$0.05 (see Table 9). This index, close to 1, is compatible with the hypothesis that the object is a rotating disk (see also Sect. 5.2). The best-fit parameters for the line and continuum modelling are presented in Table 9. The source plane model (see Fig. 10) shows that part of this disk-like profile is near the inner caustic. Significant differential magnification effects (see Sect. 6.7) can thus be expected.

In Fig. 12, we present the source-plane reconstruction of the CO and [CI] line emission. We see a clear shift of the source position with velocity. The continuum center of emission (in black) coincides with velocities nearest the systemic velocity. The eastern region is blue shifted, while the western region near the caustic is redshifted. The gradient is similar for both lines. Overall, the center of emission regions at the various velocities are well aligned. This provides an additional evidence that the source could a rotating disk.

For SPT2357-51, the decomposition of the spectrum along the various line of sight was suggestive of a two-component system (see Sect. 5.3). However, a single source-plane Gaussian component was sufficient to model the data and gave low residuals (1.72$\,\sigma$), likely due to the low spatial resolution of our data and the compact lensing geometry. The best-fit parameters of the continuum lens model are presented in Table 10. In Fig. 10 (lower panels), we present the results of our best-fit continuum model.

The same source model is used in each velocity bin. The largest residuals for each of these bins are $<$ 5$\,\sigma$. In Fig. 13, we plot the position of the source for different velocities. The geometry of this source is more disturbed than the two other sources. While, there is a velocity gradient for the velocities corresponding to the blue component, the geometry associated to the red component ($>$0 km s^-1) is much more disturbed. The kinematics from our source-plane reconstruction are only marginally consistent with the expectation for rotation, though it is hard to constrain the emission region for the red component because the signal is emitted very close to caustics. Similar to the image-plane kinematic analysis (Sect. 5.3), we also find that the continuum is offset from the main region of gas emission. Our source-plane reconstruction thus favors the scenario of a disturbed kinematics possibly linked to a merger.

From the results of the lens modelling for our sources, we find differences in the magnification between the continuum and the line emission as well as across the CO and [CI] line profiles. Figure 14 shows the magnification as a function of velocity for both CO(7-6) and the [CI](2-1) along with the continuum magnification (horizontal blue dotted line). We note that the systemic redshifts of the sources were measured from the apparent spectra (not the intrinsic spectra). We also plot the CO(7-6) and [CI](2-1) observed spectra along with the magnifications to differentiate the intrinsic features from those arising due to magnification effects.

For SPT0103-45, the agreement between the observed line profile (dashed lines) and the magnification profile (points and dotted lines) in Fig. 14, top panel, demonstrate that the velocity-dependent magnification is largely responsible for the line shape in this object.

For both lines, the magnification of SPT2147-50 increases toward redshifted velocities. The line profile of SPT2147-50 is similar to SPT0103-45 with respect to the magnification versus velocity, except that in this case the redshifted side is more magnified. The asymmetry of the integrated line profile is again likely due to differential magnification.

The red component of SPT2357-51, which was brighter in the continuum (see Sect. 5.3 and Fig. 9), is more magnified than the blue component. Since they have similar apparent line fluxes, the red component is thus intrinsically fainter.

For each source, we calculated the effective magnification ($\mu_{\rm eff}$) of each line. It is computed as the ratio between the sum of the apparent flux in all velocity bins and the sum of the intrinsic fluxes in the same bins:

where $\mu$ is the magnification in each of the bins, $\Delta\upsilon$ is the width of the bin and $\rm S_{\nu}$ is the intrinsic flux density. In our lens model, the magnification and the intrinsic flux are correlated because the total apparent flux density is very well constrained by the data. We thus do not estimate uncertainties on $\mu_{\rm eff}$ by combining the marginalised uncertainties on the input quantities. We instead calculate the effective magnification for all MCMC steps and compute the median and the standard deviation of the obtained values. The effective magnifications are tabulated in Table 5.

For SPT2147-50 and SPT2357-51, the CO(7-6), [CI](2-1) and the continuum magnifications are similar with variations $<$15%. There is thus no evidence for strong differential magnification on integrated flux measurements. For SPT0103-45, we find a larger difference between the [CI](2-1) and continuum magnification ($\sim$24%) and between the CO(7-6) and the continuum ($\sim$16%). Traditionally, luminosity ratios are used to mitigate magnification effects in lensed sources. This assumes implicitly that the luminosity of all the lines and the continuum are magnified in the same way. Our analysis suggests that it is a reasonable assumption in general, but can introduce small biases depending on the lensing configuration of the individual sources. In contrast, our results show that spectral line profiles can be significantly distorted by differential magnification. Differential magnification effects are more severe across the line profiles than between the dust continuum and molecular gas tracers.

Analysing the resolved line and continuum ratios can help to understand the ISM properties of our sources. The bolometric infrared luminosity between 8 and 1000 $\mu$m (L_IR) is a well-known tracer of the SFR of galaxies (Kennicutt, 1998). The [CI](2-1) line can be used as a tracer of the total molecular gas (e.g., Papadopoulos et al., 2004) and the high-J CO lines are a tracer of the dense molecular gas (e.g., Weiß et al., 2007). Combining these three tracers, we can map both the resolved star formation and gas properties in our objects. We continue our ratio analysis in the image plane. This avoids the introduction of any artifacts in the line ratio due to our parameterised models for the source-plane emission. As gravitational lensing conserves surface brightness, an image-plane analysis is also equivalent to the same analysis in the source plane, with the caveat that multiply-imaged regions of the source plane appear multiple times in the image plane.

To compute the continuum to line ratios accurately and avoid resolution effects, our continuum and line maps must have comparable resolutions. In Table 3, we notice that the continuum has a slightly smaller synthesised beam size than the lines. This is a result of the frequency coverage of our data, since most data used for the continuum imaging comes from the other sideband, 12 GHz higher in frequency, that samples the $uv$ plane at slightly higher spatial frequencies. Hence, we produce another set of continuum maps using only the line-free continuum channels from the same sideband as the lines to more closely match the synthesised beam size between lines and continuum.

To obtain physically meaningful ratios, we convert both continuum and line fluxes in each beam into apparent luminosities. For the continuum maps, we use the following conversion to derive the apparent infrared luminosity per beam ($\mu L_{\rm IR(beam)}$):

where $S_{\nu,\textrm{(beam})}$ is the continuum flux density per beam measured in any pixel of the continuum map, $S_{\nu\rm(total)}$ is the total continuum flux density integrated over the entire source and $\mu\,L_{IR}$ is the total apparent infrared luminosity determined by Reuter et al. (2020) from Herschel and ground-based photometry.

This conversion implicitly assumes that the shape of the spectral energy distribution (SED) is the same everywhere in the galaxy. This is an approximation. Using simulations, Cochrane et al. (2019) estimated the dust temperature ($T_{\rm dust}$) gradient in star-forming galaxies ($\gamma$ in $T_{\rm dust}=T_{\rm center}\,r^{\gamma}$,⁵⁵5Note that the factor relating to the spatial distribution of dust is referred to as $\beta$ in Cochrane et al. (2019). We chose to refer to it using the symbol $\gamma$ to avoid confusions with the dust emissivity index. where $r$ is the distance from the center in kpc). They found a range of $\gamma$ between -0.3 and -0.05. For the smallest $\gamma$ and the largest source (SPT0103-45), we obtain a temperature variation by a factor of 1.4 between the center and the half-light radius. Assuming a simple $L\propto T^{4}$ scaling, this corresponds to a factor of 3.5 in luminosity at constant flux. For $\gamma=-0.05$, the effect is much smaller, with only a 24 % luminosity difference.

For the line flux maps, we use the conversion formula from Solomon et al. (1992)

where $S_{\rm line(beam)}\Delta\upsilon$ is the value obtained from the line flux maps derived in Sect. 4.3, in the units of Jy km s^-1  beam^-1, $D_{L}$ is the luminosity distance of the source in Mpc and $\nu_{\rm obs}$ is the observed line frequency in GHz. The total line luminosities and IR luminosity corrected for magnification for all our sources are tabulated in Table 6. We use the effective magnifications of the continuum or the associated line computed in Sect. 6.7 and listed in Table 5. The uncertainties on the line luminosity ratio of SPT0103-45 are the combined uncertainties derived from the moment map luminosities and the luminosities estimated by deblending the lines (see Appendix. B).

We then create $\rm L_{[\rm CI](2-1)}/\rm L_{\rm IR}$ and $\rm L_{[\rm CI](2-1)}/\rm L_{\rm CO(7-6)}$ ratio maps, shown in Figure 15 using a 3$\,\sigma$ threshold on the IR and line luminosities. SPT0103-45 has a lower [CI]/continuum ratio and marginally lower [CI]/CO ratio at the position of the gem. In SPT2147-50, the northeast part of the main arc and its southern counter-image show a lower [CI]/continuum than the western part of the arc, while no obvious features are seen in the [CI]/CO ratio map.

The [CI]/continuum ratio of SPT2357-51 exhibits a higher value in the center than on the outskirts. This may be counter-intuitive, since galaxy sizes measured from gas tracers are usually larger than from the continuum (Spilker et al., 2015; Dong et al., 2019; Apostolovski et al., 2019; Fujimoto et al., 2020). This central region could correspond to the overlap of the two components of this candidate merger system (see Sect. 5.3). In this interacting region, the continuum luminosity could possibly be underestimated if the dust is warmer in the core of this object. We again see no clear pattern in the [CI]/CO ratio map. For this source, higher-resolution data would be necessary to make stronger conclusions, given the compact lensing configuration of this source.

We next compare the luminosity ratios of our sources with other lensed SMGs. Andreani et al. (2018) found L_CO(7-6)/L_[CI](2-1) ratios of 2.5 and 4.0⁶⁶6after converting from L^′ (K km/s pc²) to L (L_⊙) units. for their sample. Yang et al. (2017) found that the line ratio vary from 0.8 to 3.0 in their sample of Herschel-selected SMGs. These ratios are typically higher than we measure in our sources, i.e. 0.93$\pm$0.02 for SPT0103-45, 1.2$\pm$0.1 for SPT2147-50 and 1.6$\pm$0.1 for SPT2357-51. Our ratios are consistent with the ratios seen in local galaxies by Kamenetzky et al. (2014, from 0.4 to 2) and close to the value for the Milky Way center (0.91) measured by Fixsen et al. (1999). Our systems thus have less extreme line ratios than the Herschel-selected lensed galaxies of Andreani et al. (2018). Those authors claim that the L_CO(7-6)/L_[CI](2-1) ratio can be indicative of whether a system is disk-dominated (line ratio $\sim$1) or merger-driven (up to $\sim$10). According to the dynamics and the lensing reconstruction of our sources performed in Sect. 5, 6 and 8, SPT0103-45 and SPT2147-50 seem to be dominated by a rotating and extended component, consistent with their low [CI]/CO line ratios $\lesssim$1. In contrast, SPT2357-51 has a slightly higher ratio and is a possible merger candidate. These results are thus an additional indication that the L_CO(7-6)/L_[CI](2-1) ratio could be used to diagnose the kinematic nature of high-$z$ sources. However, it may be less accurate than initially suggested by Andreani et al. (2018) and larger samples with high-resolution kinematic information will be necessary to draw stronger conclusions.

To identify the physical origin of the variations of these ratios across our sources, we plot $\rm L_{[\rm CI](2-1)}/\rm L_{\rm CO(7-6)}$ against $\rm L_{[\rm CI](2-1)}/\rm L_{\rm IR}$ for each pixel of our ratio maps (Fig. 16) to explore the heterogeneity of the ISM in our sources. In Fig. 16, we compare our sources with sample compilation presented in Valentino et al. (2020).

Along with the ratios, in Fig. 16, we plot the radiation field and hydrogen volume density iso-contours from the photo dissociation region (PDR) modelling of Kaufman et al. (1999). The iso-density contours are mostly horizontal in this parameter space, indicating that the ratio between [CI](2-1) and CO(7-6) can be a proxy of the gas density. Similarly, the radiation-field iso-contours are almost vertical, indicating that the [CI](2-1)-to-$\rm L_{\rm IR}$ ratio is a proxy of the radiation field strength.

However, we should note that comparing the line ratios to PDR models does have caveats. There may be additional heating mechanisms that cannot be represented only with line ratios. The hot and dense PDRs can reproduce the line ratios but may not accurately reproduce the L_IR as a part of the L_IR may arise from the HII region. The internal properties and heating mechanisms of the sources may also be driven by the IR emission arising from the warm dust (Valentino et al., 2018).

The globally-integrated ratios of our sources suggest that they have similar radiation field and densities to the compilation presented in Valentino et al. (2020). SPT0103-45 has the lowest density and radiation field strength of our three targets, while SPT2357-51 is our most extreme target in both parameters.

Turning to the pixel-wise ratios, we see that SPT0103-45 has a large variation in radiation field strength with almost uniform density, a result of the wide range in [CI]/$\rm L_{\rm IR}$ across the source. The gem region has a slightly higher radiation field strength and density compared to the global average. This region is likely a compact star-forming region close to the center of the object (see Fig. 11). Inspection of Fig. 15 suggests that the regions closest to the gem have a stronger radiation field, possibly due to an elevated dust temperature or a high SFE producing more young stars and UV photons.

In SPT2147-50, we observe a fairly large dispersion in the gas density and the radiation field ($\sim 0.5$ dex). The kinematics and lensing reconstruction of the source (Sect. 5 and Sect. 6) suggest that the source could be a rotating disk. These variations of density could thus be caused by dynamical instabilities in this gas-rich system (Bournaud & Elmegreen, 2009; Bournaud et al., 2009; Zanella et al., 2015).

The core of SPT2357-51 dominates the flux of the source (Fig. 15) and the associated pixels are very similar to the globally-integrated ratios. In contrast, the outskirts of the system appear to have similar or lower densities, but higher radiation fields. Several explanations are possible. The center of the object could be dustier, thereby preventing UV photons from penetrating the clouds. There could also be a dust temperature gradient from the outskirts to the center, which could bias our estimate of L_IR and thus the radiation field strength (see Sect. 7.1).

Despite having a complicated morphology in the image plane (Fig. 2 and 4) and an asymmetric observed integrated line profile (Fig. 3), SPT0103-45 seems to be compatible with a rotating disk from our lensing reconstruction (Sect. 6.4), with the asymmetry of the line profiles arising from differential magnification (Sect. 6.7).

The intrinsic SFR of SPT0103-45, after correcting for the continuum magnification derived in Sect. 6, is $1900\pm 600\,\rm M_{\odot}yr^{-1}$. Furthermore, using our continuum size estimates of the extended component, we compute the SFR surface density ($\Sigma_{\rm SFR}$) and find $69\pm 20\,\rm M_{\odot}yr^{-1}kpc^{-2}$. This is similar to typical SFR surface densities observed in SMGs (e.g., Daddi et al., 2010b; Casey et al., 2014).

We can also estimate the molecular gas mass from CO(3-2) observations from Reuter et al. (2020). We assume a CO(3-2)/CO(1-0) line ratio of 0.7 $\pm$ 0.1 from Harrington et al. (2021) and use an $\alpha_{\rm CO}=3.4\,\,\rm M_{\odot}(K\,km/s\,pc^{2})^{-1}$ (where $\alpha_{\rm CO}$ is the CO-to-H₂ conversion factor) from Jarugula et al. (2021). We note that the quoted value and uncertainty on the line ratio we assume is derived from and applicable to high-redshift SMGs only, the uncertainty would increase if we also considered other galaxy populations when making this assumption. In any case, we expect the uncertainty on $\alpha_{\rm CO}$ to dominate over the uncertainty in CO(3-1)/CO(1-0) line ratio. We find M${}_{\rm gas}=(6.8\pm 2.1)\times 10^{11}\,\rm M_{\odot}$ and a gas depletion timescale $\rm t_{\rm dep}=360\pm 140$ Myr. Instead, assuming a lower $\alpha_{\rm CO}=0.8$ (Downes & Solomon, 1998), we find M${}_{\rm gas}=(1.6\pm 0.5)\times 10^{11}\,\rm M_{\odot}$ and a gas depletion time scale $\rm t_{\rm dep}=80\pm 40$ Myr. This short depletion timescale is typical of starbursts (e.g., Daddi et al., 2010b; Béthermin et al., 2015).

Furthermore, we estimate the dynamical mass of the system using a simple approach. Assuming a simple circular Keplerian motion, the enclosed dynamical mass $\rm M_{\rm dyn}$ can be computed using the following equation:

where R is the radius, which we take to be half of the distance between the lowest and highest velocity components in the source plane (see Table 8). We estimated the inclination $i$ of the system using the b/a ratio of the continuum model ($\cos i=b/a$). $\rm V_{\rm obs}$ is the projected component of the circular velocity along the line of sight and is computed as half of the velocity difference between the two extreme velocity bins. This approach is an intermediate between the crude estimate based on the line width and the size (e.g., Aravena et al., 2016; Bothwell et al., 2017; Dessauges-Zavadsky et al., 2020) and a full dynamical approach (e.g., Rizzo et al., 2020; Jones et al., 2021) done on higher resolution data or unlensed objects. Contrary to the crude approach, we can take into account the inclination of our objects and obtain the circular velocity in the outskirts of the objects instead of an integrated line width. Of course, it relies on the strong assumption that these systems are dynamically relaxed disks.

We estimate the uncertainties on the result by computing $\rm M_{\rm dyn}$ for all realizations of MCMC chains from the lens model, though we expect that systematic uncertainties in e.g. the assumed source geometry are dominant. We find a dynamical mass $\rm M_{\rm dyn}$ of $(1.2\pm 0.5)\times 10^{11}\,\rm M_{\odot}$ using the CO(7-6) radius and velocities and ($1.55\pm 0.4)\times 10^{11}\,\rm M_{\odot}$ using the [CI](2-1) radius and velocities.

We find compatible values between this dynamical mass and the molecular gas mass using the lower $\alpha_{\rm CO}=0.8$ value. Under this assumption there is only a small difference between the dynamical and gas masses, implying that the stellar mass of the source is relatively low. We note that the dynamical mass is calculated with the radius of the mid-J CO transition which requires higher excitation temperatures, typically found in the galactic centers. In Fig. 17, the depletion timescale of SPT0103-45 is similar to the bulk SMG population. However, the short depletion time and the high SFR density corresponds to a starburst-like star formation mode, which is usually not expected in isolated disks. An $\alpha_{\rm CO}=3.4$ ($\sim$4) leads to a longer depletion timescale usually found in disk galaxies, but this would also result in a strong tension with the dynamical mass. Only higher resolution data could allow us to identify a potential interaction triggering the starburst. We also cannot exclude processes linked to disk instabilities (e.g., Bournaud et al., 2007).

The dynamics of SPT2147-50 are similar to SPT0103-45, with our source plane modelling being compatible with rotation. SPT2147-50 has a lensing-corrected SFR of $830\pm 200\,\rm M_{\odot}\rm yr^{-1}$ and a $\Sigma_{\rm SFR}$ of $30\pm 7\,\rm M_{\odot}yr^{-1}kpc^{-2}$. Bothwell et al. (2017) estimated the molecular gas mass using [CI](1-0) and found $\rm M_{\rm mol}=(6.8\pm 0.5)\times 10^{10}\,\rm M_{\odot}$ (after updating their magnification by our [CI](2-1)-based measurement). Aravena et al. (2016) found a smaller molecular gas mass of $(2.3\pm 0.5)\times 10^{10}\,\rm M_{\odot}$ using CO(2-1) ($\alpha_{\rm CO}=0.8$, $\rm L^{\prime}_{\rm CO(2-1)}/\rm L^{\prime}_{\rm CO(1-0)}=0.8$ and updating the magnification using our CO(7-6)-based value). We derive a depletion timescale $\rm t_{\rm dep}=82\pm 7$ Myr using [CI](1-0) and $\rm t_{\rm dep}=28\pm 9\,\rm Myr$ using CO(2-1). This short depletion timescale is usually associated to a starburst. If instead we compute the gas mass using a Milky Way-like $\alpha_{\rm CO}=3.4$, we find a gas mass of $(9.9\pm 0.5)\times 10^{10}\,\rm M_{\odot}$ and a depletion time $\rm t_{\rm dep}=119\pm 29\,\rm Myr$.

We can also compute the dynamical mass of the system using the same method as SPT0103-45. We find $\rm M_{\rm dyn}=(2.5\pm 0.4)\times 10^{10}\,\rm M_{\odot}$ from the CO(7-6) data and $\rm M_{\rm dyn}=(3.4\pm 0.4)\times 10^{10}\,\rm M_{\odot}$ using [CI](2-1). This is very similar to a comparable dynamical estimate from Aravena et al. (2016) using only the line width and the continuum effective radius. We note that the estimate of the gas mass from CO(2-1) and using $\alpha_{\rm CO}=0.8$ seems to agree better with the dynamical mass than the [CI]-derived ones. In Fig. 17, we see that the depletion timescale is somewhat shorter compared to the bulk of the SMGs. Similar to SPT0103-45, using $\alpha_{\rm CO}=3.4$ leads to tension between the gas and dynamical masses.

The most interesting feature of this compact source is the double-peak profile of the spectrum (Fig. 3, bottom panel) seen for both lines. The PV diagram of the source (see Fig. 8), the image plane spectral modelling of the source with two components (Sect. 5.3) and the lens reconstruction (Sect. 6.6) suggest a merger between two components. The intrinsic flux ratio between the two components are 0.5:1 and 0.84:1 in CO(7-6) and [CI](2-1), respectively, with the blue component being brighter.

The blue component shows a smooth velocity gradient for both lines suggesting possible rotation or tidal streaming motions, whereas the red component exhibits a much more disturbed velocity structure. Comparing the positions of these components with the continuum emission, we notice that the red component is more centered on the continuum than the blue component. This could indicate a difference in their dust properties, with the red component being more dusty than the blue one. This could indicate a different chemical maturity of the two components which might be merging. However, it is surprising that the dustier component is the less massive one, since usually more massive galaxies are also dustier (Heinis et al., 2014; Fudamoto et al., 2020).

Using the magnification factor obtained from our continuum best-fit lens model and the L_IR from Reuter et al. (2020), we compute an intrinsic SFR of $1620\pm 260\rm\,M_{\odot}\rm yr^{-1}$. Based on the continuum source size from the lens modelling (Sect. 6.6), we estimate the SFR density $\Sigma_{\rm SFR}$ to be $59\pm 10\,\rm M_{\odot}yr^{-1}kpc^{-2}$, which is similar to the other sources and typical of SMGs.

We also estimate the gas mass of the system from the CO(3-2) flux measurement from Reuter et al. (2020). Using a CO(3-2)/CO(1-0) conversion factor of 0.7 $\pm$ 0.1 from Harrington et al. (2021) and $\alpha_{\rm CO}=3.4$, we obtain $\rm M_{\rm gas}=(1.6\pm 0.5)\times 10^{11}\,\rm M_{\odot}$ and a depletion timescale $\rm t_{\rm dep}=99\pm 29$ Myr. Alternatively, using an $\alpha_{\rm CO}=0.8$, we obtain $\rm M_{\rm gas}=(3.9\pm 1.3)\times 10^{10}\,\rm M_{\odot}$ and a very short depletion timescale $\rm t_{\rm dep}=24\pm 9$ Myr. This is much lower than the depletion timescales of the bulk of the SMG population shown in Fig. 17. Such short depletion times have also been measured in some extreme lensed SMGs but remain unusual (Béthermin et al., 2016; Ciesla et al., 2020). This reinforces the hypothesis that SPT2357-51 is an extreme starburst possibly triggered by a major merger. We do not calculate the dynamical mass of the system as we cannot assume that it is a dynamically relaxed system.