Spitzer Catalog of Herschel-selected Ultrared Dusty Star-forming Galaxies

Jingzhe Ma; Asantha Cooray; Hooshang Nayyeri; Arianna Brown; Noah Ghotbi; Rob Ivison; Ivan Oteo; Steven Duivenvoorden; Joshua Greenslade; David Clements; Julie Wardlow; Andrew Battisti; Elisabete da Cunha; Matthew L. N. Ashby; Ismael Perez-Fournon; Dominik Riechers; Seb Oliver; Stephen Eales; Mattia Negrello; Simon Dye; Loretta Dunne; Alain Omont; Douglas Scott; Pierre Cox; Stephen Serjeant; Steve Maddox; Elisabetta Valiante

doi:10.3847/1538-4365/ab4194

1. Introduction

It has become clear that observing at UV/optical wavelengths is insufficient to probe the total star formation history of the universe as a large fraction of star formation is obscured by dust (e.g., Madau & Dickinson 2014; Gruppioni et al. 2017). Wide-area infrared (IR) surveys have revolutionized our understanding of obscured star formation by discovering a large number of dusty, star-forming galaxies (DSFGs; also known as "submillimeter galaxies" or SMGs; see Casey et al. 2014 for a review), which make up the bulk of the cosmic infrared background (e.g., Dole et al. 2006). Some of the DSFGs represent the rarest and most extreme starbursts at high redshift (with star formation rates, SFRs > 10³ M_⊙ yr⁻¹, and number densities < 10⁻⁴ Mpc⁻³; Gruppioni et al. 2013), which still pose challenges to galaxy formation and evolution models (e.g., Baugh et al. 2005; Narayanan et al. 2010; Hayward et al. 2013; Béthermin et al. 2017). Discovery and detailed characterization of this population are required to understand the most extreme obscured star formation, which is only made possible now by deep and large-area surveys at far-infrared (FIR) and submillimeter/millimeter wavelengths with, e.g., the Herschel Space Observatory (Eales et al. 2010; Pilbratt et al. 2010; Oliver et al. 2012), the South Pole Telescope (SPT; Vieira et al. 2010; Carlstrom et al. 2011; Mocanu et al. 2013), the Atacama Cosmology Telescope (ACT; Marriage et al. 2011; Marsden et al. 2014; Gralla et al. 2019), and the Planck satellite (Planck Collaboration et al. 2011, 2014; Cañameras et al. 2015; Planck Collaboration et al. 2016).

The largest Herschel extragalactic surveys, the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010) and Herschel Multitiered Extragalactic Survey (HerMES; Oliver et al. 2012), covering a total area of ∼1300 deg², have revealed a large number of DSFGs. While most of them are z ∼ 1–2 starburst galaxies (e.g., Casey et al. 2012a, 2012b), selecting those with ultrared colors is extremely efficient for identifying a tail extending toward higher redshift (z > 4). A well-defined population of "ultrared" DSFGs using the Herschel SPIRE bands S₅₀₀ > S₃₅₀ > S₂₅₀ ("500 μm risers") has been established (e.g., Cox et al. 2011; Dowell et al. 2014; Asboth et al. 2016; Ivison et al. 2016). Based on their colors, these are likely to be z ≳ 4, dusty, and rapidly star-forming (>500 M_⊙ yr⁻¹) galaxies. These systems are believed to be the progenitors of massive elliptical (red and dead) galaxies identified at z ∼ 3 (e.g., Oteo et al. 2016a). Spectroscopic confirmation of a subsample of 26 sources based on CO rotational lines, an indicator of the molecular gas that fuels the prodigious star formation in these galaxies, has verified the higher redshifts compared to general DSFG samples (e.g., Cox et al. 2011; Combes et al. 2012; Riechers et al. 2013, 2017; Fudamoto et al. 2017; Donevski et al. 2018; Pavesi et al. 2018; Zavala et al. 2018). Meanwhile, relatively wide and shallow surveys with SPT have discovered a large number of gravitationally lensed DSFGs at z > 4 (Vieira et al. 2013; Spilker et al. 2016; Strandet et al. 2016), including the highest-redshift DSFG discovered so far at z = 6.9 (Strandet et al. 2017; Marrone et al. 2018).

To fully characterize the Herschel-selected ultrared DSFGs, we have been conducting a multiwavelength observational campaign to probe as many regimes of their spectral energy distributions (SEDs) as possible. Ivison et al. (2016) searched the 600 deg² H-ATLAS survey and initially selected 7961 high-redshift DSFG candidates. A subset of 109 DSFGs were further selected for follow-up observations with SCUBA-2 (Holland et al. 2013) or LABOCA (Siringo et al. 2009) at longer wavelengths (850/870 μm). Asboth et al. (2016) identified 477 "500 μm risers" from the 300 deg² HerMES Large Mode Survey (HeLMS) and 188 of the 200 brightest (S₅₀₀ > 63 mJy) sources were followed up with SCUBA-2 (Duivenvoorden et al. 2018). The addition of these longer wavelength data to the three Herschel/SPIRE bands better constrains the photometric redshifts and FIR luminosities.

The nature of the Herschel-selected ultrared sources, which can be gravitationally lensed sources, blends (multiple components including mergers at the same redshift or just projection effect from sources at different redshifts), or unlensed intrinsically bright DSFGs, requires high-resolution data to confirm. The ultrared sample provides a good opportunity for studying the gas, dust, and stellar properties of starburst galaxies at z ≳ 4 in detail out to the epoch of reionization at multiple wavelengths. In particular, ultrared sources that are not lensed, blended, or otherwise boosted are of great interest, because they may represent the most extreme galaxies in the early universe.

Gas and dust properties of DSFGs have been routinely studied with high-resolution interferometers such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the Northern Extended Millimeter Array (NOEMA). An observational campaign is being conducted with ALMA and NOEMA on a subsample (63 so far) of Herschel ultrared sources that have SCUBA-2/LABOCA data to further pinpoint their locations, reveal their morphologies, and confirm their redshifts (Fudamoto et al. 2017; Oteo et al. 2017).

Observations of the stellar populations at the optical or near-IR (NIR) are needed in order to place this population in the context of galaxy formation and evolution and provide a complete picture of their physical properties. Spitzer/IRAC is the only currently available facility that probes the rest-frame optical stellar emission of these sources and the only way to constrain their stellar masses and thus their specific SFRs (sSFRs), which are a critical diagnostic for the star formation mode of these galaxies. In this work, we present a follow-up study of 300 Herschel ultrared sources from H-ATLAS and HeLMS with Spitzer/IRAC in combination with multiwavelength ancillary data. The Spitzer data allow us to constrain the stellar masses of a statistical sample of DSFGs at z > 4 in a consistent manner for the first time and provide the first constraint on the stellar mass density and evolution of this population.

The paper is organized as follows. We describe the Spitzer/IRAC observations, source detection, and photometry in Section 2 along with the multiwavelength ancillary data. In Section 3, we introduce the cross-identification methods and describe the photometric catalog of the cross-matched Spitzer counterparts. We perform panchromatic SED modeling to derive their physical properties. The results of the SED fitting are presented in Section 4. We discuss the redshift distribution, multiplicity, unlensed fraction, SFR surface density, space density, SFR density (SFRD), and stellar mass density in Section 5. Section 6 summarizes our conclusions and future plans.

Throughout this paper, we adopt a concordance ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7, and Ω_m = 0.3. We use a Chabrier (2003) initial mass function (IMF) and AB system magnitudes.

2. Sample Selection, Spitzer Observations, and Ancillary Data

2.1. Herschel/SPIRE at 250, 350, and 500 μm

The H-ATLAS observations were performed in the parallel mode with both SPIRE (Griffin et al. 2010) and PACS (Poglitsch et al. 2010). The survey covers three equatorial fields at right ascensions (R.A.) of 9, 12, and 15 hr (namely the GAMA09, GAMA12, and GAMA15 fields), the North Galactic Pole (NGP) field in the north, and the South Galactic Pole (SGP) field in the south. The total survey area is about 600 deg². For SPIRE maps, sources were extracted and flux densities were measured based on a matched-filter approach, which mitigates the effects of confusion (e.g., Chapin et al. 2011). The depth of the PACS data (1σ ∼ 45 mJy) is insufficient to detect our FIR-rising sources. The detailed descriptions of the H-ATLAS observations and source extraction can be found in Valiante et al. (2016), Ivison et al. (2016), and Maddox et al. (2018). Ivison et al. (2016) selected the ultrared sample based on the following criteria: 3.5σ detection threshold at S₅₀₀ > 30 mJy and the color selection with S₅₀₀/S₂₅₀ ≥ 1.5 and S₅₀₀/S₃₅₀ ≥ 0.85.

The HeLMS field covers an effective area of ∼274 deg² of the equatorial region spanning 23^h14^m < R.A. < 1^h16^m and −9° < decl. < +9°. The observations were conducted with the SPIRE instrument (Griffin et al. 2010). Sources were extracted using a map-based search method (Dowell et al. 2014; Asboth et al. 2016) that combines the information in the 250, 350, and 500 μm maps simultaneously. Herschel/SPIRE maps have pixel scales of 6'', 8'' (8 farcs 333 for HeLMS maps), and 12'' at 250, 350, and 500 μm, which are one-third of the FWHM beam sizes of 18'', 25'', and 36'', respectively. We refer the reader to Asboth et al. (2016) and Duivenvoorden et al. (2018) for more information on the HeLMS observations and detailed description of the source extraction and photometry. The HeLMS ultrared sample was selected with flux densities S₅₀₀ > S₃₅₀ > S₂₅₀ and a 5σ cutoff S₅₀₀ > 52 mJy (Asboth et al. 2016).

Figure 1 demonstrates the distributions of the 300 Herschel ultrared sources on the all-sky map for which we have obtained Spitzer data. These include all of the H-ATLAS sources in Ivison et al. (2016) and an additional 31 sources from R. Ivison et al. (2019, in preparation), and 161 HeLMS sources (with a flux density cut S₅₀₀ > 64 mJy) from Asboth et al. (2016). The Herschel/SPIRE photometry of these sources at 250, 350, and 500 μm are shown in Table 1.

Table 1. Catalog of Herschel-selected Ultrared Sources

Column	Parameter	Description	Units
1	IAU_name	Survey name + Herschel source name
2	ID	Nickname or short name in Ivison et al. (2016) and Duivenvoorden et al. (2018)
3	z	Photometric redshift from FIR SED fitting (two decimal points; Ivison et al. 2016 and
		Duivenvoorden et al. 2018) or spectroscopic redshift (three decimal points; see Table 3)
4	Lensed?	Lensed (y), unlensed (n), yn (both lensed and unlensed components), unknown (x)
5	RA_IRAC	Spitzer/IRAC counterpart position: R.A. (J2000)	degree
6	Dec_IRAC	Spitzer/IRAC counterpart position: decl. (J2000)	degree
7	MAG_AUTO_3.6μm	Kron-like elliptical aperture magnitude at 3.6 μm	AB mag
8	MAGERR_AUTO_3.6μm	AUTO magnitude uncertainty at 3.6 μm	AB mag
9	MAG_AUTO_4.5μm	Kron-like elliptical aperture magnitude at 4.5 μm	AB mag
10	MAGERR_AUTO_4.5μm	AUTO magnitude uncertainty at 4.5 μm	AB mag
11	S_3.6μm	Flux density at 3.6 μm	Jy
12	S_3.6μm_err	Flux density uncertainty at 3.6 μm	Jy
13	S_4.5μm	Flux density at 4.5 μm	Jy
14	S_4.5μm_err	Flux density uncertainty at 4.5 μm	Jy
15	RA_H	Herschel/SPIRE position: R.A. (J2000)	degree
16	Dec_H	Herschel/SPIRE position: decl. (J2000)	degree
17	S_250μm	Flux density at 250 μm	Jy
18	S_250μm_err	Flux density uncertainty at 250 μm	Jy
19	S_350μm	Flux density at 350 μm	Jy
20	S_350μm_err	Flux density uncertainty at 350 μm	Jy
21	S_500μm	Flux density at 500 μm	Jy
22	S_500μm_err	Flux density uncertainty at 500 μm	Jy
23	RA_S	SCUBA-2 position: R.A. (J2000)	degree
24	Dec_S	SCUBA-2 position: decl. (J2000)	degree
25	S_850μm	Flux density at 850 μm (SCUBA-2) or at 870 μm (LABOCA)	Jy
26	S_850μm_err	Flux density uncertainty at 850 μm (SCUBA-2) or at 870 μm (LABOCA)	Jy

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

