# Multiphase ISM in the z = 5.7 Hyperluminous Starburst SPT 0346-52

Katrina C. Litke, Daniel P. Marrone, Manuel Aravena, Matthieu Béthermin, Scott C. Chapman, Chenxing Dong, Christopher C. Hayward, Ryley Hill, Sreevani Jarugula, Matthew A. Malkan, et al.

► To cite this version:

Katrina C. Litke, Daniel P. Marrone, Manuel Aravena, Matthieu Béthermin, Scott C. Chapman, et al. Multiphase ISM in the z = 5.7 Hyperluminous Starburst SPT 0346-52. The Astrophysical Journal, 2022, 928, 10.3847/1538-4357/ac58f9. insu-03667440

HAL Id: insu-03667440

https://insu.hal.science/insu-03667440

Submitted on 14 May 2022

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
OPEN ACCESS

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

CrossMark

# Multiphase ISM in the $z=5.7$ Hyperluminous Starburst SPT 0346-52

Katrina C. Litke<fnref n="1" />, Daniel P. Marrone<fnref n="1" />, Manuel Aravena<fnref n="2" />, Matthieu Béthermin<fnref n="3" />, Scott C. Chapman<fnref n="4" />,<fnref n="5" />, Chenxing Dong<fnref n="6" />, Christpher C. Hayward<fnref n="7" />, Ryley Hill<fnref n="8" />, Sreevani Jarugula<fnref n="9" />, Matthew A. Malkan<fnref n="10" />, Desika Narayanan<fnref n="6" />,<fnref n="11" />,<fnref n="12" />, Cassie A. Reuter<fnref n="9" />, Justin S. Spilker<fnref n="13" />, Nikolaus Sulzenauer<fnref n="14" />, Joaquin D. Vieira<fnref n="9" />, and Axel Weiß<fnref n="14" />

*Received 2021 April 28; revised 2022 February 15; accepted 2022 February 24; published 2022 April 7*

## Abstract

With $\Sigma_{SFR} \sim 4200 M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}$, SPT 0346-52 ($z=5.7$) is the most intensely star-forming galaxy discovered by the South Pole Telescope. In this paper, we expand on previous spatially resolved studies, using ALMA observations of dust continuum, [N II] 205 µm, [C II] 158 µm, [O I] 146 µm, and undetected [N II] 122 µm and [O I] 63 µm emission to study the multiphase interstellar medium (ISM) in SPT 0346-52. We use pixelated, visibility-based lens modeling to reconstruct the source-plane emission. We also model the source-plane emission using the photoionization code CLOUDY and find a supersolar metallicity system. We calculate $T_{\text{dust}} = 48.3$ K and $\lambda_{\text{peak}} = 80$ µm and see line deficits in all five lines. The ionized gas is less dense than comparable galaxies, with $n_e < 32 \text{ cm}^{-3}$, while ~20% of the [C II] 158 µm emission originates from the ionized phase of the ISM. We also calculate the masses of several phases of the ISM. We find that molecular gas dominates the mass of the ISM in SPT 0346-52, with the molecular gas mass ~4× higher than the neutral atomic gas mass and ~100× higher than the ionized gas mass.

*Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Interstellar medium (847); Starburst galaxies (1570)*

## 1. Introduction

The interstellar medium (ISM) of high-redshift galaxies is difficult to study directly due to cosmological dimming and angular resolution limitations. Observations of rest-frame optical and ultraviolet wavelengths also suffer from significant dust extinction, and some phases of the ISM lack suitable tracers at these wavelengths. However, for very distant objects, far-infrared (FIR) continuum and line emission that is normally obscured by the Earth's atmosphere redshifts into the submillimeter window. The angular resolution and sensitivity afforded by the Atacama Large Millimeter/submillimeter Array (ALMA) is providing new opportunities to explore the physical conditions in early galaxies through their rest-frame FIR emission.

Many recent high-redshift studies (e.g., Gullberg et al. 2015; Béthermin et al. 2020; Le Fèvre et al. 2020) have focused on the 158 µm fine-structure line of singly ionized carbon (hereafter [C II] 158 µm) because it is one of the brightest cooling lines of the ISM (Hollenbach et al. 1991). However, this line can be difficult to interpret because it can originate from both ionized

and neutral gas in photodissociation regions (PDRs). The 122 and 205 µm lines of ionized nitrogen ([N II]) arise from the ionized phase of the ISM because nitrogen has a higher ionization energy than hydrogen. Since [N II] 122 µm and [N II] 205 µm trace the ionized ISM, comparing [C II] 158 µm to [N II] 205 µm emission makes it possible to determine what fraction of the [C II] 158 µm emission originates from PDRs, with values typically in the 60%-90% range (e.g., Pavesi et al. 2016; Díaz-Santos et al. 2017; Herrera-Camus et al. 2018; Cormier et al. 2019). On the other hand, [O I] 63 µm and [O I] 146 µm originate from warm, neutral gas (Tielens & Hollenbach 1985; Hollenbach et al. 1991). Where there is more [O I] 146 µm emission compared to [C II] 158 µm emission, we would expect more dense, neutral gas in those regions (De Breuck et al. 2019).

Recently, (mostly) spatially unresolved multiline surveys of high-$z$ galaxies, including [N II] 205 µm, [C II] 158 µm, [O I] 146 µm, and [N II] 122 µm, have been conducted in individual systems. De Breuck et al. (2019) and Novak et al. (2019) and found highly enriched ISM with approximately solar metallicity in J1342+0928 and SPT 0418-47. De Breuck et al. (2019) and Lee et al. (2021) also found evidence for a dense gas-dominated ISM using the ratio of [O I] 146 µm to [C II] 158 µm in the first detections of [O I] 146 µm at $z > 1$. Rybak et al. (2020) also recently published the first [O I] 63 µm

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA; kclitke@email.arizona.edu</footnote>

<footnote n="2">²Núcleo de Astronomia, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile</footnote>

<footnote n="3">³Aix Marseille Univ., Centre National de la Recherche Scientifique, Laboratoire d'astrophysique de Marseille, Marseille, France</footnote>

<footnote n="4">⁴Eureka Scientific, Inc., 2452 Delmer Street, Suite 100, Oakland, CA 94602-3017, USA</footnote>

<footnote n="5">⁵Dalhousie University, Halifax, NS, Canada</footnote>

<footnote n="6">⁶Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="7">⁷Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA</footnote>

<footnote n="8">⁸Department of Physics and Astronomy, University of British Columbia, 6225 Agricultural Road, Vancouver, V6T 1Z1, Canada</footnote>

<footnote n="9">⁹Department of Astronomy and Department of Physics, University of Illinois, 1002 West Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="10">¹⁰Department of Physics and Astronomy, University of California, Los Angeles, CA 90095-1547, USA</footnote>

<footnote n="11">¹¹University of Florida Informatics Institute, 432 Newell Drive, CISE Bldg E251, Gainesville, FL 32611, USA</footnote>

<footnote n="12">¹²Cosmic Dawn Center (DAWN), DTU-Space, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark; Niels Bohr Institute, University of Copenhagen, Juliane Maries vej 30, DK-2100 Copenhagen, Denmark</footnote>

<footnote n="13">¹³Department of Astronomy, University of Texas at Austin, 2515 Speedway Stop C1400, Austin, TX 78712, USA</footnote>

<footnote n="14">¹⁴Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69 D-53121 Bonn, Germany</footnote>
<table><caption>**Table 1**
ALMA Observations of SPT 0346–52</caption><thead><tr><th>Line</th><th>Date</th><th>Frequency<sup>a</sup><br>(GHz)</th><th>No. Ant.</th><th>Resolution<br>(arcsec)</th><th>Flux<br>Calibrator</th><th>Phase<br>Calibrator</th><th>PWV<sup>b</sup><br>(mm)</th><th>t<sub>int</sub><sup>c</sup><br>(minutes)</th><th>Noise Level<sup>d</sup><br>(mJy beam<sup>-1</sup>)</th><th>Project ID<sup>e</sup></th></tr></thead><tbody><tr><td colspan="11"><strong>Band 6</strong></td></tr><tr><td>[N II] 205 µm</td><td>2015 Aug 30</td><td>227.518</td><td>35</td><td>0.19 × 0.25</td><td>J0334–4008</td><td>J0334–4008</td><td>1.4</td><td>44.3</td><td>0.10</td><td>2013.1.01231</td></tr><tr><td colspan="11"><strong>Band 7</strong></td></tr><tr><td>[C II] 158 µm</td><td>2014 Sep 2</td><td>291.533</td><td>34</td><td>0.22 × 0.26</td><td>J0334–4008</td><td>J0334–4008</td><td>0.9</td><td>5.3</td><td>0.25</td><td>2013.1.01231</td></tr><tr><td>[O I] 146 µm</td><td>2015 Jun 28</td><td>291.536</td><td>41</td><td>0.15 × 0.17</td><td>J0334–4008</td><td>J0334–4008</td><td>1.3</td><td>5.2</td><td>0.23</td><td></td></tr><tr><td>[N II] 122 µm</td><td>2014 Sep 2</td><td>304.136</td><td>34</td><td>0.22 × 0.27</td><td>J2357–5311</td><td>J0334–4008</td><td>0.9</td><td>23.8</td><td>0.20</td><td></td></tr><tr><td>[O I] 146 µm</td><td>2016 Jun 30</td><td>364.434</td><td>40</td><td>0.28 × 0.35</td><td>Ceres</td><td>J0253–5441</td><td>0.6</td><td>24.6</td><td>0.21</td><td>2015.1.01580</td></tr><tr><td>[N II] 122 µm</td><td>2018 Aug 28</td><td>364.431</td><td>48</td><td>0.29 × 0.37</td><td>J0519–4546</td><td>J0253–5441</td><td>0.3</td><td>35.3</td><td>0.27</td><td>2016.1.01565</td></tr><tr><td colspan="11"><strong>Band 9</strong></td></tr><tr><td>[O I] 63 µm</td><td>2018 Aug 17</td><td>711.038</td><td>46</td><td>0.23 × 0.33</td><td>J0522–3627</td><td>J0210–5101</td><td>0.4</td><td>19.2</td><td>1.7</td><td>2016.1.01565</td></tr><tr><td></td><td>2018 Aug 19</td><td>711.038</td><td>45</td><td>0.23 × 0.28</td><td>J0522–3627</td><td>J0210–5101</td><td>0.4</td><td>19.2</td><td>1.5</td><td></td></tr><tr><td></td><td>2018 Sep 03</td><td>711.034</td><td>46</td><td>0.15 × 0.19</td><td>J0522–3627</td><td>J0210–5101</td><td>0.4</td><td>18.5</td><td>1.2</td><td></td></tr><tr><td></td><td>2018 Sep 03</td><td>711.034</td><td>46</td><td>0.16 × 0.20</td><td>J0522–3627</td><td>J0210–5101</td><td>0.3</td><td>19.2</td><td>1.4</td><td></td></tr></tbody></table>

**Notes.**

<sup>a</sup> First local oscillator frequency.

<sup>b</sup> Precipitable water vapor at zenith.

<sup>c</sup> On-source integration time.

<sup>d</sup> rms noise level in continuum image.

<sup>e</sup> PI: Marrone.

detection at $z > 3$ in a dusty galaxy at $z \sim 6$ and determined that [O I] 63 µm was the main neutral gas coolant in G09.83808.

In this paper, we focus on the $z = 5.656$ gravitationally lensed dusty star-forming galaxy (DSFG) SPT-S J034640–5204.9 (hereafter SPT 0346–52; Weiß et al. 2013; Vieira et al. 2013). It is the most intensely star-forming galaxy from the 2500 deg$^2$ South Pole Telescope (SPT) survey (Vieira et al. 2010; Carlstrom et al. 2011; Everett et al. 2020), with an apparent $L_{FIR} = 1.1 \times 10^{14} L_{\odot}$ (Spilker et al. 2015; Reuter et al. 2020) and intrinsic star formation rate density $\Sigma_{SFR} = 4200 M_{\odot} \text{yr}^{-1} \text{ kpc}^{-2}$ (Hezaveh et al. 2013; Ma et al. 2015, 2016; Spilker et al. 2015), where $L_{FIR}$ is the emission from 42.5 to 122.5 µm (Helou et al. 1988). Based on Chandra observations, Ma et al. (2016) determined that the high $L_{FIR}$ is dominated by star formation with negligible contribution from an active galactic nucleus (AGN).

