# Spatially Resolved Water Emission from Gravitationally Lensed Dusty Star-forming Galaxies at $z \sim 3$

Sreevani Jarugula, Joaquin D. Vieira, Justin S. Spilker, Yordanka Apostolovski, Manuel Aravena, Matthieu Béthermin, Carlos de Breuck, Chian-Chou Chen, Daniel J. M. Cunningham, Chenxing Dong, et al.

► To cite this version:

Sreevani Jarugula, Joaquin D. Vieira, Justin S. Spilker, Yordanka Apostolovski, Manuel Aravena, et al.. Spatially Resolved Water Emission from Gravitationally Lensed Dusty Star-forming Galaxies at $z \sim 3$. The Astrophysical Journal, 2019, 880, 10.3847/1538-4357/ab290d. insu-03666737

HAL Id: insu-03666737

https://insu.hal.science/insu-03666737v1

Submitted on 12 Aug 2025

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

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

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

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

# Spatially Resolved Water Emission from Gravitationally Lensed Dusty Star-forming Galaxies at $z \sim 3$

Sreevani Jarugula¹, Joaquin D. Vieira¹,²,³, Justin S. Spilker⁴, Yordanka Apostolovski⁵, Manuel Aravena⁵, Matthieu Béthermin⁶, Carlos de Breuck⁷, Chian-Chou Chen⁷, Daniel J. M. Cunningham⁸,⁹, Chenxing Dong¹⁰, Thomas Greve¹¹,¹², Christopher C. Hayward¹³, Yashar Hezaveh¹⁴, Katrina C. Litke¹⁵, Amelia C Mangian¹, Desika Narayanan¹⁰,¹²,¹⁶, Kedar Phadke¹, Cassie A. Reuter¹, Paul Van der Werf¹⁷, and Axel Weiss¹⁸

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

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

³National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, 1205 W. Clark St., Urbana, IL 61801, USA

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

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

⁶Aix Marseille Univ., Centre National de la Recherche Scientifique, Laboratoire d'Astrophysique de Marseille, Marseille, France

⁷European Southern Observatory, Karl Schwarzschild Straße 2, D-85748 Garching, Germany

⁸Department of Astronomy and Physics, Saint Mary's University, Halifax, NS, B3H 3C3, Canada

⁹Department of Physics and Atmospheric Science, Dalhousie University, Halifax, NS, B3H 4R2, Canada

¹⁰Department of Astronomy, University of Florida, Gainesville, FL 32611, USA

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

¹²Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen, Juliane Maries vej 30, DK-2100 Copenhagen, Denmark; DTU-Space, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark

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

¹⁴Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94305, USA

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

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

¹⁷Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands

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

Received 2019 March 14; revised 2019 May 31; accepted 2019 June 10; published 2019 July 30

## Abstract

Water (H₂O), one of the most ubiquitous molecules in the universe, has bright millimeter-wave emission lines that are easily observed at high redshift with the current generation of instruments. The low-excitation transition of H₂O, p-H₂O(2₀,₂ − l₁,₁) (ν<sub>rest</sub> = 987.927 GHz), is known to trace the far-infrared (FIR) radiation field independent of the presence of active galactic nuclei (AGNs) over many orders of magnitude in FIR luminosity (L<sub>FIR</sub>). This indicates that this transition arises mainly due to star formation. In this paper, we present spatially (~0.′′5 corresponding to ~1 kiloparsec) and spectrally resolved (~100 kms<sup>-1</sup>) observations of p-H₂O(2₀,₂ − l₁,₁) in a sample of four strong gravitationally lensed high-redshift galaxies with the Atacama Large Millimeter/ submillimeter Array. In addition to increasing the sample of luminous (>10<sup>12</sup> L<sub>☉</sub>) galaxies observed with H₂O, this paper examines the L<sub>H₂O</sub>/L<sub>FIR</sub> relation on resolved scales for the first time at high redshift. We find that L<sub>H₂O</sub> is correlated with L<sub>FIR</sub> on both global and resolved kiloparsec scales within the galaxy in starbursts and AGN with average L<sub>H₂O</sub>/L<sub>FIR</sub> = 2.76<sup>+2.15</sup><sub>-1.21</sub> × 10<sup>-5</sup>. We find that the scatter in the observed L<sub>H₂O</sub>/L<sub>FIR</sub> relation does not obviously correlate with the effective temperature of the dust spectral energy distribution or the molecular gas surface density. This is a first step in developing p-H₂O(2₀,₂ − l₁,₁) as a resolved star formation rate calibrator.

**Key words:** galaxies: high-redshift – galaxies: ISM – galaxies: starburst – ISM: molecules

## 1. Introduction

Studies of molecules play a prominent role in explaining the physical, chemical, and kinematic properties of the interstellar medium (ISM) in galaxies (Omont 2007; Tielens 2013). One such molecule is H₂O, the third most abundant molecule in the warm dense ISM after H₂ and CO (Neufeld et al. 1995). As an asymmetric rotor with a large electric dipole moment, H₂O has a rich and complex spectrum giving rise to emission and absorption lines mainly in the submillimeter and far-infrared (FIR) regimes of the electromagnetic spectrum. Observations from local galaxies (van der Werf et al. 2010; Weiß et al. 2010; Rangwala et al. 2011; Yang et al. 2013), high-redshift ultra-luminous infrared galaxies (ULIRGs; Omont et al. 2013; Yang et al. 2016), and active galactic nuclei (AGNs; van der Werf et al. 2011) have shown H₂O emission to be ubiquitous with intensities as bright as CO lines. Modeling has shown that, in addition to infrared pumping where H₂O is excited by FIR

photons, collisions also contribute to the intensities of low-excitation transitions (e.g., González-Alfonso et al. 2010, 2012). This is best represented in Figure 3 from Liu et al. (2017), which shows prominent H₂O lines in different ISM components. The low-excitation lines become weaker or completely disappear in the warm and hot regions (>40 K) where infrared pumping dominates over collisions. The higher excitation transitions that require strong far-infrared radiation density are mainly found in the hotter regions (100–200 K) of the galaxy. The cascading emission lines, p-H₂O(2₀,₂ − l₁,₁) (E<sub>up</sub> = 100.8 K, ν<sub>rest</sub> = 987.927 GHz), p-H₂O(2₁,₁ − 2₀,₂) (E<sub>up</sub> = 137 K, ν<sub>rest</sub> = 752.033 GHz), and p-H₂O (2₂,₀ − 2₁,₁) (E<sub>up</sub> = 196 K, ν<sub>rest</sub> = 1228.789 GHz) are pumped by 101 µm photons from the base l₁,₁ level and are primarily excited in the warm regions of the galaxy. The collisional excitation of the low-lying levels (l₁,₁ and 2₀,₂) in optically thin or high-density hot regions might also contribute to the emission of the p-H₂O(2₀,₂ − l₁,₁)
line. Hence, H₂O transitions probe the infrared radiation field density and physical properties of the ISM such as gas density and kinetic temperature (e.g., Weiß et al. 2010; González-Alfonso et al. 2014; Liu et al. 2017).

Because of water vapor in the Earth's atmosphere, ground-based observations of H₂O in the local universe are nearly always impossible. The Herschel Space Observatory opened the window to multiple H₂O transitions in the local universe (e.g., Weiß et al. 2010) and Yang et al. (2013) demonstrated that the luminosity of submillimeter H₂O lines ($L_{H_2O}$) is linearly correlated with the total infrared luminosity ($L_{IR}$, integrated over 8–1000 µm) over three orders of magnitude in multiple transitions. This suggests that the H₂O transitions, especially p-H₂O(2₀,₂ − 1₁,₁), which is not affected by the presence of AGNs (Yang et al. 2013), trace the far-infrared field in star-forming regions. At high redshift, H₂O has been detected using the current generation of ground-based telescopes such as the CSO, PdBI, and ALMA, as the transitions are redshifted into the transparent millimeter atmospheric windows. Strong gravitational lensing, which acts as a cosmic microscope, further boosts the flux from high-redshift sources, making their detections possible. Several detections of H₂O have been reported in the literature from such lensed galaxies (e.g., Bradford et al. 2009; Omont et al. 2011; van der Werf et al. 2011; Combes et al. 2012; Bothwell et al. 2013a; Omont et al. 2013; Weiß et al. 2013; Spilker et al. 2014; Yang et al. 2016).

