COLDz: A High Space Density of Massive Dusty Starburst Galaxies ∼1 Billion Years after the Big Bang

CO(J = 2 → 1) line emission (ν_rest = 230.5380 GHz) toward HDF 850.1, AzTEC-3, and GN10 was detected in molecular line scans with the VLA within the COLDz survey data (project IDs: 13A-398 and 14A-214; PI: Riechers). A detailed description of the data is given by P18. In brief, COLDz targeted two regions in the COSMOS and GOODS-North survey fields at 35 and 34 GHz (corresponding to ∼8.7 mm), covering areas of ∼9 and 51 arcmin² in 7- and 57-point mosaics with the VLA, respectively. Observations were carried out under good Ka band weather conditions for a total of 324 hr between 2013 January 26 and 2015 December 18 in the D and DnC array configurations, as well as reconfigurations between the C, DnC, and D arrays. The correlator was set up with two intermediate frequency bands (IFs) of 4 GHz bandwidth (dual polarization) each in 3-bit mode, centered at the frequencies indicated above. Gaps between individual 128 MHz sub-bands were mitigated through frequency switching. The radio quasars J1041+0610 and J1302+5748 were observed for complex gain calibration and regular pointing corrections in the COSMOS and GOODS-North fields, respectively. The absolute flux scale was derived based on observations of 3C 286.

Data reduction and imaging were performed with the casa 4.1 package,²⁶ using the data pipeline version 1.2.0. Imaging the data with natural baseline weighting yields typical clean beam sizes of 25, with variations between individual pointings and across the large bandwidth. Typical rms noise levels are 60 and 100–200 μJy beam⁻¹ per 4 MHz (∼35 km s⁻¹) binned channel in the COSMOS and GOODS-North fields, respectively. At the positions of GN10 and AzTEC-3 (for which maps based on these data are shown in the following sections), the rms noise is 51 and 18 μJy beam⁻¹ per 76 and 60 MHz (∼623 and 491 km s⁻¹) binned channel at beam sizes of 195 × 167 and 246 × 226, respectively.

We observed CO(J = 1 → 0) line emission toward GN10 using the VLA (project ID: 16A-015; PI: Riechers). Observations were carried out under good Ku band weather conditions for a total of 11 hr during five tracks in the C array between 2016 February 2 and March 6. Two IFs with 1 GHz bandwidth (dual polarization) each in 8-bit mode were centered at 13.977 and 17.837 GHz (corresponding to 2.1 and 1.7 cm, respectively) to cover the redshifted HCN, HCO⁺, and HNC(J = 1 → 0) and CO(J = 1 → 0) lines in GN10 (ν_rest = 88.6318, 89.1885, 90.6636, and 115.2712 GHz). Observations were carried out in short cycles, spending between ∼330 and ∼470 s on source, bracketed by scans spending ∼75 s on the gain calibrator J1302+5748. Pointing was performed on the gain calibrator approximately once per hour. The absolute flux scale was derived based on observations of 3C 286.

Data reduction and imaging were performed with the casa 5.4.2 package. Imaging the data with natural baseline weighting yields a clean beam size of 175 × 128 at an rms noise level of 21 μJy beam⁻¹ over 40 MHz (656 km s⁻¹) at the CO(J = 1 → 0) line frequency. The rms noise near the HCN(J = 1 → 0) frequency is ∼27 μJy beam⁻¹ over 4 MHz (85 km s⁻¹). Imaging the data over the entire line-free bandwidth of 2.012 GHz yields an rms noise level of 1.49 μJy beam⁻¹ at a beam size of 210 × 154.

We observed continuum emission at 22.8649 GHz (K band) and 45.6851 GHz (Q band) toward GN10 (corresponding to 1.3 cm and 6.6 mm, respectively), using the VLA (project ID: AR693; PI: Riechers).²⁷ Observations were carried out under good K and Q band observing conditions for a total of 44 hr between 2009 July 19 and 2010 January 5. K-band observations were conducted for four tracks in the C array, totaling 28 hr, and Q-band observations were conducted for two tracks in the D array, totaling 16 hr. All observations used the previous generation correlator, covering two IFs of 50 MHz bandwidth (dual polarization) each at the tuning frequency and 300 MHz (K band) or 50 MHz (Q band) above, respectively, in quasi-continuum mode. Observations in the C (D) array were carried out in short cycles, spending 150 s (200 s) on source, bracketed by scans spending 60 s on the gain calibrator J13028+57486. Pointing was performed on the gain calibrator approximately once per hour. The absolute flux scale was derived based on observations of 3C 286.

Data reduction and imaging were performed with the aips package.²⁸ Imaging the data with natural baseline weighting yields clean beam sizes of 115 × 101 and 182 × 168 at rms noise levels of 44 and 58 μJy beam⁻¹ over 100 MHz in the K and Q bands, respectively.

We observed CO(J = 2 → 1) line emission toward HDF 850.1 at higher spatial resolution using the VLA (project ID: 16A-014; PI: Riechers). Observations were carried out under good Ka band weather conditions for a total of 8.8 hr during four tracks in the C array between 2016 February 3 and 20. Two IFs with 4 GHz bandwidth (dual polarization) each in 3-bit mode were centered at 34 GHz (corresponding to 8.8 mm) to cover the same frequency range as the D array observations of the main survey. Gaps between sub-bands were mitigated through frequency switching, using the same two setups with a relative shift of 16 MHz, as in the D array. Observations were carried out in short cycles, spending ∼300 s on source, bracketed by scans spending ∼100 s on the gain calibrator J1302+5748. Pointing was performed on the gain calibrator approximately once per hour. The absolute flux scale was derived based on observations of 3C 286.

