Gas and star formation from HD and dust emission in a strongly lensed galaxy

ABSTRACT

The molecular gas content of high-redshift galaxies is a highly sought-after property. However, H₂ is not directly observable in most environments, so its mass is probed through other emission lines (e.g. CO, [C i], [C ii]), or through a gas-to-dust ratio. Each of these methods depends on several assumptions, and are best used in parallel. In this work, we extend an additional molecular gas tracer to high-redshift studies by observing hydrogen deuteride (HD) emission in the strongly lensed z = 5.656 galaxy SPT0346−52 with ALMA. While no HD(1–0) emission is detected, we are able to place an upper limit on the gas mass of |$\rm M_{H_2}\lt 6.4\times 10^{11}\, M_{\odot }$|⁠. This is used to find a limit on the |$L^{\prime }_{\mathrm{ CO}}$| conversion factor of      |$\rm \alpha _{CO}\lt 5.8\,M_{\odot}(K\,km\,s^{-1}\,pc^2)^{-1}$|⁠. In addition, we construct the most complete spectral energy distribution of this source to date, and fit it with a single-temperature modified blackbody using the nested sampling code multinest, yielding a best-fitting dust mass M_dust = 10^{8.92 ± 0.02} M_⊙, dust temperature 78.6 ± 0.5 K, dust emissivity spectral index β = 1.81 ± 0.03, and star formation rate SFR = 3800 ± 100 M_⊙ yr⁻¹. Using the continuum flux densities to estimate the total gas mass of the source, we find   |$\rm M_{H_2}\lt 2.4\times 10^{11}\, M_{\odot }$| , assuming subsolar metallicity. This implies a CO conversion factor of α_CO < 2.2, which is between the standard values for MW-like galaxies and starbursts. These properties confirm that SPT0346−52 is a heavily starbursting, gas-rich galaxy.

ISM: molecules, galaxies: high-redshift, galaxies: starburst

1 INTRODUCTION

Molecular hydrogen (H₂) is both the most abundant molecule in the Universe and the main fuel for star formation (e.g. Kennicutt & Evans 2012). However, direct observations of H₂ are difficult, as the symmetry inherent in its structure dictates that line emission only originates from energetic environments. To work around this limitation, emission from other molecules may be used as a tracer.

The main alternative is the second most abundant molecule, carbon monoxide (CO), which is 10⁴ times less abundant than H₂ (e.g. Bolatto, Wolfire & Leroy 2013). The mass of molecular hydrogen in a source may be estimated by observing CO (J = 1–0) emission and assuming a conversion factor |$\alpha _{\mathrm{ CO}}=M_{\mathrm{ H_2}}/L^{\prime }_{\mathrm{ CO}}$|[ = ]M_⊙(K km s⁻¹ pc²)⁻¹, which is assumed to be ∼0.8 for starburst galaxies and ∼4.6 for relatively quiescent, Milky Way-like galaxies (e.g. Daddi et al. 2010). However, these conversion factors are strongly metallicity dependent (e.g. Narayanan et al. 2012), and observations of higher-J lines must be converted to the J = 1–0 line by the assumption of an uncertain factor (Carilli & Walter 2013). Alternatively, one may observe emission from warm dust, convert this luminosity to a dust mass (e.g. Lamperti et al. 2019), and convert this to a gas mass by assuming a gas-to-dust conversion factor (δ_gd, e.g. Saintonge et al. 2013).

A third alternative, which is still largely unexploited, is to observe hydrogen deuteride (HD), which may be converted directly into an H₂ mass using the HD/H₂ ratio at the observed redshift (e.g. Bergin et al. 2013). In regions where the gas is mainly molecular, i.e. when hydrogen is in molecular form, deuterium is mainly locked in HD. HD is rapidly photodissociated in diffuse clouds (Spitzer et al. 1973) and at the edge of photodissociation regions (PDRs), where the visual extinction falls below ∼2–3 mag (Jansen et al. 1995). However, virtually all of the deuterium is expected to be contained within HD in dense molecular clouds in our Galaxy (Tielens 2005) and beyond, as it can efficiently form both in the gas phase via H₂ + D⁺ at low metallicities and on the surface of dust grains even at dust temperatures above 150 K (Watson 1973; Cazaux & Spaans 2009).

HD has been regularly detected via absorption of electronic transitions towards stars in the Milky Way (e.g. Spitzer et al. 1973; Lacour et al. 2005) and towards other galaxies (e.g. Noterdaeme et al. 2008, 2010; Balashev, Ivanchik & Varshalovich 2010; Ivanchik et al. 2010; Oliveira et al. 2014; Daprà et al. 2017). Rotational transitions of HD, which can provide mass estimates, have not yet been observed in external galaxies. Only in the Milky Way there have been detections of the HD J = 1–0 ground state rotational transition at 112.2 μm towards the PDR Orion Bar using the Infrared Space Observatory ISO (e.g. Wright et al. 1999), and towards protoplanetary discs using the Herschel Space Observatory (Bergin et al. 2013; McClure et al. 2016). Neufeld et al. (2006) also detected the excited HD J = 4–3 and 5–4 transitions (both good pressure probes) at 28.5  and 23.1 μm, respectively, towards a supernova remnant using the Spitzer Space Telescope.

Bergin et al. (2013) used HD(1–0) to derive an important lower limit of the TW Hya protoplanetary disc mass, taking into account the fact that some atomic D could be hidden in molecular ices instead of HD; they found significantly larger masses than previously derived with CO (and rare CO isotopologue) lines and with dust continuum emission. It is therefore important to extend the study of HD rotational transitions to external galaxies and compare the new derived values of the gas masses with those obtained with the more classical CO and dust emission methods. High-redshift galaxies allow us to search for the lowest J transitions of HD using sensitive ground based telescopes, such as ALMA. As clearly shown by Bergin et al. (2013), the bulk molecular gas reservoir of molecular hydrogen can be directly assessed using HD(1–0), once the physical structure (volume density and kinetic temperature) is known. Unlike CO, HD is expected to have a constant abundance relative to H₂ throughout the dense molecular material and emission is expected as long as the gas temperature is above ∼10–15 K. Therefore, with a good knowledge of the (average) density and temperature towards a high-redshift galaxy, HD(1–0) provides a good alternative method to estimate the gas mass in these distant objects. To test this, we search for emission from the fundamental rotational transition HD(J = 1–0) from the energetic, strongly lensed galaxy SPT 0346−52.