Litke et al. (2019) performed pixelated, interferometric lens modeling of [C II] 158 µm emission in SPT 0346–52. The gas in SPT 0346–52 was found to be globally unstable, with Toomre Q instability parameters $\ll 1$ throughout the system. In addition, they found two components separated by $\sim 1$ kpc and $\sim 500 \text{ km s}^{-1}$ that appear to be merging, which is likely driving the intense star formation in SPT 0346–52. More recently, Jones et al. (2019) suggested that a rotating disk galaxy is a better explanation for a water absorption line.

In this paper, we extend the work of Litke et al. (2019), expanding their [C II] 158 µm analysis to a survey of fine-structure lines. Using [N II] 205 µm, [C II] 158 µm, [O I] 146 µm, [N II] 122 µm, [O I] 63 µm, and the underlying dust continuum emission, we conduct a multiphase study of the ISM in SPT 0346–52. This represents one of the first multiline, spatially resolved studies of the ISM at high-$z$.

We describe ALMA observations of the five fine-structure lines in Section 2. The lensing reconstruction process and results are discussed in Section 3. In Section 4, we describe the CLOUDY modeling of the ISM in SPT 0346–52. We describe the results and various line and continuum diagnostics in Section 5 and summarize the results in Section 6. We adopt the cosmology

of Planck Collaboration et al. (2016; $\Omega_m = 0.309$, $\Omega_\Lambda = 0.691$, and $H_0 = 67.7 \text{ km s}^{-1}$). At $z = 5.656$, $1'' = 6.035$ kpc.

## 2. ALMA Observations

From 2014 September through 2018 September, SPT 0346–52 was observed in ALMA Bands 6, 7, and 9 (Project IDs 2013.1.01231, 2015.1.01580, and 2016.1.01565; PI: Marrone). The [C II] 158 µm, [N II] 122 µm, and [O I] 63 µm lines were all observed on multiple dates at different resolutions, while [N II] 205 µm and [O I] 146 µm were each observed once. The details of these observations, including dates, observing frequencies, flux and phase calibrators, and resolutions, are listed in Table 1.

The data were processed using various pipeline versions of the Common Astronomy Software Applications package (CASA; McMullin et al. 2007; Petry & CASA Development Team 2012). The [C II] 158 µm, [N II] 205 µm, and [O I] 146 µm were all processed using CASA pipeline version 4.2.2, [N II] 122 µm was processed with CASA pipeline version 4.7.1, and [O I] 63 µm was processed with CASA pipeline version 5.4.0. These were the accepted pipeline versions for the cycles in which each data set was observed.

Continuum images at all five frequencies were created with the task TCLEAN in CASA version 5.4.0 using Briggs weighting (robust = 0.5) and the AUTO-MULTITHRESH masking option. The continuum images are shown in the top row of Figure 1 and as contours in the bottom row. For the line emission, the continuum was subtracted from the line cube using the CASA task UVCONTSUB with a first-order polynomial representing the continuum.

The line emission was imaged in the same manner as the continuum but tapered to $300 \text{ k}\lambda$ and integrated from $-300$ to $+300 \text{ km s}^{-1}$ (Figure 1, bottom row). This tapering was also used in the lensing reconstructions.

To evaluate the overall significance of our line detections, we construct source-integrated spectra from the visibility data.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Observed emission from SPT 0346–52. Top: continuum images. Bottom: line images with continuum contours. Continuum images were untapered and used Briggs weighting (robust = 0.5). Line images were made integrating from −300 to +300 km s⁻¹ with Briggs weighting (robust = 0.5) and tapered to 300kλ. The contours represent the observed continuum emission at 10σ, 30σ, and 50σ. From left to right: 205, 158, 146, 122, and 63 µm.</figcaption></figure>

Because SPT 0346–52 is gravitationally lensed, we use the spatial structure of the continuum emission to provide a spatial template for the line emission (see the Appendix of Litke et al. 2019). Visibilities of the line emission are weighted by a gravitational lensing model of the continuum emission to emphasize the visibilities that best sample the source structure, yielding a channelized flux density ($F_v$) determined by

$$F_v = \frac{\sum_i \frac{\tilde{v}_{v,i}}{\tilde{m}_i} |\tilde{m}_i|^2}{\sum_i |\tilde{m}_i|^2}. \qquad (1)$$

Here $\tilde{v}_{v,i}$ is the complex line data visibility, and $\tilde{m}_i$ is the complex model visibility for that data set. The model visibilities are obtained from our lensing reconstructions, described in Section 3. The observed spectrum for each line is shown in Figure 2. To obtain the uncertainties, visibilities were randomly drawn from the distribution of visibilities for each channel and line 500 times. The random spectra were then calculated using Equation (1). The uncertainties were then determined by taking the standard deviation of the 500 random noise trials.

The [N II] 122 µm line is not significantly detected in our observations. It is redshifted to 369.5 GHz, where the atmospheric transmission declines due to a strong atmospheric O₂ line at 368.5 GHz. The [O I] 63 µm line is also not significantly detected in our observations. This line is redshifted to 712.9 GHz, the high-frequency end of ALMA Band 9, where a strong atmospheric O₂ line at 715.4 GHz and the 752 GHz water line that separates the 650 and 850 GHz atmospheric windows (Bands 9 and 10) combine to produce a sharp decline in atmospheric transmission toward higher frequency (bluer velocity). To obtain the upper limits for [N II] 122 µm and [O I] 63 µm listed in Table 2, the 1σ uncertainty on a single 600 km s⁻¹ channel was calculated using the method described above. This value was multiplied by 3 to obtain the 3σ upper limit. It was then divided by $\mu = 5.6$ (Spilker et al. 2016) to correct for magnification from gravitational lensing. The nondetections of [N II] 122 µm and [O I] 63 µm are discussed in Sections 5.5 and 5.6, respectively.

### 3. Lens Modeling

Gravitational lensing is a powerful tool for studying galaxies at high redshift. Because lensing spreads the source emission over a larger solid angle on the sky while preserving the surface brightness, resolving detail in the lensed galaxy can be done with a more compact array configuration than is possible in unlensed sources, and more compact arrays have better surface brightness sensitivity. However, the image distortion introduced by the gravitational lensing makes it difficult to study the spatially resolved physical structure of the galaxy in a straightforward manner.

In order to study the source-plane structure of SPT 0346–52, we turn to lensing reconstruction models. We use a pixelated, interferometric lensing reconstruction code, RIPPLES (Hezaveh et al. 2016). Additional information in the general framework for pixelated lens modeling is described by Warren & Dye (2003) and Suyu et al. (2006). RIPPLES uses a Markov Chain Monte Carlo (MCMC) method to model the mass distribution of the foreground galaxy, as well as the background source emission. It also takes into account observational effects from the primary beam. In addition, a regularization factor is introduced that minimizes large gradients between adjacent pixels, which prevents overfitting of the data.

As mentioned above, RIPPLES models the complex visibilities observed by ALMA directly, rather than modeling CLEAN images. By modeling the complex visibilities, we use all of the data available from the ALMA observations. Because RIPPLES is a pixelated code, we do not assume a source-plane structure and can model more complex structures.

We model the ALMA observations using the same procedure as Litke et al. (2019). The mass distribution of the foreground lensing galaxy is modeled as a singular isothermal ellipsoid (SIE) at $z = 0.9$ with an external shear component. Previous lens modeling of SPT 0346–52 by Hezaveh et al. (2013), Spilker et al. (2015), and Litke et al. (2019) were used to obtain the initial parameters. The 205, 158, 146, 122, and 63 µm rest-frame continuum data were all fit independently. The best-fit lens model derived from each continuum data set
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Observed spectra of the targeted lines. Top row (left to right): [N II] 205 µm, [C II] 158 µm, and [O I] 146 µm. Bottom row (left and middle): [N II] 122 µm and [O I] 63 µm. The bottom right panel shows all five spectra overlaid on the same vertical scale. The gray regions in the [N II] 122 µm and [O I] 63 µm spectra represent velocities between spectral windows. The spectra were obtained using observed and model visibilities, as described in the Appendix of Litke et al. (2019). Typical uncertainties are plotted in the upper left corners of the [N II] 205 µm, [C II] 158 µm, [O I] 145 µm, [N II] 122 µm, and [O I] 63 µm spectra.</figcaption></figure>

<table><caption>**Table 2**
SPT 0346–52 Continuum Lens Modeling Results</caption><thead><tr><th>Parameter</th><th>205 µm</th><th>158 µm</th><th>146 µm</th><th>122 µm</th><th>63 µm</th></tr></thead><tbody><tr><td colspan="6"><strong>Lens Parameters</strong></td></tr><tr><td>log mass [M<sub>☉</sub>]</td><td>11.46 ± 0.01</td><td>11.46 ± 0.01</td><td>11.46 ± 0.02</td><td>11.50 ± 0.01</td><td>11.47 ± 0.02</td></tr><tr><td>Ellipticity x-component, e<sub>x</sub></td><td>-0.17 ± 0.01</td><td>-0.17 ± 0.01</td><td>-0.18 ± 0.02</td><td>-0.16 ± 0.01</td><td>-0.16 ± 0.01</td></tr><tr><td>Ellipticity y-component, e<sub>y</sub></td><td>+0.43 ± 0.03</td><td>+0.41 ± 0.04</td><td>+0.39 ± 0.07</td><td>+0.22 ± 0.01</td><td>+0.39 ± 0.09</td></tr><tr><td>Ellipticity, e<sup>a,c</sup></td><td>0.45</td><td>0.45</td><td>0.42</td><td>0.25</td><td>0.45</td></tr><tr><td>Position angle, φ<sub>e</sub> (E of N)<sup>b,c</sup></td><td>69°</td><td>68°</td><td>65°</td><td>55°</td><td>68°</td></tr><tr><td>Shear x-component, γ<sub>x</sub></td><td>+0.06 ± 0.01</td><td>+0.06 ± 0.01</td><td>+0.09 ± 0.02</td><td>+0.15 ± 0.01</td><td>+0.07 ± 0.03</td></tr><tr><td>Shear y-component, γ<sub>y</sub></td><td>-0.10 ± 0.01</td><td>-0.11 ± 0.01</td><td>-0.10 ± 0.01</td><td>-0.08 ± 0.01</td><td>-0.10 ± 0.01</td></tr><tr><td>Shear amplitude, γ<sup>a,c</sup></td><td>0.12</td><td>0.12</td><td>0.13</td><td>0.17</td><td>0.12</td></tr><tr><td>Shear position angle, φ<sub>γ</sub> (E of N)<sup>b,c</sup></td><td>120°</td><td>120°</td><td>127°</td><td>115°</td><td>119°</td></tr><tr><td>Lens x-position, x</td><td>0''04 ± 0''01</td><td>0''08 ± 0''01</td><td>0''00 ± 0''01</td><td>0''01 ± 0''01</td><td>0''05 ± 0''01</td></tr><tr><td>Lens y-position, y</td><td>0''38 ± 0''01</td><td>0''32 ± 0''01</td><td>0''40 ± 0''01</td><td>0''39 ± 0''01</td><td>0''35 ± 0''01</td></tr><tr><td colspan="6"><strong>Source-plane Fluxes</strong></td></tr><tr><td>Continuum flux (mJy)</td><td>6.8 ± 1.1</td><td>11.9 ± 1.2</td><td>13.2 ± 1.3</td><td>24.9 ± 2.5</td><td>28.1 ± 2.8</td></tr><tr><td>Line luminosity (10<sup>8</sup> L<sub>☉</sub>)</td><td>1.2 ± 0.2</td><td>34 ± 5</td><td>3.2 ± 0.5</td><td>&lt;3.5<sup>d</sup></td><td>&lt;53<sup>d</sup></td></tr></tbody></table>

**Notes.**

was then applied to the line data corresponding to that continuum. For example, the 158 µm model was applied to the [C II] 158 µm line data, while the 205 µm model was applied to the [N II] 205 µm line data. The lines were integrated from -300 km s⁻¹ to +300 km s⁻¹ for the reconstructions. Table 2 gives the best-fit lensing parameters for all five models, as well as the source-plane continuum fluxes and line luminosities. The parameter covariance plot from the MCMC for all five continuum sets is shown in Figure 3.

## 3.1. Comparing Models

As seen in Figure 3, the five continuum models do not have identical lens parameters, though they are very similar in most parameters.<fnref n="15" /> In order to explore the effect of the differing

<footnote><sup>a</sup> α = $\sqrt{\alpha_x^2 + \alpha_y^2}$, where α = e or γ.</footnote>

<footnote><sup>b</sup> φ<sub>α</sub> = arctan(-α<sub>y</sub>/α<sub>α</sub>), where α = e or γ.</footnote>

<footnote><sup>c</sup> Derived from best-fit parameters.</footnote>