Data reduction and imaging were performed with the casa 5.4.2 package. Imaging the data with natural baseline weighting yields a clean beam size of 071 × 068 at an rms noise level of 32.4 μJy beam⁻¹ over 66 MHz (530 km s⁻¹) at the CO(J = 2 → 1) line frequency.

We observed CO(J = 6 → 5) line emission (ν_rest = 691.4731 GHz) toward GN10 using NOEMA (project ID: X–5; PI: Riechers). Observations were carried out under good 3 mm weather conditions for four tracks in the A configuration between 2014 February 18 and 24, using six antennas (baseline range: 67–760 m). This yielded a total time of 13.8 hr (16500 visibilities) on source. Receivers were tuned to 109.7037 GHz (corresponding to 2.7 mm). The correlator was set up with a bandwidth of 3.6 GHz (dual polarization).

Data reduction and imaging were performed with the gildas package. Imaging the data with natural or uniform baseline weighting yields clean beam sizes of 082 × 071 or 063 × 059 at rms noise levels of 73 or 90 μJy beam⁻¹ over 700 MHz (1912 km s⁻¹), respectively. Imaging the data with natural weighting over the entire line-free bandwidth of 2.9 GHz yields an rms noise level of 31.6 μJy beam⁻¹.

We observed [C ii] 158 μm line emission (ν_rest = 1900.5369 GHz) toward GN10 using NOEMA (project ID: W14FH; PI: Riechers). Observations were carried out under good 0.9 mm weather conditions for one track in the C configuration on 2015 April 15, using six antennas (baseline range: 21–172 m). This yielded a total time of 1.9 hr (2249 visibilities) on source. Receivers were tuned to 301.524 GHz (corresponding to 1.0 mm). The correlator was set up with a bandwidth of 3.6 GHz (dual polarization).

Data reduction and imaging were performed with the gildas package. Imaging the data with natural or uniform baseline weighting yields clean beam sizes of 101 × 084 or 081 × 076 at rms noise levels of 0.62 or 0.71 mJy beam⁻¹ over 800 MHz (795 km s⁻¹), respectively. Imaging the data with natural or uniform weighting over the entire line-free bandwidth of 2.31 GHz yields rms noise levels of 324 or 365 μJy beam⁻¹.

We observed continuum emission at 137.057 GHz (2.2 mm)²⁹ and 250.5 GHz (1.2 mm) toward GN10 using NOEMA (project IDs: T047 and T0B7; PI: Riechers). Observations were carried out under good weather conditions for three tracks between 2009 June 4 and September 21 in the D configuration with five antennas (baseline range: 19–94 m) at 2 mm and for two tracks on 2011 January 23 and 24 in the A configuration with six antennas (baseline range: 51–665 m) at 1.2 mm. This yielded a total of 7.1 and 3.7 hr (17,040 and 8940 visibilities) six-antenna-equivalent on-source time at 2 and 1.2 mm, respectively. Observations at 2 mm were carried out with the previous generation correlator with a bandwidth of 1 GHz (dual polarization). Observations at 1.2 mm were carried out with a bandwidth of 3.6 GHz (dual polarization).

Data reduction and imaging were performed with the gildas package. Imaging the 2 mm data with natural baseline weighting yields a clean beam size of 37 × 32 at an rms noise level of 95 μJy beam⁻¹ over 1 GHz bandwidth. Imaging the 1.2 mm data with natural or uniform baseline weighting yields clean beam sizes of 045 × 038 or 038 × 033 at rms noise levels of 0.35 or 0.43 mJy beam⁻¹ over 3.6 GHz, respectively.

Serendipitous CO(J = 5 → 4) line emission (ν_rest =576.2679 GHz) was observed toward GN10 using NOEMA. These observations, taken from Daddi et al. (2009a), did not target GN10, which was offset by 97 or 197 from the phase center for different tracks, yielding primary beam attenuation factors of 1.08 or 1.30, respectively. Observations were carried out under good 3 mm weather conditions for four tracks in the B, C, and D configurations between 2008 May 4 and 2009 January 5, using six antennas (baseline range: 15–411 m). This yielded a total time of 14.6 hr (25186 visibilities) on source. Receivers were tuned to 91.375 GHz (corresponding to 3.3 mm). The correlator was set up with a bandwidth of 1 GHz (dual polarization).

We adopted the data reduction performed by Daddi et al. (2009a) but re-imaged the data with the gildas package. Imaging the data with natural baseline weighting yields a clean beam size of 264 × 190 at an rms noise level of 66.6 μJy beam⁻¹ over 365.753 MHz (1200 km s⁻¹) at the phase center (96 μJy beam⁻¹ at the position of GN10). Imaging the data over the line-free bandwidth yields an rms noise level of 55 μJy beam⁻¹.

We observed CO(J = 5 → 4) line emission (ν_rest = 576.2679 GHz) toward AzTEC-3 using NOEMA (project ID: U0D0; PI: Riechers). Observations were carried out under good 3 mm weather conditions for two tracks in the A configuration between 2011 January 19 and February 4, using six antennas (baseline range: 100–760 m). We also used previous observations (project ID: T–F; PI: Riechers) carried out for one track in the C configuration on 2010 April 1, using six antennas (baseline range: 15–176 m; see Riechers et al. 2010 for additional details). This yielded a total time of 8.8 hr (10560 visibilities; 5340 in A configuration) on source. Receivers were tuned to 91.558 GHz (corresponding to 3.3 mm). The correlator was set up with a bandwidth of 3.6 GHz (dual polarization).