SPT-S J034640−5204.9 (hereafter SPT 0346−52) is a strongly lensed dusty star-forming galaxy at z = 5.656, first studied in the ALMA survey of Weiß et al. (2013) and Vieira et al. (2013). Detailed gravitational lens modelling shows that the galaxy is magnified by a factor μ = 5.6 ± 0.1 (Spilker et al. 2016), while a source-plane reconstruction strongly resembles either an ongoing major merger or rotating disc (Spilker et al. 2015; Dong et al. 2019; Litke et al. 2019). The detection of a central outflow coincident with an extreme starburst rejects the merger scenario, adding additional credence to the rotator scenario (Jones et al. 2019). Spectral energy distribution (SED) modelling yields a massive star formation rate (SFR ∼ 4000–5000 M_⊙ yr⁻¹, Ma et al. 2015) and further observations have found substantial [C ii] 158 μm (⁠|$L_{[\mathrm{ C\, II}]}=(5.0\pm 0.7)\times 10^{10}\, \mathrm{ L}_{\odot }$|⁠, Gullberg et al. 2015) and CO (⁠|$L_{\mathrm{ CO}(2-1)}=(2.4\pm 0.2)\times 10^{8}\, \mathrm{ L}_{\odot }$|⁠, Aravena et al. 2016) emission. However, X-ray and radio observations have not revealed the presence of any AGN and suggest that this source is mainly powered by star formation (Ma et al. 2016).

In this work, we present new ALMA observations in band 8, resulting in an upper limit on the intensity of HD(1–0) emission. In addition, we add our new continuum flux densities to those previously published, and fit the resulting SED with a blackbody submm to far-infrared (FIR) model, yielding new estimates on the star formation rate, dust temperature, dust mass, and several limits on the molecular gas mass, resulting in limits for the CO-to-H₂ conversion factor. We will assume (Ω_Λ, Ω_m, h) = (0.692, 0.308, 0.678) (Planck Collaboration et al. 2016) throughout.

2 OBSERVATIONS AND DATA REDUCTION

Using band 8 of ALMA, we observed SPT 0346−52 from 2016 October 24 to November 1, using 41–46 antennas. Out of a total observation time of 12.9 h, 6.4 h were on source. Our complex gain, bandpass, and flux calibrators were J0253−5441, J0522−3627, and J0519−4546, respectively. The data were originally calibrated by ALMA staff, but excessive flagging of edge channels resulted in a gap in the resulting frequency coverage. To correct this, we recalibrated the data, flagging only three channels at each edge of the spectral windows (SPWs), rather than the default of 10 channels.

Our frequency range was covered by two sidebands, each composed of two SPWs, made of 128 channels of width 15.625 MHz. The lower sideband (387.562–391.031 GHz) was expected not to contain prominent emission or absorption lines, while the redshifted frequency of the HD(1–0) transition (ν_rest = 2674.986094 GHz, ν_obs = 401.891 GHz) falls in the upper sideband (399.757–403.351 GHz).

Since the goal of these observations is to detect HD(1–0) emission, we first create full (i.e. line and continuum) data cubes using the casa task tclean, natural weighting, and a clean threshold of 3σ, where σ is the rms noise level per channel of the dirty image. We explore both the native channel width and 2, 3, 4, and 5-channel averaging. None of these data cubes shows an obvious HD(1–0) signal, so we proceed with the native channel width.

Since the emission line is either weak or non-detected, we perform continuum subtraction, in order to isolate any low-level signal. This process is non-trivial, as the atmospheric transmission across our sidebands varies dramatically¹ and the velocity width of the HD(1–0) line is unknown. While the lower (continuum-only) sideband would ideally be used to estimate the continuum level, it shows only |$30{{\ \rm per\ cent}}$| transmission (assuming 2.0 mm PWV), while the upper (HD) sideband has a more favorable transmission (⁠|$\sim 50{{\ \rm per\ cent}}$|⁠), but is plagued by multiple atmospheric absorption lines.

In order to separate the possible line and continuum emission in these data, we explored multiple continuum-subtraction techniques. First, we used the casa task UVCONTSUB to fit and subtract a first-order polynomial continuum model directly to the visibilities of both sidebands, resulting in purely line emission. Since this task fits the continuum emission directly in the visibilities, it is independent of user-provided imaging parameters. When performing this fitting, we avoided the three frequency ranges of atmospheric absorption in the upper sideband (i.e. ∼400.0 GHz, ∼401.3 GHz, ∼402.4 GHz) and all channels that could include HD line emission. Since the linewidth of HD is unknown, we adopted a conservative line channel range of the expected HD frequency ±750 km s⁻¹, based on FWHM_CO = 613 ± 30 km s⁻¹ (Aravena et al. 2016), resulting in theoretical line channels of 400.888–402.899 GHz. Using the casa task tclean, the resulting continuum-subtracted visibilities were then imaged using natural weighting and a clean threshold of 3σ. This resulted in an obvious undersubtraction of the continuum, so this continuum-subtraction method was disregarded.

As an alternative to visibility-space continuum subtraction, we also explored the use of image-space subtraction by applying the casa task imcontsub (with both zero- and first-order polynomial fits) to a full data cube spanning both sidebands. This resulted in a nearly identical data cube as the above UVCONTSUB approach and was not used.

The above tests reveal that the lower sideband is unsuitable for continuum fitting, and that our ‘conservative’ channel exclusion is too rigorous. In order to correct this, we re-examine the imcontsub approach (both zero and first-order) for the upper sideband, but only excluding ±250 km s⁻¹ on either side of the HD(1–0) line. The first-order approach returns a flatter spectrum, so we proceed with this continuum-subtracted cube.

Our final continuum-subtracted cube has a synthesized beam of 0.19 arcsec × 0.17 arcsec at a position angle of 57°, channels of width 15.625 MHz (∼12 km s⁻¹), and an rms noise level per channel of 0.25 mJy beam⁻¹.

A continuum image was created by applying tclean to all line- and atmospheric feature-free channels of both sidebands (i.e. the ‘conservative’ approach of above), natural weighting, multifrequency synthesis, and a clean threshold of 3σ, where σ is the rms noise level per channel of the dirty image, resulting in a final rms noise level in the cleaned image of 0.11 mJy beam⁻¹.

3 RESULTS

3.1 Continuum

The resulting continuum image is shown in Fig. 1. The continuum is well detected, with a maximum significance of ∼185σ. The total continuum flux density is 171.5 ± 1.0 mJy. Note that in this section, we will state the observed (i.e. image-plane) continuum parameters, rather than the delensed (or source-plane) values.