<footnote><sup>d</sup> 3σ upper limit from observations, corrected for lensing.</footnote>

<footnote n="15">¹⁵ The exception lies with the lens x- and y-positions. These positions are measured relative to the phase centers of the observations. The variations in the modeled positions of the lens from different data sets are consistent with the variations seen in the position of the astrometric test sources.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Triangle plot with the model lens parameters for SPT 0346-52 computed for different continuum wavelengths. Red: 205 µm continuum model. Orange: 158 µm continuum model. Green: 146 µm continuum model. Blue: 122 µm continuum model. Purple: 63 µm continuum model. Here *M* is the lens mass enclosed within 10 kpc and is measured in *M*⊙, *e*<sub>x</sub> and *e*<sub>y</sub> are the two components of the lens galaxy's ellipticity, γ<sub>x</sub> and γ<sub>y</sub> are the two components of the shear, and x and y are the offset of the lens center from the phase center in arcseconds.</figcaption></figure>

models on the source-plane reconstructions, we applied each
model to each of the other continuum data sets. We then
reconstructed the source-plane continuum emission. The
derived continuum fluxes were consistent, independent of the
model used. The differences between the best-fit models are

within the errors of the MCMC fit and most likely result from
degeneracies between the ellipticity and shear parameters.
These differences do not affect our results in the source plane.

To find the uncertainty on the flux in each pixel, we created
500 sets of random visibilities from the distribution of
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Reconstructed source-plane emission. From left to right: 205, 158, 146, 122, and 63 µm. Top: continuum maps of the source plane. Color bar units are mJy. Bottom row: moment zero maps (velocity-integrated line flux). Color bar units are mJy km s⁻¹. The ellipses represent the typical resolution where there is emission in the source plane. The black crosses mark the center of the source plane. A representative caustic is shown as black lines in the top right panel (the 63 µm continuum map). All source-plane images are 8 kpc on each side.</figcaption></figure>

uncertainties in the visibilities. We then reconstructed these 500 random noise data sets and took the standard deviation in each pixel. The total error used is the standard deviation of the random noise reconstructions added in quadrature with 10% of the flux.

In order to determine the effective resolution of the reconstructed maps, we follow the method used by Litke et al. (2019). We define the effective resolution as the inferred source-plane size when reconstructing a lensed point source at that source-plane position. This effective resolution will vary depending on the position of the source relative to the caustic in the source plane, as well as the signal-to-noise ratio of the input data. For each set of visibilities, we create a point source with the flux and position of the emission at that wavelength. We then apply the corresponding lens model in Table 2 to these point sources to create a lensed set of visibilities. Next, we make source-plane reconstructions of these lensed point sources using RIPPLES. Finally, we fit a 2D Gaussian to the reconstructed image to find the effective resolution.

It is simpler to compare the different line and continuum reconstructions if they have comparable resolutions. Therefore, we tapered the visibilities in each data set to 300kλ before performing the reconstructions. The resolutions are shown as colored ellipses in the lower left corner of the continuum and moment zero maps in Figure 4. The effective resolution is ~0.08 × 0.15, which gives us ~700 pc resolution. The models, residuals, and error maps are shown in Figure 14 in Appendix A.

## 3.2. Lens Modeling Results

Figure 4 shows the reconstructed continuum maps and the moment zero (integrated flux) maps of the line reconstructions. The continuum maps mostly show similar morphologies but with different fluxes. As expected, the 63 µm continuum emission is the brightest, and the 205 µm emission is the weakest. We see differing morphologies in the line emission. The [C II] 158 µm and [N II] 205 µm lines show their brightest

emission offset from the [O I] 146 µm emission. Because the Earth's atmosphere limits access to these lines except at the highest redshifts, the sources with spatially resolved maps of FIR fine-structure lines are mostly very nearby galaxies. Parkin et al. (2013) and Herrera-Camus et al. (2018) found in M51 and NGC 1068, respectively, that the [C II] 158 µm and [N II] lines had similar morphologies. A similar offset between [O I] 146 µm and [C II] 158 µm emission was seen by Herrera-Camus et al. (2018) in NGC 1068. In their maps, the brightest [O I] 146 µm emission was associated with the central AGN in NGC 1068, while the [C II] 158 µm emission was associated with the peak in CO emission. Parkin et al. (2013) also found that the [O I] 63 µm and [O I] 146 µm emission peaked in the center of M51, most likely associated with the central AGN, offset from where the [C II] 158 µm and [N II] 122 µm emission was strongest. In the SMC, the [O I] 63 µm emission was associated with Hα emission and therefore recent massive star formation (Jameson et al. 2018). As seen in these other systems, the [O I] 146 µm and [O I] 63 µm emission has been associated with regions where one would expect dust heating, whether the dust is being heated by star formation or a central AGN. For SPT 0346-52, the [O I] 146 µm emission is concentrated closer to where the dust continuum is strongest. As any potential AGN contribution to the IR emission in SPT 0346-52 is negligible (Ma et al. 2016), the [O I] 146 µm emission appears to coincide with the most intense star formation.

Litke et al. (2019) determined that the intense star formation in SPT 0346-52 was driven by a major merger of two components. These components are centered at −310 and +160 km s⁻¹ and have similar widths (~300 km s⁻¹). The center of these two components lies near −75 km s⁻¹. The merger status was determined using position–velocity diagrams of the [C II] 158 µm emission, which is the highest signal-to-noise ratio line of those explored here. Litke et al. (2019) used higher-resolution reconstructions than are shown in this work, where we decrease the resolution of the [C II] 158 µm
reconstruction to match that of the lowest-resolution lines. More recently, Jones et al. (2019) claimed that SPT 0346–52 is better described as a disk galaxy with a molecular outflow, based on visual inspection of the image-plane structure. However, the candidate H₂O outflow Jones et al. (2019) detected is kinematically similar to the blueshifted component found by Litke et al. (2019). The mass-loading factor of the possible outflow in SPT 0346–52 is well below unity, unlike other outflows with mass-loading factors near or greater than unity. Spilker et al. (2020) did not see broad wings in [C II] 158 µm that would be indicative of outflows in any DSFGs with confirmed molecular outflows traced by blueshifted OH absorption. Thus, SPT 0346–52 is unlikely to host outflow activity. The disk structure described by Jones et al. (2019) may also have resulted from a recent major merger (e.g., Hopkins et al. 2009).

## 4. CLOUDY Modeling

In order to understand the physical conditions that explain our observations, we turn to the photoionization code CLOUDY (version 17.01; Ferland et al. 2017). CLOUDY simulates the microphysics within a cloud of gas and dust that is heated by a central source. It predicts the physical conditions throughout the cloud, including temperatures, densities, and metallicities, while also computing a predicted observed spectrum.

### 4.1. CLOUDY Parameters

We model our system as an open, or slablike, geometry. We adopt an inner radius (the distance between the central heating source and the inner face of the cloud) of 100 pc. This distance is chosen to be smaller than the size of an individual pixel.

We use the ISM gas-phase elemental abundances and grain-size distributions included with CLOUDY, which represent the average warm and cold phase abundances of the ISM in the Milky Way (Cowie & Songaila 1986; Savage & Sembach 1996). We also include polycyclic aromatic hydrocarbons (PAHs) in the simulation. Small grains such as PAHs are an important contributor of grain heating mechanisms and far-UV (FUV) radiative transfer effects (Hollenbach & Tielens 1999). For the equation of state of the ISM, we assume constant pressure. The components balanced to achieve constant-pressure gas are

$$P_{\text{tot}} = P_{\text{gas}} + P_{\text{turb}} + P_{\text{lines}} + \Delta P_{\text{rad}}, \quad (2)$$

where $P_{\text{gas}}$ is the thermal gas pressure, $P_{\text{turb}}$ is the turbulent pressure, $P_{\text{lines}}$ is the radiation pressure due to trapped emission lines, and $\Delta P_{\text{rad}}$ is the pressure from the attenuation of the incident radiation field. In order to simplify the model, we do not include a magnetic field component. Following Cormier et al. (2019), we assume a constant microturbulent velocity of 1.5 km s⁻¹. This is consistent with the microturbulent velocities of individual PDRs (Tielens & Hollenbach 1985; Kaufman et al. 1999).

The gas cloud is heated by a single-burst stellar population. The starburst spectral energy distribution (SED) was compiled by Byler et al. (2017) for use in CLOUDY using the Flexible Stellar Population Synthesis code (Conroy et al. 2009; Conroy & Gunn 2010). We use the ionizing spectrum from Byler et al. (2017) produced by the MESA Isochrones and Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016). The MIST stellar evolution tracks differ from other models in that they include

stellar rotation, which results in harder ionizing spectra and higher luminosities (Byler et al. 2017). Byler et al. (2017) compared nebular emission ionized by MIST models to Padova (for low-mass stars; Bertelli et al. 1994; Girardi et al. 2000; Marigo et al. 2008) and Geneva (for high-mass stars; Schaller et al. 1992; Meynet & Maeder 2000) evolutionary tracks (Levesque et al. 2010) and found that the MIST models can match the observed line ratios better as the starburst ages past a few megayears. We fix the stellar metallicity to log Z/Z☉ = 0 and the age to the stellar age calculated by Ma et al. (2015) using SED fitting, ~30 Myr. Allowing the stellar metallicity to scale with the gas metallicity did not change our results. The input starburst SED sets the shape of the ionizing spectrum in CLOUDY. The intensity of the ionizing radiation is determined by the ionization parameter, described below.

We also include the cosmic microwave background (CMB) at z = 5.7 and cosmic rays to contribute to the gas heating. The CMB spectrum was then subtracted from the modeled continuum SED before comparing to the observations, which are measurements of excess above the CMB spectrum.

We create a grid of models varying the ionization parameter (U), hydrogen density ($n_H$) at the face of the cloud, as gas metallicity (Z). The ionization parameter is defined as the ratio of hydrogen-ionizing photons to the total hydrogen density, or, more specifically,

$$U \equiv \frac{Q(\text{H})}{4\pi r_0^2 n_H c}, \quad (3)$$

where $r_0$ is the distance between the center of the starburst and the inner surface of the cloud (100 pc), $n_H$ is the total hydrogen density, Q(H) is the number of hydrogen-ionizing photons, and c is the speed of light. We vary log U from -4.5 to -0.5 in steps of 0.25. The total hydrogen density includes molecular, atomic, and ionized hydrogen components. The $n_H$ is defined at the inner face of the cloud for the grid values 0.5 ≤ log $n_H$ ≤ 3.5 in steps of 0.25. We vary the gas metallicity from -2.0 ≤ log Z/Z☉ ≤ 2.0 in steps of 0.25. We stop the CLOUDY simulation at visual extinction $A_V$ = 100, following Abel et al. (2009).

To find the best-fit model, we compare the CLOUDY outputs to a combination of continuum ratios (log 63 µm/122 µm, log 63 µm/146 µm, log 63 µm/158 µm, log 63 µm/205 µm), line ratios (log[C II] 158 µm/[N II] 205 µm, log[C II] 158 µm/[O I] 146 µm), and a line-to-continuum ratio (log[O I] 146 µm/ 146 µm) to constrain the relative contributions of gas and dust. The best-fit model is chosen to be the model with the lowest reduced χ² value, $\chi_r^2$. Continuum values are in mJy, and line values are in L☉. The [O I] 146 µm/146 µm has units of 10⁷ L☉ mJy⁻¹. The factor of 10⁷ is used so the line-to-continuum ratio is the same order of magnitude as the continuum and line ratios. Table 6 in the Appendix lists the values used to compare to the CLOUDY models. We consider both spatially resolved and galaxy-integrated emission models.

We tested whether the inclusion of the continuum ratios affected our results. When the continuum ratios were not included, the line ratios were still well fit, but the continuum ratios from CLOUDY did not match our observations. The ionization parameter was the variable most sensitive to the inclusion of the continuum ratios. This is consistent with previous CLOUDY modeling, where Abel et al. (2009) found that the 60 µm/100 µm continuum ratio was strongly dependent on the ionization parameter.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of observed continuum and line ratios to CLOUDY best-fit galaxy-integrated ratios. From left to right: log 63 µm/122 µm, log 63 µm/146 µm, log 63 µm/158 µm, log 63 µm/205 µm, log [O I] 146 µm/146 µm, log[C II] 158 µm/[N II] 205 µm, and log[C II] 158 µm/[O I] 146 µm. Black circles represent the observed ratios and associated uncertainties. Purple squares represent the best-fit ratios from the CLOUDY modeling and associated uncertainties.</figcaption></figure>

## 4.2. Cloudy Modeling Results