Data reduction and imaging were performed with the gildas package. Imaging the combined data with natural baseline weighting yields a clean beam size of 221 × 143 at an rms noise level of 64 μJy beam⁻¹ over 280 MHz (917 km s⁻¹). Imaging the A configuration data only with uniform baseline weighting yields a clean beam size of 139 × 085 at an rms noise level of 96 μJy beam⁻¹ over 260 MHz (852 km s⁻¹). Imaging the A (A+C) configuration data with natural weighting over the entire line-free bandwidth of 3.34 GHz yields an rms noise level of 24.8 (18.8) μJy beam⁻¹.

We observed continuum emission at 250.0 GHz (1.2 mm) toward AzTEC-3 using NOEMA (project ID: U0D0; PI: Riechers). Observations were carried out under good weather conditions for three tracks between 2011 January 25 and February 03 in the A configuration with six antennas (baseline range: 100–760 m) at 1.2 mm. This yielded a total of 9.3 hr (11160 visibilities) on source. Observations were carried out with a bandwidth of 3.6 GHz (dual polarization).

Data reduction and imaging were performed with the gildas package. Imaging the data with natural or uniform baseline weighting yields clean beam sizes of 062 × 025 or 053 × 024 at rms noise levels of 162 or 187 μJy beam⁻¹ over 3.6 GHz, respectively.

We detect strong continuum emission toward GN10 at 1.2 and 1.0 mm and weak emission between 2.2 mm and 1.9 cm (see Figure 4 and Table 2). The flux keeps decreasing between 0.9 and 1.9 cm and is >4–8 times lower than at 21 cm. This suggests that the emission detected up to 0.9 cm (i.e., rest frame 1.4 mm) likely still corresponds to thermal emission, but nonthermal emission may start to significantly contribute at 1.9 cm (i.e., rest frame 3.0 mm; see Figure 8). The continuum emission is spatially resolved along the major axis at 1.2 mm by our observations with a synthesized beam size of ∼035. By fitting two-dimensional Gaussian profiles to the emission in the visibility plane, we find a size of (025 ± 007) × (010 ± 011), which corresponds to (1.6 ± 0.4) × (0.6 ± 0.6) kpc². A circular Gaussian fit provides a full width at half power (FWHP) diameter of 018 ± 005, or (1.1 ± 0.3) kpc. Due to the agreement of the 1.2 mm flux with a previous measurement at lower spatial resolution at a close wavelength (Dannerbauer et al. 2008; see Table 2), and given the baseline coverage down to ∼50 m, it appears unlikely that the 1.2 mm size measurement is biased toward low values due to missing emission. Fits to the lower-resolution (∼075 beam size) data at 1.0 mm however suggest a size of (058 ± 012) × (050 ± 010), corresponding to (3.6 ± 0.7) × (3.1 ± 0.6) kpc², or a circular FWHP diameter of 053 ± 008 (3.3 ± 0.5 kpc²). The uncertainties for this measurement may be limited by interferometric seeing due to phase noise (which is not factored into the fitting errors), such that we treat the 1.0 mm size measurement as an upper limit only in the following sections. We however note that, in principle, the dust emission at shorter wavelengths could appear more extended due to an increasing dust optical depth as well, as discussed further in Section 4. Also, the source shape could significantly deviate from a Gaussian shape (such as a higher-index Sérsic profile; see, e.g., discussion by Hodge et al. 2016), such that more complex fitting procedures could yield different findings.

Zoom In Zoom Out Reset image size

Despite its strong dust continuum emission at sub/millimeter wavelengths, GN10 remains undetected up to an observed-frame wavelength of 2.2 μm, rendering it a "K-band dropout" (also see discussion by Wang et al. 2009; Daddi et al. 2009a, and references therein). Even sensitive space-based imaging up to 1.6 μm with the WFC3 camera on the Hubble Space Telescope (HST) from the CANDELS survey (Grogin et al. 2011) yields no hint of emission due to dust obscuration (Figure 4), but mid-infrared continuum emission is detected with Spitzer/IRAC longward of 3.6 μm (corresponding to ∼5700 Å in the rest frame; Dickinson & GOODS Team 2004; Giavalisco et al. 2004) at the position of the millimeter-wave dust continuum emission (Figure 4). Thus, the stellar light in the "optically dark" galaxy GN10 is not entirely obscured by dust.³¹

We successfully detect CO(J = 1 → 0), CO(J = 2 → 1), CO(J = 5 → 4), CO(J = 6 → 5), and [C ii] emission toward GN10 at 3.3σ, 8.6σ, 6.9σ, 7.3σ, and 18.1σ peak significance, respectively (Figures 2 and 3 and Table 1; see Appendix A for further details). In combination, these lines provide an unambiguous redshift identification. We extract the parameters of all emission lines from Gaussian fitting to the line profiles. All CO lines are fit with single Gaussian functions (see Appendix A for further details). We find line peak fluxes of S_line = (74 ± 13), (544 ± 63), (1046 ± 205), and (719 ± 144) μJy at FWHM line widths of dv_FWHM = (687 ± 144), (512 ± 72), (772 ± 220), and (681 ±173) km s⁻¹ for the CO J = 1 → 0, 2 → 1, 5 → 4, and 6 → 5 lines, respectively. The [C ii] line is fit with two Gaussian components, yielding S_line = (24.7 ± 1.6) and (6.0 ± 4.2) mJy and dv_FWHM = (617 ± 67) and (227 ± 243) km s⁻¹ for the red- and blueshifted components, respectively. The parameters of the secondary component thus are only poorly constrained, and we report [C ii] fluxes including and excluding this component in the following. It may correspond to a gas outflow, or a close-by minor companion galaxy to the DSFG.³² Considering only the main [C ii] component, the widths of all lines are consistent within the relative uncertainties.