2.2. Spitzer Observations and Source Photometry

A total of 300 Herschel ultrared sources were followed up by the Spitzer snapshot imaging program PID 13042 (PI: A. Cooray) using IRAC (Fazio et al. 2004). Images were taken at 3.6 and 4.5 μm with 30 s exposure per frame and a 36 position dither pattern in each band for each source (one source per astronomical observation request (AOR)), totaling 1080 s integration per band. We obtained the post-basic calibrated data (pBCDs) from the Spitzer Science Center (pipeline version S19.2), which have been reduced with the Mosaicing and Point-source Extraction (mopex¹⁷ ) package, including background matching, overlap correction, and mosaicking. The pBCDs, in most cases, are of good quality for our purpose of source detection and photometry. When necessary, we also downloaded and reprocessed the artifact-corrected BCDs with MOPEX (version 18.5.0; Makovoz et al. 2006) to generate improved final mosaics. The IRAC mosaics (one mosaic per source/AOR) have a resampled pixel scale of 0 farcs 6 pixel⁻¹ and an angular resolution of ∼1 farcs 9.

We use SExtractor (version 2.8.6; Bertin & Arnouts 1996) in dual-image mode to perform source detection and photometry. SExtractor is well suited for this task, because the relatively sparse mosaics have sufficient source-free pixels available for robust sky background estimation. We use the coadded 3.6 and 4.5 μm image as the detection image. The dual-image approach ensures that the photometry is measured in identical areas in both bands, yielding accurate source colors. We obtained the Kron-like elliptical aperture magnitudes, MAG_AUTO, and the corresponding flux densities for the following analyses. To estimate the accuracy of the astrometry, we compared the positions of bright IRAC sources with their 2MASS counterparts (Skrutskie et al. 2006). The astrometric discrepancy is small, about one-fourth (0 farcs 15) of the mosaic pixel.

2.3. SCUBA-2/LABOCA at 850/870 μm

For the H-ATLAS ultrared DSFGs, Ivison et al. (2016) obtained 850 μm continuum imaging with SCUBA-2 and/or 870 μm continuum imaging with LABOCA. The FWHM of the main beam of SCUBA-2 is 13 farcs 0 at 850 μm, and the astrometry accuracy σ is 2''–3''. LABOCA has an FWHM resolution of 19 farcs 2, and the positional uncertainty is estimated to be about 1''–2''. They measured the 850 or 870 μm flux densities via several methods (brightest pixel values and aperture photometry) and performed template fitting together with the three Herschel/SPIRE bands to derive photometric redshifts and FIR luminosities. Duivenvoorden et al. (2018) presented the SCUBA-2 observations of the HeLMS sources and extracted flux densities by taking the brightest pixel values within the 3σ of the positional uncertainty in the SCUBA-2 map. Photometric redshifts and integrated properties (e.g., FIR luminosities) are derived from the eazy code (Brammer et al. 2008) using representative FIR/submillimeter templates. We refer the reader to Ivison et al. (2016) and Duivenvoorden et al. (2018) for detailed descriptions of the SCUBA-2/LABOCA observations, flux density measurements, and FIR SED fitting. We adopt the peak-value photometry and photometric redshifts from both works (listed in Table 1). A total of 261 out of the 300 Spitzer follow-up sources have SCUBA-2/LABOCA data, which are used in the following SED modeling.

2.4. High-resolution Submillimeter/millimeter Data

2.4.1. Continuum Data at 870 μm, 1.1, 1.3, and 3 mm

A subset (21) of the H-ATLAS ultrared sources in Ivison et al. (2016) was selected for continuum imaging with NOEMA (PIs: R.J. Ivison, M. Krips; PIDs: W05A, X0C6, W15ET) at 1.3 and/or 3 mm and ALMA (PI: A. Conley; PID: 2013.1.00499.S) at 3 mm as described in detail in Fudamoto et al. (2017). They observed 17 with NOEMA and four with ALMA, based on their accessibility and high photometric redshifts. The synthesized FWHM beam sizes are 1''–1 farcs 5 for NOEMA and 0 farcs 6–1 farcs 2 for ALMA. Precise positions were determined from the continuum images for 18 ultrared sources, and the remaining three sources lack secure detection.

Oteo et al. (2017) presented high-resolution (∼0 farcs 12) ALMA continuum imaging at 870 μm (PI: R.J. Ivison; PIDs: 2013.1.00001.S, 2016.1.00139.S) for a sample of 44 equatorial and southern ultrared sources from both the H-ATLAS and HerMES surveys. Thirty-one of them are in our Spitzer sample and thus included in our following analysis. An additional 18 H-ATLAS ultrared sources whose ALMA data were released after Oteo et al. (2017) are also included in this work.

Five HeLMS ultrared sources (HELMS_RED_1, 2, 4, 10, 13) were observed with SMA at 1.1 mm in the compact array configuration (PI: D. Clements; PID: 2013A-S005). The reduced maps have an average rms noise level of 2.2 mJy beam⁻¹ and the beam FWHM sizes are typically 2 farcs 5. The flux densities at 1.1 mm are given in Duivenvoorden et al. (2018) and included in the SED fitting below. The details of the observations and data reduction will be presented in J. Greenslade et al. (2019, in preparation). Additional data with MUSIC/CSO (Sayers et al. 2014) and ACT (Su et al. 2017) were obtained for five HeLMS sources (HELMS_RED_1, 3, 4, 6, 7). We list the sources with submillimeter/millimeter flux density measurements from SMA, MUSIC, and ACT in the Appendix B (Table 6). Duivenvoorden et al. (2018) summarized the observations obtained so far as part of the still ongoing observational campaign.

We analyze in this paper a total of 63 Herschel-selected ultrared sources from H-ATLAS and HeLMS that have high-resolution positions from various observations as described above. The positions and flux densities of the high-resolution subsample are listed in Table 2. The original 63 ultrared sources are resolved into 86 individual submillimeter/millimeter sources as seen at high resolution (see further discussion in Sections 3 and 5.2). These sources are classified as lensed or unlensed based on the submillimeter/millimeter morphology and the presence or absence of low-redshift foreground galaxies (further discussion in Sections 3 and 5.3).

Table 2. High-resolution Sample Followed-up with ALMA, NOEMA, and SMA

Column	Parameter	Description	Units
1	ID	Source name
2	Highres_data	High-resolution data from ALMA (1), NOEMA (2), SMA (3), etc.
3	Lensed?	Lensed (y), unlensed (n), or unknown (x) according to the high-resolution interferometry data
4	RA_IRAC	Spitzer/IRAC counterpart position: R.A. (J2000)	degree
5	Dec_IRAC	Spitzer/IRAC counterpart position: decl. (J2000)	degree
6	MAG_AUTO_3.6μm	Kron-like elliptical aperture magnitude at 3.6 μm	AB mag
7	MAGERR_AUTO_3.6μm	AUTO magnitude uncertainty at 3.6 μm	AB mag
8	MAG_AUTO_4.5μm	Kron-like elliptical aperture magnitude at 4.5 μm	AB mag
9	MAGERR_AUTO_4.5μm	AUTO magnitude uncertainty at 4.5 μm	AB mag
10	S_3.6μm	Flux density at 3.6 μm	Jy
11	S_3.6μm_err	Flux density uncertainty at 3.6 μm	Jy
12	S_4.5μm	Flux density at 4.5 μm	Jy
13	S_4.5μm_err	Flux density uncertainty at 4.5 μm	Jy
14	S_250μm	Flux density at 250 μm; deblended if multiple components	Jy
15	S_250μm_err	Flux density uncertainty at 250 μm	Jy
16	S_350μm	Flux density at 350 μm; deblended if multiple components	Jy
17	S_350μm_err	Flux density uncertainty at 350 μm	Jy
18	S_500μm	Flux density at 500 μm; deblended if multiple components	Jy
19	S_500μm_err	Flux density uncertainty at 500 μm	Jy
20	S_850μm	Flux density at 850 μm (SCUBA-2) or at 870 μm (LABOCA); deblended	Jy
		if multiple components
21	S_850μm_err	Flux density uncertainty at 850 μm (SCUBA-2) or at 870 μm (LABOCA)	Jy

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

2.4.2. Spectroscopic Redshifts from Spectral Scans

Fudamoto et al. (2017) conducted spectral scans at the 3 mm atmospheric window for 21 H-ATLAS ultrared sources with NOEMA (co-PIs: R.J. Ivison, M. Krips; PIDs: W05A, X0C6) and ALMA (PI: A. Conley; PID: 2013.1.00499.S). They obtained eight secure redshifts via detections of multiple CO lines and three redshifts via a single CO line detection. One of the SGP sources, SGP-354388, was confirmed to be the core of a protocluster that lies at z = 4.002 via detections of CO, [C i], and H₂O lines with ALMA and ATCA (Oteo et al. 2018). Four HeLMS ultrared sources have spectroscopic redshifts confirmed with ALMA and CARMA (Asboth et al. 2016; Duivenvoorden et al. 2018). Table 3 lists all of the spectroscopically confirmed Herschel-selected ultrared DSFGs at z ≳ 4 in this work or in the literature.

Table 3. Spectroscopically Confirmed Herschel-selected Ultrared DSFGs at z ≳ 4

IAU name	Nickname	z_spec	L_IR	Lensed?	Reference
			(10¹³ L_⊙)
HATLAS J084937.0+001455	G09-81106	4.531	2.7	No	This work; (a), (g)
HATLAS J090045.4+004125	G09-83808	6.027	3.2	Yes	This work; (a)
HATLAS J133337.6+241541	NGP-190387	4.420	3.1	Yes	This work; (a)
HATLAS J134114.2+335934	NGP-246114	3.847	2.0	No	This work; (a)
HATLAS J133251.5+332339	NGP-284357	4.894	2.5	No	This work; (a)
HATLAS J000306.9-330248	SGP-196076	4.425	2.5	No	This work; (a), (b)
HATLAS J000607.6-322639	SGP-261206	4.242	4.4	Yes	This work; (a)
HATLAS J004223.5-334340	SGP-354388	4.002	4.8	No	This work; (c)
HerMES J004409.9+011823	HELMS_RED_1	4.163	8.9	Yes	This work; (d)
HerMES J005258.9+061319	HELMS_RED_2	4.373	8.3	Yes	This work; (d)
HerMES J002220.8-015521	HELMS_RED_4	5.161	5.8	Yes	This work; (e)
HerMES J002737.4-020801	HELMS_RED_31	3.798	3.4	Yes	This work; (e)

HATLAS J142413.9+022304	ID141	4.243	8.5	Yes	(h)
HerMES J043657.5-543809	ADFS-27	5.655	2.4	No	(f)
HLS J091828.6+514223	⋯	5.243	11	Yes	(i)
HerMES J104050.6+560654	LSW102	5.29	⋯	Yes	(j), (k), (m)
HerMES J170647.7+584623	HFLS3	6.337	2.9	Yes	(k), (l)
HerMES J170817.3+582844	HFLS1	4.286	5.6	No	(j), (k)
HerMES J172049.0+594623	HFLS5	4.44	2.8	No	(j), (k)