Multi-wavelength observations ranging from the UV to radio have improved our understanding of interstellar physics and the star formation rate (SFR) calibration. Average scaling relations from single observables are often used to estimate global SFR. Obtaining resolved SFR maps is challenging due to the difficulty of observing individual star-forming regions over multiple wavelengths. The far-infrared luminosity of galaxies ($L_{FIR}$, integrated over 42.5–122.5 µm) is often used to infer SFR, as it has some advantages over other indicators such as UV luminosity and recombination lines, which are widely discussed in Kennicutt (1998) and Kennicutt & Evans (2012). The UV emission from young stars is a direct tracer of star formation but is highly sensitive to interstellar dust attenuation. Recombination lines such as Hα and FIR cooling lines (e.g., [CII] 158 µm) originate in the ionized regions surrounding stars and are good tracers of star formation. However, these lines are affected either by dust attenuation (e.g., Casey et al. 2017) or the scatter in the estimated SFR is large (e.g., Narayanan & Krumholz 2017; Lagache et al. 2018). In contrast, $L_{FIR}$ is a good tracer of SFR at high optical depth, such as starburst galaxies where most of the UV light is re-emitted as infrared radiation. Although it is widely used as an SFR calibrator in high-redshift starburst galaxies (see review by Casey et al. (2014)), the spectral energy distribution (SED) has to be fully sampled over the SED peak at $\lambda_{rest} \sim 100$ µm to estimate $L_{FIR}$, which is observationally expensive. However, one further caveat is that infrared emission does not necessarily trace only the unobscured star formation. For instance, $L_{FIR}$ may overestimate the SFR in regions where there are other sources of dust heating such as evolved older stars or an obscured AGN (e.g., Kennicutt et al. 2009; Murphy et al. 2011; Hayward et al. 2014). Longer-wavelength spectral features such as p-H₂O(2₀,₂ − 1₁,₁), which are very well correlated with $L_{FIR}$ and observable with current generation telescopes, can be used instead of $L_{FIR}$ to estimate SFR (in environments where

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Ratio of line to far-infrared luminosity of p-H₂O(2₀,₂ − 1₁,₁) (988 GHz), CO(1–0) (115 GHz), CO(6–5) (691 GHz), and [CII] (1900 GHz) as a function of $L_{FIR}$. The SPT sources are shown as diamonds. The p-H₂O(2₀,₂ − 1₁,₁) emission is described in detail in Figure 3. The H₂O emission in the local galaxies is described in Yang et al. (2013) and the emission from high-redshift galaxies is taken from van der Werf et al. (2011), Omont et al. (2013), Yang et al. (2016), Apostolovski et al. (2019), and this paper. CO(1–0) emission from local ULIRGs is given in Solomon et al. (1997) and the ATCA observations of CO(1–0) in SPT sources (green diamonds) are described in detail in Aravena et al. (2016). CO(6–5) line emission from local luminous infrared galaxies (LIRGs) and the SPT sources (golden-yellow diamonds) are from Lu et al. (2017) and Dong et al. (2019), respectively. The [CII] sample of LIRGs is taken from Díaz-Santos et al. (2014) and the [CII] SPT sample represented by purple diamonds is from Gullberg et al. (2015). As seen in the plot, the luminosities of CO(1–0), CO(6–5) and [CII] are sublinearly correlated with $L_{FIR}$, while p-H₂O(2₀,₂ − 1₁,₁) is almost linearly correlated with $L_{FIR}$ especially for $L_{FIR} > 10^{11.5} L_{⊙}$.</figcaption></figure>

the $L_{FIR}$ based calibration holds true). Figure 1 summarizes some of the measurements from the literature and shows that p-H₂O(2₀,₂ − 1₁,₁) is almost linearly correlated with $L_{FIR}$ with Pearson's correlation coefficient of ~0.96. Among the CO transitions, it has been observed that mid to high-J CO transitions (e.g., CO(6–5) and CO(7–6)) are a good tracer of $L_{IR}$ both in local and high redshift (U)LIRGs (e.g., Lu et al. 2015; Yang et al. 2017). However, sublinear slopes in the $L_{FIR}-L_{CO}$ correlation arising possibly from shocks/turbulence and detached from star formation have also been discussed in high-J CO lines (e.g., Greve et al. 2014, see Section 4.4). While CO is collisionally excited by H₂ molecules, p-H₂O(2₀,₂ − 1₁,₁) is excited by FIR photons, which makes H₂O a more direct tracer of star formation. In nearby luminous galaxies, dense gas tracers such as HCN and CS are shown to be tightly correlated with $L_{IR}$, while HCO+ has a slightly super-linear correlation (e.g., Gao & Solomon 2004; Zhang et al. 2014). p-H₂O(2₀,₂ − 1₁,₁) is a bright emission line (compared to HCN/HCO+) that is easily observable both in local and high-redshift galaxies. The linear correlation between $L_{H_2O}$ and $L_{FIR}$ from Figure 1 suggests that it is a better tracer of $L_{FIR}$ compared to other commonly observed lines such as CO(1–0), CO(6–5), and [CII]. While the correlation is tight on global integrated scales, it is unclear if this correlation breaks down on resolved scales.