Assuming the same widths as for the CO(J = 2 → 1) line, we find 3σ upper limits of <0.017 Jy km s⁻¹ for the HCN, HCO⁺, and HNC(J = 1 → 0) lines.³³ This implies HCN, HCO⁺, and HNC to CO line luminosity ratio limits of the order of <40%, which are only modestly constraining given the expectation of a few percent to ∼20% ratios for distant starburst galaxies (see, e.g., Riechers et al. 2007; Oteo et al. 2017b, and references therein).

The integrated line fluxes and line luminosities derived from these measurements are summarized in Table 1, together with those of AzTEC-3 and HDF 850.1, including values from the literature. The CO(J = 2 → 1) flux of GN10 reported here is somewhat lower than that found by P18 using a different extraction method. Here, we adopt our updated value for consistency. Given the more complex velocity profile of the [C ii] line in GN10, we adopt the CO(J = 2 → 1)-based redshift of z = 5.3031 ± 0.0006 in the following. This measurement agrees within 1σ with the CO(J = 1 → 0), CO(J = 5 → 4), and CO(J = 6 → 5)-based measurements.³⁴ The systemic velocity centroid of the [C ii] line is slightly blueshifted, but emission is detected across the same velocity range as in the CO J = 1 → 0 to 6 → 5 lines. The line peak offset thus is likely mainly due to internal variations in the [C ii]/CO ratio.

By fitting two-dimensional Gaussian profiles to the [C ii] emission, we find a size of (104 ± 030) × (019 ± 010). This corresponds to (6.5 ± 1.9) × (1.2 ± 0.6) kpc². Attempting to fit the CO(J = 6 → 5) emission observed at a comparable spatial resolution yields a size of 042 × < 025 (2.6 × < 1.6 kpc²), but the fit does not converge well due to the moderate signal-to-noise ratio of the detection. A circular fit is consistent with a point source within the uncertainties. This suggests that the [C ii] emission is resolved at least along the major axis and that it appears to be more spatially extended than the mid-J CO and dust emission, which are consistent with having a comparable spatial extent within the current uncertainties. The [C ii] emission shows a significant spatial velocity gradient across the line emission in the main component alone (see Figure 3).

We detect strong continuum emission toward AzTEC-3 at 1.2 mm and weak emission at 3.3 mm at the same peak position (Figure 5). We measure a continuum flux of (118 ± 25) μJy at 3.3 mm (i.e., rest frame 520 μm) from the A configuration data, which is consistent with the 3σ upper limit of <0.12 mJy obtained from the C configuration data (Riechers et al. 2010). The emission appears unresolved in the A configuration data. Combining both data sets, we find a final 3.3 mm flux of (126 ± 19) μJy. The continuum emission is spatially resolved along the major axis at 1.2 mm by our observations with a synthesized beam size of ∼025. By fitting two-dimensional Gaussian profiles in the image plane, we find a size of (045 ± 014) × ($0\buildrel{\prime\prime}\over{.} {05}_{-0\buildrel{\prime\prime}\over{.} 05}^{+0\buildrel{\prime\prime}\over{.} 21}$), which corresponds to (2.8 ± 0.9) × (${0.3}_{-0.3}^{+1.3}$) kpc², and we measure a 1.2 mm flux of (3.67 ± 0.56) mJy. This is consistent with the size of the 1.0 mm continuum emission of (040 ± 004) × ($0\buildrel{\prime\prime}\over{.} {17}_{-0\buildrel{\prime\prime}\over{.} 17}^{+0\buildrel{\prime\prime}\over{.} 08}$) and the dust spectral energy distribution shape found by Riechers et al. (2014a). A circular Gaussian fit in the visibility plane (which gives more weight to the longer baseline data) yields an FWHP diameter of 014 ± 004, or (0.9 ± 0.2) kpc, but a flux of only (2.75 ± 0.29) mJy. Taken at face value, this appears to suggest that ∼75% of the emission emerges from a compact region within the 2.5–3 kpc diameter dust reservoir.³⁵

Zoom In Zoom Out Reset image size