References. (a) Fudamoto et al. (2017), (b) Oteo et al. (2016a), (c) Oteo et al. (2018), (d) Duivenvoorden et al. (2018), (e) Asboth et al. (2016), (f) Riechers et al. (2017), (g) Zavala et al. (2018), (h) Cox et al. (2011), (i) Combes et al. (2012), (j) Dowell et al. (2014), (k) Riechers et al. (2013), (l) Cooray et al. (2014), (m). J. Wardlow et al. (2019, in preparation)

A machine-readable version of the table is available.

Download table as: Data Typeset image

3. Cross-identification and Spitzer Catalogs

Before we cross-identify the Spitzer counterparts to the Herschel ultrared sources, we cross-match the IRAC data with shallow and medium-depth large-area surveys at optical/NIR wavelengths, e.g., the Sloan Digital Sky Survey (SDSS; York et al. 2000) Data Release 12 (DR12; Alam et al. 2015) and the VISTA Kilo-degree Infrared Galaxy (VIKING) Survey (Edge et al. 2013), to identify low-redshift galaxies that can potentially magnify a higher redshift ultrared source via gravitational lensing. Figure 10 in the Appendix A shows the 60'' × 60'' Spitzer/IRAC cutouts centered on the Herschel positions for all of the ultrared sources in our Spitzer program. The 3'' radius cyan circles denote the positions of the SDSS sources in the field, and the yellow circles show the positions of the VIKING sources.

For the ultrared sources with high-resolution submillimeter/millimeter detections in Table 2, we can use these data to pinpoint the locations and securely identify the Spitzer/IRAC counterparts. Within this high-resolution subsample of 63 original ultrared sources, 23 of them are gravitationally lensed sources with clear lensing signatures like rings or arcs as shown in Oteo et al. (2017), and the rest of them that do not show any lensing features are likely unlensed, intrinsically bright ultrared DSFGs (more discussion on the lensed/unlensed fraction in Section 5.3). For sources that are classified as unlensed at submillimeter and successfully cross-identified in Spitzer/IRAC, we can directly extract the Spitzer/IRAC source flux densities. We list the flux densities at 3.6 and 4.5 μm as well as deblended SPIRE, SCUBA-2/LABOCA, and ALMA flux densities of the unlensed sample in Table 4, which will be used in multiwavelength SED fitting to derive their physical properties. The unlensed sample is the focus of our SED analysis in Section 4.

Table 4. Unlensed Sample

Source name	z_FIR	z_MAGPHYS	S_3.6μm	S_4.5μm	S_250μm	S_350μm	S_500μm	S_850μm^{SCUBA − 2}	S_870μm^ALMA
			(μJy)	(μJy)	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)
