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

The largest Herschel extragalactic surveys, H-ATLAS and HerMES, have selected a sample of “ultrared” dusty star-forming galaxies (DSFGs) with rising SPIRE flux densities (S500 > S350 > S250; the so-called “500 μm risers”) as an efficient way for identifying DSFGs at higher redshift (z > 4). In this paper, we present a large Spitzer follow-up program of 300 Herschel ultrared DSFGs. We have obtained high-resolution Atacama Large Millimeter/submillimeter Array, Northern Extended Millimeter Array, and SMA data for 63 of them, which allow us to securely identify the Spitzer/IRAC counterparts and classify them as gravitationally lensed or unlensed. Within the 63 ultrared sources with high-resolution data, ∼65% appear to be unlensed and ∼27% are resolved into multiple components. We focus on analyzing the unlensed sample by directly performing multiwavelength spectral energy distribution modeling to derive their physical properties and compare with the more numerous z ∼ 2 DSFG population. The ultrared sample has a median redshift of 3.3, stellar mass of 3.7 × 1011 M⊙, star formation rate (SFR) of 730 M⊙ yr−1, total dust luminosity of 9.0 × 1012 L⊙, dust mass of 2.8 × 109 M⊙, and V-band extinction of 4.0, which are all higher than those of the ALESS DSFGs. Based on the space density, SFR density, and stellar mass density estimates, we conclude that our ultrared sample cannot account for the majority of the star-forming progenitors of the massive, quiescent galaxies found in infrared surveys. Our sample contains the rarer, intrinsically most dusty, luminous, and massive galaxies in the early universe that will help us understand the physical drivers of extreme star formation.


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 3 M e yr −1 , and number densities < 10 −4 Mpc −3 ; 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. 2011Cañ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 2 , have revealed a large number of DSFGs. While most of them are z ∼ 1-2 starburst galaxies (e.g., Casey et al. 2012aCasey et al. , 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 500 >S 350 >S 250 ("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 starforming (>500 M e yr −1 ) 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. 2013Riechers et al. , 2017Fudamoto 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 Spilker et al. 2016;Strandet et al. 2016), including the highestredshift 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 2 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 2 HerMES Large Mode Survey (HeLMS) and 188 of the 200 brightest (S 500 > 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.

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 2 . 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 500 >30 mJy and the color selection with S 500 /S 250 1.5 and S 500 /S 350 0.85.
The HeLMS field covers an effective area of ∼274 deg 2 of the equatorial region spanning 23 h 14 m <R.A.<1 h 16 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 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 500 >S 350 >S 250 and a 5σ cutoff S 500 > 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 500 > 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.

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 17 ) 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 6 pixel −1 and an angular resolution of ∼1 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 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 0 at 850 μm, and the astrometry accuracy σ is 2″-3″. LABOCA has an FWHM resolution of 19 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. 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 5 for NOEMA and 0 6-1 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 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 −1 and the beam FWHM sizes are typically 2 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 highresolution 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).  (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.

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. 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 . 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 . 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-tonoise 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 Note. z 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 FIR whenever available. z 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.
(This table is available in machine-readable form.) 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.

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 crossidentified 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 wellstudied 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 twocomponent 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 selfconsistent 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). 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.
We also compare the physical properties of the unlensed ultrared sample with those of the ALESS DSFGs below. 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.

Photometric Redshifts
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  (This table is available in its entirety in machine-readable form.) 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. 19 Given the additional data and constraining power in the NIR along with the FIR photometry, the MAGPHYSderived 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 FIRonly 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).

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 e ) < 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.
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 11 M e and a median SFR of 730± 30 M e yr −1 (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 11 M e and 180-1800 M e yr −1 . Almost all of the ultrared DSFGs have sSFRs higher than 1 Gyr −1 . The open blue circles mark the single-component ultrared sources, which are mostly hyperluminous IR galaxies (HyLIRGs; L IR 10 13 L e ) 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 . 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. 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. 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 12 L e , a dust mass of (2.8 ± 0.6) × 10 9 M e , 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. 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.

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 highredshift sources than the earlier work. We also overlay the spectroscopic redshift distribution of the SPT DSFGs from Strandet et al. (2016Strandet et al. ( , 2017also 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.

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 highresolution data for a significant subset of ultrared sources, we are able to investigate multiplicity rates and fractional contribution from the brightest components. Within our highresolution 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%  . Best-fit SEDs in the rest frame. The gray curves show the individual best-fit SEDs from MAGPHYS+photo-z of the 48 Herschel unlensed ultrared DSFGs, and the red curve and orange curve are the median and mean SED of this sample. The average (mean) SED of the ALESS DSFGs at z∼2.5 is overlaid in blue. Our Herschel ultrared DSFGs are on average more dust obscured and more luminous in the FIR than the z∼2.5 ALESS DSFGs. We do not attempt to further compare the detailed shape of the SED due to the limited filter bands, hence uncertainties in our best-fit SEDs, e.g., the rest-frame UV region is not constrained for our ultrared DSFGs. 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.

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. Figure 8 shows the SFR as a function of dust continuum size for our Herschel ultrared DSFGs as well as other highredshift (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 a R 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 . The dashed lines denote the constant SFR surface densities at 10, 100, and 1000 M e yr −1 kpc −2 . Most of our sources are above the Σ SFR =100 M e yr −1 kpc −2 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  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.

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 −7 Mpc −3 after completeness and duty-cycle corrections (Ivison et al. 2016), while this number is about an order of magnitude smaller (7×10 −8 Mpc −3 ) for the HeLMS sample (Duivenvoorden et al. 2018), due to the higher flux cut (i.e., HeLMS S 500 > 63 mJy and H-ATLAS S 500 > 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 masslimited (log(M * /M e ) > 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 −5 Mpc −3 ). 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 −4 to 10 −1 M e yr −1 Mpc −3 (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 −4 M e yr −1 Mpc −3 .
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 dustobscured 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 20 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 UVselected 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 e ) > 10.6) ZFOURGE quiescent galaxies . 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 e ) > 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 ) 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.
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.

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 highresolution 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 FIRderived 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 singlecomponent ultrared subsample from stacked probability distributions. The unlensed ultrared sample has a median stellar mass of (3.7 ± 0.2) × 10 11 M e , an SFR of 730±30 M e yr −1 , a total dust luminosity of (9.0 ± 2.0) × 10 12 L e , a dust mass of (2.8 ± 0.6) × 10 9 M e , 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 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, masslimited (log(M * /M e ) > 10.6) quiescent galaxies in the ZFOURGE survey ). The green diamond shows the K s -selected DSFGs at 3<z<4 from the ZFOURGE sample . 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 e ) > 10.6.