We detect CO(J = 5 → 4) emission toward AzTEC-3 at 15.0σ peak significance (Figure 5; CO J = 2 → 1 and 6 → 5 are also shown for reference). From a circular Gaussian fit in the visibility plane to the moment-0 map, we find a line flux of I_{CO(5 − 4)} = (0.97 ± 0.09) Jy km s⁻¹, which is consistent with a previous measurement by Riechers et al. (2010). From Gaussian fitting to the line profile (Figure 6; CO J = 2 → 1 and 6 → 5 and [C ii] are also shown for reference), we obtain a peak flux of S_{CO(5 − 4)} = (1.88 ± 0.14) mJy and an FWHM of dv_FWHM = (396 ± 37) km s⁻¹, which is consistent with the previously measured values and those found for the CO(J = 2 → 1) and [C ii] lines (Riechers et al. 2014a). The line may show an excess in its red wing that is not captured well by the Gaussian fit, but the significance of this feature is only moderate. A circular Gaussian fit in the visibility plane over the entire width of the emission of 917 km s⁻¹ (280 MHz) suggests an FWHP source diameter of 044 ± 017, which corresponds to (2.8 ± 1.1) kpc. This is comparable to the size of the 1.0 and 1.2 mm continuum emission. Fitting the source over 524 km s⁻¹ (160 MHz) to capture only the main component of the emission, we find 057 ± 015, which corresponds to (3.5 ± 0.9) kpc. This is comparable to the size of the [C ii] emission of (063 ± 009) × ($0\buildrel{\prime\prime}\over{.} {34}_{-0\buildrel{\prime\prime}\over{.} 15}^{+0\buildrel{\prime\prime}\over{.} 10}$) found by Riechers et al. (2014a) over a similar velocity range (466 km s⁻¹). There appears to be a small spatial offset between the peaks of the CO(J = 5 → 4) and [C ii] emission (which is consistent with the dust continuum peak; Figure 5). However, this offset becomes insignificant (≲1σ) when excluding the red CO line wing, i.e., when considering comparable velocity ranges, which results in a shift of the peak position by ∼01 (≲2σ significance shift) relative to the broader velocity range. If confirmed, the emission in the red line wing thus may correspond to a spatially offset kinematic component, but additional data are required to investigate this finding in more detail.

Zoom In Zoom Out Reset image size

To extract physical parameters from the spectral energy distribution (SED) of GN10, we have followed two approaches, as summarized in Figure 8 and Table 3. First, we have fit the far-infrared (FIR) peak of the SED with a modified blackbody (MBB) routine, where the MBB is joined to a smooth power law with slope α toward observed-frame mid-infrared wavelengths (e.g., Blain et al. 2003, and references therein). The dust temperature T_dust, dust emissivity parameter β_IR, and the wavelength λ₀, where the optical depth becomes unity, are used as fitting parameters. In addition, the flux f_158μm^rest at rest frame 158 μm is used as a normalization parameter. The Markov Chain Monte Carlo (MCMC) and nested sampling package emcee (Foreman-Mackey et al. 2013) is used to explore the posterior parameter distribution. By integrating the fitted functions, we obtain far-infrared (L_FIR) and total infrared (L_IR) luminosities, which we then use to determine dust-obscured star formation rates, SFR_IR, under the assumption of a Kennicutt (1998) conversion for a Chabrier (2003) stellar initial mass function. By adopting a mass absorption coefficient of κ_ν = 2.64 m² kg⁻¹ at 125 μm (Dunne et al. 2003), we also estimate a dust mass M_dust from the fits, finding a high value in excess of 10⁹ M_⊙ (see Table 3).

Zoom In Zoom Out Reset image size

Table 3.  GN10 MBB and cigale SED Modeling Parameters

Fit Parameter	Unit	Valuea
Tdust	K	${50.1}_{-9.1}^{+9.1}$
βIR		${2.98}_{-0.17}^{+0.18}$
${\lambda }_{0}^{\mathrm{rest}}$	μm	${168}_{-25}^{+22}$
α		${2.06}_{-0.11}^{+0.15}$
${f}_{158\mu {\rm{m}}}^{\mathrm{rest}}$ b	mJy	${8.3}_{-0.5}^{+0.5}$
Mdust	109 M⊙	${1.11}_{-0.27}^{+0.44}$
LFIRc	1013 L⊙	${0.58}_{-0.09}^{+0.11}$
LIRd	1013 L⊙	${1.03}_{-0.15}^{+0.19}$
SFRIR	M⊙ yr−1	${1030}_{-150}^{+190}$
M⋆e	1011 M⊙	1.19(±0.06)

Notes.

^aValues as stated are the 50th percentiles. The lower and upper error bars are stated as the 16th and 84th percentiles, respectively. ^bFit normalization parameter. ^cIntegrated between rest frame 42.5 and 122.5 μm. ^dIntegrated between rest frame 8 and 1000 μm. ^eParameter adopted from cigale. Quoted fitting uncertainties underestimate systematic effects and, thus, are not adopted to make any conclusions.

Download table as:  ASCIITypeset image

We also fit the full optical to radio wavelength photometry of GN10 using the Code Investigating GALaxy Emission (cigale; Noll et al. 2009; Serra et al. 2011), using a broad range of star formation histories and metallicities and a standard dust attenuation law (Calzetti 2001). We find parameters that are broadly consistent with the MBB fitting where applicable. Due to its high level of dust obscuration, GN10 remains undetected shortward of an observed-frame of 3.6 μm (rest frame ∼5700 Å), which leaves some parameters only poorly constrained. Thus, we only adopt the stellar mass M_⋆ from the cigale fit in the following sections. We find a high value in excess of 10¹¹ M_⊙, i.e., ∼100 times the dust mass, which we consider to be reliable to within a factor of a few (see Table 3).³⁶

From our MBB fits, we find that the dust turns optically thick at ∼170 μm (i.e., between an observed-frame wavelength of 1.0 and 1.2 mm; see Table 3), similar to the values found for other z > 4 DSFGs (see, e.g., Riechers et al. 2013, 2014a; Simpson et al. 2017), such that dust extinction may impact the observed [C ii] line luminosity at 158 μm. This would be in agreement with a larger apparent dust emission size at 1.0 mm than at 1.2 mm, as found above due to dust optical depth effects, but higher-resolution 1.0 mm imaging is required to further investigate this finding. We also find a moderately high dust temperature of (50 ± 9) K.³⁷