Figure 5 shows the observed continuum and line ratios with their associated errors, as well as the ratios calculated for the best-fit CLOUDY model for SPT 0346–52. Table 3 lists the best-fit starburst age, *U*, *n*<sub>H</sub>, and *Z* values, as well as $\chi_r^2$, for the best-fit models for the global fits shown in Figure 5. Overall, the best-fit models agree with the observed global line and continuum ratios in SPT 0346–52.

It is a crude approximation to model an entire galaxy as a single PDR. Indeed, Katz et al. (2019) found in their simulations of emission lines in high-z galaxies that there can be very large ranges in metallicity and ionization parameter across a single galaxy. In the case of SPT 0346–52, we have spatially resolved information about the line and continuum emission and can explore the distribution of properties across the galaxy by fitting our models on a pixel-by-pixel basis. We use the same technique as for the galaxy-integrated fits for each pixel. Figure 6 shows the range of values observed and the range of the best-fit models, as well as the mean values, for each set of continuum and line ratios, while Figure 7 shows the best-fit *U*, *n*<sub>H</sub>, and *Z* values in each pixel. For maps of the observed continuum and line ratios, as well as the best-fit model ratios, see Figure 15 in the Appendix.

As with the fits to the global line and continuum ratios, the range of ratios in the pixelated best-fit models mostly agree with the range of ratios observed. The modeled [C II] 158 µm/[N II] 205 µm ratios have higher maximum values than the observed ratios, though their average values are more similar. We also see a higher average [C II] 158 µm/[O I] 146 µm model ratio than observed, as we did with the galaxy-integrated fit. This may result from the CLOUDY model not probing far enough from the ionizing source to fully recover line emission from the neutral gas component of the ISM.

The best-fit ionization parameters are around $log U = -2.75$ and relatively uniform. These values are comparable to those found for the DSFG SPT 0418–47 (using [N II] 205, [C II] 158, [O I] 146, [N II] 122, and [O III] 88; De Breuck et al. 2019). A popular technique is to report the intensity of the FUV, *G*, relative to the interstellar radiation field, $G_0$, instead of the ionization parameter. To compare our *U* values to models using *G*, we turn to Figure 1 of Abel et al. (2009); for a starburst SED, $log U \approx log G - 6$, and $log U = -3.25$ corresponds to $G = 10^{3.25}G_0$. Rybak et al. (2019, 2020) found $log G/G_0 \approx 4$ for $z \sim 3$ DSFGs and the $z = 6$ DSFG G09.83808, respectively,

<table><caption>**Table 3**

Best-fit Global CLOUDY Models<sup>a</sup></caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>log U</td><td>-2.75<sup>+1.3</sup><sub>-0.1</sub></td></tr><tr><td>log n<sub>H</sub> [cm<sup>-3</sup>]</td><td>1.75<sup>+0.1</sup><sub>-1.1</sub></td></tr><tr><td>log Z/Z<sub>☉</sub></td><td>0.75<sup>+0.5</sup><sub>-0.1</sub></td></tr><tr><td>χ<sub>r</sub><sup>2</sup></td><td>6.6</td></tr></tbody></table>

**Note.**

<sup>a</sup> The likelihood distribution is calculated by summing $e^{-\chi_r^2}$. The best-fit value is the peak of the likelihood distribution, and the uncertainties are where the likelihood distribution is within $1\sigma$ of the peak.

while Novak et al. (2019) found $log G/G_0 > 3$ for the $z = 7.5$ quasar host galaxy J1342+0928. Compared to other high-z DSFGs, SPT 0346–52 has a lower ionization parameter and FUV field strength.

We find densities around $log n_H \sim 2$ throughout this system. This is lower than the inferred densities of the DSFGs SPT 0418–47 at $z = 4.225$ ($log n \sim 4.3$; De Breuck et al. 2019) and G09.83808 at $z = 6.027$ ($log n \sim 4$; Rybak et al. 2020).

From the CLOUDY modeling, SPT 0346–52 appears to have a supersolar metallicity ($log Z/Z_☉ = 0.75$). This is higher than other high-z sources that have been found to have metallicities near solar ($log Z/Z_☉ \sim 0.1$ for quasar host galaxy J1342+0928 at $z = 7.54$ and $-0.5 < log Z/Z_☉ < 0.1$ for DSFG SPT 0418–47) using [O III] 88/[N II] 122 (De Breuck et al. 2019; Novak et al. 2019). However, the [O III] 88/[N II] 122 metallicity diagnostic is highly dependent on the ionization parameter (Pereira-Santaella et al. 2017). Using the mass–metallicity relation of elliptical galaxies, Tan et al. (2014) found the metallicity of DSFGs in the protocluster GN20 ($z = 4.05$) to be $log Z/Z_☉ \sim 0.5 \pm 0.2$. The gas-to-dust ratio, $\delta_{GDR}$, has been observed to be approximately inversely proportional to the metallicity in galaxies (e.g., Magdis et al. 2012; Leroy et al. 2011). Leroy et al. (2011) determined the relation between the gas-to-dust ratio and metallicity to be $log \delta = (9.4 \pm 1.1) - (0.85 \pm 0.13) \times (12 + \log([O/H]))$. Using the gas and dust masses from Aravena et al. (2016), $\delta_{GDR} = 41 \pm 13$ in SPT 0346–52, giving us $log Z/Z_☉ = 0.7 \pm 0.3$, consistent with the CLOUDY value.

## 5. Discussion

In this section, we explore various diagnostics of the different phases of the ISM in SPT 0346–52. We begin with dust emission, then move on to line deficits and the electron density. We then explore the prevalence of ionized versus neutral gas and the nondetection of [O I] 63. Finally, we look into several gas mass estimates for different phases of the ISM. For reference, Table 4 lists the transitions explored and their excitation properties. While SPT 0346–52 has a negligible AGN contribution to its $L_{FIR}$, it is possible that there is a highly dust-obscured AGN that prevents X-ray emission from being observed (Ma et al. 2016). The DSFGs can contain AGNs; for example, Wang et al. (2013) found that $17_{-6}^{+16}$% of DSFGs in the ALMA LABOCA E-CDF-S Submm Survey sample contain AGNs. Also, SPT 0346–52 may evolve to contain an AGN (e.g., Toft et al. 2014). We therefore consider both DSFGs and quasar host galaxies for comparison.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison of observed continuum and line ratios to CLOUDY best-fit pixelated ratios. From left to right: log 63 µm/122 µm, log 63 µm/146 µm, log 63 µm/158 µm, log 63 µm/205 µm, log[O I] 146 µm/146 µm, log[C II] 158 µm/[N II] 205 µm, and log[C II] 158 µm/[O I] 146 µm. Black circles represent the means of the observed pixel ratios, and the error bars represent the maximum and minimum observed ratios. Purple squares represent the means of the best-fit pixel ratios from the CLOUDY modeling, with the maximum and minimum best-fit ratios.</figcaption></figure>

## 5.1. Dust Temperatures

We characterize the FIR continuum and dust temperatures throughout the source plane using a modified blackbody function. We use the form from Spilker et al. (2016),

$$S_{\nu_r} \propto (B_{\nu_r}(T_D) - B_{\nu_r}(T_{\text{CMB}}))(1 - e^{-\tau_{\nu_r}}). \quad (4)$$

Here $B_{\nu_r}(T)$ is the Planck function at rest-frame frequency $\nu_r$ and temperature $T$. The blackbody is modified by the dust optical depth, $\tau_{\nu_r}$, which, at long wavelengths, can be parameterized by $\tau_{\nu} = (\nu/\nu_0)^\beta = (\lambda_0/\lambda)^\beta$. In this parameterization, the optical depth reaches unity at wavelength $\lambda_0$, which, together with $T_D$, determines the peak wavelength and the width of the peak of the dust emission. The slope of the Rayleigh-Jeans tail of the dust emission is $\beta$. Typically, $\beta \sim 1.5-2$ and $\lambda_0 \sim 100-200$ µm in the rest frame (e.g., Casey et al. 2014).

This modified blackbody is used to fit for the dust temperature in each pixel of the source-plane reconstruction, as well as the global dust temperature. Five continuum bands (63, 122, 146, 158, and 205 µm) are used to fit the SEDs. Because of the limited number of photometric points available for these fits, we follow Greve et al. (2012) by fixing $\beta = 2.0$ and $\lambda_0 = 100$ µm. There are two free parameters in each SED fit: the dust temperature, $T_D$, and the normalization.

The top panel of Figure 8 shows the global SED fit (black line) to the continuum (colored stars). For SPT 0346-52, we calculate $T_D = 48.3 \pm 4.0$ K. The gray lines show the best-fit SEDs for the individual pixels. The best-fit dust temperatures are shown in the bottom middle panel of Figure 8, while the right panel shows the error on the best-fit dust temperatures in each pixel. The left panel shows the FIR surface density ($\Sigma_{\text{FIR}}$). The $L_{\text{FIR}}$ is calculated here as the integral of the best-fit modified blackbody SED from 42.5 to 122.5 µm (Helou et al. 1988). To get $\Sigma_{\text{FIR}}$, we divide by the source-plane pixel area in square kiloparsecs. Dust temperature values in individual pixels have a mean $T_D = 56$ K.

Dust temperatures can vary significantly depending on the fitting form used (e.g., Hayward et al. 2012). Casey et al. (2014) showed how, for a given $\lambda_{\text{peak}}$, the dust temperature can vary by up to 40 K, depending on what assumptions are made about the opacity, as well as the value of $\beta$ used. The dust temperatures we calculate are lower than those calculated by Reuter et al. (2020) using Equation (4) ($T_D = 79 \pm 15$ K for SPT 0346-52). However, this temperature difference is a result of the fitting procedures used in this work and by

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Best-fit spatially resolved CLOUDY parameters. Left: log U. Middle: log nH/cm⁻³. Right: log Z/Z☉.</figcaption></figure>

Reuter et al. (2020). Reuter et al. (2020) used the relation between $T_D$ and $\lambda_0$ found by Spilker et al. (2016). Using this relation tends to increase the fitted dust temperature by ~20% (Reuter et al. 2020). Jones et al. (2020) also calculated $T_D = 79 \pm 0.5$ K for SPT 0346-52 using a modified blackbody distribution with effects from the CMB. The small uncertainty is largely due to the unrealistically small photometric errors claimed for the data used in the fit, many of which are well below 1%. Using the model from Jones et al. (2020) and the data from this paper, we calculate $T_D = 71 \pm 3$ K. Apostolovski et al. (2019) calculated $T_D = 29 \pm 1$ K for SPT 0346-52 by using a radiative transfer model. The significant difference from the results of the continuum SED fits presented above undoubtedly results from the difference in methodology and the impact of including CO excitation as a constraint on the dust temperature.

## 5.2. Line Deficits

A commonly observed phenomenon is the so-called “[C II] deficit,” where $L_{[C II]}/L_{\text{FIR}}$ falls as $\Sigma_{\text{FIR}}$ increases (e.g., Luhman et al. 1998; Sargsyan et al. 2012; Farrah et al. 2013; Oteo et al. 2016; Spilker et al. 2016; Gullberg et al. 2018). Here we explore possible deficits in other FIR lines that trace different components of the ISM. Figure 9 shows the line-to-FIR luminosity ratios versus $\Sigma_{\text{FIR}}$ for the five lines in this work. In SPT 0346-52, we see deficits, i.e., apparent trends with $\Sigma_{\text{FIR}}$, in all five lines. Additionally, we see spatially resolved deficits (apparent trends with $\Sigma_{\text{FIR}}$ in individual pixels) in all three detected lines.