Figure 1.

The observed continuum (⁠|$\rm \lambda _{rest}\sim 114\, \mu m$|⁠), created using all line-free channels in both sidebands. Contours are shown in a geometric sequence: [10, 20, 40, 80, 160] × σ, where σ = 0.11 mJy beam⁻¹. The synthesized beam (0.19 arcsec × 0.17 arcsec, with major axis position angle  = 55°) is shown by the solid white ellipse to the lower left.

A 2D Gaussian fit to the northern component in the continuum image returns a deconvolved size of (0.39 ± 0.02) arcsec × (0.17 ± 0.01) arcsec at (72 ± 2)°, an integrated flux density of 56.3 ± 2.7 mJy, and a peak flux density of 17.9 ± 0.7 mJy beam⁻¹. A similar fit to the southeastern component returns a deconvolved size of (0.28 ± 0.02) arcsec × (0.16 ± 0.02) arcsec at (27 ± 6)°, an integrated flux density of 33.6 ± 1.9 mJy, and a peak flux density of 14.0 ± 0.6 mJy beam⁻¹. Although the southwestern component is extended in an arc, a Gaussian fit returns a deconvolved size of (0.68 ± 0.05) arcsec × (0.19 ± 0.02) arcsec at (117 ± 1)°, an integrated flux density of 56.7 ± 3.5 mJy, and a peak flux density of 9.9 ± 0.5 mJy beam⁻¹. Due to lensing effects, these fits are not trivially translatable to physical parameters (e.g. size, surface brightness). However, since the integrated flux density of each fit is higher than its peak, and the spatial scales have small errors, we may state that these sources are well resolved.

Using the 2σ contour as a spatial mask, we determine the total |$\lambda _{\rm rest}\sim 114\, \mu$|m continuum flux density of this source to be 171.5 ± 1.0 mJy. This value is greater than the sum of each integrated flux density (147 ± 5 mJy), suggesting that the diffuse emission between the components (especially in the south) is significant, or that these resolved components are poorly described by Gaussian fits.

3.2 HD(J = 1–0) emission

Since SPT 0346−52 is strongly detected in multiple gas tracers (i.e. CO(2–1), Aravena et al. 2016; FIR continuum emission), it is expected to contain a substantial amount of molecular gas (∼10¹¹ M_⊙), and should therefore be observable in HD(1–0) emission (e.g. Bergin et al. 2013). To explore this possibility, we examine a continuum-subtracted cube (see Section 2 for details of cube creation), searching for significant line emission at the expected redshift. However, the systemic velocity, spatial position, velocity width, and spatial extent of HD(1–0) are not known a priori, making this search non-trivial. As a further complication, multiple atmospheric absorption features are present in the observed frequency range (see Section 2), resulting in a suboptimal continuum subtraction and a varying rms noise level. By searching the cube, two tentative line features are detected, but they are believed to be noise. For details of these tentative detections (‘Tentative 1’ and ‘Tentative 2’), see Appendix  A.

This exploration of the data cube yielded no believable detections, so we turn to the possibility that the HD(1–0) emission is spread over many channels (i.e. broad-line width) with a low amplitude. Indeed, previous observations of line emission in SPT 0346−52 found FWHM values of ∼500–700 km s⁻¹ (e.g. Aravena et al. 2016; Apostolovski et al. 2019; Dong et al. 2019), which correspond to ∼40–60 channels in our data cube. If a weak emission feature is distributed over many channels, it is entirely possible that it would not be detected in a channel-by-channel inspection. To test this, we use the casa task immoments to collapse the channels corresponding to |$\rm -350\lt v\lt 350$| km s⁻¹ and search for significant features. This collapsed image shows a 3σ feature that is coincident with the northern FIR continuum image of SPT 0346−52, but this feature is comparable to noise peaks in the field of view, and is thus likely not significant (‘Tentative 3’, Appendix  A).

Instead, it is plausible that the underlying emission is very spatially extended, on a similar scale as the FIR continuum emission. As shown in Fig. 1, the FIR continuum is well resolved by our observations, so the emission is spread over multiple observing beams. If HD(1–0) is weak and as extended as the continuum emission, then it is possible that it may only be detected by integrating the emission from a wide area. To test this, we extract a spectrum from the continuum-subtracted cube using the 2σ contour of the continuum map as a mask, resulting in the spectrum shown in Fig. 2. A weak feature is detected, but it is at high velocity, and is thus unlikely to be real (‘Tentative 4’, Appendix  A).

Figure 2.

Atmospheric transmission curve for 2.0 mm PWV (top) and spectrum taken over the 2σ continuum contour (bottom). Grey shaded area in lower panel shows rms noise level. Velocity scale is shown relative to the expected frequency of HD(1–0). No lensing correction has been applied.

Since no definite emission is detected through inspecting the data cube, collapsing a wide channel range, or extracting a spectrum from a large area, we conclude that HD(1–0) is undetected in these observations. These observations may be used to place a conservative upper limit on the HD(1–0) flux from SPT 0346−52. Assuming the same spatial distribution as the FIR continuum emission, we extract a spectrum from our continuum-subtracted cube (Fig. 2), finding an rms noise level of ∼2.15 mJy at the expected frequency. Next, we assume a wide Gaussian line profile (FWHM_HD = 613 km s⁻¹ ∼ 0.8 GHz; Aravena et al. 2016) with an amplitude limit of <2σ, resulting in a conservative upper limit on the integrated intensity of HD(1–0) of F_L < 3.8 × 10⁻²⁰ W m⁻². The implications of this non-detection will be discussed in Section 4.2.1.

4 ANALYSIS

4.1 SED modelling

Since SPT 0346−52 is one of the highest-redshift sources detected in the SPT-S survey, and exhibits one of the lowest magnifications in the SPT sample, its intrinsic luminosity is substantial. Due to its extreme characteristics, it has been observed in continuum emission by a host of instruments (e.g. HST, ACT, SPIRE). We compile these observations in Table 1.

Table 1.

Continuum flux densities measured for SPT 0346−52, with no magnification correction. Includes previously published values, values from the ALMA archive, and new values from this paper. Bold entries were used in our SED model. (1) Ma et al. (2015), (2) This work, (3) 2015.1.01580.S, (4) Spilker et al. 2016, (5) Jones et al. 2019, (6) 2013.1.01231.S, (7) 2016.1.00654.S, (8) Marsden et al. 2014, (9) Apostolovski et al. 2019, (10) Dong et al. 2019, (11) Aravena et al. 2016, (12) Ma et al. 2016.