Based on the 1.2 mm size of GN10, we find a source surface area of (0.79 ± 0.44) kpc² (0.99 ± 0.30 kpc²; second quoted values indicate results from circular Gaussian fits). Its flux thus corresponds to a source-averaged rest-frame brightness temperature of T_b = (24.9 ± 8.1) K (20.0 ± 2.9 K) at rest frame 190 μm, or 50% ± 18% (40% ± 9%) of the dust temperature obtained from SED fitting. This suggests that the dust emission is likely at least moderately optically thick and that it fills the bulk of the source surface area within its inferred size. From our SED fit, we determine a dust optical depth of τ_{190 μm} = 0.69 ± 0.29 (i.e., 1−exp(−τ_{190 μm}) = $50{ \% }_{-17 \% }^{+13 \% }$), which is consistent with this finding. Using the value of L_IR derived in the previous subsection, this size corresponds to an IR luminosity surface density of Σ_IR = (7.5 ± 4.4) × 10¹² L_⊙ kpc⁻² (6.0 ± 2.0 × 10¹² L_⊙ kpc⁻²) or an SFR surface density of Σ_SFR = (750 ± 440) M_⊙ yr⁻¹ kpc⁻² (600 ± 210 M_⊙ yr⁻¹ kpc⁻²).

To obtain a dynamical mass estimate from our resolved line emission map, we have fitted visibility-plane dynamical models of a rotating disk to the [C ii] emission from GN10 (Figure 9; see Pavesi et al. 2018a for further details on the modeling approach). Rotating circular disk models of the [C ii] emission are generated through the KinMSpy code (Davis et al. 2013), super-imposed on continuum emission, which is fitted by a circular two-dimensional Gaussian function.³⁸ The fitting parameters are the disk center position, systemic velocity, gas dispersion, FWHM size of the spatial light profile of the Gaussian disk, maximum velocity, velocity scale length, inclination, position angle, and line flux. The continuum flux and FWHM size are determined based on the emission in the line-free channels. The fitting method employs the MCMC and nested sampling techniques as implemented in emcee (Foreman-Mackey et al. 2013) and MultiNest for python (Feroz et al. 2009; Buchner et al. 2014) to sample the posterior distribution of the model parameters and to calculate the model evidence. To optimize the parameter sampling, non-constraining, uniform priors are chosen for additive parameters, logarithmic priors for scale parameters, and a sine prior for the inclination angle. The data are fitted well by the disk model, leaving no significant residuals in the moment-0 map or spectrum. The results for all parameters are given in Table 4. We find a median dynamical mass of M_dyn = ${8.6}_{-2.8}^{+3.6}$(${4.5}_{-1.5}^{+2.1}$) × 10¹⁰ M_⊙ out to a radial distance of 5.5 (3.7) kpc from the center. Given the FWHM diameter of the [C ii] emission of ${6.4}_{-0.7}^{+0.7}$ kpc, the derived values are barely sufficient to host the estimated stellar mass of ∼1.2 ×10¹¹ M_⊙ if the stellar component has a similar extent as the [C ii] emission. This could either indicate that the stellar mass in GN10 is overestimated, e.g., due to the model fitting a solution that has a stellar population that is too old or the contribution of a dust-obscured active galactic nucleus (AGN) to the mid-infrared emission (see, e.g., Riechers et al. 2014b); or perhaps, that the kinematic structure of GN10 is not dominated by rotation (such that the dynamical mass is underestimated). Observations of the stellar and gas components at higher spatial resolution are required to distinguish between these scenarios.

Zoom In Zoom Out Reset image size

Table 4.  GN10 [C ii] Dynamical Modeling Parameters

Fit Parameter	Unit	Valuea
[C ii] FWHM diameter	arcsec	${1.03}_{-0.10}^{+0.11}$
	kpc	${6.4}_{-0.7}^{+0.7}$
Velocity scale length	arcsec	${0.25}_{-0.18}^{+0.29}$
	kpc	${1.6}_{-1.1}^{+1.8}$
Disk inclination	degrees	${80}_{-8}^{+7}$
Position angle	degrees	${206}_{-5}^{+5}$
Maximum Velocity	km s−1	${320}_{-80}^{+100}$
Gas dispersion	km s−1	${220}_{-25}^{+25}$
Dust FWHM diameter	arcsec	${0.56}_{-0.05}^{+0.05}$
	kpc	${3.5}_{-0.3}^{+0.3}$
Mdyn (r = 3.66 kpc)	1010 M⊙	${4.5}_{-1.5}^{+2.1}$
Mdyn (r = 5.49 kpc)	1010 M⊙	${8.6}_{-2.8}^{+3.6}$

Note.

^aValues as stated are the 50th percentiles. The lower and upper error bars are stated as the 16th and 84th percentiles, respectively.

Download table as:  ASCIITypeset image

Assuming optically thick emission, the integrated CO(J = 2 → 1) to CO(J = 1 → 0) brightness temperature ratio r₂₁ = 1.37 ± 0.45 (compared to 1.84 ± 0.38 in line profile peak brightness temperature) toward GN10 suggests that the CO(J = 1 → 0) line could be underluminous compared to expectations for thermalized or sub-thermal gas excitation. Given the comparable spatial resolution of the CO(J = 2 → 1) and CO(J = 1 → 0) observations, this is unlikely to be due to resolution effects. If significant, this effect may be due to the high cosmic microwave background (CMB) temperature at z = 5.3 (T_CMB ≃ 17.2 K), compared to the excitation potential above the ground level of the CO(J = 1 → 0) transition (which corresponds to an excitation temperature of T_ex = 5.5 K). The CMB acts as both a source of excitation and as a background for the CO emission.