In this work, we show that $L_{H_2O}$ traces far-infrared radiation not just at the integrated global scale (Omont et al. 2013; Yang et al. 2013, 2016) but also at resolved scales within galaxies at high redshift. The resolution of the observations is ~0.′′5, which corresponds to ~1 kiloparsec given the magnification
<table><caption>**Table 1**
Source Properties and Observations</caption><thead><tr><th>Source name</th><th>IAU name</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>z</th><th>ν<sup>line</sup><sub>obs</sub><br>[GHz]</th><th>ν<sup>cont</sup><sub>obs</sub><br>[GHz]</th><th>Beam<br>["]</th><th>σ<sub>cont</sub><br>[mJy/beam]</th><th>SNR<sub>line</sub></th></tr></thead><tbody><tr><td>SPT0529-54</td><td>SPT-S J052903-5436.6</td><td>05:29:03.37</td><td>-54:36:40.30</td><td>3.3689</td><td>226.13</td><td>224.87</td><td>1.10 × 1.06</td><td>0.17</td><td>7.5</td></tr><tr><td>SPT0532-50</td><td>SPT-S J053250-5047.1</td><td>05:32:51.27</td><td>-50:47:09.50</td><td>3.3988</td><td>224.59</td><td>223.80</td><td>0.73 × 0.56</td><td>0.18</td><td>24.4</td></tr><tr><td>SPT0538-50</td><td>SPT-S J053816-5030.8</td><td>05:38:16.83</td><td>-50:30:52.00</td><td>2.7817</td><td>260.98</td><td>259.32</td><td>0.63 × 0.51</td><td>0.22</td><td>18.6</td></tr><tr><td>Cloverleaf</td><td>H 1413+117</td><td>14:15:46.24</td><td>11:29:43.68</td><td>2.5579</td><td>277.67</td><td>278.83</td><td>0.54 × 0.50</td><td>1.07</td><td>11.8</td></tr></tbody></table>

Note. The position (R.A., decl.) and redshift (z) of the SPT sources are taken from ALMA 870 µm imaging in Spilker et al. (2016) and Weiß et al. (2013), respectively. The redshift of the Cloverleaf quasar is found in Solomon et al. (2003). $\nu_{\text{obs}}^{\text{line}}$ is the observed frequency of the $p$-H$_2$O($2_{0,2} - l_{1,1}$) transition at a 987.927 GHz rest frequency. $\nu_{\text{obs}}^{\text{cont}}$ is the frequency of the continuum. The beam size and the sensitivity per beam in the continuum map ($\sigma_{\text{cont}}$) at $\nu_{\text{obs}}$ are shown in columns 8 and 9. The peak signal-to-noise of the H$_2$O line (SNR$_{\text{line}}$), with a 50 kms$^{-1}$ channel resolution in all sources except SPT0529-54, with 100 kms$^{-1}$ channel width, is shown in the last column.

and redshift of the sources from Spilker et al. (2014) (the beam resolution and the physical scale for each source are given in Tables 1 and 3, respectively). This physical scale is only an approximation, as we do not perform lens modeling in this analysis and adopt magnification values obtained from 870 µm imaging. We have selected a sample of strong gravitationally lensed dusty star-forming galaxies (DSFGs) discovered in the South Pole Telescope (SPT) survey (Vieira et al. 2010; Carlstrom et al. 2011; Mocanu et al. 2013). DSFGs host intense star formation with SFR > 10$–1000 M_{\odot} \text{yr}^{-1}$ (e.g., Casey et al. 2014; Narayanan et al. 2015). These galaxies are bright in submillimeter wavelengths, as the ultraviolet (UV) radiation from young stars is absorbed and re-radiated by the dust in FIR. Long-wavelength dust continuum observations of such galaxies have the advantage of “negative K correction” (Blain & Longair 1993), where the decrease in flux due to increase in cosmological distance is compensated by the rising flux on the Rayleigh-Jeans side of the SED. Thus, sources of a given luminosity can be detected largely independent of redshift. This, in addition to gravitational lensing and the power of ALMA, provides enough sensitivity and resolution to investigate the correlation between $L_{\text{H}_2\text{O}}$ and $L_{\text{FIR}}$ at resolved scales in star-forming galaxies, which we present in this paper.

In Section 2, we summarize the ALMA observations taken over two cycles and the data reduction procedure. In Section 3, we present the results estimating the infrared luminosity and line properties. In Section 4, we analyze the results for the $L_{\text{H}_2\text{O}}-L_{\text{FIR}}$ correlation and the effect of physical properties on this correlation. We conclude with a summary in Section 5. Hereafter, H$_2$O refers to $p$-H$_2$O($2_{0,2} - l_{1,1}$) at 987.927 GHz and $L_{\text{H}_2\text{O}}/L_{\text{FIR}}$ refers to $L_{\text{H}_2\text{O}(2_{0,2}-l_{1,1})}/L_{\text{FIR}}$. We use a Planck 2015 flat $\Lambda$CDM cosmology where $h = 0.677$, $\Omega_m = 0.307$, and $\Omega_{\Lambda} = 0.693$ (Planck Collaboration et al. 2016). We estimate the total infrared luminosity ($L_{\text{IR}}$) as flux-integrated from 8 to 1000 µm and total far-infrared luminosity ($L_{\text{FIR}}$) from 42.5 to 122.5 µm in the rest frame (Helou et al. 1985).

## 2. Observations and Data Analysis

We choose $p$-H$_2$O($2_{0,2} - l_{1,1}$) as it is one of the brightest H$_2$O transitions and has been observed to be well correlated with $L_{\text{IR}}$ (e.g., Yang et al. 2013; Liu et al. 2017). This line also falls in the transparent ALMA Band 6 for the given redshift range of the sources ($z \sim 2.78–3.37$). We observed the $p$-H$_2$O($2_{0,2} - l_{1,1}$) 987.927 GHz line in SPT0529-54 ($z = 3.369$), SPT0532-50 ($z = 3.399$), and SPT0538-50 ($z = 2.782$) with ALMA. We also

include archival data on the Cloverleaf quasar, a strongly lensed AGN at $z = 2.558$ in this analysis. The source properties are listed in Table 1.

### 2.1. Sample Selection

The three SPT targets were selected such that they are at a similar redshift and within 10° of each other on the sky. This selection was chosen to observe the same line transition in the three galaxies and to make observations efficient for resolved ALMA Band 9 continuum observations, which were A-rated in Cycle 5, but not yet observed. All three sources have ALMA 870 µm imaging and lens models (Spilker et al. 2016). SPT0538-50 is a possible ongoing major merger as seen from dust continuum models (Bothwell et al. 2013b) and has resolved CO(1–0) and CO(3–2) ATCA observations (Aravena et al. 2013; Spilker et al. 2015). SPT0529-54 and SPT0532-50 have resolved CO(6–5) observations from ALMA (Dong et al. 2019) which we make use of in this work. The Cloverleaf quasar (also known as H1413+117 or QSO J1415+1129) is an extensively studied strongly lensed AGN at high redshift (e.g., Solomon et al. 2003; Weiß et al. 2003) that can be compared against the star-forming galaxies in this sample.

### 2.2. Data Reduction and Imaging

For SPT0532-50 and SPT0538-50, the observations are taken with ALMA in Cycle 3 (2015.1.01578.S) and Cycle 4 (2016.1.01554.S). SPT0529-54 was observed only in Cycle 4. The Cloverleaf quasar archival data are from ALMA 2012.1.00175.S (PI: Van Der Werf). In both Cycle 3 and Cycle 4, the SPT sources are observed in two basebands each with 2 GHz bandwidth and 240 channels with a channel resolution of 1.875 MHz. The $p$-H$_2$O($2_{0,2} - l_{1,1}$) line and continuum in the Cloverleaf were taken over a continuous bandwidth with two basebands each 2 GHz wide, 128 channels, and a channel resolution of 15.625 MHz. The observations are summarized in Table 2.

All the data are calibrated using the ALMA pipeline for the respective cycles. We inspected the quality of the reduction manually and found no major problems. The data were reduced and imaged using the Common Astronomy Software Application package CASA (McMullin et al. 2007). We imaged the continuum by combining data from all the spectral windows and by excluding the line emission using the task CLEAN. The frequency of the continuum image is given in Table 1. An outer taper of 1"0 and 0"5 is applied to SPT0529-54 and the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top: the continuum is shown as a background in log scale with the minimum pixel value as $3\sigma_{cont}$ of the continuum map and the moment 0 contours of the H$_2$O emission are overlaid in white. The contours are at [3,5,10,15 ...] × σ, where σ is the rms noise in the H$_2$O moment 0 map. The synthesized beam of the continuum image is shown in the lower left corner and the spatial scale bar of 1.0" is shown in the lower right. Bottom: the spatially integrated spectrum of the p-H$_2$O(2$_{0,2}$ − 1$_{1,1}$) transition, with a 50 kms$^{-1}$ spectral resolution in all sources except SPT0529-54 with a 100 kms$^{-1}$ resolution. The colored region shows (line center − 3σ$_{ν}$) kms$^{-1}$ to (line center + 3σ$_{ν}$) kms$^{-1}$ (definition given in the text.)</figcaption></figure>

<table><caption>**Table 2**
ALMA Observations</caption><thead><tr><th>Cycle</th><th>Proposal ID</th><th>Source</th><th>Date</th><th>Time on source (hr)</th><th>Antennas</th><th>Baseline (m)</th><th>PWV (mm)</th></tr></thead><tbody><tr><td rowspan="3">4</td><td rowspan="3">2016.1.01554.S</td><td>SPT0529-54</td><td>2016 Dec 5</td><td>1.01</td><td>41</td><td>15.1-650.3</td><td>1.57-1.84</td></tr><tr><td>SPT0532-50</td><td>2016 Dec 3</td><td>0.38</td><td>40</td><td>15.1-704.1</td><td>0.53-0.56</td></tr><tr><td>SPT0538-50</td><td>2016 Dec 5</td><td>0.90</td><td>41</td><td>15.1-650.3</td><td>1.57-1.78</td></tr><tr><td rowspan="2">3</td><td rowspan="2">2015.1.01578.S</td><td>SPT0532-50</td><td>2016 Jun 22</td><td>0.41</td><td>38</td><td>15.1-704.1</td><td>1.27-1.32</td></tr><tr><td>SPT0538-50</td><td>2016 Jun 22</td><td>0.74</td><td>38</td><td>15.1-704.1</td><td>1.23-1.28</td></tr><tr><td>2</td><td>2012.1.00175.S</td><td>Cloverleaf</td><td>30 June, 2015</td><td>0.16</td><td>42</td><td></td><td></td></tr></tbody></table>

**Note.** The baseline given in the table is the minimum and maximum. The last column shows the range of Precipitable Water Vapor (PWV) over the course of the observations.

Cloverleaf, respectively, such that the visibilities at shorter baselines are weighted more. This increases the signal-to-noise at the expense of the resolution of the image. A natural weighting is applied to all the sources. All the pixels in the image plane are correlated in the interferometric data. In order to have a minimum number of correlated pixels, we choose to have a few pixels in a beam (~3 to 5). The continuum images are shown in Figure 2. We create a mask such that only the pixels with signal-to-noise ≥3 in the continuum are selected. We use this mask for all the resolved analysis in this paper.

To get the spectral cubes, we use natural weighting for all the sources and 1"/0 and 0"/5 out tapers to SPT0529-54 and the Cloverleaf, respectively, the same as the continuum map. We use 50 kms$^{-1}$ velocity averaging for SPT0532-50, SPT0538-50 and the Cloverleaf and 100 kms$^{-1}$ in SPT0529-54. For the velocity-integrated intensity map (moment 0), in order to increase the signal-to-noise ratio, we re-image the data to create a single wide channel that contains most of the line flux. The width of this channel is ~2 × FWHM of the line (see Table 3 for FWHM and integrated line flux values). All the moment 0 maps are imaged

similar to the continuum map and the cube. The moment 0 contours are overlaid on the continuum image in Figure 2.

It has been shown by many studies that the CO gas sizes can be larger than those of infrared emission (e.g., Spilker et al. 2015; Tadaki et al. 2017; Calistro Rivera et al. 2018; Dong et al. 2019). As H$_2$O and CO are observed to have similar line profiles (e.g., Omont et al. 2013; Yang et al. 2016, and Appendix A3 (C)), it is likely that H$_2$O and CO are both tracing similar regions. Hence, H$_2$O could also be more extended than dust. In our analysis, the mask selected from the continuum includes 95%–100% of the total H$_2$O flux (depending on the source) and any possible additional extended emission would be small and not affect our results. Moreover, selecting a mask based on continuum is less biased because of the high signal-to-noise in every pixel in the continuum, unlike the moment 0 map.

## 3. Results

### 3.1. Estimating $L_{FIR}$

To fit the SED (and estimate the total $L_{FIR}$) in the SPT sources, we use the unresolved millimeter and submillimeter
<table><caption>**Table 3**
Observed Continuum and Line Properties</caption><thead><tr><th>Source</th><th>μ</th><th>μ<sup>a</sup><br>(kpc)</th><th>S<sub>cont</sub><br>(mJy)</th><th>L<sub>FIR</sub><br>(10<sup>12</sup>L<sub>☉</sub>)</th><th>S<sup>Peak</sup><sub>H<sub>2</sub>O</sub><br>(mJy)</th><th>ΔV<sub>H<sub>2</sub>O</sub><br>(km s<sup>-1</sup>)</th><th>I<sub>H<sub>2</sub>O</sub><br>(Jy km s<sup>-1</sup>)</th><th>L<sub>H<sub>2</sub>O</sub><br>(10<sup>8</sup>L<sub>☉</sub>)</th><th>L<sub>H<sub>2</sub>O</sub>/L<sub>FIR</sub><br>(10<sup>-5</sup>)</th></tr></thead><tbody><tr><td>SPT0529-54</td><td>13.2 ± 0.8</td><td>2.3</td><td>32.9 ± 0.1</td><td>30.84 ± 4.96</td><td>7.7 ± 1.0</td><td>675 ± 235</td><td>2.8 ± 0.6</td><td>5.94 ± 1.18</td><td>1.92 ± 0.49</td></tr><tr><td>SPT0532-50</td><td>10.0 ± 0.6</td><td>1.7</td><td>42.3 ± 0.1</td><td>68.33 ± 9.81</td><td>17.6 ± 0.6</td><td>475 ± 44</td><td>10.4 ± 0.2</td><td>22.25 ± 0.38</td><td>3.26 ± 0.47</td></tr><tr><td>SPT0538-50</td><td>20.1 ± 1.8</td><td>1.1</td><td>58.3 ± 0.1</td><td>63.91 ± 10.55</td><td>27.7 ± 1.5</td><td>463 ± 144</td><td>8.8 ± 1.3</td><td>13.59 ± 2.08</td><td>2.13 ± 0.48</td></tr><tr><td>Cloverleaf</td><td>11</td><td>0.4</td><td>20.1 ± 0.2</td><td>70.91 ± 15.47</td><td>49.9 ± 3.7</td><td>434 ± 156</td><td>17.0 ± 2.4</td><td>22.63 ± 3.17</td><td>3.19 ± 0.83</td></tr></tbody></table>

**Note.** Magnification (μ) of the SPT sources is taken from ALMA 870 μm lens models from Spilker et al. (2016). Magnification of the Cloverleaf quasar is based on Venturini & Solomon (2003). a) *r* is the effective resolution in kiloparsecs achieved per beam (beam size is given in Table 1) in the source plane. S<sub>cont</sub> is the spatially integrated continuum flux (≥3σ<sub>cont</sub>) and L<sub>FIR</sub> is the observed far-infrared luminosity. S<sup>Peak</sup><sub>H<sub>2</sub>O</sub> is the peak line flux and ΔV<sub>H<sub>2</sub>O</sub> is the FWHM of the spatially integrated spectrum, which is derived using the non-parametric method of estimating line properties (Section 3.3). H<sub>2</sub>O is velocity-integrated line flux under the colored region shown in Figure 2. L<sub>H<sub>2</sub>O</sub>/L<sub>FIR</sub> takes into account both the error on luminosities and magnification.

