# Molecular Line Observations in Two Dusty Star-forming Galaxies at z = 6.9

Sreevani Jarugula, Joaquin D. Vieira, Axel Weiss, Justin S. Spilker, Manuel

Aravena, Melanie Archipley, Matthieu Béthermin, Scott C. Chapman,

Chenxing Dong, Thomas R. Greve, et al.

► To cite this version:

Sreevani Jarugula, Joaquin D. Vieira, Axel Weiss, Justin S. Spilker, Manuel Aravena, et al. Molecular Line Observations in Two Dusty Star-forming Galaxies at z = 6.9. The Astrophysical Journal, 2021, 921 (1), pp.97. 10.3847/1538-4357/ac21db. hal-03577235

HAL Id: hal-03577235

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

Submitted on 19 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Molecular Line Observations in Two Dusty Star-forming Galaxies at z = 6.9

Sreevani Jarugula¹, Joaquin D. Vieira¹,²,³, Axel Weiss⁴, Justin S. Spilker⁵,<fnref n="21" />, Manuel Aravena⁶, Melanie Archipley¹,³, Matthieu Béthermin⁷, Scott C. Chapman⁸, Chenxing Dong⁹, Thomas R. Greve¹⁰,¹¹, Kevin Harrington¹², Christopher C. Hayward¹³, Yashar Hezaveh¹³,¹⁴, Ryley Hill¹⁵, Katrina C. Litke¹⁶, Matthew A. Malkan¹⁷, Daniel P. Marrone¹⁶, Desika Narayanan¹⁰,¹⁸,¹⁹, Kedar A. Phadke¹, Cassie Reuter¹, and Kaja M. Rotermund²⁰

¹Department of Astronomy, University of Illinois, 1002 West Green St., Urbana, IL 61801, USA; jarugul2@illinois.edu

²Department of Physics, University of Illinois at Urbana-Champaign, 1110 W Green St Loomis Laboratory, Urbana, IL 61801, USA

³Center for AstroPhysical Surveys, National Center for Supercomputing Applications, Urbana, IL, 61801, USA

⁴Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69 D-53121 Bonn, Germany

⁵Department of Astronomy, University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712, USA

⁶Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

⁷Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France

⁸Eureka Scientific, Inc. 2452 Delmer Street Suite 100, Oakland, CA 94602-3017, USA

⁹Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

¹⁰Cosmic Dawn Center (DAWN), DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs. Lyngby, Denmark

¹¹Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK

¹²European Southern Observatory, Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile

¹³Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA

¹⁴Département de Physique, Université de Montréal, Montreal, Quebec, H3T 1J4, Canada

¹⁵Department of Physics and Astronomy, University of British Columbia, 6225 Agricultural Rd., Vancouver, V6T 1Z1, Canada

¹⁶Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

¹⁷Department of Physics and Astronomy, University of California, Los Angeles, CA 90095-1547, USA

¹⁸Department of Astronomy, University of Florida, 211 Bryant Space Sciences Center, Gainesville, FL 32611, USA

¹⁹University of Florida Informatics Institute, 432 Newell Drive, CISE Bldg E251, Gainesville, FL 32611, USA

²⁰Department of Physics and Atmospheric Science, Dalhousie University, Halifax, Nova Scotia, Canada

Received 2021 May 25; revised 2021 August 11; accepted 2021 August 24; published 2021 November 3

Abstract

SPT0311-58 is the most massive infrared luminous system discovered so far during the Epoch of Reionization (EoR). In this paper, we present a detailed analysis of the molecular interstellar medium at z = 6.9, through high resolution observations of the CO(6–5), CO(7–6), CO(10–9), [C I](2–1), and p-H₂O(2₁,₁−2₀,₂) lines and dust continuum emissions with the Atacama Large Millimeter/submillimeter Array. The system consists of a pair of intensely star-forming, gravitationally lensed galaxies (labeled West and East). The intrinsic far-infrared luminosity is (16 ± 4) × 10¹² L⊙ in West and (27 ± 4) × 10¹¹ L⊙ in East. We model the dust, CO, and [C I] using non-local thermodynamic equilibrium radiative transfer models and estimate the intrinsic gas mass to be (5.4 ± 3.4) × 10¹¹ M⊙ in West and (3.1 ± 2.7) × 10¹⁰ M⊙ in East. We find that the CO spectral line energy distribution in West and East are typical of high-redshift submillimeter galaxies (SMGs). The CO-to-H₂ conversion factor (αCO) and the gas depletion timescales estimated from the model are consistent with the high-redshift SMGs in the literature within the uncertainties. We find no evidence of evolution of depletion time with redshift in SMGs at z > 3. This is the most detailed study of molecular gas content of a galaxy in the EoR to date, with the most distant detection of H₂O in a galaxy without any evidence for active galactic nuclei in the literature.

*Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Starburst galaxies (1570); Ultraluminous infrared galaxies (1735)*

1. Introduction

Dusty star-forming galaxies (DSFGs) are dust-enshrouded, intensely star-forming galaxies with tens to thousands of stars forming per year (e.g., Casey et al. 2014). These starburst galaxies are bright at submillimeter wavelengths as the ultraviolet (UV) photons from the young stars are absorbed and reradiated by dust into the far-infrared (FIR) wavelengths, giving rise to total infrared luminosities ($L_{IR}$) greater than $10^{12} L_{⊙}$. Observations of these galaxies have the advantage of the negative-K correction (Blain & Longair 1993) at long wavelengths, enabling galaxy selection that is independent of redshift. DSFGs contribute significantly to the cosmic star formation history at high redshift (Casey et al. 2014) and are thought to be the progenitors of present-day massive

ellipticals (e.g., Simpson et al. 2014). However, the theoretical understanding of the DSFG population has been challenging (e.g., Narayanan et al. 2015; Hayward et al. 2021), especially during the Epoch of Reionization (EoR), between z ~ 15–6 (Madau et al. 1997).