Graciá-Carpio et al. (2011) found deficits in [C II] 158 µm, [N II] 122 µm, [O I] 146 µm, and [O I] 63 µm, concluding that line deficits occurred in both ionized and neutral gas in galaxies with a variety of redshifts and optical classifications. They explained the deficits as resulting from an increase in the ionization parameter at $L_{\text{FIR}}/M_{\text{H}_2} > 80$ L☉/M☉ as highly compressed, more efficient star formation leads to enhanced ionization parameters. Zhao et al. (2016) found that there was only a [N II] 205 deficit in LIRGs that had warm ($f_{70}/f_{160} > 0.6$) colors. SPT 0346-52 falls into this warm-color regime and does indeed exhibit an [N II] 205 deficit. On the other hand, several individual sources, ranging from $z = 1.5$ main-sequence galaxies to $z = 6$ quasar host galaxies, did not show lower [O I] 146 µm/LFIR (Li et al. 2020) and [O I] 63 µm/LFIR ratios (Sturm et al. 2010; Coppin et al. 2012; Wagg et al. 2020). In the spatially resolved SHINING galaxies, a sample of nearby galaxies that includes star-forming galaxies, AGN host galaxies, and LIRGs, Herrera-Camus et al. (2018) found the trend of decreasing line/LFIR strongest for singly ionized lines like [N II] 122 µm and [C II] 158 µm and weakest for neutral [O I] 146 µm and [O I] 63 µm.
<table><caption>**Table 4**
Targeted Fine-structure Lines</caption><thead><tr><th>Line</th><th>Transition<sup>a</sup></th><th>ν<sub>0</sub> (GHz)<sup>a</sup></th><th>E<sub>ion</sub><sup>b</sup> (eV)</th><th>T<sub>e</sub> (K)<sup>b</sup></th><th>n<sub>crit,H</sub> (cm<sup>-3</sup>)<sup>b</sup></th><th>n<sub>crit,e</sub> (cm<sup>-3</sup>)<sup>a</sup></th><th>A (s<sup>-1</sup>)<sup>a</sup></th><th>g<sub>u</sub><sup>a</sup></th></tr></thead><tbody><tr><td>[N II] 205 µm</td><td><sup>3</sup>P<sub>1</sub>-<sup>3</sup>P<sub>0</sub></td><td>1461</td><td>14.53</td><td>70</td><td>1.76 × 10<sup>2</sup></td><td>48</td><td>2.1 × 10<sup>-6</sup></td><td>3</td></tr><tr><td>[C II] 158 µm</td><td><sup>2</sup>P<sub>3/2</sub>-<sup>2</sup>P<sub>1/2</sub></td><td>1901</td><td>11.26</td><td>91</td><td>4.93 × 10<sup>1</sup></td><td>50</td><td>2.1 × 10<sup>-6</sup></td><td>4</td></tr><tr><td>[O I] 146 µm</td><td><sup>3</sup>P<sub>0</sub>-<sup>3</sup>P<sub>1</sub></td><td>2060</td><td>...</td><td>327</td><td>7.65 × 10<sup>3</sup></td><td>...</td><td>1.7 × 10<sup>-5</sup></td><td>1</td></tr><tr><td>[N II] 122 µm</td><td><sup>3</sup>P<sub>2</sub>-<sup>3</sup>P<sub>1</sub></td><td>2459</td><td>14.53</td><td>188</td><td>3.86 × 10<sup>2</sup></td><td>310</td><td>7.5 × 10<sup>-6</sup></td><td>5</td></tr><tr><td>[O I] 63 µm</td><td><sup>3</sup>P<sub>1</sub>-<sup>3</sup>P<sub>2</sub></td><td>4745</td><td>...</td><td>228</td><td>3.14 × 10<sup>4</sup></td><td>...</td><td>9.0 × 10<sup>-5</sup></td><td>3</td></tr></tbody></table>

Notes. The first two columns list the targeted line and the fine-structure transition that emits that line. Here $\nu_0$ is the emitted frequency of the line, $E_{ion}$ is the ionization energy needed to remove an electron, $T_e$ is the excitation temperature needed to populate the transition level, $n_{crit,H}$ is the critical density for collisions with hydrogen at $T = 100$ K, $n_{crit,e}$ is the critical density for collisions with electrons at $T = 10,000$ K, $A$ is the Einstein A coefficient, and $g_u$ is the statistical weight of the upper level.

$^a$ Stacey (2011).

$^b$ Cormier et al. (2019).

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

**Figure 8.** Top: modified blackbody fits to galaxy-integrated (black) and pixel-by-pixel (gray) continuum flux values. The colored stars are the measured global continuum flux values. The best-fit SED has a temperature of $T_D = 48.3$ K. Bottom left: $\Sigma_{FIR}$, in $L_{\odot} \text{ kpc}^{-2}$ across the source. Bottom middle: dust temperature throughout the system. Bottom right: error in dust temperature. The mapped $T_D$ values are the temperatures from the best-fit modified blackbody SEDs in each pixel. The $\Sigma_{FIR}$ is calculated by integrating the best-fit modified blackbody SED in each pixel from 42.5 to 122.5 µm and dividing by the pixel area.

### 5.3. [C II] 158/[N II] 205

Carbon has an ionization potential of 11.26 eV, which is slightly lower than that of hydrogen. This makes interpretation of [C II] 158 µm emission difficult because it can originate from both ionized and neutral regions of the ISM. By comparing the [C II] 158 µm to the [N II] 205 µm emission, which arises only in ionized regions, we can infer how the [C II] 158 µm emission is divided between neutral and ionized gas.

We can calculate the fraction of [C II] 158 µm emission originating from neutral gas by comparing the observed [C II]

158 µm/[N II] 205 µm ratio to the expected ratio. From Abdullah et al. (2017),

$$R_{\text{ion}} = \frac{I_{[C \ II]}}{I_{[N \ II]}} = \frac{N_{[C \ II]u}}{N_{[N \ II]u}} \times \frac{E_{[C \ II]ul}}{E_{[N \ II]ul}} \times \frac{A_{[C \ II]ul}}{A_{[N \ II]ul}}, \quad (5)$$

where $R_{ion}$ is the expected line intensity ratio, $I$ is the expected line intensity of the transition, $N$ is the upper-level population, $E$ is the energy of the transition, and $A$ is the Einstein coefficient for that transition. This relation assumes that the ionic abundance ratio, [C II] 158 µm/[N II] 205 µm, is equal to the elemental abundance ratio, $C/N$. The $R_{ion}$ also depends on $n_e$. As seen in Table 4, [C II] 158 µm and [N II] 205 µm have similar critical electron densities (50 and 48 cm$^{-3}$). Therefore, if we compute the expected [C II] 158 µm/[N II] 205 µm ratio, there is little density dependence.

Croxall et al. (2017) calculated $R_{ion} = 4.0$ using collision rates from Tayal (2008; [C II] 158) and Tayal (2011; [N II] 205) and assuming Galactic gas-phase abundances for both elements, while Díaz-Santos et al. (2017) used $R_{ion} \approx 3.0 \pm 0.5$, based on photoionization models by Oberst et al. (2006). We adopt an intermediate value of $R_{ion} = 3.5$. With this expected [C II] 158 µm/[N II] 205 µm ratio, we can calculate the fraction of [C II] 158 µm from neutral gas, $f_{[C II], neutral}$, using

$$f_{[C \ II], \text{neutral}} = \frac{[C \ II] 158 \mu\text{m} - R_{\text{ion}} \times [N \ II] 205 \mu\text{m}}{[C \ II] 158 \mu\text{m}}. \quad (6)$$

For SPT 0346–52, we calculate $f_{[C II], neutral} = 0.84 \pm 0.04$. Nearby galaxies ranging from low-metallicity dwarf galaxies (Cormier et al. 2019) to star-forming galaxies (Herrera-Camus et al. 2018; Sutter et al. 2019) and (U)LIRGs (Díaz-Santos et al. 2017) have $f_{[C II], neutral} \sim 60\%–90\%$. As shown in Figure 10, SPT 0346–52 has a more comparable $f_{[C II], neutral}$ to what has been observed in other high-z sources and DSFGs ($f_{[C II], neutral} \sim 85\%$; e.g., Pavesi et al. 2018; Zhang et al. 2018; De Breuck et al. 2019; Li et al. 2020). This fraction has been observed to be higher in active star-forming regions (Herrera-Camus et al. 2018) and LIRGs with warmer $S_{63}/S_{158}$ colors (Díaz-Santos et al. 2017) like SPT 0346–52.

### 5.4. [C II] 158 µm/[O I] 146 µm

The [C II] 158 µm emission originates from both ionized and neutral gas, while [O I] 146 µm emission arises from only neutral regions. Therefore, more [O I] 146 µm emission would indicate the presence of more dense, neutral gas. The ratio of these lines has therefore been used in the literature as an
<figure><img src="image_12.png" /></figure>

**Figure 9.** Line luminosity/$L_{FIR}$ vs. $\Sigma_{FIR}$ for the observed lines. From top to bottom: [N II] 205 µm, [C II] 158 µm, [O I] 146 µm, [N II] 122 µm, and [O I] 63 µm. The right column is the same as the left column but focused on the regions where our pixels lie (the parameter space indicated by purple boxes in the left column). Colored diamonds are individual pixels, while colored stars are the global value for each line. For [N II] 122 µm and [O I] 63 µm, the upper limits are shown. Brown triangles are galaxies from SPT (Gullberg et al. 2015; De Breuck et al. 2019; Cunningham et al. 2020; Reuter et al. 2020). Gray dots are galaxy-integrated values from SHINING (Herrera-Camus et al. 2018), and black dots are from GOALS (Lutz et al. 2016; Díaz-Santos et al. 2017; Lu et al. 2017). Tan dots are a selection of high-z line detections: J1342+0928 (Venemans et al. 2017; Bañados et al. 2019; Novak et al. 2019), SPT 0311–58 (Marrone et al. 2018), PJ 231–20 (Pensabene et al. 2021), HFLS3 (Riechers 2013), PJ 308–21 (Decarli et al. 2019; Pensabene et al. 2021), G09.83808 (Zavala et al. 2018; Rybak et al. 2020), J2310+1855 (Shao et al. 2019; Li et al. 2020), BR 1202–0725 (Decarli et al. 2014; Lu et al. 2017; Lee et al. 2019, 2021), SMM J02399 (Weiß et al. 2007; Ivison et al. 2010; Ferkinhoff et al. 2011), Cloverleaf (Weiß et al. 2003; Ferkinhoff et al. 2011), SDP.11 (Lamarche et al. 2018), and MIPS J1428 (Iono et al. 2006; Hailey-Dunsheath et al. 2010).

indicator of the prevalence of dense gas (De Breuck et al. 2019; Li et al. 2020).

Figure 11 plots $L_{[C II]158\mu m}/L_{[O I]146\mu m}$ as a function of $\Sigma_{FIR}$ for SPT 0346–52, both integrated and spatially resolved, and values from the literature. Object SPT 0346–52 has a similar [C II] 158 µm/[O I] 146 µm ratio compared to galaxies in the SHINING sample (Herrera-Camus et al. 2018). It is higher than SPT 0418–47, a $z = 4.2$ lensed DSFG. De Breuck et al. (2019)

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

**Figure 10.** Fraction of [C II] 158 µm emission originating from neutral gas, $f_{[C II],neutral}$, calculated from the [C II] 158 µm/[N II] 205 µm ratio. Left: SPT 0346–52 (pink star) is compared to galaxies from the SPT (brown triangles; Gullberg et al. 2015; Cunningham et al. 2020; Reuter et al. 2020), high-z sources (tan dots; Decarli et al. 2014; Lu et al. 2017; Pensabene et al. 2021), and GOALS (black dots; Lutz et al. 2016; Díaz-Santos et al. 2017; Lu et al. 2017). Right: fraction of [C II] 158 µm originating from neutral gas in SPT 0346–52, with the galaxy-integrated value listed in the upper left corner ($f_{[C II],neutral} \sim 0.85$).

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

**Figure 11.** Shown is [C II] 158 µm/[O I] 146 µm. Left: $L_{[C II]158\mu m}/L_{[O I]146\mu m}$ vs. $\Sigma_{FIR}$. Pink diamonds are the individual pixels in SPT 0346–52, while the pink star represents SPT 0346–52. Comparison samples are taken from De Breuck et al. (2019; SPT 0418–47, brown triangle), Lee et al. (2021), Decarli et al. (2014), Novak et al. (2019; BR 1202–0725 and J1342+0928, tan dots), and Herrera-Camus et al. (2018; SHINING, gray dots). Right: mapped ratio of [C II] 158 µm/[O I] 146 µm. Lower $L_{[C II]158\mu m}/L_{[O I]146\mu m}$ values (more [O I] 146 µm) indicate more dense, neutral gas.

determined that SPT 0418–47 had a [C II] 158 µm/[O I] 146 µm ratio ~5× lower than local galaxies, leading them to conclude that the ISM in SPT 0418–47 is dominated by dense gas. Li et al. (2020) also found that [C II] 158 µm/[O I] 146 µm in their $z \sim 6$ quasar was comparable to the lowest values in ULIRGs, implying that SDSS J2310+1855 has warmer and denser gas compared to local galaxies. The higher [C II] 158 µm/[O I] 146 µm ratio in SPT 0346–52, implying a smaller dense gas component, is consistent with the lower hydrogen gas densities found using CLOUDY in Section 4.2.

## 5.5. Nondetection of [N II] 122 µm

