Full of Orions: a 200-pc mapping of the interstellar medium in the redshift-3 lensed dusty star-forming galaxy SDP.81

We present a sub-kpc resolved study of the interstellar medium properties in SDP.81, a z=3.042 strongly gravitationally lensed dusty star-forming galaxy, based on high-resolution, multi-band ALMA observations of the FIR continuum, CO ladder and the [CII] line. Using a visibility-plane lens modelling code, we achieve a median source-plane resolution of ~200 pc. We use photon-dominated region (PDR) models to infer the physical conditions - far-UV field strength, density, and PDR surface temperature - of the star-forming gas on 200-pc scales, finding a FUV field strength of ~10^3-10^4 G0, gas density of ~10^5 cm^-3 and cloud surface temperatures up to 1500 K, similar to those in the Orion Trapezium region. The [CII] emission is significantly more extended than that FIR continuum: ~50 per cent of [CII] emission arises outside the FIR-bright region. The resolved [CII]/FIR ratio varies by almost 2 dex across the source, down to ~2x10^-4 in the star-forming clumps. The observed [CII]/FIR deficit trend is consistent with thermal saturation of the C+ fine-structure level occupancy at high gas temperatures. We make the source-plane reconstructions of all emission lines and continuum data publicly available.


INTRODUCTION
Dusty star-forming galaxies (DSFGs) with star-formation rates (SFR) of 100-1000 M yr −1 play a key role in the epoch of peak star-forming activity of the Universe, in the z = 2 − 4 redshift range (e.g., Casey, Narayanan & Cooray 2014;Swinbank et al., 2014). The rest-frame far-UV radiation from the young, massive stars is often completely obscured by their large dust reservoirs and re-radiated at restframe far-infrared (FIR) wavelengths.
The key to understanding the intense star-formation in DSFGs is their interstellar medium (ISM) -dense gas and dust in their (giant) molecular clouds (GMCs). In the canonical GMC picture (e.g., Tielens & Hollenbach 1985), the surface of the clouds is exposed to the intense far-UV radiation from the new-born massive O/B-type stars. This results in a E-mail: mrybak@strw.leidenuniv.nl thermal and chemical stratification of the GMC: while deep in the cloud, the gas is in the molecular form (H2, CO); the ionizing flux and gas temperature increase towards the surface, causing CO to dissociate into C and O, with C being ionized into C + in the photon-dominated region (PDRs) at the cloud surface. Crucially, as the depth of individual layers -and the strength of the emission lines associated with individual species -is largely determined by the surface the FUV field strength (G) and the cloud density (n(H)), observing emission from multiple chemical species (e.g., C + , C and O fine-structure and CO rotational transitions, dust continuum) allows the PDR properties to be inferred, typically via forward modelling.
Previous studies of gas and dust properties in DS-FGs have been limited to spatially-integrated quantities, due to their faintness (S850 µm of few mJy) and compact sizes (few arcsec), compared to the resolution available at FIR/mm-wavelengths. Stacey et al. (2010) used unresolved [C ii], CO(2-1)/(1-0) and FIR continuum observations of a sample of twelve z = 1 − 2 galaxies to infer typical G = 10 2.5 − 10 4 G0, n(H) 10 3 − 10 4.5 cm −3 . More recently, the large samples of strongly-lensed DSFGs discovered in FIR/mm-wave surveys by Herschel -H-ATLAS (Eales et al. 2010;Negrello et al. 2010Negrello et al. , 2017 and HerMES (Oliver et al. 2012;Wardlow et al. 2012) -and the South Pole Telescope (Vieira et al. 2013;Spilker et al. 2016) have provided an opportunity to obtain robust continuum and line detections in a fraction of time compared to non-lensed sources. For example, Gullberg et al. (2015) used ground-based FIR, [C ii] and low-J CO emission to infer the ISM properties in a sample of strongly lensed galaxies from the South Pole Telescope sample. Brisbin et al. (2015) used Herschel [O i] 63-µm and ground-based [C ii] and low-J CO observations to infer n(H)= 10 3 − 10 2 cm −3 , G = 10 1 − 10 3 G0 in a sample of eight z=1.5-2.0 DSFGs. Wardlow et al. (2017) used Herschel fine-structure line spectroscopy of strongly lensed DSFGs, obtaining n(H)= 10 1 − 10 3 cm −3 , G = 10 2.2 − 10 4.5 G0. Similarly, Zhang et al. (2018)  The main limitation of deriving PDR properties from spatially-integrated observations is the implicit assumption that the PDR conditions are uniform across the source, i.e., that the FIR continuum and different emission lines are cospatial. However, resolved observations of DSFGs found that low-J CO emission is significantly more extended than the FIR continuum (e.g., Ivison et al. 2011;Riechers et al. 2011;Spilker et al. 2015;Calistro Rivera et al. 2018), and that mid-/high-J CO emission is more compact than low-J CO emission (e.g., Ivison et al. 2011, R15b). Additionally, line ratios in strongly lensed sources can be affected by differential magnification, introducing a significant bias in the inferred source properties, particularly for highly-magnified galaxies (e.g., Serjeant 2012).
Consequently, resolved multi-tracer observations are necessary to robustly infer the PDR conditions in DSFGs. For example, Rybak et al. (2019) have used resolved [C ii], CO(3-2) and FIR continuum ALMA observations to study the properties of the star-forming ISM in the central regions (R ≤ 2 kpc) of two z ∼ 3 DSFGs from the ALESS sample (Hodge et al. 2013;Karim et al. 2013), finding FUV fields in excess of 10 4 G0 and moderately high densities n(H)=10 3.5 −10 4.5 cm −3 . In this regime, the surface temperature of the PDR regions is of the order of few hundred K and the [C ii] emission becomes thermally saturated, resulting in a pronounced [C ii]/FIR deficit (Muñoz & Oh 2016).
In this work, we use high-resolution, multi-Band ALMA observations to map the ISM properties in the strongly lensed z = 3.042 DSFG SDP.81 at 200-pc resolution to address the following questions: • How are different tracers of ISM -FIR continuum, CO and [C ii] emission spatially distributed?
• What is the cooling budget of the molecular gas, and does it vary across the source?
• What are the conditions of the star-forming gas in DS-FGs at z ∼ 3 -in particular, FUV field strength, density and surface temperature of the PDR regions? Do they vary significantly across the source?
This paper is structured as follows: in Section 2, we summarize the ALMA observations analyzed in this work; Section 3 describes the synthesis imaging and lens-modelling of this data. In Section 4, we use resolved FIR continuum, CO and [C ii] maps to investigate the dust properties of SDP.81 (Section 4.1), CO spectral energy distribution and the gas cooling budget (Section 4.2) and the [C ii]/FIR ratio and deficit (Section 4.3). Section 5 then describes the PDR modelling procedure and results, and discusses the resolved ISM conditions in SDP.81 and the physical mechanisms driving the [C ii]/FIR deficit. Finally, Section 6 summarizes the results of this work.
Throughout this paper we use a flat ΛCDM cosmology from Planck Collaboration et al. (2016). Adopting the spectroscopic redshift of z = 3.042 (ALMA Partnership et al. 2015), this translates to a luminosity distance of 26480 Mpc, and an angular scale of 1 arcsec = 7.857 kpc (Wright 2006).

ALMA Band 3 observations
The CO(3-2) line was observed in ALMA Cycle 5, (Project #2016.1.00633). The observations were carried out on 2017 September 9 in ALMA Band 3. While the original project requested a total of 4 hours on-source with a resolution of 0.2 arcsec, the total observations amount to only ∼45 min. The array configuration consisted of 40 12-meter antennas, with baselines extending from 39 to 7550 m, covering angular scales between ∼0.10 and 19 arcsec. The measured precipitable water vapour was 1.5 mm. The frequency setup consisted of four spectral windows (SPWs), centered at 97. 730, 99.626, 85.813 and 87.614 GHz, respectively. Two SPWs were configured in the line mode with 1920 channels of 976.6 kHz width per channels, for a total bandwidth of 1.875 GHz per SPW. The other two SPWs were configured in the continuum mode (channel width 15.625 MHz), with a bandwidth of 2.0 GHz per SPW.

ALMA Band 8 observations
The [C ii] line was observed in ALMA Cycle 4 (Project #2016.1.01093.S). The observations were carried out on 2016 November 14 with 44 12-meter antennas with baselines extending from 13 to 850 m, corresponding to angular scales between ∼0.16 and 10.5 arcsec. The precipitable water vapour ranged between 0.52-0.64 mm. The frequency setup consisted of four SPWs, centered on 468. 35, 458.24, 456.35 and 470.24 GHz, respectively. SPWs 0-2 were configured in a continuum mode with a bandwidth of 2.0 GHz per SPW, while SPW 3 was split into 240 7.812-kHz-wide channels, for a total bandwidth of 1.875 GHz.

ALMA 2014 Long Baseline Campaign observations
As the ALMA Long Baseline Campaign (LBC) observations of SDP.81 were described in detail in ALMA Partnership et al. (2015), here we provide a only brief overview of the data The observations comprised of 5.9, 4.4 and 5.6 hours on-source in Bands 4, 6 and 7 respectively, taken in October 2014. The frequency setup covered the CO(5-4) line in Band 4, CO(8-7) and H2O (202 − 111) lines in Band 6 and CO(10-9) line in Band 7. The frequency setup provided 7.5 GHz bandwidth per Band, with a channel width of 0.976 -1.95 MHz for line-containing spectral windows (SPWs), and 15.6 MHz for the continuum SPWs. The antenna configuration consisted of 22-36 12-meter antennas, with baselines ranging from ∼15 m to ∼15 km, resulting in a synthesized beam size of ∼ 30 mas for Band 7 to ∼55 mas for Band 4 (Briggs weighing, robust = 1). For the CO lines, a strong drop in S/N was observed for baselines longer than 1 Mλ (ALMA Partnership et al. 2015). This resulted in restricting the CO line analysis to baselines ≤ 1 Mλ; here we adopt this uv-distance cut for all the ALMA datasets.

Data preparation and imaging
All the data processing was done with Casa versions 4.7 and 5.0 (McMullin et al. 2007). Figure 1 shows the Cleaned frequency-integrated maps of the Band 8, [C ii], CO(3-2) and CO(10-9) emission, obtained using natural weighting and an outer Gaussian taper of 0.2 arcsec.

ALMA Band 3
The data were calibrated and flagged using the standard ALMA Cycle 4 pipeline; in addition, antennas DA53 and DV22 were manually flagged due to a systematically high antenna temperature. The CO(3-2) line was detected between 85.509 and 85.571 GHz. Due to the low S/N, we frequency-averaged the line-containing channels into a single frequency bin. A time-averaging of 20 s was applied, as in R15a,b. This reduces the size of the problem by a factor of a few without compromising the surface-brightness sensitivity. The Cleaned images ( Figure 1) show CO(3-2) emission along the main arc and in the counterimage, with a peak S/N of 5.6σ. The Cleaned Band 3 continuum map has an rms noise of 10 µJy beam −1 (natural weighting). We do not detect any Band 3 continuum emission from the lensed source at 3σ significance. We detect the active galactic nucleus (AGN) synchrotron emission in the lensing galaxy (rest-frame frequency 111 GHz) with a total flux of 52±10 µJy. The observed flux is consistent with the Tamura et al. (2015) synchrotron spectrum for the foreground AGN within 2σ.

ALMA Band 8
After applying the standard ALMA Cycle 4 calibration, we performed manual flagging in the time-domain to remove several integrations which showed elevated noise levels; a total of 6 minutes of data was flagged. We also flagged channels affected by atmospheric lines at 468.0 and 457.3 GHz.
The The [C ii] absorption from diffuse ionized gas was detected towards some Galactic star-forming regions (e.g., Gerin et al. 2015). At high redshift, Nesvadba et al. (2016) detected a [C ii] absorption towards to z = 3.4 DSFG, which they interpreted as evidence for infalling gas. Looking for a potential [C ii] absorption in the Cleaned cube with a velocity resolution of 40 km s −1 , we do not find any evidence for [C ii] absorption against the 160-µm continuum. Campaign: Bands 6 and 7 continuum, CO(10-9) line

ALMA Long Baseline
To ensure that all the tracers are compared on the same spatial scales, we revisit the ALMA Partnership et al. (2015) Bands 6 and 7 continuum, applying a uv-distance cut of ≤1 Mλ as opposed to the 2-Mλ cut in R15a. For details of the continuum data preparation, we refer the reader to R15a. Finally, we include the CO(10-9) line, which was observed in ALMA Band 7 (see ALMA Partnership et al. 2015 for detailed description and imaging). The CO(10-9) line is detected between 284.733-285.52 GHz; due to the low S/N ratio (peak image-plane SN 5, Fig. 1) we frequencyaveraged the visibilities into a single channel. Although we obtain an integrated CO(10-9) line flux (measured from the Cleaned images) ∼50 per cent higher than reported by ALMA Partnership et al. (2015) (derived from spectral fitting), the two measurements are consistent within 2σ.

Lens modelling
We reconstruct the source-plane surface brightness distribution of individual tracers via lens-modelling, thereby minimizing the differential magnification bias. The additional complication in case of interferometric data arises from the fact that radio-interferometers do not measure directly the sky-plane surface brightness distributions, but the visibility function, which is essentially a sparsely-sampled Fourier transform of the sky. Furthermore, the high angular resolution provided by ALMA requires the source to be reconstructed on a pixellated grid, as simple analytic models for the source (e.g., Sérsic profiles) are no longer sufficient to fully capture the complex structure of the source. Consequently, several visibility-fitting lens-modelling codes aimed primarily at ALMA observations, have been developed (R15a, Hezaveh et al. 2016;Dye et al. 2018).
We perform the lens-modelling using the method first presented in Rybak et al. (2015a), which is an extension of the Bayesian technique of Vegetti & Koopmans (2009) to the radio-interferometric data. To summarize, this technique uses a parametric lens model, while the source is reconstructed on an adaptive triangular grid obtained via the Delaunay tessellation. This setup provides an increased sourceplane resolution in the high-magnification regions. We minimize the following penalty function, which combines χ 2 calculated in the visibility-space with a source-plane regularization to constrain the problem and prevent noise-fitting: where F and L are Fourier-transform (including the primary beam response) and lensing operators, λs and Rs are the regularisation constant and the regularisation matrix for the source surface brightness distribution and C −1 d is the covariance matrix of the observed visibilities. For a detailed discussion of the lens-modelling technique, we refer the reader to Koopmans (2005) and Vegetti & Koopmans (2009).
We adopt the lens model of R15a derived from 0.1-arcsec resolution ALMA Band 6 and 7 continuum observations, which follows a singular isothermal ellipsoid profile with an external shear component. Keeping the lens mass model fixed, we re-optimize for the offset between the lens center and phase-tracking centre (which differ for different observations) and the source regularization parameter λS.
The sky-plane pixel size is set to 25 mas; as the adaptive grid is obtained by casting back every second pixel, the synthesized beam is sub-sampled by a factor of ∼ 4. The primary beam response is approximated by an elliptical Gaussian. The rms noise for each visibility is estimated as the 1σ scatter of the real/imaginary visibilities for a given baseline over a 10-minute interval, rather than from the Sigma column of the ALMA Measurement Set files. We estimate the uncertainty on the reconstructed source-plane surface brightness distribution by drawing 1000 random realizations of the noise map given the σrms for each baseline, solving for the source using the best-fitting lens model and λS, and calculating the scatter per grid element from the resulting source reconstructions.
As the source-plane surface brightness sensitivity depends on the spatial position, we estimate the spatiallyintegrated magnifications for each tracer as follows: we take all the source-plane pixels with S/N ≥1, 2 and 3, forwardlens them into the image-plane and calculate the corresponding magnification and uncertainty for each S/N threshold; and finally estimate the overall magnification as the mean of the magnification factors for different S/N thresholds.

Source-plane reconstructions
We now discuss the source-plane reconstructions for the ALMA Band 8 continuum and the [C ii] and CO(3-2) emission. Given the varying S/N of individual line and continuum observations, we further assess the reliability of the sourceplane reconstruction by performing our lens modelling analysis on mock datasets in Appendix B; these confirm the robustness of our source-plane reconstruction, with the mean surface brightness recovered within 10-15 per cent for most datasets (∼20 per cent for CO(3-2)), comparable to the ALMA flux calibration uncertainty. We therefore consider our reconstructions to be robust. Fig. 2 presents the best-fitting models of the [C ii], 160µm continuum and CO(3-2) emission; the full output of the lens modelling procedure are provided in Appendix A. The source-plane models reproduce the full range of image-plane structures (given the S/N), as demonstrated by the lack of significant residuals in the dirty-image plane. The median source-plane resolution is ∼30 mas, which corresponds to ∼200 parsec; the resolution in the vicinity of the caustic is as high as ∼3 pc. To allow the comparison on a uniform spacial scale, for the analysis in Section 4 we re-sample the reconstructed source onto a Cartesian grid with 200×200 pc pixels.

ALMA Band 8 continuum
The Band 8 continuum emission is observed at the highest S/N; consequently, the small-scale structure can be robustly reconstructed. Similarly to the Bands 6 and 7 continuum, the Band 8 continuum surface-brightness distribution is concentrated into two prominent clumps, corresponding roughly to  Table 1. Observed [C ii], CO(3-2) and (10-9) line and Band 8 flux densities, measured from the naturally-weighted, velocity-averaged Cleaned images. Individual columns list the observed central frequency ν obs , synthesized beam FWHM size and position angle (east of north), rms noise of the Cleaned images σrms, total flux-density Sν and maximum surface brightness S peak , line width at zero intensity ∆v, and line flux Sν ∆v (if applicable). the northern and central parts of the source, surrounded by a fainter emission approximately 3×2 kpc in extent (R15a, Dye et al. 2015). The northern clump is quadruply imaged, which results in a superb S/N and reconstruction fidelity, while the central clump is only doubly imaged. The relative brightness of the clumps is roughly constant across Bands 6, 7 and 8, indicating only a limited change in the dust spectral energy distribution (Section 4.1).

[CII] emission
Thanks to the high S/N of the [C ii] observations, we are able to model the full velocity structure of the [C ii] emission, rather than just the frequency-averaged line. We frequencyaveraged the line into channels 62.50 MHz (∼ 40 km s −1 ) wide, which provide a good trade-off between S/N per channel (necessary for a robust source reconstruction) and resolving the velocity structure; we then optimize for the source regularization parameter λS for each channel separately. The velocity resolution matches that of the CO(5-4) and CO(8-7) reconstructions from R15b. We calculate the velocity moment-zero (surface brightness) and moment-one maps.
For the moment-one map, we mask pixels with SNR≤1 in each channel before calculating the velocity map.
The [C ii] emission is considerably more extended than both the FIR continuum and the CO(5-4) and CO(8-7) emission ( Fig. 2). Although the [C ii] luminosity increases in the FIR-bright region of the source, at the position of the peak FIR surface brightness, the [C ii] emission shows a significant drop, indicating an extremely low [C ii]/FIR ratio. A similar drop in [C ii]/FIR ratio at the continuum peak was recently found in a z = 1.7 lensed DSFG SDP.11 (Lamarche et al. 2018), in which the regions with the highest ΣSFR are essentially undetected in the [C ii] line emission. We will address the [C ii]/FIR ratio in SDP.81 in more detail in Section 4.3. , and (10-9) line emission (units mJy km s −1 kpc −2 ), and ALMA Bands 6, 7 and 8 continuum (units mJy kpc −2 ), with associated 1σ uncertainty maps; the CO(5-4) and CO(8-7) data are adopted from R15b. In the source reconstruction maps, we masked regions with S/N≤3. Grey lines indicate the lensing caustics. The scattered emission on the edges of the field is a modelling artifact. Note the large extent of the CO and [C ii] emission compared to the FIR-bright source. The source-plane reconstructions and the associated errors are available online.
The reconstructed [C ii] emission map reveals a low surface-brightness emission (Σ [CII] ∼ 3 × 10 7 L kpc −2 ) between -20 and +170 km s −1 , stretching out to ∼10 kpc north of SDP.81 (Fig. 4). The velocity maps show that this excess emission is not co-rotating with the main [C ii]/CO(5-4)/CO(8-7) component. This double-imaged feature is present in the dirty images of the corresponding velocity channels (Fig. A2), clearly offset from the main Einstein arc. The total L [CII] luminosity of this component is ∼ 0.6 × 10 9 L , less than 10 per cent of the total sourceplane [C ii] luminosity. We do not detect any FIR continuum or CO emission from the region outside the dashed box in Fig. 4, despite the high S/N of the FIR and CO(5-4)/CO(8-7) data. This second [C ii] source is not co-spatial with the optically-bright substructure detected in the HST WFC3 160W imaging (Dye et al. 2014, R15b). Given the perturbed velocity/velocity-dispersion fields (R15b), we in-terpret this morphology as a potential low-mass companion; assuming the Mgas scales with L [CII] , the gas mass ratio of the two components is 10:1. Assuming this [C ii] emission is due to star-formation, we use the Herrera-Camus et al.

CO(3-2) emission
The CO(3-2) emission is only detected in the northern part of the FIR-bright region in SDP.81, at ∼ 5σ significance. As the CO(3-2) has much lower S/N (∼ 10) than the remaining lines/continuum data, it must be interpreted carefully. A key check on the relative source-plane distribution of individual tracers is to compare their image-plane surface brightness distributions. As seen in Fig. 1, in the image  plane, the CO(3-2) emission peaks in a different region compared to the Band 8 continuum and [C ii]. We confirm that this is not an artifact of a sparse uv-plane coverage by reconstructing mock ALMA observations matching the S/N and uv-plane coverage of the CO(3-2) data (Appendix B): for a uniformly distributed CO(3-2) emission, emission from the southern part of the counterimage would be detected. This confirms that the CO(3-2) emission is concentrated to the north of the source. In the sky-plane, the CO(1-0) emission peaks to the south of the FIR source (Valtchanov et al. 2011, R15b); which implies a significant offset between the CO(3-2) and CO(1-0) surface brightness peaks.

CO(10-9) emission
The CO(10-9) emission is concentrated between the diamond caustics, and to the north, continuing the trend seen in the CO(8-7) emission. The lower S/N makes the interpretation of the faint extended structure (surface brightness <1.2 mJy km s −1 kpc −2 ) challenging: as we show in Appendix C, this structure can be an artifact of the lensing reconstruction.
Finally, Fig. 3 shows the relative distribution of the FIR continuum, [C ii] and CO emission. This highlights the compact nature of the FIR continuum with respect to the CO and [C ii] lines, and the extended, low-surface brightness [C ii] emission.

Global dust SED
First, we consider the dust thermal spectral energy distribution (SED). Namely we use the ALMA Bands 4/6/7/8 and Herschel PACS/SPIRE observations listed in Table 2 and an optically-thin modified black-body SED: where M dust is the dust mass, h is the Planck constant, k is the Boltzman constant, T dust the dust temperature and β the emissivity index. As the continuum magnification varies only marginally between ALMA Bands 6/7/8, we assume that the (unresolved) Herschel PACS/SPIRE continuum follows the same source-plane surface brightness distribution as ALMA Band 8 continuum, and hence has the same magnification factor. We assume a uniform T dust and β, and adopt a uniform magnification for the 160-µm to 1.3-mm continuum. We obtain T dust = 35 ± 4 K and β = 1.8 ± 0.3 and a (magnification-corrected) FIR luminosity (8-1000 µm) LFIR = (2.4 ± 0.5) × 10 12 L . These are in agreement with previous estimates (Bussmann et al. 2013;Swinbank et al. 2015;Sharda et al. 2018; note that Sharda et al. 2018, do not explicitly fit for β). The inferred T dust , β are in line with those derived for the general high-redshift DSFG population. For example, for the 99 DSFGs from the ALESS sample (Hodge et al. 2013) Swinbank et al. (2014) derived a median T dust = 35 ± 1 K using Herschel SPIRE and ALMA Band 7 photometry (approach directly comparable to ours); while da Cunha et al. (2015) found a median T dust = 43 ± 10 K by including additional UV-to radio-wavelength photometry. The dust emissivity index β is within the range of β = 1.5 − 2.0 inferred for ALESS galaxies (Swinbank et al. 2014).
Although the increasing CMB temperature at high redshifts can significantly bias the inferred T dust and LFIR (da Cunha et al. 2013;Zhang et al. 2016), we consider the CMB effects to be negligible for SDP.81, Namely, adopting the da Cunha et al. (2013) corrections, the CMB effect on T dust is ≤ 0.1 K, while the inferred LFIR is biased by ≤0.1 dex.
We derive the global star-formation rate and the resolved star-formation rate surface density ΣSFR (Fig. 5) using the Kennicutt (1998) SFR-LFIR relation for the Salpeter initial mass function (see also Casey et al. 2014): which yields SFR = 410±90 M yr −1 . Note that a part of the FIR luminosity in intensely star-forming systems might be due to FUV heating contribution from older stellar populations (Narayanan et al. 2015); however, this effect only becomes important at SFRs much higher than that of SDP.81. Given the compact size of the FIR emission in SDP.81 compared to the CO and [C ii] emission, we will refer to the region with ΣSFR ≥ 50 M yr −1 kpc −2 as the FIR-bright region.

Resolved dust SED
Having derived the spatially-integrated T dust and β, we now consider the spatial variation in the dust continuum SED using the source-plane reconstructions of the ALMA Bands 6, 7 and 8 continuum (rest-frame wavelength 160-350 µm). These are well within the Rayleigh-Jeans tail of the modified black-body profile, and hence can not robustly constrain both T dust and β at the same time. We therefore adopt β = 1.8 derived from the spatially-integrated SED fit, and fit for T dust only. However, fitting the modified black-body SED to the resolved Bands 6, 7 and 8 data, we do not find any significant spatial variation in T dust ; the variation in ΣFIR is fully accounted by a varying M dust . Repeating the SED-fitting Table 2. FIR and mm-wave photometry of SDP.81 used for the dust SED modelling, with Herschel PACS and SPIRE fluxdensities adopted from Zhang et al. (2018) and ALMA Bands 6 and 7 flux-densities adopted from R15a. The Herschel flux measurements are given in the image-plane; for our analysis, we assume all the 160 to 500-µm emission to follow the Band 8 surface brightness distribution.

Band
[mJy] on the Cleaned Bands 6, 7 and 8 images (thus eliminating any lens-modelling systematics) yields results consistent with the source-plane analysis. Therefore, for the rest of this paper, we adopt a uniform T dust approximation for the FIR luminosity calculation. The uniform dust temperature in SDP.81 contrasts with the expectation of T dust increasing towards the centre of DS-FGs; for example, Calistro Rivera et al. (2018) invoked a radial decrease in T dust and dust optical depth to explain the different FIR and CO(3-2) sizes in the stacked observations of the ALESS sample, with T dust ≥ 80 K in the central R ≤ 1 kpc region (A. Weiß, private communication).
We therefore consider whether our data might support high T dust (≥60 K) in the central FIR-bright clumps. We do this by fixing T dust to 60 and 80 K, respectively, and optimizing for the normalization (∝ M dust ) only. We find T dust =60 K to be inconsistent with the data at ≥ 3σ significance anywhere except in the vicinity of the FIR brightness maximum (Fig. 6); T dust =80 K is excluded at ≥ 3σ confidence level for the entire source. These results are therefore in tension with the high central T dust invoked by Calistro Rivera et al. (2018).
Finally, we note three fundamental limitations of inferring T dust from the rest-frame 160-350 µm continuum. First, at 200-pc scales probed here, the high-T dust regions might not dominate the continuum SED, resulted in a bias towards low T dust temperature estimates due to the mass-weighting in Equation (2). Second, our resolved ALMA imaging is still on the Rayleigh-Jeans tail of the dust continuum SED; higher-frequency, resolved observations (e.g., with ALMA Bands 9 or 10) will be necessary to properly constrain T dust and β on sub-kpc scales. Finally, the optically thin approximation make break down in the central regions, biasing the inferred T dust .

Emission lines
With the resolved maps of the [C ii] line and the almost complete CO SLED in hand, we can compare the total gas cooling via the [C ii], [O i] and all the CO lines combined. In particular, for the strong FUV fields and high gas (surface) densities expected in DSFGs, as the outer C + layer becomes progressively thinner, the main cooling channel of the neutral gas can switch from the [C ii] line to the [O i] Figure 6. Constraints on dust temperature in SDP.81: reduced χ 2 for T dust = 60 K (right) and 80 K (left). In contrast to temperature gradient models of Calistro Rivera et al. (2018), the FIR-bright clumps are inconsistent with T dust ≥ 80 K (reduced χ 2 ≥ 5).

[CII] dominates the gas cooling budget
We first compare the source-averaged [C ii], [O i] and CO SLED cooling rates. The [O i] 63-µm line is detected in Herschel SPIRE FTS spectra (Valtchanov et al. 2011;Zhang et al. 2018) with a line flux of (2.1 ± 0.7) × 10 −18 W m −2 , while the 3σ upper limit on the [O i] 145-µm line flux is Table 3. Spatially integrated [C ii] and CO line luminosities, with corresponding lensing magnifications µ. We list the source-plane values inferred in this work and Rybak et al., (2015b;R15b), along with CO(1-0), (7-6) and (9-8) luminosities from Valtchanov et al. (2011) and Lupu et al. (2012). To infer the source-plane line luminosities, we assume that CO(1-0) has the same spatial distribution as CO(3-2), CO(7-6) as CO (8-7), and CO(9-8) as CO(10-9). To compare the total gas cooling via the CO ladder and the [C ii], we first obtain a first-order estimate of the total cooling via the CO Jupp = 1 − 10 lines (in the absence of CO(2-1), CO(4-3), CO(6-5) and CO(7-6) observations) assuming that the CO line ratios in SDP.81 match CO ratios of the Class I ULIRGs from Rosenberg et al. (2015). Given the low CO(10-9) luminosity, the Jupp > 10 CO rotational lines are expected to contribute only marginally to the gas cooling. The total CO luminosity is LCO = Σ 10 J=1 L CO(Jupp) 0.8 × 10 9 L , ∼20 per cent of the total L [CII] and ∼10 per cent of the total gas cooling (L [CII] +LCO +L [OI] ). Therefore, the global gas cooling in SDP.81 is dominated by the [C ii] line. Fig. 7 shows the resolved CO/[C ii] gas cooling ratio. The [C ii] line dominates over the CO cooling across the entire source, with the highest CO/[C ii] cooling ratio (∼ 20 per cent) found at the position of the northern star-forming clump. The low CO/[C ii] cooling ratios indicate molecular cloud surface density ≤ 10 3 M pc −2 (c.f., Narayanan & Krumholz 2017).

CO spectral energy distribution
The global source-plane CO SLED in SDP.81 peaks between Jupp = 5 and 9 (in units of L ), with a sharp drop at Jupp = 10 (Fig. 8, see also Yang et al. 2017, for a recent sky-plane analysis). Compared to the Rosenberg et al. (2015) study of CO excitation in present-day ULIRGS, the sharp drop in LCO for Jupp ≥ 9 and minor CO contribution to the global neutral gas cooling makes SDP.81 akin to Rosenberg et al. (2015) Class I sources, which have similarly low CO cooling contributions (≤ 10 per cent). In addition to the heating by FUV radiation, mechanical heating by dissipating supernovae shocks can contribute significantly to the gas energy budget in starburst regions (e.g., Meijerink et al. 2011;Kazandjian et al. 2012Kazandjian et al. , 2015, pumping the Jupp ≥ 10 CO transitions, as well as high-J HCN, HNC and 13 CO lines. Significant mechanical heating is required to explain the high-J CO excitation in some of the archetypal present-day ULIRGs, such as e.g., Arp 220 (Rangwala et al. 2011), NGC 253 (Rosenberg et al. 2014) and NGC 6420 (Meijerink et al. 2013). Given the absence of CO Jupp > 10, HCN or 13 CO observations in SDP.81, we do not attempt to directly constrain the mechanical heating contribution. As the SDP.81 CO ladder is well-reproduced (χ 2 ≤ 1) by the PDRToolbox models which do not include mechanical heating and as the SDP.81-like Rosenberg et al. (2015) Class I ULIRGs are fully compatible with UV heating, the observed line ratios in SDP.81 are consistent with negligible mechanical heating. Finally, we consider how closely individual CO lines trace the FIR emission. We calculate the Pearson's correlation coefficient R for each line and estimate the corresponding uncertainty by drawing 1000 random realizations of CO/[C ii] and FIR surface brightness (with the mean and standard deviation from Fig. 2) for each pixel with SNR≥3 in both tracers. We find R CO(3−2)/FIR = 0.401 ± 0.010, R CO(5−4)/FIR = 0.847 ± 0.006, R CO(8−7)/FIR = 0.821 ± 0.006 and R CO(10−9)/FIR = 0.664 ± 0.09. However, the stronger regularization of CO reconstructions compared to the FIR continuum might artificially decrease the CO-FIR correlation by smoothing the small-scale CO features, especially at high ΣFIR. For comparison, the correlation coefficient for [C ii] R [CII]/FIR = 0.712 ± 0.006

Gas depletion time across the source
Finally, we combine the resolved ΣSFR maps and CO emission line maps to estimate the molecular gas depletion time, t dep ∼ ΣH 2 /ΣSFR. We estimate the molecular gas surface density ΣH 2 using the CO(3-2) line map. For Jupp > 1 CO transitions, the L CO is first down-converted to L CO(1−0) via the RJ1 = L CO(J−J−1) /L CO(1−0) factor, which is then converted to Mgas via Mgas = αCOL CO(1−0) . We adopt R31 from Sharon et al. (2016) who found a mean R31 = 0.78 ± 0.27 for a sample of z ∼ 2 sub-millimeter galaxies; and αCO = 1.0 inferred by Calistro Rivera et al. (2018) from CO kinematic studies of z = 2 − 3 DSFGs (a similar value has been derived by Bothwell et al. 2013). These are consistent with the dynamical mass constraints from the CO(5-4) and (8-7) velocity maps (R15b, Swinbank et al. 2015). Fig. 9 shows the resolved t dep map in SDP.81. With the CO(3-2) emission concentrated to the north, we find a strong variation of t dep across the FIR bright source, from ∼ 1 Myr in the south to ∼ 20 Myr in the north. By adopting the Sharon et al. (2016) r 3/1 value, the relative uncertainty is ∼33 per cent; the flux calibration uncertainty on the CO(3-2)/FIR ratio is ∼ 15 per cent.
Due to the low S/N and structure in the reconstructed CO(3-2) emission, we refrain from investigating the Kennicutt-Schmidt relation, particularly the slope of ΣH 2 − ΣSFR relation. However, as the concentration of the CO(3-2) emission to the north (and the low CO(3-2) surface brightness to the south) is robustly recovered by the lens modelling (Appendix B), the general t dep trend is robust. As a further check, using the CO(5-4) line as a molecular gas tracer 4 and adopting R51 = 0.32 ± 0.05 from the Bothwell et al.  2017) found a comparably low t dep (10-100 Myr) in an extreme strongly lensed z ∼ 3 starburst PLCK G244.8+54.9, albeit at much higher star-formation surface densities (ΣSFR ∼ 2 × 10 3 M yr −1 kpc −2 ). Given the extremely short depletion time in the centre of SDP.81, sustaining the high ΣSFR requires an efficient loss of gas angular momentum to feed the star formation.

[CII]/FIR ratio and deficit
Studies of nearby star-forming galaxies (e.g., Malhotra et al. 2001;Luhman et al. 2003) have revealed that the [C ii]/FIR ratio decreases with increasing FIR luminosity (i.e., SFR)the so-called [C ii]/FIR deficit. This trend was confirmed by Herschel FIR surveys of local main-sequence galaxies (e.g., Smith et al. 2017) and ULIRGs (e.g., Díaz-Santos et al. 2013;Díaz-Santos et al. 2017), as well as by unresolved observations of high-redshift DSFGs (e.g., Stacey et al. 2010;Gullberg et al. 2015;Spilker et al. 2016). In particular, resolved observations at z ∼ 0 have found a tight correlation Based on the source-plane FIR and [C ii] luminosities (Table 3), we infer a spatially integrated [C ii]/FIR luminosity ratio of (1.5 ± 0.3) × 10 −3 . This is consistent with low [C ii]/FIR ratios seen in other high-redshift sources (Fig. 11), as well as the Herschel observations of ULIRGs from the GOALS sample (Díaz-Santos et al. 2017). However, as the [C ii] emission is significantly more extended than the FIR continuum (in line with Gullberg et al. 2018;Rybak et al. 2019), the source-averaged [C ii/FIR measurement provides only an upper limit on the [C ii]/FIR ratio in the central, FIR-bright region. The large extent of the [C ii] emission compared to the FIR continuum implies that the spatiallyaveraged [C ii]/FIR measurements in high-redshift sourcesas inferred from unresolved observations -mask the potentially very strong [C ii]/FIR deficit in the actual star-forming regions.
Considering the resolved [C ii]/FIR maps in Fig. 10, the [C ii]/FIR ratio varies by more than 1 dex across the source, down to ∼ (2 − 3) × 10 −4 in the two FIR-bright clumps -almost an order of magnitude lower than the spatiallyintegrated [C ii]/FIR ratio. We find significant (≥ 3σ) gradients in the [C ii]/FIR ratio on sub-kpc scales. Namely, on the outskirts, the [C ii]/FIR ratio is consistent with expected gas heating efficiency via the photoelectric effect (0.1-1.0 per cent, e.g., Bakes & Tielens 1994 Fig. 7, the reduced efficiency of [C ii] cooling is partially compensated by the increase in CO/[C ii] luminosity ratio, as expected from radiative transfer arguments (e.g., Narayanan & Krumholz 2017). We note that several isolated pixels at the very outskirts show very loiw [C ii/FIR ratios -these are likely lens-modelling artifacts.

Modelling setup
We adopt the line and FIR continuum predictions from the PDRToolbox models (Kaufman et al. 1999(Kaufman et al. , 2006Pound & Wolfire 2008). The PDRToolbox models are calculated for a semi-infinite slab geometry, assuming a solar metallicity and stopping the calculation at the depth corresponding to a visual extinction AV = 10. The PDRToolbox models have been widely used in interpreting atomic and molecular line observations from both z ∼ 0 ULIRGs (e.g., Díaz-Santos et al. 2017) and high-redshift DSFGs (e.g., Valtchanov et al. 2011;Gullberg et al. 2015;Wardlow et al. 2017;Rybak et al. 2019); we now apply the same approach to SDP.81. As the FIR studies of metallicity tracers in z ≥ 1 DSFGs indicate metallicities ≥ 1 Z (Wardlow et al. 2017), we consider the Z = 1 Z default PDRToolbox model as appropriate for dust-rich ISM of SDP.81.
The default PDRToolbox models are sampled onto a 0.25 dex grid which is too coarse to estimate the uncertainties on the inferred properties. We re-sample them onto a finer grid (steps of 0.05 dex), using a degree 2 spline for interpolation.
For the comparison with observations, we consider the following independent line ratios: [C ii]/FIR, [C ii]/CO(5-4), [C ii]/CO(3-2), CO(8-7)/CO(5-4) and CO(10-9)/CO(5-4). The unresolved lines -[O ii] and CO(7-6) are excluded from this comparison, as their magnification factors are unknown. Finally it should be noted that the [O i] 63-µm line can undergo significant self-absorption (e.g., Rosenberg et al. 2015). The line ratios are calculated using line luminosities in units of L . In addition to the surface-brightness uncertainty due to lens modelling, we include an additional 10 per cent flux calibration error for each tracer, as appropriate for high-S/N ALMA observations. In the line-ratio maps, we consider pixels with SNR≥3 as robust detections; for pixels with SNR< 3, we adopt 3σ upper limits. We then calculate the χ 2 for each G, n model. Following Sawicki (2012), we include the 3σ upper limits as where R data i , R model i are the measured and predicted values of the i-th line ratio, the indices j, j run over the line ratios with SNR≥ 3 (detections) and <3 (upper limits), respectively, and σj, σ j denote corresponding uncertainties. Considering only the SNR≥3 detections (i.e. dropping the second term in Equation 4) changes the best-fitting G, n only marginally (≤0.1 dex).
Several corrections need to be applied to the data before a direct comparison with the PDR models. First, the observed [C ii] luminosity has to be corrected for the contribution from non-PDR sources (i.e. ionized gas). The con-tribution of ionized gas to [C ii] emission can be estimated using [N ii] 122/205-µm lines (Croxall et al. 2017;Herrera-Camus et al. 2018b Second, the PDRToolbox predictions need to be adjusted in case a ratio of an optically thin and an optically thick tracer is considered. Namely, the emergent line intensities are predicted for a semi-infinite slab of gas illuminated only from the face side; for the extended starburst environment in SDP.81, the clouds are likely exposed to FUV radiation from multiple sides. While for the optically-thick tracers, only the emission from the face side of the cloud is observed, for the optically-thin tracers, both the sides facing to and away from the observer will contribute to the observed flux, and the predicted fluxes need to be doubled. The FIR continuum is typically optically thin down to rest-frame 70-100µm (e.g., Casey et al. 2014). We assume that the [C ii] emission is optically thin, although moderate optical depths (τ ∼ 1) have been proposed for some Galactic PDRs (e.g., Graf et al. 2012;Sandell et al. 2015) and high-redshift sources (Gullberg et al. 2015). The CO rotational lines are all optically thick, with τ ≥ 10 at ΣSFR > 1 M yr −1 kpc −2 (Narayanan & Krumholz 2014).
As the results of PDR modelling depend directly on the reconstructed source-plane surface brightness distribution for each tracer, the robustness of the inferred PDR properties against potential bias in line ratios needs to be assessed. In Appendix C, we investigate the effect of an extreme (order unity) bias in the measured surface brightness ratios for all the tracers considered; note that this is much larger than the surface brightness bias due to lens modelling for even the lowest-S/N tracer considered. We find that the inferred G, n(H) values shift by ≤ 0.5 dex even in such an extreme scenario, confirming that the PDR results are robust to (reasonable) bias in measured surface brightness ratios.
Finally, we consider the CMB temperature effects at z = 3 on the observed line luminosities. Adopting the da Cunha et al. (2013) models (for n(H)∼ 10 4 cm −3 , see results below), the CO line luminosity is biased by ≤10 per cent even for the most affected line (CO(3-2)), below the flux calibration uncertainty. We therefore consider the CMB effect on the inferred PDR models to be negligible.

Global PDR model
We first apply the PDR modelling to the spatially-integrated source-plane line and continuum luminosities. Fig. 12 shows and CO(1-0) lines and FIR continuum), our G, n estimates are ∼1 dex higher. This is mainly due to two factors; (1) the lower [C ii] line luminosity measured by ALMA and the re-reduced Herschel spectroscopy (Section 3); (2) we are considering CO Jupp ≥ 3lines, which trace denser molecular gas that CO(1-0) in Valtchanov et al. (2011) (see Section 5.3 for a detailed discussion of systematics).
5 Due to the shape of the [C ii]/FIR isocontours in the G 0 /n(H) plane, a pure χ 2 optimization might yield very low G 0 values for the pixels in the densest regions. This would imply extremely dense regions without any star-formation; we discount these solutions as unphysical. The same issue has been encountered by Díaz-Santos et al. (2017) in the analysis of nearby ULIRGs; they likewise discard the low G 0 /n solutions.
To assess the impact of differential magnification on the global PDR model, we apply the same analysis to the skyplane (i.e., not de-lensed) line and FIR luminosities (Fig. 12). The best fitting G, n sky-plane model is only marginally (∼ 0.1 dex) offset from the source-plane one. This is due to the large spatial extent (compared to the beam size) of all tracers considered; this significantly reduces the differential magnification bias.

Resolved PDR model
We now perform the resolved PDR modelling, using the line ratios and upper limits for each 200-pc pixel. For this part of the analysis, we adopt conservative uncertainties estimates. The total error budget consists of: • Flux calibration uncertainty, assumed to be 10 per cent in each Band.
• Pixel-by-pixel surface brightness uncertainty from the lens modelling for a fixed source regularization parameter λS, estimated using the scatter in source realizations in Section 3.
• Source-averaged surface-brightness uncertainty estimated from reconstructions of mock sources (Appendix B), which account for potential bias in the maximum a posteriori λS. These vary between 6 per cent (FIR continuum) to Figure 12. Global PDR models for SDP.81, using spatially-integrated source-plane luminosities in the sky plane (right) and source plane (left), with 1σ confidence regions denoted by the colour bands. The differential magnification effect is relatively limited due to the extended nature of all tracers; the best-fitting models differ by ∼0.1 dex in G and n. 30 per cent for CO(10-9) and 65 per cent for CO(3-2) (Table B1). For CO(3-2) and CO(10-9), this error dominates the error budget.
During the fitting procedure, we discard the solutions with G0/n(H)≤ 10 −3 , which would imply an unphysically large source (c.f. Díaz-Santos et al. 2017;Wardlow et al. 2017). Fig. 13 shows the G0 and n(H) values inferred for individual 200-pc regions in SDP.81, given the surface brightness distributions from Fig. 2. In the central, FIR-bright region, we find strong FUV fields with G = 10 3.5 −10 4.0 G0, and high density n(H)=10 5.0 − 10 5.5 cm −3 , with a typical uncertainty of ∼ 0.2 dex. These agree well with the global PDR model. We find a significant (> 3σ) gradient in G in the East-West direction: G increases from 10 3.0±0.2 to 10 4.0±0.1 G0. While the star-burst region might exhibit increased gas density, we do not find a significant spatial variation in n(H). The most direct interpretation is that the star-formation is constrained to the central, disc-like part of the gas reservoir.
The PDR surface temperature TPDR varies from about 150 K in the outskirts to 1500 K in the FIR-bright clumps. However, with TPDR S/N ≤5 across the source, we can not robustly measure the spatial variation of TPDR. The full implications of the high TPDR and its local variations for the [C ii]/FIR deficit are discussed in Section 4.3.
How do these values compare to measurements for other star-forming systems? Fig. 14 compares the resolved G, n value in SDP.81 to unresolved PDR models of z = 2 − 5 DSFGs from Gullberg et al. (2015) and stacked DSFGs from Wardlow et al. (2017). Compared to both samples, SDP.81 shows a similar FUV field strength, with density ≥ 1 dex above the Wardlow et al. (2017) range, and comparable to the densest Gullberg et al. (2015) sources. Using resolved [C ii] and CO(3-2) and FIR ALMA imaging, Rybak et al. (2019) inferred G ∼ 10 4 G0, n ∼ 10 4 − 10 5 cm −3 in the central regions of two z = 2.94 DSFGs from the ALESS sample; SDP.81 shows slightly higher G and lower n(H). show densities comparable to large molecular clouds such as Sgr B2 and 30 Dor, albeit with FUV fields ∼ 1 dex stronger; these might indicate a closer distance between the stars providing FUV radiation and the PDR surfaces, compared to 11-80 pc in 30 Dor (Chevance et al. 2016). Indeed, close associations of the O/B stars and gas clouds have been invoked by Herrera-Camus et al. (2018a) to explain the [C ii]/FIR deficit in z ∼ 0 star-forming galaxies and ULIRGs.
On sub-pc scales, we compare the inferred PDR properties with the resolved (∼ 0.05 pc) study of Orion by Goicoechea et al. (2015), who used [C ii], FIR continuum and CO(2-1) observations to infer the PDR conditions in the OMC 1. Using the Stacey et al. (2010) diagnostic diagram (based on Kaufman et al. 1999Kaufman et al. , 2006 PDR models used in our analysis), Goicoechea et al. (2015) obtain G 3 × 10 4 G0, n ≥ 10 5 cm −3 for the region closest (R ≤0.1 kpc) to the Trapezium cluster, decreasing down to G 10 3 G0, n 10 4 cm −3 . These are, respectively, on the upper and lower end of the conditions seen in the FIRbright region of SDP.81; however, in SDP.81, these describe physical conditions averaged over 200-pc scales (which likely include a number of star-forming regions and the voids between them), underlining the extreme nature of the starforming regions in DSFGs. Figure 13. Left column: FUV field strength G 0 , density n(H) and PDR surface temperature T PDR inferred from a comparison with the PDRToolbox models (Kaufman et al. 1999(Kaufman et al. , 2006Pound & Wolfire 2008), with the corresponding 1σ uncertainties (middle column), reduced χ 2 per pixel and the number of line ratios used to constrain the fit for each pixel (right column). The pixel size is 200×200 pc. In the FIR-bright region (as denoted by black contours), G 10 3.5 − 10 4.0 G 0 , n(H) 10 4.5 − 10 5.0 cm −3 . The FIR-bright regions of the source shows T PDR much higher than the [C ii] transition temperature (92 K). The pixels for which the inferred G 0 and n have S/N<3 and/or with reduced χ 2 ≥ 5 have been masked for clarity.

Caveats and limitations
Our PDR analysis -both resolved and unresolved -is based on semi-infinite slab, uniform-density PDR models, which are compared with the FIR continuum (160-330 µm restframe), [C ii] and CO (Jupp ≥ 3) observations. Our choice of the modelling suite -the PDRToolbox models -was motivated by the satisfactory performance of the models given the lines studied here and corresponding the uncertainties, and for consistency with earlier, source-averaged studies of ISM properties in DSFGs. The PDRToolbox models reproduce the observed line ratios well (reduced χ 2 ≤ 5), at least given the spatial resolution, diagnostics considered and S/N of the data at hand. Here, we summarize the limitations of this approach: • Large surface-brightness uncertainties for CO(3-2) and CO(10-9). As outlined at the beginning of this Section, we adopt conservative surface-brightness uncertainty estimates to avoid over-interpreting the data. For the CO(3-2) and CO(10-9) lines, these are dominated by the error on source reconstruction due to their low S/N. Namely, excluding the uncertainties from Appendix B from the error budget decreases the quality of the best-fitting models in the FIRbright clumps(reduced χ 2 = 10 − 25). This suggests that an additional heating mechanism or a multi-phase ISM model would be more appropriate in high-ΣSFR clumps. However, as shown in Appendix B, the S/N of the CO(3-2) and (10-9) data at hand is currently insufficient to fully constrain the relative importance of different heating mechanisms.
• Geometry effects. Our semi-infinite, uniform-density slab models (even after adjusting for multi-side illumination) are likely too simplistic given the complex, clumpy Figure 14. Resolved G, n distribution in SDP.81 (grey, contours drawn at the 25th, 50th and 75th percentile), compared to the global model (Section 5.2), unresolved values inferred for DSFGs samples of Brisbin et al. (2015), Gullberg et al. (2015) and Wardlow et al. (2017), and central regions of two ALESS z ∼ 3 DSFGs (Rybak et al. 2019), z ∼ 0 ULIRGs (Díaz-Santos et al. 2017), nearby star-forming regions (adopted from Oberst et al. 2011) and resolved measurements in the vicinity of the Trapezium cluster in Orion (Goicoechea et al. 2015). Black errorbars indicate the median 1σ uncertainties for the SDP.81 data. The diagonal dashed line indicates the boundary of the region where the dust drift velocity exceeds the average turbulent velocity assumed by the PDRToolbox models (Kaufman et al. 1999); above these line, the PDR model becomes internally inconsistent.
star-forming environment expected in SDP.81, where strong gradients in ISM density and composition along the line-ofsight might be expected. These can be properly accounted for by 3D radiative transfer codes that account for the full spatial stratification of the (multi-phase) ISM, and a realistic distribution of radiation sources.
• Isochoric models. The self-gravity and external pressure on the molecular clouds will result in a density stratification, which is not included in our models. The effects of assumed cloud density profiles were recently studied by Popping et al. (2019), who found that different assumed density profiles result in [C ii] and CO (Jupp ≤ 5) luminosities changing by up to ∼0.5 dex. This effect might be even more dramatic for the CO(8-7) and CO(10-9) lines. Alternatively, isobaric models -e.g., as in a recent application of the (Le Petit et al. 2006) code to the resolved observations of 30 Doradus by Chevance et al. 2016 -might provide a better description of molecular clouds in intensely star-forming environment. However, Chevance et al. (2016) do not find significant differences between gas properies inferred using the isochoric and isobaric assumptions.
• Reliance on mid-to high-J CO lines. Our PDR modelling is largely driven by the CO(5-4) and CO(8-7) lines, however, the predictions for the emergent line intensities can vary significantly between different PDR codes. This is because the high-J CO emission is confined to the outer cloud layer; therefore, differences in CO abundance will result in a different predictions (e.g., Röllig et al. 2007). Therefore, the inferred gas properties might depend strongly on the underlying PDR model.
• Limited sensitivity to cold, low-density gas. By including only Jupp ≥ 3 CO lines (critical density ncrit ≥ 8 × 10 3 cm −3 , we effectively discount the extended cold, low-density (n ≤ 10 3 cm −3 ) gas phase traced by CO(1-0) and CO(2-1) lines. Indeed, studies of the CO ladders including the Jupp = 1 − 2 lines have found significant cold gas components both in z ∼ 0 ULIRGs (e.g., Kamenetzky et al. 2014;Greve et al. 2014) and high-redshift DSFGs (e.g., Bothwell et al. 2013;Yang et al. 2017). While the cold gas component can contribute to the CO(3-2) and CO(5-4) luminosity, the Yang et al. (2017) source-averaged CO SLED model for SDP.81 -which includes both a cold and a warm component -shows the cold-gas contribution to the CO(3-2) and CO(5-4) luminosities to be ≤30 and <10 per cent, respectively. Therefore, our results should not be significantly affected by the cold gas component (Appendix C).
• Mechanical and cosmic/X-ray heating. In the high-ΣSFR regime, additional mechanism not included in the PDRToolbox models -e.g., mechanical heating by shocks, X-ray radiation and cosmic rays -can contribute significantly to the total gas heating. Although the line ratios in SDP.81 are well-reproduce by models that include only far-UV heating by stars; some of the key diagnostics of e.g., mechanical heating, are not available for SDP.81 (Section 4.2.1). Future observations of additional tracers (e.g., high-J CO lines, mid-J HCN/HCO + lines, CO isotopologues) will be necessary to properly assess the relative importance of different heating mechanisms.

What drives the [CII]/FIR deficit in SDP.81?
With the inferred G, n and TPDR maps in hand, we now address the origin of the [C ii]/FIR deficit in SDP.81.
Recently, Rybak et al. (2019) used resolved (0.15-0.5 arcsec) observations of [C ii], CO(3-2) and FIR continuum emission in two unlensed z ∼ 3 DSFGs to show that the [C ii]/FIR deficit in their central regions is due to thermal saturation. We apply the analysis of Rybak et al. (2019) to the 200-pc resolution maps of SDP.81. Namely, we consider the following mechanisms: (1) thermal saturation of the [C ii] fine-structure levels; (2) reduced gas heating due to positive grain charging; (3) AGN contribution.

Reduced photoelectric heating
We calculate the photoelectric heating rate per hydrogen atom following Wolfire et al. (2003): where G = 1.7 × G0, Z d is the dust-to-gas ratio (normalized to the Galactic value), ne the electron density and φPAH a factor associated with the PAH molecules. Assuming ne = 1.1 × 10 −4 n and φPAH = 0.5, we calculate ΓPE for each resolution 200×200 pc resolution element. The ΓPE varies by only ∼50 per cent across the source; this variation is insufficient to explain the ∼1.5 dex variation in the observed [C ii]/FIR ratio.

AGN effects
At sub-kpc scales, a strong [C ii]/FIR deficit can be driven by AGNs, by simultaneously suppressing the [C ii] emission by X-ray-driven ionization of C + to C 2+ etc., and increasing the FIR luminosity due to extra FUV flux from the accretion disk. Resolved [C ii]/FIR studies of local active galaxies (Smith et al. 2017;Herrera-Camus et al. 2018a) have revealed a strong [C ii]/FIR drop within ∼ 500-pc radius of the AGN in some galaxies. While the AGN influence will be largely negligible for kpc-scale deficits (c.f., Rybak et al. 2019), the AGN sphere of influence is comparable to the 200-pc resolution of our data. XMM observations of the SDP.81 field detect 0.5-8 keV flux of (49 ± 7) × 10 −15 erg s −1 cm −2 , with a hardness ratio of -0.4 at the position of SDP.81 (Ranalli et al. 2015). However, it is unclear whether this flux is associated with the lensed DSFG, or the AGN in the z = 0.299 lens, which is detected in radio and mm-wave imaging. Assuming the X-ray flux is all due to a hypothetical AGN at the position of the northern FIR-bright clump (lensing magnification µ 20), we find a source-plane Xray luminosity L 0.5−8.0 keV 2 × 10 44 erg s −1 (uncorrected for obscuration). This is similar to X-ray luminosities observed in DSFGs from the ALESS sample (L 0.5−8.0 keV = 23 × 10 42 − 3 × 10 44 erg s −1 , Wang et al. 2013) 7 . Even assuming the highest column density for the ALESS sources from Wang et al. (2013), (NH 5 × 10 24 cm −3 ), a comparison with carbon ionization model in the vicinity of Xray bright AGNs of Langer & Pineda (2015) shows that the X-ray flux in SDP.81 is insufficient to drive the [C ii]/FIR deficit on scales beyond few hundred pc (c.f., Rybak et al. 2019). Finally, circumstantial evidence also disfavours the AGN-driven [C ii]/FIR deficit: (1) a single AGN is hard to reconcile with the two similarly deep [C ii]/FIR minima separated by ∼1 kpc (Fig. 11; although an "Arp 220"-like merging configuration would account for this) and (2) although the CO SLED is more excited in the northern part of the FIR-bright region (as evidenced by the higher CO(10-9) surface brightness, Fig. 2), the CO excitation is fully consistent with star-formation, without a need for additional AGNassociated X-ray/cosmic-ray heating.

C + thermal saturation
At high gas temperature Tgas, the relative occupancy of the upper/lower C + fine-structure levels becomes thermally saturated (Muñoz & Oh 2016) at temperatures higher than the Figure 15. Comparison of the [C ii]/FIR ratio in SDP.81, corrected for f [CII] , the fraction of gas traced by the [C ii] (upper), and without the correction (lower), with the thermal saturation model of Muñoz & Oh (2016) (dashed line). The solid line and the shaded region indicate the best-fitting power-law function to the data with the associated 1σ confidence region. The f [CII] correction significantly reduces the tension between the observed [C ii]/FIR ratio and the thermal saturation model.

CONCLUSIONS
We have presented matched-resolution source-plane reconstructions of the [C ii], CO (3-2), CO(10-9) and the ALMA Bands 6, 7 and 8 (rest-frame 160-320 µm) continuum emission in a z = 3.042 gravitationally lensed, dusty star-forming galaxy SDP.81. Using a visibility-fitting lens modelling code with an adaptive source-plane grid (Rybak et al. 2015a), we achieve a median source-plane resolution of ∼200 pc. Combining these with the CO(5-4) and CO(8-7) data from Rybak et al. (2015b), we use PDR models of Kaufman et al. (1999Kaufman et al. ( , 2006; Pound & Wolfire (2008), to map the physical conditions of the star-forming ISM -in particular, the FUV field strength, gas density and PDR surface temperature -on sub-kpc scales. Compared to other studies of the [C ii]/FIR deficit at high redshift, we are able to leverage the resolved CO and FIR continuum observations to characterize the ISM properties in SDP.81 on sub-kpc scales. Our main results are: • The ALMA Band 8 continuum source-plane surface brightness distribution is consistent with the previouslypublished Band 6 and Band 7 continuum reconstructions (R15a, Dye et al. 2015;Hezaveh et al. 2016). Fitting for the dust temperature on 200-pc scales, we do not find evidence for a strong dust temperature gradient, and T dust ≥ 80 K, in tension with models proposed by Calistro Rivera et al. (2018).
• The [C ii] emission is significantly more extended than the FIR continuum and CO(Jupp ≥ 5) emission: ∼50 per cent of [C ii] emission originates outside the FIR-bright region. The [C ii] maps reveal an extended, low surfacebrightness emission extending to ∼4 kpc from SDP.81, providing further evidence to SDP.81 being a merger-induced starburst.
• The CO(3-2) emission is concentrated to the north of the source, with the CO(3-2)/FIR ratio varying rapidly across the source (by ≥ 1 dex). Using the CO(3-2) as a molecular gas tracer, we find an extremely short gas depletion timescale of t dep = 10 6 − 2 × 10 7 Myr, significantly shorter than in most other resolved t dep studies in highredshift DSFGs.
• The CO SLED in SDP.81 is comparable to that of Class I ULIRGs from Rosenberg et al. (2015). Given the data at hand, the CO SLED is consistent with PDR models, without a need for additional energy source (e.g., mechanical heating). Although globally, the CO lines contribute only marginally (few per cent) to the total gas cooling, the CO/[C ii] ratio can be as high as 25 per cent for the FIRbright clumps.
• Comparing the source-integrated (unresolved) delensed [C ii], CO and FIR emission with PDR models of Kaufman et al. (1999Kaufman et al. ( , 2006; Pound & Wolfire (2008), we find FUV field strength G = 10 3 G0 and gas density n(H)=10 5 cm −3 , significantly higher than the earlier model by Valtchanov et al. (2011). The source-averaged PDR properties inferred from the de-lensed (source-plane) and observed (sky-plane) luminosities are consistent with each other; the differential lensing does not bias the global PDR model significantly.
• Modelling the PDR conditions on 200-pc scales, we find the FUV field strength to vary considerably across the source (G = 10 2.5 − 10 4.0 G0). In contrast, n(H) varies only marginally (10 5.0 −10 5.5 cm −3 ). Compared to the unresolved models, the resolved PDR model shows G up to 1 dex higher; this is due to the strong variation of [C ii]/FIR ratio (which determines G) on sub-kpc scales, the FIR continuum being much more compact than [C ii]. The physical conditions averaged over 200-pc scales are comparable to the vicinity of the Trapezium cluster in Orion (Goicoechea et al. 2015), underlying the extreme nature of DSFGs, and highlighting a need for more sophisticated radiative transfer models for interpreting the observed line ratios.
• The [C ii]/FIR ratio in SDP.81 varies by ∼1.5 dex across the source, with a strong [C ii]/FIR deficit in the central FIR-bright clumps (L [CII] /LFIR = 3 × 10 −4 ). Although the [C ii] surface brightness generally increases towards the FIRbright part of the source, it drops at the position of the FIR-bright clumps.
• Given the resolved PDR properties, the most likely cause of the [C ii]/FIR deficit in SDP.81 is thermal saturation of the C + fine-level structure due to strong FUV fields and high PDR surface temperature. After correcting the resolved [C ii]/FIR for the fraction of gas traced by the [C ii] emission, we find the [C ii]/FIR trend to be consistent with the thermal saturation model of Muñoz & Oh (2016).
We emphasise that the agreement between the observed line ratios and the simplistic PDR models is largely due to significant systematic uncertainties on the CO(3-2) and CO(10-9) source reconstructions. High-sensitivity observations of these two lines, and additional fainter diagnostics (e.g., HCN and HCO + lines, CO isotopologues) are necessary to discriminate between e.g., different heating mechanism on sub-kpc scales.
This paper presents a first look into the PDR properties of high-redshift, dusty star-forming galaxies on scales of few hundred pc, and a major improvement on the unresolved studies of DSFGs. Our results reveal the complex nature of DSFGs, with strong spatial variations in [C ii]/FIR ratio, the FUV field strength and the PDR surface temperature. By studying the properties of the photon-dominated regions on 200 pc scales, this work bridges the high-redshift starbursts and local studies of large star-forming regions. We make the reconstructed source-plane maps of the CO ladder and [C ii] line and the ALMA Bands 6, 7 and 8 continuum (Figures 2,  A2) publicly available at sdp81.strw.leidenuniv.nl. We hope this legacy dataset will facilitate independent analysis of the conditions in high-redshift DSFGs on sub-kpc scales. First, we estimate the systematic error on the reconstructed source-plane surface brightness distribution, which we use in Section 5.3. Namely, we compare the input surface brightness for the uniform-brightness model (Fig. B1) to the mean inferred surface brightness, integrated over the extent of the input source. The input and inferred surface brightness values, and the corresponding uncertainites are listed in Table B1. For the continuum data, [C ii] and CO(5-4) and (8-7) line, the resulting uncertainty is ∼10 per cent, i.e. comparable to the flux calibration uncertainty. However, for the CO(3-2) and (10-9) lines, which are observed at lower S/N, the fractional uncertainty is 65 and 30 per cent, respectively.
Second, we consider the reliability of the inferred source morphology: • FIR continuum: the reconstructed Bands 6, 7 and 8 data recover the input source structure reliably, including the positions and separation of FIR-bright clumps in Fig. B2.
• CO(3-2): the low S/N and curvature regularization cause the source to be overly extended in the east-west (i.e., radial) direction; the magnification in this direction is ∼ 1. The mean surface brightness is recovered within 20 per cent. However, neither of our mock reconstructions shows the increase in the CO(3-2) luminosity towards the north of the source, seen in the real data (Fig. 2). We therefore conclude that the observed variations in the CO(3-2) luminosity (and ΣH 2 in the depletion time analysis) is physical.
• CO(5-4) and (8-7): both the uniform and FIR-like source morphologies are reliably recovered. The surface Figure A1. Lens modelling results: maximum a posteriori models for (top to bottom): the Bands 6, 7 and 8 continuum, the CO(3-2) and the CO(10-9) lines (lower ). Individual columns are ordered as follows: (1) dirty image of the data; (2) dirty image of the maximum a posteriori (MAP) model (normalized to the peak surface brightness of the data dirty image); (3) dirty-image residuals, with contours starting at ±2σrms and increasing by steps of 1σrms; (4) MAP sky model (normalized to the peak surface brightness); (5) MAP reconstructed source, with contours drawn at 25, 50 and 75 per cent of the surface brightness maximum. Source-plane caustics are marked by the grey line. brightness bias is somewhat larger for the CO(8-7) line due to its lower S/N (15 per cent, compared to 10 per cent for CO(5-4). This confirms the robustness of the source-plane structures from Rybak et al. (2015b), namely the extent of the CO(5-4)/(8-7) emission, and their offsets with respect to each other, and the FIR continuum (Fig. 2).
• CO(10-9): the reconstructed source shows spurious structure for both mock source realizations. However, the observed concentration of the CO(10-9) emission to the north of the source is not (spuriously) reproduced in either of the mock source reconstructions; we thus consider it to be a physical. On the other hand, the low-surface brightness CO(10-9) to the south-west in Fig. 2 is reminiscent of the spurious structures seen in the mock reconstructions. As this faint component contributes negligibly to the CO SLED and PDR modelling results in this paper, we do not exclude it from our analysis. However, we note that given the S/N of the data at hand, this could be a modelling artifact. C ii : both the uniform and FIR-like source are reliably recovered. In particular, for the FIR-based source, the positions of the FIR surface brightness maxima are recovered in the reconstructed maps. This confirms that the decline in [C ii] surface brightness in these regions seen in the real data is physical, rather than a modelling artifact.     Table B1. Systematic error on the mean source surface brightness δ estimated from the mock source reconstructions in Fig. B1. The columns list the input (uniform) surface brightness (S in ), the measured mean surface brightness with 1σ uncertainties ( Sout ), and the fractional systematic error δ.  Figure C1. Sensitivity of G and n(H) values to inferred line ratios. The isocontours correspond to typical line ratios measured in SDP.81 (solid lines), and increased by a factor of 2 (dashed lines). The shaded regions correspond to an assumed 30 per cent uncertainty on line ratios.