photometry from ALMA (3 mm), SPT (2.0 and 1.4 mm), LABOCA (870 μm), and Herschel (500, 350, 250, 160, and 100 μm) (Vieira et al. 2013; Strandet et al. 2016). For the Cloverleaf, we use the unresolved photometry from Weiß et al. (2003) (references therein). We fit the modified blackbody function given by Equation (1) with the Markov Chain Monte Carlo (MCMC) algorithm using the emcee (Foreman-Mackey et al. 2013) package to sample the posterior probability function:

$$S_{\nu} = \frac{\Omega}{(1 + z)^3} (B_{\nu}(T_d) - B_{\nu}(T_{\text{CMB}}))(1 - e^{-\tau}), \quad (1)$$

where *ν* is the rest frequency, Ω is the source solid angle, B<sub>ν</sub>(T<sub>d</sub>) is the Planck function estimated at dust temperature T<sub>d</sub>, and τ is the optical depth. At long wavelengths, τ is given by τ = (λ/λ<sub>o</sub>)<sup>β</sup> with λ<sub>o</sub> being the wavelength at which the optical depth is unity (e.g., Draine 2006) and β is the spectral index that determines the slope of the Rayleigh-Jeans tail of the blackbody. This is the same method used previously in, e.g., Greve et al. (2012) and Spilker et al. (2016). It should be noted that this simplistic modified blackbody function applies only with the assumption that the entire source has a single T<sub>d</sub> and that the source is uniform (e.g., Hayward et al. 2012). Moreover, β, τ, and temperature distribution are degenerate (Papadopoulos et al. 2010).

In this analysis we consider photometry points ≥50 μm rest frame. The β value converges to ~2 for the SPT sample, which has a well-sampled SED. For the SPT sources and the Cloverleaf, we fix β to 2 and let the amplitude (Ω/(1+z)<sup>3</sup>), T<sub>d</sub>, and λ<sub>o</sub> vary. To investigate L<sub>H<sub>2</sub>O</sub>/L<sub>FIR</sub> at resolved scales, the L<sub>FIR</sub> for each pixel is obtained by scaling the L<sub>FIR</sub> of the entire source using the flux contribution of each pixel to the total continuum flux, i.e.,

$$L_{\text{FIR}}^{i} = \frac{S^{i}}{S} L_{\text{FIR}}, \quad (2)$$

where *i* denotes *i*th pixel and *S* is the continuum flux obtained by combining all the spectral windows from the observations excluding the line. It should be noted that the resolved L<sub>FIR</sub> is estimated by scaling the continuum that is not in the FIR regime, i.e., not around the peak of the SED (the frequencies of the continuum are given in Table 1). Hence, the variations in dust temperature, optical depth, etc., across the source are not taken into account in this analysis, and we assume that the galaxies have uniform temperature and opacity distribution.

Improvement on this assumption would require spatially resolved continuum observations that sample the peak of the dust SED at rest frame ~100 μm (which as of this date have been approved, but not observed).

## 3.2. Literature Sample

We draw our sample of sources detected at both low and high redshift, in p-H<sub>2</sub>O(2<sub>0,2</sub> − 1<sub>1,1</sub>), from the literature. The local sample is drawn from the Herschel Science Archive (Yang et al. 2013) and the AGNs in this sample are identified in Koss et al. (2013). The FIR luminosity in the local galaxies is estimated by fitting 60 and 100 μm photometry (Sanders et al. 2003) with a modified blackbody as discussed in Section 3.1. We fix β = 2.0 to be consistent with the SPT sources and λ<sub>o</sub> to 100 μm (Draine 2006), as there are only two photometry points available.

The high-redshift ULIRGs include SPT0346-52 (Apostolovski et al. 2019), SPT0125-47 (Appendix A2 B), HFLS3 (Riechers et al. 2013), APM08279+5255 (van der Werf et al. 2011), G12. v2.30, NBV1.78, SDP17b, NAV1.195, and SDP11 (Omont et al. 2013; Yang et al. 2016). L<sub>FIR</sub> in HFLS3 is estimated using photometry taken from Riechers et al. (2013) and magnification from Cooray et al. (2014). In the quasar APM08279+5255, L<sub>FIR</sub> is taken from Beelen et al. (2006), Weiß et al. (2007) and is magnification from Riechers et al. (2009). L<sub>FIR</sub> in all the other sources (except SPT sources) are estimated with the same process used for the local galaxies (fixing β and λ<sub>o</sub>) using 250, 350, 500, and 880 μm photometry and magnification from Bussmann et al. (2013). All the high-redshift ULIRG values are given in Table 5.

## 3.3. Spectral Analysis and L<sub>H<sub>2</sub>O</sub>

To obtain the spatially integrated flux density in each of the 50 or 100 kms<sup>-1</sup> wide channels, we apply the continuum mask to each velocity bin of the spectral cube and integrate within the region selected. Selecting a mask based on the continuum flux is less biased because of the high signal-to-noise ratio in every pixel in the continuum, unlike the moment 0 map. This mask includes 95%–100% of the total water emission. The resultant spectra are shown in Figure 2. We use the standard deviation of flux density in line-free channels as the error in each velocity bin. To obtain the line properties, we use the non-parametric estimation of line width described in Bothwell et al. (2013b), which is a more preferred method than fitting a simple Gaussian profile, especially if the line profile is asymmetric. The intensity-weighted second moment of the spectrum is
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left (A): global $L_{H_2O}/L_{FIR}$ plotted as a function of spatially integrated intrinsic $L_{FIR}$ (corrected for magnification in high-redshift lensed galaxies). $L_{H_2O}$ of the low-redshift LIRGs and ULIRGs (with mild and strong AGNs) are from Yang et al. (2013). $L_{H_2O}$ of the high-redshift ULIRGs and AGNs are taken from Omont et al. (2013), Yang et al. (2016), and van der Werf et al. (2011). The three SPT sources, SPT0529-54, SPT0532-50, and SPT0538-54 and the Cloverleaf are presented in this paper. SPT0346-52 is taken from Apostolovski et al. (2019) and SPT0125-47 is presented in Appendix B. The dotted-dashed line is the best fit to all the sources by allowing the slope as a free parameter. The best fit by fixing the slope to zero is shown as a thick black line and the gray region corresponds to the error on the fit. The dashed line is a fit to the high-redshift sources and the dotted line fits the low-redshift galaxies with a fixed slope of zero. Right (B): resolved $L_{H_2O}/L_{FIR}$ plotted as a function of surface brightness in units of $L_{\odot} \text{ kpc}^{-2}$. Each data point is the value of pixels binned within $0.05 \times 10^{12} L_{\odot} \text{ kpc}^{-2}$. The black data points are obtained by combining all the pixels from the five sources and the fit to these points by fixing the slope to zero is shown as a thick black line with the $1\sigma$ uncertainty on the fit shown in gray. The dotted line is the best fit by allowing the slope to vary and it is within the gray error region. As shown in the plots, $L_{H_2O}$ is strongly correlated with $L_{FIR}$ both at global and resolved scales within the galaxy.</figcaption></figure>

given by

$$ \sigma_{\nu} = \frac{\int (\nu - \bar{\nu})^2 S_{\nu} d\nu}{\int S_{\nu} d\nu}, \qquad (3) $$

where $\bar{\nu}$ is the intensity-weighted frequency centroid, $S_{\nu}$ is the integrated flux at frequency $\nu$, and the FWHM of the line is estimated as FWHM ~2.35$\sigma_{nu}$. The line properties estimated with this method are listed in Table 3 along with the velocity-integrated line flux ($I_{H_2O}$ in Jy kms$^{-1}$). We obtain the line luminosity from the relation given in Solomon & Vanden Bout (2005):

$$ L_{H_2O} = (1.04 \times 10^{-3}) I_{H_2O} \nu_{\text{rest}} D_L^2 (1+z)^{-1}, \qquad (4) $$

where $L_{H_2O}$ is the total line luminosity in units of $L_{\odot}$, $\nu_{\text{rest}}$ is the rest frequency of the line in GHz (987.927 GHz for the p-H$_2$O(2$_{0,2}$ - 1$_{1,1}$) transition) and $D_L$ is the luminosity distance to the source at a redshift $z$ in Mpc. We estimate $L_{H_2O}$ in each pixel using the above equation with $I_{H_2O}^i$ (the integrated line flux of the $i$th pixel taken from the moment 0 image).

# 4. Analysis and Discussion

## 4.1. $L_{H_2O}-L_{FIR}$ Correlation and SFR Calibration

Using the estimated intrinsic luminosities (corrected for magnification) of FIR and H$_2$O, we plot $L_{H_2O}/L_{FIR}$ as a function of $L_{FIR}$ in Figure 3 in both globally integrated scales and spatially resolved scales (values for high-redshift galaxies are given in Tables 3 and 5 and local galaxies are discussed in Yang et al. 2013). In our analysis, we assume that H$_2$O and FIR are spatially lensed, hence the issue of differential lensing where the lensing magnification varies across the source (Blain 1999; Hezaveh et al. 2012) is not significant.

From Figure 3(A), the correlation between global $L_{H_2O}$ and $L_{FIR}$ is slightly super-linear and shown as a dotted-dashed black line. A fit (including all sources from the literature in the MCMC) to $\log_{10}(L_{H_2O})$ and $\log_{10}(L_{FIR})$ gives

$$ L_{H_2O} \propto L_{FIR}^{(1.17 \pm 0.21)}. \qquad (5) $$