G09-47693	${3.12}_{-0.33}^{+0.39}$	2.73 ${}_{-0.71}^{+1.01}$	35.21 ± 9.06	45.30 ± 10.24	27.4 ± 7.3	34.4 ± 8.1	45.4 ± 8.6	12.5 ± 4.0	5.5 ± 0.4
G09-47693a	⋯	2.51 ${}_{-0.69}^{+0.91}$	13.66 ± 5.64	17.05 ± 6.28	11.9 ± 8.1	20.8 ± 9.0	38.0 ± 10.1	5.1 ± 1.6	2.24 ± 0.21
G09-47693b	⋯	2.79 ${}_{-0.96}^{+1.33}$	21.55 ± 7.09	28.25 ± 8.09	11.3 ± 7.9	8.1 ± 6.0	7.9 ± 6.0	7.4 ± 2.4	3.25 ± 0.29
G09-51190	${3.83}_{-0.48}^{+0.58}$	${3.31}_{-0.81}^{+1.32}$	28.96 ± 8.25	33.18 ± 8.83	28.5 ± 7.6	39.5 ± 8.1	46.6 ± 8.6	28.3 ± 7.3	9.7 ± 0.5
G09-51190a	⋯	${2.73}_{-0.69}^{+1.02}$	28.96 ± 8.25	33.18 ± 8.83	21.9 ± 6.4	27.6 ± 7.8	41.1 ± 7.3	10.4 ± 2.9	3.55 ± 0.21
G09-51190b	⋯	${4.69}_{-1.54}^{+1.58}$	<7.84	<6.50	11.4 ± 6.2	14.4 ± 7.7	6.2 ± 5.9	17.9 ± 4.6	6.15 ± 0.48
G09-59393	${3.70}_{-0.26}^{+0.35}$	${3.42}_{-0.86}^{+1.20}$	<15.07	10.19 ± 5.08	24.1 ± 7.0	43.8 ± 8.3	46.8 ± 8.6	23.7 ± 3.5	12.4 ± 0.4
G09-59393a	⋯	${3.15}_{-0.87}^{+1.01}$	<7.30	10.19 ± 5.08	23.7 ± 4.8	41.7 ± 5.2	36.3 ± 6.8	16.4 ± 2.4	7.33 ± 0.48
G09-59393b	⋯	${5.07}_{-1.28}^{+1.43}$	<7.77	<8.55	3.6 ± 3.3	2.7 ± 2.9	6.1 ± 5.4	7.3 ± 1.1	3.25 ± 0.64
G09-62610	${3.70}_{-0.26}^{+0.44}$	${3.30}_{-0.85}^{+1.23}$	12.98 ± 5.55	12.51 ± 5.51	18.6 ± 5.4	37.3 ± 7.4	44.3 ± 7.8	19.5 ± 4.9	13.6 ± 0.7
G09-62610a	⋯	${3.11}_{-0.84}^{+1.08}$	2.51 ± 2.46	2.47 ± 2.44	9.1 ± 4.4	18.0 ± 6.1	23.8 ± 8.5	6.2 ± 1.6	4.31 ± 0.35
G09-62610b	⋯	${4.30}_{-1.21}^{+1.47}$	<7.02	3.75 ± 3.05	3.6 ± 3.4	7.8 ± 6.9	9.5 ± 8.5	8.0 ± 2.0	5.57 ± 0.28
G09-62610c	⋯	3.55 ${}_{-1.01}^{+1.45}$	10.47 ± 4.97	6.29 ± 3.89	3.8 ± 3.4	8.6 ± 6.4	8.0 ± 7.5	5.3 ± 1.3	3.72 ± 0.50
G09-64889	${3.48}_{-0.40}^{+0.48}$	${3.12}_{-0.81}^{+1.18}$	11.24 ± 5.21	10.57 ± 5.09	20.2 ± 5.9	30.4 ± 7.7	34.7 ± 8.1	15.1 ± 4.3	7.91 ± 0.35
G09-79552	${3.59}_{-0.26}^{+0.34}$	${3.24}_{-0.85}^{+1.14}$	7.37 ± 4.16	10.47 ± 4.92	16.6 ± 6.2	38.1 ± 8.1	42.8 ± 8.5	17.0 ± 3.6	12.7 ± 0.6
G09-80620	${4.01}_{-0.78}^{+0.22}$	${3.17}_{-0.83}^{+1.16}$	11.97 ± 5.50	19.46 ± 6.92	13.5 ± 5.0	25.3 ± 7.4	28.4 ± 7.7	13.2 ± 4.3	8.4 ± 0.7
G09-80620a	⋯	${1.99}_{-0.67}^{+1.10}$	7.85 ± 4.43	14.00 ± 5.83	4.5 ± 3.9	7.3 ± 5.8	7.5 ± 6.8	10.3 ± 3.4	5.08 ± 0.97
G09-80620b	⋯	${2.88}_{-0.85}^{+1.19}$	4.12 ± 3.26	5.46 ± 3.73	5.5 ± 4.5	6.7 ± 5.7	12.8 ± 7.7	2.9 ± 0.9	1.45 ± 0.29
G09-80658	${4.07}_{-0.72}^{+0.09}$	${3.19}_{-0.79}^{+1.14}$	20.00 ± 6.85	31.62 ± 8.57	17.8 ± 6.4	31.6 ± 8.3	39.5 ± 8.8	17.6 ± 4.1	10.7 ± 0.7
G09-80658a	⋯	${3.75}_{-0.92}^{+1.36}$	20.00 ± 6.85	31.62 ± 8.57	5.9 ± 4.7	12.0 ± 7.2	20.3 ± 10.2	13.2 ± 3.1	4.08 ± 0.40
G09-80658b	⋯	${4.26}_{-1.17}^{+1.45}$	<5.94	<10.02	3.0 ± 2.3	3.4 ± 2.5	8.0 ± 5.5	4.4 ± 1.0	1.35 ± 0.13
G09-81106	4.531	${4.23}_{-1.06}^{+1.32}$	3.22 ± 2.81	4.68 ± 3.35	14.0 ± 6.0	30.9 ± 8.2	47.5 ± 8.8	30.2 ± 5.2	28.4 ± 0.8
G09-81271	${4.62}_{-0.38}^{+0.46}$	${4.32}_{-1.02}^{+1.39}$	4.52 ± 3.36	3.54 ± 3.00	15.0 ± 6.1	30.5 ± 8.2	42.3 ± 8.6	29.7 ± 3.7	20.5 ± 0.7
G09-87123	${4.28}_{-0.34}^{+0.52}$	${3.85}_{-0.97}^{+1.25}$	7.78 ± 4.28	13.36 ± 5.57	10.4 ± 5.8	25.3 ± 8.2	39.2 ± 8.7	20.7 ± 4.6	6.63 ± 0.47
G09-100369	${3.79}_{-0.46}^{+0.61}$	${3.42}_{-0.89}^{+1.23}$	3.75 ± 3.10	6.48 ± 3.95	15.4 ± 5.5	17.3 ± 7.6	32.3 ± 8.0	13.2 ± 3.6	3.74 ± 0.25
G09-101355	${4.20}_{-0.39}^{+0.70}$	${3.74}_{-1.02}^{+1.32}$	<10.86	10.22 ± 4.90	9.5 ± 5.5	14.6 ± 7.9	33.4 ± 8.3	13.5 ± 4.9	7.6 ± 0.7
G09-101355a	⋯	${3.99}_{-1.21}^{+1.50}$	<4.68	4.24 ± 3.16	3.6 ± 3.1	7.2 ± 4.4	15.6 ± 7.9	8.5 ± 3.1	4.78 ± 0.60
G09-101355b	⋯	${2.92}_{-0.86}^{+1.13}$	<6.18	5.99 ± 3.75	9.0 ± 4.2	11.1 ± 4.7	12.7 ± 8.0	5.0 ± 1.8	2.80 ± 0.42
NGP-101333	${3.53}_{-0.27}^{+0.34}$	${3.23}_{-0.80}^{+1.24}$	12.75 ± 5.26	13.11 ± 5.66	32.4 ± 7.5	46.5 ± 8.2	52.8 ± 9.0	24.6 ± 3.8	⋯
NGP-111912	${3.27}_{-0.26}^{+0.36}$	${2.81}_{-0.72}^{+1.03}$	36.92 ± 9.27	47.05 ± 10.42	25.2 ± 6.5	41.5 ± 7.6	50.2 ± 8.0	14.9 ± 3.9	⋯
NGP-136156	${3.95}_{-0.57}^{+0.06}$	${3.25}_{-0.79}^{+1.21}$	7.92 ± 4.15	10.23 ± 4.95	29.3 ± 7.4	41.9 ± 8.3	57.5 ± 9.2	23.4 ± 3.4	⋯
NGP-246114	3.847	${3.91}_{-0.94}^{+1.43}$	17.12 ± 5.80	17.88 ± 5.86	17.3 ± 6.5	30.4 ± 8.1	33.9 ± 8.5	25.9 ± 4.6	⋯
NGP-252305	4.34 ${}_{-0.38}^{+0.43}$	${3.95}_{-0.93}^{+1.35}$	6.32 ± 3.89	7.90 ± 4.33	15.3 ± 6.1	27.7 ± 8.1	40.0 ± 9.4	24.0 ± 3.5	⋯
NGP-284357	4.894	${4.51}_{-1.07}^{+1.33}$	4.62 ± 3.32	3.43 ± 2.87	12.6 ± 5.3	20.4 ± 7.8	42.4 ± 8.3	28.9 ± 4.3	⋯
SGP-72464	${3.06}_{-0.19}^{+0.21}$	${2.74}_{-0.73}^{+1.00}$	26.76 ± 8.31	28.51 ± 8.52	43.4 ± 7.6	67.0 ± 8.0	72.6 ± 8.9	20.0 ± 4.2	16.9 ± 0.4
SGP-93302b	${3.91}_{-0.22}^{+0.27}$	${2.84}_{-0.79}^{+0.90}$	<8.13	4.52 ± 3.66	18.3 ± 8.1	33.4 ± 14.9	19.5 ± 16.3	9.4 ± 0.9	9.95 ± 0.90
SGP-135338	${3.06}_{-0.26}^{+0.33}$	${3.14}_{-0.84}^{+0.94}$	<3.30	<6.94	32.9 ± 7.3	43.6 ± 8.1	53.3 ± 8.8	14.7 ± 3.8	6.1 ± 0.4
SGP-196076	4.425	${4.17}_{-0.96}^{+1.36}$	12.66 ± 5.19	10.17 ± 4.91	28.6 ± 7.3	28.6 ± 8.2	46.2 ± 8.6	32.5 ± 4.1	34.6 ± 2.3
SGP-196076a	4.425	${4.92}_{-1.19}^{+1.40}$	5.94 ± 3.11	5.93 ± 1.81	8.6 ± 6.0	11.2 ± 9.4	16.4 ± 12.2	21.3 ± 2.7	17.58 ± 1.03
SGP-196076b	4.425	${4.35}_{-1.24}^{+1.44}$	6.72 ± 2.08	4.24 ± 1.57	9.6 ± 6.5	8.4 ± 6.1	12.0 ± 8.6	9.6 ± 1.2	7.90 ± 0.60
SGP-196076c	4.425	${2.78}_{-0.80}^{+1.17}$	<6.24	<6.34	5.1 ± 4.0	8.2 ± 5.9	11.5 ± 8.3	1.6 ± 0.2	1.33 ± 0.17
SGP-208073	${3.48}_{-0.28}^{+0.40}$	3.10 ${}_{-0.75}^{+1.19}$	42.70 ± 9.59	50.51 ± 10.77	28.0 ± 7.4	33.2 ± 8.1	44.3 ± 8.5	19.4 ± 2.9	13.9 ± 0.9
SGP-208073a	⋯	${4.37}_{-1.13}^{+1.38}$	2.84 ± 2.10	3.19 ± 2.18	3.6 ± 2.7	9.8 ± 8.3	12.0 ± 8.6	8.7 ± 1.3	6.24 ± 0.56
SGP-208073b	⋯	3.61 ${}_{-0.94}^{+1.42}$	28.13 ± 7.76	31.03 ± 8.58	15.9 ± 9.9	9.4 ± 6.7	10.9 ± 8.1	9.6 ± 1.4	6.91 ± 0.63
SGP-208073c	⋯	1.87 ${}_{-0.52}^{+0.71}$	11.73 ± 5.23	16.29 ± 6.14	9.1 ± 8.0	12.2 ± 8.0	11.8 ± 8.2	1.2 ± 0.2	0.83 ± 0.21
SGP-213813	${3.49}_{-0.32}^{+0.40}$	2.98 ${}_{-0.71}^{+1.39}$	75.90 ± 13.25	69.53 ± 12.67	23.9 ± 6.3	35.1 ± 7.6	35.9 ± 8.2	18.1 ± 3.6	13.9 ± 0.7
SGP-219197	${2.94}_{-0.24}^{+0.25}$	${2.61}_{-0.70}^{+0.92}$	29.18 ± 8.31	36.62 ± 9.25	27.6 ± 7.4	51.3 ± 8.1	43.6 ± 8.4	12.2 ± 3.7	9.47 ± 0.40
SGP-317726	${3.69}_{-0.30}^{+0.39}$	${3.56}_{-0.92}^{+1.20}$	<10.37	<8.99	20.4 ± 6.0	35.1 ± 7.7	39.5 ± 8.0	19.4 ± 3.2	26.9 ± 2.9
SGP-354388	4.002	${4.89}_{-1.16}^{+1.45}$	5.35 ± 3.43	8.76 ± 4.58	26.6 ± 8.0	39.8 ± 8.9	53.5 ± 9.8	64.1 ± 10.9	24.1 ± 1.7
SGP-354388a	4.002	${5.60}_{-1.22}^{+1.32}$	<5.35	<8.76	6.2 ± 4.6	16.0 ± 12.0	18.8 ± 13.1	36.7 ± 6.2	9.64 ± 0.33
SGP-354388b	4.002	${4.46}_{-1.16}^{+1.40}$	<5.35	<8.76	6.3 ± 4.5	16.8 ± 11.7	19.3 ± 16.1	13.7 ± 2.3	3.61 ± 0.25
SGP-354388c	4.002	${4.67}_{-1.18}^{+1.40}$	<5.35	<8.76	6.2 ± 4.6	11.7 ± 8.4	13.6 ± 10.1	13.6 ± 2.3	3.58 ± 0.16
SGP-380990	${2.84}_{-0.21}^{+0.22}$	${2.68}_{-0.76}^{+0.92}$	4.99 ± 3.45	8.59 ± 4.48	14.4 ± 5.9	45.6 ± 8.2	40.6 ± 8.5	7.7 ± 1.8	8.93 ± 0.36
SGP-381615	${2.98}_{-0.29}^{+0.29}$	${2.97}_{-0.81}^{+0.98}$	<4.09	<5.41	19.40 ± 6.6	39.1 ± 8.1	34.7 ± 8.5	8.5 ± 3.6	6.20 ± 0.29
SGP-381637	${3.30}_{-0.25}^{+0.28}$	${3.04}_{-0.83}^{+1.02}$	3.53 ± 3.01	7.49 ± 4.27	18.7 ± 6.8	41.5 ± 8.4	49.3 ± 8.6	12.6 ± 3.7	4.45 ± 0.31
SGP-382394	${2.96}_{-0.26}^{+0.29}$	2.73 ${}_{-0.76}^{+0.98}$	9.58 ± 5.00	11.94 ± 5.46	15.7 ± 5.9	35.6 ± 8.1	35.9 ± 8.6	8.0 ± 2.4	2.11 ± 0.29
SGP-385891	${3.70}_{-0.24}^{+0.29}$	${3.32}_{-0.89}^{+1.13}$	8.44 ± 4.26	13.29 ± 5.59	13.0 ± 8.2	45.6 ± 9.8	59.6 ± 1.2	20.5 ± 3.6	11.1 ± 0.7
SGP-385891a	⋯	${5.04}_{-1.25}^{+1.41}$	<7.31	<7.81	5.3 ± 4.7	7.6 ± 5.2	20.2 ± 13.9	16.3 ± 2.9	5.9 ± 0.3
SGP-385891b	⋯	${2.43}_{-0.68}^{+0.83}$	8.44 ± 4.26	13.29 ± 5.59	12.1 ± 6.9	27.6 ± 9.0	24.7 ± 14.8	4.2 ± 0.7	1.5 ± 0.2
SGP-386447	${4.89}_{-0.73}^{+0.78}$	${4.18}_{-1.10}^{+1.35}$	14.81 ± 5.75	21.52 ± 7.24	10.5 ± 6.0	33.6 ± 8.4	34.5 ± 8.6	34.3 ± 8.4	6.5 ± 0.6
SGP-392029	${3.42}_{-0.32}^{+0.47}$	${3.70}_{-1.27}^{+1.67}$	<8.01	3.29 ± 2.91	18.3 ± 6.5	30.5 ± 8.3	35.3 ± 8.4	13.8 ± 3.5	10.8 ± 0.8
SGP-392029a	⋯	${3.81}_{-1.27}^{+1.68}$	<4.74	<7.44	10.5 ± 5.8	7.7 ± 5.8	5.2 ± 3.8	5.2 ± 1.3	3.11 ± 0.24
SGP-392029b	⋯	${3.28}_{-0.92}^{+1.03}$	<3.27	<4.56	9.8 ± 5.4	18.7 ± 6.7	28.9 ± 7.0	8.6 ± 2.2	5.15 ± 0.28
SGP-499646	${4.68}_{-0.34}^{+0.49}$	${4.37}_{-1.10}^{+1.32}$	<4.55	4.76 ± 3.37	5.8 ± 5.9	10.8 ± 8.1	41.4 ± 8.6	18.7 ± 3.0	14.1 ± 1.6
UR72S	3.35	2.84 ${}_{-0.88}^{+1.33}$	10.78 ± 4.84	12.20 ± 5.46	35.4 ± 7.3	37.0 ± 8.2	56.0 ± 8.6	⋯	8.3 ± 0.9
UR72Sa	⋯	${2.42}_{-1.00}^{+1.69}$	10.78 ± 4.84	12.20 ± 5.46	5.9 ± 4.3	7.8 ± 5.6	31.1 ± 14.4	⋯	3.57 ± 0.23
UR72Sb	⋯	${2.39}_{-0.77}^{+1.06}$	<14.53	<16.37	24.1 ± 7.8	28.5 ± 8.9	17.5 ± 12.0	⋯	0.96 ± 0.21
UR72Sc	⋯	${2.22}_{-0.85}^{+1.42}$	<14.53	<16.37	11.0 ± 6.3	5.1 ± 3.4	7.3 ± 5.3	⋯	3.72 ± 0.85
UR73S	3.45	${2.53}_{-0.82}^{+1.25}$	48.96 ± 10.13	56.22 ± 11.44	38.6 ± 7.7	39.6 ± 8.7	63.7 ± 9.1	⋯	7.6 ± 0.5
UR73Sa	⋯	${3.60}_{-0.98}^{+1.35}$	48.96 ± 10.13	56.22 ± 11.44	42.3 ± 5.2	37.8 ± 6.8	62.4 ± 7.5	⋯	2.95 ± 0.29
UR73Sb	⋯	${3.12}_{-1.11}^{+1.54}$	<9.88	<10.02	4.7 ± 3.1	19.3 ± 6.7	11.7 ± 8.2	⋯	4.65 ± 0.45
HELMS_RED_10	${4.62}_{-0.63}^{+0.75}$	${3.17}_{-0.63}^{+1.36}$	43.91 ± 9.55	48.19 ± 10.54	33.6 ± 5.7	53.9 ± 6.5	86.5 ± 6.9	37.9 ± 4.4	42.7 ± 0.9
HELMS_RED_10a	⋯	${2.37}_{-0.69}^{+0.83}$	<18.64	<21.62	16.8 ± 9.0	17.3 ± 12.4	58.0 ± 17.7	3.1 ± 0.2	2.00 ± 0.35
HELMS_RED_10b	⋯	${4.62}_{-1.14}^{+1.31}$	43.91 ± 9.55	48.19 ± 10.54	14.5 ± 8.7	36.1 ± 13.3	22.5 ± 16.2	34.8 ± 2.6	22.5 ± 1.1
HELMS_RED_23	4.20 ${}_{-0.61}^{+0.63}$	${3.15}_{-0.81}^{+1.19}$	30.65 ± 8.94	35.90 ± 9.18	48.2 ± 6.7	87.6 ± 6.3	97.2 ± 7.4	42.1 ± 4.9	47.7 ± 0.9
HELMS_RED_23a	⋯	${4.67}_{-1.16}^{+1.39}$	<9.32	<15.90	9.3 ± 6.6	6.8 ± 5.1	20.5 ± 15.0	11.5 ± 1.3	9.9 ± 0.2
HELMS_RED_23b	⋯	${4.92}_{-1.13}^{+1.40}$	9.13 ± 5.35	17.54 ± 6.99	12.8 ± 7.7	5.8 ± 4.1	13.2 ± 8.9	17.1 ± 2.0	14.7 ± 0.8
HELMS_RED_23c	⋯	${2.59}_{-0.73}^{+0.88}$	21.52 ± 7.16	35.9 ± 9.18	25.9 ± 11.3	61.6 ± 11.0	47.4 ± 22.3	13.4 ± 1.6	11.5 ± 0.7
HELMS_RED_68	${3.60}_{-0.64}^{+0.63}$	${3.11}_{-0.79}^{+1.08}$	10.33 ± 5.02	15.13 ± 6.02	55.4 ± 5.6	73.9 ± 6.1	76.1 ± 6.5	32.7 ± 3.8	24.1 ± 3.1