Studies at z ∼ 2–3 have found evidence for enhanced CO(J = 1 → 0) line widths and line strengths in some DSFGs due to the presence of low-density, low surface brightness gas, for which a low-T_ex line like CO(J = 1 → 0) is an ideal tracer (e.g., Ivison et al. 2011; Riechers et al. 2011c). The CO brightness temperature itself is measured as a contrast to the CMB. Thus, low surface brightness emission as traced by CO(J = 1 → 0) may be disproportionately affected by CMB effects toward higher redshifts (e.g., by heating colder gas to T_CMB, thereby reducing the observable brightness temperature contrast), such that a preferential impact toward weakened CO(J = 1 → 0) line emission may be expected (e.g., da Cunha et al. 2013; Zhang et al. 2016).³⁹ Thus, a reduced CO(J = 1 → 0) line flux compared to higher-J lines due to the CMB appears plausible to explain the observed line spectra toward GN10. However, while not affected as strongly, some impacts on the CO(J = 2 → 1) line (having an excitation potential above the ground level of the corresponding T_ex = 16.6 K, i.e., ∼T_CMB(z = 5.3)) would also be expected in this scenario, yielding a reduced impact on r₂₁.

Apart from effects due to the CMB, another possible scenario is that the low-J CO line emission may not be optically thick in some regions. Given the higher expected optical depth in the CO(J = 2 → 1) line compared to CO(J = 1 → 0), r₂₁ > 1 would be possible in this scenario. In principle, self absorption due to cold gas in the foreground of the warmer gas located in the star-forming regions could also disproportionately affect the strength of low-J CO lines in some source geometries. While this finding is intriguing, higher signal-to-noise measurements in the future are desirable to further investigate this effect and its origin based on a detailed comparison of the line profiles.

Using the new and updated 1–3 mm photometry from this work, Pavesi et al. (2016) and Magnelli et al. (2019), we have re-fit the SED of AzTEC-3 with the same routine as used by Riechers et al. (2014a), which is similar to that used for GN10 above. We find T_dust = ${92}_{-16}^{+15}$ K, β_IR = 2.09 ± 0.21, ${\lambda }_{0}^{\mathrm{rest}}$ = ${181}_{-34}^{+33}$, α = ${10.65}_{-6.42}^{+6.69}$, and L_FIR = (1.12 ± 0.16) × 10¹³ L_⊙. These values are consistent with those found by Riechers et al. (2014a). We also find a total infrared luminosity of L_IR = (${2.55}_{-0.74}^{+0.73}$) × 10¹³ L_⊙. The relatively large uncertainties compared to L_FIR are due to the limited reliability of the Herschel/SPIRE 250–500 μm photometry near the peak of the SED. Taken at face value, this suggests that SFR_IR = (2500 ± 700) M_⊙ yr⁻¹.⁴⁰

Based on the 1.2 mm size of AzTEC-3, we find a source surface area of <(2.95 ± 0.45) kpc² (0.62 ± 0.17 kpc²; second quoted values indicate the results from circular Gaussian fits, which account for the compact component only). Its flux thus corresponds to a source-averaged rest-frame brightness temperature of T_b > (4.7 ± 0.7) K (16.8 ± 2.2 K) at rest frame 190 μm, or >5% ± 1% (18% ± 4%) of the dust temperature obtained from SED fitting. This suggests that AzTEC-3 has significant structure on scales significantly smaller than the ∼025 beam of our 1.2 mm observations. Using the value of L_IR derived in the previous section, this size corresponds to Σ_IR > (5.0 ± 1.6) × 10¹² L_⊙ kpc⁻² (18.0 ± 7.1 ×10¹² L_⊙ kpc⁻²), or Σ_SFR > (500 ± 160)M_⊙ yr⁻¹kpc⁻² (1800 ± 700 M_⊙ yr⁻¹ kpc⁻²).⁴¹

For HDF 850.1, the size of the gas reservoir measured from the high-resolution CO(J = 2 → 1) observations at its observed ${L}_{\mathrm{CO}(2-1)}^{{\prime} }$ implies a CO luminosity surface density of Σ_{CO(2 − 1)} = (4.8 ± 1.8) × 10⁵ L_⊙ kpc⁻². We also find a rest-frame CO(J = 2 → 1) peak brightness temperature of ${T}_{{\rm{b}}}^{\mathrm{CO}}$ = 1.6 ± 0.4 K at the ∼07 resolution of our observations, which agrees to within ∼7% with the source-averaged value. This modest value is consistent with the fact that the source is resolved over less than two beams in our current data.

Adopting the apparent L_IR = (8.7 ± 1.0) × 10¹² L_⊙ reported by Walter et al. (2012) and the rest-frame 158 μm dust continuum size of (09 ± 01) × (03 ± 01) reported by Neri et al. (2014), we find an apparent physical source size of (5.7 ± 0.6) × (1.9 ± 0.6) kpc² and a source-averaged Σ_IR = 6.0 ± 0.9 × 10¹¹ L_⊙ kpc⁻² for HDF 850.1. This corresponds to Σ_SFR ∼ (60 ± 10) M_⊙ yr⁻¹ kpc⁻².⁴² We also find a source-averaged rest-frame dust brightness temperature of T_b = 1.4 ± 0.3 K at rest frame 158 μm, or ∼4% ± 1% of its dust temperature. This suggests that the dust has significant substructure on scales much smaller than the ∼035 beam of the observations presented by Neri et al. (2014). It also is comparable to the CO brightness temperature on comparable scales.

