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

ABSTRACT

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, multiband ALMA observations of the far-infrared (FIR) continuum, CO ladder, and the [C ii] 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-ultraviolet (FUV) field strength, density, and PDR surface temperature – of the star-forming gas on 200-pc scales, finding a FUV field strength of ∼10³−10⁴G₀, gas density of ∼10⁵ cm⁻³, and cloud surface temperatures up to 1500 K, similar to those in the Orion Trapezium region. The [C ii] emission is significantly more extended than that FIR continuum: ∼50 per cent of [C ii] emission arises outside the FIR-bright region. The resolved [C ii]/FIR ratio varies by almost 2 dex across the source, down to ∼2 × 10⁻⁴ in the star-forming clumps. The observed [C ii]/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.

1 INTRODUCTION

Dusty star-forming galaxies (DSFGs) with star formation rates (SFRs) of 100–1000 M_⊙ yr⁻¹ 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-ultraviolet (FUV) radiation from the young, massive stars is often completely obscured by their large dust reservoirs and re-radiated at rest-frame 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 FUV radiation from the new-born massive O/B-type stars. This results in a thermal and chemical stratification of the GMC: while deep in the cloud, the gas is in the molecular form (H₂, 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 (PDR) 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.

The PDR properties in present-day (ultra)luminous infrared galaxies [(U)LIRGs] have been extensively studied using FIR and mm-wave spectroscopy. An early study by Stacey et al. (1991) – combining Kuiper Airborne Observatory [C ii] and FIR continuum observations with ground-based CO(1–0) observations – found a range of densities n(H)=10³−10⁶ cm⁻³, with G in excess of 10³ G₀.¹ These studies were revolutionized by Herschel FIR spectroscopy, opening doors to large systematic surveys of FIR emission lines down to ∼500-pc scales in nearby star-forming galaxies and (U)LIRGs (e.g. Graciá-Carpio et al. 2011; Díaz-Santos et al. 2013, 2017; Rosenberg et al. 2015; Smith et al. 2017; Herrera-Camus et al. 2018b).

Previous studies of gas and dust properties in DSFGs have been limited to spatially-integrated quantities, due to their faintness (S_850 μ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 12 |$z$| = 1−2 galaxies to infer typical G = 10^2.5−10⁴ G₀, n(H)≃ 10³−10^4.5 cm⁻³. 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. 2010, 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-|$\mu$|m and ground-based [C ii] and low-J CO observations to infer n(H)=10³−10² cm⁻³, G = 10¹−10³ G₀ 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¹−10³ cm⁻³, G = 10^2.2−10^4.5 G₀. Similarly, Zhang et al. (2018) used Herschel [N ii] and [O i] observations for the H-ATLAS DSFGs to infer typical densities of n(H)=10³−10⁴ cm⁻³.

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 co-spatial. 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⁴ G₀ and moderately high densities n(H) = 10^3.5−10^4.5 cm⁻³. 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, multiband 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 the 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 DSFGs 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 analysed 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 Lambda cold dark matter cosmology from Planck Collaboration XIII (2016). Adopting the spectroscopic redshift of |$z$| = 3.042 (ALMA Partnership et al. 2015), this translates into a luminosity distance of 26480 Mpc, and an angular scale of 1 arcsec = 7.857 kpc (Wright 2006).

2 OBSERVATIONS

2.1 Target description

SDP.81 (J2000 09:03:11.6 + 00:39:06; Bussmann et al. 2013) is a |$z$|_S = 3.042 DSFG, strongly lensed by a foreground early-type galaxy with |$z$|_L = 0.299. Identified in the H-ATLAS survey (Eales et al. 2010; Negrello et al. 2010, 2017), SDP.81 was targeted by the ALMA 2014 Long Baseline Campaign (LBC; ALMA Partnership et al. 2015). Ancillary observations include Hubble Space Telescope (HST) near-infrared imaging (Dye et al. 2014), HerschelFIR spectroscopy (Valtchanov et al. 2011; Zhang et al. 2018), and radio/mm-wave low-resolution spectroscopy and imaging (Valtchanov et al. 2011; Lupu et al. 2012). These extensive high-resolution, high signal-to-noise ratio (S/N) ancillary data make SDP.81 perhaps the best-studied DSFG.

Rybak et al. (2015a,b, henceforth R15a and R15b, respectively) presented the reconstructed maps of Bands 6 and 7 dust continuum, and CO(5–4) and (8–7) emission with a typical resolution of 50–100 pc, alongside the reconstruction of HST near-IR imaging. These revealed a significantly disturbed morphology of stars, dust, and gas in SDP.81, with an ordered rotation in a central dusty disc some 2–3 kpc; a similar conclusion has been reached by independent analysis both in the uv plane (Hezaveh et al. 2016) and Cleaned images (Dye et al. 2015; Swinbank et al. 2015). Based on the kinematic analysis of the CO(5–4) and (8–7) maps, SDP.81 is classified as a post-coalescence merger (Dye et al. 2015; R15b; Swinbank et al. 2015). In this paper, we adopt the CO(5–4) and CO(8–7) from R15b, and obtain matched-resolution reconstructions of Band 6 and Band 7 continuum (Section 3.1.3).

2.2 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 h on-source with a resolution of 0.2 arcsec, the total observations amount to only ∼45 min. The array configuration consisted of 40 12-m 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), centred 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.

2.3 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 total on-source time was 35.8 min.The precipitable water vapour ranged between 0.52 and 0.64 mm. The frequency setup consisted of four SPWs, centred 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.

2.4 ALMA 2014 LBC observations

As the ALMA 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 h 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 H₂O (2₀₂ − 1₁₁) 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 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 data sets.

3 DATA PROCESSING AND ANALYSIS

3.1 Data preparation and imaging

All the data processing was done with casa versions 4.7 and 5.0 (McMullin et al. 2007). Fig. 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.

3.1.1 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 (Fig. 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⁻¹ (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 Tamura et al.'s (2015) synchrotron spectrum for the foreground AGN within 2σ.

3.1.2 ALMA Band 8

After applying the standard ALMA Cycle 4 calibration, we performed manual flagging in the time domain to remove several integrations that showed elevated noise levels; a total of 6 min of data was flagged. We also flagged channels affected by atmospheric lines at 468.0 and 457.3 GHz.

The [C ii] line is detected between 469.791 and 470.253 GHz. The observed image-plane [C ii] integrated line flux of 170 ± 20 Jy km s¹ ((2.7 ± 0.3) × 10⁻¹⁸ W m⁻²) is consistent with the Herschel Spire FTS measurement of (2.9 ± 0.6) × 10⁻¹⁸ W m⁻² (Zhang et al. 2018) within 1σ. This confirms that ALMA Band 8 observations are not resolving out extended [C ii] emission.²

The [C ii] absorption from diffuse ionized gas was detected towards some Galactic star-forming regions (e.g. Gerin et al. 2015). At a high redshift, Nesvadba et al. (2016) detected a [C ii] absorption towards to |$z$| = 3.4 DSFG, which they interpreted as evidence for the infalling gas. Looking for a potential [C ii] absorption in the Cleaned cube with a velocity resolution of 40 km s⁻¹, we do not find any evidence for [C ii] absorption against the 160-|$\mu$|m continuum.

3.1.3 ALMA LBC: Bands 6 and 7 continuum, CO(10–9) line

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 observations, 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 and 285.52 GHz; due to the low S/N ratio (peak image-plane S/N≃ 5, Fig. 1), we frequency-averaged 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σ.

3.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 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 centre 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 spatially integrated magnifications for each tracer as follows: We take all the source-plane pixels with S/N ≥1, 2, and 3, forward-lens 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.

3.3 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 source-plane reconstruction by performing our lens modelling analysis on mock data sets 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 data sets [∼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-|$\mu$|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 pc; 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 on to a Cartesian grid with 200 × 200 pc pixels.

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

3.3.2 [C ii] 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 frequency-averaged the line into channels 62.50 MHz (∼40 km s⁻¹) 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 S/N ≤ 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.

The reconstructed [C ii] emission map reveals a low-surface-brightness emission (⁠|$\Sigma _\mathrm{[C\, \rm {\small{II}}]}\sim 3\times 10^7\,{\rm L}_\odot$| kpc⁻²) between −20 and +170 km s⁻¹, stretching out to ∼10 kpc north of SDP.81 (Fig. 3). 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_\mathrm{[C\, \rm {\small{II}}]}$| luminosity of this component is ∼0.6 × 10⁹ L⊙, less than 10 per cent of the total source-plane [C ii] luminosity. We do not detect any FIR continuum or CO emission from the region outside the dashed box in Fig. 3, 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 sub-structure detected in the HST WFC3 160W imaging (Dye et al. 2014; R15b). Given the perturbed velocity/velocity-dispersion fields (R15b), we interpret this morphology as a potential low-mass companion; assuming the M_gas scales with |$L_\mathrm{[C\, \rm {\small{II}}]}$|⁠, the gas mass ratio of the two components is 10:1. Assuming this [C ii] emission is due to star formation, we use Herrera-Camus et al.'s (2015) SFR-[C ii] relation to derive a typical Σ_SFR = 1 M_⊙ yr⁻¹. It total, ∼50 per cent of the total [C ii] luminosity arises from regions with Σ_SFR ≤ 10 M_⊙ yr⁻¹.

3.3.3 CO(3–2) emission

The CO(3–2) emission is detected only 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 artefact 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.

3.3.4 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⁻¹ kpc⁻²) challenging: as we show in Appendix  C, this structure can be an artefact of the lensing reconstruction.

Finally, Fig. 4 shows the relative distribution of the FIR continuum, [C ii] and CO emission. This highlights the compact extent of the FIR continuum with respect to the CO and [C ii] lines, and the extended, low-surface brightness [C ii] emission.

4 RESULTS AND DISCUSSION

4.1 Dust continuum

4.1.1 Global dust SED

Although the increasing cosmic microwave background (CMB) temperature at high redshifts can significantly bias the inferred T_dust and L_FIR (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 L_FIR is biased by ≤0.1 dex.

which yields SFR = 410 ± 90 M_⊙ yr⁻¹. 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⁻¹ kpc⁻² as the FIR-bright region.

4.1.2 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 |$\mu$|m). These are well within the Rayleigh–Jeans tail of the modified blackbody 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, after fitting the modified blackbody 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 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 DSFGs; 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 |$\mu$|m continuum. First, at 200-pc scales probed here, the high-T_dust regions might not dominate the continuum SED, resulting in a bias towards low T_dust temperature estimates due to the mass-weighting in equation (2). Secondly, 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.

4.2 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] 63/145 |$\mu$|m fine-structure lines, or even CO rotational lines (Kaufman et al. 1999; Kaufman, Wolfire & Hollenbach 2006; Narayanan & Krumholz 2017).

Table 3 lists the inferred source-plane CO and [C ii] line luminosities, integrated over the entire source. We complemented the resolved ALMA observations by archival observations of the CO(1–0) (VLA; Valtchanov et al. 2011) and CO(7–6) lines (CARMA; Valtchanov et al. 2011), as well as an upper limit³ on the CO(9–8) lines (Z-Spec on the Caltech Submillimeter Observatory; Lupu et al. 2012). The unresolved/low-resolution nature of the archival data prevents a source-plane reconstruction. Therefore, we assume that the CO(1–0), (7–6), and (9–8) lines follow the same spatial distribution as the CO(3–2), (8–7), and (10–9) lines, respectively.

4.2.1 [C ii] dominates the gas cooling budget

We first compare the source-averaged [C ii], [O i], and CO SLED cooling rates. The [O i] 63-|$\mu$|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⁻¹⁸ W m⁻², while the 3σ upper limit on the [O i] 145-|$\mu$|m line flux is 2.3 × 10⁻¹⁸ W m⁻² (Zhang et al. 2018). This translates into image-plane line luminosities of μL_[OI]63 = (45 ± 12) × 10⁹ L⊙ and μL_[OI]145 < 5 × 10¹⁰ L_⊙. Assuming both lines follow the FIR continuum surface brightness distribution, we obtain the de-lensed source-plane luminosities of L_[OI]63 = (2.5 ± 0.7) × 10⁹ L⊙ and L_[OI]145 < 3 × 10⁹ L_⊙, respectively. The inferred [O i]/[C ii] luminosity ratio is 0.6 ± 0.2 for the [O i] 63-|$\mu$|m and ≤0.7 (but likely much lower) for the [O i] 145-|$\mu$|m line.

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 J_upp = 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 J_upp > 10 CO rotational lines are expected to contribute only marginally to the gas cooling. The total CO luminosity is |$L_\mathrm{CO}=\Sigma _{J=1}^{10} L_\mathrm{CO(J_\mathrm{upp})}\simeq 0.8\times 10^9\,{\rm L}_\odot$|⁠, ∼20 per cent of the total |$L_\mathrm{[C\, \rm {\small{II}}]}$|⁠, and ∼10 per cent of the total gas cooling (⁠|$L_\mathrm{[C\, \rm {\small{II}}]}+L_\mathrm{CO}+L_\mathrm{[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³ M_⊙ pc⁻² (cf. Narayanan & Krumholz 2017).

4.2.2 CO spectral energy distribution

The global source-plane CO SLED in SDP.81 peaks between J_upp = 5 and 9 (in units of L_⊙), with a sharp drop at J_upp = 10 (Fig. 8, see also Yang et al. 2017, for a recent sky-plane analysis). Compared to Rosenberg et al.'s (2015) study of CO excitation in present-day ULIRGs, the sharp drop in L_CO for J_upp ≥ 9 and minor CO contribution to the global neutral gas cooling makes SDP.81 akin to Rosenberg et al.'s (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. 2012, 2015), pumping the J_upp ≥ 10 CO transitions, as well as high-J HCN, HNC, and ¹³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 J_upp > 10, HCN or ¹³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 (χ² ≤ 1) by the PDRToolbox models, which do not include mechanical heating, and as the SDP.81-like Rosenberg et al.'s (2015) Class I ULIRGs are fully compatible with UV-only 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 S/N ≥ 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_\mathrm{[C\, \rm {\small{II}}]/FIR}=0.712\pm 0.006$|⁠.

4.2.3 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_\mathrm{dep}\sim \Sigma _\mathrm{H_2}/\Sigma _\mathrm{SFR}$|⁠. We estimate the molecular gas surface density |$\Sigma _\mathrm{H_2}$| using the CO(3–2) line map. For J_upp > 1 CO transitions, the |$L^{\prime }_\mathrm{CO}$| is first down-converted into |$L^{\prime }_\mathrm{CO(1-0)}$| via the |$R_{J1} = L^{\prime }_\mathrm{CO(J-J-1)}/L^{\prime }_\mathrm{CO(1-0)}$| factor, which is then converted into M_gas via |$M_\mathrm{gas}=\alpha _\mathrm{CO} L^{\prime }_\mathrm{CO(1-0)}$|⁠. We adopt R₃₁ from Sharon et al. (2016), who found a mean R₃₁ = 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 Sharon et al.'s (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 the |$\Sigma _\mathrm{H_2}$|–|$\Sigma _\mathrm{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⁴ and adopting R₅₁ = 0.32 ± 0.05 from the Bothwell et al. (2013) CO SLED survey of a sample of 40 DSFGs at |$z$| = 1.2−4.1, we find a t_dep = 1−30 Myr, confirming the extremely short gas depletion time.

Resolved t_dep measurements have been carried out for only a handful of DSFGs (see Sharon et al. 2019 for a recent compilation), which generally show similarly strong (∼1 dex) variation in t_dep. The depletion times in SDP.81 are in the extreme starburst regime, an order of magnitude lower than seen in some other high-redshift DSFGs – e.g. GN20 (t_dep = 50−200 Myr; Hodge et al. 2015), J14011 (50–100 Myr;, Sharon et al. 2013), and J0911 (200–1000 Myr; Sharon et al. 2019). Cañameras 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³ M_⊙ yr⁻¹ kpc⁻²). 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.

4.3 [C ii]/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, 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 with the FIR (SFR) surface density (Herrera-Camus et al. 2015). Recently, resolved [C ii] studies in |$z$| > 1 DSFGs have been presented by Gullberg et al. (2018), Lamarche et al. (2018), Litke et al. (2019), and Rybak et al. (2019), revealing a pronounced [C ii]/FIR deficit down to |$L_\mathrm{[C\, \rm {\small{II}}]}/L_\mathrm{FIR}\simeq 10^{-4}$|⁠. Similarly pronounced [C ii]/FIR deficits have been found in resolved ALMA observations of |$z$| ≥ 5 quasar hosts (e.g. Neeleman et al. 2019; Venemans et al. 2019; Wang et al. 2019). Here, we consider both the source-averaged and resolved [C ii]/FIR ratio in SDP.81, and compare it to local and high-redshift studies.

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⁻³. 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 spatially averaged [C ii]/FIR measurements in high-redshift sources – as 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⁻⁴ in the two FIR-bright clumps – almost an order of magnitude lower than the spatially integrated [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). However, the [C ii]/FIR ratio decreases steeply towards the FIR-bright clumps, down to ∼3 × 10⁻⁴, indicating a strong [C ii]/FIR deficit. Similar radial trends have been observed in other resolved [C ii]/FIR observations at a high redshift (Gullberg et al. 2015; Lamarche et al. 2018; Litke et al. 2019; Rybak et al. 2019; Wang et al. 2019). As shown in 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 low [C ii/FIR ratios – these are likely lens-modelling artifacts.

Fig. 11 compares the [C ii]/FIR ratio in SDP.81 as a function of the SFR surface density Σ_SFR to other local and high-redshift sources, and the empirical relation of Smith et al. (2017). Resolved [C ii]/FIR deficit measurements in high-redshift DSFGs (Gullberg et al. 2018; Lagache, Cousin & Chatzikos 2018; Litke et al. 2019; Rybak et al. 2019) show a ≥1 dex scatter in [C ii]/FIR for a given Σ_SFR. While this might partially stem from systematic uncertainties in Σ_SFR estimates and the varying fraction of the [C ii] emission from ionized gas (cf.  Díaz-Santos et al. 2017), the observed [C ii]/FIR dispersion is likely due to the intrinsic scatter in galactic properties, such as gas and dust conditions. From the SDP.81 data, we estimate the best-fitting log Σ_SFR–log ([C ii]/FIR) slope of −0.75 ± 0.01, using orthogonal distance regression, similar to the slope of −0.7 derived by Lamarche et al. (2018) for SDP.11. This is much steeper than Smith et al.'s (2017) trend (−0.21 ± 0.03). We address the [C ii]/FIR slope when correcting for the fraction of gas in [C ii] and the physical origin of the [C ii]/FIR deficit in more detail in Section 5.4.

5 PDR MODELLING

5.1 Modelling setup

We adopt the line and FIR continuum predictions from the PDRToolbox models (Kaufman et al. 1999, 2006; Pound & 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 A_V = 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 on to a 0.25 dex grid, which is too coarse to estimate the uncertainties on the inferred properties. We re-sample them on to a finer grid (steps of 0.05 dex), using a degree 2 spline for interpolation.

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 contribution of ionized gas to [C ii] emission can be estimated using [N ii] 122/205-|$\mu$|m lines (Croxall et al. 2017; Herrera-Camus et al. 2018b); however, [N ii] 122-|$\mu$|m emission is undetected in the Herschel spectra with only weak upper limits (≤2.6 × 10⁻¹⁸ W m⁻²; Zhang et al. 2018, giving a PDR fraction of ≥0.6), while [N ii] 205-|$\mu$|m line falls outside the Herschel spectral coverage. Alternatively, using the empirical relation for the GOALS sample, the PDR contribution to the [C ii] can be estimated using the rest-frame 63- and 158-|$\mu$|m continuum (equation 4 of Díaz-Santos et al. 2017): for SDP.81, these correspond to the SPIRE 250-|$\mu$|m and ALMA Band 8 continuum (Table 2), yielding |$f_\mathrm{[C\, \rm {\small{II}}]}^\mathrm{PDR}=0.90\pm 0.10$|⁠. Given the available constraints, we adopt a conservative estimate of the |$f_\mathrm{[C\, \rm {\small{II}}]}^\mathrm{PDR}=0.8$|⁠. We discuss the sensitivity of the inferred PDR properties to this choice below.

Secondly, 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|$\mu$|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⁻¹ kpc⁻² (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⁴ cm⁻³, 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.

5.2 Global PDR model

We first apply the PDR modelling to the spatially integrated source-plane line and continuum luminosities. Fig. 12 shows the G and n values traces by individual line ratios.⁵ The maximum a posteriori model gives G = 10^3.1 ± 0.1 G₀, n = 10^4.9 ± 0.1 cm⁻³, with a reduced χ² = 0.96. The derived values are largely insensitive to [C ii] optical depth: for an optically thick [C ii] emission, the inferred G, n values change by + 0.5 and −0.1 dex, respectively. Increasing the fraction of [C ii] emission originating in PDR regions from 80 to 100 per cent results in G, n changing by ≤0.1 dex.

Compared to Valtchanov et al.'s (2011) global PDR model for SDP.81 [inferred using the same PDRToolbox models based on [O i], [C ii], 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); and (2) we are considering CO J_upp ≥ 3 lines, which trace denser molecular gas that CO(1–0) in Valtchanov et al. (2011; see Section 5.3 for a detailed discussion of systematics).

To assess the impact of differential magnification on the global PDR model, we apply the same analysis to the sky-plane (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.

5.3 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 the following:

• 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 (FIR continuum) and 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 G₀/n(H)≤10⁻³, which would imply an unphysically large source (cf. Díaz-Santos et al. 2017; Wardlow et al. 2017). Fig. 13 shows the G₀ 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 G₀, and high density n(H) = 10^5.0−10^5.5 cm⁻³, 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 G₀. While the starburst 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 T_PDR varies from about 150 K in the outskirts to 1500 K in the FIR-bright clumps. However, with T_PDR S/N ≤ 5 across the source, we can not robustly measure the spatial variation of T_PDR. The full implications of the high T_PDR 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⁴ G₀, n ∼ 10⁴−10⁵ cm⁻³ in the central regions of two |$z$| = 2.94 DSFGs from the ALESS sample; SDP.81 shows slightly higher G and lower n(H).

Compared to the PDR properties of the present-day ULIRGs (GOALS sample; Díaz-Santos et al. 2017) – inferred from the [O i] 63-|$\mu$|m and [C ii] lines and FIR continuum – both the global and resolved PDR properties of SDP.81 show density almost ∼2 dex higher, with somewhat (∼0.5 dex) stronger FUV fields.

Finally, we consider star-forming regions in the Milky Way and the Large Magellanic Cloud⁶ from Oberst et al.'s (2011) compilation. The resolved PDR properties of SDP.81 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.'s 1999, 2006 PDR models used in our analysis), Goicoechea et al. (2015) obtain G ≃ 3 × 10⁴ G₀, n ≥ 10⁵ cm⁻³ for the region closest (R ≤0.1 kpc) to the Trapezium cluster, decreasing down to G ≃ 10³ G₀, n ≃ 10⁴ cm⁻³. These are, respectively, on the upper and lower end of the conditions seen in the FIR-bright 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 star-forming regions in DSFGs.

5.3.1 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 |$\mu$|m rest-frame), [C ii], and CO (J_upp ≥ 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 χ² ≤ 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 overinterpreting 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 FIR-bright clumps(reduced χ² = 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 star-forming environment expected in SDP.81, where strong gradients in ISM density and composition along the line of sight might be expected. These can be properly accounted for by three-dimensional 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 (J_upp ≤ 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 Le Petit et al.'s (2006) code to the resolved observations of 30 Doradus by Chevance et al. 2016 – might provide a better description of molecular clouds in an 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 J_upp ≥ 3 CO lines (critical density n_crit ≥ 8 × 10³ cm⁻³), we effectively discount the extended cold, low-density (n ≤ 10³ cm⁻³) gas phase traced by CO(1–0) and CO(2–1) lines. Indeed, studies of the CO ladders including the J_upp = 1−2 lines have found significant cold gas components both in |$z$| ∼ 0 ULIRGs (e.g. Kamene(e.g. Greve et al. 2014; Kamenetzky 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, Yang et al.'s (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 FUV 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.

5.4 What drives the [C ii]/FIR deficit in SDP.81?

With the inferred G, n, and T_PDR 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, and (3) AGN contribution.

5.4.1 Reduced photoelectric heating

5.4.2 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 disc. 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 (cf. 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⁻¹⁵ erg s⁻¹ cm⁻², 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 X-ray luminosity L_{0.5–8.0 keV} ≃ 2 × 10⁴⁴ erg s⁻¹ (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⁴²−3 × 10⁴⁴ erg s⁻¹; Wang et al. 2013).⁷ Even assuming the highest column density for the ALESS sources from Wang et al. (2013) (N_H ≃ 5 × 10²⁴ cm⁻³), a comparison with the carbon ionization model in the vicinity of X-ray 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 (cf. 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 AGN-associated X-ray/cosmic ray heating.

5.4.3 C⁺ thermal saturation

Given the high PDR surface temperature inferred from PDR modelling (Section 5), the good match between the [C ii] abundance-corrected [C ii]/FIR slope with Muñoz & Oh's (2016) predictions, and the limited effects of the positive grain charging and AGN heating, we consider the thermal saturation of the C⁺ fine-structure levels to be the most likely cause of the strong [C ii]/FIR deficit in SDP.81.

6 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 |$\mu$|m) continuum emission in a |$z$| = 3.042 gravitationally lensed DSFG SDP.81. Using a visibility-fitting lens modelling code with an adaptive source-plane grid (R15a), we achieve a median source-plane resolution of ∼200 pc. Combining these with the CO(5–4) and CO(8–7) data from R15b, we use PDR models of Kaufman et al. (1999, 2006) and 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 as follows:

• The ALMA Band 8 continuum source-plane surface brightness distribution is consistent with the previously published Band 6 and Band 7 continuum reconstructions (Dye et al. 2015; R15a; 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(J_upp ≥ 5) emission: ∼50 per cent of [C ii] emission originates outside the FIR-bright region. The [C ii] maps reveal an extended, low-surface-brightness 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 time-scale of t_dep = 10⁶−2 × 10⁷ Myr, significantly shorter than in most other resolved t_dep studies in high-redshift 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 FIR-bright clumps.

• Comparing the source-integrated (unresolved) de-lensed [C ii], CO, and FIR emission with PDR models of Kaufman et al. (1999, 2006) and Pound & Wolfire (2008), we find FUV field strength G = 10³ G₀ and gas density n(H) = 10⁵ cm⁻³, 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 G₀). In contrast, n(H) varies only marginally (10^5.0−10^5.5 cm⁻³). 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_\mathrm{[C\, \rm {\small{II}}]}/L_\mathrm{FIR}=3\times 10^{-4}$|⁠). Although the [C ii] surface brightness generally increases towards the FIR-bright 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 emphasize 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, for example, different heating mechanisms on sub-kpc scales.

This paper presents a first look into the PDR properties of high-redshift DSFGs 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 PDRs 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 (Figs 2 and A2) publicly available at sdp81.strw.leidenuniv.nl. We hope this legacy data set will facilitate independent analysis of the conditions in high-redshift DSFGs on sub-kpc scales.

ACKNOWLEDGEMENTS

The authors thank the anonymous referee for their thorough comments and M. Kazadjian and F. Israel for helpful discussions on the radiative transfer modelling, Ana Lopez-Sepulcre for help with ALMA observing setup, and C. Lamarche and K. Litke for sharing their [C ii]/FIR data. MR thanks MPA Garching, and S. D. M. White in particular, for the access to computational facilities.

MR and JH acknowledge the support of the VIDI research programme with project number 639.042.611, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO). LG acknowledges support from the Amaldi Research Center funded by the MIUR program ‘Dipartimento di Eccellenza’ (CUP:B81I18001170001). This paper makes use of the following ALMA data: ADS/JAO.ALMA #2011.0.00016.SV, #2016.1.00633.S, and #2016.1.01093.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. Fig. 14 utilizes the corner package by Foreman-Mackey (2016).

Footnotes

REFERENCES

APPENDIX A: LENSING RECONSTRUCTIONS

Figs A1 and A2 show the lens modelling results for the Band 8 continuum, velocity-integrated CO(3–2) and CO(10–9) lines, and the reconstructions for individual 40 km s⁻¹ [C ii] velocity channels. The dirty images are calculated using natural weighting (i.e. each visibility is weighted by 1/σ²).

APPENDIX B: ESTIMATING THE SYSTEMATIC ERROR ON THE SOURCE RECONSTRUCTIONS

We reconstructed the lensed source using source-plane regularization to impose a spatial correlation in the source plane and prevent noise fitting. However, for low-S/N data, the optimization process might yield high levels of regularization, overly smoothing the source. Here, we assess the systematic uncertainty on the reconstructed source by creating mock ALMA observations for each tracer.

We consider two different source-plane surface brightness distributions: (1) a uniform surface-brightness elliptical source with major/minor axis of 3.4 × 2.4 kpc (corresponding to the FIR-bright source size at 45° inclination) and (2) a surface-brightness distribution matching the FIR luminosity map. These sources are then forward-lensed using the lens model for SDP.81. We then create mock ALMA observations for each tracer by setting the sky-plane flux density to the observed value (see Table 1), using the same uv-plane coverage and adding the observed noise per baseline. Finally, we perform the full lens-modelling procedure for the each data set. The reconstructed sources are shown in Figs B1 (uniform input source) and B2 (FIR-like input source).

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

Secondly, 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 |$\Sigma _\mathrm{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 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 R15b, 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 artefact.

• [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 artefact.

APPENDIX C: SENSITIVITY OF G, n TO INFERRED LINE RATIOS

We now assess the robustness of G and n inferred from the PDR modelling to the bias in the reconstructed source-plane brightness distribution. In other words, if the lens modelling introduces a bias in the reconstructed source, how much will this affect the inferred G and n? Fig. C1 shows the G–n space isocontours corresponding to one and two times of the mean line ratios in SDP.81, using the PDRToolbox models (Kaufman et al. 1999). As the systematic error the inferred line ratios is less than unity (Appendix  B), a change by a factor of 2 thus overestimates the actual uncertainty. As the line ratios considered are either between different species or non-neighbouring CO rotational transitions, a change by a factor of 2 results in relatively minor (≤0.5 dex) shifts in the G/n(H) space. The only exception is the CO(5–4)/(3–2) ratio; however, given the relatively large associated uncertainty (Fig. 12), we do not expect the inferred PDR properties to be significantly biased.