Discovery of a Dark, Massive, ALMA-only Galaxy at z ∼ 5–6 in a Tiny 3 mm Survey

, , , , , , , , , and

Published 2019 October 22 © 2019. The American Astronomical Society. All rights reserved.
, , Citation Christina C. Williams et al 2019 ApJ 884 154 DOI 10.3847/1538-4357/ab44aa

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/884/2/154

Abstract

We report the serendipitous detection of two 3 mm continuum sources found in deep ALMA Band 3 observations to study intermediate-redshift galaxies in the COSMOS field. One is near a foreground galaxy at 1farcs3, but is a previously unknown dust-obscured star-forming galaxy (DSFG) at probable zCO = 3.329, illustrating the risk of misidentifying shorter wavelength counterparts. The optical-to-millimeter spectral energy distribution (SED) favors a gray λ−0.4 attenuation curve and results in significantly larger stellar mass and SFR compared to a Calzetti starburst law, suggesting caution when relating progenitors and descendants based on these quantities. The other source is missing from all previous optical/near-infrared/submillimeter/radio catalogs ("ALMA-only"), and remains undetected even in stacked ultradeep optical (>29.6 AB) and near-infrared (>27.9 AB) images. Using the ALMA position as a prior reveals faint signal-to-noise ratio ∼ 3 measurements in stacked IRAC 3.6+4.5, ultradeep SCUBA2 850 μm, and VLA 3 GHz, indicating the source is real. The SED is robustly reproduced by a massive M* = 1010.8M and Mgas = 1011M, highly obscured AV ∼ 4, star-forming SFR ∼ 300 M yr−1 galaxy at redshift z = 5.5 ± 1.1. The ultrasmall 8 arcmin2 survey area implies a large yet uncertain contribution to the cosmic star formation rate density CSFRD(z = 5) ∼ 0.9 × 10−2 M yr−1 Mpc−3, comparable to all ultraviolet-selected galaxies combined. These results indicate the existence of a prominent population of DSFGs at z > 4, below the typical detection limit of bright galaxies found in single-dish submillimeter surveys, but with larger space densities ∼3 × 10−5 Mpc−3, higher duty cycles of 50%–100%, contributing more to the CSFRD, and potentially dominating the high-mass galaxy stellar mass function.

Export citation and abstract BibTeX RIS

1. Introduction

In past decades, single-dish submillimeter surveys have identified populations of massive, dusty star-forming galaxies at z > 1 (e.g., Casey et al. 2014). While these galaxies are rare even at Cosmic Noon (1 < z < 3) when the star formation activity in the universe peaks, their contribution to the cosmic star formation rate density (CSFRD) equals that of all optical and near-infrared selected galaxies combined (Madau & Dickinson 2014). However, at z > 3 the situation is much less certain. A tail of submillimeter-selected galaxies have been confirmed beyond z > 4, but they trace only the very tip of the star formation rate (SFR) distribution at early times (e.g., Cooray et al. 2014; Strandet et al. 2017; Marrone et al. 2018). The total contribution of dust-obscured star formation, and therefore the census of star formation in the early universe, is unknown. Despite the strongly negative k-correction allowing sources to be found to z = 10, the overwhelming majority of (sub)millimeter-selected galaxies continue to be confirmed at z < 3 with Atacama Large Millimeter/submillimeter Array (ALMA) spectroscopy (Brisbin et al. 2017; Danielson et al. 2017). While some dusty galaxies have been discovered beyond z > 5 through gravitational lensing (Spilker et al. 2016; Zavala et al. 2018b), the lensing correction and selection effects make it challenging to establish their contribution to the CSFRD. Progress is hampered by the limited sensitivity and low spatial resolution of single-dish submillimeter observations and the difficulty of associating detections with counterparts in the optical–near-infrared (NIR). Ultradeep SCUBA surveys over moderate ∼100 arcmin2 are now pushing into the range of "normal" SFRs (several 100 M yr−1, main-sequence galaxies; e.g., Koprowski et al. 2016; Cowie et al. 2017, 2018) and extending to z > 4, but the analysis is often limited by the ability to identify counterparts at other wavelengths and derive accurate redshifts.