λobs [μm]	λrest[μm]	Sν [mJy]	Instrument	Reference
1.1	0.17	<3.8 × 10−4	HST/WFC3	1
1.6	0.24	<9.1 × 10−4	HST/WFC3	1
3.6	0.54	<2.4 × 10−3	Spitzer/IRAC	1
4.5	0.68	<3.6 × 10−3	Spitzer/IRAC	1
100	15	<6.0	Herschel/PACS	1
160	24	33 ± 9	Herschel/PACS	1
250	38	122 ± 11	Herschel/SPIRE	1
350	53	181 ± 14	Herschel/SPIRE	1
500	75	204 ± 15	Herschel/SPIRE	1
758	114	171.5 ± 1.0	ALMA	2
823	124	123.9 ± 1.0	ALMA	3
870	131	123.0 ± 13.0	ALMA	4
894	134	82.3 ± 1.7	ALMA	5
1028	155	79.7 ± 1.2	ALMA	6
1320	198	43.9 ± 0.4	ALMA	7
1375	207	47.1 ± 4.3	ACT	8
1400	210	46.0 ± 6.8	SPT	1
2000	301	8.80 ± 1.35	ALMA	9
2026	304	16.7 ± 2.6	ACT	8
2064	310	14.75 ± 0.03	ALMA	10
3000	451	3.06 ± 0.05	ALMA	9
8081	1214	0.16 ± 0.02	ATCA	11
54 508	8189	<0.114	ATCA	12
142 758	21 448	<0.213	ATCA	12

In an effort to fill in the |$\lambda _{\rm obs}\sim 1\,$| mm regime, we have also compiled multiple continuum detections from the ALMA data archive. For each of these measurements, we use the ALMA-staff calibrated continuum images (i.e. no further calibration has been performed). These continuum measurements may be combined to create a SED, which can be fit with a dust model.

4.1.1 Model description

We choose a modified blackbody (MBB) (also known as a greybody; e.g. Greve et al. 2012):

$$\begin{eqnarray*} S_{\mathrm{ obs}}(\nu _{\mathrm{ obs}})=\frac{\Omega }{(1+z)^3}B(\nu ,T_{\mathrm{ dust}})(1-\mathrm{ e}^{-\tau _{\nu }}), \end{eqnarray*}$$

(1)

where S_obs(ν_obs) is the observed flux density at ν_obs = ν_rest/(1 + z), |$\Omega =(1+z)^4A_{\mathrm{ gal}}/D_\mathrm{ L}^2$| is the solid angle of the galaxy with area A_gal at luminosity distance D_L, B(ν, T_dust) is the blackbody function with a dust temperature T_dust, and τ_ν is the optical depth:

$$\begin{eqnarray*} \tau _{\nu }=\Sigma _{\mathrm{ dust}}\kappa _{\nu }=\frac{M_{\mathrm{ dust}}}{A_{\mathrm{ gal}}}\kappa _o(\nu /\nu _o)^{\beta }, \end{eqnarray*}$$

(2)

where we assume the dust absorption coefficient κ_o = 0.04 m² kg⁻¹ at ν_o = 250 GHz (Beelen et al. 2006). It is often assumed that the emission is optically thin (τ_ν < <1), which allows equation (1) to be greatly reduced (e.g. Casey 2012). However, we find that this assumption is not applicable for our data (see Section 4.1.3), and proceed with the general form.

Due to the high redshift of this source, we include the effects of the cosmic microwave background (CMB) on the observed dust continuum, as outlined in da Cunha et al. (2013) and implemented in Carniani et al. (2019). First, the higher temperature CMB will heat the dust, resulting in a correction on our derived dust temperature:

$$\begin{eqnarray*} T_{\mathrm{ dust}}^{\prime }=\left[T_{\mathrm{ dust}}^{4+\beta }+T_o^{4+\beta }\left((1+z)^{4+\beta }-1\right)\right]^{1/(4+\beta)}, \end{eqnarray*}$$

(3)

where T_dust is the true dust temperature, |$T_{\mathrm{ dust}}^{\prime }$| is an effective dust temperature, T_o = 2.73 K is the CMB temperature at z = 0, and β is the dust emissivity spectral index. In addition to this effect, the hotter CMB provides a background against which we observe the dust emission. This contribution must be removed, resulting in an MBB function:

$$\begin{eqnarray*} B^{\prime }(\nu ,T_{\mathrm{ dust}})&=&B(\nu ,T_{\mathrm{ dust}}^{\prime })-B(\nu ,T_{\mathrm{ CMB}})\nonumber \\ &=& \frac{2h\nu ^3}{c^2}\left[\frac{1}{\mathrm{ e}^{h\nu /k_BT_{\mathrm{ dust}}^{\prime }}-1}-\frac{1}{\mathrm{ e}^{h\nu /k_BT_{\mathrm{ CMB}}}-1} \right], \end{eqnarray*}$$

(4)

where T_CMB = (1 + z)T_o.

Combining these equations, we find the following equation for the observed flux density:

$$\begin{eqnarray*} S_{\mathrm{ obs}}(\nu _{\mathrm{ obs}})=\frac{(1+z)\pi R_{\mathrm{ gal}}^2}{D_\mathrm{ L}^2}B^{\prime }(\nu ,T_{\mathrm{ dust}})\left(1-\mathrm{ e}^{\frac{-M_{\mathrm{ dust}}\kappa _o(\nu /\nu _o)^{\beta }}{\pi R_{\mathrm{ gal}}^2}} \right). \end{eqnarray*}$$

(5)

To approximate the radius of emission, we average the source-plane λ_obs = 2.0  and 3.0 mm effective radii of Apostolovski et al. (2019), resulting in R_gal ∼ 0.76 kpc.

Using equation (5), it is possible to fit the full FIR dust SED with only three free parameters: the dust emissivity spectral index (β), dust temperature (T_dust), and dust mass (M_dust).

4.1.2 Model implementation

Models were fit to the SED using the Bayesian inference code multinest (Feroz & Hobson 2008) and its python wrapper (pymultinest; Buchner et al. 2014). This code returns the most likely parameter values, parameter value probability distributions, and covariance distributions for each parameter pair. We assume uniform priors for each variable, with log₁₀(M_dust/M_⊙) = 8 − 10, β = 1 − 3, and T_dust = 10–200 K.

However, we must first consider what subset of our sample we may fit with this model. Since we are not including contributions from synchrotron or thermal free–free emission (e.g. Yun & Carilli 2002), we do not include the low-frequency λ_obs = 8 mm data point.