Removing all strongly lensed and cluster-lensed galaxies from the full z > 5 sample to investigate potential differential lensing bias, we find a median value of (52.7 ± 6.7) K based on eight galaxies. This value is higher than the full sample median and, thus, does not provide direct evidence for preferential magnification of higher T_dust regions in the lensed subsample. Removing all galaxies selected at wavelengths shorter than 850 μm to investigate potential SED peak selection bias, we find (46.0 ± 4.7) K based on 10 galaxies. We thus do not find significant evidence that the different selection methods strongly bias the median T_dust for this z > 5 sample. Comparing sources that used the same SED fitting method as GN10 and AzTEC-3 (eight galaxies; see Table 6 caption for details) to the SPT sample (six galaxies), which was fitted with a common but somewhat different method (using graybody-only fits with fixed β_IR = 2 and λ₀ = 100 μm; e.g., Weiß et al. 2013), we find median values of (55.6 ± 4.5) and (46.0 ± 4.2) K, respectively. The former subsample is selected at shorter wavelengths (except one serendipitous discovery), while the latter contains a higher fraction of strongly and cluster-lensed systems (25% versus 83%). Thus, the different median T_dust cannot be directly attributed to differences in the fitting methods, but our study would be consistent with it playing a role. Overall, we thus find that the median value appears not to be significantly biased toward high values based on the heterogenous composition of the sample, while keeping in mind that even the combination of all current selection methods could entirely miss z > 5 DSFGs in undersampled regions of the parameter space.

Using the median redshift of the sample of 5.66, we find expected median values of T_dust of 37.3 and 54.6 K when applying the T_dust–z relations by Magnelli et al. (2014), their Equation (4), and Schreiber et al. (2018), their Equation (15),⁵⁷ for galaxies on the star-forming main sequence (Figure 11, top left panel). At face value, the Magnelli et al. relation appears consistent with the lowest T_dust sources, but the z > 5 sample would need to exhibit median star formation rates that are a factor of ∼90 higher (∼240 when only considering unlensed and weakly lensed systems) than those of main-sequence galaxies at the same redshift with the same M_⋆ for this relation to agree with the median T_dust value of our sample. While most of the current sample is likely to be in excess of the main sequence, there is no evidence for a difference of more than a factor of several in SFR compared to the main sequence. As such, current observations of z > 5 DSFGs would appear to be in favor of a stronger evolution in T_dust with redshift than indicated by this relation when assuming that redshift is the main property responsible for explaining the observed differences, unless systematic effects in the sample selection and SED fitting methods are more significant than currently known. On the other hand, a significant fraction of the z > 5 sample appears to be in reasonable agreement with expectations based on the Schreiber et al. relation, but the predicted median value is higher than that observed for the sample as a whole.⁵⁸ At least half of the sample has a lower T_dust than expected from this relation.

Both of the proposed T_dust–z relations have to be extrapolated significantly in redshift to match our sample, such that the observed mismatches are likely consistent with the true uncertainties of these relations. In particular, the T_dust of main-sequence galaxies at z > 5 is currently only poorly constrained by observations (e.g., Pavesi et al. 2016), such that it remains unclear that an offset in SFR from the main sequence can currently be meaningfully translated to expected differences in T_dust for dusty starbursts at these redshifts. Also, the increasing CMB temperature toward z > 5 provides a natural cutoff at low T_dust, both due to a reduced brightness temperature contrast and due to a higher contribution to dust heating (e.g., da Cunha et al. 2013). This effect leads to an overall change in SED shape and increase in the measured T_dust with redshift but is not accounted for in current extrapolations of the T_dust–z relations. On the other hand, it has recently been suggested that previously proposed T_dust–z relations could be an observational artifact due to selection effects (Dudzevičiūtė et al. 2020). If there were to be no trend in T_dust with redshift, another explanation would be required to explain the high median T_dust of the z > 5 sample.

We find a trend between T_dust and L_FIR in our sample (Figure 11), which spans about an order of magnitude in intrinsic L_FIR and a factor of 2.8 in T_dust. The data are consistent with an overall increase of L_FIR with T_dust. For z > 5 DSFGs with T_dust < 50.1 K (i.e., below the sample median), we find a median L_FIR = (7.3 ± 3.9) × 10¹² L_⊙, but for those with >50.1 K, we find (11.0 ± 5.0) × 10¹² L_⊙, as is consistent with the general trends found for lower-redshift DSFG samples (e.g., da Cunha et al. 2015; Simpson et al. 2017; Dudzevičiūtė et al. 2020). Taken at face value, this trend may at least partially, and perhaps entirely, explain the high median T_dust of the current z > 5 DSFG sample as being due to their relatively high median L_FIR. Indeed, the three intrinsically least-luminous sources in the sample with an average L_FIR = (3.4 ± 0.7) × 10¹² L_⊙ have an average T_dust = (38 ± 4) K, which is much closer to lower-redshift samples with comparable L_FIR, like ALESS and AS2UDS.⁵⁹ In the simplest scenario, the higher T_dust in the z > 5 DSFG sample thus could be due to their high median SFRs. This subject thus warrants further study based on the detailed dust properties of larger galaxy samples in the future.