COLDz: A High Space Density of Massive Dusty Starburst Galaxies ~1 Billion Years after the Big Bang
Dominik A. Riechers, Jacqueline A. Hodge, Riccardo Pavesi, Emanuele Daddi, Roberto Decarli, Rob J. Ivison, Chelsea E. Sharon, Ian Smail, Fabian Walter, Manuel Aravena, Peter L. Capak, Christopher L. Carilli, Pierre Cox, Elisabete da Cunha, Helmut Dannerbauer, Mark Dickinson, Roberto Neri, Jeff Wagg
DD RAFT VERSION M AY
7, 2020
Preprint typeset using L A TEX style AASTeX6 v. 1.0
COLDZ: A HIGH SPACE DENSITY OF MASSIVE DUSTY STARBURST GALAXIES ∼ D OMINIK
A. R
IECHERS , J
ACQUELINE
A. H
ODGE , R ICCARDO P AVESI , E MANUELE D ADDI , R OBERTO D ECARLI ,R OB J. I
VISON , C HELSEA
E. S
HARON , I AN S MAIL , F ABIAN W ALTER , M
ANUEL A RAVENA , P ETER
L. C
APAK ,C HRISTOPHER
L. C
ARILLI , P
IERRE C OX , E LISABETE DA C UNHA , H
ELMUT D ANNERBAUER ,M ARK D ICKINSON , R OBERTO N ERI , AND J EFF W AGG (Received 20/01/2020; Revised 17/04/2020; Accepted 21/04/2020) Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany Leiden Observatory, Leiden University, P.O. Box 9513, NL2300 RA Leiden, The Netherlands Laboratoire AIM, CEA/DSM-CNRS-Univ. Paris Diderot, Irfu/Service d’Astrophysique, CEA Saclay, Orme des Merisiers, F-91191 Gif-sur-Yvette cedex, France INAF - Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129, Bologna, Italy European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748 Garching, Germany Yale-NUS College, Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK National Radio Astronomy Observatory, Pete V. Domenici Array Science Center, P.O. Box O, Socorro, NM 87801, USA Núcleo de Astronomía, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile Spitzer Science Center, California Institute of Technology, MC 220-6, 1200 East California Boulevard, Pasadena, CA 91125, USA Cavendish Astrophysics Group, University of Cambridge, Cambridge, CB3 0HE, UK Sorbonne Université, UPMC Université Paris 6 and CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98bis boulevard Arago, F-75014 Paris, France Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia International Centre for Radio Astronomy Research, University of Western Australia, 35 Stirling Hwy, Crawley, WA 6009, Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain Universidad de La Laguna, Departamento de Astrofísica, E-38206 La Laguna, Tenerife, Spain NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA Institut de RadioAstronomie Millimétrique, 300 Rue de la Piscine, Domaine Universitaire, F-38406 Saint Martin d’Héres, France SKA Organization, Lower Withington, Macclesfield, Cheshire SK11 9DL, UK
ABSTRACTWe report the detection of CO( J =2 →
1) emission from three massive dusty starburst galaxies at z > z =5.183) and AzTEC-3 ( z =5.298),were previously known. We revise a previous redshift estimate for the third source GN10 ( z =5.303), which wehave independently confirmed through detections of CO J =1 →
0, 5 →
4, 6 →
5, and [C II ] 158 µ m emission withthe VLA and the NOrthern Extended Milllimeter Array (NOEMA). We find that two currently independentlyconfirmed CO sources in COLDz are “optically dark”, and that three of them are dust-obscured galaxies at z > ∼
60 arcmin , our results appear to imply a ∼ z > II ] emission compared to the dust continuum and higher [C II ]-to-far-infrared lu-minosity ratios in those galaxies with lower gas excitation. We thus find substantial variations in the conditionsfor star formation between z > ∼ ±
25% warmerthan starbursts at z =2–3 due to their enhanced star formation activity. Keywords: cosmology: observations — galaxies: active — galaxies: formation — galaxies: high-redshift —galaxies: starburst — radio lines: galaxies [email protected] a r X i v : . [ a s t r o - ph . GA ] M a y R IECHERS ET AL . INTRODUCTIONLuminous dusty star-forming galaxies (DSFGs) representthe most intense episodes of star formation throughout cos-mic history (see, e.g., Blain et al. 2002; Casey et al. 2014 forreviews). While the bulk of the population likely existed atredshifts z ∼ z > z > has now been detected at z > z ≥ z > z < z > In this work, galaxies with infrared luminosities of 10 < L IR < L (cid:12) are considered to be moderately luminous, and those above as luminous. nous DSFGs appear to contain large molecular gas reservoirsthat fuel their star formation, and to be significantly metal-enriched, leading to bright CO line emission. As such, thesesystems are preferentially picked up by panoramic molecularline scan surveys, and they may even dominate among detec-tions at the highest redshifts, where most other galaxy pop-ulations may exhibit only weak CO emission (e.g., Pavesi etal. 2019) due to a combination of lower characteristic galaxymasses at earlier epochs, lower metallicity (which is thoughtto lead to an increase in the α CO conversion factor, i.e., alower CO luminosity per unit molecular gas mass; see Bo-latto et al. 2013 for a review), and possibly lower CO lineexcitation (e.g., Daddi et al. 2015). Support for this idea wasprovided by the detection of the “optically-dark” z =5.183DSFG HDF 850.1 in a molecular line scan in the Hubble
Deep Field North (Walter et al. 2012), but the survey area of ∼ was too small to make more quantitative state-ments regarding the broader properties and space density ofsuch sources. To build upon this encouraging finding, and to provide abetter understanding of the true space densities of the mostdistant DSFGs, we here study the properties of dusty star-bursts at z > ∼
60 arcmin VLA COLDz survey data (Pavesi et al.2018b; Riechers et al. 2019, hereafter P18, R19). We re-port the detection of CO( J =2 →
1) emission from three sys-tems initially detected in 850 µ m and 1.1 mm continuumsurveys, HDF 850.1, AzTEC-3, and GN10 (Hughes et al.1998; Pope et al. 2005; Scott et al. 2008), two of whichhad previous correct redshift identifications through CO mea-surements (AzTEC-3 and HDF 850.1; Riechers et al. 2010;Capak et al. 2011; Walter et al. 2012). We also reporthigher-resolution observations of AzTEC-3 and HDF 850.1,and detailed follow-up observations of the third system, the“optically-dark” galaxy GN10, as well as CO line excita-tion modeling for the sample. Our analysis is used to con-strain the evolution of dust temperature with redshift, andthe space density of z > z > z > See, e.g., Calabro et al. (2019) and references therein for a more detaileddiscussion of the nature of such sources at lower redshift. Also, while its redshift was not known at the time, the telescope point-ing was chosen to include HDF 850.1. As such, this measurement did notconstitute an unbiased discovery of an “optically-dark” source. See http://coldz.astro.cornell.edu for additional infor-mation.
OLD Z : D USTY S TARBURSTS AT z > z > Λ CDM cosmology throughout, with H = 69.6 km s − Mpc − , Ω M = 0.286, and Ω Λ = 0.714 (Ben-nett et al. 2014). DATA2.1.
Very Large Array
COLDz survey data
CO( J =2 →
1) line emission ( ν rest =230.5380 GHz) towardHDF 850.1, AzTEC-3, and GN10 was detected in molecu-lar line scans with the VLA within the COLDz survey data(project IDs: 13A-398 and 14A-214; PI: Riechers). A de-tailed 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 ∼ ∼ in 7- and57-point mosaics with the VLA, respectively. Observationswere carried out under good Ka band weather conditions fora total of 324 hr between 2013 January 26 and 2015 Decem-ber 18 in the D and DnC array configurations, as well asreconfigurations between C, DnC, and D array. The cor-relator was set up with two intermediate frequency bands(IFs) of 4 GHz bandwidth (dual polarization) each in 3-bitmode, centered at the frequencies indicated above. Gaps be-tween individual 128 MHz sub-bands were mitigated throughfrequency switching. The radio quasars J1041+0610 andJ1302+5748 were observed for complex gain calibration andregular pointing corrections in the COSMOS and GOODS-North fields, respectively. The absolute flux scale was de-rived based on observations of 3C 286.Data reduction and imaging was performed with the CASA using the data pipeline version 1.2.0. Imagingthe data with natural baseline weighting yields typical cleanbeam sizes of 2.5 (cid:48)(cid:48) , with variations between individual point-ings and across the large bandwidth. Typical rms noise levelsare 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 (forwhich maps based on these data are shown in the following),the rms noise is 51 and 18 µ Jy beam − per 76 and 60 MHz( ∼
623 and 491 km s − ) binned channel at beam sizes of1.95 (cid:48)(cid:48) × (cid:48)(cid:48) and 2.46 (cid:48)(cid:48) × (cid:48)(cid:48) , respectively.2.1.2. GN10 CO (J=1 →
0) follow-up
We observed CO( J =1 →
0) line emission toward GN10 us-ing the VLA (project ID: 16A-015; PI: Riechers). Obser-vations were carried out under good Ku band weather con-ditions for a total of 11 hr during 5 tracks in C array be-tween 2016 February 02 and March 06. Two IFs with https://casa.nrao.edu + , and HNC( J =1 →
0) and CO( J =1 →
0) lines in GN10( ν rest =88.6318, 89.1885, 90.6636, and 115.2712 GHz). Ob-servations 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 per-formed on the gain calibrator approximately once per hour.The absolute flux scale was derived based on observations of3C 286.Data reduction and imaging was performed with the
CASA (cid:48)(cid:48) × (cid:48)(cid:48) at anrms noise level of 21 µ Jy beam − over 40 MHz (656 km s − )at the CO( J =1 →
0) line frequency. The rms noise near theHCN( J =1 →
0) frequency is ∼ µ Jy beam − over 4 MHz(85 km s − ). Imaging the data over the entire line-freebandwidth of 2.012 GHz yields an rms noise level of1.49 µ Jy beam − at a beam size of 2.10 (cid:48)(cid:48) × (cid:48)(cid:48) .2.1.3. GN10 1.3 cm and 6.6 mm continuum follow-up
We observed continuum emission at 22.8649 GHz (Kband) and 45.6851 GHz (Q band) toward GN10 (correspond-ing to 1.3 cm and 6.6 mm, respectively), using the VLA(project ID: AR693; PI: Riechers). Observations were car-ried out under good K and Q band observing conditions fora total of 44 hr between 2009 July 19 and 2010 January 05.K-band observations were conducted for 4 tracks in C array,totaling 28 hr, and Q-band observations were conducted for2 tracks in D array, totaling 16 hr. All observations used theprevious generation correlator, covering two IFs of 50 MHzbandwidth (dual polarization) each at the tuning frequencyand 300 MHz (K band) or 50 MHz (Q band) above, respec-tively, in quasi-continuum mode. Observations in C (D) ar-ray were carried out in short cycles, spending 150 s (200 s)on source, bracketed by scans spending 60 s on the gain cal-ibrator J13028+57486. Pointing was performed on the gaincalibrator approximately once per hour. The absolute fluxscale was derived based on observations of 3C 286.Data reduction and imaging was performed with the
AIPS package. Imaging the data with natural baseline weightingyields clean beam sizes of 1.15 (cid:48)(cid:48) × (cid:48)(cid:48) and 1.82 (cid:48)(cid:48) × (cid:48)(cid:48) atrms noise levels of 44 and 58 µ Jy beam − over 100 MHz inK and Q band, respectively.2.1.4. HDF 850.1 CO (J=2 →
1) high-resolution follow-up
We observed CO( J =2 →
1) line emission towardHDF 850.1 at higher spatial resolution using the VLA(project ID: 16A-014; PI: Riechers). Observations were car-ried out under good Ka band weather conditions for a total of These observations were tuned to the CO( J =1 →
0) and CO( J =2 → R IECHERS ET AL .8.8 hr during 4 tracks in C array between 2016 February 03and 20. Two IFs with 4 GHz bandwidth (dual polarization)each in 3-bit mode were centered at 34 GHz (correspondingto 8.8 mm) to cover the same frequency range as the D arrayobservations of the main survey. Gaps between sub-bandswere mitigated through frequency switching, using the sametwo setups with a relative shift of 16 MHz as in D array.Observations were carried out in short cycles, spending ∼
300 s on source, bracketed by scans spending ∼
100 s onthe gain calibrator J1302+5748. Pointing was performedon the gain calibrator approximately once per hour. Theabsolute flux scale was derived based on observations of3C 286.Data reduction and imaging was performed with the
CASA (cid:48)(cid:48) × (cid:48)(cid:48) at an rmsnoise level of 32.4 µ Jy beam − over 66 MHz (530 km s − ) atthe CO( J =2 →
1) line frequency.2.2.
NOEMA
GN10 CO (J=6 →
5) follow-up
We observed CO( J =6 →
5) line emission( ν rest =691.4731 GHz) toward GN10 using NOEMA (projectID: X–5; PI: Riechers). Observations were carried outunder good 3 mm weather conditions for 4 tracks in the Aconfiguration between 2014 February 18 and 24, using 6antennas (baseline range: 67–760 m). This yielded a totaltime of 13.8 hr (16500 visibilities) on source. Receiverswere tuned to 109.7037 GHz (corresponding to 2.7 mm).The correlator was set up with a bandwidth of 3.6 GHz (dualpolarization).Data reduction and imaging was performed with the GILDAS package. Imaging the data with natural or uniformbaseline weighting yields clean beam sizes of 0.82 (cid:48)(cid:48) × (cid:48)(cid:48) or 0.63 (cid:48)(cid:48) × (cid:48)(cid:48) at rms noise levels of 73 or 90 µ Jy beam − over 700 MHz (1912 km s − ), respectively. Imaging the datawith natural weighting over the entire line-free bandwidth of2.9 GHz yields an rms noise level of 31.6 µ Jy beam − .2.2.2. GN10 [C II ]( P / → P / ) follow-up We observed [C II ] 158 µ m line emission( ν rest =1900.5369 GHz) toward GN10 using NOEMA(project ID: W14FH; PI: Riechers). Observations were car-ried out under good 0.9 mm weather conditions for 1 trackin the C configuration on 2015 April 15, using 6 antennas(baseline range: 21–172 m). This yielded a total time of1.9 hr (2249 visibilities) on source. Receivers were tuned to301.524 GHz (corresponding to 1.0 mm). The correlator wasset up with a bandwidth of 3.6 GHz (dual polarization).Data reduction and imaging was performed with the GILDAS package. Imaging the data with natural or uniform baseline weighting yields clean beam sizes of1.01 (cid:48)(cid:48) × (cid:48)(cid:48) or 0.81 (cid:48)(cid:48) × (cid:48)(cid:48) at rms noise levels of 0.62 or0.71 mJy beam − over 800 MHz (795 km s − ), respectively.Imaging the data with natural or uniform weighting over theentire line-free bandwidth of 2.31 GHz yields rms noise lev-els of 324 or 365 µ Jy beam − .2.2.3. GN10 1.2 and 2 mm continuum follow-up
We observed continuum emission at 137.057 GHz(2.2 mm) and 250.5 GHz (1.2 mm) toward GN10 usingNOEMA (project IDs: T047 and T0B7; PI: Riechers). Ob-servations were carried out under good weather conditionsfor 3 tracks between 2009 June 4 and September 21 in theD configuration with 5 antennas (baseline range: 19–94 m)at 2 mm and for 2 tracks on 2011 January 23 and 24 in theA configuration with 6 antennas (baseline range: 51–665 m)at 1.2 mm. This yielded a total of 7.1 and 3.7 hr (17040 and8940 visibilities) 6 antenna-equivalent on source time at 2and 1.2 mm, respectively. Observations at 2 mm were carriedout with the previous generation correlator with a bandwidthof 1 GHz (dual polarization). Observations at 1.2 mm werecarried out with a bandwidth of 3.6 GHz (dual polarization).Data reduction and imaging was performed with the GILDAS package. Imaging the 2 mm data with natural base-line weighting yields a clean beam size of 3.7 (cid:48)(cid:48) × (cid:48)(cid:48) at anrms noise level of 95 µ Jy beam − over 1 GHz bandwidth.Imaging the 1.2 mm data with natural or uniform base-line weighting yields clean beam sizes of 0.45 (cid:48)(cid:48) × (cid:48)(cid:48) or0.38 (cid:48)(cid:48) × (cid:48)(cid:48) at rms noise levels of 0.35 or 0.43 mJy beam − over 3.6 GHz, respectively.2.2.4. Archival: GN10 CO (J=5 → Serendipitous CO( J =5 →
4) line emission( ν rest =576.2679 GHz) was observed toward GN10 usingNOEMA. These observations, taken from Daddi et al.(2009a), did not target GN10, which was offset by 9.7 (cid:48)(cid:48) or19.7 (cid:48)(cid:48) from the phase center for different tracks, yielding pri-mary beam attenuation factors of 1.08 or 1.30 respectively.Observations were carried out under good 3 mm weatherconditions for 4 tracks in the B, C, and D configurationsbetween 2008 May 04 and 2009 January 05, using 6 antennas(baseline range: 15–411 m). This yielded a total time of14.6 hr (25186 visibilities) on source. Receivers were tunedto 91.375 GHz (corresponding to 3.3 mm). The correlatorwas 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 aclean beam size of 2.64 (cid:48)(cid:48) × (cid:48)(cid:48) at an rms noise level of66.6 µ Jy beam − over 365.753 MHz (1200 km s − ) at the These observations were tuned to the CO( J =6 →
5) emission line at theprevious, incorrect redshift estimate.
OLD Z : D USTY S TARBURSTS AT z > Figure 1 . COLDz molecular line scan CO( J =2 →
1) spectra (histograms) of z > ∼
130 km s − ) spectral resolution.Line emission in GN10, AzTEC-3 and HDF 850.1 is detected at 8.6 σ , 14.7 σ , and 5.3 σ significance, respectively. phase center (96 µ Jy beam − at the position of GN10). Imag-ing the data over the line-free bandwidth yields an rms noiselevel of 55 µ Jy beam − .2.2.5. AzTEC-3 CO (J=5 →
4) high-resolution follow-up
We observed CO( J =5 →
4) line emission( ν rest =576.2679 GHz) toward AzTEC-3 using NOEMA(project ID: U0D0; PI: Riechers). Observations were carriedout under good 3 mm weather conditions for 2 tracks inthe A configuration between 2011 January 19 and February04, using 6 antennas (baseline range: 100–760 m). We alsoused previous observations (project ID: T–F; PI: Riechers)carried out for 1 track in the C configuration on 2010 April1, using 6 antennas (baseline range: 15–176 m; see Riecherset al. 2010 for additional details). This yielded a total timeof 8.8 hr (10560 visibilities; 5340 in A configuration) onsource. Receivers were tuned to 91.558 GHz (correspondingto 3.3 mm). The correlator was set up with a bandwidth of3.6 GHz (dual polarization).Data reduction and imaging was performed with the GILDAS package. Imaging the combined data with naturalbaseline weighting yields a clean beam size of 2.21 (cid:48)(cid:48) × (cid:48)(cid:48) at an rms noise level of 64 µ Jy beam − over 280 MHz(917 km s − ). Imaging the A configuration data only withuniform baseline weighting yields a clean beam size of1.39 (cid:48)(cid:48) × (cid:48)(cid:48) at an rms noise level of 96 µ Jy beam − over260 MHz (852 km s − ). Imaging the A (A+C) configura-tion 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 − .2.2.6. AzTEC-3 1.2 mm continuum follow-up
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 weatherconditions for 3 tracks between 2011 January 25 and Febru-ary 03 in the A configuration with 6 antennas (baseline range:100–760 m) at 1.2 mm. This yielded a total of 9.3 hr (11160visibilities) on source. Observations were carried out with abandwidth of 3.6 GHz (dual polarization).Data reduction and imaging was performed with the
GILDAS package. Imaging the data with natural or uniformbaseline weighting yields clean beam sizes of 0.62 (cid:48)(cid:48) × (cid:48)(cid:48) or 0.53 (cid:48)(cid:48) × (cid:48)(cid:48) at rms noise levels of 162 or 187 µ Jy beam − over 3.6 GHz, respectively. RESULTS3.1.
COLDz Molecular Line Scan CO (J=2 →
1) Detections
Our CO search in the COSMOS and GOODS-Northfields carried out as part of the COLDz molecular linescan survey (P18, R19) yielded four matches with mas-sive dusty star-forming galaxies initially selected in single- One match was found in the COSMOS field, and three matches werefound in the ∼ R IECHERS ET AL . Table 1 . Line fluxes and line luminosities for COLDz z > Line GN10 AzTEC-3 HDF 850.1 References
COLDz.GN.0 COLDz.COS.0 (cid:63)
COLDz.GN.31 (J2000.0) a (12:36:33.45, +62:14:08.85) (10:00:20.70, +02:35:20.50) (12:36:52.07, +62:12:26.49) I line L (cid:48) line I line L (cid:48) line I line L (cid:48) line [Jy km s − ] [10 K km s − pc ] [Jy km s − ] [10 K km s − pc ] [Jy km s − ] [10 K km s − pc ]CO( J =1 →
0) 0.054 ± b ± < < J =2 →
1) 0.295 ± ± ± ± ± ± ± ± ± ±
4, 5CO( J =5 →
4) 0.86 ± ± ± ± ± ± ± ± J =6 →
5) 0.52 ± ± ± ± ± ± J =7 →
6) 0.35 ± ± J =16 → < < Π / J =3/2 → ± ± P → P ) 0.14 ± ± P / → P / ) 17.6 ± ± ± ± ± ± ± c ± ± ± ± ±
1, 10, 5[NII]( P → P ) 0.46 ± ± J =2 →
1) centroid positions adopted from P18.b A Gaussian fit to the line profile formally suggests 0.054 ± − , but we consider these uncertainties to be somewhat optimistic due to the increasing noise leveltowards the blue edge of the bandpass. We thus adopt more conservative error bars based on the signal-to-noise ratio of the detection in the moment-0 map.c Main component only. (cid:63) Alternative ID: AS2COS0059.1 (Simpson et al. 2020).
References —[1] this work; [2] Wagg et al. (2007); [3] Pavesi et al. (2018b); [4] Riechers et al. (2010); [5] Walter et al. (2012); [6] Daddi et al. (2009a); [7] Decarli et al.(2014); [8] Riechers et al. (2014a); [9] Neri et al. (2014); [10] Pavesi et al. (2016). dish bolometer surveys with the JCMT at 850 µ m or 1.1 mm.One of the matches at 6.1 σ significance corresponds toCO( J =1 →
0) emission associated with the z =2.488 DSFGGN19 (Pope et al. 2005; Riechers et al. 2011c; Ivisonet al. 2011; see P18), and will not be discussed furtherhere. Two of the matches correspond to CO( J =2 →
1) emis-sion in the z =5.183 and z =5.298 DSFGs HDF 850.1 andAzTEC-3 (Fig. 1 and Table 1), which are detected at 5.3 σ and 14.7 σ significance, respectively. From Gaussian fits tothe line spectra and moment-0 maps, we find CO( J =2 → v FWHM =(490 ± ±
44) km s − for HDF 850.1 and AzTEC-3, yielding line fluxes of I CO ( − ) =(0.148 ± ± − , re-spectively. These flux levels are consistent with previous,lower-significance detections within the relative uncertainties(Riechers et al. 2010; Walter et al. 2012).Unexpectedly at the time of observation, we also detect anemission line at 8.6 σ significance toward the DSFG GN10(Fig. 1), previously thought to be at z =4.042 based on a singleline detection at 3 mm and photometric redshift information(Daddi et al. 2009a). We identify this line with CO( J =2 → z =5.303, which implies that the line detected byDaddi et al. (2009a) corresponds to CO( J =5 →
4) emission, rather than CO( J =4 →
3) emission. This identification wasconfirmed through the successful detection of CO( J =1 → J =6 → II ] emission at the same redshift, as de-scribed in detail below (see Figs. 2 and 3). This explains whyour earlier attempts to detect CO( J =1 → J =2 → J =6 →
5) emission at z =4.042 (see Sect. 2) were unsuc-cessful. 3.2. GN10 Follow-Up
Continuum Emission
We detect strong continuum emission toward GN10 at1.2 and 1.0 mm, and weak emission between 2.2 mm and1.9 cm (see Fig. 4, and Table 2). The flux keeps decreas-ing between 0.9 and 1.9 cm, and is > ∼ (cid:48)(cid:48) . By fitting two-dimensionalGaussian profiles to the emission in the visibility plane, wefind a size of (0.25 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.10 (cid:48)(cid:48) ± (cid:48)(cid:48) ), which corre-OLD Z : D USTY S TARBURSTS AT z > Figure 2 . VLA and NOEMA line spectra of GN10 ( z =5.3031; his-tograms) and Gaussian fits to the line profiles (black curves). Spec-tra are shown at resolutions of 66, 66, 75, 109, and 40 km s − (4, 8,23, 40, and 40 MHz; top to bottom ), respectively. Continuum emis-sion (9.55 ± II ] spectrumonly. sponds to (1.6 ± × (0.6 ± . A circular Gaussianfit provides a full width at half power (FWHP) diameter of0.18 (cid:48)(cid:48) ± (cid:48)(cid:48) , or (1.1 ± ∼
50 m,it appears unlikely that the 1.2 mm size measurement is bi-ased towards low values due to missing emission. Fits tothe lower-resolution ( ∼ (cid:48)(cid:48) beam size) data at 1.0 mm how-ever suggest a size of (0.58 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.50 (cid:48)(cid:48) ± (cid:48)(cid:48) ), corre-sponding to (3.6 ± × (3.1 ± , or a circular FWHPdiameter of 0.53 (cid:48)(cid:48) ± (cid:48)(cid:48) (3.3 ± ). The uncertaintiesfor this measurement may be limited by interferometric see-ing due to phase noise (which is not factored into the fittingerrors), such that we treat the 1.0 mm size measurement as anupper limit only in the following. We however note that, inprinciple, the dust emission at shorter wavelengths could ap-pear more extended due to an increasing dust optical depth aswell, as discussed further in Section 4. Also, the source shapecould significantly deviate from a Gaussian shape (such as ahigher-index Sérsic profile; see, e.g., discussion by Hodge etal. 2016), such that more complex fitting procedures couldyield different findings.Despite its strong dust continuum emission atsub/millimeter wavelengths, GN10 remains undetected Table 2 . GN10 continuum photometry.
Wavelength Flux density a Telescope Reference( µ m) (mJy)0.435 < × − HST /ACS 10.606 < × − HST /ACS 10.775 < × − HST /ACS 10.850 < × − HST /ACS 11.25 < × − Subaru/MOIRCS 11.60 < × − HST /NICMOS 12.15 < × − Subaru/MOIRCS 13.6 (1.29 ± × − Spitzer /IRAC 24.5 (2.07 ± × − Spitzer /IRAC 25.8 (2.96 ± × − Spitzer /IRAC 28.0 (5.30 ± × − Spitzer /IRAC 216 (17.5 ± × − Spitzer /IRS 224 (33.4 ± × − Spitzer /MIPS 270 < Spitzer /MIPS 3110 < Herschel /PACS 2160 < Herschel /PACS 2160 < Spitzer /MIPS 3250 b ± Herschel /SPIRE 2350 b ± Herschel /SPIRE 2500 b ± Herschel /SPIRE 2850 12.9 ± ± ± ± ± ± ± ± < < ± × − VLA 413100 < ± × − VLA 4210000 (35.8 ± × − VLA 2, 3a Limits are 3 σ .b De-blended fluxes. Uncertainties do not account for confusion noise, whichformally is 5.9, 6.3, and 6.8 mJy (1 σ ) at 250, 350, and 500 µ m, respectively(Nguyen et al. 2010), but also is reduced though the de-blending process. References —[1] Wang et al. (2009) and references therein; [2] Liu et al.(2018); [3] Dannerbauer et al. (2008) and references therein; [4] this work;[5] Daddi et al. (2009a). up to observed-frame 2.2 µ m, rendering it a “K-banddropout” (also see discussion by Wang et al. 2009; Daddiet al. 2009a, and references therein). Even sensitive space-based imaging up to 1.6 µ m with the WFC3 camera onthe Hubble Space Telescope ( HST ) from the CANDELSsurvey (Grogin et al. 2011) yields no hint of emission due R
IECHERS ET AL . [CII] velocity[-550;+330] Figure 3 . Velocity-integrated line contour maps overlaid on a
HST /WFC3 F160W continuum image from the CANDELS survey toward GN10.Contour maps are averaged over 656, 623, 1199, 1913, and 795 km s − , respectively. Contours start at ±
2, 3, 2, 3, and 4 σ , and are shown insteps of 1 σ =21, 51, 96, 73, and 710 µ Jy beam − for the first five panels, respectively. The synthesized beam size is indicated in the bottom leftcorner of each panel where applicable. Last panel shows CO( J =6 →
5) (cyan) and [C II ] (magenta) contours and is zoomed-in by a factor of 1.8compared to all other panels. The inset in the last panel shows a velocity map over the central 880 km s − of the [C II ] emission (created from80 km s − velocity channels and adopting a detection threshold of 9 mJy beam − ). Velocity contours are shown in steps of 50 km s − . to dust obscuration (Fig. 4), but mid-infrared continuumemission is detected with Spitzer /IRAC longward of 3.6 µ m(corresponding to ∼ Line Emission
We successfully detect CO( J =1 → J =2 → J =5 → J =6 → II ] emission towardGN10 at 3.3 σ , 8.6 σ , 6.9 σ , 7.3 σ , and 18.1 σ peak sig-nificance, respectively (Figs. 2 and 3 and Table 1; seeAppendix A for further details). In combination, theselines provide an unambiguous redshift identification.We extract the parameters of all emission lines fromGaussian fitting to the line profiles. All CO lines arefit with single Gaussian functions (see Appendix A forfurther details). We find line peak fluxes of S line =(74 ± ± ± ± µ Jy at FWHM line Contributions to the rest-frame optical light by a dust-obscured activegalactic nucleus in GN10 cannot be ruled out; see discussion in Sect. 4. widths of d v FWHM =(687 ± ± ± ± − for the CO J =1 →
0, 2 →
1, 5 →
4, and 6 → II ] line is fit with two Gaussiancomponents, yielding S line =(24.7 ± ± v FWHM =(617 ±
67) and (227 ± − for the red-and blue-shifted components, respectively. The parametersof the secondary component thus are only poorly con-strained, and we report [C II ] fluxes including and excludingthis component in the following. It may correspond toa gas outflow, or a close-by minor companion galaxy tothe DSFG. Considering only the main [C II ] component,the widths of all lines are consistent within the relativeuncertainties.Assuming the same widths as for the CO( J =2 →
1) line,we find 3 σ upper limits of < − for the HCN,HCO + , and HNC( J =1 →
0) lines. This implies HCN,HCO + , and HNC to CO line luminosity ratio limits of order The CO( J =6 →
5) spectrum shows excess positive signal at comparablevelocities as the secondary [C II ] component, but its significance is too lowto permit further analysis. The observations also cover the CCH( N =1 →
0) line, which is not de-tected down to a comparable limit (see Appendix).
OLD Z : D USTY S TARBURSTS AT z > ×
10 kpc × (NA)
10 kpc × (UN)
10 kpc D e c li n a t i o n ( J ) CO & 1.2 mm on HST/WFC3 F160W ×
10 kpc [CII] & 1.0 mm on HST/WFC3 F160W ×
10 kpc +62:13:54.014:00.006.012.018.024.012:36:32.033.034.035.0 µ m
10 kpc µ m Figure 4 . Rest-frame far-infrared continuum contour maps at observed-frame 1.0 ( top left ) and 1.2 mm ( top middle and right ; imaged usingnatural and uniform baseline weighting, respectively), overlaid on a
HST /WFC3 F160W continuum image toward GN10. Contours start at ± σ , and are shown in steps of 1 σ =365, 352, and 421 µ Jy beam − , respectively. The synthesized beam size is indicated in the bottomleft corner of each panel where applicable. The 4 σ peaks south of GN10 in the top left panel are due to sidelobe residuals given imperfectionsin the calibration. Bottom panels show overlays of 1.2 and 1.0 mm continuum (yellow) with CO( J =6 →
5) and [C II ] emission ( bottom left and middle ; same contours as in Fig. 3), and 1.2 mm contours (natural weighting; shown in steps of ± σ ) on Spitzer /IRAC 3.6 (inset) and 4.5 µ mimages ( bottom right ). The dashed gray box in the last panel indicates the same area as shown in the other panels. < ∼
20% ratios for distant starburstgalaxies (see, e.g., Riechers et al. 2007; Oteo et al. 2017a,and references therein).The integrated line fluxes and line luminosities derivedfrom these measurements are summarized in Table 1, to-gether with those of AzTEC-3 and HDF 850.1, includ-ing values from the literature. The CO( J =2 →
1) flux ofGN10 reported here is somewhat lower than that found byP18 using a different extraction method. Here we adoptour updated value for consistency. Given the more com-plex velocity profile of the [C II ] line in GN10, we adoptthe CO( J =2 → z =5.3031 ± σ with theCO( J =1 → J =5 → J =6 → The systemic velocity centroid of the [C II ] lineis slightly blueshifted, but emission is detected across thesame velocity range as in the CO J =1 → → The redshift uncertainties for the CO( J =1 → J =2 → J =5 → J =6 → II ] lines from fitting Gaussian functionsto the line profiles are 60, 30, 71, 70, and 24 km s − , respectively. line peak offset thus is likely mainly due to internal variationsin the [C II ]/CO ratio.By fitting two-dimensional Gaussian profiles to the [C II ]emission, we find a size of (1.04 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.19 (cid:48)(cid:48) ± (cid:48)(cid:48) ).This corresponds to (6.5 ± × (1.2 ± . Attempt-ing to fit the CO( J =6 →
5) emission observed at compa-rable spatial resolution yields a size of 0.42 (cid:48)(cid:48) × < (cid:48)(cid:48) (2.6 × < ), but the fit does not converge well due tothe moderate signal-to-noise ratio of the detection. A circu-lar fit is consistent with a point source within the uncertain-ties. This suggests that the [C II ] emission is resolved at leastalong the major axis, and that it appears to be more spatiallyextended than the mid- J CO and dust emission, which areconsistent with having comparable spatial extent within thecurrent uncertainties. The [C II ] emission shows a significantspatial velocity gradient across the line emission in the maincomponent alone (see Fig. 3).3.3. AzTEC-3 Follow-Up
Continuum Emission
We detect strong continuum emission toward AzTEC-3 at 1.2 mm, and weak emission at 3.3 mm at the samepeak position (Fig. 5). We measure a continuum flux of0 R
IECHERS ET AL . CO( J =2 −
1) on HST/ACS F814W ×
10 kpc
CO( J =5 −
4) on HST/ACS F814W × (AC configuration, NA)
10 kpc
CO( J =6 −
5) on HST/ACS F814W ×
10 kpc D e c li n a t i o n ( J ) CO( J =5 −
4) & [CII] 158 µ m on HST/ACS F814W × (A configuration, UN)0.63 ×
10 kpc × (A configuration, NA)0.62 × (NA)
10 kpc +2:35:19.020.021.022.010:00:20.620.720.8 × (UN) 0.63 ×
10 kpc
Figure 5 . Velocity-integrated line and rest-frame far-infrared continuum contour maps overlaid on a
HST /ACS F814W continuum image fromthe COSMOS survey toward AzTEC-3. CO( J =6 → II ], and 1.0 mm continuum data are adopted from Riechers et al. (2010; 2014a). Dataare imaged with natural (NA) baseline weighting unless mentioned otherwise. Line contour maps in the top row include all available dataand are averaged over 491, 917, 874 km s − ( left to right ), respectively. Contours start at ± σ , and are shown in steps of 1 σ =18, 64, and224 µ Jy beam − , respectively. Line contour map in the bottom left panel includes only A configuration CO( J =5 →
4) data (cyan; imaged withuniform weighting). Contours are averaged over 852 (CO), and 466 km s − (magenta; [C II ]), start at ± σ , and are shown in steps of 1 and4 σ , where 1 σ =96 and 200 µ Jy beam − , respectively. Continuum contour maps in the bottom middle panel include only A configuration 3.3 mmcontinuum data (white). Contours start at ± σ , and are shown in steps of 1 σ =24.8 (3.3 mm) and 162 µ Jy beam − (1.2 mm; gray), respectively.Contours in the bottom right panel start at ± σ , and are shown in steps of 1 (1.2 mm; white, imaged with uniform weighting) and5 σ (1.0 mm; gray), where 1 σ =187 and 58 µ Jy beam − , respectively. The synthesized beam size is indicated in the bottom left (white or cyancontours) or right (magenta or gray) corner of each panel. Bottom panels are zoomed-in by a factor of 3.3, except last panel, which is zoomed-inby a factor of ∼ (118 ± µ Jy at 3.3 mm (i.e., rest-frame 520 µ m) from theA configuration data, which is consistent with the 3 σ up-per limit of < ± µ Jy. The con-tinuum emission is spatially resolved along the major axis at1.2 mm by our observations with a synthesized beam size of ∼ (cid:48)(cid:48) . By fitting two-dimensional Gaussian profiles in theimage plane, we find a size of (0.45 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.05 (cid:48)(cid:48) + . (cid:48)(cid:48) − . (cid:48)(cid:48) ),which corresponds to (2.8 ± × (0.3 + . − . ) kpc , and wemeasure a 1.2 mm flux of (3.67 ± (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.17 (cid:48)(cid:48) + . (cid:48)(cid:48) − . (cid:48)(cid:48) ) and the dust spectral energydistribution shape found by Riechers et al. (2014a). A cir-cular Gaussian fit in the visibility plane (which gives moreweight to the longer baseline data) yields a FWHP diam- eter of 0.14 (cid:48)(cid:48) ± (cid:48)(cid:48) , or (0.9 ± ± ∼
75% of the emission emerge from a compact re-gion within the 2.5–3 kpc diameter dust reservoir. Line Emission
We detect CO( J =5 →
4) emission toward AzTEC-3 at15.0 σ peak significance (Fig. 5; CO J =2 → → I CO ( − ) =(0.97 ± − , which is con-sistent with a previous measurement by Riechers et al.(2010). From Gaussian fitting to the line profile (Fig. 6; CO J =2 → → II ] are also shown for reference), We caution that the source shape could significantly deviate from aGaussian shape, such that the structure and size of the dust emission couldbe more complex than indicated by these findings.
OLD Z : D USTY S TARBURSTS AT z > Figure 6 . VLA, NOEMA, and ALMA line spectra of AzTEC-3( z =5.2979; histograms) and Gaussian fits to the line profiles (blackcurves). CO( J =6 →
5) and [C II ] data are adopted from Riechers etal. (2010; 2014a). Spectra are shown at resolutions of 66, 33 55,and 20 km s − (8, 10, 20, and 20 MHz; top to bottom ), respectively. we obtain a peak flux of S CO ( − ) =(1.88 ± v FWHM =(396 ±
37) km s − , which is consistentwith the previously measured values, and those found for theCO( J =2 →
1) and [C II ] lines (Riechers et al. 2014a). The linemay show an excess in its red wing that is not captured wellby the Gaussian fit, but the significance of this feature is onlymoderate. A circular Gaussian fit in the visibility plane overthe entire width of the emission of 917 km s − (280 MHz)suggests a FWHP source diameter of 0.44 (cid:48)(cid:48) ± (cid:48)(cid:48) , whichcorresponds to (2.8 ± − (160 MHz) to capture only the maincomponent of the emission, we find 0.57 (cid:48)(cid:48) ± (cid:48)(cid:48) , whichcorresponds to (3.5 ± II ] emission of (0.63 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.34 (cid:48)(cid:48) + . (cid:48)(cid:48) − . (cid:48)(cid:48) )found by Riechers et al. (2014a) over a similar velocity range(466 km s − ). There appears to be a small spatial offsetbetween the peaks of the CO( J =5 →
4) and [C II ] emission(which is consistent with the dust continuum peak; Fig. 5).However, this offset becomes insignificant ( (cid:46) σ ) when ex-cluding the red CO line wing, i.e., when considering compa- rable velocity ranges, which results in a shift of the peak posi-tion by ∼ (cid:48)(cid:48) ( (cid:46) σ significance shift) relative to the broadervelocity range. If confirmed, the emission in the red line wingthus may correspond to a spatially offset kinematic compo-nent, but additional data are required to investigate this find-ing in more detail. D e c li n a t i o n ( J ) CO( J =2 −
1) on HST/WFC3 F160W ×
10 kpc
Figure 7 . Velocity-integrated CO( J =2 →
1) line contour map (VLAC array data only) overlaid on a
HST /WFC3 F160W continuum im-age from the CANDELS survey toward HDF 850.1. Contour mapis averaged over 530 km s − . Contours start at ± σ , and are shownin steps of 1 σ =32.4 µ Jy beam − . The synthesized beam size is in-dicated in the bottom left corner. HDF 850.1 Follow-Up
We have imaged the CO( J =2 →
1) emission in HDF 850.1at ∼ ∼ σ peak significance, and also spa-tially resolved (Fig. 7). A two-dimensional Gaussian fitto the emission in the image plane suggests a deconvolvedsize of (1.06 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.49 (cid:48)(cid:48) ± (cid:48)(cid:48) ), which correspondsto (6.7 ± × (3.1 ± . Given the moderate signal-to-noise ratio of the detection, the real uncertainty on the minoraxis extent of the emission is likely of order 50%. The orien-tation and size of the emission are consistent with that seenin the rest-frame 158 µ m dust continuum emission (Neri etal. 2014), and the peak positions agree to within one syn-thesized beam size. No stellar light is detected at the po-sition of the CO emission even in deep HST /WFC3 imag-ing at observed-frame 1.6 µ m due to dust obscuration, whichindependently confirms HDF 850.1 as an “optically-dark”galaxy. It is strongly blended with a bright foreground el-liptical galaxy in Spitzer /IRAC imaging longward of 3.6 µ m,such that only moderately constraining upper limits can beobtained (see also discussion by Pope et al. 2006).2 R IECHERS ET AL . ° Observed wavelength ∏ obs [ µ m] ° ° ° ° ° ° O b s e r v e d flu x F ∫ [ m J y ] CIGALE ModelMedian MBB FitArp 220M82GN10 data ° Rest wavelength ∏ rest [ µ m] ° ° ° ° ° ° O b s e r v e d flu x F ∫ [ m J y ] CIGALE ModelALESS ( z ª z =4.05)HFLS3 ( z =6.34)GN10 data Figure 8 . Spectral energy distribution of GN10 compared to well-known starbursts and a sample composite from the literature, showing itsunusual rest-frame optical/infrared colors even when compared to other dust-obscured galaxies.
Left:
Continuum photometry (points), overlaidwith median modified black body fit (MBB; black line) and
CIGALE (red line) models. SED fits for the nearby starbursts M82 (dotted orangeline) and Arp 220 (dashed magenta line) are shown for comparison (Silva et al. 1998), normalized to the 500 µ m flux of GN10. Right:
Samepoints and red line, but also showing SED fits for the z =4.06 and z =6.34 DSFGs GN20 (dashed pink line) and HFLS3 (dotted green line), as wellas a composite fit to DSFGs in the ALESS survey (dash-dotted gray line) for comparison (Magdis et al. 2011; Riechers et al. 2013; Swinbanket al. 2014), normalized in the same way.4. ANALYSIS OF INDIVIDUAL SOURCES4.1.
Properties of GN10
Given the new redshift identification of GN10, we here de-scribe in detail the determination of its key physical proper-ties. 4.1.1.
Spectral Energy Distribution
To extract physical parameters from the spectral energydistribution (SED) of GN10, we have followed two ap-proaches, as summarized in Fig. 8 and Table 3. First, we havefit the far-infrared peak of the SED with a modified blackbody (MBB) routine, where the MBB is joined to a smoothpower law with slope α toward observed-frame mid-infraredwavelengths (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 unityare used as fitting parameters. In addition, the flux f rest158 µ m at rest-frame 158 µ m is used as a normalization parameter.The Markov Chain Monte Carlo (MCMC) and Nested Sam-pling package EMCEE (Foreman-Mackey et al. 2013) is usedto explore the posterior parameter distribution. By integrat-ing the fitted functions, we obtain far-infrared ( L FIR ) and to-tal infrared ( L IR ) luminosities, which we then use to deter-mine dust-obscured star-formation rates SFR IR under the as-sumption of a Kennicutt (1998) conversion for a Chabrier(2003) stellar initial mass function. By adopting a mass ab-sorption coefficient of κ ν =2.64 m kg − at 125 µ m (Dunneet al. 2003), we also estimate a dust mass M dust from the fits, finding a high value in excess of 10 M (cid:12) (see Table 3).We also fit the full optical to radio wavelength photome-try of GN10 using CIGALE , the Code Investigating GALaxyEmission (Noll et al. 2009; Serra et al. 2011), using a broadrange in star-formation histories and metallicities and a stan-dard dust attenuation law (Calzetti 2001). We find parametersthat are broadly consistent with the MBB fitting where appli-cable. Due to its high level of dust obscuration, GN10 re-mains undetected shortward of observed-frame 3.6 µ m (rest-frame ∼ M (cid:63) fromthe CIGALE fit in the following. We find a high value in ex-cess of 10 M (cid:12) , i.e., ∼
100 times the dust mass, which weconsider to be reliable to within a factor of a few (see Ta-ble 3). From our MBB fits, we find that the dust turns opticallythick at ∼ µ m (i.e., between observed-frame 1.0 and1.2 mm; see Table 3), similar to the values found for other z > II ] line luminosity at 158 µ m. This would bein agreement with a larger apparent dust emission size at See Appendix B for an alternative M (cid:63) value obtained with the MAG - PHYS code, which we consider to be consistent within the uncertainties. Seealso Liu et al. (2018) for a discussion of the potential uncertainties associ-ated with the determination of M (cid:63) for the most distant DSFGs, and Simpsonet al. (2014; 2017) for a more detailed discussion of the uncertainties of M (cid:63) estimates for DSFGs in general. OLD Z : D USTY S TARBURSTS AT z > Table 3 . GN10 MBB and
CIGALE
SED modeling parameters.
Fit Parameter unit value a T dust K 50.1 + . − . β IR + . − . λ rest0 µ m 168 + − α + . − . f rest158 µ mb mJy 8.3 + . − . M dust M (cid:12) + . − . L FIRc L (cid:12) + . − . L IRd L (cid:12) + . − . SFR IR M (cid:12) yr − + − M (cid:63) e M (cid:12) ± th percentiles. Lower and upper error bars are stated as 16 th and 84 th percentiles, respectively.b Fit normalization parameter.c Integrated between rest-frame 42.5 and 122.5 µ m.d Integrated between rest-frame 8 and 1000 µ m.e Parameter adopted from CIGALE . Quoted fitting uncertainties underestimate sys-tematic effects, and thus, are not adopted to make any conclusions in the following.
Table 4 . GN10 [C II ] dynamical modelingparameters. Fit Parameter unit value a [C II ] FWHM diameter arcsec 1.03 + . − . kpc 6.4 + . − . Velocity scale length arcsec 0.25 + . − . kpc 1.6 + . − . Disk inclination degrees 80 + − Position angle degrees 206 + − Maximum Velocity km s − + − Gas dispersion km s − + − Dust FWHM diameter arcsec 0.56 + . − . kpc 3.5 + . − . M dyn ( r =3.66 kpc) 10 M (cid:12) + . − . M dyn ( r =5.49 kpc) 10 M (cid:12) + . − . a Values as stated are 50 th percentiles. Lower and up-per error bars are stated as 16 th and 84 th percentiles,respectively. ±
9) K. See Appendix for an alternative, luminosity-averaged T dust value ob-tained by fitting multiple dust components with the MAGPHYS code.
Star-Formation Rate Surface Density
Based on the 1.2 mm size of GN10, we find a sourcesurface area of (0.79 ± (0.99 ± ; secondquoted values indicate results from circular Gaussian fits).Its flux thus corresponds to a source-averaged rest-framebrightness temperature of T b =(24.9 ± ± µ m, or 50% ±
18% (40% ± τ µ m =0.69 ± − exp( − τ µ m )=50% + − ), which is consistent with thisfinding. Using the L IR derived in the previous subsection,this size corresponds to an IR luminosity surface density of Σ IR =(7.5 ± × L (cid:12) kpc − (6.0 ± × L (cid:12) kpc − ),or a SFR surface density of Σ SFR =(750 ± M (cid:12) yr − kpc − (600 ± M (cid:12) yr − kpc − ).4.1.3. Dynamical Mass Estimate
To obtain a dynamical mass estimate from our resolvedline emission map, we have fitted visibility-plane dynamicalmodels of a rotating disk to the [C II ] emission from GN10(Fig. 9; see Pavesi et al. 2018a for further details on the mod-eling approach). Rotating circular disk models of the [C II ]emission are generated through the KinMSpy code (Davis etal. 2013), super-imposed on continuum emission which is fit-ted by a circular two-dimensional Gaussian function. Thefitting parameters are the disk center position, systemic ve-locity, gas dispersion, FWHM size of the spatial light pro-file of the Gaussian disk, maximum velocity, velocity scalelength, inclination, position angle, and line flux. The contin-uum flux and FWHM size are determined based on the emis-sion in the line-free channels. The fitting method employsMCMC and Nested Sampling techniques as implemented in
EMCEE (Foreman-Mackey et al. 2013) and M
ULTI N EST for python (Feroz et al. 2009; Buchner et al. 2014) to sam-ple the posterior distribution of the model parameters andto calculate the model evidence. To optimize the param-eter sampling, non-constraining, uniform priors are chosenfor additive parameters, logarithmic priors for scale parame-ters, and a sine prior for the inclination angle. The data arefitted well by the disk model, leaving no significant resid-uals in the moment-0 map or spectrum. The results for allparameters are given in Table 4. We find a median dynam-ical mass of M dyn =8.6 + . − . (4.5 + . − . ) × M (cid:12) out to a radialdistance of 5.5 (3.7) kpc from the center. Given the FWHMdiameter of the [C II ] emission of 6.4 + . − . kpc, the derived val-ues are barely sufficient to host the estimated stellar mass of ∼ × M (cid:12) if the stellar component has a similar extent We include the continuum emission in the fitting to properly accountfor uncertainties associated with continuum modeling and subtraction.
IECHERS ET AL . F l u x d e n s i t y [ m J y ] DATA F l u x d e n s i t y [ m J y ] MODEL F l u x d e n s i t y [ m J y ] RESIDUAL
DATAmoment 0 = 12.5 kpc DATAmoment 1 = 12.5 kpc DATAmoment 2 = 12.5 kpc MODELmoment 0 = 12.5 kpc MODELmoment 1 = 12.5 kpc MODELmoment 2 = 12.5 kpc RESIDUALmoment 0 = 12.5 kpc RESIDUALmoment 1 = 12.5 kpc RESIDUALmoment 2 = 12.5 kpc Figure 9 . Visibility space dynamical modeling results for the [C II ] emission toward GN10. Intensity (moment-0), velocity (moment-1), andvelocity dispersion (moment-2) maps ( left ) for the data ( top ), model ( middle ), and “data–model” residuals ( bottom ), and spectra (histograms; right ). Median model fits and residuals are shown. Moment-0 contours are overlaid on all panels, and are shown in steps of ± σ . Data areimaged with “natural” baseline weighting. Beam size is shown in the bottom left corner of each map panel. Continuum emission was includedas part of the additional fitting parameters for the rotating disk model. to the [C II ] emission. This could either indicate that the stel-lar mass in GN10 is overestimated, e.g., due to the modelfitting a solution which has too old a stellar population orthe contribution of a dust-obscured active galactic nucleus(AGN) to the mid-infrared emission (see, e.g., Riechers et al.2014b), or that the kinematic structure of GN10 is not dom-inated by rotation (such that the dynamical mass is underes-timated). Observations of the stellar and gas components athigher spatial resolution are required to distinguish betweenthese scenarios.4.1.4. “Underluminous” CO (J=1 →
0) Emission?
Assuming optically thick emission, the integratedCO( J =2 →
1) to CO( J =1 →
0) brightness temperature ratio r =1.37 ± ± J =1 →
0) line could be underluminous compared toexpectations 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 toresolution effects. If significant, this effect may be due to thehigh cosmic microwave background (CMB) temperature at z =5.3 ( T CMB (cid:39) J =1 →
0) transition (which corre-sponds to an excitation temperature of T ex =5.5 K). The CMBacts as both a source of excitation and as a background forthe CO emission.Studies at z ∼ J =1 →
0) line widths and line strengths in some DSFGsdue to the presence of low density, low surface brightnessgas, for which a low- T ex line like CO( J =1 →
0) is an idealtracer (e.g., Riechers et al. 2011a; Ivison et al. 2011). TheCO brightness temperature itself is measured as a contrast tothe CMB. Thus, low surface brightness emission as tracedby CO( J =1 →
0) may be disproportionately affected by CMBeffects toward higher redshifts (e.g., by heating colder gas to T CMB , thereby reducing the observable brightness tempera-OLD Z : D USTY S TARBURSTS AT z > J =1 →
0) line emission may be expected (e.g., daCunha et al. 2013; Zhang et al. 2016). Thus, a reducedCO( J =1 →
0) line flux compared to higher- J lines due to theCMB appears plausible to explain the observed line spec-tra toward GN10. However, while not affected as strongly,some impact on the CO( J =2 →
1) line (having an excita-tion potential above ground corresponding to 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 sce-nario is that the low- J CO line emission may not be opti-cally thick in some regions. Given the higher expected opti-cal depth in the CO( J =2 →
1) line compared to CO( J =1 → r > J CO lines in some sourcegeometries. While this finding is intriguing, higher signal-to-noise measurements in the future are desirable to furtherinvestigate this effect and its origin based on a detailed com-parison of the line profiles.4.2.
Properties of AzTEC-3
We here update some of the key properties on AzTEC-3,based on the new information available in this work.4.2.1.
Spectral Energy Distribution
Using the new and updated 1–3 mm photometry fromthis work, Pavesi et al. (2016), and Magnelli et al.(2019), we have re-fit the SED of AzTEC-3 with thesame routine as used by Riechers et al. (2014a), whichis similar to that used for GN10 above. We find T dust =92 + − K, β IR =2.09 ± λ rest0 =181 + − , α =10.65 + . − . ,and L FIR =(1.12 ± ) × L (cid:12) . These values are consistentwith those found by Riechers et al. (2014a). We also finda total infrared luminosity of L IR =(2.55 + . − . ) × L (cid:12) . Therelatively large uncertainties compared to L FIR are due to thelimited reliability of the
Herschel /SPIRE 250–500 µ m pho-tometry near the peak of the SED. Taken at face value, thissuggests SFR IR =(2500 ± M (cid:12) yr − . Star-Formation Rate Surface Density
Based on the 1.2 mm size of AzTEC-3, we find a sourcesurface area of < (2.95 ± (0.62 ± ; secondquoted values indicate results from circular Gaussian fits, In this scenario, a disproportional impact on emission from cold dustdue to a reduced contrast toward and increased heating by the CMB wouldalso be expected, resulting in a higher apparent dust temperature due tochanges in the SED shape. This would be consistent with the relativelyhigh measured T dust of GN10. Given the uncertainties on L IR , previous works adopted L FIR to deter-mine the SFR of AzTEC-3. We adopt this updated value here for internalconsistency of the analysis presented in this work. which account for the compact component only). Its fluxthus corresponds to a source-averaged rest-frame brightnesstemperature of T b > (4.7 ± ± µ m, or > ±
1% (18% ± ∼ (cid:48)(cid:48) beam of our 1.2 mm observations. Using the L IR derived in the previous section, this size corresponds to Σ IR > (5.0 ± × L (cid:12) kpc − (18.0 ± × L (cid:12) kpc − ),or Σ SFR > (500 ± M (cid:12) yr − kpc − (1800 ± M (cid:12) yr − kpc − ). Properties of HDF 850.1
The new high-resolution CO( J =2 →
1) data and a combina-tion of constraints from the literature allow us to determinesome additional key properties of HDF 850.1, as detailed inthe following.4.3.1.
CO Luminosity Surface Density
For HDF 850.1, the size of the gas reservoir measuredfrom the high-resolution CO( J =2 →
1) observations at itsobserved L (cid:48) CO ( − ) implies a CO luminosity surface den-sity of Σ CO ( − ) =(4.8 ± × L (cid:12) kpc − . We also finda rest-frame CO( J =2 →
1) peak brightness temperature of T COb =1.6 ± ∼ (cid:48)(cid:48) resolution of our observations,which agrees to within ∼
7% with the source-averaged value.This modest value is consistent with the fact that the sourceis resolved over less than two beams in our current data.4.3.2.
Star-Formation Rate Surface Density
Adopting the apparent L IR =(8.7 ± × L (cid:12) reportedby Walter et al. (2012) and the rest-frame 158 µ m dustcontinuum size of (0.9 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.3 (cid:48)(cid:48) ± (cid:48)(cid:48) ) reported byNeri et al. (2014), we find an apparent physical sourcesize of (5.7 ± × (1.9 ± and a source-averaged Σ IR =6.0 ± × L (cid:12) kpc − for HDF 850.1. This corre-sponds to Σ SFR ∼ (60 ± M (cid:12) yr − kpc − . We also finda source-averaged rest-frame dust brightness temperatureof T b =1.4 ± µ m, or ∼ ±
1% of itsdust temperature. This suggests that the dust has significantsubstructure on scales much smaller than the ∼ (cid:48)(cid:48) beamof the observations presented by Neri et al. (2014). It also iscomparable to the CO brightness temperature on comparablescales. We here assume that the fraction of the flux contained by the compactcomponent at 1.2 mm is constant with wavelength, which may yield a lowerlimit on Σ IR if this component is warmer than the rest of the source, or anupper limit if it has a higher optical depth. The lensing magnification factor cancels out of the surface density cal-culations, such that the resulting properties are conserved under lensing, bar-ring potential differential lensing effects.
IECHERS ET AL . DISCUSSION OF THE COLDZ Z > z > z > Star-Formation Rate Surface Densities
The SFR surface densities of Σ SFR =(750 ± ± M (cid:12) yr − kpc − found above for GN10 and thecentral region of AzTEC-3 are even higher than the val-ues found for some other compact starbursts like ADFS-27 ( z =5.66; 430 ± M (cid:12) yr − kpc − ) and HFLS3 ( z =6.34;480 ± M (cid:12) yr − kpc − ; Riechers et al. 2013; 2017). Thesource-averaged value of > (500 ± M (cid:12) yr − kpc − forAzTEC-3 is comparable to these sources. Like these sys-tems, GN10 and AzTEC-3 thus show the properties ex-pected for so-called “maximum starbursts”. At face value,the peak Σ SFR in AzTEC-3 may slightly exceed the ex-pected Eddington limit for starburst disks that are sup-ported by radiation pressure (e.g., Thompson et al. 2005;Andrews & Thompson 2011), but is potentially consistentunder the assumption of a more complex source geome-try. On the other hand, the high Σ SFR value of GN10could be boosted by an obscured AGN contribution to thedust heating. GN10 exhibits strong 0.5–8 keV X-ray emis-sion. Using Equation 15 of Ranalli et al. (2003), its ob-served (absorption-corrected) rest-frame 2–10 keV X-rayluminosity of L X =5.6(12.5) × erg s − (Laird et al. 2010)corresponds to a SFR X of 1100(2500) M (cid:12) yr − . Given itsSFR IR =1030 + − M (cid:12) yr − , its L X remains consistent with in-tense star formation, but a contribution from a modestly lu-minous obscured AGN cannot be ruled out, depending onthe (relatively uncertain) absorption correction required. Thiswould also be consistent with a possible excess mid-infraredemission due to an obscured AGN, which may be favored bysome of the SED fits.In contrast, the source-averaged Σ SFR in HDF 850.1 isonly ∼ (60 ± M (cid:12) yr − kpc − , i.e., ∼
12 times lower thanin GN10, and ∼
30 times lower than the peak value inAzTEC-3. On the other hand, the values for GN10 andHDF 850.1 would become comparable when assuming thelarger 1.0 mm continuum size limit for GN10 instead of thesmaller 1.2 mm continuum size adopted above. As such,the source-averaged Σ SFR may be more similar when ac- AzTEC-3 is not detected at X-ray wavelengths. From fitting a Galactic absorption plus power-law model, an effectivephoton index Γ =1.40 + . − . was found for GN10, but it was not possible tosimultaneously fit the absorbing column N H and Γ due to the limited pho-ton counts. The data are consistent with no intrinsic absorption within theuncertainties, such that the absorption-corrected L X could be considered anupper limit (Laird et al. 2010) – or, alternatively, a lower limit in case it isheavily absorbed. Figure 10 . CO excitation ladders (points) and LVG modeling (lines)of the three COLDz z > J =2 →
1) line in AzTEC-3. Sources are slightly offset from each other in J upper to improveclarity. The models (solid lines) of AzTEC-3 and HDF 850.1 consistof two gas components each (dashed lines). counting for potentially extended dust emission if present,and thus, significantly lower than in AzTEC-3 on aver-age. Such more modest, “sub-Eddington” Σ SFR on few kpcscales would be consistent with what is found from high-resolution studies of lower- z DSFG samples on average (e.g.,Simpson et al. 2015; Hodge et al. 2016; 2019). For ref-erence, using the sizes and SFRs for z > Σ SFR =(200 ± M (cid:12) yr − kpc − ,where the error corresponds to the bootstrap uncertainty onthe median. Extending the sample down to z > Σ SFR =(170 ± M (cid:12) yr − kpc − . However, a more precisemeasurement of the 1.0 mm continuum morphology of GN10is required to make firm statements about extended dustemission. The intriguingly high peak Σ SFR of AzTEC-3 willbe explored further in future work.5.2.
CO Large Velocity Gradient Modeling
To constrain the line radiative transfer properties of thethree COLDz z > T kin and the gas density ρ gas as free parameters (Fig. 10).For all calculations, the H ortho-to-para ratio was fixed to3:1, the CMB temperature was set to the values appropriatefor our targets (16.85–17.18 K at z =5.183 to 5.303), and theFlower et al. (2001) CO collision rates were adopted. We alsoused a ratio between the CO abundance and velocity gradientof 10 − pc ( km s − ) − (see, e.g., Weiß et al. 2005b; 2007;Riechers et al. 2006).OLD Z : D USTY S TARBURSTS AT z > T kin =40 K and ρ gas =10 . cm − , suggesting the presenceof gas at moderate excitation. This model would imply thatthe CO emission should fill the bulk of the [C II ]-emitting re-gion, which perhaps suggests that the warm, compact nucleustraced by the 190 µ m dust emission is embedded in a moreextended cold gas reservoir. For AzTEC-3, we adopt the LVG model discussed byRiechers et al. (2010). This model consists of a low-excitation component with T kin =30 K and ρ gas =10 . cm − and a high-excitation component with T kin =45 K and ρ gas =10 . cm − . Based on the updated CO( J =2 →
1) flux,we reduce the area filling factor of the lower-temperaturecomponent. This component now contributes only ∼
20% tothe CO( J =2 →
1) line flux. This model suggests an expectedCO( J =1 →
0) flux of 0.057 Jy km s − , or a CO( J =2 →
1) toCO( J =1 →
0) line ratio of r =0.91. If the high-excitation gaswere distributed over the same 3.9 × area as the [C II ]emission imaged by Riechers et al. (2014a), the LVG modelwould imply an area filling factor of close to 50%. Assumingthe extent of the main CO( J =5 →
4) component determinedabove instead would imply a filling factor of just above 30%.The low-excitation gas, on the other hand, would need tobe distributed over an area at least twice as large in linearextent as the [C II ] emission, corresponding to an area ap-proximately half as large as the current observed limit on theCO( J =2 →
1) diameter of ∼ J =1 → − ),but the model does not reproduce the ratio between theCO( J =6 →
5) and CO( J =7 →
6) lines well. Thus, to improvethe model fit, it is necessary to add a second, high-excitationgas component with T kin =70 K and ρ gas =10 . cm − . Themoderate-and high-excitation components in this model fill ∼
20% and ∼ II ]-emitting area imaged by Neriet al. (2014), respectively. Assuming the CO( J =2 →
1) sizedetermined above instead yields ∼ The model is consistent with high optical depth in all CO lines. Both radial variations in temperature and optical depth may contributeto the observed effect; see also discussion by Calistro Rivera et al. (2018). Since no detections above J upper =7 exist, we caution that the parametersof this second component are only poorly constrained by the data. tion between z > J =5 →
4) – FIR luminosity relation for star-forming galaxies (Fig. 11; e.g., Daddi et al. 2015; Liu et al.2015), which shows that their level of star-formation activity(as traced by L FIR ) is consistent with what is expected basedon the properties of the warm, dense molecular gas (as tracedby CO J =5 → Gas Masses, Gas Surface Densities, Gas Fractions,Gas-to-Dust Ratios, Gas Depletion Times, andConversion Factor
Adopting the CO( J =1 →
0) fluxes from the LVG model-ing and α CO =1 M (cid:12) (K km s − pc ) − to convert the result-ing CO( J =1 →
0) luminosities to molecular gas masses M gas ,we find M gas =(7.1 ± ± ± × M (cid:12) for GN10, AzTEC-3, and HDF 850.1, respectively, wherethe overall uncertainties are dominated by systematic uncer-tainties in α CO . For GN10, this corresponds to 83% + − of the dynamical mass estimate, which is representative un-der the assumption that the CO( J =1 →
0) emission is as ex-tended as the [C II ] emission on which the dynamical mod-eling was carried out. Adopting the dynamical mass esti-mate of M dyn =(9.7 ± × M (cid:12) found for AzTEC-3 byRiechers et al. (2014a), its updated gas mass correspondsto 59% ±
11% of the dynamical mass under the same as-sumptions. Using an isotropic virial estimator (e.g., En-gel et al. 2010) and the [C II ] sizes and line widths ofits two components derived by Neri et al. (2014), wefind a combined dynamical mass of M dyn =7.5 × M (cid:12) for HDF850.1. As such, its gas mass corresponds to29% ±
18% of the dynamical mass under the same assump-tions as for GN10. Conversely, based on their observed L (cid:48) CO ,the dynamical masses of GN10, AzTEC-3, and HDF 850.1at face value imply upper limits of α CO < < < M (cid:12) (K km s − pc ) − , respectively. For HDF 850.1,the size of the gas reservoir measured from the high-resolution CO( J =2 →
1) observations implies a molecular gasmass surface density of Σ gas =(1.3 ± × M (cid:12) pc − . Thissuggests that HDF 850.1 follows the spatially-resolved star Gas and dust masses and sizes for HDF 850.1 are corrected by a fac-tor of µ L =1.6 to account for gravitational lensing magnification (Neri et al.2014). Statistical uncertainties in r from the LVG modeling likely contributeonly at the 10%–15% level. The uncertainties of this estimate are dominated by systematic effectsdue to the choice of fitting method, which corresponds to at least 50%.
IECHERS ET AL . Table 5 . Derived properties for COLDz z > Quantity unit GN10 AzTEC-3 HDF 850.1 D dust (major × minor axis) kpc (1.6 ± × (0.6 ± ± × (0.3 + . − . ) (5.7 ± × (1.9 ± IR M (cid:12) yr − + − ±
700 870 ± Σ SFR M (cid:12) yr − kpc − ± > ±
160 60 ± Σ peakSFR M (cid:12) yr − kpc − ±
440 1800 ± M gas M (cid:12) ± ± ± M dust M (cid:12) + . − . + . − . ± M gas / M dust + − + − ± M dyn M (cid:12) + . − . ± ± f gas = M gas / M dyn + − ±
11% 29% ± α dynCO M (cid:12) (K km s − pc ) − < < < τ dep = M gas /SFR IR Myr 70 ±
15 22 ± ± L CII / L FIR − ± ± ± L CII / L CO ( − ) ±
650 2400 ±
300 4500 ± formation law, which has been determined at lower redshift(e.g., Hodge et al. 2015).Given the dust masses of (11.1 + . − . ), (2.66 + . − . ), and(1.72 ± × M (cid:12) for GN10, AzTEC-3, and HDF 850.1(this work; Riechers et al. 2014a; Walter et al. 2012), the gasmasses correspond to gas-to-dust ratios of ∼ + − , 215 + − ,and 130 ±
50, respectively, which are in the range of valuesfound for nearby infrared-luminous galaxies (e.g., Wilson etal. 2008). Given the SFRs of (1030 + − ), (2500 ± ± M (cid:12) yr − , (this work; Riechers et al. 2014a; Wal-ter et al. 2012; Neri et al. 2014), they also yield gas de-pletion times of τ dep =(70 ± ± ±
15) Myr,respectively. This is consistent with short, (cid:46)
100 Myr star-burst phases, as also typically found for lower-redshift DS-FGs (e.g., Simpson et al. 2014; Dudzeviciute et al. 2020),with the lowest value found for the source showing the high-est gas excitation.5.4. [C II ]/FIR and [C II ]/CO Luminosity Ratios To further investigate the difference in conditionsfor star formation among the three COLDz z > II ]( P / → P / )/FIR and[C II ]( P / → P / )/CO( J =1 →
0) luminosity ratios (e.g.,Stacey et al. 1991; 2010). We adopt the CO( J =1 → II ]/FIR ratiosof L CII / L FIR =(2.5 ± ± ± × − forGN10, AzTEC-3, and HDF 850.1, respectively, i.e., a factorof ∼ II ]/FIR ratioswith increasing CO excitation. In GN10, the CO( J =6 → II ] emission. Thus, the source-averaged CO excitation may be comparatively low, but it maybe significantly higher in the nuclear region with the highest Σ SFR . Given the larger spatial extent of the [C II ] emission,the [C II ]/FIR ratio likely also is substantially lower in thisregion than the source average, consistent with the apparentglobal trend between [C II ]/FIR and CO excitation.We also find [C II ]/CO ratios of L CII / L CO ( − ) (cid:39) ± ± ± ∼ The lowestratio is found for AzTEC-3, i.e., the source with the highestgas excitation (but the two lower-excitation sources showcomparable values), the highest dust temperature (seeTable 6), and the highest peak Σ SFR . This may suggest areduced [C II ] line strength due to a stronger, more intenseinterstellar radiation field, but it could also be due to dustextinction of the [C II ] line emission or a low brightnesstemperature contrast between the [C II ] and 158 µ m dustemission in the nuclear starburst region.5.5. Dust Temperatures
We find that the COLDz z > z =5.18) has a compara-tively low T dust =(35 ±
5) K (Walter et al. 2012), while GN10 Adopting the measured CO( J =1 →
0) luminosity instead of that inferredfrom the LVG modeling would yield a ratio of (cid:39) We caution that a direct comparison in terms of physical propertiesto lower redshifts (as typically done with photon dominated region modelscalculated at z =0) is not straight forward due to the reduction in the intensityof low- J CO line emission at z > OLD Z : D USTY S TARBURSTS AT z > Table 6 . Properties of known z > Name redshift lensed? µ La selection b S S T dustd µ L L FIRe L FIR references[mJy] [mJy] [K] [10 L (cid:12) ] [10 L (cid:12) ]HXMM-30 5.094 strongly ... 250–500 µ m 55 ± ± µ m 116.3 ± ± ± + . − . ... 3, 1 HDF 850.1 f µ m <
14 7.8 ± ± ± ± ± µ m 212 ±
15 125 ± ± ± ± ± µ m 140 ± ± ± ± ± ± ± ± ± ± + . − . ± AzTEC-3 ± ± g + − ± ± GN10 µ m 12.4 ± ± ± + . − . ± ± ± ± ± µ m 24.0 ± ± + . − . ± (cid:46) ± ± ± ± ± + . − . ± ± h ± ± + . − . ± ± ± ± ± ± ± ± ± ± ± ± + . − . ... 9ID 85001929 5.847 no? ... SED template i ± ± + . − . + . − . ± µ m 97 ± k ± k ... ... ... 1, 2G09-83808 6.0269 strongly 8.2 ± µ m 44.0 ± ± ± ± ± ± l µ m 47.3 ± ± + . − . ± ± f ± ± ± + . − . ± OTE —Uncertainties in T dust are the statistical uncertainties from the fitting procedures adopted by different authors and do not account for systematic uncertainties dueto differences between the procedures (see references provided for additional details). The method adopted here for GN10 is virtually identical to those used for fittingHELMS_RED_4, HLock-102, AzTEC-3, ADFS-27, CRLE, ID 85001929, and HFLS3.a Lensing magnification factor. Modeled with two source components when multiple values are given. The uncertainties of full SED-based quantities for strongly-lensedsources (i.e., those with µ L > µ m bands were typically preferentially followed up if they showed high 850 µ m–1.3 mm fluxes.c S or S are used where S is not available.d Dust temperatures with small uncertainties typically are due to keeping β IR , λ , or both fixed in the fitting process. HLS0918 was fitted with an optically-thin modelonly, which may underestimate T dust . HDF 850.1 was fitted with a MAGPHYS -based model, and its dust peak is only constrained by upper limits shortward of 850 µ m.e Apparent far-infrared luminosity, i.e., not corrected for lensing magnification where applicable.f No uncertainties are quoted in the original works. Here we assume 20% uncertainty for the calculation of derived quantities.g Simpson et al. (2020) report an updated value of 8.3 ± µ m (D. Riechers et al., in prep.).i Source was selected at 850 µ m (originally reported at 1.2 mm; Bertoldi et al. 2007), but followed up based on a photometric redshift obtained by fitting a template SEDbased on HFLS3 (Riechers et al. 2013).j Values obtained by fitting the de-blended fluxes at 100 µ m–3 mm using a method virtually identical to that adopted for GN10. The T dust is compatible with the value of61 ± z source that contributes ∼
30% of the flux at 870 µ m, and likely a higher fraction at 500 µ m.l Updated value based on visibility-plane lens modeling of the dust and gas emission. References —[1] D. Riechers et al., in prep.; [2] Oteo et al. (2017b); [3] Asboth et al. (2016); [4] Walter et al. (2012); [5] this work; [6] Rawle et al. (2014); [7] Riechers etal. (2013); [8] Dowell et al. (2014); [9] Strandet et al. (2016); [10] Riechers et al. (2014a); [11] Riechers et al. (2017); [12] Pavesi et al. (2018a); [13] Jin et al. (2019);[14] Fudamoto et al. (2017); [15] Zavala et al. (2018); [16] Strandet et al. (2017).
IECHERS ET AL . . . . . . redshift z L F I R [ L Ø ] COLDz
COLDz z > DSFGs (apparent)COLDz z > DSFGs (intrinsic)Literature z > DSFGs (apparent)Literature z > DSFGs (intrinsic)
30 40 50 60 70 80 90 100 110 T dust [K] L F I R [ L Ø ] intrinsic L / T COLDz z > DSFGsLiterature z > DSFGs L CO(5 ° [K kms ° pc ] L F I R [ L Ø ] Low- z galaxies (L15) z = 1 . ° . DSFGs (CW13)COLDz z > DSFGsLiterature z > DSFGs . . . . . redshift z T du s t [ K ] GN10AzTEC-3HDF 850.1 HFLS3ALESS z ª . DSFGsArp 220 z > median S r e l a t i o n M M S r e l a t i o n M
14 5 x SS F R M S COLDz z > DSFGsLiterature z > DSFGs
Figure 11 . Dust temperature ( top left ) and far infrared luminosity ( top right ) of COLDz z > z > topright panel; see Table 6 for references), with the exception of one source discovered recently from de-blended single-dish catalogs (Jin etal. 2019). There is a weak trend between dust temperature and far infrared luminosity within the current sample, which is largely consistentwith a standard L ∝ T scaling (dashed line in bottom left panel, shown scaled to the median values and with ± top left panel show representative dust temperatures and uncertainty rangesfor ALESS z ∼ z > × higher than the MS, respectively.The COLDz z > J =5 →
4) – FIR luminosity relation of nearby galaxies ( bottom right ; green symbols, upperlimit arrows without symbols, and dashed line show the
Herschel /SPIRE spectroscopy sample and relation by Liu et al. 2015), suggesting thatthe properties of the dense, warm gas in their star-forming regions are as expected for starburst galaxies. Magenta stars show a lower redshiftcomparison DSFG sample, updated from the compilation of Carilli & Walter (2013), featuring data from Andreani et al. (2000), Weiß et al.(2005a; 2009), Carilli et al. (2010), Riechers et al. (2011c; 2011b), Cox et al. (2011), Danielson et al. (2011), and Salome et al. (2012). Valuesare corrected for gravitational magnification unless mentioned otherwise.
OLD Z : D USTY S TARBURSTS AT z > z =5.30) shows a moderate T dust =(50 ±
9) K, especially whencompared to the relatively high T dust =(92 + − ) K displayed byAzTEC-3 ( z =5.30; Riechers et al. 2014a). It is interesting to place these galaxies into the more gen-eral context of all currently known z > Thefull z > T dust =33–46 Ksuch as CRLE ( z =5.67; Pavesi et al. 2018a) and the grav-itationally lensed HLS0918 ( z =5.24; Rawle et al. 2014),SPT 2319–55, 2353–50, 0243–49, and 0311–58 ( z =5.29,5.58, 5.70, and 6.90; Strandet et al. 2016; 2017), and G09-83808 ( z =6.03; Fudamoto et al. 2017), sources with higher T dust =50–59 K like GN10 such as ADFS-27 ( z =5.66; Riech-ers et al. 2017), ID 85001929 ( z =5.85; Jin et al. 2019),HFLS3 ( z =6.34; Riechers et al. 2013), and the gravitation-ally lensed HLock-102 ( z =5.29; Riechers et al. 2013; Dow-ell et al. 2014) and SPT 0346–52 and 2351–57 ( z =5.66 and5.81; Strandet et al. 2016), and sources with high T dust >
60 Kcloser to AzTEC-3 like the strongly-lensed HELMS_RED_4( z =5.16; T dust =67 ± z > ± with an average value of 50.3 K) in comparison to thebulk of the population at z ∼ ± ∼ ± Given the heteroge-nous selection and some differences in the fitting methodsused, we investigate in the following to what degree theycould be responsible for this difference, but we find no trendsthat would be sufficient to explain this difference. We thenexplore if the higher T dust in the z > L FIR , and find that the latter is likely sufficientto explain the observed differences to lower-redshift samples.5.5.1.
Potential Biases due to Gravitational Lensing, SelectionWavelength, or SED Fitting Method
Removing all strongly-lensed and cluster-lensed galaxiesfrom the full z > ± Note that an optically-thin fitting procedure suggests a more modestdust temperature of (53 ±
5) K compared to the general fitting result adoptedhere, but it also yields a worse fit to the data. Jin et al. (2019) report ID 20010161, a source with a candidate redshiftof z =5.051 based on a single emission line. While the identification is plau-sible, we do not include it in the current data compilation until the redshiftis confirmed. Uncertainties are given as the median absolute deviation. Luminosity-averaged T dust from energy balance modeling for ALESSgalaxies are higher, with an average of (43 ±
2) K (da Cunha et al. 2015).We here adopt those from Swinbank et al., since their derivation is moreconsistent with the methods used for the z > median, and thus, does not provide direct evidence for pref-erential magnification of higher T dust regions in the lensedsubsample. Removing all galaxies selected at wavelengthsshorter than 850 µ m to investigate potential SED peak selec-tion bias, we find (46.0 ± T dust for this z > β IR =2 and λ =100 µ m; e.g.,Weiß et al. 2013), we find median values of (55.6 ± ± T dust cannot be directly attributed to differences inthe fitting methods, but our study would be consistent with itplaying a role. Overall, we thus find that the median value ap-pears not to be significantly biased towards high values basedon the heterogenous composition of the sample, while keep-ing in mind that even the combination of all current selectionmethods could entirely miss z > Comparison to Proposed T dust –z Relations
Using the median redshift of the sample of 5.66, we findexpected median values of T dust of 37.3 and 54.6 K when ap-plying the T dust – z relations by Magnelli et al. (2014; theirEquation 4) and Schreiber et al. (2018; their Equation 15) for galaxies on the star-forming “main sequence” (Fig. 11,top left). At face value, the Magnelli et al. relation appearsconsistent with the lowest T dust sources, but the z > ∼
90 higher ( ∼
240 when only considering un-lensed and weakly-lensed systems) than those of “main se-quence” galaxies at the same redshift with the same M (cid:63) forthis relation to agree with the median T dust value of our sam-ple. While most of the current sample is likely to be in ex-cess of the “main sequence”, there is no evidence for a dif-ference by more than a factor of several in SFR compared tothe “main sequence”. As such, current observations of z > T dust with redshift than indicated by this relation when assum-ing that redshift is the main property responsible for explain-ing the observed differences, unless systematic effects in thesample selection and SED fitting methods are more signifi-cant than currently known. On the other hand, a significantfraction of the z > Value obtained after scaling by their Equation 6 to obtain modifiedblack-body temperatures.
IECHERS ET AL .ment with expectations based on the Schreiber et al. relation,but the predicted median value is higher than observed forthe sample as a whole. At least half the sample have lower T dust than expected from this relation.Both of the proposed T dust – z relations have to be extrap-olated significantly in redshift to match our sample, suchthat the observed mismatches are likely consistent with thetrue uncertainties of these relations. In particular, the T dust of “main sequence” galaxies at z > T dust for dusty starbursts at these red-shifts. Also, the increasing CMB temperature toward z > T dust , both due to a reducedbrightness temperature contrast and due to a higher contribu-tion to dust heating (e.g., da Cunha et al. 2013). This effectleads to an overall change in SED shape and increase in themeasured T dust with redshift, but is not taken into accountin current extrapolations of the T dust – z relations. On the otherhand, it has recently been suggested that previously proposed T dust – z relations could be an observational artifact due to se-lection effects (Dudzeviciute et al. 2020). If there were to beno trend in T dust with redshift, another explanation would berequired to explain the high median T dust of the z > Investigation of T dust –L FIR
Relations
We find a trend between T dust and L FIR in our sample(Fig. 11), which spans about an order of magnitude in in-trinsic L FIR , and a factor of 2.8 in T dust . The data are consis-tent with an overall increase of L FIR with T dust . For z > T dust < L FIR =(7.3 ± × L (cid:12) , but for those with > ± × L (cid:12) , as is consistent withthe general trends found for lower-redshift DSFG samples(e.g., da Cunha et al. 2015; Simpson et al. 2017; Dudze-viciute et al. 2020). Taken at face value, this trend mayat least partially, and perhaps entirely explain the high me-dian T dust of the current z > L FIR . Indeed, the three intrinsi-cally least luminous sources in the sample with an average L FIR =(3.4 ± × L (cid:12) have an average T dust =(38 ±
4) K,which is much closer to lower-redshift samples with compa-rable L FIR like ALESS and AS2UDS. In the simplest sce-nario, the higher T dust in the z > The relation also suggests 38.7 K at z ∼ ∼
20% higher thanthe ALESS sample median. For reference, the z ∼ L FIR =(3.0 ± × L (cid:12) at T dust =(32 ±
1) K (Swinbank et al. 2014).
Number Densities and Space Densities
The overall population of DSFGs was found to be suffi-cient to explain the early formation of most luminous localand intermediate redshift elliptical galaxies, under the as-sumption that DSFGs are their intensely star-forming pro-genitors (e.g., Simpson et al. 2014). On the other hand, thehighest-redshift luminous DSFGs are the likely progenitorsof some of the most massive compact “quiescent” galaxies at z > z > J =2 →
1) emission in an areaof ∼
60 arcmin , compared to four independently-confirmedsources detected in CO( J =1 →
0) at z =2.0–2.8 (one of whichis a DSFG) in the same area. On the one hand, the comov-ing volume covered by the survey is ∼ J =2 →
1) emission at z =4.9–6.7 compared tothose of CO( J =1 → L (cid:48) CO limit (P18, R19). Also, one of the fields was chosento include AzTEC-3 at z =5.298, which thus needs to be re-moved from comparisons to avoid potential bias. As such,there are significantly fewer serendipitously discovered CO-rich galaxies per unit volume in the higher-redshift bin. Onthe other hand, once AzTEC-3 is removed from further con-sideration, the two other sources alone still imply a signif-icant excess in the CO luminosity function at z ∼ z > COLDz Number Density and Space Density — Conservativelyrestricting the analysis to the GOODS-North field (i.e.,excluding AzTEC-3 and its environment), the two de-tected sources alone correspond to a number density of ∼ (150 ± − for z > ± × − Mpc − , within the volume probedin CO( J =2 →
1) emission by COLDz ( ∼ in The redshift of HDF 850.1 was known at the time the COLDz surveywas carried out, but the field was chosen based on the availability of the deep
HST /WFC3-IR CANDELS data, not the presence of this source. Thus, it isnot excluded from the space density discussion, since it should not introducea bias.
OLD Z : D USTY S TARBURSTS AT z > S limit [mJy] − − N u m b e r d e n s i t y [ d e g − ] µ m-selected K-band dropoutsCOLDz z> z> z> z> z> z> z> z> z> z> z> z> z> Figure 12 . Number densities of high-redshift dusty galaxy samplesdown to a given 500 µ m flux limit, showing that the values foundfor COLDz z > z DSFG samples. Greenand gray squares show red
Herschel source samples identified byRiechers et al. (2013), Dowell et al. (2014), Asboth et al. (2016),and Ivison et al. (2016). Open symbols show the full samples,and filled symbols show estimated fractions or lower limits of z > z >
6) sources where applicable. Red symbols show lower limits at z > ) and the HerMESfields from which sources were selected (325 deg total). The dot-ted black line shows “K-band dropout” galaxies identified at 870 µ mfor reference (Dudzeviciute et al. 2020), many of which are likely at z ∼ z > z > z > S > µ m, 850, and 850 µ m, respectively(Simpson et al. 2014; 2020; Miettinen et al. 2017; Dudzeviciute etal. 2020). The COLDz and COSMOS samples are not selected at500 µ m; the symbols are shown at the lowest detected 500 µ m fluxin the samples instead of a formal flux limit. GOODS-North alone; R19). Comparison to Bright Samples from Large-Area Surveys — Atbright flux levels (i.e., down to a limit of S >
52 mJy at500 µ m), Asboth et al. (2016) find 477 candidate z > based on the selection of red Herschel sources (i.e., sources with S < S < S ). Basedon a power law fit to their sample, they find a number den-sity of 2.8 deg − . Follow-up observations of a subsampleof 188 sources at longer wavelengths suggest that at least21, or ∼ z > > − . Applying a similar selection down to Quoted uncertainties here and in the following are the standard devia-tion of a Poisson distribution given the number of sources detected. a limit of S >
30 mJy over an area of 21 deg , Dowell et al.(2014) find a number density of 3.3 deg − based on 38 suchsources, of which at least ∼
11% are spectroscopically con-firmed to be at z >
4. This corresponds to a number densityof > − . Based on the z =6.34 DSFG HFLS3 alone,Riechers et al. (2013) estimate that the number density ofsuch sources may drop to values as low as 0.014 deg − at z >
6. For a sample of 109 out of 708 sources found withsimilar selection criteria ( S / S ≥ S / S ≥ < S <
73 mJy (with the exception of one source at S =118 mJy) selected over ∼
600 deg , Ivison et al. (2016)find that about one-third of their candidates (or ∼ − when correcting for completeness) are likely at z >
4. Thisleads them to infer a space density of ∼ × − Mpc − for z > < − for z > S >
28 mJy (Fig. 12). Indeed, thenumber density of spectroscopically confirmed z > Her-schel sources in the HerMES fields provides a lower limit ofonly > (0.021 ± − down to S >
24 mJy. In ad-dition, some simulations appear to suggest that the numberdensities of red
Herschel sources could be biased towardshigh values due to selection effects introduced by the noiselevel and limited resolution of the
Herschel data (e.g., Bether-min et al. 2017). If significant, this would render currentspace density estimates upper limits.With S =(12.4 ± ± <
14 mJy forGN10, AzTEC-3, and HDF 850.1 (Table 6; Walter et al.2012; Riechers et al. 2014a; Liu et al. 2018), respectively,all of the COLDz z > > S (Fig. 12; seealso Fig. 11 for a comparison to other z > L FIR ). However, the more than two orders of magnitude differencein the implied space and number densities likely either sug-gests that such sources are overrepresented in the COLDzsurvey area due to large scale structure even after neglect-ing AzTEC-3, that the counts fall very steeply with flux inthe 10 (cid:46) S (cid:46)
30 mJy regime at z >
5, or that current selec-tion methods underlying these previous estimates are moreincomplete than currently thought. Substantially larger sur-vey areas than accessible with surveys like COLDz wouldbe required to distinguish between these possibilities. Un-til such surveys are available, further insight may be gainedthrough a comparison to samples selected over smaller ar-eas, but down to lower flux levels, than the bright samplesdiscussed so far.
Comparison to Faint Samples — Targeted studies of DSFGpopulations down to S (cid:38) z > These surveys however may contain strongly-lensed analogs of COLDzsources (Fig. 11).
IECHERS ET AL .with S / =(12.0 ± ± ± > × − Mpc − ; Cooke et al. 2018) for 4 < z < This lowerlimit provides a closer match to the space density implied bythe COLDz data at face value, but it likely still falls shortgiven the lower redshift and lower flux limit of this compar-ison sample. On the other hand, interferometric follow-upsurveys of flux-limited DSFG samples such as ALESS andAS2UDS find number densities of “K-band dropout” sourcesof ∼ − , or ∼
70 deg − (Fig. 12; e.g., Simpson etal. 2014; 2017; Dudzeviciute et al. 2020; also see discussionof similar sources, e.g., by Dannerbauer et al. 2004; Frayeret al. 2004; Franco et al. 2018). These sources may be con-sidered analogs of HDF 850.1 and GN10 in some respects,but they typically have lower 870 µ m fluxes and currentlymostly lack spectroscopic confirmation and molecular gasmass measurements. Also, many of them are consistent with z ∼ Based on largely photometricredshifts, the ALMA-COSMOS, ALESS, and AS2UDS sam-ples selected at 1.1 mm, 870 µ m, and 850 µ m provide spacedensity estimates of (24 ± ± ±
7) deg − for candidate z > ± ± ±
4) deg − for z > ±
1) and(2.7 ± down to S > z > z >
5, re-spectively, based on the AS2COSMOS survey (Simpson etal. 2020). These samples suggest that the space densities ofDSFGs decline by factors of ∼ z =4 to 5, and that atmost ∼ µ m–1.1 mm selected DSFGs at the cur-rently achieved flux limits appear to be at z >
5. Moreover, Also see Aravena et al. (2010) for comparable estimates at z =3–5 basedon single-dish 1.2 mm-selected sources. The Ivison et al. (2016) sample has single-dish fluxes of S =8–71 mJy, including upper limits consistent with this range. Thus, it reachesdown to the long wavelength fluxes of the COLDz z > µ m. Simpson et al. (2017) estimate a comoving space density of ∼ − Mpc − for AS2UDS DSFGs with a median S =(8.0 ± S to the COLDz z > z > The Simpson et al. (2017) AS2UDS sample contains two “K-banddropout” sources with S > There appears to be a trend that the brightest S -selected DSFGs tendto be found at higher redshifts (e.g., Younger et al. 2007; Simpson et al.2020). spectroscopic follow-up of two of the three z > < z <
5, reducingthe true ALESS z > z > > (1.5 ± − , i.e., (cid:46) z > ∼ S (cid:38) z > Potential Role of Selection Effects — A possible concern for thecomparisons between COLDz and continuum surveys (be-yond the increased impact of gravitational lensing on thebrighter populations) is that samples selected at a singlesub/millimeter wavelength may be incomplete at a given L FIR due to selection effects, in particular those related to dusttemperature or optical depth, which could lead to an underes-timate in the volume density of luminous z > µ m (i.e., (cid:46) µ m rest-frame at z >
5) and below, couldmiss sources with low dust temperatures or high dust opti-cal depths at the highest redshifts, for which the dust SEDpeak shifts to significantly longer wavelengths than 500 µ m.An example of such sources are “870 µ m risers”, which arerequired to either be extremely luminous or strongly gravita-tionally lensed to remain detectable at ≤ µ m (Riechers etal. 2017). Indeed, despite its moderately high dust tempera-ture, GN10 has a 870 µ m flux that is just below the “870 µ mriser” selection criterion. On the other hand, samples selectedat long wavelengths, e.g., 2 mm (i.e., (cid:46) µ m rest-frame at z >
5) and above, may be anticipated to be more complete atthe highest redshifts (see, e.g., Staguhn et al. 2014, or similardiscussion by Casey et al. 2018), but they could miss sourcesat high dust temperatures. As an example, a recent sensitive2 mm survey in COSMOS finds several dusty sources whichare claimed to be at a relatively high median redshift com-pared to shorter wavelength selected samples as expected, butit missed AzTEC-3 (i.e., the highest significance CO detec-tion in the COLDz survey and the most distant DSFG cur-rently known in the area surveyed at 2 mm) due to its rela-tively high dust temperature (Magnelli et al. 2019). Sur- Values after accounting for gravitational magnification of HDF 850.1. GN10 and HDF 850.1 however are solidly and tentatively detected ina 2 mm survey in GOODS-North, consistent with their lower dust tempera-tures (Staguhn et al. 2014).
OLD Z : D USTY S TARBURSTS AT z > µ m (i.e., (cid:46) µ m in the rest frame) are lessincomplete due to variations in T dust at a given L FIR at z > z > Potential Contribution of Cosmic Variance due to Clustering — Inprinciple, the COLDz survey may have led to the identifi-cation of an unexpectedly high number of z > z (cid:39) ∆ v (cid:39) − ), whichcorresponds to a comoving distance of 61.4 Mpc. As such,there already is evidence for large scale structure, which,in principle, could include both DSFGs in GOODS-North.On the other hand, the redshift difference between GN10and AzTEC-3 is only d z (cid:39) ∆ v (cid:39)
240 km s − , or2.5 Mpc comoving distance), and thus, much smaller, butthey are located in widely different parts of the sky. As such,a physical connection between these sources is not consid-ered to be possible. Thus, the modest redshift differencebetween GN10 and HDF 850.1 does not necessarily implya physical connection (and indeed, would require a ratherlarge correlation length). Moreover, it remains plausible toinfer that modest overdensities such as the one associatedwith HDF 850.1 may simply be highlighted by the presenceof DSFGs, and thus, a common feature of the environmentsof such sources at z > z >
4, but perhaps are much less common; e.g., Oteo etal. 2018; Miller et al. 2018, but also see Robson et al. 2014).From these considerations, it remains unclear that the fieldselection alone is sufficient to explain the apparent excess inCO-bright galaxies at z > SUMMARY AND CONCLUSIONSWe have detected CO( J =2 →
1) emission toward three z > ∼
60 arcmin VLA COLDz survey data (see P18, R19, for acomplete description), including a new secure redshift iden-tification of the “optically-dark” source GN10 at z =5.303.Despite star-formation rates of ∼ M (cid:12) yr − , two ofthe sources (which are separated by only ∼ (cid:48) on the sky) remain undetected at rest-frame ultraviolet to optical wave-lengths, below observed-frame 3.6 µ m, due to dust obscura-tion. Molecular line scans such as COLDz thus are an idealmethod to determine the redshifts for such sources, whichwill remain challenging to study in their stellar light at leastuntil the launch of the James Webb Space Telescope ( JWST ).By carrying out a multi-wavelength analysis includingnew NOEMA and VLA observations of GN10, AzTEC-3,and HDF 850.1, we find a broad range of physical proper-ties among this CO-selected z > ∼ T dust =35–92 K), a range of a factor of ∼ M gas =2.2–7.1 × M (cid:12) ) and gas-to-dust ratios (65–215), a factor ofup to ∼
30 difference in SFR surface densities ( Σ SFR =60–1800 M (cid:12) yr − kpc − ), factors of ∼ ∼ L CII / L FIR and L CII / L CO ( − ) ratios (0.6–2.5 × − and 2400–4500, respectively), and significant differences in CO lineexcitation and the implied gas densities and kinetic tempera-tures. In particular, we find a trend that appears to suggest adecrease in L CII / L FIR with increasing CO excitation. We alsofind that the gas depletion times vary by a factor of ∼ τ dep =22–70 Myr), consistent with short starburstphases. Given the high inferred Σ SFR , we cannot rule outa contribution of heavily obscured AGN to the dust heatingand/or gas excitation in these compact systems.At the resolution of the current follow-up data, GN10appears to consist of a compact, ∼ ∼ ∼ ∼
75% of the dust luminosity, em-bedded in a more extended, ∼ ∼ × lower SFR than AzTEC-3 spread across a ∼ z > z > ∼ z ∼ T dust – z relations proposed in the litera-ture were to hold (e.g., Magnelli et al. 2014; Schreiber etal. 2018), but the level of increase in T dust with redshift is notcaptured well by these relations (which require significant ex-trapolation in redshift). We investigate potential biases dueto preferential gravitational magnification, sample selectionwavelength, and SED fitting methods, but find no clear evi-dence that the median T dust is biased towards high values dueto the heterogenous composition of the parent sample. At the6 R IECHERS ET AL .same time, it cannot be ruled out that currently employed se-lection methods miss z > L FIR towards higher T dust .Given their typically very high dust luminosities, this trendis consistent with findings for lower- z DSFG samples (e.g.,da Cunha et al. 2015; Simpson et al. 2017) without requiringany redshift evolution. This perhaps suggests that the ob-served trends in T dust could be mostly, if not entirely due tothe relatively high median L FIR (and thus, SFR) of the current z > z > z > A. GN10 LINE PARAMETERSHere we provide a table including the full line parameters from Gaussian fitting to the line profiles for GN10 (Table A1). Wealso provide a figure showing the upper limit spectrum for the HCN, HCO + , and HNC J =1 → B. GN10 MAGPHYS SED FITIn addition to
CIGALE , we have also used the
MAGPHYS code (da Cunha et al. 2015) to fit the full optical to radio wavelengthphotometry of GN10. The fit to the spectral energy distribution is shown in Fig. B2, and the resulting physical parameters areprovided in Table B2.
MAGPHYS suggests a dust luminosity L dust that is ∼
25% higher than the L IR found from MBB fitting, butconsistent within the uncertainties. It also finds a total SFR that is comparable to the SFR IR found from the MBB fit, consistentwith the expectation that dust-obscured star formation dominates the SFR of GN10. M dust is about half the value found fromthe MBB fit, and T dust is ∼
20% lower (but consistent within the uncertainties). These differences are likely due to the fact that
MAGPHYS uses multiple dust components in the fitting. In particular, T dust corresponds to a luminosity-averaged value, calculatedover multiple dust components. We adopt the values from the MBB fit in the main text to enable a more straight forward Figure A1 . VLA line limit spectrum of GN10 ( z =5.3031; histogram). Spectrum is shown at a resolution of 85 km s − (4 MHz), referenced tothe expected redshift of the HCN( J =1 →
0) line. Velocities where the peaks of the HCN, HCO + , and HNC J =1 → N =1 →
0) line.
OLD Z : D USTY S TARBURSTS AT z > Table A1 . GN10 line parameters.
Line S ν d v FWHM v I line [ µ Jy] [ km s − ] [ km s − ] [Jy km s − ]CO( J =1 →
0) 74 ±
13 687 ± − ±
60 0.054 ± J =2 →
1) 544 ±
63 512 ±
72 0 ±
30 0.295 ± J =5 →
4) 1046 ±
205 772 ± − ±
71 0.86 ± J =6 →
5) 719 ±
144 681 ± − ±
70 0.52 ± P / → P / ) (A) (24.7 ± × ± − ±
24 17.6 ± b (B) (6.0 ± × ± − ± J =1 → (512) c < + ( J =1 → (512) c < J =1 → (512) c < N =1 → (512) c < d a Velocity offset relative to CO( J =2 →
1) redshift and uncertainty from Gaussian fitting.b Summed over both components. Component (A) alone is (16.2 ± − .c Fixed to CO( J =2 →
1) line width.d We conservatively assume equal strength of the hyperfine-strucutre transitions to obtain thislimit. Assuming that the three strongest components dominate would yield a 3 σ limit of < − . Figure B2 . Spectral energy distribution of GN10 ( top; red symbols), overlaid with best-fit
MAGPHYS model (black line) and unattenuated stellarlight emission spectrum before dust reprocessing (blue line), and residuals after subtracting the best-fit model ( bottom ). comparison to other z > MAGPHYS (e.g., da Cunha et al. 2015; Dudzeviciute et al. 2020; Simpson et al.2020).
MAGPHYS suggests approximately two times the M (cid:63) found by CIGALE . We adopt the value determined using
CIGALE
IECHERS ET AL . Table B2 . GN10
MAGPHYS
SED modeling parameters.
Fit Parameter unit value a T dust K 41.7 + . − . M dust M (cid:12) + . − . L dust L (cid:12) + . − . SFR total M (cid:12) yr − + − M (cid:63) e M (cid:12) + . − . a Median values are given. Lower and upper error bars are stated as 16 th and 84 th percentiles, respectively. in the main text, as CIGALE provides a better fit to the break between 2.2 and 3.6 µ m and to the 16–24 µ m photometry. Weconsider the two values to be consistent within the expected uncertainties in determining M (cid:63) for highly-obscured z > C. z =4.0508)and GN20.2b ( z =4.0563), two member galaxies of the GN20 protocluster environment (e.g., Daddi et al. 2009b; Hodge et al.2013; Tan et al. 2014) in track sharing, leading to nearly the same amount of on source time and u − v coverage (8936 visibili-ties). The pointing was centered between the two galaxies. From circular Gaussian fitting to the visibility data, we find primarybeam-corrected 1.2 mm fluxes of (3.85 ± ± (cid:48)(cid:48) ± (cid:48)(cid:48) and 0.36 (cid:48)(cid:48) ± (cid:48)(cid:48) forGN20.2a and b, corresponding to surface areas of (1.74 ± ± , respectively. The 1.2 mm fluxes thuscorrespond to source-averaged rest-frame brightness temperatures of T b =(8.5 ± ± Σ IR =(2.6 ± ± × L (cid:12) kpc − and SFR surface densities of Σ SFR =(260 ± ± M (cid:12) yr − kpc − for GN20.2a and b,respectively. This suggests that the star-formation activity in GN20.2b is almost as intense as that in the central region of GN20( Σ SFR =(120 ± M (cid:12) yr − kpc − ; Hodge et al. 2015), while GN20.2a appears to be a more intense starburst, approaching theactivity level of “maximum starbursts”.The dust continuum emission in GN20.2a appears to be more compact than the CO( J =2 →
1) emission imaged by Hodge et al.(2013), which has an extent of (0.7 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.4 (cid:48)(cid:48) ± (cid:48)(cid:48) ). The dust and cold molecular gas emission appear to peak at the sameposition, such that the rest-frame 237 µ m luminosity associated with the intense starburst is likely dominantly emerging from theregions containing the highest-density gas. The dust emission in GN20.2b also appears to be more compact than the CO( J =2 → (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.7 (cid:48)(cid:48) ± (cid:48)(cid:48) ) as measured by Hodge et al. (2013). Interestingly, the gas and dustemission appear to be spatially offset by < (cid:48)(cid:48) , with the dust emission peaking much closer to the likely near-infrared counterpartof the dusty galaxy. This may indicate the presence of multiple galaxy components, where the brightest CO-emitting componentidentified by Hodge et al. (2013) is not the same as that dominating the dust emission and stellar light. Another, perhaps lesslikely possibility is that the dust-emitting component is not at the redshift of the GN20 protocluster. On the other hand, theCO( J =6 →
5) emission appears to peak at a position that is more consistent with the 1.2 mm dust continuum peak and the stellarlight, albeit observed at about two times lower linear spatial resolution (Tan et al. 2014). More sensitive CO observations arerequired to further investigate the nature of this offset. We adopt the de-blended photometry throughout, but we caution thatuncertainties due to deblending are significant for photometry in the latterwavelength range.
JWST will be critical to overcome these uncertainties. We also observed GN20 ( z =4.0553) as part of this project in a secondsetup. Results from these data were reported by Hodge et al. (2015). Two-dimensional Gaussian fits suggest that GN20.2a is not resolvedalong its minor axis ( < (cid:48)(cid:48) , or < (cid:48)(cid:48) (2.1 kpc), but the fit does not converge well. As such, thisestimate is considered to be a weak constraint at best. For GN20.2b, a two-dimensional Gaussian fit suggests a size of(0.42 (cid:48)(cid:48) ± (cid:48)(cid:48) ) × (0.28 (cid:48)(cid:48) ± (cid:48)(cid:48) ), or (3.0 ± × (2.0 ± . OLD Z : D USTY S TARBURSTS AT z > J =2 → × (NA)
10 kpc × (UN)
10 kpc D e c li n a t i o n ( J ) J =2 → × (NA)
10 kpc × (UN)
10 kpc
Figure C3 . Rest-frame far-infrared continuum contour maps at observed-frame 1.2 mm overlaid on a
HST /WFC3 F160W continuum imagetoward GN20.2a (top) and GN20.2b (bottom). Contours start at ± σ , and are shown in steps of 1 σ =0.35 (left; imaged using natural baselineweighting) and 0.41 mJy beam − (right; imaged using uniform weighting), respectively. In the left panels, magenta CO( J =2 →
1) contourstapered to 0.38 (cid:48)(cid:48) (GN20.2a) and 0.77 (cid:48)(cid:48) (GN20.2b; Hodge et al. 2013) resolution are shown for comparison. Contour steps are the same, where1 σ =20 and 28 µ Jy beam − for GN20.2a and b, respectively. The synthesized beam size for the 1.2 mm observations is indicated in the bottomleft corner of each panel. Facilities:
VLA, NOEMA REFERENCES
Andreani, P., Cimatti, A., Loinard, L., & Röttgering, H. 2000, A&A, 354,L1Andrews, B. H., & Thompson, T. A. 2011, ApJ, 727, 97Aravena, M., Younger, J. D., Fazio, G. G., et al. 2010, ApJ, 719, L15Asboth, V., Conley, A., Sayers, J., et al. 2016, MNRAS, 462, 1989Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ApJ, 794,135Bertoldi, F., Carilli, C., Aravena, M., et al. 2007, ApJS, 172, 132Béthermin, M., De Breuck, C., Sargent, M., & Daddi, E. 2015, A&A, 576,L9Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89Blain, A. W., Barnard, V. E., & Chapman, S. C. 2003, MNRAS, 338, 733Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002,Phys. Rep., 369, 111Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207Bolatto, A. D., Chatterjee, S., Casey, C. M., et al. 2017, arXiv e-prints, arXiv:1711.09960 [astro-ph.IM]
Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047Brisbin, D., Miettinen, O., Aravena, M., et al. 2017, A&A, 608, A15 Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125Calabrò, A., Daddi, E., Puglisi, A., et al. 2019, A&A, 623, A64Calzetti, D. 2001, PASP, 113, 1449Capak, P., Carilli, C. L., Lee, N., et al. 2008, ApJ, 681, L53Capak, P. L., Riechers, D., Scoville, N. Z., et al. 2011, Nature, 470, 233Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105Carilli, C. L., Daddi, E., Riechers, D., et al. 2010, ApJ, 714, 1407Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45Casey, C. M., Zavala, J. A., Spilker, J., et al. 2018, ApJ, 862, 77Chabrier, G. 2003, PASP, 115, 763Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772Cooke, E. A., Smail, I., Swinbank, A. M., et al. 2018, ApJ, 861, 100Coppin, K. E. K., Chapman, S. C., Smail, I., et al. 2010, MNRAS, 407,L103Cowie, L. L., Barger, A. J., Wang, W.-H., & Williams, J. P. 2009, ApJ, 697,L122Cox, P., Krips, M., Neri, R., et al. 2011, ApJ, 740, 63da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110
IECHERS ET AL . Daddi, E., Dannerbauer, H., Krips, M., et al. 2009a, ApJ, 695, L176Daddi, E., Dannerbauer, H., Stern, D., et al. 2009b, ApJ, 694, 1517Daddi, E., Dannerbauer, H., Liu, D., et al. 2015, A&A, 577, A46Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2011, MNRAS, 410,1687—. 2017, ApJ, 840, 78Dannerbauer, H., Lehnert, M. D., Lutz, D., et al. 2002, ApJ, 573, 473—. 2004, ApJ, 606, 664Dannerbauer, H., Walter, F., & Morrison, G. 2008, ApJ, 673, L127Davis, T. A., Alatalo, K., Bureau, M., et al. 2013, MNRAS, 429, 534Decarli, R., Walter, F., Carilli, C., et al. 2014, ApJ, 782, 78Dickinson, M., & GOODS Team. 2004, in Bulletin of the AmericanAstronomical Society, Vol. 36, American Astronomical Society MeetingAbstracts, 1614Dowell, C. D., Conley, A., Glenn, J., et al. 2014, ApJ, 780, 75Downes, D., Neri, R., Greve, A., et al. 1999, A&A, 347, 809Dudzeviˇci¯ut˙e, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, arXiv:1910.07524 [astro-ph.GA]
Duivenvoorden, S., Oliver, S., Scudder, J. M., et al. 2018, MNRAS, 477,1099Dunne, L., Eales, S. A., & Edmunds, M. G. 2003, MNRAS, 341, 589Engel, H., Tacconi, L. J., Davies, R. I., et al. 2010, ApJ, 724, 233Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601Flower, D. R. 2001, Journal of Physics B: Atomic, Molecular and OpticalPhysics, 34, 2731Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP,125, 306Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152Frayer, D. T., Reddy, N. A., Armus, L., et al. 2004, AJ, 127, 728Fudamoto, Y., Ivison, R. J., Oteo, I., et al. 2017, MNRAS, 472, 2028Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600,L93Graciá-Carpio, J., Sturm, E., Hailey-Dunsheath, S., et al. 2011, ApJ, 728,L7Greve, T. R., Bertoldi, F., Smail, I., et al. 2005, MNRAS, 359, 1165Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35Gullberg, B., Smail, I., Swinbank, A. M., et al. 2019, MNRAS, 490, 4956Güsten, R., Philipp, S. D., Weiß, A., & Klein, B. 2006, A&A, 454, L115Hodge, J. A., Carilli, C. L., Walter, F., Daddi, E., & Riechers, D. 2013, ApJ,776, 22Hodge, J. A., Riechers, D., Decarli, R., et al. 2015, ApJ, 798, L18Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103Hodge, J. A., Smail, I., Walter, F., et al. 2019, ApJ, 876, 130Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Nature, 394, 241Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913Ivison, R. J., Lewis, A. J. R., Weiss, A., et al. 2016, ApJ, 832, 78Jin, S., Daddi, E., Magdis, G. E., et al. 2019, ApJ, 887, 144Kennicutt, Robert C., J. 1998, ApJ, 498, 541Laird, E. S., Nandra, K., Pope, A., & Scott, D. 2010, MNRAS, 401, 2763Liu, D., Gao, Y., Isaak, K., et al. 2015, ApJ, 810, L14Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172Magdis, G. E., Daddi, E., Elbaz, D., et al. 2011, ApJ, 740, L15Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86Magnelli, B., Karim, A., Staguhn, J., et al. 2019, ApJ, 877, 45Miettinen, O., Delvecchio, I., Smolˇci´c, V., et al. 2017, A&A, 606, A17Miller, T. B., Chapman, S. C., Aravena, M., et al. 2018, Nature, 556, 469Neri, R., Downes, D., Cox, P., & Walter, F. 2014, A&A, 562, A35Nguyen, H. T., Schulz, B., Levenson, L., et al. 2010, A&A, 518, L5Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793Oteo, I., Zhang, Z.-Y., Yang, C., et al. 2017a, ApJ, 850, 170Oteo, I., Ivison, R. J., Negrello, M., et al. 2017b, arXiv e-prints, arXiv:1709.04191arXiv:1709.04191