On the high-frequency side of the model, we choose to include all of the Herschel data (λ_rest = 24–75 μm), which trace the peak of the dust emission. Previous implementations of a similar model to fit dust SEDs of HyLIRGs (Yun & Carilli 2002; Wagg et al. 2014) used data from beyond the thermal dust peak. In addition, the SEDs of Ma et al. (2016) and Apostolovski et al. (2019) included these points, which were well fitted by dust-only models.

Each observed flux density was corrected for magnification by assuming a constant lensing magnification of μ = 5.6, based on the detailed modelling of |$\lambda _{\mathrm{ obs}}=870\, \mu$|m continuum emission (Spilker et al. 2016).

4.1.3 Model results

Our best-fitting model is shown in Fig. 3, the resulting values are listed in Table 2, and the associated covariance plots are shown in Fig. 4. Fitting an MBB model to our combined data set yields best-fitting values of T_d = 78.6 ± 0.5 K, β = 1.81 ± 0.03, and log(M_dust/M_⊙) = 8.92 ± 0.02. This model exhibits an FIR luminosity (i.e. 42.5–122.5 μm) of L|$_{\rm FIR}=(2.2\pm 0.1)\times 10^{13}\,$| L_⊙. Using the scaling relation from Kennicutt (1998b), this corresponds to an SFR = 3800 ± 100 M_⊙ yr⁻¹.

Figure 3.

The submm-FIR SED for SPT0346−52, showing points from Table 1. Data used in the fit are shown in black, while the one excluded low-frequency point and upper limits from literature (which are not used in the fit) are shown by open red circle and triangles, respectively. The flux densities were corrected for lensing magnification by assuming μ = 5.6 ± 0.1 (Spilker et al. 2016). The best-fitting MBB model to the sample is shown in black. The frequency range used to measure L_FIR is denoted by the vertical black dashed lines.

Figure 4.

The covariance distributions of the three parameters (M_dust, β, and T_d) used to fit our SED. The top plot of each column is the overall posterior probability distribution for each parameter.

Table 2.

Best-fitting values from an MBB fit to the FIR SED (top) and the resulting FIR luminosity and SFR (below).

Tdust [K]	78.6 ± 0.5
β	1.81 ± 0.03
log10(Mdust/M⊙)	8.92 ± 0.02
LFIR [L⊙]	(2.2 ± 0.1) × 1013
SFR [M⊙ yr−1]	3800 ± 100

Using the code cigale (Noll et al. 2009), Ma et al. (2015) fit a mid-infrared (MIR) to FIR (without the ALMA, ACT, or ATCA flux densities, but including the high-frequency upper limits) SED of this source, and derived an SFR of |$4840^{+1090}_{-890}$| M_⊙ yr⁻¹. The same study used a derived IR luminosity and the conversions of Kennicutt (1998a) and Murphy et al. (2011) to find SFR_IR = 3830–5340 M_⊙ yr⁻¹. Our value of SFR is only slightly (i.e. <2σ) lower than these estimates.

However, our value of β is >3σ lower than the regularly assumed value of β = 2.0 (Weiß et al. 2013; Gullberg et al. 2015; Ma et al. 2015). Using β = 2.0, previous investigations have yielded dust temperatures of 53 ± 5 K (Weiß et al. 2013) and 52 ± 2 K (Gullberg et al. 2015), which are lower than our value. Since both of these investigations fit an FIR SED with a similar MBB (i.e. with no assumption of optical thinness) to our model (Greve et al. 2012), this discrepancy in T_dust is due to their assumption of β ≡ 2.0.

In order to test whether these previously determined values of β and T_d could agree with our results, we fixed each parameter to their literature value, and fit a model using only two or one free variables. In the case where β was fixed to the previously determined value of 2, a nearly identical fit was returned, although with a lower dust mass and temperature, and a worse goodness of fit (i.e. a lower Bayesian evidence). If the dust temperature is instead fixed to 52 K, then we are unable to fit a reliable model to the data, as all models greatly underpredict the high-frequency flux density.

We may also examine the results of erroneously assuming optically thin dust emission. In this case, we find a significantly larger dust mass (⁠|$\rm \mathit{ M}_{dust}=9.36\pm 0.01\, M_{\odot }$|⁠), lower dust temperature (⁠|$\rm \mathit{ T}_{dust}=48.3\pm 0.7$| K), and lower dust emissivity constant (β = 1.44 ± 0.02), as well as a comparable FIR luminosity (L|$_{\rm FIR}=(2.3\pm 0.2)\times 10^{13}\,$| L_⊙) and SFR (4000 ± 400 M_⊙ yr⁻¹). This low dust temperature is similar to previously determined values, but the best-fitting value of β is greatly discrepant. This fit underpredicts the high-frequency data (⁠|$\lambda _{\rm obs}=160\!-\!500\, \mu$|m), suggesting a need for an additional MIR power law (e.g. Casey 2012). However, we find that our source is not optically thin over the examined frequency range (see Section 4.1.4), so this fit is non-physical.

4.1.4 Discussion of uncertainties

Here, we discuss possibly detrimental assumptions in our model, as well as the ways in which we attempt to minimize this uncertainty.

Modified Blackbody: The primary source of unaccounted uncertainty is the use of a single-temperature MBB model. Because of this, our ‘T_dust’ is a single luminosity-weighted dust temperature, rather than a mass-weighted dust temperature (e.g. Scoville et al. 2016) or a distribution of temperatures (e.g. Zhang et al. 2018). This model also assumes no contributions from separate MIR (e.g. Casey 2012), radio (e.g. Yun & Carilli 2002), or higher frequency (e.g. Leitherer et al. 1999) components. In particular, some studies have found that the addition of an MIR power law to an MBB results in a better fit to the FIR SEDs of some galaxies (e.g. Ibar et al. 2015; Faisst et al. 2020; Reuter et al. 2020).

Even though the MBB function is indeed simple, it has been found to fit dust SEDs very well (e.g. Bianchi 2013; Jiménez-Andrade et al. 2018; Carniani et al. 2019; Crocker et al. 2019; Lamperti et al. 2019). The frequency domain of our model is also controlled to ensure that the contributions from higher and lower frequency components are negligible. If an MIR component was required, then the fit shown in Fig. 3 would show a deviation at the high-frequency edge of the dust peak (see fig. 1D of Casey 2012). Since this frequency range is well fitted, we acknowledge that such a component may be present, but minimal.