The [N II] 122 µm is expected to be brighter than [N II] 205 µm. However, due to atmospheric $O_2$ at 368.5 GHz, just ~800 km s⁻¹ from the expected center of the [N II] 122 µm emission, [N II] 122 µm is not detected in SPT 0346–52. However, we can use the upper limit obtained in Section 2 to place constraints on the ISM conditions.

Nitrogen ions are only expected to be found in ionized regions of the ISM. In this regime, the [N II] 122 µm and [N II] 205 µm fine-structure lines would be excited mostly through collisions with electrons (Goldsmith et al. 2015). Therefore, the relative intensity of [N II] 122 µm compared to [N II] 205 µm will depend on the electron density of the ISM and can be used to calculate this density.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Theoretical and calculated electron densities for SPT 0346–52. The gray line is the theoretical relation between the [N II] 122 µm/[N II] 205 µm ratio and electron density. The purple star represents the global [N II] 122 µm/[N II] 205 µm ratio and corresponding $n_e$ in SPT0346-52. The black dashed line represents the median value from the GOALS sample, and the dotted lines are the maximum and minimum values (Díaz-Santos et al. 2017).</figcaption></figure>

We calculate the theoretical relation between the [N II] 122 µm/[N II] 205 µm flux ratio and the electron density, $n_e$, following Goldsmith et al. (2015) and using the collision rate coefficients from Tayal (2011). This relation is shown as the gray line in Figure 12. This figure also shows the upper limit to the ratios and corresponding $n_e$ values for SPT 0346–52.

In SPT 0346–52, we find $n_e < 32 \text{ cm}^{-3}$. The densities observed in SPT 0346–52 are lower than the densities calculated in other, comparable systems using the [N II] 122/[N II] 205 ratio. For example, Díaz-Santos et al. (2017) found a median $n_e = 41 \text{ cm}^{-3}$ for local LIRGs, with densities in the range $20 \text{ cm}^{-3} < n_e < 100 \text{ cm}^{-3}$ (black dashed and dotted lines in Figure 12), comparable to those found by Zhao et al. (2016), while De Breuck et al. (2019) calculated $n_e \sim 50 \text{ cm}^{-3}$ for the lensed DSFG SPT 0418–47. The nondetection of [N II] 122 indicates low densities in the ionized phase of the ISM.

## 5.6. Nondetection of [O I] 63 µm

The [C II] 158 µm and [O I] 63 µm transitions are both major coolants of the ISM. As one transitions to high-density ($n > 10^{3–4} \text{ cm}^{-3}$), high-temperature ($T > 10^4$ K), high-radiation ($G_0 > 10^3$) regimes, [O I] 63 becomes the dominant coolant over [C II] 158 (Tielens & Hollenbach 1985; Hollenbach et al. 1991). One would expect to find bright [O I] 63 µm emission comparable to or greater than the [C II] 158 µm emission in FIR-bright systems like starbursts, where the ISM is heated by strong FUV radiation and warmer, denser gas is expected (Tielens & Hollenbach 1985). In the first $z > 3$ detection of [O I] 63 µm, Rybak et al. (2020) found that [O I] 63 µm was ~four times brighter than [C II] 158 in the $z \sim 6$ DSFG, G09.83808. Given the intense star formation in SPT 0346–52, we could reasonably expect that [O I] 63 µm might be significantly more luminous than [C II] 158 µm. However, we are unable to detect this line in our observations, finding a luminosity ratio of $L_{[C II]158}/L_{[O I]63} > 0.7$. Our inability to place a tighter detection limit reflects the impact of atmospheric O$_2$ absorption centered at 715.4 GHz, just ~1000 km s$^{-1}$ from the redshifted [O I] 63 µm line center.

The [O I] 63 µm emission may be intrinsically weak. Spinoglio & Malkan (1992) also found lower [O I] 63 µm intensities with their CLOUDY modeling of starburst regions. In addition, Abel et al. (2009) found that low values of the

ionization parameter, $U$, were associated with a lower [O I] 63 µm/[C II] 158 µm ratio in their CLOUDY models of ULIRGs. For SPT 0346–52, [O I] 63 µm/[C II] 158 µm < 1.5, which is lower than [O I] 63 µm/[C II] 158 µm ~ 4 in the $z = 6$ DSFG G09.83808 (Rybak et al. 2020). Rybak et al. (2020) found $G = 10^4 G_0$ (corresponding to $\log U \approx -2$), which is higher than the values found in SPT 0346–52. Our CLOUDY modeling of SPT 0346–52 indicates lower ionization parameters, consistent with the lower [O I] 63 µm/[C II] 158 µm ratio.

Often, [O I] 63 µm is optically thick (Tielens & Hollenbach 1985; Kaufman et al. 1999; Liseau et al. 2006). It is also easily self-absorbed; small amounts of cold foreground gas can absorb [O I] 63 µm while leaving [O I] 146 µm unaffected (Liseau et al. 2006). This effect has been measured to reduce the [O I] 63 µm emission by factors of $1.3 \pm 1.8$ (Kramer et al. 2020) up to $2.9 \pm 1.6$ (Liseau et al. 2006). If [O I] 63 µm/[O I] 146 µm < 10, [O I] 63 µm is likely self-absorbed (Tielens & Hollenbach 1985; Cormier et al. 2015; Díaz-Santos et al. 2017). Taking the $3\sigma$ upper limit and correcting for the lensing magnification, [O I] 63 µm/[O I] 146 µm < 14 in SPT 0346–52. The [O I] 63 µm/[O I] 146 µm intensity ratio may be <10, so [O I] 63 µm may be self-absorbed.

The [O I] 63 µm/[C II] 158 µm ratio may also be influenced by the presence of an AGN. With its high critical density (Table 4), [O I] 63 µm is produced primarily in dense, neutral gas. In models of PDRs and X-ray-dominated regions (XDRs), the [O I] 63 µm/[C II] 158 µm ratio is higher in XDRs (which are expected near AGNs) than in PDRs (Maloney et al. 1996; Meijerink et al. 2007). Ma et al. (2016) found no evidence of AGN activity in this source, so we do not expect an AGN-based enhancement of [O I] 63 µm emission.

## 5.7. Gas Mass Estimates

In this section, we estimate the ionized, neutral, and molecular gas masses using the various fine-structure lines observed in SPT 0346–52. For ease of comparison, and because of the higher metallicity expected by the gas-to-dust mass ratio in SPT 0346–52, the ionic abundances are assumed to be the same as the global abundances in H II regions from Savage & Sembach (1996) (i.e., $\chi([C II]) = C/H$, $\chi([N II]) = N/H$, and $\chi([O I]) = O/H$). Table 5 lists a summary of the different masses calculated using the various methods described below. Figure 13 shows the various masses calculated for SPT 0346–52 and several other high-redshift sources, normalized by the molecular gas mass calculated using CO.

### 5.7.1. Molecular Gas Mass from $\alpha_{[C II]}$

In their study of $z \sim 2$ main-sequence galaxies, Zanella et al. (2018) found that the [C II] 158 µm luminosity and the molecular gas mass of a galaxy were correlated. With a standard deviation of 0.3 dex, $\alpha_{[C II]} = 31 M_\odot/L_\odot$ is mostly independent of depletion time, metallicity, and redshift.

Using our observed [C II] 158 µm luminosity, we find a molecular gas mass of $M_{\text{gas},\alpha_{[C II]}} = 1.1 \pm 0.3 \times 10^{11} M_\odot$. This is between the molecular gas mass calculated by Aravena et al. (2016) using CO ($8.2 \pm 0.6 \times 10^{10} M_\odot$) and the molecular gas mass calculated by Apostolovski et al. (2019) using radiative transfer modeling ($3.9 \pm 2.2 \times 10^{11} M_\odot$).

Zanella et al. (2018) calibrated $\alpha_{[C II]}$ using main-sequence galaxies at $z \sim 2$. However, for DSFGs and quasar host galaxies, this method tends to result in higher molecular gas
<table><caption>**Table 5**
Summary of Mass Estimates</caption><thead><tr><th>Mass Type</th><th>Line or Reference</th><th>SPT 0346-52<br>(×10<sup>9</sup> M<sub>☉</sub>)</th></tr></thead><tbody><tr><td rowspan="2">Dust</td><td>a</td><td>2.1 ± 0.3</td></tr><tr><td>b</td><td>2.0 ± 0.6</td></tr><tr><td rowspan="2">Molecular</td><td>b</td><td>82 ± 6</td></tr><tr><td>c</td><td>390 ± 220</td></tr><tr><td>Stellar</td><td>d</td><td>&lt;310</td></tr><tr><td>Molecular</td><td>[C II] 158 µm</td><td>106 ± 15</td></tr><tr><td>PDR</td><td>[C II] 158 µm</td><td>4.1 ± 0.6</td></tr><tr><td rowspan="2">Neutral</td><td>[C II] 158 µm</td><td>24 ± 4</td></tr><tr><td>[O I] 146 µm</td><td>18 ± 3</td></tr><tr><td>Ionized</td><td>[N II] 205 µm</td><td>0.8 ± 0.1</td></tr></tbody></table>

**Notes.**

a Spilker et al. (2015).

b Aravena et al. (2016).

c Apostolovski et al. (2019).

d Ma et al. (2015).

masses than those obtained using CO, as shown by the green points in Figure 13. The molecular masses calculated using $\alpha_{[C \text{II}]}$ are also higher than the neutral gas masses calculated using the method in Section 5.7.3 in both this work and previous studies (Decarli et al. 2019; Novak et al. 2019).

### 5.7.2. Neutral Gas Mass from [C II] 158 µm PDR Modeling

Hailey-Dunsheath et al. (2010) calculated the atomic mass associated with PDRs in the $z = 1.3$ hyperluminous starburst galaxy MIPS J1428 using the [C II] 158 µm luminosity and PDR models from Kaufman et al. (1999). Assuming that the [C II] 158 µm emission is optically thin and that a single temperature characterizes the [C II] 158 µm-emitting region, the PDR mass is given by

$$ \frac{M_{\text{PDR}}}{M_{\odot}} = 0.77 \left( \frac{f_{[C \text{II}],\text{neutral}} L_{[C \text{II}]}}{L_{\odot}} \right) \left( \frac{1.4 \times 10^{-4}}{\chi(C^{+})} \right) \times \frac{1 + 2 \exp\left(\frac{-91K}{T}\right) + \frac{n_{\text{crit}}}{n}}{2 \exp\left(\frac{-91K}{T}\right)}. \quad (7) $$

Following Hailey-Dunsheath et al. (2010), we assume that the gas temperature is the surface temperature of their modeled PDR, $T \approx 230$ K, $n = 10^{4.2} \text{ cm}^{-3}$, $n_{\text{crit}} = 2.7 \times 10^{3} \text{ cm}^{-3}$ (Launay & Roueff 1977), and $C^{+}/H = 1.4 \times 10^{-4}$ (Savage & Sembach 1996). We use our calculated value of $f_{[C \text{II}],\text{neutral}} = 0.84$. We calculate a PDR mass of $M_{\text{PDR}} = 4.1 \pm 0.6 \times 10^{9} M_{\odot}$ for SPT 0346–52. This is ~5% of the total gas mass from Aravena et al. (2016).

The PDR mass fraction calculated for SPT 0346–52 is much lower than that found by Hailey-Dunsheath et al. (2010) for MIPS J1428 (~55%). We also find a similarly lower PDR mass compared to the molecular gas mass in SPT 0346–52 than in HFLS3 (20%; Riechers 2013) and SDP.11 (23%; Lamarche et al. 2018). This method assumes a gas density $\log n/\text{cm}^{-3} = 4.2$, ~100× higher than the density found using CLOUDY (Section 4.2). If we instead use $\log n/\text{cm}^{-3} = 2$, we calculate an atomic PDR mass fraction of 58%, which is more comparable to other high-z sources. This discrepancy could

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Mass estimates for SPT 0346–52 and other high-redshift sources, normalized by the molecular gas mass from CO measurements, vs. $L_{\text{FIR}}$. The stars represent the mass estimates for SPT 0346–52. The masses for J1342 +0928 ($z \approx 7.5$; Venemans et al. 2017; Novak et al. 2019), PJ 308–21 ($z \approx 6.2$; Decarli et al. 2019; Pensabene et al. 2021), SDP.11 ($z \approx 1.8$; Lamarche et al. 2018), MIPS J1428 ($z \approx 1.3$; Iono et al. 2006; Hailey-Dunsheath et al. 2010), Cloverleaf ($z \approx 2.5$; Weiß et al. 2003; Ferkinhoff et al. 2011), J2310+1855 ($z \approx 6.0$; Shao et al. 2019; Li et al. 2020), HFLS3 ($z \approx 6.3$; Riechers 2013; Cooray et al. 2014), and SPT 0311–58 (Marrone et al. 2018; Jarugula et al. 2021, $z \approx 7$;) are reported in the literature. The dust and CO molecular masses for the SPT galaxies, including SPT 0346–52, are taken from Aravena et al. (2016) and Marrone et al. (2018), and the stellar masses are taken from Ma et al. (2015) and Marrone et al. (2018); the rest of the masses are calculated using the methods described in Section 5.7 using data from Gullberg et al. (2015), Cunningham et al. (2020), and Reuter et al. (2020). Tan represents dust masses, green is molecular gas estimates using [C II] 158 µm (Molecular$_{[C \text{II}]}$; see Section 5.7.1), silver is a PDR model-based gas mass estimate (Section 5.7.2), blues are neutral gas mass estimates (Section 5.7.3), pinks/reds are ionized gas mass estimates (Section 5.7.4), and gold represents stellar mass estimates. The black dashed line represents $M = M_{\text{CO}}$.</figcaption></figure>

