1,1 (ν rest =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 FIR ). 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 -1 ) observations of - ( )

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

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 2 O, the third most abundant molecule in the warm dense ISM after H 2 and CO (Neufeld et al. 1995). As an asymmetric rotor with a large electric dipole moment, H 2 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 ultraluminous 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 2 O emission to be ubiquitous with intensities as bright as CO lines. Modeling has shown that, in addition to infrared pumping where H 2 O is excited by FIR photons, collisions also contribute to the intensities of lowexcitation 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 2 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 1 2 0,2 1,1

(E up =100.8 K,ν rest =987.927 GHz), - ( ) p H O 2 2 2 1,1 0,2 (E up = 137 K, ν rest =752.033 GHz), and p-H 2 O (2 2,0 -2 1,1 ) (E up = 196 K, ν rest =1228.789 GHz) are pumped by 101 μm photons from the base 1 1,1 level and are primarily excited in the warm regions of the galaxy. The collisional excitation of the low-lying levels (1 1,1 and 2 0,2 ) in optically thin or high-density hot regions might also contribute to the emission of the - ( )

The Astrophysical Journal, 880:92 (13pp), 2019 August 1 https://doi.org/10.3847/1538-4357/ab290d line. Hence, H 2 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, groundbased observations of H 2 O in the local universe are nearly always impossible. The Herschel Space Observatory opened the window to multiple H 2 O transitions in the local universe (e.g., Weiß et al. 2010) and Yang et al. (2013) demonstrated that the luminosity of submillimeter H 2 O lines (L H O 2 ) 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 2 O transitions, especially - ( ) p H O 2 1 2 0,2 1,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 2 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 2 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 λ rest ∼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

1,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 theL FIR based calibration holds true). Figure 1 summarizes some of the measurements from the literature and shows that -

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

1,1 is excited by FIR photons, which makes H 2 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)

1,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 O 2 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), 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 O 2 traces far-infrared radiation not just at the integrated global scale (Omont et al. 2013;Yang et al. 2013Yang et al. , 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

1 , 1 emission is described in detail in Figure 3. The H 2 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 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 -

1 , 1 is almost linearly correlated with L FIR especially for L FIR >10 11.5 L e . 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 and3, 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 e 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 H O 2 and L 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 H O 2 -L FIR correlation and the effect of physical properties on this correlation. We conclude with a summary in Section 5.

. We use a Planck 2015 flat ΛCDM cosmology where h=0.677, Ω m =0.307, and Ω Λ =0.693 (Planck Collaboration et al. 2016). We estimate the total infrared luminosity (L IR ) as flux-integrated from 8 to 1000 μm and total far-infrared luminosity (L FIR ) from 42.5 to 122.5 μm in the rest frame (Helou et al. 1985).

We choose -

1,1 as it is one of the brightest H 2 O transitions and has been observed to be well correlated with L 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∼2. 78-3.37). We observed the 1.

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.

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 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). n obs line is the observed frequency of the -

1 , 1 transition at a 987.927 GHz rest frequency. n obs cont is the frequency of the continuum. The beam size and the sensitivity per beam in the continuum map (σ cont ) at ν obs are shown in columns 8 and 9. The peak signal-to-noise of the H 2 O line (SNR 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.

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 velocityintegrated 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.

To fit the SED (and estimate the total L FIR ) in the SPT sources, we use the unresolved millimeter and submillimeter

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.) 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:

where ν is the rest frequency, Ω is the source solid angle, B ν (T d ) is the Planck function estimated at dust temperature T d , and τ is the optical depth. At long wavelengths, τ is given by

with λ o 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 d 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) 3 ), T d , and λ o vary. To investigate L L H O FIR 2 at resolved scales, the L FIR for each pixel is obtained by scaling the L FIR of the entire source using the flux contribution of each pixel to the total continuum flux, i.e.,

where i denotes ith 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 FIR 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).

We draw our sample of sources detected at both low and high redshift, in -

1,1 , 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 λ 0 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 FIR 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 FIR is taken from Beelen et al. (2006), Weiß et al. (2007) and is magnification from Riechers et al. (2009). L FIR in all the other sources (except SPT sources) are estimated with the same process used for the local galaxies (fixing β and λ 0 ) 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.

To obtain the spatially integrated flux density in each of the 50 or 100 kms -1 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 nonparametric 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 given by

where n is the intensity-weighted frequency centroid, S ν is the integrated flux at frequency ν, and the FWHM of the line is estimated as FWHM ∼2.35σ nu . The line properties estimated with this method are listed in Table 3 along with the velocityintegrated line flux (I H O 2 in Jy kms -1 ). We obtain the line luminosity from the relation given in Solomon & Vanden Bout (2005):

where L H O 2 is the total line luminosity in units of L e , ν rest is the rest frequency of the line in GHz (987.927 GHz for the

1,1 transition) and D L is the luminosity distance to the source at a redshift z in Mpc. We estimate L H O 2 in each pixel using the above equation with I i H O 2 (the integrated line flux of the ith pixel taken from the moment 0 image).

-L L

2

Using the estimated intrinsic luminosities (corrected for magnification) of FIR and H 2 O, we plot L L