Dust Absorption Coefficient: While most of the values in our model are allowed to vary, we assume a fixed value for the dust absorption coefficient κ_o. This practice has been widely adopted for radio-submm SED modelling (e.g. Wang et al. 2013; Wagg et al. 2014), and we adopt a relatively recent dust absorption coefficient (Beelen et al. 2006, rather than Hildebrand 1983) that falls into the range of values found by observational studies (Alton et al. 2004).

Source Radius: For the radius of SPT 0346−52, we average two magnification-corrected effective radii for continuum data taken at ν_rest ∼ 666 and 1000 GHz (0.79 ± 0.02  and 0.73 ± 0.03 kpc, respectively; Apostolovski et al. 2019), resulting in a well-constrained value of r_eff = 0.76 ± 0.02 kpc. These continuum observations were taken within the frequency range of our fit, supporting their use. While our model does not take the uncertainty of R_eff into account, we find that perturbing this value by 1σ results in minimal (i.e. <3σ) variations in the best-fitting values of each parameter.

Optical Thickness: As discussed in Section 4.1.1, we make no assumption on the optical thickness of our source, and therefore use a generalized MBB (similar to e.g. Leech, Metcalfe & Altieri 2001; Cortzen et al. 2020), rather than one that assumes τ_ν < <1 (e.g. Carniani et al. 2019; Lamperti et al. 2019; Valentino et al. 2020). To test whether this source is optically thin, we apply equation (2) to the data in Table 1, using our best-fitting dust mass and emissivity index. We find that τ_ν ranges from the marginally optically thin regime (0.22 at λ_obs = 3.0 mm) to optically thick (43.99 at λ_obs = 160 μm), with τ_ν ∼ 1 at λ_obs ∼ 200 μm. This a posteriori calculation verifies that the optically thin assumption cannot be applied to this source.

Uniformity: Our model assumes a circular source with a constant radius, luminosity-weighted dust temperature, lensing magnification, and set of dust emission properties (i.e. emission and absorption) across all FIR frequencies. Of course, SPT 0346−52 is not perfectly circular (e.g. Litke et al. 2019) and the magnification factor and intrinsic size of its dust continuum emission vary slightly with frequency (Apostolovski et al. 2019). Specifically, the dust continuum magnification factor has been found to decrease with increasing wavelength: μ_870μm = 5.6 ± 0.1, μ_2mm = 5.04 ± 0.09, and μ_3mm = 4.63 ± 0.03 (Spilker et al. 2016; Apostolovski et al. 2019). Therefore, our global assumption of μ = 5.6 may slightly skew the intrinsic SED. The details of dust heating and emission are also quite complex (e.g. Draine et al. 2007). Therefore, we note that these assumptions have a detrimental effect on the physicality of our best-fitting values.

Fitting Method: We note that the uncertainties in the best-fitting parameter values presented in Table 2 are simply the standard deviations of the posterior deviations for each parameter as found by multinest, and thus do not account for the additional sources of uncertainty discussed here. However, the use of multinest allows us to robustly explore the parameter space and to inspect the fits for model degeneracies through the construction of covariance plots.

4.2 Molecular hydrogen mass

4.2.1 HD-based estimate

Using our estimates of the integrated intensity of HD(J = 1–0) emission from this source, we may place limits on the amount of HD and H₂ present in SPT 0346−52. We begin by using our HD(1–0) flux density and equation (2) of Bergin et al. (2013):

$$\begin{eqnarray*} F_L=\int S_{\nu }\mathrm{ d}\nu =\frac{N_{\mathrm{ HD}}A_{10}h\nu f_u}{4\pi D_\mathrm{ L}^2}, \end{eqnarray*}$$

(6)

where N_HD is the number of HD molecules, A₁₀ is the Einstein coefficient for spontaneous emission (5.44 × 10⁻⁸ s⁻¹; Müller et al. 2005), h is Planck’s constant, ν is the rest frequency of the transition (2.674986 × 10¹² Hz), D_L is the luminosity distance (1.70 × 10²⁷ m), and f_u is the fraction of HD that is in J = 1:

$$\begin{eqnarray*} f_u=\frac{3\mathrm{ e}^{-128.5/T}}{Q(T)}, \end{eqnarray*}$$

(7)

where T is our temperature in Kelvin and Q(T) is the partition function.

Next, we assume that (HD/H₂)/(D/H) = 1, based on observations of z ∼ 2.5 galaxies with high metallicity and H i column density (Liszt 2015). We also adopt the primordial value of D/H = (2.6 ± 0.1) × 10⁻⁵ (Coc et al. 2004; Planck Collaboration et al. 2016; Novosyadlyj, Sergijenko & Shulga 2017), as the Universe is only 1 Gyr old at z = 5.656, so the deuterium fraction is still nearly primordial (van de Voort et al. 2018). This results in a simple equation for the molecular gas mass:

$$\begin{eqnarray*} M_{\mathrm{ H}_2}=\frac{4\pi m_{\mathrm{ H}_2} D_\mathrm{ L}^2 F_\mathrm{ L}}{A_{10}h \nu f_u D/H}, \end{eqnarray*}$$

(8)

where |$m_{\mathrm{ H}_2}$| is the mass of a single H₂ molecule. Note that while Bergin et al. (2013) consider the total gas mass, we do not consider helium and heavy elements and focus on H₂, the main gas component. We assume a gas temperature of T ∼ 60 K, based on Large Velocity Gradient modelling of CO(8–7) and CO(2–1) data (Dong et al. 2019).

Using our 2σ upper limit on the integrated flux of HD(1–0) (⁠|$\rm \mathit{ F}_L\lt 3.8\times 10^{-20}\, W\, m^{-2}$|⁠), this results in a magnification-corrected upper limit on the molecular gas mass of |$\rm \mathit{ M}_{H_2}\lt (6.4\times 10^{11})(5.6/\mu)$| M_⊙, where μ is the magnification factor (e.g. Spilker et al. 2015). This is in agreement with the CO(2–1) based gas mass of Aravena et al. (2016): |$\rm \mathit{ M}_{H_2}=(8.2\pm 0.6)\times 10^{10}(5.6/\mu)(\alpha _{CO}/0.8)$| M_⊙.

