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

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'' (8333 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.

Zoom In Zoom Out Reset image size

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 06 pixel⁻¹ and an angular resolution of ∼19.

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 (015) of the mosaic pixel.

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 130 at 850 μm, and the astrometry accuracy σ is 2''–3''. LABOCA has an FWHM resolution of 192, 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.1. Continuum Data at 870 μm, 1.1, 1.3, and 3 mm

2.4.2. Spectroscopic Redshifts from Spectral Scans

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).

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

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 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.