Note. ${z}_{\mathrm{FIR}}$ is the photometric redshift derived from FIR SED fitting by Ivison et al. (2016) and Duivenvoorden et al. (2018). We list the spectroscopic redshift (three decimals) instead of ${z}_{\mathrm{FIR}}$ whenever available. ${z}_{\mathrm{MAGPHYS}}$ is the photometric redshift derived from magphys+photo-z SED fitting in Section 4. For nondetections in IRAC, we quote the 3σ upper limits. The ALMA flux densities are from Oteo et al. (2017) and this work. The deblended SPIRE and SCUBA-2/LABOCA flux densities for the individual components are described in Section 3.

A machine-readable version of the table is available.

Download table as: DataTypeset images: 1 2

For lensed sources, however, a significant fraction of the emission seen at Spitzer/IRAC is likely due to the foreground lensing galaxies, and higher-resolution optical/NIR imaging is required to perform source/lens deblending. Therefore, we can only place an upper limit on the IRAC photometry for the background lensed ultrared sources until we are able to deblend the source/lens photometry. For a few lensed ultrared sources, high-resolution Hubble Space Telescope (HST) imaging is being acquired, and the deblending results will be presented in A. Brown et al. (2019, in preparation).

Within the high-resolution sample of the 63 original Herschel-selected ultrared sources, 17 ultrared sources are resolved into multiple components. A total of 86 individual submillimeter sources were identified from the 63 ultrared sources as a result of the multiple components (more discussion on the multiplicity in Section 5.2). Corresponding Spitzer/IRAC counterparts are identified based on the high-resolution positions from ALMA, NOEMA, and SMA. We use the probabilistic deblender XID+ (Hurley et al. 2017), which is a new prior-based source extraction tool in confusion-dominated maps, with the positional priors from the high-resolution interferometry data to disentangle the SPIRE flux densities over the subcomponents. XID+ is developed and has been tested on SPIRE maps using a probabilistic Bayesian framework, which includes prior information and uses the Bayesian inference to obtain the full posterior probability distribution function on flux estimates. The deblended flux densities are shown in Table 2. The SCUBA-2 or LABOCA flux densities are split among subcomponents assuming the same relative ratios between the individual components derived from the ALMA 870 μm flux densities (Oteo et al. 2017). We caution though that there is inconsistency between the ALMA and SCUBA-2 photometry that could be partly due to the small field of view of ALMA Band 7 (Oteo et al. 2017). We use the SCUBA-2/LABOCA flux densities in the following SED fitting in Section 4.

For the ultrared sources that currently lack high-resolution interferometry data, we use the SCUBA-2 or Herschel positions and cross-match with any IRAC sources that are within 2σ positional uncertainties of SCUBA-2 or Herschel. For sources with robust SCUBA-2 detections (i.e., signal-to-noise ratio, S/N > 3), we combine a statistical positional accuracy of σ_pos = 0.6 × FWHM/(S/N) (Ivison et al. 2007) and the JCMT pointing accuracy of 2''–3''. For sources with low S/N or no SCUBA-2 data, we use the Herschel positions for cross-matching (positional uncertainty σ_H ∼ 6''; Asboth et al. 2016). If the SCUBA-2/Herschel positions are on top of or in close proximity to an IRAC source that is identified as an SDSS or VIKING low-redshift object, the ultrared source is likely lensed by the foreground galaxy and most often the background emission in IRAC is blended with that from the foreground lens. We provide a catalog of the IRAC counterpart candidates identified as the closest IRAC source in search radius in Table 1 but defer counterpart confirmation and robust assessment of lensing until we obtain more high-resolution interferometry data.

4. SED Fitting

In this work, we aim to model the multiwavelength observed SEDs of the cross-identified ultrared sources and derive their physical properties in order to place this population in the context of galaxy formation and evolution. Here we only show the results of the 41 unlensed Herschel ultrared sources (i.e., 63 individual DSFGs at high resolution), which have been cross-identified at multiple wavelengths. To facilitate comparison (i.e., avoid systematic uncertainties in SED fitting due to different choices of SED codes, SED models, and assumptions) with the more abundant DSFGs at z ∼ 2, particularly the well-studied ALMA-LESS (ALESS; ALMA follow-up of the LABOCA submillimeter survey in the Extended Chandra Deep Field South) SMGs at z(median) ∼ 2.5 (Hodge et al. 2013), we use the same SED code, i.e., the updated version of magphys as in da Cunha et al. (2015) and Danielson et al. (2017). magphys relies on a self-consistent energy balance argument to combine stellar emission with dust attenuation and dust emission in galaxies (da Cunha et al. 2008). The updated version extends the model parameter space to account for properties that are more likely applicable to heavily obscured galaxies at high redshifts. The major SED model components include the stellar population synthesis models of Bruzual & Charlot (2003), the delayed-τ star formation histories, the two-component dust attenuation model (Charlot & Fall 2000), dust emission, and radio emission based on the radio-to-FIR correlation.

We use the spectroscopic redshift as the input redshift whenever available. magphys has also been tested and used as a photometric redshift code by leaving the redshift as a free parameter and being simultaneously constrained with all other physical parameters (da Cunha et al. 2015; Battisti et al. 2019). The key feature of the photo-z version (referred to as "magphys+photo-z" hereafter; Battisti et al. 2019) is the self-consistent incorporation of photometric uncertainty into the uncertainty of all the derived properties. By accounting for this effect, magphys+photo-z therefore provides more realistic uncertainties for the physical properties.

We run magphys+photo-z with all of the multiwavelength data and compare with the photo-z's derived from the FIR SEDs (Ivison et al. 2016; Duivenvoorden et al. 2018). The resultant best-fit parameters and associated uncertainties are determined from the posterior probability distributions. We take the median value as the best-fit parameter and the 16th or 84th percentile range as the 1σ uncertainty. In order to analyze the overall properties of the unlensed sample and take into account the associated uncertainties, we stack individual probability distributions of the key physical parameters in Figure 2, including the photometric redshift, stellar mass, SFR, sSFR, mass-weighted stellar population age, V-band dust extinction, total dust luminosity, luminosity-weighted dust temperature, and dust mass. The average (median) properties and associated 16th–84th percentile ranges are listed in Table 5 for each parameter. We display and compare the stacked probability distributions for the following three subsamples:

1.
All unlensed DSFGs: this is the whole unlensed sample containing 63 individual DSFGs (Table 4).
2.
Unlensed ultrared DSFGs: this subsample includes all the unlensed DSFGs that satisfy the ultrared selection criteria¹⁸ (48 DSFGs) including multicomponent sources. At high resolution, some subcomponents no longer meet the ultrared selection criteria so we remove those from this subsample.
3.
Unlensed single-component ultrared DSFGs: this subsample (31 DSFGs) contains the single-component unlensed sources that have FIR photo-z's or the components with spectroscopic redshifts (Table 4). This is the subset on which we can run magphys using the fixed-z version (i.e., FIR photo-z or spectroscopic redshift as input) and compare with the results from magphys+photo-z. As will be shown later, this subsample contains the intrinsically most FIR-luminous and massive DSFGs in our sample, which are likely the best candidates for real z > 4 DSFGs.

**Figure 2.** Normalized stacked posterior probability distributions of the key physical parameters of the whole unlensed sample (63 DSFGs; black), the unlensed ultrared subsample (48 DSFGs; red), and the single-component ultrared subsample (31 DSFGs; blue). The median values of the distributions are indicated by the vertical lines with the corresponding colors (three lines are overlapping if only one line is seen). We also plot the prior distributions of z, age_M, A_V, and T_dust in gray for comparison.
Download figure:
Standard image High-resolution image

Table 5. Average Properties of the Herschel Unlensed Sample Derived from magphys+photo-z SED Fitting

Parameter	All Unlensed	Unlensed Ultrared	Single-comp. Ultrared	ALESS Sample
	(63 DSFGs)	(48 DSFGs)	(31 DSFGs)	(99 DSFGs)
z_phot	${3.3}_{-1.0}^{+1.8}$	${3.3}_{-1.0}^{+1.6}$	${3.5}_{-1.0}^{+1.6}$	${2.7}_{-1.1}^{+1.4}$
log(M_*/M_⊙)	11.45 ${}_{-0.5}^{+0.4}$	${11.55}_{-0.5}^{+0.4}$	11.65 ${}_{-0.4}^{+0.3}$	${10.95}_{-0.8}^{+0.6}$
log(SFR/M_⊙ yr⁻¹)	${2.75}_{-0.5}^{+0.5}$	${2.85}_{-0.6}^{+0.4}$	2.95 ${}_{-0.4}^{+0.4}$	${2.45}_{-0.5}^{+0.4}$
log(sSFR/Gyr⁻¹)	${0.35}_{-0.5}^{+0.5}$	${0.35}_{-0.5}^{+0.5}$	0.35 ${}_{-0.4}^{+0.5}$	0.45 ${}_{-0.6}^{+0.6}$
log(age_M/yr⁻¹)	${8.45}_{-0.3}^{+0.5}$	8.45 ${}_{-0.3}^{+0.5}$	${8.45}_{-0.4}^{+0.5}$	8.35 ${}_{-0.6}^{+0.5}$
A_V	3.6 ${}_{-1.8}^{+2.8}$	${3.9}_{-1.8}^{+2.8}$	${4.1}_{-1.8}^{+2.5}$	1.9 ${}_{-1.0}^{+1.2}$
log( ${M}_{* }/{L}_{H}$ / ${M}_{\odot }/{L}_{\odot }$ )	0.28 ${}_{-0.5}^{+0.6}$	${0.33}_{-0.5}^{+0.6}$	0.33 ${}_{-0.4}^{+0.6}$	−0.13 ${}_{-0.4}^{+0.4}$
log( ${L}_{\mathrm{dust}}$ /L_⊙)	${12.85}_{-0.4}^{+0.4}$	12.95 ${}_{-0.5}^{+0.4}$	13.05 ${}_{-0.3}^{+0.3}$	${12.55}_{-0.5}^{+0.3}$
${T}_{\mathrm{dust}}$ /K	37 ${}_{-6}^{+12}$	${37}_{-6}^{+10}$	${37}_{-6}^{+10}$	43 ${}_{-10}^{+10}$
log(M_dust/M_⊙)	9.35 ${}_{-0.5}^{+0.4}$	${9.45}_{-0.5}^{+0.3}$	9.55 ${}_{-0.4}^{+0.3}$	${8.75}_{-0.4}^{+0.3}$

Note. A machine-readable table of the physical parameters of individual sources is available. ${z}_{\mathrm{phot}}$ : photometric redshift; M_*: stellar mass; SFR: current SFR defined as the average of the star formation history over the last 10 Myr; sSFR: specific SFR; age_M: mass-weighted stellar population age; A_V: V-band extinction; M_*/L_H: stellar mass to H-band luminosity ratio; L_dust: total dust luminosity = IR luminosity; T_dust: luminosity-weighted dust temperature; M_dust: dust mass. The median values and 16th–84th percentile ranges are determined from the stacked posterior probability distributions shown in Figure 2.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as: Data Typeset image