This gas mass estimate may be used to place a constraint on the CO luminosity–H₂ mass conversion ratio |$\alpha _{\mathrm{ CO}}=M_{\mathrm{ H}_2}/L^{\prime }_{\mathrm{ CO}}$|⁠. We convert the CO(2–1) integrated flux density of Aravena et al. (2016) to a CO(1–0) luminosity using |$L^{\prime }_{\mathrm{ CO}(2-1)}/L^{\prime }_{\mathrm{ CO}(1-0)}=0.85$| (appropriate for starburst galaxies) from table 2 of Carilli & Walter (2013), yielding an observed μL’ _CO(1–0) = (6.1 ± 0.4) × 10¹¹  K km s ⁻¹ pc², which we then correct for lensing (μ = 5.6 ± 0.1; Spilker et al. 2016). This results in a limit of α_CO < 5.8.

The above values assume that the kinetic temperature of HD is equal to that of CO, as derived by Dong et al. (2019). As an alternative estimate of the gas temperature, we briefly consider the possibility that the gas and dust are thermally coupled due to high densities (n > 10^4.5 cm⁻³; Goldsmith 2001). With this coupling, we may consider the dust temperature derived through our SED fit (78.6 ± 0.05 K), resulting in a slightly less conservative upper limit on the gas mass: |$\rm \mathit{ M}_{H_2}\lt 4.4\times 10^{11}\, M_{\odot }$|⁠. On the other hand, if the suggested mass-weighted dust temperature of Scoville et al. (2016) is assumed (25 K), then we find a more conservative gas mass limit of |$\rm \mathit{ M}_{H_2}\lt 9.9\times 10^{12}\, M_{\odot }$|⁠.

4.2.2 Continuum-based estimates

We may also use our derived dust mass (10^9.38 M_⊙) and a standard gas-to-dust ratio (GDR) of |$\rm \delta _{GDR}=\mathit{ M}_{H_2}/\mathit{ M}_{dust}=100\times 10^{Z-Z_{\odot }}$| (Draine et al. 2007) to determine a separate estimate on |$\rm M_{H_2}$|⁠. Since no metallicity information is available, we assume subsolar metallicity for this early galaxy, resulting in M|$_{\mathrm{ H}_2,\mathrm{ GDR}}\lt 2.4\times 10^{11}$| M_⊙.

Alternatively, equation (A14) of Scoville et al. (2016) may be applied to our four continuum flux density values that satisfy their λ_rest > 250 μm criterion (λ_obs > 1664 μm). Using our derived dust temperature, this yields a magnification-corrected average value of |$M_{\mathrm{ H}_2}=(5.1\pm 0.3)\times 10^{11}$| M_⊙. It should be noted that this equation was derived using a calibration sample of z ∼ 0 − 3 galaxies with high stellar masses, and thus high metallicities. Its applicability to higher-redshift objects with possibly lower metallicities like SPT 0346−52 is uncertain.

These two continuum-based gas mass estimates may also be used to place constraints on the CO luminosity–H₂ mass conversion ratio |$\alpha _{\mathrm{ CO}}=M_{\mathrm{ H}_2}/L^{\prime }_{\mathrm{ CO}}$|⁠. Both the Scoville et al. (2016) and GDR approaches result in values between the Milky Way (4.6) and starburst (0.8) limits (e.g. Daddi et al. 2010). These gas mass estimates and their resulting |$\rm \alpha _{CO}$| values are listed in Table 3.

Table 3.

Derived molecular gas masses using different estimators, and the corresponding α_CO values based on the L’_CO(2–1) value of Aravena et al. (2016). For reference, the CO(2–1)-based gas mass of Aravena et al. (2016) (assuming α_CO = 0.8) is also included.

Method	M|$_{\rm H_2}$| [ M⊙]	αCO
HD	<6.4 × 1011	<5.8
GDR	<2.4 × 1011	<2.2
Scoville + 16	(5.1 ± 0.3) × 1011	4.0 ± 0.4
Aravena + 16	(8.2 ± 0.6) × 1010	⋅⋅⋅

5 CONCLUSIONS

In this work, we have presented the first upper limit on the luminosity of HD(J = 1–0) in SPT 0346−52, fitted an MBB model to an FIR SED of the source, and derived several estimates of the mass of molecular gas in the system. All of these results confirm the extreme starburst nature of this galaxy, but there is variation in some derived values.

No significant HD(J = 1–0) emission is detected, implying a molecular gas mass of |$\rm \mathit{ M}_{H_2}\lt 6.4\times 10^{11}\, M_{\odot }$|⁠. This estimate is in agreement with a previous estimate of |$\rm \mathit{ M}_{H_2}=(8.2\pm 0.6)\times 10^{10}\, M_{\odot }$| based on CO(2–1) emission (Aravena et al. 2016).

Using archival ALMA observations, we are able to fill in the red side of the dust SED, resulting in a better constraint in the dust emissivity spectral index β and new estimates on the dust mass and dust temperature. Our best-fitting luminosity-weighted dust temperature (T = 78.6 ± 0.5 K) is greater than both the extreme mass-weighted temperature case of 25 K (Scoville et al. 2016) and the previous estimate of 53 K (Weiß et al. 2013). This discrepancy is possibly due to the extreme optical depth of our source (i.e. τ_o reaches 43.99 at |$\lambda _{\mathrm{ obs}}=160\, \mu$|m), or our different method of deriving the dust emissivity spectral index β. Our derived value of β = 1.81 ± 0.03 is reasonable considering the standard range of β = 1.5–2.0 for high redshift, dusty galaxies (e.g. Chapman et al. 2005; Casey et al. 2011; Chapin et al. 2011). We present a thorough discussion of the assumptions of our models, and state possible sources of unstated uncertainty in our best-fitting models.

The molecular mass of the entire system was also estimated with two continuum-based methods. First, the Scoville et al. (2016) S(ν_rest < 1.2 THz) estimator was applied to four of our SED values, resulting in an average value of |$\rm \mathit{ M}_{H_2}=(5.1\pm 0.3)\times 10^{11}$| M_⊙. In addition, our dust mass was converted to a gas mass by assuming subsolar metallicity and a possible dust-to-gas ratio, yielding an upper limit of |$\rm \mathit{ M}_{H_2}\lt 2.4\times 10^{11}$| M_⊙. The difference between these values may suggest that the dust temperature is underestimated, the galaxy has a high metallicity, or simply that this source does not follow the Scoville et al. (2016) relation, which was created using local starbursts.