ALMA has opened an avenue for addressing this issue through surveys at superior sensitivity and spatial resolution. ALMA deep fields at ∼1 mm (Aravena et al. 2016; Dunlop et al. 2017) have probed to extremely deep flux limits over small areas (<5 arcmin2). Progress has still been limited likely because dust-obscured star formation preferentially occurs in massive galaxies (Whitaker et al. 2017), which are clustered and relatively rare (∼0.1 arcmin−2 at $\mathrm{log}({\text{}}M/{M}_{\odot })\gt 10.8$ and z ∼ 4; Davidzon et al. 2017). Wider (10's arcmin2) and shallower (∼100–200 μJy) ALMA surveys at ∼1 mm (Franco et al. 2018; Hatsukade et al. 2018) are now approaching large enough areas to identify tentative massive candidates at z > 4 (e.g., Yamaguchi et al. 2019).

A promising development is to select at longer wavelengths (>2 mm), which optimizes the selection to dusty star formation at redshift z > 4 (Béthermin et al. 2015; Casey et al. 2018). However, the current state-of-the-art ALMA Spectroscopic Survey (ASPECS) at 3 mm (González-López et al. 2019), still only covered ∼5 arcmin2, and identified 6 continuum sources, all at z < 3. Larger archival studies of ALMA 3 mm observations to find high-redshift candidates report some spectroscopic confirmations, but like the 1 mm redshift distribution, the majority lie at z < 3 (Zavala et al. 2018a). Recent advances with IRAM/GISMO provide evidence that 2 mm surveys favor selecting higher-redshift sources (Magnelli et al. 2019), although counterpart identification continues to be problematic due to large beam sizes. Overall, the number of strong candidates for dust-obscured sources at z > 4 remains small and as a result, the contribution of dust-obscured star formation in the early universe is poorly constrained.

Here we report the serendipitous discovery of 2 previously unknown sources in deep ALMA 3 mm observations in the COSMOS field, and assess the implications for dust-obscured star formation at z > 3. We assume a ΛCDM cosmology with H0 = 70 km s−1 Mpc−1, ΩM = 0.3, ΩΛ = 0.7, and a Kroupa (2001) initial mass function (IMF).

2. Methods

2.1. ALMA Millimeter Interferometry

The ALMA observations are part of a program targeting CO(2-1) line emission in z ∼ 1.5 galaxies (2018.1.01739.S, PI: Williams; 2015.1.00853.S, PI: Bezanson) in the COSMOS field. The data presented here include ALMA maps of five galaxies, and results on these targets are presented elsewhere (Bezanson et al. 2019; C. C. Williams et al. 2019, in preparation). In one of these pointings we identified two serendipitous continuum sources, unrelated to the original target, which we will now describe in detail.

ALMA Band 3 observations were carried out in two observing blocks on 2018 December 23 and 24 under program 2018.1.01739.S (PI: Williams). One 1.875 GHz spectral window was centered at the sky frequency (94.92 GHz) with 7.8 MHz (∼24 km s−1) channelization. Three additional 1.875 GHz bandwidth spectral windows were placed at sky frequencies 96.8, 106.9, and 108.8 GHz for continuum observations, each with 15.6 MHz channelization. A total of 43 antennas were active, reaching maximum baselines of 500 m, for an angular resolution of ∼2 arcsec (synthesized beam major axis = 2farcs2, minor axis = 1farcs7 under natural weighting). The total time on-source was 97 minutes.

The data were reduced using the standard ALMA Cycle 6 pipeline, and no problems with the pipeline calibration were found. The imaged data reach a continuum sensitivity of 5.7 μJy/beam at 101.9 GHz, and a typical line sensitivity of 55–65 μJy/beam per 100 km s−1 channel. We imaged the data using natural weighting, creating a 101 GHz continuum image of the field from all four spectral windows. We imaged the data using 0farcs2 pixels and created a 500 × 500 pixel image, yielding images 100'' on a side. Given the ALMA primary beam at this frequency (∼57'' FWHM), these images extend to approximately the 0.05 response point of the primary beam.

2.2. ALMA Source Detection

Interferometric maps without correction for the primary beam response have uniform, normally distributed noise properties across the field, and source detection significance is straightforward to measure from such maps. Two blind 3 mm continuum sources were apparent in this map, located 24farcs6 and 38farcs2 from the phase center, corresponding to primary beam response levels of 0.57 and 0.29, respectively. Each source is detected at a peak signal-to-noise ratio (S/N) of ∼8; the probability of finding a Gaussian noise fluctuation of this magnitude given the number of independent beams in the images is exceedingly low (<10−9). Both sources thus are real. The 3 mm flux densities of these sources, corrected for the primary beam response, are 155 ± 20 μJy and 75 ± 10, hereafter referred to as 3MM-1 and 3MM-2, respectively. Neither source is spatially resolved, based on a comparison of the peak pixel values and the integrated flux densities.

After finding both continuum sources in the combined map, we reimaged the upper and lower sidebands of the data separately in order to determine the spectral index of each source at these frequencies. Thermal dust emission on the Rayleigh–Jeans tail has a very steep power-law index with Sν ∝ ν2+β and β ∼ 1.5–2, while non-thermal synchrotron emission typically exhibits a negative spectral index, Sν ∝ ν−0.8. Given their respective frequencies, we expect dust emission to be ∼50%–60% brighter in the upper sideband than in the lower sideband. We find spectral indices for both continuum sources in excellent agreement with the expectation for thermal dust emission.

We also reimaged the data of both sidebands to search for blindly detected emission lines in each of the two continuum sources, using channel widths ranging from 100 to 400 km s−1. 3MM-2 contains a serendipitous emission line centered at 106.5 GHz, with an integrated flux density of 0.66 ±0.1 Jy km s−1. The spectrum is shown in Figure 2. A Gaussian fit indicates a width of 630 ± 70 km s−1. Assuming the line is a transition of carbon monoxide, the possible redshifts are z = [0.08, 1.16, 2.25, 3.33, 4.41, 5.49]. We find no significant emission lines in the spectrum of 3MM-1, and no evidence for a line at the same frequency of 3MM-2. A Gaussian fit restricted to the same frequency and width results in an integrated flux density of 0.17 ± 0.2 Jy km s−1, indicating no evidence for a line at that location.

Both 3MM-1 and 3MM-2 were also contained within the field of view of an additional ALMA Band 3 program, 2015.1.00861.S. Both sources were again far out in the primary beam of these data, at approximately the 0.25 and 0.1 response points, respectively. These data are described in more detail in Silverman et al. (2018), and have non-overlapping frequency coverage with our own data. We downloaded and imaged these data following the same procedure as for our own data. The new images reach a continuum sensitivity at a phase center of 10 μJy/beam at 93.5 GHz and a line sensitivity of 120–170 μJy/beam per 100 km s−1 channel, approximately a factor of two higher than that in our data. Neither source is detected in the continuum in these data, nor was expected to be detected given the sensitivity, effective frequency, and position of our sources within the ALMA primary beam. We additionally searched these data for blindly detected CO lines as in our own data, but found no significant emission lines. The limited (primary beam-corrected) sensitivity of these data preclude us from drawing strong conclusions about the redshifts of either source.

For convenience, we define an equivalent survey area as the total area across all five ALMA maps at which a source with the same S3 mm as 3MM-1 would be detected at >5σ. Taking into account the small variations in the central frequency of each map (tuned to the specific redshift of the target z ∼ 1.5 galaxies), which change the primary beam response shape and the detection threshold (25–31 μJy) due to the steep spectral index of dust emission, we derive a total survey area of 8.0 arcmin2 to 155 μJy (5σ limit).

2.3. Multiwavelength Photometry

3MM-1 has extremely deep coverage at all optical-to-submillimeter wavelengths, yet it has no counterpart within a radius of 3farcs3 in the deepest published multiwavelength catalogs in COSMOS to date, from 0.6 to 1100 μm (Le Floc'h et al. 2009; Aretxaga et al. 2011; Muzzin et al. 2013; Laigle et al. 2016; Geach et al. 2017; Hurley et al. 2017) or 3–1.4 GHz (Schinnerer et al. 2010; Smolčić et al. 2017). The astrometry of multiwavelength catalogs in COSMOS are excellent owing to the registration between VLA observations (e.g., Schinnerer et al. 2007), ground-based optical and space-based facilities, with an astrometric accuracy of 5 mas (Koekemoer et al. 2007). Similar astrometric accuracy was found between ALMA and COSMOS multiwavelength data sets (Schreiber et al. 2017). There is no apparent flux at optical-to-Spitzer/IRAC wavelengths at the ALMA position (Figure 1). 3MM-2 is also missing from these catalogs, likely because it is blended with a bright neighboring galaxy ∼1farcs3 to its north, with a photometric redshift 0.95 (see Figure 1), but appears detected in Ks and IRAC. It is possible that the ALMA source is simply a highly obscured region in the low-redshift galaxy, which can be ruled out by spatially deblended SED analysis. We therefore proceeded to perform deblended photometry on both ALMA positions using the following data sets.

Figure 1.

Figure 1. Cutouts (25'' × 25'') centered at a 3 mm position of 3MM-1 (red circle; 3'' diameter). 3MM-1 was not previously detected (>3σ) at any shorter wavelength, including deep optical and near-IR stacks, Spitzer, Herschel, and S2COSMOS SCUBA2 850 μm. Remeasuring with the ALMA position as a prior reveals marginal 2σ–3σ measurements in IRAC 3.6+4.5 and 850 μm, consistent with heavy dust obscuration at z > 4. It is faintly detected at 3 GHz (4σ), indicative of a moderate radio excess due to a possible active galactic nucleus (AGN). The 1.4 GHz image is excluded because neither source is significantly detected. 3MM-2 is also identified (blue circle), and is blended with a foreground galaxy at z = 0farcs95, 1farcs3 north (Muzzin et al. 2013; Laigle et al. 2016).

Standard image High-resolution image

2.3.1. Optical, Near-infrared, and Spitzer/IRAC

The optical data consist of the Subaru/Suprime-cam Bj, Vj, g+, r+, i+ and z+-imaging (Taniguchi et al. 2007), with 5σ limits of ∼25–27.4 mag in 1farcs2 apertures, and Subaru HyperSuprimeCam (HSC) g, r, i, z, and y (Aihara et al. 2018a, 2018b) imaging (∼25–26.8 mag, 5σ). Ultradeep NIR coverage is provided by the 4th data release of the UltraVISTA survey (McCracken et al. 2012), thanks to mosaics in the Y, J, H, and Ks filters to ∼25 mag (AB, 5σ). Remarkably, the coverage in the Ks band from DR4 is ∼0.9 mag deeper than that from DR3, allowing us to place strong constraints on the flux density of 3MM-1 at NIR wavelengths. Stacked images were constructed using the optical imaging 0.4–0.8 μm and the NIR imaging 0.9–1.6 μm. We use Spitzer/IRAC 3.6 μm and 4.5 μm mosaics that combine data from the S-COSMOS (Sanders et al. 2007) and the Spitzer Large Area Survey with HSC (SPLASH, PI: Capak) programs (∼24.5 mag, 5σ in 1farcs8 apertures), and the 5.8 μm and 8.0 μm from the S-COSMOS program (∼20.7 mag, 5σ in 1farcs8 apertures).

Figure 2.

Figure 2. Portion of the observed ALMA Band 3 spectrum covering the detected CO line in 3MM-2 at 106.5 GHz (upper sideband). The CO solution corresponding to CO(4-3) at z = 3.329 is in excellent agreement with the photometric redshift measured in Section 2.4. The line flux is 0.66 ± 0.1 Jy km s−1 and a Gaussian fit produces a width of 650 km s−1. No line is found at the same frequency in the spectrum of 3MM-1 (lower panel) with formal SNR = 0.8. No lines are detected in the lower sideband from 94 < ν < 96.8 GHz.

Standard image High-resolution image

We measured flux densities for 3MM-1 in the optical and UltraVISTA bands in 1farcs2-diameter apertures after subtracting the neighbors using mophongo (Labbé et al. 2013, 2015). This procedure carefully models the light profiles of the sources using a higher-resolution image as a prior, minimizing potential contamination by bright nearby objects (see Figure 3). In our analysis we adopted the HSC z-band image as a prior for 3MM-1, as it provides the best compromise between depth and resolution and the F814W-band image for 3MM-2, given the nearby bright neighbor at ∼1farcs3 distance. Because of the broader PSF, we adopted 1farcs8 apertures for our estimates in the IRAC bands. Total flux densities were then estimated from the spatial profile and the relevant PSF-correction kernel.

Figure 3.

Figure 3. Illustration of the deblending procedure using the mophongo software (Labbé et al. 2015). Deblending is performed by simultaneously fitting the pixels of all sources using the deep optical images and the ALMA positions as priors and accounting for differences in the PSFs. Top row: deblending results for the 4.5 μm band centered on the ALMA position of 3MM-2 (12'' × 12''). The left panel shows the original 4.5 μm image, the middle panel shows 3MM-2 after other modeled sources have been subtracted, and the right panel shows a 3.6–4.5 μm color image clearly indicating a vastly different IRAC color at the location of the ALMA source. Bottom row: higher contrast, zoomed-in panels centered on the ALMA position of 3MM-1, showing a faint AB ∼ 25.2 IRAC source after subtracting the PSF wings of a bright foreground neighbor.

Standard image High-resolution image

2.3.2. Far-infrared to Submillimeter

Far-infrared and submillimeter Herschel fluxes were measured for both galaxies by simultaneously fitting Gaussian profiles at fixed prior locations in the Herschel images specified by the ALMA locations and augmented by the MIPS positions of all neighboring objects from the S-COSMOS data (Sanders et al. 2007; Le Floc'h et al. 2009). The FWHMs of the Gaussians were 7farcs7, 12'' (PACS 100, 160 μm) and 18, 25, and 37' (SPIRE 250, 350, 500 μm), respectively. Uncertainties were computed by fitting Gaussians at random locations within a 2 arcmin radius and computing the rms. Flux calibration was performed by comparing to the 24 μm prior catalog of Herschel DR4 (Oliver et al. 2012; Hurley et al. 2017). Comparison of the measured SPIRE fluxes with those published showed a scatter of 30%; this was added in quadrature to the flux uncertainties for those bands. The correction is minor considering the S/N of ≲1 of the SPIRE fluxes. For SPIRE 500 μm only, the flux of a neighboring zspec = 1.45 galaxy 5'' to the south was subtracted separately prior to deblending, using its predicted IR emission based on the infrared SED of Wuyts et al. (2011) and its SFR24,IR = 65 M yr−1.

The procedure for measuring photometry at 850 μm follows that of Herschel/PACS and SPIRE (250 and 350 μm) using the ALMA and MIPS positions as priors. The sources lie in a region of shallower (σ850 ∼ 3 mJy/beam) coverage in the inhomogeneous SCUBA-2 Cosmology Legacy Survey (S2CLS; Geach et al. 2017) and are undetected in this map. We therefore use the state-of-the-art S2COSMOS observations, which achieved a σ850 ∼ 1.2 mJy over the entire field (Simpson et al. 2019).

2.3.3. Radio VLA 3GHz and 1.4 GHz

Neither source has a counterpart in either the VLA 3 GHz 5σ source catalog (Smolčić et al. 2017) or the 1.4 GHz deep survey 5σ source catalog (Schinnerer et al. 2010). Photometry on the 3 GHz map using the ALMA position as a prior reveals a 9.98 ± 2.39 μJy/beam (4σ) point source for 3MM-1, below the detection limit of the Smolčić et al. (2017) catalog. Blindly detected sources at this flux density have a high probability (>50%) of being spurious (Smolčić et al. 2017), therefore requiring the ALMA prior in order to be considered real. 3MM-2 shows no significant flux to 3σ limits of ∼7 μJy/beam. Within 10 arcmin2, the 1.4 GHz map has an rms of ∼17 μJy/beam and no significant flux (>3σ) is detected from either galaxy.

2.4. SED Modeling

The deep photometry from λ = 0.4–3000 μm places strong constraints on the SEDs and allows us to model their stellar populations and redshift. It should be noted that even if a source is not formally detected (>3σ) at most wavelengths, the absence of flux can still provide useful constraints, in particular on the redshift. All fitting is performed in linear fluxes and uncertainties and no upper limits are enforced on measurements with low S/Ns.

We use the Bayesian Analysis of Galaxies for Physical Inference and Parameter EStimation (bagpipes) code (Carnall et al. 2018), which assumes the stellar population synthesis models of Bruzual & Charlot (2003) and implements nebular emission lines following the methodology of Byler et al. (2017) using the cloudy photoionization code (Ferland et al. 2017). We select the flexible dust absorption model of Charlot & Fall (2000), with the exponent of the effective absorption ∝λn as a free parameter and adopt the Draine & Li (2007) dust emission model under the assumption of energy balance, such that dust-absorbed light is re-radiated in the far-infrared. All stellar populations have this effective absorption, while the youngest stars (defined as those with age <10 Myr) have an extra factor of attenuation applied (η) to account for dusty birthclouds.

The dust emission model is parameterized using the starlight intensity U incident on the dust (translating into a distribution of dust temperatures), the amount of PAH emission, and the fraction of dust at the coldest temperature. We further assume a delayed exponential star formation history with τ and age left free and metallicities ranging from 0.2 to 2.5Z. For each parameter, uniform, diffuse priors are assumed (see Table 1). In general, we do not expect to constrain most free parameters, but we will marginalize over those and limit the discussion to the more important parameters including redshift, stellar mass, and SFR. Where relevant, we will also quote results based on assumptions often used in the literature, such as assuming a fixed Calzetti et al. (2000) dust attenuation model.

Table 1.  Stellar Population Properties of 3 mm Continuum Sources

Parameter 3MM-1a 3MM-2 3MM-2 Range (Uniform Prior)
Attenuation Curve CF00 CF00 Calzetti
Coordinates 10:02:36.82 + 02:08:40.60 10:02:36.30 + 2:08:49.55  
Redshift ${5.5}_{-1.1}^{+1.2}$ zCO = 3.329 zCO =3.329 (0.0, 10.0)
SFR [M yr−1] ${309}_{-149}^{+241}$ ${247}_{-76}^{+77}$ ${112}_{-38}^{+27}$  
Log10 Stellar Mass [M] 10.8${}_{-0.4}^{+0.4}$ 11.1${}_{-0.2}^{+0.2}$ 10.2${}_{-0.2}^{+0.2}$ (6.0, 14.0)
Mass-weighted age [Gyr] ${0.2}_{-0.1}^{+0.1}$ ${0.5}_{-0.2}^{+0.1}$ ${0.1}_{-0.1}^{+0.2}$  
AV ${4.0}_{-1.0}^{+1.4}$ ${2.7}_{-0.3}^{+0.3}$ ${1.4}_{-0.1}^{+0.1}$ (0.0, 10.0)
ηb ${2.0}_{-0.6}^{+0.6}$ ${2.1}_{-0.8}^{+0.6}$ ${2.5}_{-0.5}^{+0.3}$ (1, 3.0)
Uminc 13.0${}_{-7.3}^{+7.7}$ ${7.6}_{-4.0}^{+6.3}$ ${4.8}_{-2.4}^{+3.4}$ (1, 25.0)
gammad ${0.5}_{-0.3}^{+0.3}$ ${0.5}_{-0.3}^{+0.3}$ ${0.5}_{-0.3}^{+0.3}$ (0.01, 0.99)
ne ${0.9}_{-0.4}^{+0.6}$ ${0.4}_{-0.1}^{+0.1}$ (0.1, 2.0)
qpahf ${2.3}_{-1.1}^{+1.1}$ ${3.3}_{-1.0}^{+0.5}$ ${3.0}_{-1.1}^{+0.7}$ (0.5, 4.0)
Log10 LIRg [L] 12.6 12.4 12.1
Gas mass [M] 0.5–1.5 × 1011 4–9 × 1010    

Notes.

aFitting priors are uniform with the range as defined in the last column. bMultiplicative factor producing extra attenuation for young stars. cStarlight intensity on dust grains, related to dust temperature as Tdust ∝ U1/6 (Draine & Li 2007). dFraction of stars at Umin. ePower-law slope of the dust extinction law (for Charlot & Fall 2000). fPAH mass fraction. gTotal infrared luminosity (8–1000 μm) in units from integrating the median posterior spectrum in Figure 4.

Download table as:  ASCIITypeset image

3. Results

The photometry for both galaxies is presented in the Appendix (Table 2).

3.1. 3MM-1

3MM-1 is undetected (<3σ) in individual optical (>27AB), near-infrared (>25.2), and Spitzer/IRAC (>25.2) bands, and remains undetected in the stacked ultradeep optical (>29.6 AB) and near-IR (>27.9 AB) data (1σ limits). The source is faintly detected (25.2 AB ∼3σ) in deep stacked IRAC 3.6+4.5 observations, indicating the source is likely real with extremely red colors. Owing to the shallower depth, the source is undetected in deep Spitzer/MIPS 24 μm and Herschel PACS+SPIRE 160–500 μm. The flux measured with SCUBA is S850 = 3.5 ± 1.1 mJy (SNR = 3.1). This flux density is fainter than the depth typically achieved by deep single-dish submillimeter surveys (>3.5–5 mJy; e.g., Weiß et al. 2009; Aretxaga et al. 2011; Cowie et al. 2017; Geach et al. 2017) at robust detection thresholds (≳4σ). The observed SED is shown in Figure 4 (left panel). For clarity, photometry with S/N < 1 is indicated with arrows at the 1σ rms level, and S/N > 1 with 1σ error bars. The deep non-detections at all wavelengths between 24 and 500 μm and the extreme flux ratios (S24/S850 < 4 × 10−3, S24/S3 mm < 9 × 10−2, S4.5/S850 <9 × 10−5, S4.5/S3 mm < 2 × 10−3; see, e.g., Cowie et al. 2018; Yamaguchi et al. 2019) demand that the dust peak be highly redshifted, z > 4.

Figure 4.

Figure 4. Left: observed photometry of 3MM-1 (points). For display purposes, we show 1σ upper limits (downward arrows) for data with S/N < 1 and stacked optical and NIR fluxes (horizontal bars) instead of the individual optical/NIR non-detections. Any photometry with SNR > 1 is plotted with 1σ uncertainties but does not necessarily indicate a significant detection. The light blue squares indicate photometry with S/N < 3, and the dark blue circles indicate detections with S/N ≥ 3. Shown are the median posterior spectrum (dark orange) and 16–84th percentile range (light orange) from Bayesian SED fitting with bagpipes assuming the Charlot & Fall (2000) dust model. Insets contain the posterior redshift distribution and the 16, 50, and 84th percentiles (black lines). The deep FIR photometric limits at 24 < λ < 250 μm favor high redshift (90% probability at z > 4.1). The dotted line indicates the radio spectrum ∝ν−0.8 expected from the total LIR (Tisanić et al. 2019), suggesting a 3 GHz radio excess. Right: SED of 3MM-2. As in the left panel, all photometry with S/N > 1 are plotted as error bars with 1σ uncertainties. The photometric redshift agrees with CO(4–3) at zCO = 3.329 (red line). The SED fit (orange) assumes the Charlot & Fall (2000) dust model and the redshift is fixed to zCO. Also shown is a fit assuming a Calzetti et al. (2000) starburst attenuation law (blue curve), a common assumption in high-redshift studies that can drastically change the estimated stellar mass, LIR, and SFR. The blue curve has a 8× lower stellar mass and 2× lower SFR while still adequately fitting the observations.

Standard image High-resolution image

Fitting stellar population models, we find the observations are well reproduced by a massive 1010.8±0.4M, star-forming $\mathrm{SFR}={309}_{-149}^{+241}$M yr−1, highly obscured ${A}_{V}\sim {4}_{-1.0}^{+1.4}$ galaxy at very high redshift $z\sim {5.5}_{-1.1}^{+1.2}$. The SED-fitting results are shown in Figure 4 and the posterior values and priors are listed in Table 1. The posterior distributions of all parameters are presented in Figure 7. Adopting a classical Calzetti dust law would increase the stellar mass, SFR, and redshift by 60%, 10%, and 5%, respectively. These changes are within the uncertainties, so we elect to adopt our more conservative values measured assuming Charlot & Fall (2000) as fiducial. We also experimented with fitting only the mid-infrared to submillimeter SED, finding the same results. In all cases, the solutions require the galaxy to be at high redshift, with z > 4.1 at 90% confidence and a total infrared luminosity LIR = 4 × 1012L (based on integrating the median posterior spectrum). Adopting an SED typical for the most obscured AV > 3.5 SMGs (da Cunha et al. 2015) produces similar redshift and LIR, as does adopting the Arp 220 SED. Converting to SFR using the Kennicutt & Evans (2012) conversion, which does not implicitly assume a SFH, results in a larger estimated SFR = 540M yr−1, but consistent within the 1σ uncertainties from SED fitting.

Calculating redshift using the radio-to-submillimeter spectral index method (using 850 μm and the upper limit on 1.4 GHz) implies a lower limit to the photometric redshift of ∼4.7 (e.g., Carilli & Yun 1999). We note that this is very similar to the massive dusty galaxy HDF850.1, which had a redshift 4.1 according to this relation (Dunlop et al. 2004), and was later spectroscopically confirmed at z = 5.2 (Walter et al. 2012). The SED of 3MM-1 is very similar to that of HDF850.1; scaling its best-fit SED to the 850 μm and 3 mm fluxes of 3MM-1 predicts its observed 3 GHz radio flux of 10 μm (Figure 5).

Figure 5.

Figure 5. Comparison between the SED of 3MM-1 with that of both HDF850.1 and Arp220. The SED-fitting results for 3MM-1 are the same as those in the left panel of Figure 4. The SED of Arp220 scaled to the 850 μm and 3 mm photometry of 3MM-1, at the photometric redshift of 3MM-1, is remarkably similar to the 3MM-1 SED. Forcing Arp220 to the redshift of 3MM-2 violates the Herschel FIR+sub-mm constraints. The SED of 3MM-1 is very similar to that of the bright zspec = 5.18 dusty galaxy HDF850.1 (using the SWIRE SED of IRAS 20551-4250 as in Serjeant & Marchetti 2014) scaled to match the average observed 850 μm and 3 mm flux).

Standard image High-resolution image

The observed 3 GHz radio flux (4σ) is in excess of that expected from empirical SED templates for obscured star-forming galaxies (e.g., SMGs; da Cunha et al. 2015) and recent calibrations of the redshift-dependent ratio of total infrared luminosity to rest-frame 1.4 GHz luminosity density (Delhaize et al. 2017; Tisanić et al. 2019). For typical assumptions of the radio spectral slope Sν ∝ ν−0.8 (dashed line in Figure 4), the radio emission of 3MM-1 is a factor of ∼6 ± 1.5 higher than expectations for star formation, in excess of the commonly adopted ×3 threshold for AGN activity (Daddi et al. 2009). No evidence for redshift evolution in the spectral slope, or dependence on sSFR or distance from the main sequence has been reported (Magnelli et al. 2015).

The radio excess therefore indicates plausible evidence for the presence of an AGN. We note that no X-ray counterpart exists within 5'' (Civano et al. 2016), although this is not surprising given the high redshift of 3MM-1. Given the possible presence of an AGN, it is worth noting that inferred parameters from SED fitting could be biased by AGN emission, in particular in the mid-infrared (e.g., Leja et al. 2018). Therefore, we investigate any possible contamination using the empirical AGN-dominated template published by Kirkpatrick et al. (2015). We use their mid-infrared template based on galaxies with the largest mid-infrared luminosity contribution from AGN emission scaled to the 24 μm flux of 3MM-1 and subtract the predicted AGN contribution from the observed SED. This represents the maximum AGN contribution that could be accommodated by the data without violating the observed photometry. We then remeasure the SED fitted parameters. Subtracting this empirical AGN template reduces the SFR by ∼15%, well within the uncertainty, but does not change the inferred stellar mass significantly. We conclude that any AGN contribution does not significantly impact our SED-fitting results.

Finally, the 3 mm flux density enables an estimate of the molecular gas mass for 3MM-1 following the calibration of submillimeter flux to gas mass (Scoville et al. 2016). Exploring the range of the photometric redshift PDF, and typical assumptions about the dust emissivity, temperature, and dust-to-gas ratio, we find that 3MM-1 likely has Mgas ∼ (0.5–1.5) ×1011 M. This is independent evidence that the galaxy is massive, with a high inferred gas fraction (∼60%).

3.2. 3MM-2

3MM-2 is optically faint, but detected in Ks,AB = 24.2 and relatively bright in IRAC ∼23 mag, with red optical-to-IRAC colors. There is an apparent break between the Js and H-band, consistent with a Balmer/4000 Å break at z ∼ 3. The source is undetected at 24–870 μm, but with less extreme flux ratios (S24/S3 mm < 3.5 × 10−1, S4.5/S3 mm < 3 × 10−2) compared to 3MM-1, indicating a lower redshift.

We follow the same procedure as before to fit the SED, finding a well-constrained photometric redshift of z = 3.3 ± 0.2 (99% probability between 2.7 < z < 3.7), ruling out that this source is simply a highly obscured region in the nearby z = 0.95 foreground galaxy to the north. Considering the CO emission line detection (consistent with redshifts 2.25, 3.33, and 4.41), we determine the line is likely CO(4–3) at z = 3.329, in agreement with the best-fit photometric redshift. Fixing the redshift to the spectroscopic redshift, the observations are best reproduced with a massive 1011.1M, star-forming SFR = 250 M yr−1, highly obscured AV ∼ 2.7 galaxy. Using the Kennicutt (1998), Kennicutt & Evans (2012) conversion from LIR to SFR yields ∼375 M yr−1. The derived gas mass using the CO(4-3) line flux (assuming the average CO excitation from Bothwell et al. 2013) is in the range ∼(1–4) × 1011 M, given the factor of four scatter in measured excitation for high-redshift dust-obscured galaxies (e.g., Figure 1 in Narayanan & Krumholz 2014). The 3 mm derived gas mass is (4–9) × 1010 M. Both estimates indicate gas fractions in the range 30%–70%. The SED properties for 3MM-2 are listed in Table 1, the SED fit is shown in Figure 4, and posterior distributions for the parameters are presented in Figure 8.

The best-fit power-law index for the Charlot & Fall (2000) dust model is n = 0.4 ± 0.1, flatter than the n = 0.7 appropriate for typical nearby starburst galaxies. The flat spectral index means that the attenuation curve is "grayer" than a Calzetti attenuation law (e.g., Chevallard et al. 2013), resulting in larger attenuation for the same amount of reddening E(BV) = 0.3. The assumed attenuation model can have a large impact on the derived stellar mass. If instead of a flexible attenuation model a classical Calzetti starburst law is assumed then the fits result in a factor ${8}_{-3.9}^{+7.1}$ lower stellar mass (1010.2 M) and factor 2.2 ± 0.7× lower SFR (same answer if the infrared-to-millimeter constraints are also removed). The large difference in stellar mass is remarkable. The reason is that the steeper reddening curve of Calzetti reproduces the data with an intrinsically blue, young, low M/L ratio stellar population, driving the lower masses. The Charlot & Fall (2000) model requires an older, redder, higher M/L population, while also implying larger attenuation at all wavelengths, including the near-infrared. We note that these results all assume energy balance and a parametric delayed–τ SFH. We explore the dependence on SFH by also fitting the SED with prospector (Leja et al. 2017), which is capable of fitting a non-parametric SFH in logarithmic bins of age (Leja et al. 2019). These fits produce qualitatively similar results, with Calzetti yielding ∼2× lower SFR and stellar mass.

An independent estimate of mass can be derived by estimating the dynamical mass from the 630 km s−1 width of the CO emission line at 106.5 GHz (Figure 2). Following the procedure used by Wang et al. (2013), we compute Vcirc =0.75FWHM(CO)/sin(i) and ${M}_{\mathrm{dyn}}=1.16\times {10}^{5}{V}_{\mathrm{circ}}^{2}D$, where we note that inclination angle i and disk diameter D (in kpc) are unknown. Adopting the mean size ${R}_{e}={3}_{-0.9}^{+1.8}\,\mathrm{kpc}$ expected for massive 1011M star-forming galaxies at z = 3 (van der Wel et al. 2014), we find ${M}_{\mathrm{dyn}}{\sin }^{2}(i)={1.5}_{-0.5}^{+1.0}\times {10}^{11}$M. Assuming a disk geometry with an inclination angle of i = 45° would imply Mdyn ∼ 3 × 1011M. Overall, the high dynamical mass appears to agree better with the high stellar mass inferred from the Charlot & Fall (2000) dust model with a flat n than with the Calzetti-based stellar mass, but the large uncertainties in dynamical and stellar mass prevent firmer conclusions.

3.3. Redshifts

Given the close r ∼ 15'' separation on the sky of the two ALMA sources we consider the possibility that they are at the same redshift. The redshift of 3MM-2 (z = 3.329) is strongly disfavored for 3MM-1 (1% probability z < 3.3). Alternatively, our identification of zCO(4–3) = 3.329 for 3MM-2 is erroneous and the galaxy is at zCO(5–4) = 4.41, more consistent with 3MM-1. This is unlikely based on the SED fit of 3MM-2, which is well constrained by the presence of both a Lyman Break and Balmer break. We find no significant emission line (S/N > 3) in the spectrum of 3MM-1, which suggests it is unlikely at 3.20 < z < 3.35 and 3.72 < z < 3.90. A Gaussian fit to the spectrum of 3MM-1 at fixed 106.5 GHz indicates SNR = 0.8, providing no evidence for 3MM-1 and 3MM-2 occupying the same dark matter halo (DMH). Additionally, we scan the nearby velocity space in case both galaxies are in a large-scale structure filament, but we find no emission line of S/N > 2 within Δv = ±2000 km s−1. Finally, we consider the odds of finding two 3 mm sources within r ∼ 15'' on the sky. Using the 3 mm number counts of Zavala et al. (2018a), we simulate uniformly random distributions of two sources, which indicate a ∼10% chance of finding a 75 μJy source within r < 15'' of a 155 μJy source. These odds are low, but not negligible. Overall, there is no conclusive evidence that the two sources are at the same redshift. We therefore proceed and take the zphot = 5.5 ± 1.1 for 3MM-1 at face value.

4. Discussion

4.1. Implications from Full Optical-to-millimeter Spectrum SED Fitting

We consider the SED-fitting analysis using a stellar population model with self-consistent dust absorption and emission under the assumption of energy balance and constrained by high-quality optical-to-millimeter observations. In the case of the z ∼ 5.5 source 3MM-1, which lacks strong detections at any wavelength other than 3 mm, it is remarkable that the Bayesian posterior probability distributions for key parameters such as redshift, stellar mass, dust attenuation, and SFR are as constrained as they are (see Table 1 and Figure 7). Closer inspection indicates that the results are mostly driven by the combination of the high S/N ALMA measurement with deep photometric limits in the short-wavelength infrared (Spitzer/MIPS and Herschel/PACS 100–160 μm), which demand that the dust emission peak is highly redshifted (z > 4) and the SFR is high. It is likely that the uncertainties on the stellar mass and other stellar population parameters are underestimated, because of the tight prior imposed by the parametric SFH (Carnall et al. 2019; Leja et al. 2019). We note that our choice of a rising SFH produces a relatively conservative (lower) estimate of stellar mass compared to constant or declining SFH.

There is accumulating evidence that the attenuation law in very dusty galaxies can be flatter than the Calzetti law (Salmon et al. 2016; Leja et al. 2017; Salim et al. 2018), possibly caused by a more uniformly mixed geometry of both old and young stars with dust (Narayanan et al. 2018). The flat inferred attenuation law λ−0.4 for 3MM-2 is not surprising in that regard. It is notable and sobering, however, that modeling the optical-to-millimeter with a classical Calzetti starburst law instead (keeping everything else the same) results in significantly lower SFR and stellar mass (1010.2), in apparent tension with the high dynamical mass ((1.5–3) × 1011M). Detailed analyses of ULIRGs at z ∼ 2 with photometric constraints on far-infrared dust emission (e.g., Lo Faro et al. 2017) also found that stellar masses inferred using a Calzetti law are systematically lower because of the smaller amount of reddening at near-infrared wavelengths. Clearly, caution should be exercised when applying a single locally calibrated attenuation law at high redshift. Overall, the results highlight the need for deep (sub-)millimeter measurements to determine bolometric luminosities and provide high-redshift empirical constraints on the dust law/energy balance (e.g., Hodge et al. 2016).

For the remainder of the paper we emphasize discussing the implications of the ALMA-only source 3MM-1, which likely represents a population that is absent from all current optical-IRAC selected galaxy studies, is below the nominal detection threshold of deep single-dish sub-mm surveys, and therefore deserves more scrutiny.

4.2. Number Counts

A recent unbiased ALMA 3 mm archival search for continuum sources over 130 independent pointings (Zavala et al. 2018a) enabled the first 3 mm number counts. For archival searches, however, detection limits and effective search area are necessarily very inhomogeneous due to the strong variation of the ALMA primary beam response, complicating a straightforward comparison of results. Restricting the comparison to 3MM-1, we apply their best-fit power law to the cumulative number counts, which considers effective selection area and incompleteness, predicting one source brighter than S3mm > 155μJy per ${16}_{-8}^{+16}$ arcmin2. This is consistent with our finding one source in our effective survey area of 8 arcmin2 (see Section 2.2).

We note that given the estimated stellar mass of ∼1010.8M at z = 5.5, we expect 3MM-1 to occupy a massive DMH (∼1012M) and to be strongly clustered. While our results are dominated by Poissonian uncertainties, this may imply that the uncertainties in published number counts are underestimates. Clustering may affect source counts from studies over very small areas such as the ASPECS 4.6 arcmin2 survey in the HUDF (González-López et al. 2019), and even the counts in the multiple pointings of Zavala et al. (2018a), as these observations were generally targeting moderately sized deep fields (COSMOS, CDFS, UDS). Clustering effects may be exacerbated by the fact that dust obscuration preferentially occurs in massive galaxies (Whitaker et al. 2017). We note that our results are unlikely to be biased by the original primary targets, which were all at low redshift (z ∼ 1.5). Clearly, larger blind surveys are needed to constrain number densities at the bright end.

4.3. Unbiased ALMA 3 mm Selection of High-redshift Galaxies

Simulations of dust-obscured galaxies in the early universe predict 3 mm continuum selection optimizes the sensitivity to DSFGs at redshift z > 4 (Casey et al. 2018). Current evidence is still mixed: continuum-detected faint galaxies in the small ultradeep ALMA ASPECS field are at average $\bar{z}=2.3$ (González-López et al. 2019), while limited ground-based spectroscopic evidence for sources found in wider-area archival data indicates $\bar{z}=3.1$ (Zavala et al. 2018a). Note, however, that the luminosity of the sources and volume covered could impact the redshift selection function (Strandet et al. 2016; Brisbin et al. 2017), and follow up ground-based spectroscopy in the optical and near-IR is possibly biased against highly obscured sources at z > 3, due to the faintness of Lyα, nebular optical lines, and a lack of wavelength coverage >2.4 μm. Overall, the redshifts of the two continuum sources in this study, a probable zCO = 3.329 (based on a single CO line and congruous zphot = 3.25 ± 0.15) and zphot = 5.5 ± 1.1 ($\bar{z}=4.4$), are consistent with 3 mm favoring higher redshifts than the $\bar{z}=2.2$ typical for 870 μm selected galaxies (Simpson et al. 2014).

A challenge in determining the selection function is the difficulty of identifying counterparts and determining redshifts. Neither of our two 3 mm sources had counterparts in previous deep optical-to-radio selected catalogs, raising some concern for analyses where this is a critical step (e.g., photometric redshifts or spectroscopic follow up in optical/NIR). This is particularly problematic for single-dish sub-mm observations, due to the large beam size, but here we find it to be challenging even with the high spatial resolution FWHM = 2'' and accurate astrometry offered by ALMA. In the case of 3MM-2 the source is very close to a bright foreground galaxy (1farcs3 to the north), which was initially mistaken as the counterpart. The optical/NIR faintness and low resolution of the Spitzer/IRAC data caused it to be missing or blended in existing multiwavelength catalogs. In the case of 3MM-1 the combination of high redshift and high obscuration resulted in it being extremely faint and undetected at all optical-infrared wavelengths. It is therefore possible that the identification of the highest-redshift galaxies in existing (sub-)mm selected samples is still incomplete. In addition, obtaining spectroscopic redshifts in the optical/NIR is unfeasible for the faintest sources. ALMA spectral scans targeting CO and [C ii] are likely the only recourse until other facilities (e.g., JWST/NIRSpec, LMT Redshift Search Receiver) become available.

4.4. Comparison to Other "dark" Galaxies

The observed SED of 3MM-1 is different from that of other optical/near-IR/IRAC-faint populations. The ALESS survey identified 9 IRAC-faint SMGs to 1 magnitude brighter IRAC limits (Simpson et al. 2014), but these sources were much brighter in the far-infrared (peaking at 250–350 μm to ∼10 mJy), with a significantly bluer dust peak consistent with z ∼ 2−3. Such an infrared SED is ruled out by our observations. 3MM-1 is fainter at all wavelengths than the so-called H/K−dropout galaxies (e.g., Caputi et al. 2012; Wang et al. 2016; Schreiber et al. 2017), which are generally selected on H − 4.5 color, >10× brighter in IRAC, and estimated to be at lower redshift $\bar{z}=3.7$.

More recently, SCUBA surveys have been carried out to unprecedented depth over CANDELS fields (∼1.5 mJy detection limit; e.g., Cowie et al. 2017, 2018; Koprowski et al. 2017; Zavala et al. 2017; Simpson et al. 2019). Such surveys are deep enough to discover S870 = 2–4 mJy objects like 3MM-1 at z > 4. Indeed, the deepest SCUBA surveys have reported submillimeter selected objects that are likely at z > 4 and that are either extremely faint at optical-IRAC and MIPS wavelengths (Cowie et al. 2018) or likely initially misidentified and determined to be at high redshift based on the long-wavelength IR SED (e.g., Koprowski et al. 2016). With optical-IRAC magnitudes in the range 29–25 AB, no MIPS detection, and only a marginal 3 GHz radio detection, 3MM-1-like objects would have likely defied secure identification, and its probable z = 5–6 redshift could have only been determined from the IR SED.

More similar are recently reported ALMA 1.2 mm selected sources from the ASAGAO survey (Yamaguchi et al. 2019) in GOODS-S. These sources also lack obvious counterparts at shorter wavelengths, show extremely small S4.5/S1.2mm and S24/S850 flux ratios (indicative of high redshift), but are generally fainter S1.2 mm = 0.44–0.8 mJy than 3MM-1 (S1.2 mm ∼ 2.3). Additional candidates that lack deep multiwavelength ancillary data may exist outside legacy fields (Wardlow et al. 2018).

4.5. Contribution to CSFRD

The fact that 3MM-1 was identified in a survey of such small area suggests that similar galaxies are common in the early universe. Using the effective area of 8 arcmin2 as derived in Section 2.2, our finding implies a source density of ${0.13}_{-0.10}^{+0.30}$ arcmin−2, an order of magnitude higher than the rare submillimeter-selected starbursts at z > 4 (0.01–0.02 arcmin−2; e.g., Danielson et al. 2017; Marrone et al. 2018). If our results are representative, galaxies like 3MM-1 could represent the "iceberg under the tip" of the known extreme dust-obscured star-forming galaxies in the early universe.

Estimating the space densities and contribution to the CSFRD requires estimating a selection volume, which is difficult as we do not know the expected properties or true abundance of DSFGs at z > 4. Instead, we derive a simple estimate based on the derived properties of 3MM-1. Encouraged by the strong redshift constraints from the long-wavelength photometry, we set the lower bound of the selection volume to the lower 10th percentile redshift probability (z > 4.1). The upper bound is chosen to be the redshift beyond which we expect only 0.1 halos of sufficient mass based on the cumulative halo mass function. Using the halo mass function calculator HMF published by Murray et al. (2013) and assuming the halo mass function of Behroozi et al. (2013b) and a high 30% baryon conversion into stars and (converting M* = 1010.8M into a conservative Mhalo ∼1012M) we compute this upper bound to be z = 5.7.

Using this selection volume we determine a space density of ${2.9}_{-2.4}^{+6.5}\times {10}^{-5}\,{\mathrm{Mpc}}^{-3}$. We find that the contribution to the CSFRD by 3MM-1 is ρSFR ${0.9}_{-0.7}^{+2.0}\times {10}^{-2}$ M yr−1 Mpc−3 (converted to a Chabrier 2003 IMF for comparison with literature measurements). Assuming instead the cumulative space density of massive halos (>1012) at z ∼ 5 using the Behroozi et al. (2013b) halo function, produces a very similar ρSFR 0.6 × 10−2 M yr−1 Mpc−3 for our derived SFR.

With only a single object, the Poissonian uncertainty is large and dominant (Gehrels 1986). Cosmic variance is only on the order of ∼30% based on the calculator by Trenti & Stiavelli (2008), and is therefore not further included. No completeness corrections are applied, because the true distribution is unknown. Figure 6 shows comparisons to various results at 0 < z < 10 from literature that report dust-uncorrected UV-derived SFR and the dust-obscured SFR (IR-to-millimeter derived).

Figure 6.

Figure 6. The cosmic star formation history of the universe. Blue shades are UV (dust uncorrected) and red shades are IR-to-millimeter derived SFRs. We add to the compilation of Madau & Dickinson (2014; blue circles) more recent z > 4 measurements by Finkelstein et al. (2015), Bouwens et al. (2016), McLeod et al. (2016), and Oesch et al. (2018). Red circles are the dust-obscured IR compilation of Madau & Dickinson (2014), to which we add recent z > 2 measurements from Swinbank et al. (2014), Koprowski et al. (2017), Magnelli et al. (2019), Cowie et al. (2018), Dunlop et al. (2017), and Liu et al. (2018). The contribution of 3MM-1 is indicated by the black star and shaded region, where the redshift range indicates the estimated selection volume as discussed in Section 4.5.

Standard image High-resolution image

The contribution of 3MM-1 is higher than inferred for the two near-infrared dark ALMA 1.2 mm sources in Yamaguchi et al. (2019), mostly owing to their smaller implied total infrared luminosities. This study does not report formal uncertainties or derive redshifts, which prohibits a more quantitative comparison. Bright SMGs beyond z > 4 contribute ∼1 × 10−3 M yr−1 Mpc−3 (Swinbank et al. 2014; Michałowski et al. 2017), about an order or magnitude lower than our best estimate. Results from fainter SMGs found in the deepest SCUBA surveys, with luminosities similar to 3MM-1 (e.g., Koprowski et al. 2017; Cowie et al. 2018), show a declining contribution at 2 < z < 5, consistent with our estimates.

Interestingly, if 3MM-1 is representative, a population with similar properties could contribute as much to the CSFRD as all known ultraviolet-selected galaxies at similar redshifts combined. This could even imply that dust-obscured star formation continues to dominate the cosmic star formation history beyond z > 4, where current infrared-to-mm measurements are incomplete.

4.6. Contribution to Stellar Mass Density

The high stellar mass 1010.8 M and large space density ∼3 × 10−5 Mpc−3 of 3MM-1 imply a considerable cosmic stellar mass density (CSMD) in similar objects at z ∼ 5: ${\rho }^{* }={1.9}_{-1.5}^{+4.4}\times {10}^{6}$ M Mpc−3, higher than reported for bright (S850 > 4 mJy) submillimeter galaxies (≈0.5 × 106 M Mpc−3; Michałowski et al. 2017). Comparing to the estimate of the CSMD based on HST-selected galaxies ρ*(z = 5) = 6.3 ×106M Mpc−3 (Song et al. 2016), suggests that they could contribute a significant fraction (${22}_{-16}^{+25} \% $) to the total and perhaps even dominate the high-mass end. The high-mass end of the galaxy stellar mass function at high redshift z > 4 is still uncertain and susceptible to biases. Current estimates for the number density of ∼1010.8 M galaxies at z ∼ 5 are (1–5) ×10−5 Mpc−3dex−1 (Duncan et al. 2014; Grazian et al. 2015; Song et al. 2016; Davidzon et al. 2017; Stefanon et al. 2017), comparable to the space density derived for our ALMA-only galaxy 3MM-1. Given that optical-IRAC dark galaxies are missing from these previous studies, it is therefore possible that about half the stellar mass density in high-mass galaxies at z ∼ 5 remains unaccounted for.

4.7. Implications for the Formation of Massive Galaxies

Bright submillimeter-selected galaxies at z > 3 with SFR ≳ 1000 M yr−1 are often hypothesized to be progenitors of massive z ∼ 2 quiescent galaxies (e.g., Toft et al. 2014; Spilker et al. 2018). More recently, massive quiescent galaxies have been spectroscopically confirmed at 3 < z < 4 (Glazebrook et al. 2017; Schreiber et al. 2018), but finding their progenitors at even earlier times is challenging. Generally, UV-selected galaxies at z > 4 are presumed to not be abundant, massive, and star-forming enough to produce the population of the earliest known massive quiescent galaxies with N ∼ (3−5) × 10−5 Mpc−3 and Log (M/M) ≳ 10.6 (Straatman et al. 2014). Bright (>4 mJy) submillimeter galaxies are a possible avenue, but it is difficult to establish a conclusive connection. Estimated SMG number densities at z > 4 are low and uncertain ∼0.1–3 × 10−6 Mpc−3 (Ivison et al. 2016; Michałowski et al. 2017; Jin et al. 2018). While their high SFRs indicate they will rapidly form massive galaxies, the modest inferred gas masses indicate gas depletion timescales on the order of 10–100 Myr (e.g., Aravena et al. 2016; Spilker et al. 2018), and large duty cycle corrections are needed to make up for the low number densities.

The inferred space densities and stellar mass of 3MM-1, on the other hand, are already as large as those reported for the population of early quiescent galaxies at 3 < z < 4 (Nayyeri et al. 2014; Straatman et al. 2014). The large ∼1011M gas mass and modest SFR indicates much longer depletion timescales (∼200–500 Myr), half the age of the universe at this epoch, and implies a ∼50%–100% duty cycle for our adopted 4 ≲ z  ≲ 6 selection window. 3MM-1 could therefore represent a more gradual path for massive galaxy growth compared to a rapid and bursty formation, as has been found in some bright merger-induced SMGs (e.g., Marrone et al. 2018; Pavesi et al. 2018). Overall, our results provide evidence for the existence of a sustained growth mode for massive galaxies in the early universe.

Finally, the large systematic difference in SFR and stellar mass in particular depending on the assumed attenuation model for 3MM-2 suggests exercising caution when relating progenitors and descendants galaxies. Often, such links are determined based on the capability of a progenitor population to produce adequate numbers of sufficiently massive galaxies at later times (e.g., Straatman et al. 2014; Toft et al. 2014; Williams et al. 2014, 2015), or by the assumption that the rank-order on stellar mass can be reliably determined (e.g., Brammer et al. 2011; Behroozi et al. 2013a; Leja et al. 2013). Such determinations could be more uncertain than has been accounted for if the derived stellar masses of individual galaxies are off by factors of ∼8 (see Section 3.2). This is a particular concern at high redshift z > 3, where very few galaxies have high enough S/N(sub-)millimeter observations for meaningful constraints on the dust attenuation model. These results are obviously based on only a single galaxy and larger samples with full optical-to-millimeter photometry are needed to determine the scatter in stellar mass.

4.8. Future ALMA and JWST Observations

The dark nature of the 3MM-1 non-detection in very deep stacked optical and near-IR data, and the faint IRAC fluxes, suggest a prominent population of DSFGs at z > 4. The high inferred redshift may point to the efficacy of blind surveys at 2–3 mm at finding the earliest dusty star-forming galaxies (Béthermin et al. 2015; Casey et al. 2018). These galaxies are below the classical detection limit (>4 mJy) of bright galaxies found in single-dish submillimeter surveys. Deep 2 mm single-dish and wider 1–3 mm ALMA (sub)-millimeter surveys are only just starting to push into this territory at very high redshift (e.g., Magnelli et al. 2019; Yamaguchi et al. 2019). Until the launch of the James Webb Space Telescope (JWST), ALMA alone can find and study these galaxies.

Wide-area unbiased ALMA surveys covering hundreds of arcmin2 are necessary to further constrain their prominence in the early universe. Such surveys are feasible at 2–3 mm with ALMA because of the relatively large size of the primary ALMA beam and the exquisite sensitivity to high-redshift star-forming galaxies even in short integration times (Casey et al. 2019). To date, none of the ALMA-only galaxies found have been spectroscopically confirmed, but spectral line scans for CO and in particular [C ii] are efficient and feasible with ALMA. Future surveys with JWST will enable systematic studies of large samples of faint SMG galaxies. Large legacy surveys such as the JWST Advanced Deep Extragalactic Survey (JADES) will likely characterize ∼15–30 galaxies similar to 3MM-1 (based on expected number densities from this work and Zavala et al. 2018a, and observations described in Williams et al. 2018). JWST will have the capability to measure stellar population properties and redshifts, and in combination with ALMA far-infrared constraints, will enable a detailed investigation into the star formation, dust, and stellar population properties of massive galaxies in the early universe.

We thank the anonymous referee whose valuable suggestions have improved the paper. We are grateful to James Simpson for making part of S2COSMOS available for our analysis in advance of publication. We thank Jorge Zavala, Peter Behroozi, and Pieter van Dokkum for helpful discussions. C.C.W. acknowledges support from the National Science Foundation Astronomy and Astrophysics Fellowship grant AST-1701546. J.L. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1701487. This paper makes use of the following ALMA data: ADS/JAO.ALMA #2018.1.01739.S, ADS/JAO.ALMA #2015.1.00853.S, ADS/JAO.ALMA #2015.1.00861.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), 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.

Appendix

Here we present supplemental observational details for 3MM-1 and 3MM-2. In Table 2, we present the deblended photometry (measurements and 1σ uncertainties) based on the ALMA priors for the two 3-mm sources that are used in the SED-fitting analysis presented in Section 2.4. In Figures 7 and 8 we present the posterior probability distributions for the galaxy properties of both sources measured using the Bayesian SED fitting with bagpipes.

Figure 7.

Figure 7. Posterior probability distributions for the fitted parameters of 3MM-1 corresponding to the orange region in Figure 4, assuming Charlot & Fall (2000) dust attenuation with the slope left free. Mdel refers to the total stars formed (integral of the delayed star formation history), whereas M* is the stellar mass excluding remnants and including mass loss due to stellar evolution, and is the parameter referred to throughout the text. Parameters, their definitions, and their uncertainties are presented in Table 1.

Standard image High-resolution image
Figure 8.

Figure 8. Posterior probability distributions for the fitted parameters of 3MM-2 corresponding to the orange region in Figure 4, assuming Charlot & Fall (2000) dust attenuation with the slope left free, and fixed to the spectroscopic redshift. Measured parameters, their definitions, and their uncertainties are presented in Table 1.

Standard image High-resolution image

Table 2.  ALMA Prior-based Deblended Photometry

Band 3MM-1 flux density rms 3MM-2 flux density rms
  (μJy)   (μJy)  
SUBARU B [−1.3E-02] 7.8E-03
HSC g [3.6E-03] 1.5E-02 [2.3E-02] 1.1E-02
SUBARU V [−7.6E-03] 1.8E-02
HSC r [−1.2E-02] 1.4E-02 9.3E-02 9.7E-03
SUBARU rp [1.2E-03] 1.7E-02
SUBARU ip [2.0E-03] 2.6E-02
HSC i [−1.0E-02] 2.1E-02 1.2E-01 1.5E-02
HSC z [1.2E-02] 3.0E-02 1.3E-01 2.4E-02
SUBARU zp [−1.5E-01] 6.5E-02
HSC Y [2.2E-02] 7.4E-02 [4.9E-02] 5.8E-02
UltraVISTA Y [9.7E-03] 8.9E-02 [1.6E-01] 8.2E-02
UltraVISTA J [−1.1E-01] 1.0E-01 [2.0E-01] 8.9E-02
UltraVISTA H [−5.9E-02] 1.3E-01 6.3E-01 1.1E-01
UltraVISTA Ks [1.2E-01] 8.8E-02 8.4E-01 8.7E-02
Spitzer/IRAC 3.6 μm [3.2E-01] 1.2E-01 1.7E+00 9.2E-02
Spitzer/IRAC 4.5 μm [3.2E-01] 1.3E-01 3.0E+00 9.9E-02
Spitzer/IRAC 5.8 μm [5.2E+00] 3.3E+00 8.8E+00 2.2E+00
Spitzer/IRAC 8 μm [6.0E+00] 4.1E+00 [7.7E+00] 2.9E+00
Spiter/MIPS 24 μm [1.1E+01] 1.5E+01 [2.6E+01] 1.5E+01
Herschel/PACS 100 μm [7.1E+02] 1.4E+03 [1.5E+02] 1.4E+03
Herschel/PACS 160 μm [−2.4E+03] 3.2E+03 [−1.9E+03] 3.2E+03
Herschel/SPIRE 250 μm [3.7E+03] 3.9E+03 [1.0E+03] 3.9E+03
Herschel/SPIRE 350 μm [4.3E+03] 4.9E+03 [−5.0E+02] 4.9E+03
Herschel/SPIRE 500 μm [5.5E+03] 5.8E+03 [4.4E+03] 5.8E+03
SCUBA2 850 μm 3.5E+03 1.1E+03 [4.0E+02] 1.1E+03
ALMA 3 mm 1.6E+02 2.2E+01 7.5E+01 1.0E+01
VLA 3 GHz 1.0E+01 2.4E+00 [3.9E+00] 2.4E+00
VLA 1.4 GHz [1.7E+01] 1.7E+01 [4.0E+01] 1.7E+01

Note. Subaru optical photometry for 3MM-2 are not measured owing to cosmetic defects in the mosaics that cross the location of the source. Photometric upper limits as indicated by downward arrows in Figure 4 are at the rms values of photometric points in this table where S/N < 1. Non-significant photometric points (with SNR < 3) are indicated with brackets.

Download table as:  ASCIITypeset image

Please wait… references are loading.
10.3847/1538-4357/ab44aa