We also compare the physical properties of the unlensed ultrared sample with those of the ALESS DSFGs below.

4.1. Photometric Redshifts

magphys has been tested as a photometric redshift code on ALESS SMGs (as a preliminary version of magphys+photo-z in da Cunha et al. 2015) by comparing to the spectroscopic redshifts from Danielson et al. (2017). da Cunha et al. (2015) find a good agreement between the magphys-based photometric redshifts and the ALESS spectroscopic redshifts, with a small median relative difference of Δz/(1+z_spec) = −0.005. Battisti et al. (2019) further demonstrated the success of magphys+photo-z in estimating redshifts and physical properties based on over 4000 IR-selected galaxies at 0.4 < z < 6.0 in the COSMOS field with robust spectroscopic redshifts. They achieved high photo-z accuracy (median offset Δz/(1+z_spec) ≲ 0.02) and low catastrophic failure rates (η ≲ 4%; a catastrophic failure is defined as a source with Δz/(1+z_spec) > 0.15) over all redshifts. The median value and uncertainties on the photometric redshift are determined in a self-consistent manner as all the other physical parameters, as the likelihood distributions of the redshift and physical parameters are computed simultaneously. Battisti et al. (2019) also demonstrated that the choice of priors, especially the nonuniform prior for the model redshift distribution (Figure 2), does not introduce significant redshift bias in the results.

We compare the photometric redshifts from magphys+photo-z and FIR SED fitting with spectroscopic redshifts based on the five sources with available spectroscopic redshifts (Figure 3, top panel). The mean relative offset Δz/(1+z_spec) is 0.096 for the FIR method and −0.005 for magphys+photo-z, while the median relative offset is 0.076 for the FIR method and −0.047 for magphys+photo-z. The magphys-derived redshifts on average (both mean and median) are more accurate than the FIR method. We further make the comparison between the magphys-derived redshifts and the FIR-derived redshifts for the entire unlensed sample (Figure 3, bottom panel). The magphys-derived redshifts, with a median value of 3.3, are systematically lower than the FIR-derived redshifts with a median value of 3.7. Ivison et al. (2016) and Duivenvoorden et al. (2018) compared the FIR-derived redshifts with available spectroscopic redshifts of a larger sample and found a median relative offset of Δz/(1+z_spec) = 0.08. We compare the magphys redshifts with the expected true redshifts based on their tests.¹⁹ Given the additional data and constraining power in the NIR along with the FIR photometry, the magphys-derived redshifts are on average more consistent with the expected true redshifts based on this comparison (Figure 3 bottom panel), although the error bars are larger. The large error bars from magphys+photo-z are partly due to the fact that we do not have enough data to constrain the full SED but reflect more realistic uncertainties than the errors for the FIR-only photo-z's, which are based on the templates used. Obtaining more rest-frame UV/optical data would further improve the accuracy and precision (Battisti et al. 2019).

**Figure 3.** Top: we have five sources for directly comparing photometric redshifts (magphys and FIR) with available spectroscopic redshifts. The red solid line represents the 1:1 ratio. Bottom: comparison of the photometric redshifts derived from the FIR photometry only and from magphys+photo-z multiwavelength SED fitting for the unlensed sample. The blue line shows the median relative offset (Δz/(1+z) = 0.08) of the FIR-derived photometric redshifts from the available spectroscopic redshifts based on reliability tests in Ivison et al. (2016) and Duivenvoorden et al. (2018). The magphys-derived redshifts are systematically lower than the FIR-derived redshifts but are on average more consistent with the expected spectroscopic redshifts.
Download figure:
Standard image High-resolution image

4.2. Stellar Mass versus SFR

Galaxy surveys at low and high redshifts have shown that star-forming galaxies form a power-law relation between their SFR and stellar mass, known as the main sequence (MS) of star-forming galaxies (e.g., Daddi et al. 2007; Magdis et al. 2010; Speagle et al. 2014). Figure 4 shows the stellar masses versus SFRs of the unlensed ultrared DSFGs compared to those of ALESS DSFGs. The green solid line shows the star-forming MS at the median redshift of our sample, z = 3.3, from Speagle et al. (2014), and the dashed lines show three times above or below the MS. The wide spread of the stellar masses and SFRs of the ALESS DSFGs compared to the star-forming MS suggests that these DSFGs are not a homogeneous population, with some lying significantly (>3 times) above the MS thus defined as starbursts and some being consistent with the MS (<3 times) but just at the high-mass end of this relation (da Cunha et al. 2015). This bimodal distribution is also in line with theoretical predictions, e.g., simulations by Hayward et al. (2012). However, the fraction of starbursts at high redshift depends on how we define the normal star-forming MS at these redshifts. Because the sSFR of the MS predicted by Speagle et al. (2014) continues increasing with redshift, this fraction therefore would be lower than if the MS flattens at z ∼ 2 as suggested by some other studies (e.g., Weinmann et al. 2011; González et al. 2014). The star-forming MS is also dependent on the stellar mass, and the slope of the MS is found to be steeper at low masses (log(M_*/M_⊙) < 10.5) and flattens at the high-mass end out to z of 2.5 (e.g., Whitaker et al. 2014; Leja et al. 2015). A shallower slope at high masses than the one in Speagle et al. (2014) could mean a higher starburst fraction.

**Figure 4.** SFR vs. stellar mass. The red circles are our unlensed *Herschel* ultrared DSFGs while the gray circles denote the ALESS DSFGs at z ∼ 2.5 from da Cunha et al. (2015). The typical error bar derived from magphys SED fitting is shown at the lower-right corner. The green solid line represents the star-forming main sequence from Speagle et al. (2014) at the median redshift z = 3.3 of our sample, and the dashed lines are three times above or below this relation. The blue open circles mark the single-component ultrared DSFGs with parameters derived from magphys+photo-z, while the orange triangles are derived from fixing the redshift to the FIR photo-z or known spectroscopic redshift as input.
Download figure:
Standard image High-resolution image

The Herschel ultrared DSFGs on average have higher stellar masses and SFRs than the ALESS DSFGs, with a median stellar mass of (3.7 ± 0.2) × 10¹¹ M_⊙ and a median SFR of 730 ± 30 M_⊙ yr⁻¹ (modulo assumptions about the IMF; e.g., Romano et al. 2017; Zhang et al. 2018). The 16th–84th percentile ranges of the stacked stellar mass and SFR probability distributions are (1.1–8.9) × 10¹¹ M_⊙ and 180–1800 M_⊙ yr⁻¹. Almost all of the ultrared DSFGs have sSFRs higher than 1 Gyr⁻¹. The open blue circles mark the single-component ultrared sources, which are mostly hyperluminous IR galaxies (HyLIRGs; L_IR ≥ 10¹³ L_⊙) and are at the high-mass end of our sample. Because leaving redshift as a free parameter introduces an extra degree of freedom and degeneracy in SED fitting than fixing the redshift, here we also compare with the magphys results by fixing the input redshift to the FIR-derived photo-z or spectroscopic redshift for the single-component ultrared sample (open triangles). The stellar masses and SFRs from the photo-z and fixed-z versions are generally consistent with each other although the SFRs from the fixed-z version are higher, which is mostly driven by the higher FIR-derived photo-z than the magphys-derived photo-z as shown in Figure 3. Nevertheless, they are all consistent with and at the high-mass end of the star-forming MS by Speagle et al. (2014). If the MS slope flattening continues at z > 3 at the most massive end where our ultrared DSFGs reside, a shallower MS slope and lower sSFR than those in Speagle et al. (2014) would infer a higher starburst fraction, but major mergers, which trigger short-phased enhanced SFRs and thus lie significantly above the MS, are not a dominant driver of our ultrared DSFG population.

4.3. Dust Properties

Figure 5 shows the comparisons of the total dust luminosity, dust temperature, V-band extinction, and dust mass between the unlensed Herschel ultrared DSFGs and ALESS DSFGs. The Herschel ultrared DSFGs have a median total dust luminosity of (9.0 ± 2.0) × 10¹² L_⊙, a dust mass of (2.8 ± 0.6) × 10⁹ M_⊙, a luminosity-averaged dust temperature of 38 ± 2 K, and a V-band extinction of 4.0 ± 0.3. Again, we show the properties of the single-component ultrared subsample derived from both the photo-z version and the fixed-z version. As noted in da Cunha et al. (2015) and other SMG studies (e.g., Chapman et al. 2005; Wardlow et al. 2010; Magnelli et al. 2012; Swinbank et al. 2014), there exists a correlation between the total dust luminosity and the average temperature of DSFGs. The dust temperatures of the ultrared sample distribute in a narrower range across the dust luminosity although the error bars are large. This narrow range is almost entirely driven by the prior distribution as shown in Figure 2, i.e., T_dust is basically not constrained by the data. Although we have FIR–submillimeter/millimeter data to locate the dust emission peak, the degeneracy between T_dust and redshift still exists. More scatter in T_dust can be seen for the single-component subsample if we fix the input redshift to the FIR photo-z or spectroscopic redshift, i.e., break the T_dust–redshift degeneracy. The single-component subsample (fixed-z version) shifts to the higher-luminosity end, with most of them being HyLIRGs, and the luminosity–temperature correlation is very weak. The V-band extinction values are widely spread over the range from ∼1–7.5 for the whole ultrared sample and the single-component subsample and is on average higher than those of ALESS DSFGs by ∼2 magnitudes. High A_V values (A_V > 4) have also been found in other SMGs (e.g., Hopwood et al. 2011; Ma et al. 2015). The Herschel ultrared DSFGs, which have significantly higher IR luminosities, shift to the high dust-mass end. The sample selection methods factor into the differences we observe here for the two DSFG samples. The ALESS SMGs are selected based on a single flux density limit (S_870μm > 4.2 mJy; Swinbank et al. 2014), while the Herschel ultrared DSFGs are selected by their rising SPIRE flux densities and limited by the sensitivity that Herschel probes, i.e., the Herschel ultrared selection naturally chooses higher 870 μm flux density sources than the ALESS sample. We also notice that the Herschel ultrared sample distributes similarly on the A_v versus log L_dust plane as the A_v–logM_dust plot without a clear correlation. Any correlation may be diluted by the statistical errors on these parameters due to the limited sampling of the SEDs. Only by obtaining spectroscopic redshifts and a better sampling of the SEDs can we ultimately break the degeneracy between these parameters (e.g., dust temperature–redshift degeneracy) and investigate potential intrinsic correlations between the physical properties.

4.4. The Average SED of the Unlensed Ultrared Sample

Figure 6 shows the best-fit SED shifted to the rest-frame wavelength for each unlensed ultrared DSFG (gray) to highlight the intrinsic SED variations of these sources. We generate the median SED (red) and mean SED (orange) of this sample by averaging the flux densities of all the best-fit SEDs at each rest-frame wavelength in the same manner as for the ALESS DSFGs. The average (mean) SED of the ALESS DSFGs from da Cunha et al. (2015) is also overlaid for comparison. The comparison of the average SEDs shows that the Herschel ultrared DSFGs on average are more luminous at FIR–submillimeter, more dust obscured, and peak at longer wavelengths than ALESS DSFGs at z ∼ 2.5. We caution though that our current sample of the SED fitting analysis is smaller and our SEDs are not as well sampled as the ALESS DSFGs due to limited photometry (e.g., the rest-frame UV region is not constrained); therefore, we do not attempt to further compare the detailed shape of the SED.

5. Discussion

5.1. Redshift Distribution