This result is similar to the conclusion presented in the literature (Omont et al. 2013; Yang et al. 2016) where $L_{H_2O}/L_{FIR}$ is slightly higher in luminous high-redshift galaxies when compared to less luminous local galaxies. The fit to $L_{H_2O}/L_{FIR}$ with the slope fixed to zero is shown as a thick black line (the error region is shown in gray) gives $L_{H_2O}/L_{FIR} = 1.69_{-0.54}^{+0.79} \times 10^{-5}$. The fit (slope fixed to zero) to high-redshift sources is shown as a dashed line and the fit to the low-redshift galaxies is shown as a dotted line. Both the lines (with their respective errors) are within the error bar of the thick black line, which shows that the increase in $L_{H_2O}/L_{FIR}$ in this sample might not be significant. The almost linear relation between $L_{H_2O}$ and $L_{FIR}$ over more than three orders of magnitude supports the previously found conclusions in the literature that p-H$_2$O(2$_{0,2}$ - 1$_{1,1}$) traces $L_{FIR}$ both in local galaxies and in high-redshift ULIRGs irrespective of the presence of AGNs.

It should be noted that the low-redshift sample from Yang et al. (2013) is not complete and the super-linear correlation between $L_{H_2O}$ and $L_{FIR}$ seen in Figure 3(A) might be real. The deficiency of H$_2$O in less luminous galaxies ($L_{FIR} < 10^{11.5} L_{\odot}$) has also been observed in the compilation of fluxes presented by Liu et al. (2017), where a few LIRGs have $L_{H_2O}$ similar to ULIRGs while the others either have lower or no H$_2$O detection. This effect could be arising because of p-H$_2$O(2$_{0,2}$ - 1$_{1,1}$) excitation requirements. H$_2$O molecules have to be well shielded from UV radiation to avoid dissociation but also have to reside in warm gas (not UV heated) to escape into the gas phase from
<figure><img src="image_6.png" /><figcaption>Figure 4. Normalized deviation of $L_{H_2O}/L_{FIR}$ in each pixel from the median value, i.e., $\Delta \frac{L_{H_2O}}{L_{FIR}}$, for each source. The contours in black correspond to continuum emission at [3, 5, 10, 20, 40, 80 ...] $\times \sigma$, where $\sigma$ is the rms noise in the continuum map. Values around zero are closer to the median value in that source. The deviation in all the sources is within 10%.</figcaption></figure>

grain mantles. Moreover, p-H₂O(2,0,2 − 1,1) requires dense gas to populate the 1,1 level through collisions, which is the base for 101 µm excitations. Hence, H₂O couples with far-infrared radiation strongly in such warm, dense, well-shielded gas, which is prevalent in galaxies with $L_{FIR} > 10^{11.5} L_⊙$, but less so in low luminous galaxies. H₂O emission is also enhanced as a result of shocks or intense radiation fields prominent in starbursts. Shocks could increase the abundance of H₂O and strong radiation could lead to an increase in excitation (e.g., Gonzalez et al. 2010; Omont et al. 2011). This might result in the slightly super-linear correlation.

Figure 3(B) shows the correlation between $L_{FIR}$ and $L_{H_2O}$ on resolved scales. In each spatial pixel, the surface brightness ($\Sigma_{L_{FIR}}$) is estimated by dividing the observed $L_{FIR}$ in each pixel by the area of the pixel in the image plane. Because lensing conserves surface brightness, the pixel-by-pixel values of $\Sigma_{L_{FIR}}$ do not require a lensing correction. $L_{H_2O}/L_{FIR}$ per pixel also does not require a lensing correction, as we assume that the FIR and H₂O emission are spatially lensed. The data points corresponding to each source are the average value of $L_{H_2O}/L_{FIR}$ in pixels binned within $(0.05 \times 10^{12}) L_⊙$ kpc⁻² and the error on each data point corresponds to the combined error of the standard deviation of pixel values within that bin and the propagated error due to averaging. This averaging reduces the number of degenerate pixels arising from multiple images of the same region in the source due to gravitational lensing. However, pixel averaging contributes to the error on $\Sigma_{L_{FIR}}$ and $L_{H_2O}/L_{FIR}$ because of differential lensing within pixels in each bin. These errors are not shown in Figure 3(B), as it is difficult to quantify without a good lens model, but this does not affect the results discussed below. The black points are obtained by combining all the pixels from all the sources. The best fit to the black points by fixing the slope to zero is given by

$$ \frac{L_{H_2O}}{L_{FIR}} = 2.76^{+2.15}_{-1.21} \times 10^{-5}. \qquad (6) $$

This fit is shown as a thick black line and the gray shaded region in Figure 3(B), which is consistent with the global $L_{H_2O}/L_{FIR}$ (Figure 3(A)). The dotted line in Figure 3(B) is the best fit to the resolved data by allowing the slope to vary. We obtain $L_{H_2O}/L_{FIR} \propto L_{FIR}^{0.14 \pm 1.44}$, which is consistent within the gray shaded region. Comparing the two fits using a F-test, we conclude that the fit with varying slope does not provide

significantly more information than the fit with a slope fixed to zero. We obtain an F-distribution value of ~0.5 with the null hypothesis that the fit with varying slope provides more information than the fit with a slope fixed to zero and we reject this hypothesis if the F-distribution value is >0.05.

We observe that $L_{H_2O}/L_{FIR}$ remains linear in SPT0532-50, SPT0538-50, and the Cloverleaf with $\Sigma_{L_{FIR}}$. The correlation is not obvious in SPT0529-54, as the signal-to-noise ratio of the H₂O emission in this source is lower than those of other sources (peak SNR ~7.5 in 100 kms⁻¹ channels). This strong correlation between $L_{H_2O}$ and $\Sigma_{L_{FIR}}$ suggests that H₂O is tracing SFR not just at global scales (as discussed in the previous paragraph) but also at resolved scales (~1 kiloparsec) within the galaxy. Thus, this result demonstrates that we can use resolved H₂O as a resolved SFR indicator in high-redshift intense star-forming regions. However, high-resolution continuum observations in the FIR regime are required to quantify temperature (and hence $L_{IR}$) variations within the galaxy.

To understand which regions of the galaxy H₂O best traces $L_{FIR}$ in the image plane at resolved pixel scales, we plot the relative deviation of $L_{H_2O}/L_{FIR}$ in Figure 4. Here, $\Delta \frac{L_{H_2O}}{L_{FIR}}$ corresponds to [ratio¹ - ratiomed]/ratiomed where ratio¹ is $L_{H_2O}/L_{FIR}$ of the ith pixel and ratiomed is the median value of $L_{H_2O}/L_{FIR}$ in each source. The closer the $\Delta \frac{L_{H_2O}}{L_{FIR}}$ value is to 0, the closer it is to the median value, which implies that H₂O is well correlated with FIR in those pixels. As seen in the images, the correlation is stronger in regions with good signal-to-noise and the ones which deviate the most from $\Delta \frac{L_{H_2O}}{L_{FIR}}$ are at the edges with low signal-to-noise. Moreover, the bright regions are multiple images of the same region in the source (due to gravitational lensing) and as expected, we see that $L_{H_2O}/L_{FIR}$ is the same in these regions. This suggests that water emission faithfully traces the FIR luminosity on resolved scales.

The almost linear correlation between $L_{H_2O}$ and $L_{FIR}$ at the resolved scales in the galaxies (Figure 3(B) and Equation (6)) allows us to calibrate the SFR as a function of $L_{H_2O}$ for high-redshift intense star-forming regions. The caveat discussed previously in Section 1 that the SFR calibration from $L_{FIR}$ cannot be applied to certain environments, such as regions around AGNs and resolved SFR—depends on variations in dust temperature and opacity. However, we assume a uniform temperature and opacity distribution across the sources, as we
<table><caption>**Table 4**

Observed Physical Properties in High-redshift ULIRGs</caption><thead><tr><th>Source</th><th>λ<sub>max</sub><br>(μm)</th><th>M<sub>gas</sub><br>(10<sup>10</sup> L<sub>☉</sub>)</th><th>A<sub>eff</sub><br>(kpc<sup>2</sup>)</th><th>Σ<sub>gas</sub><br>(10<sup>10</sup> M<sub>☉</sub> kpc<sup>-2</sup>)</th><th>References</th></tr></thead><tbody><tr><td>SPT0529-54</td><td>108.57 ± 5.59</td><td>4.05 ± 0.90</td><td>9.17 ± 3.46</td><td>0.44 ± 0.19</td><td>Bothwell et al. (2017)</td></tr><tr><td>SPT0532-50</td><td>92.43 ± 6.83</td><td>6.05 ± 1.71</td><td>7.46 ± 1.69</td><td>0.81 ± 0.29</td><td>Bothwell et al. (2017)</td></tr><tr><td>SPT0538-50</td><td>94.91 ± 6.83</td><td>1.7 ± 0.3</td><td>15.71 ± 3.05</td><td>0.11 ± 0.03</td><td>Aravena et al. (2016)</td></tr><tr><td>SPT0125-47</td><td>84.98 ± 11.17</td><td>11.5 ± 1.0</td><td>11.93 ± 9.93</td><td>0.96 ± 0.81</td><td>Aravena et al. (2016)</td></tr><tr><td>SPT0346-52</td><td>73.80 ± 5.59</td><td>8.2 ± 0.6</td><td>2.81 ± 0.45</td><td>2.92 ± 0.52</td><td>Aravena et al. (2016)</td></tr><tr><td>Cloverleaf</td><td>71.32 ± 6.83</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

**Note.** $\lambda_{max}$, the rest-frame wavelength at which dust SED peaks, is estimated from the modified blackbody fit to photometry using the MCMC algorithm. Gas mass ($M_{gas}$) is taken from the references shown in the last column. The Sérsic area ($A_{eff}$) is calculated from the best-fit source parameters from lens modeling (Spilker et al. 2016). $\Sigma_{gas}$ is the gas surface density.