explain the difference in atomic PDR masses between SPT 0346–52 and other high-z sources. Masses calculated with this method are shown as silver points in Figure 13.

### 5.7.3. Neutral Gas Masses from [C II] 158 µm and [O I] 146 µm

We next estimate the neutral gas mass in SPT 0346–52 using [C II] 158 µm and [O I] 146 µm. Based on the work by Weiß et al. (2005) and Li et al. (2020), the mass associated with a single transition can be calculated using

$$ M_x = m_x \frac{8\pi k\nu_0^2 Q(T_{\text{ex}})}{hcA g_u} e^{T_e/T_{\text{ex}}} L', \quad (8) $$

where $m_x$ is the mass of atom x, $\nu_0$ is the emission frequency, A is the Einstein coefficient, Q is the partition function, $g_u$ is the statistical weight of the upper level, k is the Boltzmann constant, h is the Planck constant, $T_e$ is the excitation temperature needed to populate the transition level from Table 4, $T_{\text{ex}}$ is the excitation temperature of the gas, and $L'$ is the line luminosity in K km s$^{-1}$ pc$^{-2}$. The Einstein A coefficients, upper-level statistical weights, and transition temperatures can be found in Table 4.

We adopt an excitation temperature of $T_{\text{ex}} = 100$ K. Based on PDR modeling by Meijerink et al. (2007), 100 K is the temperature of a PDR cloud near the outer regions where [C II] 158 µm is primarily emitted. The cloud will cool off to ~20 K
in the inner parts of the cloud as carbon transitions from [C II] to neutral C to CO. We also use this temperature for the [O I] 146 µm–based mass calculation. This method assumes that the [C II] 158 µm and [O I] 146 µm emission is optically thin.

For [C II] 158 µm, Equation (8) becomes

$$ \frac{M_{C+}}{M_{\odot}} = 3.21 \times 10^{-4} \frac{Q(T_{ex})}{4} e^{91/T_{ex}} L'_{[C II]158\mu m}, \quad (9) $$

where $Q(T_{ex}) = 2 + 4e^{-91/T_{ex}}$, and $L'_{[C II]158\mu m}$ is multiplied by $f_{[C II],neutral}$ so that only the neutral gas contribution to the [C II] 158 emission is included. Decarli et al. (2019) argued that the mass calculated using [C II] 158 µm is a lower limit, as it does not include nonionized carbon or effects from lower-metallicity gas, suppressed [C II] 158 µm emission from collisional de-excitation, and optical depth effects.

For [O I] 146 µm, we can calculate the oxygen mass using

$$ \frac{M_O}{M_{\odot}} = 6.20 \times 10^{-5} Q(T_{ex}) e^{329/T_{ex}} L'_{[O I]146\mu m}. \quad (10) $$

In this case, $Q(T_{ex}) = 5 + 3e^{-228/T_{ex}} + e^{-329/T_{ex}}$.

To get the neutral gas mass, we divide by the C/H or O/H abundance. Using the abundance from Savage & Sembach (1996; C/H = 3.98 × 10⁻⁴ and O/H = 5.89 × 10⁻⁴), we find neutral gas masses of $M_{neutral,[CII]158\mu m} = 2.4 \pm 0.4 \times 10^{10}$ and $M_{neutral,[OII]146\mu m} = 1.8 \pm 0.3 \times 10^{10} M_{\odot}$. The [C II] 158 µm and [O I] 146 µm neutral gas masses are approximately one-fourth of the total gas mass calculated by Aravena et al. (2016). In general, the neutral gas masses calculated using [C II] 158 µm, [O I] 146 µm, and [C I] 369 µm in SPT0346-52, DSFGs, and in quasar host galaxies in the literature are within a factor of ~10 of the molecular gas mass (blue points in Figure 13).

The mass calculated using [C II] 158 µm is ~6× higher than the mass calculated using [C II] 158 µm and PDR modeling in Section 5.7.2. While two different temperatures were used (230 K versus 100 K), this difference only changes the calculated masses by ~30% and does not fully account for the large discrepancy in the mass estimates. The mass estimate in Section 5.7.2 is based on the PDR model of a cloud illuminated on one side from Kaufman et al. (1999). As discussed in Section 4, this is a simple model compared to the complexity of a galaxy. On the other hand, the method used in this section assumes optically thin emission and that the lines are in local thermodynamic equilibrium. These assumptions may not be valid for all of the [C II] 158 µm emission in SPT 0346–52 and could also account for the discrepancy between the [C II] 158 µm neutral gas mass estimates. In addition, as discussed in Section 5.7.2, the discrepancy also arises from the density of the gas in the PDR model. Using $n = 10^2 \text{ cm}^{-3}$ from our CLOUDY modeling instead of $n = 10^{4.2} \text{ cm}^{-3}$ from Hailey-Dunsheath et al. (2010) results in consistent PDR masses using the different methods.

### 5.7.4. Ionized Gas Mass

Following the method Ferkinhoff et al. (2010) used for [O III] 88 µm and Ferkinhoff et al. (2011) adapted for [N II] 122 µm, we can calculate the minimum ionized gas mass required to produce the observed [N II] 205 µm emission. This method assumes that all of the nitrogen in the H II regions is singly ionized, and that the gas is in the high-temperature limit, as would be expected in active star-forming regions. With these assumptions, we can calculate the minimum ionized gas mass

using

$$ M([H \ II])_{205} = \frac{L_{[N \ II]205\mu m}}{\frac{g_u}{Q(T_{ex})} A_{205} h \nu_{205}} \frac{m_H}{\chi([N \ II])}. \quad (11) $$

Here h is the Planck constant, and $m_H$ is the mass of the hydrogen atom. In the high-temperature limit, the partition function of [N II] 205 is $Q(T_{ex}) \to 9$. The $g_u$, A, and $\nu_{205}$ can all be found in Table 4. Using the nitrogen abundance from Savage & Sembach (1996), N/H = 7.76 × 10⁻⁵, we find a minimum ionized gas mass of $M([H II])_{205} \le 7.5 \pm 1.0 \times 10^8 M_{\odot}$.

Ferkinhoff et al. (2011) found that $M_{[H II]} / M_{mol}$ is correlated with $\Sigma_{SFR}$, where more intensely star-forming galaxies have higher fractions of ionized gas. The minimum ionized gas mass in SPT 0346–52 is ~1% of the total molecular gas mass from Aravena et al. (2016). These values are lower than the ionized gas mass fractions determined for the Cloverleaf (~8%; Ferkinhoff et al. 2011) and J1342+0928 (~4%; Novak et al. 2019). As mentioned above, this method assumes that the gas is in the high-temperature limit. This is likely not the case for SPT 0346–52. At lower temperatures, more mass will be required to produce the observed [N II] 205 µm emission. The ionized gas masses calculated here are therefore lower limits to the total ionized gas mass.

As shown in Figure 13, the molecular phase is the most significant mass component in SPT 0346–52. There is ~4× more molecular gas than neutral gas and ~100× more molecular gas than ionized gas. This is in contrast to galaxies in the nearby universe; for example, the Milky Way has a molecular gas mass (~2 × 10⁹ $M_{\odot}$) very close to its ionized gas mass ($\gtrsim 1.6 \times 10^9 M_{\odot}$) and more than twice as much atomic gas compared to molecular gas (Ferrière 2001). The large molecular gas reservoir in SPT 0346–52 fuels the large star formation rate observed in this system.

## 6. Summary and Conclusions

In this work, we present ALMA Bands 6, 7, and 9 observations of [N II] 205 µm, [C II] 158 µm, [O I] 146 µm, and undetected [N II] 122 µm and [O I] 63 µm, as well as the underlying continuum at all five wavelengths, in the $z=5.7$ lensed DSFG SPT 0346–52. We reconstruct the lensed continuum and line data using the pixelated, interferometric lens modeling code RIPPLES in order to study the source-plane structure of SPT 0346–52. We analyze both the galaxy-integrated properties and the spatially resolved properties of SPT 0346–52.

We use the photoionization code CLOUDY to model the physical conditions of the ISM in SPT 0346–52. It has lower ionization parameter ($\log U \sim -2.75$) and hydrogen density ($\log n_H / \text{cm}^{-3} \sim 2$) than other high-z DSFGs. Based on CLOUDY modeling, we find supersolar metallicity ($\log Z / Z_{\odot} = 0.75$), similar to what would be expected from the gas-to-dust ratio in SPT 0346–52.

We calculate the dust temperatures throughout SPT 0346–52 and compare the global dust temperature and the wavelength where the SED peaks to other models. We look at line deficits and find deficits in all five lines and spatially resolved deficits in all three detected lines, [N II] 205 µm, [C II] 158 µm, and [O I] 145 µm. We use the limit on the [N II] 122 µm/[N II] 205 µm ratio to find $n_e < 32 \text{ cm}^{-3}$ in SPT 0346–52, which is
lower than what is observed in ULIRGs and other DSFGs. Using [C II] 158 µm/[N II] 205 µm, we determine that ~84% of the [C II] 158 µm emission originates from neutral gas, comparable to other high-z sources and ULIRGs. Using the [C II] 158 µm/[O I] 146 µm ratio, we see that SPT 0346-52 has similar dense gas in PDRs to local galaxies.

Finally, we calculate ionized, neutral, and molecular gas masses using a variety of methods. The molecular gas mass is ~100× the ionized gas mass and ~4× the neutral atomic gas mass. The molecular ISM dominates the mass budget of SPT 0346-52, fueling the intense star formation in this system.

We thank the anonymous referee for the insightful and thorough comments. We also thank Chris Marslender for his computational support. The SPT is supported by the US National Science Foundation (NSF) through grant OPP-1852617. K.C.L., D.P.M., J.D.V., K.P., and S.J. acknowledge support from the US NSF under grants AST-1715213 and AST-1716127. K.C.L. and S.J. acknowledge support from the US NSF NRAO under grants SOSPA4-007 and SOSPA5-001, respectively. J.D.V. acknowledges support from an A. P. Sloan Foundation Fellowship. M.A. acknowledges partial support from FONDECYT grant 1211951, ANID+PCI+INSTITUTO MAX PLANCK DE ASTRONOMIA MPG 190030, ANID +PCI+REDES 190194, and ANID BASAL project FB210003. N.S. is a member of the International Max Planck Research School (IMPRS) for Astronomy and Astrophysics at the Universities of Bonn and Cologne. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This material has made use of the Ocelote high-performance computer, which is part of the

high-performance computing (HPC) resources supported by the University of Arizona TRIF, UITS, and RDI and maintained by the UA Research Technologies department. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2013.1.01231.S, #2015.1.01580.S, and #2016.1.01565.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. This research has made use of NASA's Astrophysics Data System. This research also uses the Cosmology Calculator by Wright (2006). Cloudy has been supported by NSF (1816537), NASA (ATP 17-ATP17-0141), and STScI (HST-AR-15018).