The raw photometric redshift distribution of the 300 Herschel ultrared DSFGs are shown in Figure 7 (left), which has a median redshift of 3.7 ± 0.7 based on the FIR method as we do not have magphys-derived redshifts for all of them. We compare it with the photometric redshift distribution of ALESS DSFGs from Simpson et al. (2014). Danielson et al. (2017) present spectroscopic redshifts for 52 ALESS DSFGs, and the distribution is consistent with the photometric redshift distribution for these sources. Here we use the photometric redshift distribution for comparison because it is more complete although less precise. The median redshift of ALESS DSFGs, 2.5 ± 0.2, is consistent with that of Chapman et al. (2005), which is relied on radio-wavelength counterpart identification, although Simpson et al. (2014) shows a higher fraction of high-redshift sources than the earlier work. We also overlay the spectroscopic redshift distribution of the SPT DSFGs from Strandet et al. (2016, 2017; also Weiß et al. 2013), which are almost purely gravitationally lensed sources. The observed median redshift of the SPT sample, 3.9 ± 0.4, is higher than other samples due to two major selection effects: longer selection wavelengths and gravitational lensing (Béthermin et al. 2015). After correcting for the lensing effect, the median redshift of the SPT DSFGs decreases to 3.1 ± 0.3. Figure 7 (Right) demonstrates the normalized dN/dz that scales to the total number of source in each sample. We do not attempt to correct the raw distribution for any selection effects, due to the complicated selection functions for our sample. Spectroscopic observations are required to ascertain the redshift distribution of our ultrared sample especially the ones at z > 4.

5.2. Multiplicity

Source blending or multiplicity has always been a concern for 500 μm risers due to the relatively large SPIRE beam at 500 μm. Different multiplicity rates have been reported in the literature, depending upon various selection criteria and instruments used (e.g., Cowley et al. 2015; Simpson et al. 2015; Scudder et al. 2016). Now that we have obtained high-resolution data for a significant subset of ultrared sources, we are able to investigate multiplicity rates and fractional contribution from the brightest components. Within our high-resolution sample, ∼27% (17 out of 63) of the Herschel ultrared sources are resolved into multiple components of two or more. The multiplicity rate is about 39% (16 out of 41) for the unlensed sample. The brightest components seen in ALMA contribute 41%–80% of the total ALMA flux. This fraction is lower, in the range of 15%–59%, if we use the SCUBA-2/LABOCA flux at similar wavelengths to the total flux as the ALMA observations may not recover the total flux. Based on SMA follow-up observations at 1.1 mm and/or 870 μm of thirty-six 500 μm risers from various Herschel fields, J. Greenslade et al. (2019, in preparation) found a multiplicity rate of 33% with a fractional contribution of 50%–75% due to the brightest components. The deblending results from XID+ suggest that the brightest components contribute to 35%–87% of the total 500 μm flux. Donevski et al. (2018) investigated multiplicity using simulated SPIRE maps and compared the extracted flux of the brightest galaxy to the total flux, resulting in an average brightest-galaxy fraction of 64%. The average observed brightest-galaxy fraction of our sample is consistent with this prediction.

Multiple components at the same redshift will not have a serious impact on the determination of redshift, and physical properties such as IR luminosity and SFR are determined for the combined system. Several Herschel ultrared sources have been spectroscopically confirmed as major mergers (e.g., SGP-196076 at z = 4.425 (Oteo et al. 2016a), ADFS-27 at z = 5.655 Riechers et al. 2017), or galaxies in the core of protoclusters (e.g., SGP-354388 at z = 4.002; Oteo et al. 2018). Detailed SED analysis, especially their stellar properties, with additional high-resolution optical/NIR data, for the individual merging systems or protoclusters will be published in forthcoming papers. Multiple components at different redshifts, e.g., blending with foreground objects, will produce composite SEDs that may lead to an intermediate redshift estimate. Once decomposed, the individual components might not satisfy the 500 μm riser selection criteria. Duivenvoorden et al. (2018) derived from mock observations that 60% of the detected HeLMS galaxies pass the selection threshold due to flux boosting partially caused by blending with foreground objects. Our SPIRE deblending results with XID+, which are based on high-resolution positional priors of mostly H-ATLAS galaxies, suggest that ∼20% of our sources would not pass the selection criteria of 500 μm risers if there is no blending.

5.3. Unlensed Fraction

Based on the high-resolution subsample, 65% of the ultrared sources (41 out of 63) do not have clear signatures of lensing and thus are classified as unlensed. This fraction is higher, 73% (63 out of 86), if we count the individual components. We do not exclude the possibility that a small fraction of them might be moderately lensed (a lensing magnification factor of 1 < μ < 2). This would not significantly affect our derived properties.

Donevski et al. (2018) compared the observed lensed fractions of 500 μm risers in different fields, H-ATLAS (Ivison et al. 2016), HeLMS (Asboth et al. 2016), and HerMES (Dowell et al. 2014), with predicted fractions using the Béthermin et al. (2017) model and applying the same selection criteria in each study. The observed lensed fractions are consistent with the corresponding predicted lensed fractions. Because the HeLMS ultrared sources are selected with a higher flux cut at 500 μm, the lensed fraction is the highest (75%), compared to the lowest lensed fraction (28%), in H-ATLAS. Our sample contains sources from both H-ATLAS and HeLMS, therefore the observed lensed fraction (35%) is in between as expected. This number is more toward the lower end because most of the sources in our high-resolution sample are from H-ATLAS.

We acknowledge that there could be lenses that we are missing with the current ancillary data, for example, those with small Einstein radii and faint, high-redshift lensing galaxies. High-resolution imaging (e.g., with HST/James Webb Space Telescope (JWST)/ALMA) and spectroscopic confirmation are needed to truly address the lensing fraction.

5.4. SFR Surface Density

Figure 8 shows the SFR as a function of dust continuum size for our Herschel ultrared DSFGs as well as other high-redshift (z > 3) galaxies and quasar hosts in the literature. All sources in the literature comparison sample have dust continuum size measurements. The SFR surface density is defined as Σ_SFR = SFR/2R_aR_b, where R_a = FWHM_major/2 and R_b = FWHM_minor/2 are the measured semimajor and semiminor axes. The sizes and areas are measured by carrying out 2D elliptical Gaussian fitting on the high-resolution ALMA dust continuum images at 870 μm (Oteo et al. 2017). The dashed lines denote the constant SFR surface densities at 10, 100, and 1000 M_⊙ yr⁻¹ kpc⁻². Most of our sources are above the Σ_SFR = 100 M_⊙ yr⁻¹ kpc⁻² curve, and some are at or close to the Eddington limit for radiation-pressure-supported starbursts (Thompson et al. 2005). The radiation pressure would drive dusty gaseous outflows, which have been observed in starbursting galaxies (e.g., Martin 2005; Diamond-Stanic et al. 2012; Spilker et al. 2018). Future observations are required to confirm this. The very high Σ_SFR could be explained by either high star formation efficiency, high gas fraction, or both. Gas-rich mergers are a viable mechanism for triggering the compact, enhanced star formation (Mihos & Hernquist 1996). For example, SGP-196076, one of our unlensed ultrared sources, has been spectroscopically confirmed to be at least two interacting galaxies at z = 4.425 that will eventually merge (Oteo et al. 2016a). However, major mergers do not seem to be the dominant driver of our ultrared DSFGs as discussed in Section 4.2.

5.5. Space Density, SFR Density, and Stellar Mass Density

The space density of the H-ATLAS ultrared DSFGs in the redshift range 4 < z < 6 is estimated to be ∼6 × 10⁻⁷ Mpc⁻³ after completeness and duty-cycle corrections (Ivison et al. 2016), while this number is about an order of magnitude smaller (7 × 10⁻⁸ Mpc⁻³) for the HeLMS sample (Duivenvoorden et al. 2018), due to the higher flux cut (i.e., HeLMS S₅₀₀ > 63 mJy and H-ATLAS S₅₀₀ > 30 mJy).

DSFGs at z > 4 have been proposed to be the star-forming progenitors of the population of massive, quiescent galaxies at z ∼ 3 uncovered from NIR surveys (e.g., Nayyeri et al. 2014; Simpson et al. 2014; Straatman et al. 2014; Toft et al. 2014). However, the quiescent galaxies represented by the mass-limited (log(M_*/M_⊙) > 10.6) sample at 3.4 < z < 4.2 in the ZFOURGE survey, whose star formation is predicted to occur at z ∼ 5, have a comoving space density of ∼2 × 10⁻⁵ Mpc⁻³ (Straatman et al. 2014). This is more than a factor of 30 higher than the H-ATLAS ultrared sample (Ivison et al. 2016) and more than two orders of magnitude higher than the HeLMS sample (Duivenvoorden et al. 2018). Based on the space density comparison, Ivison et al. (2016) and Duivenvoorden et al. (2018) conclude that the Herschel ultrared sample, which is limited by the flux density levels probed by Herschel, cannot account for the formation of massive, quiescent galaxies at z ∼ 3. Our Herschel-selected ultrared sample contains more FIR-luminous and thus rarer DSFGs than the star-forming progenitors of the massive, quiescent galaxies.

The SFRD of DSFGs can be estimated by summing SFRs of all the sources divided by the comoving volume contained in a redshift range. Although it has become clear that the UV/optical is insufficient to probe the total cosmic star formation at z < 3, there is no consensus on the significance of the contribution of dust-obscured star formation in DSFGs at z > 3, due to the lack of complete surveys. The SFRD estimates at 4 < z < 5 based on different samples vary by three orders of magnitude from ∼10⁻⁴ to 10⁻¹ M_⊙ yr⁻¹ Mpc⁻³ (e.g., Dowell et al. 2014; Rowan-Robinson et al. 2016; Bourne et al. 2017; Donevski et al. 2018; Duivenvoorden et al. 2018), with our HeLMS ultrared sample representing the lower limit of ∼10⁻⁴ M_⊙ yr⁻¹ Mpc⁻³.

The space density and SFRD of DSFGs at z > 4 have been widely explored as summarized above. However, hardly any estimates have been made yet on the stellar mass density contribution of DSFGs at z > 4, mainly due to the highly dust-obscured nature of these objects, and therefore the difficulty in detecting the rest-frame optical stellar emission to constrain their stellar masses. Our Spitzer follow-up sample provides the first attempt to constrain the stellar mass density of ultrared DSFGs at z > 4. Figure 9 shows the comparison of the stellar mass density (SMD) as a function of redshift for different populations. We construct the SMD of our H-ATLAS²⁰ unlensed ultrared DSFGs in seven redshift bins and have corrected for the H-ATLAS ultrared sample completeness (Ivison et al. 2016) and duty cycle (assuming a starburst duty cycle of ∼100 Myr). The total SMD curve is a simultaneous fit to the total SMD of the K_s-selected galaxies at z < 3.5 from the COSMOS/UltraVISTA survey and the SMD of the UV-selected samples at z > 3.5 (Muzzin et al. 2013; see also Davidzon et al. 2017 and references therein). The UV samples at z > 3.5 contain galaxies that are based on dropout selection methods, e.g., Lyman break galaxies (LBGs), B-dropout etc. (Stark et al. 2009; Labbé et al. 2010; González et al. 2011; Lee et al. 2012). The total SMD increases with cosmic time as galaxies build up their stellar masses through star formation. The maximum contribution of our Herschel ultrared sources to the total SMD is 1.7% at the z = 3.25 bin. The SMD of massive, quiescent galaxies also evolves with redshift, as represented by the K_s-selected, mass-limited (log(M_*/M_⊙) > 10.6) ZFOURGE quiescent galaxies (Straatman et al. 2014). The DSFGs at 3 < z < 4 in the ZFOURGE sample and the H–[4.5] color-selected "HIEROs" at z > 3.5 as defined in Wang et al. (2016), which are also massive DSFGs, have comparable SMDs to the quiescent galaxies at 3 < z < 4. However, our Herschel ultrared DSFGs at z ∼ 5 have ∼2 orders of magnitude lower SMDs than the quiescent sample at z ∼ 3 and the dusty star-forming HIEROs at similar redshifts, which cannot be reconciled even by increasing the stellar mass limit to log(M_*/M_⊙) > 11.0. Again, this suggests that our ultrared sample cannot account for the star-forming progenitors of the massive, quiescent galaxies at z ∼ 3 found in NIR surveys (although likely include the progenitors of the most extremely massive, quenched systems), while other selections, such as HIEROs (Wang et al. 2016) and S_850μm- or S_870μm-selected DSFGs at z > 4 (e.g., Oteo et al. 2016b; Michałowski et al. 2017) likely include the majority of the progenitors of the massive, quiescent galaxies at z ∼ 3.