do not have more resolved continuum observations around the peak of the SED. The SFR is generally estimated from the $L_{IR}$ scaling relations discussed in Kennicutt & Evans (2012):

$$ \text{SFR} [M_{\odot} \text{yr}^{-1}] = 1.47 \times 10^{-10} L_{\text{IR}} [L_{\odot}]. \quad (7) $$

To convert $L_{IR}$ to $L_{FIR}$, we use $\langle L_{IR}/L_{FIR} \rangle \sim 1.38$ obtained from the SEDs of SPT0529-54, SPT0532-50, SPT0538-50, and the Cloverleaf. This value is similar to $\langle L_{IR}/L_{FIR} \rangle$ in the low-redshift galaxies, which is $\sim 1.29$ (Soifer et al. 1987). Using the $L_{IR}$ to $L_{FIR}$ conversion from the SPT sources and the relation between $L_{H_2O}$ and $L_{FIR}$ given by Equation (6), we calibrate SFR using $L_{H_2O}$ at resolved galaxy scales as:

$$ \text{SFR} [M_{\odot} \text{yr}^{-1}] = 7.35^{+5.74}_{-3.22} \times 10^{-6} L_{\text{H}_2\text{O}} [L_{\odot}]. \quad (8) $$

We have shown that $L_{H_2O}$ is well correlated with the FIR continuum at resolved scales. As mentioned previously, this calibration is applicable to high-redshift intense star-forming regions assuming no spatial variations in temperature and optical depth. A similar analysis using resolved continuum observations at the peak of the SED has to be performed to obtain a more accurate SFR calibration from $L_{H_2O}$, which could then be used instead of the observationally expensive $L_{IR}$. We note that resolved FIR continuum observations at the peak of the SED (even at these redshifts) are observationally expensive, while $H_2O$ is a bright line and is easily observable in high-redshift sources with ALMA. In addition to providing an alternative to the expensive FIR continuum observations, $H_2O$ additionally provides kinematics of the star-forming regions.

## 4.2. Effect of AGNs

Aside from the SPT sources, which are dominated by star formation and show no evidence of an AGN, our sample also includes the Cloverleaf quasar, a well characterized AGN at $z = 2.56$.

The higher transitions of $H_2O$ ($E_{up} \ge 400$ K) are mainly excited by absorption of short-wavelength far-infrared photons ($\le 50$ µm) emitted by the hot dust surrounding the AGN. Modeling of a lensed quasar, APM 08279+5255 at $z \sim 3.9$ (van der Werf et al. 2011), showed that the higher $H_2O$ transitions are arising from the compact central region with $T_d \sim 200$ K. The AGN contributes less to the $J \le 3$ excitations (mainly excited by 75 µm and 101 µm photons) in the warm regions but does contribute significantly to the total $L_{IR}$. This results in APM 08279+5255 lying low on the $L_{H_2O}/L_{IR}$ correlation

(Figure 3(A), Table 5). Similar results are found in Mrk 231 where $p-H_2O(4_{2,2} - 4_{1,3})$ is detected (van der Werf et al. 2010) and $L_{H_2O}$ in $p-H_2O(2_{0,2} - 1_{1,1})$ and $p-H_2O(2_{1,1} - 2_{0,2})$ is lower than other ULIRGs without AGNs (Yang et al. 2016).

Previous work (Yang et al. 2013) has shown that although the presence of strong AGN lowers the global $L_{H_2O}/L_{IR}$ ratio, there does not appear to be a significant effect of AGN on $H_2O$ emission lines. This can be seen from Figure 3(A). Our spatially resolved analysis of the Cloverleaf quasar in Figure 3(B) indicates that $L_{H_2O}/L_{FIR}$ remains constant at resolved scales and is similar to ULIRGs even in the presence of an AGN. This suggests that the presence of an AGN has little impact on $p-H_2O(2_{0,2} - 1_{1,1})$ excitation—not just at the global scale but also down to kiloparsec scales.

## 4.3. Correlation of $L_{H_2O}/L_{FIR}$ with Physical Properties

The strength of $H_2O$ emission depends on a number of physical properties such as $T_d$, $H_2O$ column density, and continuum opacity ($\tau$), which are better constrained through modeling of multiple $H_2O$ excitations. As we observed only a single transition, we now investigate to what extent the global variations in $L_{H_2O}/L_{FIR}$ we observe correlate with other properties we constrain such as $\lambda_{max}$ (the wavelength in rest frame at which the dust SED peaks) and gas-mass density ($\Sigma_{gas}$). We estimate these by using the available photometry and values from the literature (Table 4).

The correlation of $L_{H_2O}/L_{FIR}$ with $\lambda_{max}$ can be interpreted as a correlation with the dust temperature or $L_{FIR}$ surface density. Because $p-H_2O(2_{0,2} - 1_{1,1})$ is mainly excited by the FIR radiation, we are interested in understanding the correlation with dust temperature. We use $\lambda_{max}$, as it is a more direct observable than dust temperature, which is degenerate with optical thickness (e.g., Papadopoulos et al. 2010). In Figure 5(A), there is no correlation between $L_{H_2O}/L_{FIR}$ and $\lambda_{max}$, consistent with previous results on low-redshift galaxies (Table 2 in Yang et al. 2013), where no relationship is observed between $L_{H_2O}/L_{IR}$ and $S_{60\mu m}/S_{100\mu m}$, a dust temperature indicator (these wavelengths are used to fit the SED in our analysis for local sources).

In Figure 5(B), we plot $L_{H_2O}/L_{FIR}$ as a function of $\Sigma_{gas}$ in an effort to understand whether collisions significantly affect the $H_2O$ excitation. $\Sigma_{gas}$ is calculated by dividing gas mass by the effective area of the source. The gas-mass values from Bothwell et al. (2017; estimated using [CI] observations) in SPT0529-54 and SPT0532-50 and from Aravena et al. (2016) for SPT0538-50, SPT0125-47, and SPT0346-52 (estimated using CO(1–0) observations) are used. The source properties
<figure><img src="image_7.png" /><figcaption>Figure 5. Left (A): correlation between global L<sub>H2O</sub>/L<sub>FIR</sub> and λ<sub>max</sub>, the rest-frame wavelength at which dust SED peaks. Right (B): global L<sub>H2O</sub>/L<sub>FIR</sub> is plotted as a function of gas surface density, Σ<sub>gas</sub> in units of M<sub>☉</sub> kpc<sup>-2</sup>. The gas masses for SPT0529-54 and SPT0532-50 are taken from Bothwell et al. (2017). SPT0538-50 and other two SPT sources (SPT0125-47 and SPT0346-52) are detailed in Aravena et al. (2016). The intrinsic SPT source sizes obtained from lens modeling can be found in Spilker et al. (2016). As shown in the plots, the variation in L<sub>H2O</sub>/L<sub>FIR</sub> is uncorrelated with either the effective temperature of the dust SED or the gas surface density.</figcaption></figure>

are detailed in Spilker et al. (2016), where lens modeling of the 870 µm dust continuum is performed by assuming single or multiple Sérsic source profiles. Using these values, the area ($A_{eff}$) under a Sérsic profile is calculated. This method might overestimate $\Sigma_{gas}$, but because the sizes are within a factor of ~2, the overestimated value might only be by a factor of few. Moreover, the CO (gas) sizes can be larger than that of the infrared emission (Spilker et al. 2015; Tadaki et al. 2017; Calistro Rivera et al. 2018; Dong et al. 2019). For a simple calculation, we assume that the dust and gas sizes are similar. All the values are given in Table 4. As seen in Figure 5(B), there is no observed correlation of $L_{H2O}/L_{FIR}$ with $\Sigma_{gas}$.

The dust opacity at 100 µm, the wavelength at which $p$-H$_2$O(2$_{0,2}$ – 1$_{1,1}$) is excited, could also affect the intensity of the line and could explain the spread in the $L_{H2O}-L_{FIR}$ correlation. The slightly super-linear $L_{H2O}/L_{FIR}$ correlation where the increase in H$_2$O line emission is faster than $L_{FIR}$ (e.g., Omont et al. 2013; Yang et al. 2016) could be because of the increase in $\tau_{100}$ (dust opacity at 100 µm) with the increase in $L_{FIR}$ (González-Alfonso et al. 2014), in turn enhancing $L_{H2O}$ because of photon trapping. In a high $\tau_{100}$ medium, the 100 µm photons are trapped and scattered thereby increasing the local radiation field. This amplifies the $p$-H$_2$O(2$_{0,2}$ – 1$_{1,1}$) pumping and hence the $p$-H$_2$O(2$_{0,2}$ – 1$_{1,1}$) line photons. We do a simple estimation of $\tau_{100}$ in the three SPT sources using the equation from Yang et al. (2016) where $\tau_{100}$ is given by

$$ \tau_{100} = \kappa_{100} \frac{M_{\text{dust}}}{2\pi r^2}, \qquad (9) $$

$\kappa_{100}$ is the dust absorption opacity at 100 µm and r is the radius of the source at submillimeter wavelengths. We use $\kappa_{\lambda} = 2.92 \times 10^5 (\lambda/\mu\text{m})^{-2} \text{cm}^2 \text{g}^{-1}$ (Li & Draine 2001) at rest wavelength $\lambda$ and dust mass ($M_{dust}$) given by

$$ M_{\text{dust}} = \mu^{-1} \frac{D_L^2 S_\nu}{(1 + z_s) \kappa_\lambda [B_\nu(\text{T}) - B_\nu(\text{T}_{\text{CMB}})]} \qquad (10) $$

Here, $S_\nu$ is the flux density at observed frequency, $z_s$ is the redshift of the source, $D_L$ is the luminosity distance, $B_\nu(\text{T})$ is