To understand star formation in DSFGs, it is important to study the molecular gas content of the interstellar medium (ISM) (e.g., Carilli & Walter 2013). However, cold molecular hydrogen (H₂), which is the most abundant molecule and the fuel for star formation, is difficult to observe directly because it lacks a permanent dipole moment. The most commonly used tracer of H₂ is carbon monoxide (¹²C¹⁶O, hereafter CO), which is the second most abundant molecule in the universe and also an important coolant. The observations of multiple CO rotational transitions from local galaxies (e.g., Greve et al. 2014; Rosenberg et al. 2015; Lu et al. 2017), high-redshift ultraluminous infrared galaxies ((U)

<footnote n="21">²¹NHFP Hubble Fellow.</footnote>
<table><caption>**Table 1**
ALMA observations</caption><thead><tr><th>Tuning</th><th>Time on source<br>(minutes)</th><th>ν<sub>central</sub><br>LSB, USB<br>(GHz)</th><th>Line</th><th>ν<sub>rest</sub><br>(GHz)</th><th>FWHM<sub>cont</sub><br>(")</th><th>σ<sub>cont</sub><br>(μJy beam<sup>-1</sup>)</th><th>SNR<sub>line</sub><br>W, E</th></tr></thead><tbody><tr><td>1</td><td>47</td><td>95, 107</td><td>p-H<sub>2</sub>O(2<sub>1,1</sub> - 2<sub>0,2</sub>)</td><td>752.03</td><td>0.31 × 0.22</td><td>16.19</td><td>4, ...</td></tr><tr><td rowspan="3">2</td><td rowspan="3">43</td><td>88, 90</td><td>CO(6-5)</td><td>691.50</td><td rowspan="3">0.47 × 0.36</td><td rowspan="3">8.83</td><td>10, 4</td></tr><tr><td rowspan="2">100, 102</td><td>CO(7-6)</td><td>806.65</td><td>8, 2</td></tr><tr><td>[C I](2-1)</td><td>809.34</td><td>5, ...</td></tr><tr><td rowspan="2">3</td><td rowspan="2">46</td><td rowspan="2">134, 146</td><td>CO(10-9)</td><td>1151.98</td><td rowspan="2">0.27 × 0.20</td><td rowspan="2">8.27</td><td>6, 2</td></tr><tr><td>p-H<sub>2</sub>O(3<sub>1,2</sub> - 2<sub>2,1</sub>)</td><td>1153.13</td><td>...</td></tr></tbody></table>

**Note.** $\nu_{\text{central}}$ corresponds to the central frequency in the lowerside band (LSB) and the upperside band (USB) of the observations. $\nu_{\text{rest}}^{\text{line}}$ is the rest frequency of the lines. $\text{FWHM}_{\text{cont}}$ is the resolution of the continuum image. $\sigma_{\text{cont}}$ is the rms noise in the continuum image by combining LSB and USB. $\text{SNR}_{\text{line}}$ is the signal-to-noise ratio of the lines in W and E, which is calculated by using the maximum flux density of the scaled [C II] spectrum and the standard deviation of the flux density in the line-free channels (Section 3.1). $\text{SNR}_{\text{line}} < 2$ are not shown.

LIRGs) and quasars (e.g., Weiß et al. 2007; Yang et al. 2017), and simulations (Kamenetzky et al. 2018) have shown that the CO spectral line energy distribution (SLED) can probe the physical conditions of the galaxies, such as density and temperature. The low-$J$ CO transitions emitted from the low density ($n_{H_2} \lesssim 10^3 \text{ cm}^{-3}$) diffuse the ISM (e.g., Harris et al. 2010; Ivison et al. 2011). The mid-$J$ CO transitions, such as CO(6–5) and CO(7–6), are excited in denser regions of the galaxy ($n_{H_2} \sim 10^4 \text{ cm}^{-3}$) where star formation mainly occurs (Lu et al. 2015, 2017). The high-$J$ CO transitions ($J > 10$) originate in the compact high-density regions and may be further enhanced in the presence of active galactic nuclei (AGNs) (e.g., Weiß et al. 2007; Bradford et al. 2009; Lu et al. 2017). Molecular gas mass is traditionally calculated from the CO luminosity by using the CO-to-H$_2$ gas conversion factor, $\alpha_{\text{CO}}$ (e.g., Bolatto et al. 2013). However, there is a large uncertainty in the estimation of gas mass because of factors such as the dependence of carbon abundance on optical depth, metallicity, and the destruction of CO due to UV radiation (e.g., Narayanan et al. 2011; Bolatto et al. 2013). Observations have shown that [C I](1–0) can be an independent tracer of molecular gas (e.g., Papadopoulos et al. 2004; Weiß et al. 2005a; Walter et al. 2011; Bothwell et al. 2017), given its simpler excitation level structure.

Another abundant molecule in the universe after H$_2$ and CO is water, H$_2$O (Neufeld et al. 1995). The complex level populations of H$_2$O are tightly coupled to the infrared radiation field. This is due to infrared pumping, where the molecular transitions in H$_2$O are mainly excited by the FIR photons. The different transitions of H$_2$O have been observed to be correlated with $L_{\text{IR}}$; hence, star formation over more than three orders of magnitude, both in the local- and high-redshift (U)LIRGs (e.g., Omont et al. 2013; Yang et al. 2013, 2016; Liu 2017; Jarugula et al. 2019). H$_2$O is observed to be more tightly correlated to $L_{\text{IR}}$ than CO(6–5) or [C II] (Jarugula et al. 2019).

SPT0311-58 consists of a pair of gravitationally lensed galaxies discovered in the South Pole Telescope (SPT) Survey (Vieira et al. 2013; Everett et al. 2020). The spatially unresolved observations of low- and mid-$J$ CO transitions from the Atacama Large Millimeter/ submillimeter Array (ALMA), the Australia Telescope Compact Array (ATCA), and the Atacama Pathfinder Experiment (APEX) are presented in Strandet et al. (2017). These observations confirm the redshift of the source at $z = 6.9$. High-resolution dust, [C II], and [O III] ALMA observations of SPT0311-58 are discussed in Marrone et al. (2018), and a detailed lens modeling shows a pair of

<figure><img src="image_3.png" /><figcaption>**Figure 1.** SPT0311-58 West and East. The background gray-scale image is ALMA 423 GHz Band 8 continuum with magnification of ~1.3 in East and ~2.2 in West (Marrone et al. 2018). The moment 0 contours of CO(6–5) (red), CO(10–9) (magenta), p-H$_2$O(2$_{1,1}$ – 2$_{0,2}$) (green), and [C I](2–1) (cyan) from this analysis are overlaid on the continuum image. The CO contours are at [3, 4, 5, 7, and 9] × $\sigma$ and [C I] and H$_2$O contours at [2, 3, and 4] × $\sigma$. The synthesized beams are shown in the bottom left and the 1".0 scale bar, which corresponds to ~5.4 kpc at $z = 6.9$, is at the bottom right.</figcaption></figure>

galaxies separated by a projected distance of ~8 kpc (~1".5), likely in a state of merging. In this paper, we combine the previous observations with high-resolution molecular line imaging with ALMA, including the high-$J$ CO transition, CO(10–9) and p-H$_2$O(2$_{1,1}$ – 2$_{0,2}$). We use realistic non-Local Thermodynamic Equilibrium (non-LTE) models and estimate the physical properties of the ISM in both the galaxies.

Throughout the paper, we refer to SPT0311-58 West as W and SPT0311-58 East as E. In Section 2, we present the observations used in this analysis and the data reduction. In Section 3, we detail the data analysis procedure to estimate the continuum and line luminosities. The lens modeling is also discussed in this section. The results are presented in Section 4. In the first half of Section 4, we present the results from observations and in the second half, we discuss the radiative transfer models and results. In Section 5, we continue to
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: the continuum is shown as the background in square root scale with minimum pixel value as $3\sigma_{cont}$, where $\sigma_{cont}$ is the rms of the continuum map. The moment 0 contours of the lines are overlaid at [3, 5, 7, and 9] $\times \sigma$ in W and at [3 and 4] $\times \sigma$ in E, where $\sigma$ is the rms noise in the moment 0 map. Red and blue contours correspond to CO(6–5), CO(7–6), and CO(10–9) in the W and E sources, respectively. The [C I] $3\sigma$ contour (in the W source of second panel) is shown in white and the p-H$_2$O($2_{1,1} - 2_{0,2}$) $3\sigma$ contour (in the W source of last panel) is shown in magenta. The synthesized beam of the continuum image is shown in the lower left and the $1''/0$ scale bar is at the lower right. Bottom: the spectrum of each line is shown in W (blue) and E (green), shifted down for clarity. The scaled [C II] spectrum is shown in red for CO and p-H$_2$O($2_{1,1} - 2_{0,2}$) lines and in black for [C I] and p-H$_2$O($3_{1,2} - 2_{2,1}$). The scaling procedure is described in Section 3.1.</figcaption></figure>

discuss the results, and in Section 6, we conclude with a summary. We use flat $\Lambda$CDM cosmology where $h = 0.677$, $\Omega_m = 0.307$, and $\Omega_\Lambda = 0.693$ (Planck Collaboration et al. 2016). We estimate the $L_{IR}$ as flux integrated from 8–1000 µm and total FIR luminosity ($L_{FIR}$) from 42.5–122.5 µm in the rest frame (Helou et al. 1985).

## 2. Observations

We obtained observations of the different emission lines and the continuum using the ALMA Bands 3 and 4 over three tunings (Project ID: 2017.1.01168.S, PI: Vieira). Table 1 presents the details of the observations. We detect CO(6–5), CO(7–6), and CO(10–9) in both of the components, and detect [C I] (2–1) and p-H$_2$O($2_{1,1} - 2_{0,2}$) in W. In Figure 1, we show the high-resolution ALMA Band 8 continuum image from Marrone et al. (2018) and CO(6–5), CO(10–9), [C I](1–0), and p-H$_2$O($2_{1,1} - 2_{0,2}$) observations from this analysis (see also Figure 2).

### 2.1. Data Reduction and Imaging

The data reduction and the imaging were performed using the Common Astronomy Software Application (CASA) package (McMullin et al. 2007). We use the calibrated data products from the ALMA data reduction pipeline (CASA version 5.1.1). The continuum images are produced by combining the LSB and the USB, excluding the channels containing the line emission. Briggs weighting and a robust parameter of 0.5 are used, which provides a good compromise between resolution and noise. This gives a synthesized beam of $\sim 0''$.5 at 95 GHz and $\sim 0''$.3 at 140 GHz. To generate the spectral cubes, we use the same weighting as with the continuum maps and a velocity averaging of 100 km s$^{-1}$, after subtracting the continuum using the CASA task uvcontsub with a polynomial fit of order 1.

To generate the velocity-integrated intensity maps (i.e., moment 0), we produce single-channel cubes with a width of 1000 km s$^{-1}$ in W and 500 km s$^{-1}$ in E, which is $\sim 2 \times$ full width at half maximum (FWHM) of the lines. We consider a velocity separation of $\sim 750$ km s$^{-1}$ between W and E based on the high-resolution [C II] observations (Marrone et al. 2018). The synthesized beam of the moment 0 maps for p-H$_2$O($2_{1,1} - 2_{0,2}$) is $\sim 0''$.3, CO(6–5) and CO(7–6) is $\sim 0''$.5, and CO(10–9) is $\sim 0''$.2.

The signal-to-noise ratio of the lines and the flux densities of the continuum are given in Tables 1 and 2, respectively. The continuum images with moment 0 contours overlaid on top are shown in Figure 2.

### 2.2. Ancillary Data

In addition to the molecular line observations and the continuum flux density at 95 GHz (B3) and 140 GHz (B4), we include ancillary data on SPT0311-58 from the literature in this analysis (Strandet et al. 2016, 2017; Marrone et al. 2018). The source has been observed in previous ALMA cycles (Project ID: 2015.1.00504.S, PI: Strandet and Project ID: 2016.1.01293.S, PI: Marrone) in B3, B6, B7, and B8 corresponding to 95, 240, 340, and 420 GHz, respectively. For the continuum flux density in B3, we use current observations that are at a higher resolution of $\sim 0''$.5 compared to previous observations at $\sim 3''$.5. We also include Herschel/ SPIRE photometry at 250, 350, and 500 µm (Project ID: DDT_mstrande_1; Strandet et al. 2016). The intrinsic continuum flux densities (corrected for magnification) in W and E are obtained by performing the lens modeling described in detail in Marrone et al. (2018).

We include CO(3–2) data observed with the ATCA to constrain the CO SLED. The data reduction and the line
<table><caption>**Table 2**
ALMA Band 3 and Band 4 Continuum Properties</caption><thead><tr><th>Tuning</th><th>Source</th><th>ν<sub>obs</sub><br>(GHz)</th><th>S<sub>obs</sub><br>(mJy)</th><th>μ</th><th>S<sub>int</sub><br>(mJy)</th><th>L<sub>FIR int</sub><br>(10<sup>12</sup> L<sub>☉</sub>)</th><th>L<sub>IR int</sub><br>(10<sup>12</sup> L<sub>☉</sub>)</th></tr></thead><tbody><tr><td>1 and 2</td><td>W</td><td>95</td><td>1.26 ± 0.04</td><td>2.09 ± 0.06</td><td>0.60 ± 0.09</td><td>16 ± 4</td><td>26 ± 12</td></tr><tr><td>3</td><td>W</td><td>140</td><td>4.12 ± 0.01</td><td>2.09 ± 0.10</td><td>1.97 ± 0.31</td><td></td><td></td></tr><tr><td>1 and 2</td><td>E</td><td>95</td><td>0.04 ± 0.02</td><td>1.3</td><td>0.03 ± 0.02</td><td>3.0 ± 0.4</td><td>3.5 ± 0.7</td></tr><tr><td>3</td><td>E</td><td>140</td><td>0.24 ± 0.01</td><td>1.32 ± 0.06</td><td>0.18 ± 0.03</td><td></td><td></td></tr></tbody></table>

**Note.** $\nu_{\text{obs}}$ corresponds to the observed frequency of the continuum. $S_{\text{obs}}$ is the observed flux density from the continuum image. $\mu$ is the flux weighted magnification of the continuum obtained from lens modeling. Note that, due to low signal-to-noise ratio at 95 GHz in E, we take magnification from Marrone et al. (2018). $S_{\text{int}}$ is intrinsic flux density ($S_{\text{obs}}/\mu$) with 15% additional uncertainty to account for uncertainties in the absolute flux calibration and lens modeling, added in quadrature. $L_{\text{FIR int}}$ is the intrinsic $L_{\text{FIR}}$ calculated from the spectral energy distribution (SED) for each source (flux integrated from 42.5–122.5 µm in the rest frame). $L_{\text{IR int}}$ is the intrinsic $L_{\text{IR}}$ (flux integrated from 8–1000 µm in the rest frame).

properties are present in Strandet et al. (2017). Since CO(3–2) is spatially unresolved (~5"–6") into W and E, we distribute the flux density of CO(3–2) by scaling to CO(6–5) in the two galaxies.

### 3. Data Analysis

#### 3.1. Spectral Line Decomposition

The spectra in W and E are extracted from the 100 km s⁻¹ spectral cubes using apertures of 2".5 and 1".2, respectively. The continuum subtracted spectrum is shown in Figure 2. Since the spectral lines are non-Gaussian and blended in the case of CO(7–6) and [C I](2–1) and in the case of CO(10–9) and p-H₂O(3,1,2 – 2,2,1), we use [C II] from the ALMA 240 GHz observations (Marrone et al. 2018) as a template to derive the spectral properties. [C II] is observed at a high spatial and spectral resolution with a signal-to-noise ratio of ~56 in W and ~23 in E. We reimaged the [C II] data using the same weighting and velocity averaging as described in Section 2.1. The [C II] spectrum obtained has a FWHM of 779 ± 25 km s⁻¹ in W and 371 ± 12 km s⁻¹ in E. This [C II] spectrum is scaled to the observed spectrum of the lines using a least squares fit. We adopt the standard deviation of the flux densities in line-free channels as the error in each velocity bin. The velocity-integrated observed line flux, $I_{\text{obs}}$ (Jy km s⁻¹), is obtained by adding the flux density under the scaled [C II] spectrum from line center –3$\sigma_{\nu}$ to line center +3$\sigma_{\nu}$ km s⁻¹ where $\sigma_{\nu}$ is FWHM/2.35. We estimate the line properties by repeating the scaling 3000 times with random Gaussian noise added to the flux density and taking the median value of all the fits. The uncertainty of the values is taken as the standard deviation of all the fits. In the case of blended lines, we shift the [C II] spectrum to the centers of the two lines and perform a joint fit.

Using the velocity-integrated line flux densities, we calculate the line luminosities using the equation from Solomon et al. (1997) as given below:

$$ 
\begin{aligned}
L_{\text{line}} &= (1.04 \times 10^{-3}) I_{\text{obs}} \nu_{\text{rest}} D_{L}^{2} (1+z)^{-1} \\
L'_{\text{line}} &= (3.25 \times 10^{7}) I_{\text{obs}} D_{L}^{2} (1+z)^{-3} \nu_{\text{obs}}^{-2}
\end{aligned}
\quad (1) 
$$

where $L_{\text{line}}$ is the total observed line luminosity in units of $L_{\odot}$, $L'_{\text{line}}$ is the luminosity in units of K km s⁻¹ pc⁻², $\nu_{\text{rest}}$ and $\nu_{\text{obs}}$ are the rest and observed frequencies of the line in GHz, and $D_L$ is the luminosity distance to the source at a redshift $z$, in Mpc. The line properties are given in Table 3.

<table><caption>**Table 3**
Line Properties</caption><thead><tr><th>Line</th><th>Source</th><th>I<sub>obs</sub><br>(Jy km s<sup>-1</sup>)</th><th>L<sub>obs</sub><br>(×10<sup>8</sup> L<sub>☉</sub>)</th><th>μ</th></tr></thead><tbody><tr><td rowspan="2">CO(6–5)</td><td>W</td><td>2.33 ± 0.09</td><td>10.22 ± 0.39</td><td>2.08 ± 0.21</td></tr><tr><td>E</td><td>0.18 ± 0.03</td><td>0.79 ± 0.13</td><td>1.3</td></tr><tr><td rowspan="2">CO(7–6)</td><td>W</td><td>2.15 ± 0.12</td><td>11.00 ± 0.61</td><td>2.13 ± 0.14</td></tr><tr><td>E</td><td>0.17 ± 0.04</td><td>0.87 ± 0.20</td><td>1.3</td></tr><tr><td rowspan="2">CO(10–9)</td><td>W</td><td>1.41 ± 0.10</td><td>10.30 ± 0.73</td><td>2.16 ± 0.11</td></tr><tr><td>E</td><td>0.14 ± 0.04</td><td>1.02 ± 0.29</td><td>1.3</td></tr><tr><td rowspan="2">[C I](2–1)</td><td>W</td><td>1.42 ± 0.12</td><td>7.28 ± 0.61</td><td>2.13 ± 0.14</td></tr><tr><td>E</td><td>0.04 ± 0.04</td><td>0.21 ± 0.21</td><td>1.3</td></tr><tr><td rowspan="2">p-H<sub>2</sub>O(3,1,2 – 2,2,1)</td><td>W</td><td>0.22 ± 0.10</td><td>1.63 ± 0.74</td><td>2.2</td></tr><tr><td>E</td><td>0.03 ± 0.03</td><td>0.23 ± 0.26</td><td>1.3</td></tr><tr><td rowspan="2">p-H<sub>2</sub>O(2,1,1 – 2,0,2)</td><td>W</td><td>0.97 ± 0.09</td><td>4.63 ± 0.43</td><td>2.2</td></tr><tr><td>E</td><td>0.01 ± 0.04</td><td>0.05 ± 0.19</td><td>1.3</td></tr></tbody></table>

**Note.** $I_{\text{obs}}$ is the observed (not corrected for magnification) integrated flux density in Jy km s⁻¹ (Section 3.1). $L_{\text{obs}}$ is the observed line luminosity in $L_{\odot}$, $\mu$ is the magnification of the lines. In W, the CO magnification is obtained from lens modeling using velocity-integrated measurement sets (Section 3.2). For the H₂O lines, we adopt $\mu = 2.2$ in W. For all the lines in E, we adopt $\mu = 1.3$ (Marrone et al. 2018). p-H₂O(3,1,2 – 2,2,1) is not detected in W and E. [C I](2–1) and p-H₂O(2,1,1 – 2,0,2) are not detected in E. The integrated flux densities in these lines are obtained through the template fitting procedure described in Section 3.1.

#### 3.2. Lens Modeling

SPT0311-58 is a gravitationally lensed system. To infer the intrinsic properties of the source, a lensing reconstruction has to be performed. Lens modeling using a pixellated lensing reconstruction technique (Hezaveh et al. 2016) on high-resolution (0".3) rest frame of 160, 110, and 90 µm continuum and [C II] gives a magnification of $\mu_E = 1.3$, $\mu_W = 2.2$, and $\mu_{tot} = 2.0$ (Marrone et al. 2018). In this analysis, we perform the lens modeling on the 140 and 101 GHz (~272 and 375 µm rest frame, respectively) continuum and the CO(6–5), CO(7–6), and CO(10–9) molecular line transitions in W. Since the signal-to-noise ratio is not sufficiently high for pixellated modeling, we use a parametric lens modeling code, visilens (Spilker et al. 2016). In this code, the background source is parameterized by one or more Sérsic (Sersic 1968) profiles with seven free parameters: the source position relative to the lens ($x_S$, $y_S$), the flux density ($S$), the half light radius major axis ($a_S$), the Sérsic index ($n_S$), the axis ratio ($b_S/a_S$), and the position angle ($\phi_S$). The source profile can be constrained to be circular with $n_S = 0.5$ (i.e., Gaussian profile). The lensing source is represented by one or more Singular Isothermal Ellipsoids (SIEs), and it has five free parameters: the center of the
lens relative to phase center ($x_L, y_L$), the lens mass ($M_L$) which determines the strength of lensing, the ellipticity of the lens ($e_L$), and the position angle of the major axis ($\phi_L$) of the lens in degrees east of north. This model uses a Markov Chain Monte Carlo (MCMC) algorithm (emcee package; Foreman-Mackey et al. 2013) to sample the parameter space. The analysis is performed in the visibility plane instead of the image plane to account for the correlated noise in the interferometric images.

For lens modeling in this analysis, we average the data to decrease the number of visibilities. The performance of different models is measured by Deviance Information Criterion (DIC; Spiegelhalter et al. 2002) to determine the number of sources to include in the model. As the best-fit parameters, we take the median value of the 1000 MCMC chains with $1\sigma$ uncertainty. For all the lensing reconstruction in this paper, we vary the lens parameters only within the best-fit values obtained from the high resolution 140 GHz measurement set using a single lens. The lens parameters obtained from the 140 GHz data agree with Marrone et al. (2018) model. For the lines, we fit single-channel models with the width of 1000 km s$^{-1}$, which is equivalent to the ~2 × FWHM to estimate the overall magnification. The best-fit parameters for continuum and lines in W are given in the Appendix (Tables 8, A.2.2, and A.2.3). Since we do not have a good signal-to-noise for lensing reconstruction in E, we adopt a magnification of 1.3 for all the molecular lines. For the same reason, we adopt $\mu = 2.2$ from Marrone et al. (2018) in p-H$_2$O(2$_{1,1}$ − 2$_{0,2}$) in W. The magnification values used in this analysis for continuum and molecular lines are given in Tables 2 and 3, respectively.

### 3.3. SED Fitting and Dust Mass

To estimate $L_{IR}$ and $L_{FIR}$ in the two sources, we fit a dust SED to the magnification corrected photometry given in the A.1 (Table A.1.1). A 15% absolute calibration error is added in quadrature to the statistical error to account for the uncertainty from the absolute flux calibration and the lens modeling. Here, we fit a simplistic SED model to the dust, and in Section 4.4 we perform a joint fit of dust, CO, and [C I], based on radiative transfer modeling (see Figures 6 and 7). In this section, we assume a uniform dust temperature ($T_{dust}$) in each source and fit a modified blackbody function to the photometry greater than 50 $\mu$m rest frame, and the equations are adopted from Weiß et al. (2007):

$$S_{\nu} = \frac{\pi r^2}{D_A^2} \frac{1}{(1+z)^3} \times [B_{\nu}(T_{\text{dust}}) - B_{\nu}(T_{\text{CMB}})](1 - e^{-\tau_{\nu}}) \quad (2)$$

where $B_{\nu}(T_{\text{dust}})$ is the Planck function at dust temperature, $D_A$ is the angular diameter distance, and r is the effective radius of the source. The optical depth $\tau_{\nu}$ is given by

$$\tau_{\nu} = \frac{\kappa(\nu) M_{\text{dust}}}{\pi r^2} \quad (3)$$

where $M_{\text{dust}}$ is the dust mass. The frequency dependent dust absorption coefficient is adopted from Draine (2003):

$$\kappa_{\nu} = 0.038(\nu/372.7\text{GHz})^{\beta_{T_{d}}} [\text{m}^2 \text{kg}^{-1}]. \quad (4)$$

Here, $\beta_{T_d}$ is the spectral index that determines the slope of the Rayleigh–Jeans tail of the blackbody. We fix $\beta_{T_d}$ to 2.0 (Greve et al. 2012).

In the case of SPT0311-58, we have an estimate of the effective radius of the sources from the 95 and 140 GHz lens models. We use an MCMC algorithm to sample the dust mass

and the dust temperature by giving an upper limit on the effective radius as 5000 pc in W and 1300 pc in E.

From the SED modeling in SPT0311-58, we get $T_{\text{dust}} = 69 \pm 20$ K in W and $49 \pm 9$ K in E, and an intrinsic $M_{\text{dust}} = (1.3 \pm 0.5) \times 10^9 M_{\odot}$ in W and $(1.8 \pm 0.6) \times 10^8 M_{\odot}$ in E. We estimate intrinsic $L_{\text{IR}} = (26 \pm 12) \times 10^{12} L_{\odot}$ and $L_{\text{FIR}} = (16 \pm 4) \times 10^{12} L_{\odot}$ in W, and $L_{\text{IR}} = (35 \pm 7) \times 10^{11} L_{\odot}$ and $L_{\text{FIR}} = (27 \pm 4) \times 10^{11} L_{\odot}$ in E. These values agree with the luminosities inferred from CIGALE SED fits from Marrone et al. (2018) within the uncertainties.

For comparison with SPT0311-58, we have also fitted the SEDs to the SPT-SMG sample and the other literature sources at low and high redshift. In the case of the literature sources, we do not constrain the radius in the SED fit. Instead, we adopt the $\lambda_0 - T_{\text{dust}}$ relationship from Spilker et al. (2016), where $\lambda_0$ is the wavelength at which the optical depth is unity. We then estimate the radius using $\sqrt{\kappa(\nu_0)M_{\text{dust}}}$, where $\kappa(\nu_0)$ is the dust absorption coefficient at $\nu_0 = c/\lambda_0$. The MCMC sampling of the dust mass and dust temperature is the same as those that followed for SPT0311-58.

## 4. Results

### 4.1. $L_{H_2O} - L_{FIR}$ Correlation

Water is the third most abundant molecule in the universe (Neufeld et al. 1995), and its line intensity can be as bright as mid-J CO transitions in high-redshift (U)LIRGs (e.g., Omont et al. 2013; Yang et al. 2013). Multiple H$_2$O transitions from low-redshift galaxies are detected by Herschel Space Observatory while ground-based facilities, such as Caltech Submillimeter Observatory (CSO), Northern Extended Millimeter Array (NOEMA), and the ALMA have detected redshifted H$_2$O transitions from high-redshift galaxies. From these observations, it has been shown that H$_2$O traces $L_{\text{IR}}$, both in low-redshift galaxies and high-redshift (U)LIRGs (Omont et al. 2013; Yang et al. 2013, 2016; Jarugula et al. 2019). This can be explained by the excitation mechanism of H$_2$O molecules where infrared pumping plays a major role (e.g., González-Alfonso et al. 2010, 2012): the higher transitions, such as p-H$_2$O(2$_{1,1}$ − 2$_{0,2}$), p-H$_2$O(2$_{2,0}$ − 1$_{1,1}$), and p-H$_2$O(2$_{2,0}$ − 2$_{1,1}$) are pumped by the 101 $\mu$m FIR photons from the base 1$_{1,1}$ level. These lines are mainly found in the warm, star-forming regions of galaxies.

We detect p-H$_2$O(2$_{1,1}$ − 2$_{0,2}$) in W with a signal-to-noise ratio ~4 above the continuum. In the literature, this is the most distant detection of H$_2$O in a galaxy without an AGN. H$_2$O has been detected in the host galaxy of a quasar candidate at redshift 7 (Koptelova & Hwang 2019). Since the p-H$_2$O(2$_{1,1}$ − 2$_{0,2}$) is not significantly detected in E, we give a flux density based on the template fit (see Section 3.1 for the template fitting procedure). We compare the H$_2$O luminosity ($L_{H_2O}$) in the two galaxies of SPT0311-58 with the local-redshift galaxies and other high-redshift galaxies as shown in Figure 3. For all the galaxies in the plot, we derive $L_{\text{FIR}}$ by fitting a modified blackbody function to the photometry. In the case of the local galaxies, the continuum flux density values are from Sanders et al. (2003), and $L_{H_2O}$ is from Yang et al. (2013). The photometry and the magnification in the high-redshift galaxies are from Weiß et al. (2007), Riechers et al. (2009), Riechers et al. (2013), Bussmann et al. (2013), and Cooray et al. (2014), and $L_{H_2O}$ is from Omont et al. (2013) and Yang et al. (2016). We additionally include another SPT source, SPT0346-52 from Apostolovski et al. (2019), in the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** SPT0311-58 W is shown in red, and a 1σ upper limit in E is shown in blue. The H₂O observations in local (U)LIRGs are from Yang et al. (2013), and the high-redshift ULIRGs are from Omont et al. (2013) and Yang et al. (2016). We also include SPT0346-52 observations from Apostolovski et al. (2019) shown in magenta. The dashed black line shows the slightly super-linear correlation. The fit by fixing the slope to 0.0 is shown as the thick black line, and the error is the gray shaded region.</figcaption></figure>

plot. The p-H₂O(2,1,1-2,0,2) luminosity in SPT0311-58 W and E is consistent within the scatter of the other high-redshift ULIRGs.

To get the correlation between L<sub>H₂O</sub> and L<sub>FIR</sub>, we perform two MCMC linear regressions to all the sources in log-log space as shown in Figure 3. We assume that there is no differential magnification between H₂O and the continuum emission in the gravitationally lensed sources. In Figure 3, one of the fits fixes the slope of log₁₀(L<sub>H₂O</sub>/L<sub>FIR</sub>) versus log₁₀(L<sub>FIR</sub>) to zero (thick black line), and the other allows the slope to vary (dashed line). From the second fit, we see that the relationship is slightly super-linear, which is discussed in the literature (e.g., Omont et al. 2013; Yang et al. 2016). One explanation could be the increase in the optical depth at 100 µm with increasing L<sub>FIR</sub>, which further enhances the p-H₂O(2,1,1-2,0,2) emission. A similar super-linear correlation is observed in another H₂O transition, p-H₂O(2,0,2-1,1,1), which is also excited by 100 µm photons (e.g., Yang et al. 2016; Jarugula et al. 2019). The fit with the slope fixed to zero gives:

$$ \frac{L_{\text{H}_2\text{O}}}{L_{\text{FIR}}} = (1.10 \pm 0.4) \times 10^{-5} \qquad (5) $$

The star formation rate (SFR) is traditionally calculated from L<sub>FIR</sub>, which is a good tracer of the star formation under the assumption that young stars are dust obscured. This assumption holds true in the (U)LIRGs. However, to estimate L<sub>FIR</sub>, the peak of the SED at λ<sub>rest</sub> = 100 µm has to be well-sampled, which is observationally expensive. The long-wavelength spectral lines, such as H₂O, which is bright and well-correlated with L<sub>FIR</sub>, can be used as an alternative tracer of star formation. We use the L<sub>H₂O</sub> - L<sub>FIR</sub> correlation to estimate the SFR in W and E. The SFR scaling relation from L<sub>IR</sub> is taken from Kennicutt & Evans (2012):

$$ \text{SFR}[M_{\odot}/\text{yr}] = 1.47 \times 10^{-10} L_{\text{IR}}[L_{\odot}]. \qquad (6) $$

We calibrate this SFR with L<sub>H₂O</sub> using L<sub>IR</sub>/L<sub>FIR</sub> values from the SED fit and Equation (5). We use the following equation to

estimate the SFR from p-H₂O(2,1,1-2,0,2):

$$ \text{SFR}[M_{\odot}/\text{yr}] = 2.07 \times 10^{-5} L_{\text{H}_2\text{O}}[L_{\odot}]. \qquad (7) $$

We estimate an intrinsic SFR of (4356 ± 2143) M<sub>☉</sub>/yr in W and a 1σ upper limit of 385 M<sub>☉</sub>/yr in E.

## 4.2. Gas Mass from [C I]

The molecular gas mass is traditionally estimated from CO(1–0) emission line luminosity by assuming a CO-to-H₂ conversion factor (α<sub>CO</sub>; e.g., Bolatto et al. 2013). However, α<sub>CO</sub> can vary significantly depending on the physical environments of the galaxies, such as the gas density, the temperature, and the starburst phase driven by mergers (e.g., Maloney & Black 1988; Narayanan 2011; Bolatto et al. 2013). Moreover, it is challenging to observe CO(1–0) from high-redshift galaxies where the CMB temperature is high and can dominate the CO signal. The dependence of the CO brightness temperature on the CMB temperature is discussed in the literature (e.g., da Cunha et al. 2013; Tunnard & Greve 2016).

An alternative molecular gas tracer is [C I](1–0) (492 GHz), whose luminosity is observed to be linearly correlated with CO(1–0) luminosity across a wide range of environments both in local- and high-redshift galaxies (e.g., Papadopoulos et al. 2004; Alaghband-Zadeh et al. 2013; Jiao et al. 2017; Valentino et al. 2020). In our current analysis, we have observations of [C I](2–1) (809 GHz), which is brighter than [C I](1–0) and comes for free with CO(7–6). We estimate the gas mass from [C I](2–1) using the equation from Weiß et al. (2003).<fnref n="22" /> To consider the effect of the background CMB temperature at redshift 6.9, which is ~21.5 K, we include a factor B<sub>ν</sub>(T<sub>ex</sub>)/(B<sub>ν</sub>(T<sub>ex</sub>) - B<sub>ν</sub>(T<sub>CMB</sub>)) to the equation.

$$ M_{\text{C I}} = 4.566 \times 10^{-4} Q(T_{\text{ex}}) \frac{1}{5} e^{62.5/T_{\text{ex}}} \frac{B_{\nu}(T_{\text{ex}})}{(B_{\nu}(T_{\text{ex}}) - B_{\nu}(T_{\text{CMB}}))} L'_{\text{C I}}(2-1) [M_{\odot}] \qquad (8) $$

Here, Q(T<sub>ex</sub>) is the partition function, which is given by:

$$ Q(T_{\text{ex}}) = 1 + 3e^{-T_{\text{I}}/T_{\text{ex}}} + 5e^{-T_{\text{2}}/T_{\text{ex}}} \qquad (9) $$

where T<sub>I</sub> = 23.6 K and T<sub>2</sub> = 62.5 K are the excitation energy transitions of the two transitions of [C I]. The excitation temperature, T<sub>ex</sub>, can be estimated from the [C I] line ratios (Stutzki et al. 1997). However, since we only have observations of [C I](2–1), we adopt a typical value of T<sub>ex</sub> = 30 K (Walter et al. 2011).

The gas mass from [C I](2–1) is converted into the total H₂ gas mass by assuming [C I]/[H₂] abundance of (8.4 ± 3.5) × 10⁻⁵ (Walter et al. 2011). To include the contribution from Helium, we multiply the H₂ gas mass by 1.36. We estimate an intrinsic gas mass of (24.9 ± 10.7) × 10¹⁰ M<sub>☉</sub> in W and (11.4 ± 11.4) × 10⁹ M<sub>☉</sub> in E. Note that [C I](2–1) is not detected in E, and we present the gas mass derived from the flux density obtained by the template fitting (see Section 3.1 for the template fitting procedure).

<footnote n="22">²² Note that the coefficient in Weiß et al. (2003) should be 4.556 × 10⁻⁴. This correction is mentioned in the footnote in Weiß et al. (2005a).</footnote>
<figure><img src="image_6.png" /></figure>

**Figure 4.** The CO SLED normalized by CO(6–5). In the left panel, we compare the two galaxies, W (red) and E (blue), with star-forming galaxies. The local galaxies shown in gray are from Rosenberg et al. (2015). The high-redshift SMGs are from Spilker et al. (2014), including ($z$ = 2–5.7), and Bothwell et al. (2013) is shown as ($z$ = 1.2–4.1). In the right panel, we compare W and E with the representative local- and high-redshift AGNs. Mrk 231 observations are from van der Werf et al. (2010). The high-redshift quasars are Cloverleaf (Bradford et al. 2009; Uzgil et al. 2016) and J2310+1855 (Li et al. 2020). We observe that the SED in W and E peak at transitions similar to high-redshift star-forming galaxies.

### 4.3. CO SLED and AGN Fraction

CO emission lines have traditionally been used as tracers of molecular gas in the local and high-redshift galaxies. CO(1–0) has a critical density of ~$10^3$ cm$^{-3}$ and traces the bulk of the molecular gas in the ISM. However, star formation occurs in the dense molecular clouds, and this warm, dense molecular gas is better traced by mid-$J$ CO lines, such as CO(6–5) and CO(7–6) with a critical density of ~$10^4$ cm$^{-3}$. A linear correlation has been observed between mid-$J$ CO luminosity ($4 \lesssim J_{up} \lesssim 8$) and $L_{IR}$ in local- and high-redshift (U)LIRGs (e.g., Greve et al. 2014; Rosenberg et al. 2015; Lu et al. 2017; Yang et al. 2017). The high-$J$ CO transitions ($J_{up} \gtrsim 9$) require high gas densities ($\ge 10^5$ cm$^{-3}$) and high temperatures. Such favorable conditions for high-$J$ CO excitations can arise in the presence of an AGN (van der Werf et al. 2010), warm photodissociation regions (PDRs), or shocks (Mashian et al. 2015).

The CO SLED shape provides information about the physical conditions of the molecular gas. In Figure 4, we compare the SLED shape of the two galaxies of SPT0311-58 with the local- and high-redshift starburst galaxies in the left panel and the known AGN in the right panel. All the flux densities are normalized to CO(6–5), which is detected at a high signal-to-noise ratio in SPT0311-58 W and E. The local (U)LIRGs sample includes 29 galaxies taken from Rosenberg et al. (2015), excluding NGC 6240 that is an outlier due to galaxy-wide shocks (Meijerink et al. 2013). We also compare SPT0311-58 with the average of 22 SPT-SMGs from $z$ = 2–5.7 (Spilker et al. 2014) and 32 SMGs from $z$ = 1.2–4.1 (Bothwell et al. 2013). The local starburst galaxies peak at $J_{up} \sim 4$ while the high-redshift sample peaks at mid-$J$ CO ($J_{up} \sim 6$) transitions. In both the local- and high-redshift galaxies, there is a drastic decrease in the high-$J$ CO emission at $J_{up} \gtrsim 8$. In the right panel of Figure 4, we see that Mrk231 peaks at $J_{up} \sim 5$ (Weiß et al. 2005a), higher than local starburst galaxies, and the CO emission can be explained by an XDR model detailed in van der Werf et al. (2010). The high-redshift

quasars, Cloverleaf at $z$ = 2.56 (Bradford et al. 2009) and J2310+1855 at $z$ = 6.0 (Li et al. 2020), peak at $J_{up} \sim 9$, and the CO emission is luminous even after $J_{up} \gtrsim 10$. Both the galaxies in SPT0311-58 peak at $J_{up} \sim 6$, similar to other high-redshift starburst galaxies.

We use the ratio of high-$J$ CO to mid-$J$ CO transitions (CO(10–9)/CO(6–5)) to explore the presence of heating mechanisms in addition to the photoelectric heating in SPT0311-58. The ratio of high-$J$ to mid-$J$ CO has been used in literature to characterize the starburst and AGN activity (e.g., Rosenberg et al. 2015; Mashian et al. 2015), which defines the drop-off slope of the CO SLED. In Figure 5, we explore the correlation between $L_{[C II]}/L_{FIR}$, $L_{CO(10-9)}/L_{CO(6-5)}$, and the presence of an AGN. We compare the two galaxies of SPT0311-58 with low-redshift starburst, AGN, and composite galaxies compared to high-redshift SMGs and quasars. The low-redshift [C II] sample is from Díaz-Santos et al. (2014), and the CO observations are from Rosenberg et al. (2015), Mashian et al. (2015). The AGN fraction in the local sample is detailed in Rosenberg et al. (2015). The CO and [C II] observations in the high-redshift galaxies are from Danielson et al. (2011), Frayer et al. (2011), Yang et al. (2017), Andreani et al. (2018), Zhang et al. (2018), Li et al. (2019), Shao et al. (2019), Yang et al. (2019), Wang et al. (2019), and Rybak et al. (2020).

[C II] is a dominant cooling line in the ISM and has been widely studied. It has been observed that $L_{[C II]}/L_{FIR}$ decreases with increasing $L_{FIR}$ (e.g., Malhotra et al. 2001; Brauher et al. 2008; Díaz-Santos et al. 2013; Gullberg et al. 2015; Litke et al. 2019). Several mechanisms have been proposed to explain this [C II] deficit such as positively charged dust grains in the presence of ionizing UV photons, which reduce the photoelectric heating (i.e., UV heating) efficiency, saturation of [C II] emission in the PDR regions, and self absorption (e.g., Narayanan & Krumholz 2014; Muñoz & Oh 2016; Narayanan & Krumholz 2017). If the [C II] deficit is due to the reduced photoelectric heating, this effect will also result in reduced $L_{CO(10-9)}/L_{CO(6-5)}$ since the main collision partner, H$_2$, is
<figure><img src="image_7.png" /><figcaption>Figure 5. The bottom right corner of the plot with enhanced higher-$J$ CO emission and lower [C II] emission, might correspond to heating from X-rays and/or mechanical heating from processes such as AGN activity, stellar winds, or supernovae. We compare the two galaxies, W (red) and E (blue), with starburst, composite, and AGN at local and high redshift. Local galaxies with an AGN fraction >0.3 are indicated by red circles and other Class II or Class III objects (i.e., a high ratio of high-$J$ CO to mid-$J$ CO emission) are shown by gray circles (Rosenberg et al. 2015). The local galaxy observations are detailed in Díaz-Santos et al. (2014), Mashian et al. (2015), and Rosenberg et al. (2015) while the high-redshift sample is from Danielson et al. (2011), Frayer et al. (2011), Yang et al. (2017), Andreani et al. (2018), Zhang et al. (2018), Li et al. (2019), Shao et al. (2019), Yang et al. (2019), Wang et al. (2019), and Rybak et al. (2020). The position of SPT0311-58 E and W on the plot indicates that there might be additional heating mechanisms beyond photoelectric heating.</figcaption></figure>

heated by electrons deep within the molecular clouds in a pure PDR context. However, in Figure 5, we observe that $L_{CO(10-9)}/L_{CO(6-5)}$ is increasing with decreasing $L_{[C II]}/L_{FIR}$. This indicates that other heating mechanisms (e.g., mechanical processes), in addition to the photoelectric effect, result in the observed trend or that the [C II] deficit is probably not due to the effect of reduced photoelectric heating efficiency. The presence of an AGN can also decrease [C II] emission due to possible destruction of dust grains (Smith et al. 2017). There could also be an apparent deficit in [C II]/FIR due to contribution to the total infrared emission from the AGN. However, it has been observed that this effect is not significant (Díaz-Santos et al. 2013). As discussed in Rosenberg et al. (2015), heating mechanisms in addition to photoelectric heating are required to explain the high ratio of high-$J$ to mid-$J$ CO emission.

In Figure 5, reduced [C II] emission and increased CO(10-9) might indicate X-ray heating from the AGN and/or mechanical heating from the AGN or other sources, such as stellar winds, mergers, or supernovae explosions in addition to photoelectric heating. We observe that, in Figure 5, the known AGN and other Class II or III composite galaxies (galaxies with a high or low AGN fraction and a high ratio of high-$J$ CO to mid-$J$ CO emission; Rosenberg et al. 2015) are mostly separated from the starburst galaxies. Both the galaxies of SPT0311-58 have a higher $L_{CO(10-9)}/L_{CO(6-5)}$ comparable to the AGN of Class II or III galaxies. W has a lower $L_{[C II]}/L_{FIR}$ ratio possibly due to an increased star formation density or a higher metallicity than E

(Marrone et al. 2018). We conclude that there is a possibility of other heating mechanisms in addition to photoelectric heating in SPT0311-58, but we cannot infer the presence of an AGN. We further discuss mechanical heating from shocks and stellar feedback in detail in Section 5.4.

## 4.4. Radiative Transfer Modeling

To model the dust emission, observed CO, and [C I] flux densities in SPT0311-58 W and E, we use two non-local thermodynamic equilibrium (non-LTE) LVG radiative transfer modeling methods: N-component and Turbulence. The first model is described in detail in Weiß et al. (2007), and an updated version is in Strandet et al. (2017). The second modeling method, based on the equations from Weiß et al. (2007), is summarized in Harrington et al. (2021).

In the radiative transfer models, we model the excitation of the dust continuum, CO, and [C I] flux densities simultaneously.

We model the dust continuum using the modified blackbody function shown in Equation (2). The dust optical depth (Equation (3)) depends on the dust mass and the dust absorption coefficient. In the radiative transfer models, we adopt $\kappa_{\nu}/\text{m}^2\text{kg}^{-1} = 0.04 \times (\nu/250 \text{ GHz})^{\beta_{T_0}}$. The dust mass is obtained from the gas mass by using the gas-to-dust mass ratio that is a free parameter in the models.

For the CO and [C I] line emission modeling, the infrared radiation is considered as background radiation in addition to the CMB radiation at $z = 6.9$. The line flux density is modeled by the following:

$$S_{\text{CO}/[\text{C I}]} = T_b 2k \nu_{\text{obs}}^2 \Omega/[c^2(1+z)] \quad (10)$$

where c is the speed of light, $\nu_{\text{obs}}$ is the observed frequency of the CO or [C I] lines, and $\Omega$ is the source solid angle given by $\pi r^2/D_A^2$. The brightness temperature of the line, $T_b$, is calculated using the excitation temperature of the line, background temperature, dust temperature, and the dust opacity. These are further dependent on the gas number density, the kinetic temperature, the gas phase abundance per velocity gradient, and the gas-to-dust mass ratio. All these parameters are discussed in detail in this section.

In this section, we also highlight the differences between the two models and present the best-fit and derived parameters obtained from the radiative transfer modeling.

### 4.4.1. N-component and Turbulence Models

Both the N-component and Turbulence codes model the dust, CO, and [C I] flux densities simultaneously. The N-component model is the more basic of these two codes and can describe the ISM using N components. However, due to the increased degeneracies with increased number of free parameters, we use two components in this analysis. These models will henceforth be referred to as 1-component and 2-component, which have 10 and 20 input parameters, respectively (see Table A.3.1 in the Appendix). The Turbulence model realistically models the ISM with 12 input parameters (see Table 4).

One of the main differences between the Turbulence and N-component models is that, in the Turbulence code, the ISM is modeled with the dependence of the source solid angle ($\Omega$) on the gas volume density ($n_{H_2}$, normalized to the mean density).
<table><caption>**Table 4**
Turbulence Model Parameters</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Range</th><th>West</th><th>West (Z<sub>☉</sub>)</th><th>East</th><th>East (Z<sub>☉</sub>)</th></tr></thead><tbody><tr><td colspan="7"><strong>Model input</strong></td></tr><tr><td>log(n<sub>H₂</sub>)</td><td>cm<sup>-3</sup></td><td>1-7</td><td>5.3 ± 1.3</td><td>5.1 ± 0.9</td><td>4.7 ± 2.4</td><td>4.4 ± 2.2</td></tr><tr><td>T<sub>K</sub></td><td>K</td><td>21-600</td><td>116 ± 42</td><td>90 ± 24</td><td>166 ± 65</td><td>136 ± 57</td></tr><tr><td>T<sub>K</sub>/T<sub>dust</sub></td><td>...</td><td>0.5-6.0</td><td>2.0 ± 0.9</td><td>1.6 ± 0.4</td><td>3.4 ± 1.3</td><td>2.7 ± 1.1</td></tr><tr><td>β<sub>T<sub>K</sub></sub></td><td>...</td><td>-0.5-0.05</td><td>-0.1 ± 0.0</td><td>-0.1 ± 0.0</td><td>-0.1 ± 0.0</td><td>-0.1 ± 0.0</td></tr><tr><td>β<sub>T<sub>dust</sub></sub></td><td>...</td><td>1.5-2.0</td><td>1.8 ± 0.1</td><td>1.9 ± 0.1</td><td>2.0 ± 0.0</td><td>2.0 ± 0.0</td></tr><tr><td>R<sub>eff</sub></td><td>pc</td><td>0.1-5000</td><td>3450 ± 750</td><td>3342 ± 628</td><td>2548 ± 954</td><td>2055 ± 870</td></tr><tr><td>κ<sub>vir</sub></td><td>km s<sup>-1</sup> pc<sup>-1</sup> cm<sup>3/2</sup></td><td>1-3</td><td>1.4 ± 0.4</td><td>1.3 ± 0.4</td><td>1.7 ± 0.6</td><td>1.5 ± 0.6</td></tr><tr><td>dv<sub>turb</sub></td><td>km s<sup>-1</sup></td><td>5-200</td><td>134 ± 46</td><td>127 ± 43</td><td>55 ± 61</td><td>71 ± 61</td></tr><tr><td>GDMR</td><td>...</td><td>90-1100</td><td>124 ± 62</td><td>...</td><td>138 ± 95</td><td>...</td></tr><tr><td>GDMR(Z<sub>☉</sub>)</td><td>...</td><td>120-150</td><td>...</td><td>129 ± 9</td><td>...</td><td>130 ± 10</td></tr><tr><td>[CO/H<sub>2</sub>]</td><td>...</td><td>1.0 × 10<sup>-7</sup>-2.0 × 10<sup>-4</sup></td><td>(7 ± 5) × 10<sup>-5</sup></td><td>...</td><td>(7 ± 5) × 10<sup>-5</sup></td><td>...</td></tr><tr><td>[CO/H<sub>2</sub>] (Z<sub>☉</sub>)</td><td>...</td><td>1.0 × 10<sup>-4</sup>-2.0 × 10<sup>-4</sup></td><td>...</td><td>(14 ± 3) × 10<sup>-5</sup></td><td>...</td><td>(14 ± 3) × 10<sup>-5</sup></td></tr><tr><td>[C I/H<sub>2</sub>]</td><td>...</td><td>1.0 × 10<sup>-7</sup>-1.0 × 10<sup>-4</sup></td><td>(3 ± 2) × 10<sup>-5</sup></td><td>(5 ± 2) × 10<sup>-5</sup></td><td>(3 ± 2) × 10<sup>-5</sup></td><td>(4 ± 2) × 10<sup>-5</sup></td></tr><tr><td>β<sub>[C I]</sub></td><td>...</td><td>-5-0</td><td>-0.8 ± 1.1</td><td>-0.7 ± 0.9</td><td>-2.5 ± 1.5</td><td>-2.3 ± 1.6</td></tr><tr><td colspan="7"><strong>Estimated within the model</strong></td></tr><tr><td>T<sub>dust</sub></td><td>K</td><td></td><td>52 ± 5</td><td>55 ± 5</td><td>48 ± 4</td><td>51 ± 4</td></tr><tr><td>M<sub>gas</sub></td><td>M<sub>☉</sub></td><td></td><td>(5.4 ± 3.4) × 10<sup>11</sup></td><td>(4.5 ± 1.8) × 10<sup>11</sup></td><td>(3.1 ± 2.7) × 10<sup>10</sup></td><td>(2.6 ± 0.7) × 10<sup>10</sup></td></tr><tr><td>M<sub>dust</sub></td><td>M<sub>☉</sub></td><td></td><td>(4.3 ± 3.5) × 10<sup>9</sup></td><td>(3.5 ± 1.4) × 10<sup>9</sup></td><td>(2.2 ± 2.5) × 10<sup>8</sup></td><td>(2.0 ± 0.5) × 10<sup>8</sup></td></tr><tr><td colspan="7"><strong>Derived from the model</strong></td></tr><tr><td>L'<sub>CO</sub></td><td>K km s<sup>-1</sup> pc<sup>2</sup></td><td></td><td>(8 ± 1) × 10<sup>10</sup></td><td>(9 ± 1) × 10<sup>10</sup></td><td>(9 ± 2) × 10<sup>9</sup></td><td>(11 ± 2) × 10<sup>9</sup></td></tr><tr><td>α<sub>CO</sub></td><td>M<sub>☉</sub>/K km s<sup>-1</sup> pc<sup>2</sup></td><td></td><td>7.1 ± 5.3</td><td>5.3 ± 2.7</td><td>3.9 ± 4.4</td><td>2.5 ± 0.8</td></tr><tr><td>L<sub>FIR</sub></td><td>L<sub>☉</sub></td><td></td><td>(19 ± 1) × 10<sup>12</sup></td><td>(19 ± 1) × 10<sup>12</sup></td><td>(31 ± 1) × 10<sup>11</sup></td><td>(31 ± 2) × 10<sup>11</sup></td></tr><tr><td>SFR</td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td></td><td>5046 ± 944</td><td>5043 ± 949</td><td>701 ± 151</td><td>708 ± 151</td></tr><tr><td>t<sub>dep</sub></td><td>Myr</td><td></td><td>107 ± 70</td><td>90 ± 40</td><td>44 ± 40</td><td>36 ± 12</td></tr></tbody></table>

**Note.** The input and derived parameters of the *Turbulence* model with units and the explored range. The model outputs are the intrinsic source properties as we use the magnification-corrected photometry and line flux densities for modeling. There are 12 input parameters in the model that are explored within the given range. We present two Turbulence models: the first one is modeled by assuming less than solar metallicity and varying the gas-to-dust mass ratio (GDMR) and the CO abundance; in the second model, to get better constraints on the parameters, we assume solar metallicity in the galaxies ([Z<sub>☉</sub>]) and fix GDMR and CO abundance accordingly. The best-fit parameters from the Turbulence model in SPT0311-58 West and East are shown with and without assuming solar metallicity.

This is given by the following:

$$ \Omega = \Omega_s \times dp(n', dv_{\text{turb}}) \quad (11) $$

where $\Omega_s$ is the solid angle of the total emission region. $dp(n', dv_{\text{turb}})$ is the log-normal probability distribution function (PDF) in supersonically isothermal turbulent gas as given in Krumholz & McKee (2005), which depends on the normalized gas density ($n'$) and the turbulence line width ($dv_{\text{turb}}$). The normalized gas density is given by $n' = n_{H_2}/<n_{H_2}>$ where $<n_{H_2}>$ is the mean gas density. The PDF is given by the following:

$$ dp(n', dv_{\text{turb}}) = \frac{1}{\sqrt{(2\pi\sigma^2)}} \times \frac{1}{n'} e^{-[\log(n')/2\sigma^2]} \quad (12) $$

where, the spread of the PDF is sigma, defined as $\sigma = [\log(1 + \frac{3}{4}\text{Mach}^2)]^{1/2}$ and Mach = $dv_{\text{turb}}/[\sqrt{(k T_K)/2 m_H}]$, where $m_H$ is proton mass, and $T_K$ is the kinetic temperature. For each value of $<n_{H_2}>$, we sample 50 gas densities within $10-10^{10} \text{ cm}^{-3}$. The final dust and gas SED, which describes the observations, is the sum of all 50 SEDs.

In the Turbulence model, the kinetic temperature ($T_K$) is coupled to the molecular gas density as $T_K \propto (n_{H_2})^{\beta_{T_K}}$ where $\beta_{T_K}$ is a negative power-law index. Chemical modeling and simulations have shown that at gas densities $\lesssim 10^5 \text{ cm}^{-3}$, the temperature slightly decreases with increasing densities. This correlation is because the photoelectric and cosmic ray heating has a linear density dependence, but cooling due to CO and [C II] has a super linear relationship (Larson 2005; Meijerink et al. 2007; Krumholz 2014), which leads to an overall cooling per unit of mass in this regime. In addition to $T_K$, the [C I] abundance relative to $H_2$ is also coupled to the gas density as a power law with a negative index, $\beta_{[C I]}$. At gas densities $>10^4 \text{ cm}^{-3}$, the [C I] line ratio, [C I] (2–1)/[C I](1–0), is shown to increase in the PDR models (e.g., Meijerink et al. 2007). However, to reproduce the subthermally excited [C I] gas, i.e., lower [C I](2–1)/[C I](1–0), as typically observed, it is therefore required to reduce the [C I] abundance in the dense gas. Since, we do not have observations of both the [C I] lines, this is of no importance in the current analysis. The N-component code, on the other hand, considers $\Omega$ and $n_{H_2}$ to be independent, and there is no explicit dependence of $T_K$ and [C I] abundance on the gas density.

To constrain the 2-component model, which has more free parameters than the 1-component and Turbulence models, we
consider the maximum value of CO(13–12) from Rosenberg et al. (2015) scaled to CO(10–9) in W and E as upper limits. We also include these upper limits in the 1-component and Turbulence models, although the results do not change significantly when the limits are not included.

The models are optimized to fit the observed flux densities of dust and gas using the Monte Carlo Bees Algorithm (Pham & Castellani 2009). To briefly summarize, the algorithm explores the parameter space within the given ranges by estimating $\chi^2$ for a few different models referred to as bees. Extra bees are assigned to regions of best $\chi^2$ while the rest of the space continues to be explored by the other bees. We run the code 30 times to avoid artificially narrow PDFs for the solutions and hence have 30 final solutions. The best solution is the one with the best $\chi^2$. We evaluate ~10⁷ models, and the best-fit Turbulence model parameter values shown in Table 4 are obtained by taking the mean and standard deviations from all the model runs. The best-fit 1-component and 2-component model parameters are shown in Table A.3.1 in the Appendix.

### 4.4.2. Model Parameters

The free parameters and the corresponding ranges given as inputs in the Turbulence model are shown in Table 4. We use the same ranges for all the parameters in the N-component model, except $\beta_{T_K}$ and $\beta_{[C\,I]}$ since there is no coupling between $T_K$ and [C I] abundance with the gas density in this model. In this analysis, we run two models. The first one is exploring a range of gas-to-dust mass ratio (GDMR) and CO abundance in SPT0311-58 W and E, corresponding to metallicity less than solar metallicity. The second model is by assuming solar metallicity and constraining the GDMR and CO abundance accordingly. In starburst galaxies at high redshift, it is possible for the ISM to be enriched to solar metallicity (Novak et al. 2019; De Breuck et al. 2019). The parameters which are given as inputs to the models are discussed below:

**Gas volume density (log(nH₂))**: We consider a wide range of gas densities, from 10¹–10¹⁰ cm⁻³ in the Turbulence model to sample the density PDF. The mean density of the PDF is sampled from 10¹–10⁷ cm⁻³, and we use the same range in the N-component model.

**Gas kinetic temperature ($T_K$) and Dust temperature ($T_{dust}$)**: The gas kinetic temperature is coupled to the gas volume density in the Turbulence model as mentioned in Section 4.4.1. We explore $T_K$ with the CMB temperature at redshift 6.9 (~21 K) as the lower limit. The dust temperature is weakly coupled to the kinetic temperature through a free parameter, $T_K/T_{dust}$, which allows for additional heating mechanisms. We limit this parameter to 0.5–6 in the models.

**Kinetic temperature power-law index ($\beta_{T_K}$)**: As discussed in Section 4.4.1, the kinetic temperature is expected to decrease with an increase in density below a certain threshold value. From theory and observations of nearby low mass and low density galaxies, Larson (2005) gives the equation of state for this relationship where $\beta_{T_K} = -0.27$. In this analysis of high-redshift galaxies, we explore a wide range between -0.5 and -0.05.

**Dust emissivity ($\beta_{T_d}$)**: We explore the dust spectral emissivity index within a range of 1.5–2.0, which is consistent with the observations from the high-redshift DSFGs (e.g., Conley et al. 2011; Casey et al. 2014).

**Effective radius ($R_{eff}$)**: This radius defines the source solid angle of the emitting region, which normalizes the density PDF

in the Turbulence model. In W and E, we give an upper limit of 5 kpc. This is larger than the size we obtain from the lens modeling because low-J CO is more diffuse than mid- and high-J CO, giving rise to a large radius. CO gas sizes are also found to be larger than the infrared emission (e.g., Spilker et al. 2015; Dong et al. 2019; Apostolovski et al. 2019).

**Virial parameter ($\kappa_{vir}$)**: The velocity gradient, which determines the escape velocity of the gas, is coupled to the gas volume density through the virial parameter (Goldsmith 2001) as

$$dv/dr = 3.1 \kappa_{vir} \sqrt{\frac{n_{H_2}}{10^4}}. \qquad (13)$$

We explore a range of virial parameters from 1–3 where $\kappa_{vir} = 1$ corresponds to virialized gas and >1 corresponds to unbound motions (Greve et al. 2009).

**Turbulence line width ($dv_{turb}$)**: The turbulence line width along the line of sight, in addition to the other free parameters, determines the gas mass including the contribution from Helium (Weiß et al. 2007), which is given by:

$$M_{\text{gas}} = 1.36 \Omega \frac{n_{\text{H}_2}}{dv/dr} dv_{\text{turb}}. \qquad (14)$$

$dv_{turb}/(dv/dr)$ is the equivalent path-length of the molecules. The velocity profile of CO (from the lens modeling in this analysis) and [C II] (Marrone et al. 2018) show a velocity gradient across SPT0311-58 W, which could indicate either rotation of the galaxy or a complicated merger pattern. The observed line width (~1000 km s⁻¹ in all the CO and H₂O lines in W) is a combination of both galaxy rotation and random motions due to turbulence. In this analysis, we explore a range from 5–200 km s⁻¹ for the turbulence velocity line width.

**Gas-to-dust mass ratio (GDMR)**: This parameter is used to calculate the dust mass from the gas mass given in Equation (14). Observations and models suggest that GDMR has a dependence on metallicity where it increases with decreasing metallicity (e.g., Sandstrom et al. 2013; Leroy et al. 2011; Li et al. 2019). We run two models: the first one assumes metallicity less than or similar to solar metallicity and explores the GDMR in the range of 90–1100; in the second model, we use a GDMR in the range of 120–150 assuming the galaxies have enriched Milky Way metallicity (e.g., Draine & Li 2007; Elia et al. 2017). In starburst systems, dense regions can build up metals relatively early and approach solar metallicity (Cen & Ostriker 1999; Novak et al. 2019; De Breuck et al. 2019). We fix the GDMR assuming solar metallicity to better constrain the models. To estimate the dust mass in the LVG models from the gas mass and the GDMR, we adopt $\kappa_{\nu}/\text{m}^2\text{kg}^{-1} = 0.04 \times (\nu/250\,\text{GHz})^{\beta_{T_d}}$.

**CO abundance ([CO/H₂]**): We explore a range of 1.0 × 10⁻⁷–2.0 × 10⁻⁴ for the CO abundance assuming metallicity less than or similar to solar metallicity. Under the assumption of solar metallicity conditions, we use the abundance in the Milky Way and the nearby giant molecular clouds (GMCs) which is in the range of 1.0 × 10⁻⁴–2.0 × 10⁻⁴ (e.g., Blake et al. 1987; Kulesa 2002).

**[C I] abundance ([C I/H₂]**): We explore a range of 1.0–10⁻⁷ – 1.0 × 10⁻⁴. This range includes the [C I] abundance values from the GMCs (e.g., Glover et al. 2015; Fuente et al. 2019).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top: dust continuum, CO, and [C I] SED in W from the 1-component model. The observed data is shown as black points with error bars. The best-fit SED obtained from the model is shown as the thick black line and the 30 other good models in gray. The best-fit model in [C I] is shown as a black circle. Middle: the SED from the 2-component model in W. The contribution from component 1 is shown in blue, and component 2 is in red. Bottom: the SED from the Turbulence model. We present the contribution to the final SED from the ISM at five different gas number densities shown in colored dashed lines. The remaining densities (from the 50 samples) are shown as dashed gray SEDs.</figcaption></figure>

**[C I] abundance power-law index ($\beta_{[C I]}$):** As mentioned in Section 4.4.1, the [C I] abundance decreases with increasing gas density. This is modeled as a power law where $\beta_{[C I]}$ is explored in a wide range of -5 to 0.

Throughout the paper, we refer to the Turbulence, 2-component, and 1-component models as the models run by assuming that SPT0311-58 has a metallicity that is less than the solar metallicity. In the case of models run with an enriched solar metallicity assumption, we specify the assumption explicitly.

### 4.4.3. Model Outputs

The best-fit parameters from the Turbulence model in SPT0311-58 W and E assuming less than solar metallicity and solar metallicity are shown in Table 4. The best-fit value

for each parameter is taken as the mean of outputs from all the 10⁷ models weighted by $\chi^2$, and the error on the parameters is the $\chi^2$ weighted standard deviation. In the Appendix, we show the $\chi^2$ weighted parameter values from all the Turbulence models assuming less than solar metallicity conditions (Figures 16 and 17). The SEDs estimated using the best-fit parameters are shown in Figures 6 and 7 where the top panel corresponds to the dust, CO, and [C I] SEDs obtained from the 1-component model, the middle panel is from the 2-component model, and the bottom panel is from the Turbulence model. In all the models, we show the best-fit SED as the thick black line and the 30 other good $\chi^2$ models in gray. In the 2-component model, we show the relative contribution of component 1 in blue and component 2 in red. As mentioned in Section 4.4.1, in the Turbulence model, the density PDF is sampled at 50
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top: dust continuum, CO, and [C I] SED in E from the 1-component model. Middle: the SED from the 2-component model. Bottom: the SED from the Turbulence model.</figcaption></figure>

densities, which correspond to a source solid angle, and the final SED is the sum of the SEDs obtained at each of the 50 models. In the Turbulence model panel, we show the relative contribution of five densities to the final SED as different dashed colored lines and other densities as dashed gray lines. In this analysis, we focus our discussion on the Turbulence model.

The best-fit mean gas densities are $\log(n_{H_2}/\text{cm}^{-3}) = (5.3 \pm 1.3)$ in W and $(4.7 \pm 2.4)$ in E from the Turbulence model. In the Turbulence model panels of Figures 6 and 7, we see that the larger gas densities ($\log(n_{H_2}) > 5$) emitted from smaller regions of the galaxy do not contribute significantly to the overall dust and CO emission. In the 2-component fit, component 1 has a lower gas density and is emitted from a larger area than component 2 in both W and E (see Table A.3.1 in the Appendix). This indicates that component 1, associated with low-$J$ CO excitations, is primarily tracing diffuse regions of the ISM, while component 2 traces dense and compact gas where

high-$J$ CO is excited. The mean of the gas density from the two components is $\log(n_{H_2}/\text{cm}^{-3}) = (4.7 \pm 1.2)$ in W and $(4.6 \pm 1.3)$ in E. These values are consistent with the Turbulence model within the uncertainty. The kinetic temperature from the Turbulence model is $(116 \pm 42)$ K in W and $(166 \pm 65)$ K in E. This is larger than the dust temperatures, $(52 \pm 5)$ K in W and $(48 \pm 4)$ K in E, since the ratio of $T_K/T_{\text{dust}}$ from the model is greater than 1.0. This could imply that, in addition to photoelectric heating, there are other heating mechanisms such as X-rays, cosmic rays, and mechanical energy input from processes such as AGN outflows, mergers, or stellar feedback. However, it has to be noted that, while $T_K/T_{\text{dust}}$ is a function of visual extinction (Tielens & Hollenbach 1985), we consider a constant value at all extinctions.

The best-fit GDMR from the Turbulence model is $(124 \pm 62)$ in W and $(138 \pm 95)$ in E, which is similar to the enriched Milky Way value. The best-fit CO abundance is $(7 \pm 5) \times 10^{-5}$ in W and E,
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Brightness temperature ratio at each CO excitation transition. The Turbulence model SPT0311-58 values are the mean of the 30 SEDs shown in Figure 6 and Figure 7. The gravitationally lensed SPT-SMG sample is from Spilker et al. (2014), and the unlensed SMGs from $z = 1.2-4.1$ are from Bothwell et al. (2013). A combined sample of starburst galaxies and AGN from $z = 1-7$ is taken from Harrington et al. (2021) and Kirkpatrick et al. (2019). The average Turbulence model output from the 24 Planck sources (Harrington et al. 2021) are also shown.</figcaption></figure>

consistent with the canonical value of $8 \times 10^{-5}$ adopted in the literature (Frerking et al. 1982). The best-fit [C I/H$_2$] is $(3.4 \pm 1.8) \times 10^{-5}$ in W and $(3.4 \pm 2.3) \times 10^{-5}$ in E. These values are consistent with the carbon abundance in dense star-forming environments where a value of $\sim 5 \times 10^{-5}$ is reported in the center of the local starburst galaxy M82 (White et al. 1994). In high-redshift ($z > 2$) samples of submillimeter and quasar host galaxies, Walter et al. (2011) derive a carbon abundance of $(8.4 \pm 3.5) \times 10^{-5}$. Such high values are possible at high redshift where CO molecules are dissociated due to cosmic rays or interstellar radiation, increasing the [C I] abundance (Bisbas et al. 2015). In addition to cosmic rays, other factors are also found to affect the [C I] abundance such as the gas density, the temperature, and the metallicity. In hydrodynamical simulations, the carbon abundance is observed to decrease with an increasing metallicity (Glover & Clark 2016). In the Turbulence model, the [C I] abundance decreases with increasing $n_{H_2}$ using the negative power-law index $\beta_{[C I]}$. This relationship is also manifested in the 2-component model where the [C I] abundance is lower in the denser component 2 than the diffuse component 1 by a factor of $\gtrsim 3$. However, it has to be noted that, due to the lack of [C I](1–0) observations, we cannot reliably confirm the relationship between [C I] abundance and gas density. The [C I] emission is mainly dominated by the diffuse component as seen in the 2-component model (component 1 shown in blue) in Figures 6 and 7.

### 4.4.4. CO and [C I] Excitations

From the radiative transfer modeling, we can estimate the intrinsic brightness temperature ($T_b$) ratios for all the CO transitions from $J = 1-15$ in terms of the line luminosities as following:

$$r_{J_i, J_j} = L'_{CO(J_i \rightarrow J_i-1)} / L'_{CO(J_j \rightarrow J_j-1)} \quad (15)$$

<figure><img src="image_11.png" /><figcaption>**Figure 9.** $L_{[C I](1-0)}/L_{CO(4-3)}$ representing total gas to dense gas ratio against $L_{[C I](1-0)}/L_{FIR}$. The values obtained from the Turbulence model in SPT0311-58 are compared to the observations in local- and high-redshift galaxies. The SPT-SMG data is from Bothwell et al. (2017) and the other lensed SMGs are from Alagband-Zadeh et al. (2013). Local galaxies from Herschel/SPIRE are from Kamenetzky et al. (2016). CO(4–3) and [C I](1–0) line luminosities in SPT0311-58 W and E are obtained from the Turbulence model, where we model the observed dust, CO, and [C I](2–1) transitions. In SPT0311-58 W and E, we show the mean line luminosities and the uncertainty obtained from running the Turbulence model 30 times. The direction of increase in the UV field strength and the density are shown at the top right corner of the plot. The model values in W and E agree with the observations in the SPT-SMG sample.</figcaption></figure>

where $J_j$ is usually 1. In Figure 8 and Table 5 we show the brightness temperature ratios in the two SPT0311-58 galaxies by normalizing to $J_j = 3$, which is the lowest observed CO transition. We compare the observed and model CO excitations in the two galaxies with the stacked $T_b$ from 22 gravitationally lensed SPT-SMGs from Spilker et al. (2014) and with the median values obtained from 32 SMGs from $z = 1.2-4.1$ (Bothwell et al. 2013). We also include a heterogeneous sample of star-forming galaxies and AGN from Kirkpatrick et al. (2019) combined with 24 lensed Planck selected sources from Harrington et al. (2021), which gives a range of galaxies from $z = 1-7$. The average Turbulence model output of the 24 Planck sources from Harrington et al. (2021) is also shown. The Turbulence model reproduces the observed values in SPT0311-58 W and E. Both the SPT0311-58 galaxies have a brightness temperature profile similar to the high-redshift SMG samples, and a subset of sources from Harrington et al. (2021) also show similar excitations.

The non-LTE 2-component best-fit model outputs the excitation temperature ($T_{ex}$) of [C I] in component 1 and component 2. We calculate the total $T_{ex}$ of [C I](1–0) by performing a flux density weighted sum of $T_{ex}$ in component 1 and component 2. This gives $T_{ex} \sim 30$ K and 32 K in W and E, respectively. These values are similar to the typical $T_{ex}$ of 30 K adopted in the literature (Walter et al. 2011). However, due to a lack of [C I](1–0) observations, we cannot get an accurate estimate of the excitation temperatures from modeling alone.

In Figure 9, we compare the Turbulence model [C I](1–0) and CO(4–3) in SPT0311-58 W and E to the sources from the
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparing the gas mass calculated from different methods in SPT0311-58 W and E. Data points from top to bottom: (1) The Turbulence model output from Equation (14) not constraining the GDMR and CO abundance. (2) The Turbulence model output by assuming solar metallicity ([Z<sub>☉</sub>]) and fixing the GDMR and CO abundance. (3) Using the same equation in the 2-component model and adding the gas mass from component 1 and 2. (4) Gas mass from the 2-component model assuming solar metallicity. (5) Gas mass from the 1-component model. (6) Gas mass from the 1-component model assuming solar metallicity. (7) Gas mass from the observed $L_{\text{C1(2-1)}}$ (Equation (8)), and assuming a typical excitation temperature and [C I] abundance. (8) Gas mass estimated from dust mass by assuming GDMR = 100. (9) Gas mass from Marrone et al. (2018) where CO(1–0) is converted from CO(3–2), and assuming a conversion factor $\alpha_{\text{CO}} = 1.0 M_{\odot}/\text{km s}^{-1}\text{pc}^2$.</figcaption></figure>

literature. [C I](1–0) is observed to be a good tracer of CO(1–0), and hence, it is a good tracer of the bulk of the H₂ gas mass (e.g., Jiao et al. 2017, 2019). CO(4–3) traces denser H₂ gas participating in the star formation. Alaghband-Zadeh et al. (2013) observed that the ratio $L_{[C I](1-0)}/L_{\text{CO(4-3)}}$ decreases with increasing $L_{\text{FIR}}$. In Figure 9 we compare the ratio of [C I](1–0) / CO(4–3) in SPT0311-58 with the lensed SMGs (Alaghband-Zadeh et al. 2013; Bothwell et al. 2017) and the local star-forming galaxies from Kamenetzky et al. (2016). The direction in which the UV field strength and the density of the gas increases is also shown in Figure 9 (e.g., Kaufman et al. 1999; Valentino et al. 2020). SPT0311-58 W and E have ratios similar to the ~40 lensed SPT-SMGs sample detailed in Bothwell et al. (2017). W is similar to the mean value of SPT-SMGs while E has a lower value, which might indicate that E has a higher UV radiation and is a more compact starburst than W.

### 4.4.5. Gas Mass and SFR

The gas mass in the LVG model is estimated using Equation (14). From the Turbulence model, we get $M_{\text{gas}} = (5.4 \pm 3.4) \times 10^{11} M_{\odot}$ in W and $(3.1 \pm 2.7) \times 10^{10} M_{\odot}$ in E. Using the GDMR from the Turbulence model, we derive $M_{\text{dust}} = (4.3 \pm 3.5) \times 10^9 M_{\odot}$ in W and $(2.2 \pm 2.5) \times 10^8 M_{\odot}$ in E. By assuming solar metallicity in SPT0311-58, we get $M_{\text{gas}} = (4.5 \pm 1.8) \times 10^{11} M_{\odot}$ in W and $(2.6 \pm 0.7) \times 10^{10} M_{\odot}$ in E from the Turbulence model. Because the GDMR and CO abundance is fixed for solar metallicity, the gas mass is better constrained than in the case where we do not constrain the parameters.

We compare the gas mass calculated using different methods in Figure 10. We compare the gas mass from the LVG models with the calculations from [C I](2–1), discussed in Section 4.2. We also calculate the gas mass from the dust mass (Section 3.3) by assuming the GDMR = 100 (Sandstrom et al. 2013), which includes contribution from Helium. We also show the gas

masses reported in Marrone et al. (2018) where they are obtained by scaling CO(3–2) to CO(1–0) and converting to gas mass by assuming $\alpha_{\text{CO}} = 1.0 M_{\odot}/\text{km s}^{-1}\text{pc}^2$. The total gas mass from the 2-component model is consistent with the Turbulence models within the uncertainties. The gas mass estimated from [C I](2–1) agrees with the LVG models within the uncertainties. The gas mass from [C I](2–1) in E is the upper limit due to non-detection of the line. SPT0311-58 W and E gas mass estimates reported in Marrone et al. (2018) and the estimates from dust mass include assumptions about the CO scaling, gas mass conversion factor, and the GDMR which could result in discrepancy with the LVG models.

All the methods used to estimate the gas mass yield a canonical value of a CO-to-H₂ conversion factor, $\alpha_{\text{CO}} > 0.8 M_{\odot}/\text{km s}^{-1}\text{pc}^2$ as discussed in Section 5.2. The intrinsic SFR from the Turbulence model is calculated using Equation (6), giving an SFR = $5046 \pm 944 M_{\odot}$ in W and $701 \pm 151 M_{\odot}$ in E.

## 5. Discussion

In this section we discuss the differential magnification across the different CO transitions and dust in SPT0311-58. We also compare the spatial extent of CO and dust. We further estimate the CO-to-H₂ conversion factor and gas depletion timescales in SPT0311-58 and compare them with the other high-redshift SMGs in the literature. Toward the end of this section, we briefly discuss some of the heating mechanisms contributing to the CO emission in SPT0311-58.

### 5.1. Magnification and Intrinsic Size

In Figure 11 we compare the magnification and intrinsic size of CO and the dust continuum regions in SPT0311-58 W obtained from the lens models. The magnification for the dust continuum and CO are given in Tables 2 and 3 respectively, and the lens model parameters are given in Tables A.2.2 and A.2.3.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Top: the magnification of CO and the dust emission in SPT0311-58 W. There is no significant differential magnification between the different CO transitions and the dust. Bottom: CO and the dust emission region sizes (circularized radius) in SPT0311-58 W. The higher-$J$ CO emission arises from more compact regions of the galaxy than the lower-$J$ CO transitions. CO has been observed to have a marginally larger size than the dust.</figcaption></figure>

The circularized intrinsic radius is calculated from the semimajor axis ($a_s$) and the semiminor axis ($b_s$) as $\sqrt{(a_s b_s)}$. Differential magnification can occur depending on the position of the source relative to the lensing caustic and the size of the emission region (Hezaveh et al. 2012; Spilker et al. 2015). Significant differential magnification between the CO lines, particularly between those tracing different physical conditions of the galaxy, can affect the physical properties derived from the CO SLED (Dong et al. 2019). From the top panel in Figure 11, we see that the magnification is consistent between the CO transitions and the dust continuum (i.e., we do not observe differential magnification between these components in SPT0311-58 W). From the bottom panel of the figure, we see that the CO emission area is marginally decreasing with an increasing excitation level. This trend has been observed in sources such as M82 (Weiß et al. 2005b) and SPT0346-52 (Apostolovski et al. 2019) where low-$J$ CO transitions, such as CO(1–0) and CO(2–1), are tracing the diffuse regions of the

galaxy while the high-$J$ CO is emitted from the compact regions. We also observe that the CO emission regions are larger than the dust emission, although not significantly. Our result is consistent with the previous literature where CO is observed to have a larger radial extent than the dust (e.g., Spilker et al. 2015; Apostolovski et al. 2019; Dong et al. 2019). In Calistro Rivera et al. (2018) and Dong et al. (2019), the authors discuss several possibilities for such a trend (for example, a temperature gradient across the galaxy giving rise to a lower intensity of continuum emission in the outskirts or a spatially varying gas-to-dust ratio). The compact size of the dust continuum in comparison to the cold gas due to a temperature gradient across the galaxy is also discussed in Cochrane et al. (2019). In SPT0346-52 (Apostolovski et al. 2019), it is observed that the radial extent of the dust at 3 mm is larger than the radial extent at 2 mm. In SPT0311-58, the size of the emission region of the dust at 2 mm (ALMA B4) is slightly higher than the size at 3 mm (ALMA B3), but this difference is not statistically significant (bottom panel in Figure 11).

## 5.2. Gas Mass Conversion Factor

Traditionally, the gas mass is estimated by assuming a CO-to-H$_2$ conversion factor ($\alpha_{CO}$) where $M_{gas} = \alpha_{CO} L'_{CO(1-0)}$. In the literature, $\alpha_{CO} = 0.8 M_{\odot}/\text{K km s}^{-1}\text{pc}^2$ is typically adopted for (U)LIRGs (e.g., Downes & Solomon 1998; Carilli & Walter 2013). We hereby mention $\alpha_{CO}$ values without units for convenience. There is a large uncertainty in $\alpha_{CO}$ values ranging from 0.4–6 (e.g., Ivison et al. 2011; Papadopoulos et al. 2012; Carilli & Walter 2013; Mashian et al. 2013). In SPT0311-58, we calculate $\alpha_{CO}$ using three different estimates of gas mass as shown in Figure 12. The CO(1–0) luminosity is obtained from the Turbulence model where $L'_{CO(1-0)} = (8.3 \pm 1.3) \times 10^{10} \text{ K km s}^{-1}\text{pc}^2$ in W and $(8.8 \pm 2.6) \times 10^{10} \text{ K km s}^{-1}\text{pc}^2$ in E. In the first method, we derive the conversion factor based on the intrinsic gas mass (Equation (14)) from the Turbulence model. Assuming less than solar metallicity in SPT0311-58, we obtain $\alpha_{CO} = 7.1 \pm 5.3$ in W and $3.9 \pm 4.4$ in E. Assuming the solar metallicity GDMR and CO abundance, we obtain $\alpha_{CO} = 5.3 \pm 2.7$ in W and $2.5 \pm 0.8$ in E. Since we fix two parameters under solar metallicity assumptions, the $\alpha_{CO}$ values are better constrained. This is similar to $\alpha_{CO} = 4.8 \pm 2.9$ reported in Strandet et al. (2017) derived from unresolved observations of SPT0311-58 using the 2-component model, assuming solar metallicity, and fixing the CO abundance. In the second method, we estimate the intrinsic gas mass from the dust mass (Section 3.3) by assuming the GDMR = 100. We obtain $\alpha_{CO} = 1.6 \pm 0.7$ in W and $2.1 \pm 0.9$ in E. In the third method, we use gas mass from [C I](2–1) as detailed in Section 4.2, obtaining $\alpha_{CO} = 3.0 \pm 1.4$ in W and an upper limit of 2.8 in E.

In Figure 12, we compare $\alpha_{CO}$ in SPT0311-58 W and E with the literature sample. In all the literature sources, we estimate the conversion factor using the second method, where we calculate the dust mass following the method detailed in Section 3.3. This is done to be consistent and to reduce uncertainties from different assumptions in the dust mass calculations, the GDMR values, and the gas mass calculated from dynamical mass estimates. In the left and right panels of Figure 12, we plot $\alpha_{CO}$ as a function of $L_{FIR}$ and the dust temperature ($T_{dust}$), respectively. In the literature sample, both the $L_{FIR}$ and $T_{dust}$ are obtained from the SED fitting procedure as detailed in 3.3. The main sequence (MS) galaxy sample is
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Left: the main sequence (MS) galaxies are taken from Magnelli et al. (2012). The high-redshift SMGs are from Carilli et al. (2010), Fu et al. (2012), Walter et al. (2012), Fu et al. (2013), Ivison et al. (2013), and Alaghband-Zadeh et al. (2013). The SPT-SMGs are from Aravena et al. (2016) and Reuter et al. (2020). The $\alpha_{CO}$ in the literature sample is estimated using $M_{gas}$ from $M_{dust}$ assuming a gas-to-dust mass ratio of 100. The $\alpha_{CO}$ values in SPT0311-58 W and E obtained using $M_{gas}$ from the Turbulence model assuming less than solar metallicity are shown as solid red and blue data points. The Turbulence model $\alpha_{CO}$ values under solar metallicity assumption are shown as open red and blue data points. The $\alpha_{CO}$ using $M_{gas}$ from $M_{dust}$ are shown as solid magenta and cyan data points in W and E, respectively. The $\alpha_{CO}$ using $M_{gas}$ from [C I](2–1) are shown as open magenta and cyan data points. Right: the correlation between $\alpha_{CO}$ and $T_{dust}$ is statistically significant in the MS galaxies and in the SPT-SMGs.</figcaption></figure>

<table><caption>**Table 5**
Brightness Temperature Ratios from the Turbulence Model</caption><thead><tr><th>Source</th><th>r<sub>1,3</sub></th><th>r<sub>2,3</sub></th><th>r<sub>4,3</sub></th><th>r<sub>5,3</sub></th><th>r<sub>6,3</sub></th><th>r<sub>7,3</sub></th><th>r<sub>8,3</sub></th><th>r<sub>9,3</sub></th><th>r<sub>10,3</sub></th></tr></thead><tbody><tr><td>W</td><td>1.04 ± 0.13</td><td>1.09 ± 0.08</td><td>0.87 ± 0.05</td><td>0.73 ± 0.04</td><td>0.59 ± 0.03</td><td>0.45 ± 0.02</td><td>0.33 ± 0.02</td><td>0.22 ± 0.01</td><td>0.13 ± 0.01</td></tr><tr><td>E</td><td>0.88 ± 0.16</td><td>1.05 ± 0.09</td><td>0.89 ± 0.06</td><td>0.76 ± 0.05</td><td>0.63 ± 0.04</td><td>0.50 ± 0.03</td><td>0.37 ± 0.03</td><td>0.26 ± 0.02</td><td>0.17 ± 0.01</td></tr></tbody></table>

**Note.** The brightness temperature ratios are obtained from Equation (15) using the mean of the 30 best SEDs from the Turbulence large velocity gradient (LVG) modeling.

taken from Magnelli et al. (2012). The high-redshift SMGs are from Carilli et al. (2010), Fu et al. (2012), Walter et al. (2012), Fu et al. (2013), Ivison et al. (2013), and Alaghband-Zadeh et al. (2013). The SPT-SMG sample is from Aravena et al. (2016) and Reuter et al. (2020). In Alaghband-Zadeh et al. (2013), CO(1–0) is not observed directly and is derived from CO(4–3) using a conversion from Bothwell et al. (2013). In the SPT sources where CO(1–0) is not observed, it is derived from CO(2–1) by assuming a line brightness ratio of 0.9 (Aravena et al. 2016). The errors are large in sources which have photometry data available at three or fewer wavelengths. In the figure, we also show the conversion factors typically adopted in the literature to be 4.6 for the Milky Way and 0.8 for high-redshift (U)LIRGs (Downes & Solomon 1998).

In SPT0311-58, the values from the Turbulence model are closer to the Milky Way value. The conversion factor estimated from the dust mass in the SPT sample is larger than the (U) LIRGs value of 0.8. This discrepancy in $\alpha_{CO}$ is discussed in Ivison et al. (2011) where they find that a range of $\alpha_{CO} \sim 5–10$ and 0.4–1 can both explain the gas excitations in the SMGs. Constraining $\alpha_{CO}$ in the (U)LIRGs based solely on low-$J$ CO, which traces an extended low density and a warm diffuse gas, results in a lower $\alpha_{CO}$ as it does not account for all the gas mass. A higher density and a lower $T_{kin}$ turbulent gas in the (U)LIRGs can dominate most of the gas mass and can increase $\alpha_{CO}$ to galactic values (Ivison et al. 2011; Papadopoulos et al. 2012;

Scoville et al. 2012). From our Turbulence model, we find that the dense gas in both E and W ($\log(n_{H_2}/\text{cm}^{-3}) = 4–5$) has a significant contribution to the overall gas emission, which might give galactic values of $\alpha_{CO}$ in SPT0311-58.

In the right panel of Figure 12, we explore the correlation of $\alpha_{CO}$ and $T_{dust}$. We fit a linear function to the MS galaxies and the SPT-SMGs using an MCMC and find a statistically significant correlation between $\alpha_{CO}$ and $T_{dust}$. We use an F-test to determine that the model with a negative slope is statistically better than the model with no slope or zero correlation. We obtain a p-value of <0.05, and we reject the null hypothesis that the complex model (i.e., the fit with negative slope) does not provide more information than the fit with slope fixed to zero. However, we observe no significant correlation if we include all the SMGs in the model. The strong correlation between $\alpha_{CO}$ and $T_{dust}$ in the MS galaxies has been discussed in Magnelli et al. (2012). We find a similar correlation in the SPT-SMG sample where SPT0311-58 W and E (from the dust mass method) follow a similar trend as the other SPT-SMGs. It has to be noted that, while $T_{dust}$ and the dust mass are estimated from the same photometry, the correlation between $M_{dust}$ and $T_{dust}$ is not as steep as the $\alpha_{CO}$ and $T_{dust}$ correlation. This empirical relation between $\alpha_{CO}$ and $T_{dust}$ can be used to select the appropriate $\alpha_{CO}$ value for the gas mass calculations in normal and starburst galaxies. One caveat is that the dust SED fitting, the dust mass calculation, and the gas mass from the
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Depletion time as a function of the redshift. The values in SPT0311-58 are estimated using the two methods outlined in Section 5.3. The depletion time from the Turbulence model, assuming less than solar metallicity, is shown as solid red and blue data points in W and E, respectively. Assuming solar metallicity ($Z_\odot$), the depletion time is better constrained and shown as open red and blue data. The depletion time estimated using the SFR from p-H$_2$O ($2_{1,1} - 2_{0,2}$) and the gas mass from [C I](2–1) are shown as magenta and cyan in W and E. The high-redshift gravitationally lensed SPT-SMGs are from Aravena et al. (2016) and Reuter et al. (2020), and the other high-redshift SMGs are taken from Carilli et al. (2010), Fu et al. (2012), Walter et al. (2012), Fu et al. (2013), Ivison et al. (2013), and Alaghband-Zadeh et al. (2013). The MS galaxies are described in Saintonge et al. (2013) as $t_{\text{dep}} = 1.5(1+z)^\alpha$ where $\alpha$ is from -1.5 (Davé et al. 2012) to -1.0 (Magnelli et al. 2013), which is shown as the gray shaded region. There is no clear evidence of an evolution of the depletion time with the redshift in the SMG sample above $z > 3$.</figcaption></figure>

dust mass have to be estimated using the same assumptions in all the sources to understand this correlation.

From the high-redshift SMG sample together with the SPT-SMGs, we estimate a mean $\alpha_{CO} = 3.2 \pm 2.7$.

### 5.3. Depletion Time

The gas depletion timescale ($t_{\text{dep}}$) is defined as $M_{\text{gas}}/\text{SFR}$. We calculate and compare the depletion timescales in SPT0311-58 W and E with the literature sources in Figure 13. We estimate $t_{\text{dep}}$ in SPT0311-58 using two methods. In the first method, the gas mass and SFR are from the Turbulence model (Section 4.4.5). We estimate $t_{\text{dep}} = 107 \pm 70$ Myr and $44 \pm 40$ Myr in W and E, respectively, assuming less than solar metallicity. Under the assumption of solar metallicity in SPT0311-58, we estimate $t_{\text{dep}} = 90 \pm 40$ Myr and $36 \pm 12$ Myr in W and E, respectively. The second method is using the observed lines where gas mass is calculated from the observed [C I](2–1) (see Section 4.2), and the SFR is calculated from p-H$_2$O($2_{1,1} - 2_{0,2}$) (see Section 4.1). We obtain $t_{\text{dep}} = (57 \pm 37)$ Myr in W and $(143 \pm 570)$ Myr in E. The non-detection of [C I](2–1) and p-H$_2$O($2_{1,1} - 2_{0,2}$) in E gives a large uncertainty on the depletion time estimate. The literature high-redshift SMGs are taken from Carilli et al. (2010), Fu et al. (2012), Walter et al. (2012), Fu et al. (2013), Ivison et al. (2013), and Alaghband-Zadeh et al. (2013). The SPT-SMGs are taken from Aravena et al. (2016) and Reuter et al. (2020). The gas

<table><caption>**Table 6**
Intrinsic Cooling Power and Percent of the Total Cooling Fraction</caption><thead><tr><th rowspan="2">Cooling channel</th><th colspan="2">Power<br>(erg s<sup>-1</sup>)</th><th colspan="2">Fraction relative to [C II]</th><th colspan="2">Total cooling %</th></tr><tr><th>W</th><th>E</th><th>W</th><th>E</th><th>W</th><th>E</th></tr></thead><tbody><tr><td>CO</td><td>~ 1 ×<br>10<sup>43</sup></td><td>~ 2 ×<br>10<sup>42</sup></td><td>0.36</td><td>0.12</td><td>~10%.7</td><td>~3%.7</td></tr><tr><td>[C I]</td><td>~ 2 ×<br>10<sup>42</sup></td><td>~ 8 ×<br>10<sup>40</sup></td><td>0.04</td><td>0.004</td><td>~1%.2</td><td>~0%.1</td></tr><tr><td>[C II]</td><td>~ 4 ×<br>10<sup>43</sup></td><td>~ 2 ×<br>10<sup>43</sup></td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** The intrinsic cooling power of CO is the sum $\sum_{J=1}^{J=15} L_{COJ}$, and [C I] is $\sum_{J=1}^{J=2} L_{[C I]J}$ from the Turbulence model. The [C II] intrinsic power is obtained from Marrone et al. (2018). The percent cooling fraction for CO and [C I] is obtained by considering the fraction of the cooling channel luminosity relative to the [C II] luminosity and assuming that [C II] contributes to ~30% of the total cooling (Rosenberg et al. 2015).

mass in the literature sources is calculated from the CO(1–0) luminosity by assuming the mean $\alpha_{CO} = 3.2 \pm 2.7$ from the high-redshift SMG sample and the SPT-SMG sample (Figure 12). The SFR is calculated using Equation (6) where the $L_{IR}$ is obtained from the modified blackbody SED fit described in Section 3.3. The MS galaxies are described in Saintonge et al. (2013) as $t_{\text{dep}} = 1.5(1+z)^\alpha$ where $\alpha$ is from -1.5 (Davé et al. 2012) to -1.0 (Magnelli et al. 2013), which is shown as the gray shaded region.

From Figure 13, we observe that the depletion time in the SMGs between $z = 2–3$ follows a trend similar to the MS galaxies where $t_{\text{dep}}$ decreases with an increasing redshift (e.g., Saintonge et al. 2013; Tacconi et al. 2013). However, this evolution does not seem to exist in the SMG sample at $z > 3$. In SPT0311-58 W and E, the depletion times obtained from both the methods are within the typical range for high-redshift SMGs.

### 5.4. Energy Budget

In this section, we discuss the heating and cooling budget of the neutral gas in SPT0311-58 with a focus on the dense molecular gas traced by CO. Some of the important cooling lines in the neutral gas regions include [C II] (158 µm), [C I] (369, 609 µm), [O I] (63 µm), [Si II] (35 µm), and CO. The total neutral gas cooling budget is the sum of luminosities of the lines. We combine the observations of [C II] from Marrone et al. (2018) and the total CO and [C I] luminosities from the LVG modeling to estimate the cooling power. However, since we do not have observations of the two important coolants [O I] and [Si II], which are found to contribute to more than 50% of the total cooling budget (Rosenberg et al. 2015), we cannot provide the complete picture of the neutral gas. The total power of CO ($\sum_{J=1}^{J=15} L_{COJ}$) and neutral carbon ($\sum_{J=1}^{J=2} L_{[C I]J}$) from the Turbulence model and [C II] from Marrone et al. (2018) are given in Table 6. We also show the contribution of CO and [C I] to the total cooling by assuming that [C II] contributes to ~30% of the total cooling Rosenberg et al. (2015). In the dense molecular gas regions (i.e., high visual extinction), CO transitions are the dominant cooling lines (Tielens & Hollenbach 1985), where molecular collisions play an important role over photoelectric heating. From the Turbulence
model, we derive a gas density >10⁴ cm⁻³ and a gas kinetic temperature ~100–170 K in SPT0311-58, which is the dense molecular gas regime. Over the depletion time of the molecular gas, i.e., ~100 and 40 Myr in W and E, we estimate the total energy output from CO emission to be ~5 × 10⁵⁸ erg and ~3 × 10⁵⁷ erg, respectively.

The Turbulence LVG model has a free parameter ($T_{\text{kin}}/T_{\text{dust}}$) to account for heating from sources such as X-rays, cosmic-rays, and mechanical heating, in addition to photoelectric heating by modeling $T_{\text{kin}}$ and $T_{\text{dust}}$ simultaneously. $T_{\text{kin}}/T_{\text{dust}} = 2.0 \pm 0.9$ and $3.4 \pm 1.3$ in SPT0311-58 W and E, respectively, which suggests that there are other heating mechanisms in addition to photoelectric heating (traced by $T_{\text{dust}}$). A high value of $T_{\text{kin}}/T_{\text{dust}}$ could also occur in the photodissociation regions at low extinctions or a low density (Tielens & Hollenbach 1985). However, since both galaxies have a high gas density where high-J CO is excited, we can consider heating from other processes. One caveat is that we do not consider the change in $T_{\text{kin}}/T_{\text{dust}}$ with density in the LVG model. Following Harrington et al. (2021), we calculate the contribution of turbulent kinetic energy to the CO line emission in the dense molecular gas. The turbulent energy can be calculated from the turbulence line width and the gas mass (Table 4) as $E_{\text{turb}} = 0.5 M_{\text{gas}} \times d\nu_{\text{turb}}^2$. We get a turbulence energy of ~1 × 10⁵⁹ erg in W and ~1 × 10⁵⁷ erg in E, similar to the total CO cooling energy. This shows that turbulence energy is sufficient to heat the molecular gas in SPT0311-58. There are several possible sources of turbulence-driven mechanical heating such as the presence of AGN outflows, stellar winds, or supernovae explosions.

We also estimate heating from the star formation such as stellar feedback and supernovae explosions based on the calculations from Harrington et al. (2021) and references therein. The SPT sources have a negligible contribution from AGN to the total infrared luminosity, even in the galaxy with one of the highest star formation rate densities (Ma et al. 2016). Hence, the presence of an AGN is probably not a significant source of molecular gas heating in SPT0311-58. It is also argued in Harrington et al. (2021) that the X-ray luminosity from non-AGN sources is not a major source of heating. We estimate heating from the stellar outflows and the supernovae. Stellar feedback in galaxies is primarily contributed by massive stars, such as an O-type, which evolve into core collapse supernovae, neutron stars, or black holes. Using an O-type stellar wind luminosity of ~10⁴⁹⁻⁵¹ erg in a lifetime of 5 Myr (Leitherer et al. 1999; Smith 2014; Ramachandran et al. 2019) and assuming O-type stars comprise 0.2% of the total number of stars formed (Kroupa initial mass function), we can calculate the total energy output over the depletion time of the galaxies. We estimate a total stellar feedback energy of ~10⁵⁸⁻⁶⁰ erg in W and ~7 × 10⁵⁶⁻⁵⁸ erg in E. Stars with masses in the range of ~10–40 $M_{\odot}$ collapse as supernovae (Heger et al. 2003), which emit an energy of ~10⁵¹ erg. Assuming the stars with these stellar masses comprise ~7% of the total stellar mass, we estimate a total energy output from supernovae explosions of ~3 × 10⁶⁰ erg in W and ~2 × 10⁵⁹ erg in E, over the depletion timescale. This is an upper limit as the energy input from the supernovae explosions is not a continuous process. About 1% or less of the supernovae energy goes into turbulent energy (e.g., Iffrig & Hennebelle 2015; Martizzi et al. 2016), which is consistent with the CO cooling energy.

The energy estimates show that the mechanical heating from star formation (i.e., stellar outflows and supernovae

explosions), some of which is converted into turbulent energy through a turbulent energy cascade from galaxy scales to smaller physical scales, may contribute significantly to the CO cooling budget over the molecular gas depletion timescale.

## 6. Summary and Conclusion

SPT0311-58 is comprised of two intensely star-forming galaxies, W and E, at a redshift of 6.9, in the EoR. We characterize the physical properties of the galaxies using new observations of CO(6–5), CO(7–6), CO(10–9), [C I](2–1), and p-H₂O(2₁,₁ – 2₀,₂) transitions. We perform lensing reconstruction by assuming a Sérsic source profile using *visilens* (Spilker et al. 2016). We show that there is no significant differential magnification between dust and CO. We run non-LTE LVG radiative transfer models, which consider dust, CO, and [C I] simultaneously: 1-component, 2-component, and Turbulence models. In the 1-component and 2-component models, the ISM is modeled by one and two gas densities, respectively. The Turbulence model is more sophisticated because the gas density is sampled from a log-normal PDF in turbulent gas. The main results and conclusion are given below:

1. We detect p-H₂O(2₁,₁ – 2₀,₂) in SPT0311-58 W, which is the most distant detection of water in a galaxy without any evidence of an AGN in the literature. The $L_{\text{H}_2\text{O}}/L_{\text{FIR}}$ ratio in SPT0311-58 is consistent with other high-redshift galaxies. From the $L_{\text{H}_2\text{O}}-L_{\text{FIR}}$ correlation, we estimate a SFR of 4356 ± 2143 $M_{\odot}$ yr⁻¹ in W and an upper limit of 385 $M_{\odot}$ yr⁻¹ in E. The SFR calculated from $L_{\text{IR}}$ is 5046 ± 944 $M_{\odot}$ yr⁻¹ in W and 701 ± 151 $M_{\odot}$ yr⁻¹ in E. Both the calculations give consistent values within errors in W. These measurements of H₂O are broadly consistent with the possibility that the cascade transition from the FIR absorption pathways may trace the total FIR luminosity and thus star formation.

2. The CO SLED and brightness temperature ratios in SPT0311-58 W and E are consistent with the other high-redshift starburst galaxies. We explore $L_{\text{[C II]}}/L_{\text{FIR}}$ versus $L_{\text{CO(10–9)}/L_{\text{CO(6–5)}}$ as an indicator of the presence of heating mechanisms in addition to photoelectric heating. We observe that $L_{\text{[C II]}}/L_{\text{FIR}}$ decreases with an increasing $L_{\text{CO(10–9)}/L_{\text{CO(6–5)}}$. The ratio of $L_{\text{CO(10–9)}/L_{\text{CO(6–5)}}$ in SPT0311-58 is comparable to low-redshift AGN and Class II and III galaxies (Rosenberg et al. 2015), which suggests that there are additional heating mechanisms, but we cannot confirm the presence of an AGN.

3. By comparing the radial extent of dust and CO transitions, we observe that the CO is emitted from a larger area than the dust, although the result is not significant. We also observe that the CO emission region is marginally decreasing with an increasing excitation level.

4. The mean density derived from the Turbulence LVG model is $\log(n_{\text{H}_2}/\text{cm}^{-3}) = (5.3 \pm 1.3)$ in W and $(4.7 \pm 2.4)$ in E. The gas mass from the Turbulence model is $5.4 \pm 3.4 \times 10^{11} M_{\odot}$ in W and $3.1 \pm 2.7 \times 10^{10} M_{\odot}$ in E. From the predicted LVG model CO(1–0), we estimate a gas conversion factor $\alpha_{\text{CO}} = 7.1 \pm 5.3$ and $3.9 \pm 4.4$ in W and E, respectively. This is consistent with the $\alpha_{\text{CO}}$ in the high-redshift SMGs within the uncertainties. From the high-redshift SMG sample
together with the SPT-SMGs, we estimate a mean of $\alpha_{CO} = 3.2 \pm 2.7$.

5. From the Turbulence model, we estimate a depletion timescale of 107 ± 70 Myr in W and 44 ± 40 Myr in E. We observe that there is no evidence for the evolution of depletion time with the redshift among the SMG sample at $z > 3$. The gas depletion time in SPT0311-58 W and E is within the range of the other high-redshift starburst galaxies.

6. The ratio of $T_{kin}/T_{dust}$ is >1 in both galaxies, which indicates that there are additional heating mechanisms such as X-rays, stellar outflows, and supernovae in addition to photoelectric heating in dense molecular gas. The mechanical heating from stellar outflows and supernovae explosions, some of which is converted into turbulent energy, may contribute significantly to the total CO cooling over the depletion timescale.

SPT0311-58 is one of the most well-characterized galaxies in the EoR. We observed the brightest source at a high spatial resolution of ~2–3 kpc and could detect the CO lines with a peak signal-to-noise ratio greater than 4. These resolved observations at $z \sim 7$ highlight the power of the ALMA. SPT0311-58 is the highest redshift source from the SPT-SZ survey, and it is expected that more than 100 sources at $z > 7$ can be found in the SPT-3G survey (e.g., Benson et al. 2014; Guns et al. 2021).

The SPT is supported by the NSF through grant OPP-1852617. D.P.M., J.D.V., and S.J. acknowledge support from the US NSF under grants AST-1715213 and AST-1716127. S.J. acknowledges support from the US NSF NRAO under grant SOSPA7-006. J.D.V. acknowledges support from an A. P. Sloan Foundation Fellowship. M.A. and J.D.V. acknowledge support from the Center for AstroPhysical Surveys at the National Center for Supercomputing Applications in Urbana, IL. M.A. has been supported by the grant CONICYT+PCI+REDES 190194. T.R.G. acknowledges the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation under grant No. 140. The Flatiron Institute is supported by the

Simons Foundation. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2017.1.01168.S, ADS/JAO.ALMA #2016.1.01293.S, and ADS/JAO.ALMA #2015.1.00504.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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Software: CASA (v5.1.1; McMullin et al. 2007), emcee (Foreman-Mackey et al. 2013), visilens (Spilker et al. 2016).

# Appendix

## A.1. Far-Infrared Photometry

The intrinsic FIR photometry in SPT0311-58 W and E are shown in Table A.1.1.

<table><caption>Table A.1.1
Intrinsic Far-Infrared Photometry</caption><thead><tr><th>Wavelength (µm)</th><th>East (mJy)</th><th>West (mJy)</th></tr></thead><tbody><tr><td>250</td><td>1.9 ± 0.6</td><td>12.7 ± 4.4</td></tr><tr><td>350</td><td>2.5 ± 0.6</td><td>16.6 ± 3.3</td></tr><tr><td>500</td><td>3.5 ± 0.7</td><td>22.7 ± 4.8</td></tr><tr><td>710</td><td>3.1 ± 0.4</td><td>19.9 ± 2.0</td></tr><tr><td>869</td><td>2.9 ± 0.4</td><td>15.9 ± 1.6</td></tr><tr><td>1260</td><td>1.18 ± 0.03</td><td>9.8 ± 1.0</td></tr><tr><td>2140</td><td>0.18 ± 0.03</td><td>2.0 ± 0.3</td></tr><tr><td>3150</td><td>0.03 ± 0.02</td><td>0.6 ± 0.1</td></tr></tbody></table>

Note. The continuum photometry given in the table is corrected for magnification. The photometry from 250 to 1260 µm are taken from Extended Data Table 3 in Marrone et al. (2018). The continuum flux densities at 2140 µm and 3150 µm (140 GHz and 95 GHz, respectively) are given in Table A.2.2. A 15% calibration error is added in quadrature to the statistical error to account for the uncertainty from the absolute flux calibration and the lens modeling.
## A.2. Lens Models

We present the source plane reconstruction of the continuum and spectral lines in W in Figures 14 and 15. The first two panels show the observations and the beam-convolved image from the model, respectively. The third panel is the high-resolution image

from the lens model with the lensing caustic shown. The final panel shows the source plane reconstructed image. The best-fit lens parameters from the lens modeling are given in Table A.2.1. The best-fit source plane continuum and line parameters are given in Tables A.2.2 and A.2.3, respectively.

### 95 GHz Continuum

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

### 140 GHz Continuum

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

**Figure 14.** Lens models in continuum at 95 GHz (top panel) and 140 GHz (bottom panel). The first panel is from the observations and the second panel is the high-resolution model convolved with the telescope beam. The contours are at ± [3, 5, 10, 20, 40, and 80] × σ where σ is the rms noise in the map. The residual (i.e., observed - model) contours are shown in blue at ± [2, 3, 4, and 5] × σ. The third panel is the high-resolution source plane obtained from the model with the lensing caustic shown in black. The last panel is the source plane reconstruction. The third and fourth panels are shown in logarithmic scale to emphasize the features.
CO(6-5)

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

CO(7-6)

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

CO(10-9)

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Lens models in the velocity-integrated single-channel CO(6–5) (top), CO(7–6) (middle), and CO(10–9) (bottom). The contours in the dirty images are at ± [3, 5, 7, 9, and 11] × σ. The residual description is the same as Figure 14.</figcaption></figure>

<table><caption>**Table A.2.1**

Lens parameters</caption><thead><tr><th>x<sub>L</sub><br>(")</th><th>y<sub>L</sub><br>(")</th><th>e<sub>L</sub></th><th>M<sub>L</sub><br>(10<sup>11</sup> M<sub>☉</sub>)</th><th>θ<sub>L</sub><br>(degree)</th></tr></thead><tbody><tr><td>-0.31 ± 0.02</td><td>-0.15 ± 0.01</td><td>0.68 ± 0.07</td><td>0.30 ± 0.02</td><td>56.09 ± 3.99</td></tr></tbody></table>

Note. $x_L$ and $y_L$ is the position of the lens relative to the phase center. $e_L$ is the ellipticity of the lens. $M_L$ is the mass of the lens. $\theta_L$ is the position angle of the major axis counter-clockwise from East. These parameters are consistent with the lens parameters from Marrone et al. (2018).
<table><caption>**Table A.2.2**
Continuum Source parameters in W</caption><thead><tr><th>ν<sub>obs</sub><br>(GHz)</th><th>x<sub>s</sub><br>(")</th><th>y<sub>s</sub><br>(")</th><th>S<sub>cont</sub><br>(mJy)</th><th>a<sub>s</sub><br>(")</th><th>n<sub>s</sub></th><th>b<sub>s</sub>/a<sub>s</sub></th><th>φ<sub>s</sub><br>(degree)</th></tr></thead><tbody><tr><td>95</td><td>0.24 ± 0.01</td><td>0.04 ± 0.01</td><td>0.60 ± 0.02</td><td>0.33 ± 0.01</td><td>0.67 ± 0.17</td><td>0.15 ± 0.02</td><td>113 ± 1</td></tr><tr><td rowspan="3">140</td><td>0.24 ± 0.01</td><td>0.11 ± 0.01</td><td>1.46 ± 0.06</td><td>0.29 ± 0.01</td><td>0.84 ± 0.09</td><td>0.32 ± 0.03</td><td>111 ± 2</td></tr><tr><td>0.34 ± 0.02</td><td>-0.24 ± 0.02</td><td>0.19 ± 0.07</td><td>0.08 ± 0.01</td><td>0.54 ± 0.17</td><td>0.58 ± 0.09</td><td>-72 ± 10</td></tr><tr><td>0.54 ± 0.04</td><td>-0.53 ± 0.04</td><td>0.06 ± 0.02</td><td>0.08 ± 0.02</td><td></td><td></td><td></td></tr></tbody></table>

**Note.** The best model includes one source at 95 GHz and three sources at 140 GHz. $x_s$ and $y_s$ is the position relative to the lens. $S_{cont}$ is the continuum flux density of the source. $a_s$ is the half light radius of the major axis of the Sérsic profile or radius of the Gaussian. $n_s$ is the Sérsic index (fixed to 0.5 for a Gaussian source). $b_s/a_s$ is the axis ratio. $\phi_s$ is the position of the angle counter-clockwise from East.

<table><caption>**Table A.2.3**
Single-Channel CO Source Parameters in W</caption><thead><tr><th>Line</th><th>x<sub>s</sub><br>(")</th><th>y<sub>s</sub><br>(")</th><th>S<sub>line</sub><br>(mJy)</th><th>a<sub>s</sub><br>(")</th><th>n<sub>s</sub></th><th>b<sub>s</sub>/a<sub>s</sub></th><th>φ<sub>s</sub><br>(degree)</th></tr></thead><tbody><tr><td>CO(6–5)</td><td>0.20 ± 0.02</td><td>0.08 ± 0.04</td><td>0.52 ± 0.05</td><td>0.48 ± 0.06</td><td>0.84 ± 0.47</td><td>0.22 ± 0.07</td><td>109 ± 4</td></tr><tr><td>CO(7–6)</td><td>0.21 ± 0.02</td><td>0.06 ± 0.03</td><td>0.61 ± 0.06</td><td>0.38 ± 0.05</td><td>0.96 ± 0.53</td><td>0.21 ± 0.06</td><td>113 ± 4</td></tr><tr><td>CO(10–9)</td><td>0.19 ± 0.01</td><td>0.10 ± 0.02</td><td>0.35 ± 0.03</td><td>0.38 ± 0.03</td><td>0.32 ± 0.23</td><td>0.18 ± 0.03</td><td>110 ± 3</td></tr></tbody></table>

**Note.** $x_s$ and $y_s$ is the position relative to the lens. $S_{line}$ is the line flux density of the source. $a_s$ is the half light radius of the major axis. $n_s$ is the Sérsic index. $b_s/a_s$ is the axis ratio. $\phi_s$ is the position of the angle counter-clockwise from East.

## A.3. LVG Models

The best-fit 1-component and 2-component model parameters are shown in Table A.3.1.

We also present the output parameters from the ~10<sup>7</sup> Turbulence models in W and E in Figures 16 and 17, respectively. The mean value is indicated in the parameter histograms.

<table><caption>**Table A.3.1**
1- and 2-component Model Parameters</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">1-component model</th><th colspan="4">2-component model</th></tr><tr><th>W</th><th>E</th><th colspan="2">W</th><th colspan="2">E</th></tr><tr><th></th><th></th><th></th><th>component 1</th><th>component 2</th><th>component 1</th><th>component 2</th></tr></thead><tbody><tr><td colspan="7"><strong>Model input</strong></td></tr><tr><td>log(n<sub>H2</sub>)</td><td>4.0 ± 1.6</td><td>3.9 ± 1.1</td><td>3.9 ± 0.4</td><td>5.6 ± 1.1</td><td>3.6 ± 0.4</td><td>5.5 ± 1.2</td></tr><tr><td>T<sub>K</sub></td><td>222 ± 79</td><td>225 ± 94</td><td>88 ± 65</td><td>168 ± 93</td><td>100 ± 80</td><td>142 ± 76</td></tr><tr><td>T<sub>K</sub>/T<sub>dust</sub></td><td>4.0 ± 1.3</td><td>3.9 ± 1.6</td><td>2.7 ± 1.4</td><td>1.5 ± 0.5</td><td>2.5 ± 1.5</td><td>1.7 ± 0.8</td></tr><tr><td>β<sub>Tdust</sub></td><td>1.8 ± 0.1</td><td>1.9 ± 0.1</td><td>1.9 ± 0.1</td><td>1.9 ± 0.1</td><td>2.0 ± 0.1</td><td>2.0 ± 0.1</td></tr><tr><td>R<sub>eff</sub></td><td>1676 ± 646</td><td>702 ± 844</td><td>2458 ± 920</td><td>651 ± 380</td><td>1028 ± 918</td><td>356 ± 778</td></tr><tr><td>κ<sub>vir</sub></td><td>1.4 ± 0.6</td><td>1.4 ± 0.5</td><td>1.9 ± 0.6</td><td>1.7 ± 0.4</td><td>2.0 ± 0.6</td><td>1.8 ± 0.4</td></tr><tr><td>dν<sub>turb</sub></td><td>161 ± 40</td><td>156 ± 38</td><td>134 ± 42</td><td>146 ± 32</td><td>112 ± 42</td><td>116 ± 44</td></tr><tr><td>GDMR</td><td>110 ± 56</td><td>149 ± 131</td><td>155 ± 97</td><td>155 ± 97</td><td>245 ± 148</td><td>245 ± 148</td></tr><tr><td>[CO/H<sub>2</sub>]</td><td>(6 ± 5) × 10<sup>-5</sup></td><td>(10 ± 7) × 10<sup>-5</sup></td><td>(8 ± 5) × 10<sup>-5</sup></td><td>(7 ± 6) × 10<sup>-5</sup></td><td>(12 ± 5) × 10<sup>-5</sup></td><td>(7 ± 6) × 10<sup>-5</sup></td></tr><tr><td>[C I/H<sub>2</sub>]</td><td>(6 ± 3) × 10<sup>-5</sup></td><td>(3 ± 2) × 10<sup>-5</sup></td><td>(6 ± 2) × 10<sup>-5</sup></td><td>(2 ± 2) × 10<sup>-5</sup></td><td>(2 ± 1) × 10<sup>-5</sup></td><td>(0.6 ± 0.6) × 10<sup>-5</sup></td></tr><tr><td colspan="7"><strong>Estimated<br>within the model</strong></td></tr><tr><td>T<sub>dust</sub></td><td>54 ± 4</td><td>58 ± 4</td><td>31 ± 9</td><td>106 ± 37</td><td>38 ± 16</td><td>80 ± 16</td></tr><tr><td>M<sub>gas</sub></td><td>(1.4 ± 1.2) × 10<sup>11</sup></td><td>(1.8 ± 2.5) × 10<sup>10</sup></td><td>(2.4 ± 2.4) × 10<sup>11</sup></td><td>(1.1 ± 2.0) × 10<sup>11</sup></td><td>(2.0 ± 2.3) × 10<sup>10</sup></td><td>(1.6 ± 2.2) × 10<sup>10</sup></td></tr><tr><td>M<sub>dust</sub></td><td>(1.2 ± 1.3) × 10<sup>9</sup></td><td>(1.2 ± 2.0) × 10<sup>8</sup></td><td>(1.5 ± 1.8) × 10<sup>9</sup></td><td>(0.7 ± 1.3) × 10<sup>9</sup></td><td>(0.8 ± 1.0) × 10<sup>8</sup></td><td>(0.7 ± 1.0) × 10<sup>8</sup></td></tr><tr><td colspan="7"><strong>Derived<br>from the model</strong></td></tr><tr><td>L'<sub>CO</sub></td><td>(6 ± 3) × 10<sup>10</sup></td><td>(9 ± 4) × 10<sup>9</sup></td><td>(5 ± 1) × 10<sup>10</sup></td><td>(1 ± 0.6) × 10<sup>10</sup></td><td>(2 ± 1) × 10<sup>9</sup></td><td>(8 ± 3) × 10<sup>8</sup></td></tr><tr><td>L<sub>FIR</sub></td><td>(19 ± 2) × 10<sup>12</sup></td><td>(31 ± 2) × 10<sup>11</sup></td><td>(2 ± 2) × 10<sup>12</sup></td><td>(16 ± 2) × 10<sup>12</sup></td><td>(6 ± 6) × 10<sup>11</sup></td><td>(25 ± 5) × 10<sup>11</sup></td></tr></tbody></table>

**Note.** The input and derived parameters of the 1-component and 2-component models assuming less than solar metallicity. The units and the explored range are same as the Turbulence model given in Table 4. The model outputs are the intrinsic source properties as we use the magnification corrected photometry and the line flux densities for modeling.
<figure><img src="image_21.png" /><figcaption>**Figure 16.** Turbulence model output in W. The contours are at 1, 2, and 3σ confidence levels. The mean value of each parameter is shown as a black vertical line on the histogram. Note that $\beta_{[CI]}$ is not constrained because of a lack of observations of [CI](1–0).</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 17.** Turbulence model output in E. The contours are at 1, 2, and 3σ confidence levels. The mean value of each parameter is shown as a black vertical line on the histogram. Note that $β_{[CI]}$ is not constrained because of a lack of observations of [CI](1–0) and [CI](2–1).</figcaption></figure>

## ORCID iDs

Sreevani Jarugula https://orcid.org/0000-0002-5386-7076

Joaquin D. Vieira https://orcid.org/0000-0001-7192-3871

Axel Weiss https://orcid.org/0000-0003-4678-3939

Justin S. Spilker https://orcid.org/0000-0003-3256-5615

Manuel Aravena https://orcid.org/0000-0002-6290-3198

Melanie Archipley https://orcid.org/0000-0002-0517-9842

Matthieu Béthermin https://orcid.org/0000-0002-3915-2015

Chenxing Dong https://orcid.org/0000-0002-5823-0349

Thomas R. Greve https://orcid.org/0000-0002-2554-1837

Christopher C. Hayward https://orcid.org/0000-0003-4073-3236

Yashar Hezaveh https://orcid.org/0000-0002-8669-5733

Katrina C. Litke https://orcid.org/0000-0002-4208-3532

Matthew A. Malkan https://orcid.org/0000-0001-6919-1237

Daniel P. Marrone https://orcid.org/0000-0002-2367-1080

Desika Narayanan https://orcid.org/0000-0002-7064-4309

Kedar A. Phadke https://orcid.org/0000-0001-7946-557X

Cassie Reuter https://orcid.org/0000-0001-7477-1586

## References

Alaghband-Zadeh, S., Chapman, S. C., Swinbank, A. M., et al. 2013, MNRAS, 435, 1493

Andreani, P., Retana-Montenegro, E., Zhang, Z.-Y., et al. 2018, A&A, 615, A142

Apostolovski, Y., Aravena, M., Anguita, T., et al. 2019, A&A, 628, A23

Aravena, M., Spilker, J. S., Bethermin, M., et al. 2016, MNRAS, 457, 4406

Benson, B. A., Ade, P. A. R., Ahmed, Z., et al. 2014, Proc. SPIE, 9153, 915312
Bisbas, T. G., Papadopoulos, P. P., & Viti, S. 2015, ApJ, 803, 37

Blain, A. W., & Longair, M. S. 1993, MNRAS, 264, 509

Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Bothwell, M. S., Aguirre, J. E., Aravena, M., et al. 2017, MNRAS, 466, 2825

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

Bradford, C. M., Aguirre, J. E., Aikin, R., et al. 2009, ApJ, 705, 112

Brauher, J. R., Dale, D. A., & Helou, G. 2008, ApJS, 178, 280

Bussmann, R. S., Pérez-Fournon, I., Amber, S., et al. 2013, ApJ, 779, 25

Calistro Rivera, G., Hodge, J. A., Smail, I., et al. 2018, ApJ, 863, 56

Carilli, C. L., Daddi, E., Riechers, D., et al. 2010, ApJ, 714, 1407

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

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

Cen, R., & Ostriker, J. P. 1999, ApJL, 519, L109

Cochrane, R. K., Hayward, C. C., Anglés-Alcázar, D., et al. 2019, MNRAS, 488, 1779

Conley, A., Cooray, A., Vieira, J. D., et al. 2011, ApJL, 732, L35

Cooray, A., Calanog, J., Wardlow, J. L., et al. 2014, ApJ, 790, 40

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2011, MNRAS, 410, 1687

Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98

De Breuck, C., Weiß, A., Bétheimin, M., et al. 2019, A&A, 631, A167

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2014, ApJL, 788, L17

Dong, C., Spilker, J. S., Gonzalez, A. H., et al. 2019, ApJ, 873, 50

Downes, D., & Solomon, P. M. 1998, ApJ, 507, 615

Draine, B. T. 2003, ARA&A, 41, 241

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Elia, D., Molinari, S., Schisano, E., et al. 2017, MNRAS, 471, 100

Everett, W. B., Zhang, L., Crawford, T. M., et al. 2020, ApJ, 900, 55

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Frayer, D. T., Harris, A. I., Baker, A. J., et al. 2011, ApJL, 726, L22

Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590

Fu, H., Cooray, A., Feruglio, C., et al. 2013, Natur, 498, 338

Fu, H., Jullo, E., Cooray, A., et al. 2012, ApJ, 753, 134

Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, A&A, 624, A105

Glover, S. C. O., & Clark, P. C. 2016, MNRAS, 456, 3596

Glover, S. C. O., Clark, P. C., Micic, M., & Molina, F. 2015, MNRAS, 448, 1607

Goldsmith, P. F. 2001, ApJ, 557, 736

González-Alfonso, E., Fischer, J., Graciá-Carpio, J., et al. 2012, A&A, 541, A4

González-Alfonso, E., Fischer, J., Isaak, K., et al. 2010, A&A, 518, L43

Greve, T. R., Leonidaki, I., Xilouris, E. M., et al. 2014, ApJ, 794, 142

Greve, T. R., Papadopoulos, P. P., Gao, Y., & Radford, S. J. E. 2009, ApJ, 692, 1432

Greve, T. R., Vieira, J. D., Weiß, A., et al. 2012, ApJ, 756, 101

Gullberg, B., De Breuck, C., Vieira, J. D., et al. 2015, MNRAS, 449, 2883

Guns, S., Foster, A., Daley, C., et al. 2021, ApJ, 916, 98

Harrington, K. C., Weiss, A., Yun, M. S., et al. 2021, ApJ, 908, 95

Harris, A. I., Baker, A. J., Zonak, S. G., et al. 2010, ApJ, 723, 1139

Hayward, C. C., Sparre, M., Chapman, S. C., et al. 2021, MNRAS, 502, 2922

Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288

Helou, G., Soifer, B. T., & Rowan-Robinson, M. 1985, ApJL, 298, L7

Hezaveh, Y. D., Dalal, N., Marrone, D. P., et al. 2016, ApJ, 823, 37

Hezaveh, Y. D., Marrone, D. P., & Holder, G. P. 2012, ApJ, 761, 20

Iffrig, O., & Hennebelle, P. 2015, in SF2A-2015: Proc. of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. F. Martins et al. (Toulouse: Société Française d'Astronomie et d'Astrophysique), 333

Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913

Ivison, R. J., Swinbank, A. M., Smail, I., et al. 2013, ApJ, 772, 137

Jarugula, S., Vieira, J. D., Spilker, J. S., et al. 2019, ApJ, 880, 92

Jiao, Q., Zhao, Y., Lu, N., et al. 2019, ApJ, 880, 133

Jiao, Q., Zhao, Y., Zhu, M., et al. 2017, ApJL, 840, L18

Kamenetzky, J., Privon, G. C., & Narayanan, A. 2008, ApJ, 859, 9

Kamenetzky, J., Rangwala, N., Glenn, J., Maloney, P. R., & Conley, A. 2016, ApJ, 829, 93

Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, ApJ, 527, 795

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

Kirkpatrick, A., Sharon, C., Keller, E., & Pope, A. 2019, ApJ, 879, 41

Koptelova, E., & Hwang, C.-Y. 2019, ApJL, 880, L19

Krumholz, M. R. 2014, PhR, 539, 49

Krumholz, M. R., & McKee, C. F. 2005, ApJ, 630, 250

Kulesa, C. A. 2002, PhD thesis, Univ. Arizona

Larson, R. B. 2005, MNRAS, 359, 211

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Leroy, A. K., Bolatto, A., Gordon, K., et al. 2013, ApJ, 737, 12

Li, J., Wang, R., Riechers, D., et al. 2020, ApJ, 889, 162

Li, Q., Narayanan, D., & Davé, R. 2019, MNRAS, 490, 1425

Litke, K. C., Marrone, D. P., Spilker, J. S., et al. 2019, ApJ, 870, 80

Liu, L., Weiß, A., Perez-Beaupuits, J. P., et al. 2017, ApJ, 846, 5

Lu, N., Zhao, Y., Díaz-Santos, T., et al. 2017, ApJS, 230, 1

Lu, N., Zhao, Y., Xu, C. K., et al. 2015, ApJL, 802, L11

Ma, J., Gonzalez, A. H., Vieira, J. D., et al. 2016, ApJ, 832, 114

Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429

Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132

Magnelli, B., Saintonge, A., Lutz, D., et al. 2012, A&A, 548, A22

Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766

Maloney, P., & Black, J. H. 1988, ApJ, 325, 389

Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Natur, 553, 51

Martizzi, D., Fielding, D., Faucher-Giguère, C.-A., & Quataert, E. 2016, MNRAS, 459, 2311

Mashian, N., Sternberg, A., & Loeb, A. 2013, MNRAS, 435, 2407

Mashian, N., Sturm, E., Sternberg, A., et al. 2015, ApJ, 802, 81

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XV, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJL, 762, L16

Meijerink, R., Spaans, M., & Israel, F. P. 2007, A&A, 461, 793

Muñoz, J. A., & Oh, S. P. 2016, MNRAS, 463, 2085

Narayanan, D. 2011, arXiv:1112.1073

Narayanan, D., Krumholz, M. R. 2014, MNRAS, 442, 1411

Narayanan, D., & Krumholz, M. R. 2017, MNRAS, 467, 50

Narayanan, D., Turk, M., Feldmann, R., et al. 2015, Natur, 525, 496

Neufeld, D. A., Lepp, S., & Melnick, G. J. 1995, ApJS, 100, 132

Novak, M., Bañados, E., Decarli, R., et al. 2019, ApJ, 881, 63

Omont, A., Yang, C., Cox, P., et al. 2013, A&A, 551, A115

Papadopoulos, P. P., Thi, W.-F., & Viti, S. 2004, MNRAS, 351, 147

Papadopoulos, P. P., van der Werf, P., Xilouris, E., Isaak, K. G., & Gao, Y. 2012, ApJ, 751, 10

Pham, D. T., & Castellani, M. 2009, Proc. Inst. Mech. Eng. C, 223, 2919

Planck Collaboration, Ade, P. A. R., & Aghanim, N. 2016, A&A, 594, A13

Ramachandran, V., Hamann, W. R., Oskinova, L. M., et al. 2019, A&A, 625, A104

Reuter, C., Vieira, J. D., Spilker, J. S., et al. 2020, ApJ, 902, 78

Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Natur, 496, 329

Riechers, D. A., Walter, F., Carilli, C. L., & Lewis, G. F. 2009, ApJ, 690, 463

Rosenberg, M. J. F., van der Werf, P. P., Alto, S., et al. 2015, ApJ, 801, 72

Rybak, M., Hodge, J. A., Vegetti, S., et al. 2020, MNRAS, 494, 5542

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Sanders, D. B., Mazzarella, J. M., Kim, D.-C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607

Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5

Scoville, N., Lee, N., LeFloch, E., & Sanders, D. 2012, AAS Meeting, 220, 517:05

Sersic, J. L. 1968, Atlas de galaxias australes (Cordoba, Argentina: Observatorio Astronomico)

Shao, Y., Wang, R., Carilli, C. L., et al. 2019, ApJ, 876, 99

Simpson, J. M., Swinbank, A. M., Smail, I., et al. 2014, ApJ, 788, 125

Smith, J. D. T., Croxall, K., Draine, B., et al. 2017, ApJ, 834, 5

Smith, N. 2014, ARA&A, 52, 487

Solomon, P. M., Downes, D., Radford, S. J. E., & Barrett, J. W. 1997, ApJ, 478, 144

Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. 2002, J. R. Stat. Soc. B, 64, 583

Spilker, J. S., Aravena, M., Marrone, D. P., et al. 2015, ApJ, 811, 124

Spilker, J. S., Marrone, D. P., Aguirre, J. E., et al. 2014, ApJ, 785, 149

Spilker, J. S., Marrone, D. P., Aravena, M., et al. 2016, ApJ, 826, 112

Strandet, M. L., Weiss, A., De Breuck, C., et al. 2017, ApJL, 842, L15

Strandet, M. L., Weiss, A., Vieira, J. D., et al. 2016, ApJ, 822, 80

Stutzki, J., Graf, U. U., Haas, S., et al. 1997, ApJL, 477, L33

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Tielens, A. G. G. M., & Hollenbach, D. J. 1985, ApJ, 291, 722

Tunnard, R., & Greve, T. R. 2016, ApJ, 819, 161
Uzgil, B. D., Bradford, C. M., Hailey-Dunsheath, S., Maloney, P. R., & Aguirre, J. E. 2016, ApJ, 832, 209

Valentino, F., Magdis, G. E., Daddi, E., et al. 2020, ApJ, 890, 24

van der Werf, P. P., Isaak, K. G., Meijerink, R., et al. 2010, A&A, 518, L42

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

Walter, F., Decarli, R., Carilli, C., et al. 2012, Natur, 486, 233

Walter, F., Weiß, A., Downes, D., Decarli, R., & Henkel, C. 2011, ApJ, 730, 18

Wang, F., Wang, R., Fan, X., et al. 2019, ApJ, 880, 2

Weiß, A., Downes, D., Henkel, C., & Walter, F. 2005a, A&A, 429, L25

Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955

Weiß, A., Henkel, C., Downes, D., & Walter, F. 2003, A&A, 409, L41

Weiß, A., Walter, F., & Scoville, N. Z. 2005b, A&A, 438, 533

White, G. J., Ellison, B., Claude, S., Dent, W. R. F., & Matheson, D. N. 1994, A&A, 284, L23

Yang, C., Gao, Y., Omont, A., et al. 2013, ApJL, 771, L24

Yang, C., Omont, A., Beelen, A., et al. 2016, A&A, 595, A80

Yang, C., Omont, A., Beelen, A., et al. 2017, A&A, 608, A144

Yang, J., Venemans, B., Wang, F., et al. 2019, ApJ, 880, 153

Zhang, Z.-Y., Ivison, R. J., George, R. D., et al. 2018, MNRAS, 481, 59