**Figure 9.** Stellar mass density (SMD) as a function of redshift. The black curve shows a simultaneous fit to the total SMD of the K_s-selected galaxies at z < 3.5 from the UltraVISTA survey and the SMD of the UV-selected samples at z > 3.5 (Stark et al. 2009; Labbé et al. 2010; González et al. 2011; Lee et al. 2012; Muzzin et al. 2013). The orange squares are the K_s-selected, mass-limited (log(M_*/M_⊙) > 10.6) quiescent galaxies in the ZFOURGE survey (Straatman et al. 2014). The green diamond shows the K_s-selected DSFGs at 3 < z < 4 from the ZFOURGE sample (Spitler et al. 2014). The blue circles are the H–[4.5] color-selected massive DSFGs (so-called "HIEROs" as defined in Wang et al. 2016) using the ZFOURGE catalog and applying the same stellar mass cut, i.e., log(M_*/M_⊙) > 10.6.
Download figure:
Standard image High-resolution image

Nevertheless, our Herschel ultrared sample contains the intrinsically most FIR-luminous (i.e., HyLIRGs) and massive galaxies in the early universe that are extremely interesting by themselves. These ultrared DSFGs are crucial to understanding the drivers of extreme star formation and the assembly and evolution of massive galaxies with cosmic time.

6. Summary and Conclusions

We have presented a large Spitzer follow-up program of 300 Herschel-selected ultrared DSFGs. For a significant subset, we have obtained high-resolution interferometry data such that we can pinpoint the locations and securely cross-identify Spitzer counterparts, and classify them as lensed or unlensed based on the morphology and the presence or absence of low-redshift foreground galaxies. For the rest of the sample, we have selected Spitzer counterpart candidates based on the SCUBA-2 positions. We have provided a catalog of all the cross-matched sources, including their positions, Spitzer/IRAC magnitudes and flux densities, as well as multiwavelength photometry. In this paper, we have focused on analyzing the unlensed sample by performing magphys SED modeling with the multiwavelength photometry to derive their physical properties and compare with the more abundant z ∼ 2 DSFG population. We have also estimated the stellar mass density as a function of redshift and compared with massive, quiescent galaxies at lower redshifts. Our main results are summarized as follows.

1.
Within the 63 Herschel ultrared sources that have high-resolution data, ∼65% (41 out of 63) appear to be unlensed, and about 27% (17 out of 63) are resolved into multiple components. Some of the deblended components are no longer 500 μm risers. About 20% of the original ultrared sources would not pass the selection criteria without blending with other sources at lower redshifts.
2.
We run magphys+photo-z on the unlensed sample to simultaneously constrain their redshifts and physical properties. The ultrared sample has a median redshift of 3.3, which is lower than the median value of the FIR-derived redshifts, and the 16th–84th percentile range is from 2.3 to 4.9. The magphys-based photometric redshifts are more in line with the expected true redshifts based on the test by comparing with spectroscopic redshifts.
3.
We derive the median properties of the whole unlensed sample, the unlensed ultrared subsample, and the single-component ultrared subsample from stacked probability distributions. The unlensed ultrared sample has a median stellar mass of (3.7 ± 0.2) × 10¹¹ M_⊙, an SFR of 730 ± 30 M_⊙ yr⁻¹, a total dust luminosity of (9.0 ± 2.0) × 10¹² L_⊙, a dust mass of (2.8 ± 0.6) × 10⁹ M_⊙, and a V-band extinction of 4.0 ± 0.3. These properties are all higher than those of the ALESS DSFGs.
4.
We estimate the stellar mass densities of our ultrared DSFGs as a function of redshift. The stellar mass density at z ∼ 5 is significantly lower than that of the massive, quiescent galaxies at lower redshifts from the ZFOURGE survey and HIREOs at similar redshifts. Combined with the comparison of the space density and SFRD, we conclude that our ultrared sample cannot account for the majority of the star-forming progenitors of the massive, quiescent galaxies. Our sample is limited by the flux density levels probed by Herschel and thus contains more FIR-luminous and rarer DSFGs than the progenitors of the massive, quiescent galaxies found in NIR surveys.
5.
We have identified a sample of unlensed, intrinsic HyLIRGs. These HyLIRGs are potentially extremely valuable for our understanding of galaxy evolution, because they present the most luminous, massive, and active galaxies in the early universe. Future investigations of their detailed kinematics are needed to understand the physical drivers of such extreme star formation.

This paper provides a catalog of high-redshift DSFGs for spectroscopic follow-up observations and future JWST observations to probe the mid-IR and rest-frame UV continuum. Our sample contains largely unlensed DSFGs that are especially advantageous because it avoids uncertainties in lens modeling and differential lensing, which will enable us to draw definite conclusions on the connections between stellar, gas, and dust emission components.

We thank the anonymous referee for the constructive comments that have greatly improved the manuscript. We thank Maria Strandet and Tao Wang for sharing data with us. This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This research was supported by NASA grants NNX15AQ06A, NNX16AF38G, and NSF grant AST-131331. J.L.W. acknowledges support from an STFC Ernest Rutherford Fellowship (ST/P004784/2). L.D. and S.J.M. acknowledge funding from the European Research Council Consolidator Grant: CosmicDust. D.R. acknowledges support from the National Science Foundation under grant number AST-1614213.

This paper makes use of the following ALMA data: ADS/JAO.ALMA #2013.1.00001.S and #2016.1.00139.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), and NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China, and the Chinese Academy of Sciences (grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada.

The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen), responsible for the payload module and for system testing at spacecraft level; Thales Alenia Space (Turin), responsible for the service module; and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors.

SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, University of Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK); and Caltech, JPL, NHSC, University of Colorado (USA). This development has been supported by national funding agencies CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA).

Appendix A: Spitzer/IRAC Cutout Images of Herschel-selected Ultrared Sources

Figure 10 and its figure set shows the 60'' × 60'' Spitzer/IRAC cutouts centered on the Herschel positions (green cross) for the 300 ultrared sources in Spitzer program PID 13042. The white circle shows the 36'' FWHM beam size of SPIRE at 500 μm. The 3'' radius cyan circle denotes the positions of the SDSS sources in the field, and the yellow circle shows the positions of the VIKING sources. The red circles indicate the SCUBA-2 positions and corresponding closest IRAC counterparts (red plus). The magenta circles show the high-resolution positions from ALMA, NOEMA, and/or SMA.

Appendix B

Table 6 lists the sources with additional Submilimeter/millimeter photometry from SMA, MUSIC, and ACT.

Table 6. Sources with Submillimeter/millimeter Photometry from SMA, MUSIC, and ACT

Source name	MUSIC	SMA	ACT	MUSIC	ACT	MUSIC	ACT	MUSIC
	0.92 mm	1.1 mm	1.1 mm	1.1 mm	1.4 mm	1.4 mm	2.0 mm	2.0 mm
	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)	(mJy)
HELMS_RED_1		28.6 ± 2.3	72.32 ± 6.26	35.11 ± 2.62	12.49 ± 1.74
HELMS_RED_2		33.9 ± 2.3
HELMS_RED_3	100.3 ± 52.8		35.32 ± 6.24	24.9 ± 12.3	19.50 ± 2.56	16.0 ± 6.9	6.14 ± 1.76
HELMS_RED_4	65.2 ± 57.3	21.3 ± 1.9		32.8 ± 12.4		19.4 ± 6.7		8.4 ± 7.4
HELMS_RED_6				32.3 ± 9.8		8.8 ± 4.3		13.1 ± 5.7
HELMS_RED_7	108.7 ± 35.1			62.4 ± 13.5		18.1 ± 5.1
HELMS_RED_10		13.3 ± 2.8
HELMS_RED_13		11.5 ± 1.8

Source name	NOEMA	NOEMA/ALMA
	1.3 mm	3 mm
	(mJy)	(mJy)

G09-59393	4.0 ± 0.6
G09-62610	5.2 ± 0.8	<0.18
G09-81106	9.7 ± 1.3	0.24 ± 0.04
G09-83808	19.4 ± 2.0	0.66 ± 0.12
NGP-101333	10.8 ± 1.3	<0.25
NGP-111912	4.7 ± 0.9	<0.26
NGP-113609	13.0 ± 2.3	<0.26
NGP-126191	12.3 ± 1.7	0.30 ± 0.11
NGP-136156	3.1 ± 0.8	<0.25
NGP-190387	12.2 ± 1.2	0.84 ± 0.14
NGP-206987	9.2 ± 1.8	<0.32
NGP-246114	8.0 ± 1.5	0.42 ± 0.06
NGP-252305	6.5 ± 0.7	<0.29

Download table as: ASCII Typeset image

Spitzer Catalog of Herschel-selected Ultrared Dusty Star-forming Galaxies

Article metrics

Permissions

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Sample Selection, Spitzer Observations, and Ancillary Data

2.1. Herschel/SPIRE at 250, 350, and 500 μm

2.2. Spitzer Observations and Source Photometry

2.3. SCUBA-2/LABOCA at 850/870 μm

2.4. High-resolution Submillimeter/millimeter Data

2.4.1. Continuum Data at 870 μm, 1.1, 1.3, and 3 mm

2.4.2. Spectroscopic Redshifts from Spectral Scans

3. Cross-identification and Spitzer Catalogs

4. SED Fitting

4.1. Photometric Redshifts

4.2. Stellar Mass versus SFR

4.3. Dust Properties

4.4. The Average SED of the Unlensed Ultrared Sample

5. Discussion

5.1. Redshift Distribution

5.2. Multiplicity

5.3. Unlensed Fraction

5.4. SFR Surface Density

5.5. Space Density, SFR Density, and Stellar Mass Density

6. Summary and Conclusions

Appendix A: Spitzer/IRAC Cutout Images of Herschel-selected Ultrared Sources

Appendix B

Footnotes

Spitzer Catalog of Herschel-selected Ultrared Dusty Star-forming Galaxies

Article metrics

Permissions

Share this article

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Sample Selection, Spitzer Observations, and Ancillary Data

2.1. Herschel/SPIRE at 250, 350, and 500 μm

2.2. Spitzer Observations and Source Photometry

2.3. SCUBA-2/LABOCA at 850/870 μm

2.4. High-resolution Submillimeter/millimeter Data

2.4.1. Continuum Data at 870 μm, 1.1, 1.3, and 3 mm

2.4.2. Spectroscopic Redshifts from Spectral Scans

3. Cross-identification and Spitzer Catalogs

4. SED Fitting

4.1. Photometric Redshifts

4.2. Stellar Mass versus SFR

4.3. Dust Properties

4.4. The Average SED of the Unlensed Ultrared Sample

5. Discussion

5.1. Redshift Distribution

5.2. Multiplicity

5.3. Unlensed Fraction

5.4. SFR Surface Density

5.5. Space Density, SFR Density, and Stellar Mass Density

6. Summary and Conclusions

Appendix A: Spitzer/IRAC Cutout Images of Herschel-selected Ultrared Sources

Appendix B

Footnotes