the Planck function at rest frequency (described in Section 3.1), and $\mu$ is the magnification of the source. From the lens model parameters derived in Spilker et al. (2016), we estimate $\tau_{100} \sim 0.34, 1.36$ and $0.46$ for SPT0529-54, SPT0532-50, and SPT0538-54, respectively. The higher value of $\tau_{100}$ in SPT0532-50 could be enhancing the H$_2$O luminosity above the average value. However, several physical factors like $T_d$, opacity, H$_2$O abundance, etc., can also influence the intensity of the line.

To summarize, the global variations in $L_{H2O}/L_{FIR}$ are not observed to be affected by the physical properties of the galaxy such as $\lambda_{max}$ and $\Sigma_{gas}$. Large dust opacity at 100 µm might enhance $L_{H2O}$ due to photon trapping. Modeling with multiple transitions will give a better understanding of the factors influencing the correlation between $L_{FIR}$ and $L_{H2O}$.

### 4.4. H$_2$O and CO

CO(6–5) traces relatively dense gas (with critical density of H$_2$ ~ $10^5$ cm$^{-3}$) in molecular clouds, although not as dense as HCN or HCO$^+$ (Shirley 2015; Béthermin et al. 2018). The high –J CO lines are therefore found to be correlated with the far-infrared field in these star-forming regions (Figure 1 in Liu et al. 2015). Here, we investigate this correlation in the context of the $L_{H2O}-L_{FIR}$ relation. We make use of the spatially and spectrally resolved observations of mid–J CO(6–5) in SPT0529-54, SPT0532-50 (Dong et al. 2019), and supplement these data with observations of the two other SPT sources SPT0346-52 (Apostolovski et al. 2019) and SPT1247-50 (which is not detected in $p$-H$_2$O(2$_{0,2}$ – 1$_{1,1}$)) from Dong et al. (2019). The imaging of the CO data is similar to that described in Section 2.2. The mask used to select the pixels recovers 93%–100% of the CO emission, depending on the source.

Figure 6 shows $L_{CO(6-5)}/L_{FIR}$ as a function of $L_{FIR}$ similar to Figure 3. Figure 6(A) contains the global integrated values in local luminous infrared galaxies from Lu et al. (2017) and high-redshift SPT sources. Figure 6(B) shows the resolved correlation between CO(6–5) and $L_{FIR}$ (the binning procedure
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left (A): global $L_{CO(6-5)}/L_{FIR}$ as a function of $L_{FIR}$. The local LIRGs are shown as yellow data points taken from Lu et al. (2017). The high-redshift ULIRGs are represented by the SPT sources. The thick black line is a fit to all the sources by fixing the slope to zero, with the $1\sigma$ error shown as the gray region. The dotted-dashed line is the fit from allowing the slope to vary. From the plot, we see that the correlation is almost linear. Right (B): resolved $L_{CO(6-5)}/L_{FIR}$ as a function of surface brightness in units of $L_{\odot}~kpc^{-2}$. Each data point is the value of pixels binned within $0.05 \times 10^{12} L_{\odot}~kpc^{-2}$. The fits are to the combined binned pixels shown in black. The correlation within the sources follow a similar pattern as the global values in the left plot and is also nearly linear. This plot and Figure 3 together suggest that H$_2$O is as good a tracer of the far-infrared radiation as CO(6-5).</figcaption></figure>

is similar to that described in Section 4.1), where the resolved $L_{FIR}$ is estimated using the continuum around the CO(6-5) line. CO(6-5) is observed to have an almost linear correlation with $L_{FIR}$ both at global and resolved scales, similar to H$_2$O.

The spectra of CO(6-5) and $p$-H$_2$O($2_{0,2} - 1_{1,1}$) in SPT0532-50 and SPT0346-52 (with a good detection of both the lines) show that CO has an FWHM consistent with H$_2$O within the errors (Figure 8(A)). This may indicate that both the lines are emitted from similar regions in the galaxy (See also Omont et al. 2013; Yang et al. 2016; Liu et al. 2017). It can further be seen from the spatial distribution comparison in the image plane (Figure 8(B)). This agrees with the results in Yang et al. (2019), where they find similar spatial distribution and also similar kinematic structure between CO(6-5) and $p$-H$_2$O($2_{1,1} - 2_{0,2}$) in G09v1.97.

While $p$-H$_2$O($2_{0,2} - 1_{1,1}$) excitation is due to FIR pumping mechanism and depends mainly on the radiation field density, the CO excitation is due to collisions with the H$_2$ molecules. Hence, CO intensity increases with increase in the gas density and temperature (e.g., Narayanan & Krumholz 2014). The mid and high-$J$ CO lines ($J = 6-5$ and above) are shown to have increasingly sublinear slopes with $L_{FIR}$, which suggests that denser and much warmer gas than star-forming regions (possibly due to shocks or turbulence) is needed for CO excitation (Narayanan et al. 2008; Greve et al. 2014). Thus, variations in star formation efficiency, or other physical properties within the molecular gas such as the shape of the density probability distribution function (PDF) and the median density of the gas within and between galaxies might affect $L_{CO(6-5)}/L_{FIR}$ more strongly than $L_{H_2O}/L_{FIR}$. Although the mid-$J$ $L_{CO}/L_{IR}$ ratio is not expected to be enhanced in galaxies with supernovae or stellar wind driven shocks, NGC 6240 shows a higher ratio (Lu et al. 2017). This suggests that H$_2$O is an intrinsically better tracer of the far-infrared field than CO(6-5). To confirm this result, we need a larger sample of sources across a broad range in $L_{FIR}$ to compare H$_2$O and CO (and other dense gas traces such as HCN) to determine which one is an empirically better tracer of star formation.

## 5. Summary and Conclusion

We observed the $p$-H$_2$O($2_{0,2} - 1_{1,1}$) 987.927 GHz line in SPT0529-54 ($z = 3.369$), SPT0532-50 ($z = 3.399$), and SPT0538-50 ($z = 2.782$) with ALMA. We also include the Cloverleaf quasar at $z = 2.558$ to compare with the star-forming galaxies. The observational results and conclusions from this analysis are as follows:

1. $L_{H_2O}$ is empirically correlated with $L_{FIR}$ over more than three orders of magnitude from low-redshift LIRGs to high-redshift ULIRGs.

2. The relationship between $L_{H_2O}$ and $L_{FIR}$ stays linear even at resolved scales within individual galaxies with average $L_{H_2O}/L_{FIR} = 2.76_{-1.21}^{+2.15} \times 10^{-5}$.

3. This linear correlation holds even in the presence of a strong AGN in the Cloverleaf quasar.

4. We present $p$-H$_2$O($2_{0,2} - 1_{1,1}$) as a resolved SFR calibrator for high-redshift intense star-forming regions, assuming a single temperature and opacity across the source

$$ \text{SFR} [M_{\odot} \text{yr}^{-1}] = 7.35_{-3.22}^{+5.74} \times 10^{-6} L_{H_2O} [L_{\odot}]. $$

5. There is no observed correlation of $L_{H_2O}/L_{FIR}$ with $\lambda_{max}$, the wavelength at which SED peaks or $\Sigma_{gas}$, the gas-mass surface density. The dust opacity at 100 µm ($\tau_{100}$), on the other hand, may influence $L_{H_2O}$ due to photon trapping. However, the current sample is too small to give any definite result.

6. $p$-H$_2$O($2_{0,2} - 1_{1,1}$) is intrinsically a better tracer of $L_{FIR}$ than CO(6-5). A larger sample size is needed to confirm this result.

This work shows that $p$-H$_2$O($2_{0,2} - 1_{1,1}$) traces $L_{FIR}$ at resolved ~kiloparsec scales in high-redshift galaxies with intense star-forming regions while assuming a single temperature and dust opacity across the source. In order to validate these assumptions and obtain a more accurate SFR calibration,
<table><caption>**Table 5**
Observed Properties in High-redshift ULIRGs</caption><thead><tr><th>Source</th><th>z</th><th>μ</th><th>λ<sub>max</sub><br>(μm)</th><th>L<sub>FIR</sub>/μ<br>(10<sup>12</sup> L<sub>☉</sub>)</th><th>L<sub>H<sub>2</sub>O</sub>/μ<br>(10<sup>8</sup> L<sub>☉</sub>)</th><th>L<sub>H<sub>2</sub>O</sub>/L<sub>FIR</sub><br>(10<sup>-5</sup>)</th><th>References</th></tr></thead><tbody><tr><td>SPT0125-47</td><td>2.5148</td><td>5.467 ± 0.120</td><td>84.98 ± 11.17</td><td>19.40 ± 4.37</td><td>5.12 ± 0.87</td><td>2.64 ± 0.74</td><td>Appendix A2 B</td></tr><tr><td>SPT0346-52</td><td>5.6559</td><td>5.570 ± 0.117</td><td>73.80 ± 5.59</td><td>21.50 ± 2.31</td><td>6.36 ± 0.24</td><td>2.96 ± 0.32</td><td>Apostolovski et al. (2019)</td></tr><tr><td>G12.v2.30</td><td>3.259</td><td>9.5 ± 0.6</td><td>82.49 ± 4.34</td><td>8.16 ± 1.02</td><td>1.35 ± 0.27</td><td>1.65 ± 0.39</td><td>Omont et al. (2013)</td></tr><tr><td>NAv1.195</td><td>2.951</td><td>4.1 ± 0.3</td><td>93.67 ± 4.34</td><td>10.25 ± 1.38</td><td>1.63 ± 0.27</td><td>1.59 ± 0.34</td><td>Yang et al. (2016)</td></tr><tr><td>SDP11</td><td>1.786</td><td>10.9 ± 1.3</td><td>91.18 ± 5.59</td><td>2.60 ± 0.63</td><td>0.58 ± 0.12</td><td>2.22 ± 0.71</td><td>Yang et al. (2016)</td></tr><tr><td>NBv1.78</td><td>3.111</td><td>13.0 ± 1.5</td><td>80.01 ± 4.97</td><td>4.65 ± 0.79</td><td>0.94 ± 0.21</td><td>2.02 ± 0.57</td><td>Omont et al. (2013)</td></tr><tr><td>SDP17</td><td>2.305</td><td>4.9 ± 0.7</td><td>92.43 ± 5.59</td><td>7.47 ± 1.65</td><td>1.73 ± 0.32</td><td>2.32 ± 0.67</td><td>Omont et al. (2013)</td></tr><tr><td>HFLS3</td><td>6.337</td><td>2.2 ± 0.3</td><td>76.28 ± 7.45</td><td>13.47 ± 3.74</td><td>5.51 ± 1.12</td><td>4.09 ± 1.41</td><td>Riechers et al. (2013)</td></tr><tr><td>APM08279<br>+5255</td><td>3.9</td><td>4.0</td><td></td><td>50.0 ± 12.0</td><td>6.0 ± 1.2</td><td>1.2 ± 0.3</td><td>van der Werf et al. (2011)</td></tr></tbody></table>