2 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 and5 and local galaxies are discussed in Yang et al. 2013). In our analysis, we assume that H 2 O and FIR are cospatially 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 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 e ) has also been observed in the compilation of fluxes presented by Liu et al. (2017) 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 L

2 plotted as a function of surface brightness in units of L e kpc -2 . Each data point is the value of pixels binned within 0.05 × 10 12 L e 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σ 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 O 2 is strongly correlated with L FIR both at global and resolved scales within the galaxy. grain mantles. Moreover, -

1,1 requires dense gas to populate the 1 1,1 level through collisions, which is the base for 101 μm excitations. Hence, H 2 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 e , but less so in low luminous galaxies. H 2 O emission is also enhanced as a result of shocks or intense radiation fields prominent in starbursts. Shocks could increase the abundance of H 2 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. 14±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 L

2 remains linear in SPT0532-50, SPT0538-50, and the Cloverleaf with S L FIR . The correlation is not obvious in SPT0529-54, as the signal-to-noise ratio of the H 2 O emission in this source is lower than those of other sources (peak SNR ∼7.5 in 100 kms -1 channels). This strong correlation between L H O 2 and S L FIR suggests that H 2 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 2 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 value is to 0, the closer it is to the median value, which implies that H 2 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

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 L

2 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 O 2 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 O 2 for highredshift 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 , for each source. The contours in black correspond to continuum emission at [3, 5, 10, 20, 40, 80 ...] × σ, where σ 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%.

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):

To convert L IR to L FIR , we use á L IR /L FIR ñ ~1.38 obtained from the SEDs of SPT0529-54, SPT0532-50, SPT0538-50, and the Cloverleaf. This value is similar to á L IR /L FIR ñ in the lowredshift galaxies, which is ∼1.29 (Soifer et al. 1987). Using the L IR to L FIR conversion from the SPT sources and the relation between L H O 2 and L FIR given by Equation (6), we calibrate SFR using L H O 2 at resolved galaxy scales as: 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 O 2 , 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 2 O is a bright line and is easily observable in highredshift sources with ALMA. In addition to providing an alternative to the expensive FIR continuum observations, H 2 O additionally provides kinematics of the star-forming regions.

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 2 O (E up 400 K) are mainly excited by absorption of short-wavelength far-infrared photons (50 μm) emitted by the hot dust surrounding the AGN. Modeling of a lensed quasar, APM 08279+5255 at z ∼ 3.9 (van der Werf et al. 2011), showed that the higher H 2 O transitions are arising from the compact central region with T d ∼ 200 K. The AGN contributes less to the J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 L ratio, there does not appear to be a significant effect of AGN on H 2 O emission lines. This can be seen from Figure 3(A). Our spatially resolved analysis of the Cloverleaf quasar in Figure 3

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 -

1,1 excitation-not just at the global scale but also down to kiloparsec scales.

H we observe correlate with other properties we constrain such as λ max (the wavelength in rest frame at which the dust SED peaks) and gas-mass density (Σ gas ). We estimate these by using the available photometry and values from the literature (Table 4).

The correlation of L L

2 with λ max can be interpreted as a correlation with the dust temperature or L FIR surface density. Because

is mainly excited by the FIR radiation, we are interested in understanding the correlation with dust temperature. We use λ 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 L

2 and λ max , consistent with previous results on low-redshift galaxies (Table 2 in Yang et al. 2013), where no relationship is observed between L L H O IR 2 and S 60μm /S 100μ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 L

2 as a function of Σ gas in an effort to understand whether collisions significantly affect the H 2 O excitation. Σ 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 Note. λ 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). Σ gas is the gas surface density.

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 Σ 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 Omont et al. 2013;Yang et al. 2016) could be because of the increase in τ 100 (dust opacity at 100 μm) with the increase in L FIR (González-Alfonso et al. 2014), in turn enhancing L H O 2 because of photon trapping. In a high τ 100 medium, the 100 μm photons are trapped and scattered thereby increasing the local radiation field. This amplifies the -

1,1 line photons. We do a simple estimation of τ 100 in the three SPT sources using the equation from Yang et al. (2016) where τ 100 is given by

κ 100 is the dust absorption opacity at 100 μm and r is the radius of the source at submillimeter wavelengths. We use κ λ =2.92×10 5 (λ/μm) -2 cm 2 g -1 (Li & Draine 2001) at rest wavelength λ and dust mass (M dust ) given by 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 farinfrared field in these star-forming regions (Figure 1 in Liu et al. 2015). Here, we investigate this correlation in the context of the L H O 2 -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 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 highredshift SPT sources. Figure 6(B) shows the resolved correlation between CO(6-5) and L FIR (the binning procedure is plotted as a function of gas surface density, Σ gas in units of M e kpc -2 . 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 L

2 is uncorrelated with either the effective temperature of the dust SED or the gas surface density.

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 -

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 -

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 L H O FIR 2

. 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.

We observed the - ( ) p H O 2 1 2 0,2 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 starforming galaxies. The observational results and conclusions from this analysis are as follows: 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, 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σ error shown as the gray region. The dotteddashed 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 e kpc -2 . Each data point is the value of pixels binned within 0.05 × 10 12 L e 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).

we need resolved continuum observations around the peak of the SED. We also need to perform a similar analysis on less luminous galaxies (L FIR <10 12 L e ) 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 2 O lines with other dense gas tracers.