These three gas mass estimates may be used to place constraints on α_CO, resulting in <5.8, <2.2, and 4.0 ± 0.4 for the HD, GDR, and Scoville et al. (2016) approaches, respectively. While they are disparate, they either agree with or are between the MW-like α_CO ∼ 4.6 (e.g. Daddi et al. 2010) and the starburst α_CO ∼ 0.8, which was assumed by Aravena et al. (2016). Due to the large star formation rate of SPT 0346−52 (∼4000 M_⊙ yr⁻¹), these high conversion factors may be surprising. However, a recent investigation of the submillimeter galaxy AzTEC/C159 by Jiménez-Andrade et al. (2018) found α_CO ∼ 4, despite the high star formation rate of the source (∼700 M_⊙ yr⁻¹). One of the interpretations of these results was that AzTEC/C159 is in the early stages of a starburst. It is thus possible that SPT 0346−52 is undergoing a massive starburst driven by the presence of a large amount of pristine gas.

Due to the lensed nature of the source, we are unable to present any exact SFR or molecular mass surface densities. However, as a test of our values, we consider a size of ∼0.76 kpc (Apostolovski et al. 2019) and a maximal limit of |$\rm A_{SFR}=A_{H_2}$|⁠. This would suggest a Σ_SFR ∼ 2 × 10³ M_⊙ yr⁻¹ kpc⁻² and Σ_H2 ∼ 5.5 × 10⁴ M_⊙ pc⁻². This places SPT 0346−52 squarely in the starburst region of a Kennicutt–Schmidt diagram (e.g. Kennicutt 1998b). These results confirm that SPT 0346−52 is a highly starbursting, dusty, extremely luminous galaxy.

We note that our observations were designed to maximize sensitivity to compact emission, and so the non-detection of HD(1–0) may be partly due to the extended nature of the emission. In addition, the expected line is coincident with a series of atmospheric transmission lines, resulting in a non-constant noise level. Future observations of HD(1–0) may benefit from targeting sources at higher-redshift (i.e. z > 7.5), where the line will fall into a frequency range of higher atmospheric transmission (i.e. ALMA band 7 or below). With the rising number of detections of CO emission at z > 5 (e.g. D’Odorico et al. 2018; Pavesi et al. 2018, 2019; Riechers et al. 2020), it will be of interest to have an independent tracer of molecular gas mass in the early Universe.

ACKNOWLEDGEMENTS

This paper makes use of the following ALMA data: 2015.1.01580.S, 2013.1.01231.S, 2016.1.00654.S, and 2016.1.01313.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. GCJ and RM acknowledge ERC Advanced Grant 695671 ‘QUENCH’ and support by the Science and Technology Facilities Council (STFC). SC acknowledges support from the ERC Advanced Grant INTERSTELLAR H2020/740120. We thank the anonymous referee for constructive feedback that strengthened this manuscript.

DATA AVAILABILITY

The data analysed in this work are available from the ALMA data archive (https://almascience.nrao.edu/asax/) under project code 2016.1.01313.S (Band 8).

Footnotes

1

https://almascience.eso.org/about-alma/atmosphere-model

APPENDIX A: TENTATIVE DETECTIONS

While there are no significant detections of HD(J = 1–0)112.2 μm line emission in these data, there are four emission features that could be interpreted as tentative detections (see Section 3.2). Here, we identify the velocity width of each feature, use these velocity ranges to create moment zero maps and integrated spectra, and discuss why they are disregarded as spurious noise peaks. In order to avoid including noise, spectra are extracted from a combination of two spatial masks: >2σ in FIR continuum emission and >3σ in the moment zero map created using the relevant channels. Velocities are given with respect to the expected redshifted frequency of HD(1–0) at z = 5.656, or 401.890 GHz.

The first tentative line feature is detected through an exploration of the continuum-subtracted data cube. It is primarily detected at the northern image (‘Tentative 1’, −192 to 30 km s⁻¹, 402.148 to 401.851 GHz). The moment zero map of this feature shows >3σ emission near the north continuum peak (red contours in top left panel of Fig. A1), and the emission line itself is broad and has a small velocity offset from SPT 0346−52 (Fig. A2), so this feature is marginally believable. However, multiple other 3σ noise peaks are apparent in the collapsed image, so this feature is not confidently detected.

Figure A1.

Moment zero maps of four tentative line detections (coloured contours) compared with FIR continuum (background grey-scale). Contours are shown at −3, 3, 4 × σ, where 1σ = 25.0, 18.0, 43.0, and 9.6 mJy beam⁻¹ km s⁻¹ for Tentative 1 through 4, respectively.

Figure A2.

Spectra of each tentative line feature in Fig. A1. Each is extracted using combined masks of >2σ in FIR continuum emission and >3σ in the moment zero map. The possible line feature is denoted by the yellow shaded region, while the expected redshift of HD(1–0) is shown by a vertical dashed line. No lensing correction has been applied.

This exploration of the data cube resulted in an additional feature detected near the southeastern image (‘Tentative 2’, 193 to 368 km s⁻¹, 401.632 to 401.398 GHz). The collapsed image of these channels shows >4σ emission near the southeastern image and >3σ emission near the northern image (blue contours of Fig. A1). However, the spectrum shows a significant velocity offset from SPT 0346−52(∼300 km s⁻¹) and a double-peaked profile, suggesting noise domination. At first glance, this line shape may be interpreted as evidence for rotation. However, previous studies of line emission have revealed that the FWHM of this source is ∼600 km s⁻¹ (e.g. Aravena et al. 2016), while this feature is <200 km s⁻¹ wide. Since it has a high-velocity offset, <5σ significant peaks in the collapsed image, and a non-Gaussian profile, we conclude that this feature is noise.

One possible feature is present in a wide-velocity (⁠|$\rm -350\lt v\lt 350$| km s⁻¹) moment zero map (‘Tentative 3’). It is exactly coincident with the northern image of the FIR continuum, and is detected at ∼3σ, suggesting a real detection. However, it is comparable to noise peaks in the field of view, and a spectrum extracted from its 3σ contour shows a broad velocity range of positive, low-level emission. Since the amplitude of this feature is nearly constant across this range (i.e. no central peak) and no emission is detected from the other lensed images, we believe this increased emission to be an artefact of continuum subtraction due to the very high significance of the northern continuum peak (∼185σ; Section 3.1).

When it is assumed that the spatial distributions of HD(1–0) and FIR continuum emission are equally extended, we find a single spectral feature (‘Tentative 4’). This emission is quite narrow (1247 to 1294 km s⁻¹, 400.219 to 400.156 GHz), is close to an atmospheric transmission line, and is at a large velocity offset from the expected redshift of HD(1–0). In addition, the moment zero map shows that there is only one 3σ peak present in the southwestern arc. Therefore, this feature is likely only noise.