Note. For the SPT sources, z and magnification (μ) are given in Spilker et al. (2016) and λmax (the wavelength in rest frame at which the dust SED peaks) and L<sub>FIR</sub> are estimated by fitting a modified blackbody function to the photometry by fixing β = 2.0. H<sub>2</sub>O observations of SPT0346-52 is discussed in detail in Apostolovski et al. (2019). In HFLS3, photometry is from Riechers et al. (2013) and magnification is from Cooray et al. (2014). The magnification in APM08279+5255 is taken from Riechers et al. (2009). For all other sources, μ and photometry are from Bussmann et al. (2013). L<sub>FIR</sub> and λmax are estimated by fixing β = 2.0 and λ<sub>0</sub> = 100 μm except in APM08279+5255, where L<sub>FIR</sub> is taken from (Beelen et al. 2006; Weiß et al. 2007). L<sub>H<sub>2</sub>O</sub> is taken from the references given in the last column.

we need resolved continuum observations around the peak of the SED. We also need to perform a similar analysis on less luminous galaxies (L<sub>FIR</sub> < 10<sup>12</sup> L<sub>☉</sub>) to extend the SFR calibration. Future work will involve detailed lens modeling of the sources with a pixelated lens model (Hezaveh et al. 2016). In the future, it would also be interesting to compare and model multiple resolved H<sub>2</sub>O lines with other dense gas tracers.

This paper makes use of the following ALMA data: ADS/JAO. ALMA #2015.1.01578.S, #2016.1.01554.S, #2012.1.00844.S, #2012.1.00175.S, and #2011.0.00957.S. The SPT is supported by the NSF through grant PLR-1248097, with partial support through PHY-1125897, the Kavli Foundation, and the Gordon and Betty Moore Foundation grant GBMF 947. S.J., J.D.V., D.P.M., and K.C.L. acknowledge support from the US NSF under grants AST-1715213 and AST-1716127. S.J. and K.C.L. acknowledge support from the US NSF NRAO under grants SOSPA5-001 and SOSPA4-007, respectively. J.D.V. acknowledges support from an A. P. Sloan Foundation Fellowship. D.N. was supported in part by NSF Award AST-1715206 and HST Theory Award 15043.0001. 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.

## Appendix A
## Appendix: A1

The high redshift ULIRGs sample taken from the literature is shown in Table 5.

## Appendix B
## Appendix: A2 APEX Observations of H<sub>2</sub>O in SPT0125-47

We observed p-H<sub>2</sub>O(2<sub>0.2</sub> − l<sub>1,1</sub>) (ν<sub>rest</sub> = 987.927 GHz) line in SPT0125-47 at z = 2.5148 using the APEX-2 receiver of the Swedish Heterodyne Facility Instrument (Vassilev et al. 2008) on the Atacama Pathfinder Experiment (APEX). The observations in the shared ESO+Swedish project 092.A-0467 (PI M. Aravena) were done between July and November 2013 in excellent conditions with precipitable water vapor 0.25 < PWV < 0.5 mm,

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spatially integrated spectrum of the p-H<sub>2</sub>O(2<sub>0.2</sub> − l<sub>1,1</sub>) transition with a 75 kms<sup>-1</sup> spectral resolution in SPT0125-47.</figcaption></figure>

and a total on-source integration time of 3 hours. We reduced the data using the standard procedures in the IRAM CLASS software. The line is clearly detected (Figure 7), with a line flux of 21.8 ± 3.7 Jykm s<sup>-1</sup> and a line width of ~117 km s<sup>-1</sup>. Note that the source is unresolved in the 280 GHz APEX beam of 22''/3.

## Appendix C
## Appendix: A3 Spatial distribution of CO(6–5) and H<sub>2</sub>O in SPT0532-50 and SPT0346-52

Figure 8 shows the spatial distribution of CO(6–5) and p-H<sub>2</sub>O(2<sub>0.2</sub> − l<sub>1,1</sub>) in two high-redshift ULIRGs, SPT0532-50 and SPT0346-52. From the spatial distribution comparison we can see that both the lines are tracing similar regions in the velocity space. Although the source is gravitationally lensed, the similar distribution in the image plane might indicate that they are tracing the same regions of the galaxy in the source plane.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left (A): spectra of CO(6–5) and *p*-H₂O(2₀,₂ – *l*₁,₁) in SPT0532-50 and SPT0346-52 integrated over 50 kms⁻¹ channels. Both CO and H₂O have similar FWHMs. Right (B): the channel map of *p*-H₂O(2₀,₂ – *l*₁,₁) in the background and CO(6–5) contours overlaid in black in both the sources. The contours are at [3, 5, 10, 15 ...] × σ, where σ is the rms noise in the CO cube. CO(6–5) and H₂O may be tracing similar warm dense regions of the galaxy.</figcaption></figure>

## ORCID iDs

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

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

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

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

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

Carlos de Breuck https://orcid.org/0000-0002-6637-3315

Chian-Chou Chen https://orcid.org/0000-0002-3805-0789

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

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

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

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

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

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

Paul Van der Werf https://orcid.org/0000-0001-5434-5942

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

## References

Apostolovski, Y., Aravena, M., Anguita, T., et al. 2019, arXiv:1905.12738

Aravena, M., Murphy, E. J., Aguirre, J. E., et al. 2013, MNRAS, 433, 498

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

Beelen, A., Cox, P., Benford, D. J., et al. 2006, ApJ, 642, 694

Béthermin, M., Greve, T. R., De Breuck, C., et al. 2018, A&A, 620, A115

Blain, A. W. 1999, MNRAS, 304, 669

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

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

Bothwell, M. S., Aguirre, J. E., Chapman, S. C., et al. 2013a, ApJ, 779, 67

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

Bradford, C. M., Aguirre, J. E., Aikin, R., et al. 2009, ApJ, 705, 112
Bussmann, R. S., Pérez-Fournon, I., Amber, S., et al. 2013, ApJ, 779, 25

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

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

Casey, C. M., Cooray, A., Killi, M., et al. 2017, ApJ, 840, 101

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

Combes, F., Rex, M., Rawle, T. D., et al. 2012, A&A, 538, L4

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

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

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

Draine, B. T. 2006, ApJ, 636, 1114

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

Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271

Gonzalez, A. H., Papovich, C., Bradač, M., & Jones, C. 2010, ApJ, 720, 245

González-Alfonso, E., Fischer, J., Aalto, S., & Falstad, N. 2014, A&A, 567, A91

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

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

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

Greve, T. R., 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

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

Hayward, C. C., Lanz, L., Ashby, M. L. N., et al. 2014, MNRAS, 445, 1598

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

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

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

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

Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189

Kennicutt, R. C., Jr., Hao, C.-N., Calzetti, D., et al. 2009, ApJ, 703, 1672

Koss, M., Mushotzky, R., Baumgartner, W., et al. 2013, ApJL, 765, L26

Lagache, G., Cousin, M., & Chatzikos, M. 2018, A&A, 609, A130

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Liu, D., Gao, Y., Isaak, K., et al. 2015, ApJL, 810, L14

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

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

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

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

Mocanu, L. M., Crawford, T. M., Vieira, J. D., et al. 2013, ApJ, 779, 61

Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67

Narayanan, D., Cox, T. J., Shirley, Y., et al. 2008, ApJ, 684, 996

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

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

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

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

Omont, A. 2007, RPh, 70, 1099

Omont, A., Neri, R., Cox, P., et al. 2011, A&A, 530, L3

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

Papadopoulos, P. P., van der Werf, P., Isaak, K., & Xilouris, E. M. 2010, ApJ, 715, 775

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

Rangwala, N., Maloney, P. R., Glenn, J., et al. 2011, ApJ, 743, 94

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

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

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

Shirley, Y. L. 2015, PASP, 127, 299

Soifer, B. T., Neugebauer, G., & Houck, J. R. 1987, ARA&A, 25, 187

Solomon, P., Vanden Bout, P., Carilli, C., & Guelin, M. 2003, Natur, 426, 636

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

Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677

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

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

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

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

Tadaki, K.-i., Kodama, T., Nelson, E. J., et al. 2017, ApJL, 841, L25

Tielens, A. G. G. M. 2013, RvMP, 85, 1021

van der Werf, P. P., Berciano Alba, A., Spaans, M., et al. 2011, ApJL, 741, L38

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

Vassilev, V., Meledin, D., Lapkin, I., et al. 2008, A&A, 490, 1157

Venturini, S., & Solomon, P. M. 2003, ApJ, 590, 740

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

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

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

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

Weiß, A., Requena-Torres, M. A., Güsten, R., et al. 2010, A&A, 521, L1

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

Yang, C., Gavazzi, R., Beelen, A., et al. 2019, A&A, 624, A138

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

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

Zhang, Z.-Y., Gao, Y., Henkel, C., et al. 2014, ApJL, 784, L31