*Software:* CASA (v5.4.0; McMullin et al. 2007; Petry & CASA Development Team 2012), ripples (Hezaveh et al. 2016), visilens (https://github.com/jspilker/visilens; Hezaveh et al. 2013; Spilker et al. 2016), Cloudy (v17.01; Ferland et al. 2017), FSPS (Conroy et al. 2009; Conroy & Gunn 2010), MIST (Choi et al. 2016; Dotter 2016).

## Appendix A
### Lens Modeling Details

Figure 14 shows the observed data and modeled image- and source-plane data, along with the residual image-plane emission and the uncertainty in the source plane. As described in Section 3.1, the uncertainty maps are obtained by creating 500 random sets of visibilities and taking the standard deviation in each pixel of the 500 noise reconstructions.
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Image- and source-plane maps. From left to right: dirty image of data visibilities, dirty image of model visibilities, residual map of dirty images, high-resolution (non-visibility-sampled) image-plane model, source-plane model, and source-plane uncertainty. From top to bottom: 205, 158, 145, 122, 63, [N II] 205, [C II] 158, and [O I] 145 µm. Contours indicate the observed continuum emission at each wavelength. The dirty observed image, dirty model image, and residual image are on the same color scale for each set of models.</figcaption></figure>
## Appendix B
### Spatially Resolved Best-fit CLOUDY Ratios

Table 6 lists the data values compared to the CLOUDY output to determine the best-fit models. For global values, the log of the ratios and the associated errors are listed. For pixelated values, the mean value is listed. The recorded error is the mean

of the uncertainties on the log of the ratios and represents a typical error in a pixel. In parentheses, the minimum and maximum pixel values are listed. Figure 15 shows the observed continuum and line ratios, as well as the best-fit model ratios from CLOUDY. The fitting procedures and best-fit parameters are described in Section 4.

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Maps of observed and modeled ratios from CLOUDY. From left to right: log 63 µm/122 µm, log 63 µm/146 µm, log 63 µm/158 µm, log 63 µm/205 µm, log [O I] 146 µm/146 µm, log [C II] 158 µm/[N II] 205 µm, and log [C II] 158 µm/[O I] 146 µm. Top: observed. Bottom: modeled. Each column has the same color scale.</figcaption></figure>

<table><caption>**Table 6**
Data Used to Fit to CLOUDY Models</caption><thead><tr><td>log Ratio</td><td>Integrated Value<sup>d</sup></td><td>Pixel-by-pixel Values<sup>e</sup></td></tr></thead><tbody><tr><td>63 µm/122 µm<sup>a</sup></td><td>0.05 ± 0.06</td><td>0.10 ± 0.07 (−0.64–0.88)</td></tr><tr><td>63 µm/146 µm<sup>a</sup></td><td>0.33 ± 0.06</td><td>0.36 ± 0.07 (−0.14–0.89)</td></tr><tr><td>63 µm/158 µm<sup>a</sup></td><td>0.37 ± 0.06</td><td>0.55 ± 0.07 (−0.32–1.35)</td></tr><tr><td>63 µm/205 µm<sup>a</sup></td><td>0.62 ± 0.06</td><td>0.77 ± 0.07 (−0.09–1.62)</td></tr><tr><td>[O I] 146 µm/146 µm<sup>b</sup></td><td>0.38 ± 0.09</td><td>0.42 ± 0.44 (0.04–0.99)</td></tr><tr><td>[C II] 158 µm/[N II] 205 µm<sup>c</sup></td><td>1.45 ± 0.08</td><td>1.46 ± 0.62 (1.08–1.69)</td></tr><tr><td>[C II] 158 µm/[O I] 146 µm<sup>c</sup></td><td>1.03 ± 0.09</td><td>0.71 ± 0.62 (0.48–1.11)</td></tr></tbody></table>

**Notes.**

a Continuum flux density ratio.

b Line/continuum ratio in $10^7 L_\odot$ mJy$^{-1}$ (see Section 4.1).

c Line luminosity ratio.

d Values and uncertainties used to compare to CLOUDY models for galaxy-integrated fits.

e For each set of pixelated values, we list $\bar{R} \pm \bar{\sigma}$ ($R_{min} - R_{max}$), where $\bar{R}$ is the mean value in pixels for the log of the ratio, $\bar{\sigma}$ is the mean of the uncertainties in the pixels for the log of the ratio, and $R_{min}$ and $R_{max}$ are the minimum and maximum values of the log of the ratio.
ORCID iDs

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

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

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

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

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

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

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

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

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

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

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

Nikolaus Sulzenauer https://orcid.org/0000-0002-3187-1648

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

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

References

Abdullah, A., Brandl, B. R., Groves, B., et al. 2017, ApJ, 842, 4

Abel, N. P., Dudley, C., Fischer, J., Satyapal, S., & van Hoof, P. A. M. 2009, ApJ, 701, 1147

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

Aravena, M., Spilker, J. S., Béthermin, T., et al. 2016, MNRAS, 457, 4406

Bañados, E., Novak, M., Neeleman, M., et al. 2019, ApJL, 881, L23

Bertelli, G., Bressan, A., Chiosi, C., Fagotto, F., & Nasi, E. 1994, A&AS, 106, 275

Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2

Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44

Carlstrom, J. E., Ade, P. A. R., Aird, K. A., et al. 2011, PASP, 123, 568

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

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833

Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486

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

Coppin, K. E. K., Danielson, A. L. R., Geach, J. E., et al. 2012, MNRAS, 427, 520

Cormier, D., Abel, N. P., Hony, S., et al. 2019, A&A, 626, A23

Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2015, A&A, 578, A53

Cowie, L. L., & Songaila, A. 1986, ARA&A, 24, 499

Croxall, K. V., Smith, J. D., Pellegrini, E., et al. 2017, ApJ, 845, 96

Cunningham, D. J. M., Chapman, S. C., Aravena, M., et al. 2020, MNRAS, 494, 4090

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

Decarli, R., Dotti, M., Bañados, E., et al. 2019, ApJ, 880, 157

Decarli, R., Walter, F., Carilli, C., et al. 2014, ApJL, 782, L17

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2017, ApJ, 846, 32

Dotter, A. 2016, ApJS, 222, 8

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

Farrah, D., Lebouteiller, V., Spoon, H. W. W., et al. 2013, ApJ, 776, 38

Ferkinhoff, C., Brisbin, D., Nikola, T., et al. 2011, ApJL, 740, L29

Ferkinhoff, C., Hailey-Dunsheath, S., Nikola, T., et al. 2010, ApJL, 714, L147

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMxAA, 53, 385

Ferrière, K. M. 2001, RvMP, 73, 1031

Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371

Goldsmith, P. F., Yıldız, U. A., Langer, W. D., & Pineda, J. L. 2015, ApJ, 814, 133

Graciá-Carpio, J., Sturm, E., Hailey-Dunsheath, S., et al. 2011, ApJL, 728, L7

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

Gullberg, B., Swinbank, A. M., Smail, I., et al. 2018, ApJ, 859, 12

Hailey-Dunsheath, S., Nikola, T., Stacey, G. J., et al. 2010, ApJL, 714, L162

Hayward, C. C., Jonsson, P., Kereš, D., et al. 2012, MNRAS, 424, 951

Helou, G., Khan, I. R., Malek, L., & Boehmer, L. 1988, ApJS, 68, 151

Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2018, ApJ, 861, 94

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

Hezaveh, Y. D., Marrone, D. P., Fassnacht, C. D., et al. 2013, ApJ, 767, 132

Hollenbach, D. J., Takahashi, T., & Tielens, A. G. G. M. 1991, ApJ, 377, 192

Hollenbach, D. J., & Tielens, A. G. G. M. 1999, RvMP, 71, 173

Hopkins, P. F., Cox, T. J., Younger, J. D., & Hernquist, L. 2009, ApJ, 691, 1168

Iono, D., Tamura, Y., Nakanishi, K., et al. 2006, PASJ, 58, 957

Ivison, R. J., Smail, I., Papadopoulos, P. P., et al. 2010, MNRAS, 404, 198

Jameson, K. E., Bolatto, A. D., Wolfire, M., et al. 2018, ApJ, 853, 111

Jarugula, S., Vieira, J. D., Weiss, A., et al. 2021, ApJ, 921, 97
Jones, G. C., Maiolino, R., Caselli, P., & Carniani, S. 2019, A&A, 632, L7
Jones, G. C., Maiolino, R., Caselli, P., & Carniani, S. 2020, MNRAS, 498, 4109

Katz, H., Galligan, T. P., Kimm, T., et al. 2019, MNRAS, 487, 5902

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

Kramer, C., Nikola, T., Anderl, S., et al. 2020, A&A, 639, A61

Lamarche, C., Verma, A., Vishwas, A., et al. 2018, ApJ, 867, 140

Launay, J. M., & Roueff, E. 1977, A&A, 56, 289

Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1

Lee, M. M., Nagao, T., De Breuck, C., et al. 2019, ApJL, 883, L29

Lee, M. M., Nagao, T., De Breuck, C., et al. 2021, ApJ, 913, 41

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

Levesque, E. M., Kewley, L. J., & Larson, K. L. 2010, AJ, 139, 712

Li, J., Wang, R., Cox, P., et al. 2020, ApJ, 900, 131

Liseau, R., Justtanont, K., & Tielens, A. G. G. M. 2006, A&A, 446, 561

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

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

Luhman, M. L., Satyapal, S., Fischer, J., et al. 1998, ApJL, 504, L11

Lutz, D., Berta, S., Contursi, A., et al. 2016, A&A, 591, A136

Ma, J., Gonzalez, A. H., Spilker, J. S., et al. 2015, ApJ, 812, 88

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

Madden, S. C., Rémy-Ruyer, A., Galametz, M., et al. 2013, PASP, 125, 600

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6

Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 466, 561

Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883

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

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Aolap, K. 2007, in ASP Conf. Ser. 376, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

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

Meynet, G., & Maeder, A. 2000, A&A, 361, 101

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

Oberst, T. E., Parshley, S. C., Stacey, G. J., et al. 2006, ApJL, 652, L125

Oteo, I., Ivison, R. J., Dunne, L., et al. 2016, ApJ, 827, 34

Parkin, T. J., Wilson, C. D., Schirm, M. R. P., et al. 2013, ApJ, 776, 65

Pavesi, R., Riechers, D. A., Capak, P. L., et al. 2016, ApJ, 832, 151

Pavesi, R., Riechers, D. A., Sharon, C. E., et al. 2018, ApJ, 861, 43

Pensabene, A., Decarli, R., Bañados, E., et al. 2021, A&A, 652, A66

Pereira-Santaella, D., Rigopoulou, D., Farrah, D., Lebouteiller, V., & Li, J. 2017, MNRAS, 470, 1218

Petry, D. & CASA Development Team 2012, in ASP Conf. Ser. 461, Astronomical Data Analysis Software and Systems XXI, ed. P. Ballester, D. Egret, & N. P. F. Lorente (San Francisco, CA: ASP), 849

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

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

Riechers, D. A. 2013, ApJL, 765, L31

Rybak, M., Zavala, J. A., Hodge, J. A., Casey, C. M., & Werf, P. V. d. 2020, ApJL, 889, L11

Rybak, M., Calistro Rivera, G., Hodge, J. A., et al. 2019, ApJ, 876, 112

Sargsyan, L., Lebouteiller, V., Weedman, D., et al. 2012, ApJ, 755, 171

Savage, B. D., & Sembach, K. R. 1996, ARA&A, 34, 279

Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269

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

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

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

Spilker, J. S., Phadke, K. A., Aravena, M., et al. 2020, ApJ, 905, 85

Spinoglio, L., & Malkan, M. A. 1992, ApJ, 399, 504

Stacey, G. J. 2011, ITTST, 1, 241

Sturm, E., Verma, A., Graciá-Carpio, J., et al. 2019, A&A, 518, L36

Sutter, J., Dale, D. A., Croxall, K. V., et al. 2010, ApJ, 886, 60

Suyu, S. H., Marshall, P. J., Hobson, M. P., & Blandford, R. D. 2006, MNRAS, 371, 983

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

Tayal, S. S. 2008, A&A, 486, 629

Tayal, S. S. 2011, ApJS, 195, 12

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

Toft, S., Smolčić, V., Magnelli, B., et al. 2014, ApJ, 782, 68

Venemans, B. P., Walter, F., Decarli, R., et al. 2017, ApJL, 851, L8
Vieira, J. D., Crawford, T. M., Switzer, E. R., et al. 2010, ApJ, 719, 763

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

Wagg, J., Aravena, M., Brisbin, D., et al. 2020, MNRAS, 499, 1788

Wang, S. X., Brandt, W. N., Luo, B., et al. 2013, ApJ, 778, 179

Warren, S. J., & Dye, S. 2003, ApJ, 590, 673

Weiß, A., De Breuck, C., Marrone, D. P., et al. 2013, ApJ, 767, 88

Weiß, A., Downes, D., Henkel, C., & Walter, F. 2005, 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

Wright, E. L. 2006, PASP, 118, 1711

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976

Zavala, J. A., Montaña, A., Hughes, D. H., et al. 2018, NatAs, 2, 56

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

Zhao, Y., Lu, N., Xu, C. K., et al. 2016, ApJ, 819, 69