The ALMA Spectroscopic Survey Large Program: The Infrared Excess of z = 1.5–10 UV-selected Galaxies and the Implied High-redshift Star Formation History

The principal data used are the band 6 ALMA observations from the 2016.1.00324.L program over the HUDF. Those observations were obtained through a full frequency scan in band 6 (212–272 GHz) with ALMA in its most compact configuration. The observations are distributed over 85 pointings separated by 11'' and cover an approximate area of ∼4.2 arcmin² to near uniform depth. Our construction of a continuum mosaic from ALMA data is described in González-López et al. (2020). The peak sensitivity in our 1.2 mm continuum observations is 9.3 μJy (1σ) per synthesized beam (153 × 108; González-López et al. 2020).

For HST optical Advanced Camera for Surveys (ACS)/Wide Field Camera (WFC) and near-infrared WFC3/IR observations, we make use of the eXtreme Deep Field (XDF) reductions (Illingworth et al. 2013), which incorporated all ACS+WFC3/IR data available over the HUDF in 2013. The XDF reductions are ∼0.1–0.2 mag deeper than original Beckwith et al. (2006) reductions at optical wavelengths and also provide coverage in the F814W band. The WFC3/IR reductions made available as part of the XDF release include all data from the original HUDF09 (Bouwens et al. 2011), the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011), and the HUDF12 (Ellis et al. 2013) programs. Subsequent to the XDF release, only 17 additional orbits of HST imaging data have been obtained with HST over the XDF region (five of which are in the F105W band and 12 in the F435W band). Given that this is <4% of the integration time already included in the XDF release, we elected to use the XDF release due to the effort putting into using super sky flats to optimize the sensitivity. ²²

For the 0.2–0.4 μm WFC3/UVIS data over the ASPECS field, we made use of the v2 release of the UVUDF epoch 3 data (Teplitz et al. 2013; Rafelski et al. 2015) which included imaging data in the F225W, F275W, and F336W bands. The Spitzer/Infrared Array Camera (IRAC) observations we utilize are from the ∼200 hr stacks of the IRAC observations over the HUDF from the GOODS Reionization Era wide-Area Treasury from Spitzer (GREATS) program (M. Stefanon et al. 2020, in preparation: PI: Labbé).

Photometry for sources in our samples is performed in the same way as in the Bouwens et al. (2016) analysis from the ASPECS pilot program. HST fluxes are derived using our own modified version of the SExtractor (Bertin & Arnouts 1996) software. Source detection is performed on the square root of ${\chi }^{2}$ image (Szalay et al. 1999; similar to a coadded image) constructed from the V₆₀₆, i₇₇₅, Y₁₀₅, J₁₂₅, JH₁₄₀, and H₁₆₀ images. After point-spread function (PSF)-correcting fluxes to match the H₁₆₀-band image, color measurements are made in Kron-style (1980) scalable apertures with a Kron factor of 1.6. "Total magnitude" fluxes are derived by (1) correcting up the fluxes in smaller scalable apertures to account for the additional flux seen in a larger scalable aperture (Kron factor of 2.5) seen on the square root of ${\chi }^{2}$ image and (2) correcting for the flux outside these larger scalable apertures and on the wings of the PSF using tabulations of the encircled energy, appropriate for point sources (Dressel 2012).

As in our earlier analysis and many other analyses (e.g., Shapley et al. 2005; Grazian et al. 2006; Labbé et al. 2006, 2010, 2015; Laidler et al. 2007; Merlin et al. 2015), Spitzer/IRAC photometry was performed using the HST observations as a template to model the fluxes of sources in the Spitzer/IRAC observations and thus perform photometry below the nominal confusion limit. In performing photometry, a simultaneous fit of the flux of a source of interest and its neighbors is performed, the flux from neighboring sources is subtracted, and then aperture photometry on the source of interest is performed. Photometry is performed in 18 diameter circular apertures for the Spitzer/IRAC 3.6 μm and 4.5 μm bands and 20 diameter circular apertures for the 5.8 μm and 8.0 μm bands. The observed fluxes are corrected to total based on the inferred growth curve for sources after PSF correction to the Spitzer/IRAC PSF.

A similar procedure is used to derive fluxes for sources from the deep ground-based K-band observations available from the Very Large Telescope (VLT)/Hawk-I UDS and GOODS Survey (Fontana et al. 2014), VLT/ISAAC, and Panoramic Near Infrared Camera (PANIC) observations over the HUDF (5σ depths of 26.5 mag).

The purpose of this subsection is to summarize our approach in modeling the far-IR SED of faint, UV-selected z = 1.5–10 galaxies. Having accurate constraints on the overall form of the far-IR SED for these galaxies is potentially important for interpreting far-IR continuum observations of the distant universe to quantify the dust-obscured SFRs. The goal of this subsection will be to use a variety of published observations from the literature to motivate the approach we will utilize throughout the rest of this manuscript.

As is common practice (e.g., Casey 2012), we will adopt a modified blackbody (MBB) form to model the far-IR spectral energy distributions of galaxies (e.g., Casey 2012), with a dust emissivity power-law spectral index of ${\beta }_{d}=1.6$, which is toward the center of the range of values (1.5–2.0) frequently found in the observations (Eales et al. 1989; Klaas et al. 1997). MBB SEDs have the advantage of being relatively simple in form, but are known to show less flux at mid-IR wavelengths than galaxies with a prominent mid-IR power-law component. Fortunately, the impact of such differences on the conversion factors from the 1.2mm flux densities we observe and the total IR luminosity is relatively modest (i.e., factors of ≲1.5; see, e.g., Casey et al. 2018), especially relative to other issues like the dust temperature.

Characterizing the evolution of the dust temperature as a function of redshift is challenging due to both selection bias and the significant dependence the dust temperature can show on other quantities like the bolometric luminosity, specific star formation, and the wavelength where dust becomes opaque (e.g., Magnelli et al. 2014; Liang et al. 2019; Ma et al. 2019) which are arguably larger and more significant than the impact of redshift on the dust temperature.

Nevertheless, there have been multiple studies looking at the evolution of dust temperature in galaxies with redshift for fixed values of the bolometric luminosity (e.g., Béthermin et al. 2015; Schreiber et al. 2018). One particularly comprehensive recent study on this front has been by Schreiber et al. (2018), who consider the apparent evolution in dust temperatures from z ∼ 4 to z ∼ 0 using stacks of the available Herschel observations.

In Figure 1, we present the same observations that Schreiber et al. (2018) consider, and then add earlier results from Béthermin et al. (2015) to their constraints. Finally, we also include the dust temperature measurements obtained by Pavesi et al. (2016) on a $z\sim 5.25$ galaxy, by Knudsen et al. (2017) on a $z\sim 7.5$ galaxy, by Hashimoto et al. (2019) on a z = 7.15 galaxy, by Harikane et al. (2020) on two $z\sim 6.1$ galaxies, by Bakx et al. (2020) on the Tamura et al. (2019) z = 8.31 galaxy, by Faisst et al. (2020) on four $z\sim 5.5$ galaxies, and by Béthermin et al. (2020) on stacks of z = 4–5 and z = 5–6 galaxies, as well as the median dust temperatures measured by Strandet et al. (2016) on their sample of bright South Pole Telescope (SPT) sources. Each of these temperature measurements is reported to be corrected for the impact of cosmic microwave background (CMB) radiation (da Cunha et al. 2013).

Zoom In Zoom Out Reset image size

To make the present dust temperature measurements in Figure 1 as consistent as possible, all measurements have been converted to their equivalent values using an emissivity index ${\beta }_{d}$ of 1.6 and using the light-weighted dust temperatures (converting the Schreiber et al. 2018 temperatures from the mass-weighted temperatures to light-weighted temperatures using their Equation (6)). Pursuing a joint fit to all dust temperature measurements in Figure 1, we derive the following relationship between dust temperature and redshift:

$$\begin{eqnarray}&&{T}_{d}[K]=(34.6\pm 0.3)+(3.94\pm 0.26)(z-2).\end{eqnarray} \tag{ 1 }$$

The best-fit evolution we derive for the dust temperature is higher than what Schreiber et al. (2018) derive (${T}_{d}[K]\,=(32.9\pm 2.4)+(4.60\pm 0.35)(z-2)$) due to our use of light-weighted dust temperatures where the dust temperatures are higher. Our best-fit relation for the temperature evolution does, however, evolve slightly less steeply with redshift, largely as a result of our inclusion of constraints from SPT sources, the four Faisst et al. (2020) $z\sim 5.5$ galaxies, and the new Béthermin et al. (2020) stack constraints for z = 4–6 galaxies. This best-fit evolution is also not especially dissimilar from the trends found in theoretical models such as those by Narayanan et al. (2018), Liang et al. (2019), and Ma et al. (2019). In the Narayanan et al. (2018) results, the dust temperature increases from 40–50 K in galaxies at z ∼ 2–3 galaxies to 55–70 K at z ∼ 6–7. In Liang et al. (2019) and Ma et al. (2019), the evolution in dust temperature expected on the basis of the evolution of the MASSIVEFIRE sample is ${(1+z)}^{0.36\pm 0.06}$ (their Table 2), similar to that implied by Equation (1) above.

Despite the clear evolution in temperature found here and earlier by Béthermin et al. (2015) and Schreiber et al. (2018), other recent studies find no less evolution in dust temperature with redshift. For example, Ivison et al. (2016) infer only ∼50% as much evolution in the dust temperature as we find, while other studies, e.g., Dudzevičiūtė et al. (2020), find no significant evolution in the dust temperature of galaxies with redshift when a purely luminosity-limited sample is studied (see also Strandet et al. 2016). Dudzevičiūtė et al. (2020) have argued that the apparent temperature evolution that studies such as Schreiber et al. (2018) have found is likely a consequence of luminosity variations in that study. Given this, we also consider there being less evolution of the dust temperature of galaxies with cosmic time than in our fiducial models.

Assuming that the effective dust temperature of obscured star formation in z ∼ 1.5–10 galaxies follows the same evolution as given by Equation (1), we can derive the limiting dust-obscured SFR we would be able to detect as a function of redshift from our program. Adopting a modified blackbody form for the SED shape described at the beginning of this subsection and accounting for the impact of the CMB (e.g., Section 3.1.1 of da Cunha et al. 2013), we estimate that we should be able to tentatively detect at 4σ any star-forming galaxy at z > 2 with an IR luminosity (8–1000 μm rest frame) in excess of 6.8 × 10¹⁰ L_⊙ at z ∼ 2, 9.0 × 10¹⁰ L_⊙ at z ∼ 3, and ∼11.2–28.4 × 10¹⁰ L_⊙ at z ∼ 4–10. We verified that use of potentially more realistic far-IR SED templates than a modified blackbody form, following, e.g., Álvarez-Márquez et al. (2016) with a mid-IR power law, yields similar 1.2 mm to IR luminosity conversion factors (see also Appendix A of Fudamoto et al. 2020a).

Adopting the Kennicutt (1998) conversion between IR luminosity and the SFR, these limits translate to 4σ limits on the obscured SFRs of 6.8 M_⊙ yr⁻¹, 9.0 M_⊙ yr⁻¹, and 11.2–28.4 M_⊙ yr⁻¹, respectively, at these redshifts. If we instead allow for much less evolution in the dust temperature, such that the typical dust temperature at z ∼ 4–8 is 35 K, the 4σ limits from ASPECS translates to limits on the obscured SFRs of 4–5 M_⊙ yr⁻¹.

In Table 1, we provide these limiting luminosities and SFRs in tabular form, while providing for context these limits for modified blackbody SEDs if the dust temperature is fixed at 35 K or 50 K.

Table 1. 4σ Sensitivity Limits for Our Probe of Obscured Star Formation from Individual z ≳ 1.5 Galaxies and the Dependence on SED

Far-IR	4σ Sensitivity Limits (1010 L⊙)
SED Model	z∼2	z∼3	z∼4	z∼5	z∼6	z∼7	z∼8	z∼9	z∼10
Fiducial evolving a , b	6.8	9.0	11.2	13.6	16.1	18.7	21.7	24.9	28.4
35K graybody b	7.1	6.3	5.5	5.1	4.8	4.7	4.7	4.9	5.2
50K graybody b	30.8	25.2	20.9	17.8	15.7	14.4	13.6	13.3	13.4
	4σ Limit for Probes of the Obscured SFR (M⊙ yr−1) c
SED Model	z∼2	z∼3	z∼4	z∼5	z∼6	z∼7	z∼8	z∼9	z∼10
Fiducial evolving a , b	6.8	9.0	11.2	13.6	16.1	18.7	21.7	24.9	28.4
35K graybody b	7.1	6.3	5.5	5.1	4.8	4.7	4.7	4.9	5.2
50K graybody b	30.8	25.2	20.9	17.8	15.7	14.4	13.6	13.3	13.4
	Dust Temperatures for Fiducial Evolving SED Model (K)
	34.6	38.5	42.5	46.4	50.4	54.3	58.2	62.2	66.1

Notes.

^a Using Equation (1). ^b Standard modified blackbody form (e.g., Casey 2012) with a dust emissivity power-law spectral index of ${\beta }_{d}=1.6$ (Eales et al. 1989; Klaas et al. 1997). ^c The Kennicutt (1998) conversion factor from IR luminosity to SFR is adopted.

Download table as:  ASCIITypeset image

In constructing samples of z = 1.5–10 galaxies for examination with the ASPECS ALMA data, we utilize both Lyman-break selection criteria as well as photometric redshift selection to ensure our samples are as comprehensive as possible.

To ensure consistency with earlier results from our pilot study (Bouwens et al. 2016), we have adopted essentially identical color–color and photometric redshift selection criteria to those applied in Bouwens et al. (2016). z = 1.5–3.5 sources are identified using the same Lyman-break color criteria we had earlier used in Bouwens et al. (2016) and identified by running the EAZY photometric redshift code (Brammer et al. 2008) on our own HST WFC3/UVIS, ACS, and WFC3/IR photometric catalogs. Our z ∼ 2 and z ∼ 3 color criteria are as follows:

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 2:({{\rm{UV}}}_{275}-{U}_{336}\gt 1)\wedge ({U}_{336}-{B}_{435}\lt 1)\wedge \\ & & ({V}_{606}-{Y}_{105}\lt 0.7)\wedge ({\rm{S}}/{\rm{N}}({{\rm{UV}}}_{225})\lt 1.5),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 3:({U}_{336}-{B}_{435}\gt 1)\wedge ({B}_{435}-{V}_{606}\lt 1.2)\wedge \\ & & ({i}_{775}-{Y}_{105}\lt 0.7)\wedge ({\chi }_{{{\rm{UV}}}_{225},{{\rm{UV}}}_{275}}^{2}\lt 2)\end{array}\end{eqnarray*}$$

where ∧, ∨, and S/N represent the logical AND, OR symbols, ${\chi }^{2}$ is the ${\chi }^{2}$ parameter defined in Bouwens et al. (2011), and signal-to-noise ratio in our smaller scalable apertures, respectively. We also made use of the photometric catalog of Rafelski et al. (2015) and included those sources in our samples, if not present in the other selections.

Our z = 4–8 samples are drawn from the Bouwens et al. (2015) samples and include all z = 3.5–8.5 galaxies located over the 4.2 arcmin² ASPECS region. The Bouwens et al. (2015) samples were based on the deep optical ACS and WFC3/IR observations within the HUDF. z = 4–8 samples were constructed by applying Lyman-break-like color criteria to the XDF reduction (Illingworth et al. 2013) of the HUDF. Those criteria are the following for our z ∼ 4, 5, 6, 7, and 8 selections:

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 4:({B}_{435}-{V}_{606}\gt 1)\wedge ({i}_{775}-{J}_{125}\lt 1)\wedge \\ & & ({B}_{435}-{V}_{606}\gt 1.6({i}_{775}-{J}_{125})+1),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 5:({V}_{606}-{i}_{775}\gt 1.2)\wedge ({z}_{850}-{H}_{160}\lt 1.3)\wedge \\ & & ({V}_{606}-{i}_{775}\gt 0.8({z}_{850}-{H}_{160})+1.2),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 6:({i}_{775}-{z}_{850}\gt 1.0)\wedge ({Y}_{105}-{H}_{160}\lt 1.0)\wedge \\ & & ({i}_{775}-{z}_{850}\gt 0.777({Y}_{105}-{H}_{160})+1.0),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 7:({z}_{850}-{Y}_{105}\gt 0.7)\wedge ({J}_{125}-{H}_{160}\lt 0.45)\wedge \\ & & ({z}_{850}-{Y}_{105}\gt 0.8({J}_{125}-{H}_{160})+0.7),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 8:({Y}_{105}-{J}_{125}\gt 0.45)\wedge ({J}_{125}-{H}_{160}\lt 0.5)\wedge \\ & & ({Y}_{105}-{J}_{125}\gt 0.75({J}_{125}-{H}_{160})+0.525).\end{array}\end{eqnarray*}$$

The six galaxies in our z = 9–10 samples are identified by applying the following Y₁₀₅/J₁₂₅-dropout Lyman-break color criteria to the available HST data:

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 9:(({Y}_{105}-{H}_{160})+2({J}_{125}-{{JH}}_{140})\gt 1.5)\wedge \\ & & (({Y}_{105}-{H}_{160})+2({J}_{125}-{{JH}}_{140})\gt \\ & & \,\,\,\,\,\,\,1.5+1.4({{JH}}_{140}-{H}_{160}))\wedge \\ & & ({{JH}}_{140}-{H}_{160}\lt 0.5)\wedge ({J}_{125}-{H}_{160}\lt 1.2),\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}z & \sim & 10:({J}_{125}-{H}_{160}\gt 1.2)\wedge (({H}_{160}-[3.6]\lt 1.4)\vee \\ & & ({\rm{S}}/{\rm{N}}([3.6])\lt 2)).\end{array}\end{eqnarray*}$$

Selected sources are required to be undetected (<2σ) in all HST passbands blueward of the break, both individually and in a stack. Potential stars are excluded from our selection using the measured SExtractor (Bertin & Arnouts 1996) stellarity criterion.

To ensure that our color criteria were not included, for galaxies which are likely passive we adopt a UVJ-like criterion (Williams et al. 2009) which allow us to exclude passive galaxies from our z ≳ 1.5 selection of star-forming galaxies. Specifically, we adopt the prescription given in Pannella et al. (2015):

$$\begin{eqnarray*}&&\begin{array}{l}(U-V\lt 1.3)\wedge (V-J\gt 1.6)\wedge \\ \quad (U-V\lt 0.88(V-J)+0.59),\end{array}\end{eqnarray*}$$

which is very similar to the prescription given in Williams et al. (2009). Application of this criteria to our z ∼ 1.5–10 selection results in the exclusion of just one source from our selection.

The z ∼ 2, 3, 4, 5, 6, 7, 8, 9, and 10 selections we consider over the ASPECS footprint include 447, 203, 395, 139, 94, 54, 24, 4, and 2 distant sources, respectively (Table 2). The expected contamination levels in these color-selected samples by lower-redshift galaxies (or stars) is estimated to be on the order of 3%–8% (e.g., Bouwens et al. 2015). Sources in our selection have apparent magnitude in the UV-continuum extending from 23.5–30.5 mag (Figure 2 (left panel)).

Zoom In Zoom Out Reset image size

Table 2. Number of UV-selected z ∼ 2, z ∼ 3, z ∼ 4, z ∼ 5, z ∼ 6, z ∼ 7, z ∼ 8, z ∼ 9, and z ∼ 10 Galaxies Located within Our 4.2 arcmin² ASPECS Footprint

		Number of
Redshift	Selection Criterion	Sources	References
z ∼ 2	UV275-dropout or
	$1.5\lt {z}_{\mathrm{phot}}\lt 2.5$	447	R15/This work
z ∼ 3	U336-dropout or
	$2.5\lt {z}_{\mathrm{phot}}\lt 3.5$	203	R15/This work
z ∼ 4	B435-dropout or
	$3.5\lt {z}_{\mathrm{phot}}\lt 4.5$	395	B15/This work
z ∼ 5	V606-dropout	139	B15
z ∼ 6	i775-dropout	94	B15
z ∼ 7	z850-dropout or
	$6.5\lt {z}_{\mathrm{phot}}\lt 7.5$	54	B15/This work
z ∼ 8	Y105-dropout	24	B15
z ∼ 9	Y105-dropout	4	This work
${zi}\sim 10$	J125-dropout	2	This work
	Total	1362

References. B15 = Bouwens et al. (2015), R15 = Rafelski et al. (2015).

Download table as:  ASCIITypeset image

Based on an abundance of previous work, it is well known that the IRX is correlated with the measured UV-continuum slope of galaxies (e.g., Meurer et al. 1999) and also the stellar mass (e.g., Whitaker et al. 2017).

For each of the sources over ASPECS, we derive a UV-continuum slope β fitting the HST photometry in various bands probing the UV-continuum to a power law ${f}_{1600}{(\lambda /1600\mathring{\rm A} )}^{\beta }$ to derive a mean flux at ∼1600 Å and also a spectral slope β. Flux measurements in band passes that could be impacted by IGM absorption or rest-frame optical ≳3500 Å light are excluded. The inclusion of photometric constraints on the UV-continuum even to ∼3000 Å is expected to have little impact on the derived β given the general power-law-like shape of the UV continuum (e.g., see Appendix A in Wilkins et al. 2016). Due to the limited wavelength leverage available to derive UV-continuum sources for sources at z = 8–10, we take the UV-continuum slope β to be uniformly −2.2 consistent with the results of Bouwens et al. (2014a).

As in other work (e.g., Sawicki & Yee 1998; Brinchmann & Ellis 2000; Papovich et al. 2001; Labbé et al. 2006; González et al. 2014), we estimate stellar masses for individual sources in our samples by modeling the observed photometry using stellar population libraries and considering variable (or fixed) star formation histories, metallicities, and dust content.

For z ∼ 1.5–10 sources in our catalogs, we make use of the publicly-available code FAST (Kriek et al. 2009) to perform this fitting. We assume a Chabrier (2003) IMF, a metallicity of 0.2 Z_⊙, a stellar population age from 10 Myr to the age of the universe, and allow the dust extinction in the rest-frame V to range from 0 to 2 mag, which we acknowledge may be inadequate for some especially dust-rich galaxies (e.g., Simpson et al. 2017). We assume an ${e}^{-t/\tau }$ star formation history and allow the τ parameter to have any value from 1–100 Gyr. Our fixing the fiducial metallicity to 0.2 Z_⊙ is motivated by studies of the metallicity of individual z ∼ 2–4 galaxies (Pettini et al. 2000) or as predicted from cosmological hydrodynamical simulations (Finlator et al. 2011; Wise et al. 2012). While the current choice of parameters can have a sizeable impact on inferred quantities like the age of a stellar population (changing by >0.3–0.5 dex), these choices typically do not have a major impact ($\gtrsim 0.2$ dex) on the inferred stellar masses.

In deriving the stellar masses for individual sources, use is made of flux measurements from 11 HST bands (UV₂₂₅, UV₂₇₅, U₃₃₆, B₄₃₅, V₆₀₆, i₇₇₅, z₈₅₀, Y₁₀₅, J₁₂₅, JH₁₄₀, H₁₆₀), one band in the near-IR from the ground (K_s ), and four Spitzer/IRAC bands (3.6 μm, 4.5 μm, 5.8 μm, and 8.0 μm). The HST photometry we use for estimating stellar masses is derived applying the same procedure as used for selecting our z ∼ 1.5–3.5 Lyman-break galaxies (LBG) samples (see Section 2.2).

A modest correction is made to the Spitzer/IRAC 3.6 μm and 4.5 μm photometry to account for the impact of nebular emission lines on the observed IRAC fluxes. Specifically, the 3.6 μm and 4.5 μm band fluxes of galaxies in the redshift ranges z = 3.8–5.0 and z = 5.1–6.6, respectively, are reduced by 0.32 mag and 0.35 mag, respectively, to remove the contribution of the Hα+[N ii] emission lines to the broadband fluxes. A 0.32 mag and 0.35 mag correction is appropriate for a rest-frame equivalent width (EW) of ∼500 Å and ∼540 Å , respectively, for the Hα+[N ii] emission lines, consistent with most determinations of the Hα+[N ii] emission line EW over the range z = 3.8–5.4 (Stark et al. 2013; Mármol-Queraltó et al. 2016; Faisst et al. 2016; Smit et al. 2016; Rasappu et al. 2016). For galaxies in the redshift ranges z = 5.4–7.0 and z = 7.0–9.1, the measured fluxes in the 3.6 μm and 4.5 μm bands are reduced by 0.5 mag. A 0.5 mag correction is appropriate for a rest-frame EW of ∼680 Å for the Hα+[N ii] emission lines, consistent with most determinations of the Hα+[N ii] emission line EW over the range z = 3.8–5.4 (Labbé et al. 2013; Smit et al. 2014, 2015; Faisst et al. 2016; Endsley et al. 2020). The fiducial stellar mass estimates we derive using FAST are typically ∼0.1 dex lower than using other stellar population codes like MAGPHYS and Prospector (see Appendix A).

The middle panel of Figure 2 illustrates the effective range in stellar mass probed by our z = 1.5–10 sample. Most sources from our HUDF z = 1.5–10 sample have stellar masses in the range ${10}^{7.5}$–${10}^{9.5}$ M_⊙. The most massive sources probed by our program extend to ${10}^{11.5}$ M_⊙. Beyond the stellar mass itself, Figure 2 also illustrates the range in UV-continuum slope β probed by our samples (see Section 3.1 for details on how β is derived). Since the measured β has been demonstrated to be quite effective in estimating the IRX for lower-redshift UV-selected samples (e.g., M99; Reddy et al. 2006; Daddi et al. 2007), it is useful for us to probe a broad range in β. As can be seen from Figure 2, our samples probe the range $\beta \sim -1.5$ to $\sim -2.5$ quite effectively.

Thanks to the limited evolution seen in the IRX versus stellar mass and IRX versus β results over the entire redshift range z ∼ 3 to z ∼ 0 (Reddy et al. 2006; Whitaker et al. 2017; Fudamoto et al. 2020a), we might expect these relations to be at least approximately valid to even higher redshifts.

Before looking in detail at which sources show continuum detections and what their properties are, let us briefly calculate how many sources we would expect to detect based on published IRX versus stellar mass and IRX versus UV-continuum slope β relations. Given the limited evolution in these relations, we expect the predicted results to be reasonably accurate in estimating the overall numbers from our program. For our baseline IRX–stellar mass M relation, we take the relation derived in our pilot program (Appendix A from Bouwens et al. 2016):

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\mathrm{IRX}}_{M,0}={\mathrm{log}}_{10}M-9.17.\end{eqnarray} \tag{ 3 }$$

For our baseline IRX–β relation, we make use of the consensus low-redshift relation derived in Appendix B based on the following three studies (Overzier et al. 2011; Takeuchi et al. 2012; Casey et al. 2014). The relation we derive is the following:

$$\begin{eqnarray}&&{{\rm{IRX}}}_{z=0}=1.7({10}^{0.4(1.86(\beta +1.85))}-1).\end{eqnarray} \tag{ 4 }$$

The IRX implied by the above relation are ≈0.5 times that of the Meurer et al. (1999) relation. Equivalent expressions for a Reddy (similar to Calzetti et al. 2000) and SMC-like dust law are the following:

$$\begin{eqnarray}&&{\mathrm{IRX}}_{\mathrm{Reddy}}=1.7({10}^{0.4(1.84(\beta +1.85))}-1)\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{\mathrm{IRX}}_{\mathrm{SMC}}=1.7({10}^{0.4(1.1(\beta +1.85))}-1).\end{eqnarray} \tag{ 6 }$$

Based on the above relations and observed UV fluxes, we can compute the equivalent flux at an observed wavelength of 1.26 mm adopting a modified blackbody form with a dust emissivity power-law spectral index of ${\beta }_{d}=1.6$ and dust temperature given by Equation (1). To account for the impact of the CMB at z ∼ 1.5–10 on the expected flux densities we would measure, we multiply the predicted flux (before consideration of CMB effects) by ${C}_{\nu }$

$$\begin{eqnarray}&&{C}_{\nu }=\left[1-\displaystyle \frac{{B}_{\nu }({T}_{\mathrm{CMB}}(z))}{{B}_{\nu }({T}_{d}(z))}\right],\end{eqnarray} \tag{ 7 }$$

following prescriptions given in da Cunha et al. (2013).

Using the above procedure, we calculated the expected flux for our entire sample of 1362 z = 1.5–10 galaxies identified over the 4.2 arcmin² ASPECS footprint alternatively making use of the consensus IRX–stellar mass relation from Bouwens et al. (2016), our consensus low-redshift IRX–β relation, and also a SMC-like IRX–β relation (Equations (3)–(6)). Fifteen, 28, and eight sources, respectively, are predicted to show >4σ detections in the ASPECS observations in the 1.2 mm continuum. Assuming a fixed dust temperature of 35 K, the predicted numbers would be 27, 42, and 11, respectively. Figure 3 shows the predicted IR luminosities versus redshift using either the aforementioned IRX–stellar mass relation (left) or the IRX–β relation (right) for our fiducial dust temperature model. The solid red and dotted lines show the 4σ IR luminosity limit we probe with the ASPECS data set adopting the fiducial dust temperature model given in Equation (1) (solid red line) and assuming the dust temperature remains fixed at 35 K for all of cosmic time (dotted red line).

Zoom In Zoom Out Reset image size

Examination of the 1362 z = 1.5–10 galaxies over our sensitive ASPECS mosaic shows that 18 of these galaxies are detected at >4.0σ in the 1.2 mm continuum images. We use the flux densities and uncertainties that González-López et al. (2020) derive for each source from the 1.2 mm continuum images. González-López et al. (2020) make use of flux density measurements made from the tapered images, allowing for a more complete account of the total dust-continuum flux density in sources, many of which are spatially extended. The coordinates and source properties of the continuum detected sources are provided in Table 3. 1.2 mm continuum images of the 4σ-detected sources are presented in Figure 4 and shown with respect to the HST and Spitzer/IRAC images.

Zoom In Zoom Out Reset image size

Table 3.  z ≳ 1.5 UV-selected Galaxies Showing 4σ Detections in Our Deep ALMA Continuum Observations

					log10		Measured	Inferred
			${m}_{\mathrm{UV},0}$		$M/$ M⊙		${f}_{1.2\mathrm{mm}}$	LIR
ID e	R.A.	Decl.	(mag)	z		β	(μJy) f	(1010 L⊙)	References a
XDFU-2435246390 (C06)	03:32:43.52	−27:46:39.0	27.6	2.696 c	10.92	−0.3 ± 0.4	1071 ± 46	259 ± 11	3
XDFU-2385446340 (C01)	03:32:38.54	−27:46:34.0	24.4	2.543 c	9.90	−1.2 ± 0.1	752 ± 10	226 ± 3	1, 2, 3
XDFU-2397246112 (C05)	03:32:39.72	−27:46:11.2	24.9	1.551 c	11.10	−0.4 ± 0.1	461 ± 14	112 ± 3	1, 2, 3
XDFU-2369747272 (C02)	03:32:36.97	−27:47:27.2	26.9	1.76 b	10.66	1.3 ± 0.2	432 ± 9	104 ± 2	3
XDFU-2400547554 (C10)	03:32:40.05	−27:47:55.4	23.6	1.997 c	10.83	−0.4 ± 0.1	342 ± 18	83 ± 4	3
XDFU-2410746315 (C04)	03:32:41.07	−27:46:31.5	27.0	2.454 c	9.39	−0.8 ± 0.1	316 ± 11	95 ± 3	3
XDFU-2433446471 (C11)	03:32:43.34	−27:46:47.1	28.2	2.76 b	11.00	0.5 ± 0.2	289 ± 21	87 ± 6	3
XDFU-2350746475 (C07)	03:32:35.07	−27:46:47.5	26.6	2.58 c	10.89	0.5 ± 0.2	233 ± 11	56 ± 3	3
XDFU-2416846554 (C14a)	03:32:41.68	−27:46:55.4	27.4	1.999 c	10.47	0.6 ± 0.3	185 ± 10	45 ± 2
XDFB-2380246263 (C08)	03:32:38.02	−27:46:26.3	25.4	3.711 d	10.81	2.9 ± 0.1	163 ± 10	59 ± 4	1
XDFB-2355547038 (C09)	03:32:35.55	−27:47:03.8	26.2	3.601 c	9.47	−0.8 ± 0.1	155 ± 9	56 ± 3
XDFU-2387248103 (C24)	03:32:38.72	−27:48:10.3	26.0	2.68 b	9.45	−0.5 ± 0.1	134 ± 24	40 ± 7
XDFU-2373546453 (C18)	03:32:37.35	−27:46:45.3	23.9	1.845 d	10.49	−0.7 ± 0.1	107 ± 10	26 ± 2	1, 2
XDFU4596 (C17)	03:32:38.80	−27:47:14.8	24.5	1.848 d	10.46	−0.6 ± 0.1	97 ± 9	23 ± 2
XDFU-2361746276 (C19)	03:32:36.17	−27:46:27.6	25.4	2.574 c	10.59	−0.2 ± 0.1	85 ± 12	20 ± 3	1
XDFU9838 (C26)	03:32:34.68	−27:46:44.5	25.5	1.552 d	10.31	−0.2 ± 0.1	65 ± 15	16 ± 4
XDFU-2359847256 (C21)	03:32:35.98	−27:47:25.6	25.2	2.69 b	10.24	−1.0 ± 0.1	58 ± 10	18 ± 3
XDFU-2370746171 g (C31)	03:32:37.07	−27:46:17.1	23.7	2.227 d	9.49	−1.3 ± 0.1	47 ± 11	14 ± 3	2

Notes.

^a References previously reporting continuum detections of the identified sources: [1] Aravena et al. (2016), [2] Bouwens et al. (2016), [3] Dunlop et al. (2017). ^b Photometric redshift. ^c Spectroscopic redshift from the detection of a CO line in the ASPECS ALMA data (Boogaard et al. 2019a, 2019b). ^d Spectroscopic redshift available for this source from the MUSE guaranteed time observations over the HUDF (Bacon et al. 2017). ^e The source IDs included inside the parentheses are as in González-López et al. (2020) and Aravena et al. (2020). ^f Measurements as in González-López et al. (2020). ^g This source was previously reported as a tentative 2.3σ detection in Bouwens et al. (2016).

Download table as:  ASCIITypeset image

The IR luminosities we estimated based on our far-IR SEDs and fiducial dust temperature evolution (Equation (1)) are presented in Table 3 and range from 1.4 × 10¹¹ L_⊙ to 2.6 × 10¹² L_⊙. Aravena et al. (2020), in a separate analysis of these same sources using SED fits from MAGPHYS, find the range to be 1.1 × 10¹¹ L_⊙ to 3.4 × 10¹² L_⊙. Our derived IR luminosities are just 0.01 dex higher in the mean than those employed by Aravena et al. (2020), demonstrating that the modified blackbody form we utilize here produces IR luminosities very similar to SED analyses that include a mid-IR power law.

The total number of >4σ detections in the z = 1.5–10 galaxies found over the ASPECS footprint is 18. In Section 3.1, we had predicted that 15, 28, and 8 sources would be found from this selection using the consensus IRX–stellar mass relationship, the consensus low-redshift IRX–β relationship, and an SMC-like IRX–β relationship. If in our use of the IRX–β relationship we only consider those sources with stellar masses greater than ${10}^{9.5}$ M_⊙, the predicted number of 4σ detections decreases to 16, almost identical to the observed number. As discussed in Bouwens et al. (2016; their Section 3.1.1) and McLure et al. (2018), the impact of scatter on the breadth of the UV-continuum slope β distribution is to increase the fraction of sources with redder UV-continuum slopes β, increasing the predicted number of sources expected to be detected in the dust continuum.

As in most previous works (Pannella et al. 2009 Bouwens et al. 2016; Dunlop et al. 2017), detected sources from our selection tend to be the star-forming galaxies with the highest stellar masses. In Figure 5 we present the stellar masses and redshifts inferred for the 1362 z = 1.5–10 galaxies over our ASPECS field, indicating which sources are detected in ASPECS. All 11 z ∼ 1.5–3.5 sources with high stellar masses (>10^10.0 M_⊙) and sensitive ALMA observations from ASPECS (<20 μJy beam⁻¹) are detected in our combined data set. If we repeat this exercise on sources in our z = 1.5–10 samples, 11 of 13 are detected, implying a ${85}_{-18}^{+7}$% detection fraction at >10¹⁰ M_⊙.

Zoom In Zoom Out Reset image size

In Figure 6, we present the fraction of sources detected at >4σ as a function of stellar mass. In computing this fraction, we only consider those sources (939 out of 1362) over the ASPECS field where the 1.2 mm continuum sensitivities are the highest, i.e., with 1σ rms noise <20 μJy beam⁻¹. As in previous work (e.g., Bouwens et al. 2016; Dunlop et al. 2017), it is clear that stellar mass is a useful predictor of the dust-continuum flux from star-forming galaxies.

Zoom In Zoom Out Reset image size

Figure 7 shows both the continuum and tentative detections in our sample relative to the stellar mass–β trend found for galaxies in CANDELS (see Section 3.4.1). All 4σ detected sources from ASPECS have a UV-continuum slope β of −1.3 or redder and a stellar mass of ≳${10}^{9.4}$ M_⊙. Detected sources with the largest IRXs (red circles) are distributed toward the reddest UV slopes and highest stellar masses, as expected, but with a significant amount of scatter.

Zoom In Zoom Out Reset image size

Fainter, lower-mass sources in our selections are not sufficiently bright in the dust continuum to be individually detected. It is therefore useful to stack the continuum observations from ASPECS to derive constraints on their dust-continuum properties. We consider various subdivisions of our samples in terms of the physical properties.

For sources included in the stack, the ALMA continuum maps of the relevant sources are mapped onto the same position and stacked in the image plane, weighting each in proportion to the expected 1.2 mm continuum signal divided by the noise squared (per beam). We derive a flux density from the stack based on a convolution of the image stack (33 × 33 aperture) with the primary beam. Individually undetected sources are assumed to be unresolved at the resolution of our observations.

3.3.1. Infrared Excess versus Stellar Mass

We first look at the average IRX of z = 1.5–10 galaxies as a function of stellar mass. We consider six different bins of stellar mass: >10^10.75 M_⊙, ${10}^{10.25}\mbox{--}{10}^{10.75}$ M_⊙, ${10}^{9.75}\mbox{--}{10}^{10.25}$ M_⊙, 10^9.25–10^9.75 M_⊙, ${10}^{8.75}$–10^9.25 M_⊙, and <10^8.75 M_⊙. For these stacks, we weight sources according to the inverse square of the noise [in μJy], i.e., $\sigma {({f}_{1.2\mathrm{mm}})}^{-2}$.

Our stack results are presented in Figure 8 for both our z = 1.5–3.5 and z = 3.5–10 samples, including both the individually detected and undetected sources. Galaxies in our ${10}^{9.75}\mbox{--}{10}^{10.25}$ M_⊙ mass bin are detected at 10σ, while sources in the 10^9.25–10^9.75 M_⊙ bin only show a tentative 2σ detection. Table 4 in the main text and Table 9 from Appendix C presents these results in tabular form. Our stack results for star-forming galaxies which are individually undetected (<4σ) are presented in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Inferred IRX vs. Galaxy Stellar Mass and β from ASPECS (Assuming the Dust Temperature Evolution Specified in Equation (1)) ^a

Stellar		Number of
Mass (M⊙)	β	sources	IRX b
z = 1.5–3.5
$\gt {10}^{10.75}$	All	5	${51.34}_{-21.51}^{+65.82}\pm 1.29$
${10}^{10.25}\mbox{--}{10}^{10.75}$	All	6	${26.99}_{-12.53}^{+27.18}\pm 0.64$
${10}^{9.75}\mbox{--}{10}^{10.25}$	All	11	${16.73}_{-11.72}^{+9.37}\pm 0.51$
${10}^{9.25}\mbox{--}{10}^{9.75}$	All	33	${2.23}_{-0.89}^{+1.17}\pm 0.23$
${10}^{8.75}\mbox{--}{10}^{9.25}$	All	123	${0.90}_{-0.45}^{+0.43}\pm 0.38$
$\lt {10}^{8.75}$	All	467	${0.72}_{-0.80}^{+0.77}\pm 0.66$
z = 3.5–10
$M\gt {10}^{10.25}$	All	1	${19.08}_{-0.00}^{+0.00}\pm 1.02$
${10}^{9.75}\mbox{--}{10}^{10.25}$	All	6	$-{0.22}_{-0.87}^{+0.76}\pm 1.11$
${10}^{9.25}\mbox{--}{10}^{9.75}$	All	31	${4.12}_{-2.38}^{+3.23}\pm 0.49$
${10}^{8.75}\mbox{--}{10}^{9.25}$	All	69	${0.41}_{-0.51}^{+0.50}\pm 0.61$
$\lt {10}^{8.75}$	All	594	$-{0.72}_{-0.66}^{+0.59}\pm 0.59$
z = 1.5–10
$\lt {10}^{9.25}$	All	1253	${0.50}_{-0.35}^{+0.34}\pm 0.31$
z = 1.5–3.5
>${10}^{9.5}$	−4.0 < β <  −1.75	4	${0.02}_{-0.16}^{+0.12}\pm 0.21$
	−1.75 < β <  −1.00	16	${6.54}_{-4.97}^{+4.88}\pm 0.28$
	−1.00 < β <  −0.20	14	${10.27}_{-2.21}^{+3.74}\pm 0.30$
	−0.20 < β	4	${174.57}_{-41.65}^{+104.96}\pm 3.32$
< ${10}^{9.5}$	−4.0 < β <  −1.75	369	${0.83}_{-0.52}^{+0.54}\pm 0.43$
	−1.75 < β <  −1.00	204	${0.84}_{-0.44}^{+0.39}\pm 0.36$
	−1.00 < β	34	${5.57}_{-4.73}^{+6.07}\pm 1.13$
z = 3.5–10
	−4.0 < β <  −1.75	537	$-{0.24}_{-0.48}^{+0.39}\pm 0.37$
	−1.75 < β <  −1.00	125	${0.65}_{-0.54}^{+0.62}\pm 0.56$
	−1.00 < β	32	${7.67}_{-4.98}^{+4.42}\pm 0.96$

Notes.

^a See Tables 9–10 from Appendix C for a more detailed presentation of the stack results summarized here. ^b Both the bootstrap and formal uncertainties are quoted on the result (presented first and second, respectively).

Download table as:  ASCIITypeset image

Our z = 1.5–3.5 stack results provide us with the highest S/N results to derive a dependence of the IRX on stellar mass. In quantifying the dependence, we made use of the power-law relation

$$\begin{eqnarray}&&{\mathrm{IRX}}_{M}={(M/{M}_{s})}^{\alpha }\end{eqnarray} \tag{ 8 }$$

where M_s is the characteristic stellar mass for significant IR emission (L_IR = L_UV) and α gives the power by which the IRX depends on mass. We then fit our z = 1.5–3.5 stacked IRX measurements to this relation and arrived at a best-fit value for M_s and α of ${10}^{{9.15}_{-0.16}^{+0.18}}$ M_⊙ and ${0.97}_{-0.17}^{+0.17}$, respectively. The best-fit relation is shown in both the left and right panels of Figure 8 with the light-red shaded region. Broadly, our z ∼ 1.5–3.5 results are consistent with the consensus relation that we derived in our earlier analysis based on results in the literature (Bouwens et al. 2016).

At z ∼ 3.5–10, our stack results for the IRX show a clear detection in the highest stellar mass bin and a tentative 2σ detections in the third highest stellar mass bin, i.e., 10^9.25–10^9.75 M_⊙, while at lower masses, there is still no detection in our stack results. Our new stack results for the IRXs at z = 3.5–10 seem consistent with what we derive at lower redshifts. Previously, Pannella et al. (2015) had found no strong evidence for evolution in the IRX–stellar mass relation to z ∼ 3.5, and Whitaker et al. (2017) found this same lack of evolution to z ∼ 3. From first principles, one might expect some evolution in this relationship due to the observed evolution in the mass–metallicity relation (e.g., Erb et al. 2006a); however, it is possible that a higher gas and ISM mass in $z\gtrsim 2$ galaxies compensate for the lower metal content to produce a relatively unevolving IRX–stellar mass relation (Tan et al. 2014).

However, we emphasize that this conclusion is sensitive to the dust temperature evolution we adopt. If there is no significant evolution in the dust temperatures with redshift, then the IRXs at z = 3.5–10 would be lower by ∼0.4 dex than what we infer at z = 1.5–3.5, and we would therefore infer that the IRX–stellar mass relation increases at early cosmic times. In Appendix D, we investigated the extent to which our IRX versus stellar mass relation showed a dependence on the stellar population code used to estimate the mass for individual sources, and recovered a steeper IRX–stellar mass relation using Prospector masses.

For stacks of sources with stellar masses less than ${10}^{9.25}$ M_⊙, we do not find a detection in the IR continuum. In an effort to provide a dramatic illustration of this, we include in Figure 10 three different stacks of all 1253 z = 1.5–10 sources with stellar mass estimates <10^9.25 M_⊙ over our ASPECS footprint. Our first stack weights sources by their UV flux, our second stack weights sources by their estimated stellar mass, and our third stack weights sources equally (left, center, and right panels, respectively). None of the stacks show a significant detection, and in our unweighted stack, the mean continuum flux density is −0.1 ± 0.4 ± 0.4 μJy beam⁻¹. Even weighting sources in the stack by the measured UV-continuum slope β fails to result in a significant detection. This demonstrates, rather dramatically, that faint, UV-selected galaxies show essentially no dust-continuum emission (see also Carvajal et al. 2020). Converting this flux density constraint to a SFR for a galaxy at z ∼ 4, we derive a SFR of 0.0 ± 0.1 M_⊙ yr⁻¹.

Zoom In Zoom Out Reset image size

3.3.2. Infrared Excess versus β

Stacked results of z = 1.5–3.5 and z = 3.5–10 sources over our ASPECS footprint are presented as a function of UV-continuum slope β in Figure 11 with the large solid circles and 2σ upper limits. Five different bins in β are utilized to better map out the trend with UV-continuum slope β.

Zoom In Zoom Out Reset image size

Separate stack results are presented for sources with stellar masses $\gt {10}^{9.5}$ M_⊙ (large red circles and downward arrows, respectively) and $\lt {10}^{9.5}$ M_⊙ (large green circles and downward arrows, respectively) to evaluate whether higher-mass galaxies show a different IRX–β relationship from lower-mass galaxies. This treatment also ensures that results in the redder, high-mass bins are not impacted by the inclusion of bluer, lower-mass sources (but where the measured UV-continuum slopes β are much redder than the actual slopes due to the impact of noise). Figure 12 presents our stack results for star-forming galaxies which are individually undetected (<4σ). Our IRX–β stack results are presented in Table 4 in the main text and Table 11 in Appendix C.

Zoom In Zoom Out Reset image size

For our highest-mass z ∼ 1.5–3.5 samples, our stack results lie closest to the Reddy (Calzetti-like) IRX–β relations. As in our earlier analysis of the ASPECS pilot data, we formalize this analysis by finding those parameters which best match the stacked IRX results versus β and then computing 68% confidence intervals on the derived parameters. Here we derive constraints on both ${{dA}}_{\mathrm{UV}}/d\beta$ and β as

$$\begin{eqnarray}&&{\mathrm{IRX}}_{\beta }=1.7\times {10}^{0.4({{dA}}_{\mathrm{UV}}/d\beta )(\beta -{\beta }_{\mathrm{int}})}-1,\end{eqnarray} \tag{ 9 }$$

instead of just deriving constraints on ${{dA}}_{\mathrm{UV}}/d\beta$ as in our previous analysis.

Our maximum-likelihood derived values for ${{dA}}_{\mathrm{UV}}/d\beta$ and ${\beta }_{\mathrm{int}}$ are ${1.81}_{-0.14}^{+0.18}$ and $-{1.86}_{-0.10}^{+0.14}$ and presented in Table 5. The ${{dA}}_{\mathrm{UV}}/d\beta$ we derive is similar to the Calzetti or Reddy value, i.e., 1.97 or 1.84. Meanwhile, the ${\beta }_{\mathrm{int}}=-1.86$ we derive is not only redder than the ${\beta }_{\mathrm{int}}=-2.23$ implicit in the Meurer et al. (1999) formulation, but also redder than what might be expected for dust-free galaxies with a constant SFR for 100–500 Myr (e.g., as in Reddy et al. 2018). Both the ${{dA}}_{\mathrm{UV}}/d\beta$ and ${\beta }_{\mathrm{int}}$ we derive are consistent with the consensus low-redshift values for these quantities (e.g., Equation (4)). If we instead take ${\beta }_{\mathrm{int}}=-2.23$ as has been conventional (following Meurer et al. 1999), the ${{dA}}_{\mathrm{UV}}/d\beta$ we recover is ${1.48}_{-0.11}^{+0.09}$. In our pilot study, our best-fit determination for ${{dA}}_{\mathrm{UV}}/d\beta$ is ${1.26}_{-0.36}^{+0.27}$ when taking ${\beta }_{\mathrm{int}}$ equal to −2.23. For a ${\beta }_{\mathrm{int}}=-2.30$, we recover ${{dA}}_{\mathrm{UV}}/d\beta$ equal to ${1.42}_{-0.11}^{+0.09}$.

Table 5. Present Constraints on the IRX–β Relationship

Sample	Mass Range	${{dA}}_{\mathrm{UV}}/d\beta$	${\beta }_{\mathrm{int}}$
Current Determinations
z ∼ 1.5–3.5	$\gt {10}^{9.5}$ M⊙	${1.81}_{-0.14}^{+0.18}$	$-{1.86}_{-0.10}^{+0.14}$
z ∼ 1.5–3.5	$\lt {10}^{9.5}$ M⊙	${1.12}_{-0.30}^{+0.31}$	−1.86 (fixed)
z ∼ 1.5–3.5	$\gt {10}^{9.5}$ M⊙	${1.48}_{-0.11}^{+0.09}$	−2.23 (fixed)
z ∼ 1.5–3.5	$\gt {10}^{9.5}$ M⊙	${1.42}_{-0.11}^{+0.09}$	−2.30 (fixed)
Canonical IRX–β Relations
Consensus: z ∼ 0 a		1.86	−1.87
Reddy et al. (2015): z ∼ 2		1.84	−2.43
Overzier et al. (2011): z ∼ 0		1.96	−1.96
Takeuchi et al. (2012): z ∼ 0		1.58	−1.94
Casey et al. (2014): z ∼ 0		2.04	−1.64
Meurer et al. (1999): z ∼ 0		1.99	−2.23
Dust Laws
Calzetti		1.97	⋯
SMC		∼1.10	⋯

Note.

^a Taking the median of the IRX–β relations derived by Overzier et al. (2011), Takeuchi et al. (2012), and Casey et al. (2014). See Appendix B.

Download table as:  ASCIITypeset image

For lower-mass (<10^9.5 M_⊙) z ∼ 1.5–3.5 galaxies found over ASPECS, significant ALMA continuum flux is found in two of the three β bins we consider. Fixing ${\beta }_{\mathrm{int}}$ to be the same as for the higher-mass galaxies, we find a best-fit value for ${{dA}}_{\mathrm{UV}}/d\beta$ of ${1.12}_{-0.30}^{+0.31}$. This is most consistent with an SMC-like dust curve, but is nevertheless consistent with our constraints on ${{dA}}_{\mathrm{UV}}/d\beta$ value in the higher mass >10^9.5 M_⊙ bin.

We now look at the constraints we can set on the IRX–β relationship at z ∼ 3.5–10. We focus on sources with the highest stellar masses, i.e., >10^9.25 M_⊙ to minimize the impact of intrinsically blue, lower-mass sources scattering to redder colors (see Section 3.1.1 from Bouwens et al. 2016). Our z ∼ 3.5–10 stack results for sources shows prominent detections in the reddest two β bins, one at −0.8 and 1.6. Those two detections imply very different IRX–β relationships. Fixing the value of ${\beta }_{\mathrm{int}}$ to be −2.23 and fitting to two bluest β bins plus the $\beta \sim -0.8$ bin, we derive a ${{dA}}_{\mathrm{UV}}/d\beta$ value of 2.27. By contrast, if we fit to the two bluest β bins plus the $\beta \sim 1.6$ bin, we derive a ${{dA}}_{\mathrm{UV}}/d\beta$ value of 0.63. Given how different the two relations are and the fact that there are only two significant detections at z > 3.5 we can use from ASPECS, perhaps it is best for us simply to quote our z = 3.5–10 results as the range spanned by these two relations. As this range includes both Reddy/Calzetti-like and SMC-like dust relations, the ASPECS data provides us with very little information on the IRX–β relation evolution.

3.3.3. Summary of Stack Results

3.4.1. Correlation with Stellar Mass and UV-continuum Slope β

Having looked at the correlation of the IRX with the stellar mass and UV-continuum slope β, it is interesting to try to link these relations based on both the empirical correlation of these two quantities, and on the large samples that now exist from various legacy data sets. Given the significant correlation between the dust content and the metallicity of galaxies and their stellar mass (e.g., Reddy et al. 2010; Pannella et al. 2015), one would expect a strong correlation between the UV-continuum slope of galaxies and their stellar mass, as in fact is observed (e.g., McLure et al. 2018; Carvajal et al. 2020).

For this exercise, we take all the z = 1.3–2.5 sources identified over the five CANDELS fields by the 3D-HST team (Skelton et al. 2014) and compare their UV-continuum slopes β with their stellar masses derived by prospector (Leja et al. 2017, 2019). The results are presented in Figure 7, and it is clear that for sources with stellar masses to ${10}^{8.8}$ M_⊙ the UV-continuum slopes β of galaxies generally lie in the range −2.5 to −1.8. For sources with stellar masses >10⁹, the UV-continuum slopes β show a strong correlation with stellar mass to 10¹¹ M_⊙.

Using the correlations we derive between the IRX and the stellar mass (Section 3.3.1),

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\mathrm{IRX}=\alpha {\mathrm{log}}_{10}(M/{M}_{s})\end{eqnarray} \tag{ 10 }$$

and between the IRX and the UV-continuum slope β (Section 3.3.2)

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\mathrm{IRX}={\mathrm{log}}_{10}\left({10}^{0.4\left(\tfrac{{{dA}}_{\mathrm{FUV}}}{d\beta }\left(\beta -{\beta }_{\mathrm{int}}\right)\right)}-1\right)+0.23.\end{eqnarray} \tag{ 11 }$$

This results in

$$\begin{eqnarray}&&\beta ={\beta }_{\mathrm{int}}+\displaystyle \frac{2.5}{\tfrac{{{dA}}_{\mathrm{FUV}}}{d\beta }}{\mathrm{log}}_{10}\left(\displaystyle \frac{1}{1.7}{\left(M/{M}_{s}\right)}^{\alpha }+1\right).\end{eqnarray} \tag{ 12 }$$

Fixing ${\beta }_{\mathrm{int}}=-2.3$ and taking the best-fit value we find for $\tfrac{{{dA}}_{\mathrm{FUV}}}{d\beta }$ (i.e., 1.42), we look for the optimal values of M_s and α to capture the observed relationship between stellar mass and UV-continuum slope β shown in Figure 7. In deriving this relationship, we segregate sources into those above and below the β versus M relation, determine the number of such sources in six distinct regions along the relation, compute the square of the difference in the number of sources on each side for each of the six regions, and then minimize the square of the differences. The best-fit values of M_s and α are ${10}^{9.07}$ M_⊙ and 0.92, respectively. This best-fit relation is included in Figure 7 as the blue line. For comparison, Figure 7 also shows the β versus stellar mass relationship derived by McLure et al. (2018). Encouragingly enough, the best-fit value for M_s and α are consistent (at 1σ) with the values we derive from our IRX–stellar mass analysis, i.e., ${10}^{{9.15}_{-0.16}^{+0.18}}$ M_⊙ and ${0.97}_{-0.17}^{+0.17}$, respectively, demonstrating that the IRX–β and IRX–stellar mass relations we derive are essentially equivalent.

3.4.2. Infrared Excess of a Function of Stellar Mass and UV-continuum Slope β

Having quantified the approximate relationship between the stellar mass and UV-continuum slope β of galaxies at z ∼ 1.5–2.5, we now move on to try to express the IRX as a bivariate function of the UV-continuum slope β and the stellar mass M.

One reason for pursuing such a parameterization would be to take advantage of the greater information content present in both the measured UV-continuum slope β and the inferred stellar mass of a galaxy. While the two parameters are clearly correlated (e.g., Section 3.4.1), they do provide us with independent information on sources and therefore should theoretically be able to improve our estimates of the IRX.

We use the following functional form:

$$\begin{eqnarray}&&\mathrm{IRX}{(\beta ,M)=1.7({10}^{0.4({{dA}}_{\mathrm{UV}}/d\beta )(\beta +2.3)}-1)(M/M(\beta ))}^{\alpha }\end{eqnarray} \tag{ 13 }$$

where $M(\beta )$ is as follows and gives the expected stellar mass for a given UV-continuum slope (as derived in the previous subsection):

$$\begin{eqnarray}&&M{(\beta )=({10}^{9.07}{M}_{\odot })(1.7\times {10}^{0.4(1.42)(\beta +2.3)}-1)}^{1/0.92}.\end{eqnarray} \tag{ 14 }$$

The expression we adopt for $\mathrm{IRX}(\beta ,M)$ is the standard form for the IRX–β relation, but then allows for a dependence on whether a source is more or less massive than one would expect for a given UV-continuum slope β.

Sources from ASPECS were divided in stellar mass and β in the same way as in the previous subsections, stacked using the same weighting scheme as described in Section 3.3, and then an average IRX derived for each stellar mass–β bin. The derived IRXs versus β and stellar mass were then fit using the expression given in Equation (13). The best-fit values we recovered for ${{dA}}_{\mathrm{UV}}/d\beta$ and α were 1.48 ± 0.10 and 0.67 ± 0.06. Encouragingly enough, the best-fit value for ${{dA}}_{\mathrm{UV}}/d\beta$ is very similar to what we found expressing the IRX as a function of the UV-continuum slope β alone. We do find a minor additional dependence on whether the inferred stellar mass is greater or less than given by the general correlation between stellar mass and β, but the dependence is not particularly strong. The blue lines in Figure 13 presents the suggested regions in β/M_* parameter space with IRXs of 4, 20, and 100, shown relative to the detected and undetected sources from ASPECS.

Zoom In Zoom Out Reset image size

Álvarez-Márquez et al. (2019) had previously attempted to quantify the IRX as a function of both the UV-continuum slope β and stellar mass, as ${\mathrm{log}}_{10}(\mathrm{IRX})=(0.51\pm 0.06){\beta }_{\mathrm{UV}}$ $+\,(0.37\pm 0.08){\rm{log}}({M}_{\ast }[{M}_{\odot }])-1.89\pm 0.40$. While the functional form Álvarez-Márquez et al. (2019) utilize is different from what we consider, it is interesting to try to compute the logarithmic dependence of IRX on ${\beta }_{\mathrm{UV}}$ and log₁₀ M_* to investigate how similar the results are. For simplicity, we compute the dependence at a β = 0.5 and ${\mathrm{log}}_{10}{M}_{* }$ of ${10}^{10.5}$ M_⊙. For the $\mathrm{IRX}(\beta ,M)$ function we derive, we compute a $d{\mathrm{log}}_{10}(\mathrm{IRX})/d\beta$ of 0.18 and a $d{\mathrm{log}}_{10}(\mathrm{IRX})/d{\mathrm{log}}_{10}{M}_{* }$ of 0.67 versus 0.51 ± 0.06 and 0.37 ± 0.08 found by Álvarez-Márquez et al. (2019). These relations are in reasonably good agreement, which is encouraging given the differences in approach (the Álvarez-Márquez et al. (2019) relations are based on deep Herschel stacks).

Given the strong correlation between both parameters, where ${\rm{\Delta }}\beta \sim 1.5{\rm{\Delta }}{M}_{* }$ (see Section 3.4.1), it is also interesting to reformulate the Álvarez-Márquez et al. (2019) IRX relation to be just a single function of β. We find $d{\mathrm{log}}_{10}(\mathrm{IRX})/d\beta \sim 0.63$. If we make same change to our bivariate $\mathrm{IRX}(\beta ,M)$ relation, we find $d{\mathrm{log}}_{10}(\mathrm{IRX})/d\beta \sim 0.68$. As with the previous comparison, the two dependencies are similar, which is encouraging given differences in the two approaches.

Before concluding this section, it is useful to summarize the predicted 1.2mm flux densities expected for different z ≳ 1.5 galaxies over the ASPECS footprint and compare those predictions with the observations. A compilation of the results are presented in Table 6 and include the predicted flux densities using (1) the Meurer et al. (1999) IRX–β relation (Equation (B4): Appendix B), (2) the consensus low-redshift IRX–β relation (Equation (4)) derived here in Appendix B from literature results, (3) an SMC-like IRX–β relation (Equation (6)), (4) the consensus IRX–stellar mass relation (Equation (3)) presented in our previous study (Bouwens et al. 2016), (5) our derived IRX–β relation for >10^9.5 M_⊙, z ∼ 1.5–3.5 galaxies (Equation (9): Section 3.3.2), (6) our derived IRX–stellar mass relation for z ∼ 1.5–3.5 galaxies (Equation (8): Section 3.3.1), and (7) our derived IRX(β, M) relation (Equation (13): Section 3.4.2). As one final predictor, we include a comparison against the flux density predicted taking the geometric mean of our derived z = 1.5–3.5 IRX–β relation and our IRX–stellar mass relationship, i.e., ${({\mathrm{IRX}}_{\beta }{\mathrm{IRX}}_{M})}^{1/2}$, and using Equations (9) and (8) while taking ${{dA}}_{\mathrm{UV}}/d\beta$, ${\beta }_{\mathrm{int}}$, M_s , and α to be ${1.81}_{-0.14}^{+0.18}$, $-{1.86}_{-0.10}^{+0.14}$, ${10}^{{9.15}_{-0.16}^{+0.18}}$ M_⊙, and ${0.97}_{-0.17}^{+0.17}$, respectively. This should provide for an alternate way of using both the UV-continuum slopes β and stellar masses in estimating the IRX.

Table 6. Comparisons between the Predicted and Measured 1.2 μm Flux Densities for z ≳ 1.5 UV-selected Galaxies Showing 4σ Detections ⁱ

	Predicted f1.2mm (μJy)								Measured f1.2mm
ID	${\mathrm{IRX}}_{M99}$ a	${\mathrm{IRX}}_{z=0}$ b	IRXSMC c	${\mathrm{IRX}}_{M,0}$ d	${\mathrm{IRX}}_{\beta }$ e	IRXM f	$\mathrm{IRX}(\beta$,M) g	$({{\rm{IRX}}}_{\beta }\,{{\rm{IRX}}}_{M}{)}^{1/2}$ h	(μJy)
XDFU-2435246390	60	24	7	63	23	58	50	36	1071 ± 46
XDFU-2385446340	184	66	31	111	65	110	126	85	752 ± 10
XDFU-2397246112	380	151	45	642	143	587	436	290	461 ± 14
XDFU-2369747272	1572	534	56	43	469	40	82	137	432 ± 9
XDFU-2400547554	1434	571	177	1492	541	1391	1206	868	342 ± 18
XDFU-2410746315	41	16	6	3	15	3	6	7	316 ± 11
XDFU-2433446471	173	64	11	44	58	40	47	48	289 ± 21
XDFU-2350746475	710	264	47	148	240	137	170	182	233 ± 11
XDFU-2416846554	294	109	19	21	99	20	34	44	185 ± 10
XDFB-2380246263	262940	73791	2555	514	60164	480	1716	5373	163 ± 10
XDFB-2355547038	124	49	18	11	47	12	22	23	155 ± 9
XDFU-2387248103	203	81	26	10	77	10	22	28	134 ± 24
XDFU-2373546453	655	261	91	474	250	452	451	336	107 ± 10
XDFU4596	459	183	61	259	174	248	263	208	97 ± 9
XDFU-2361746276	552	218	60	225	205	213	238	209	85 ± 12
XDFU9838	253	100	28	56	94	54	73	71	65 ± 15
XDFU-2359847256	145	56	23	123	54	119	119	80	58 ± 10
XDFU-2370746171	244	85	40	67	84	69	100	76	47 ± 11
Performance i
$({f}_{\mathrm{obs}}-{f}_{\mathrm{pred}})/{{ef}}_{\mathrm{obs}}$
25%/75% Quartiles	[−23.2,−1.7]	[−8.0,0.8]	[−2.0,4.0]	[−4.8,1.1]	[−7.2,0.8]	[−4.8,1.0]	[−5.5,0.9]	[−3.8,1.3]
${({f}_{\mathrm{obs}}-{f}_{\mathrm{pred}})/({f}_{\mathrm{pred}}^{2}+{{ef}}_{\mathrm{obs}}^{2})}^{0.5}$
25%/75% Quartiles	[−1.1,−0.4]	[−1.1,0.6]	[−1.1,1.4]	[−1.1,0.7]	[−1.1,0.7]	[−1.1,0.7]	[−1.1,0.3]	[−1.1,0.6]
${\mathrm{log}}_{10}({f}_{\mathrm{obs}}/{f}_{\mathrm{pred}})$ j
Mean/Standard Deviation	−0.42 ± 0.81	0.01 ± 0.80	0.54 ± 0.70	0.29 ± 0.69	0.03 ± 0.79	0.30 ± 0.68	0.14 ± 0.65	0.17 ± 0.67
Median	−0.59	−0.19	0.37	0.12	−0.16	0.12	−0.05	−0.02

Notes.

^a From Equation (B4), which is the Meurer et al. (1999) IRX–β relationship. ^b From Equation (4), which is the consensus low-redshift IRX–β relation derived here in Appendix B from literature results. ^c Equation (6), which gives an SMC-like IRX–β relation. ^d From Equation (3), which is the consensus IRX–stellar mass relation presented in our previous study (Bouwens et al. 2016). ^e From Equation (9), which is the IRX–β relation we derived for >10^9.5 M_⊙, z ∼ 1.5–3.5 galaxies (Section 3.3.2). ^f From Equation (8), which is the IRX–stellar mass relation we derived for z ∼ 1.5–3.5 galaxies (Section 3.3.1). ^g From Equation (13), which is the IRX(β, M) relation we derived (Section 3.4.2). ^h Geometric mean of our derived z = 1.5–3.5 IRX–β relation ${\mathrm{IRX}}_{\beta }$ and our IRX–stellar mass relationship IRX_M . ⁱ See Section 3.5 for a discussion. ^j Only for those 25 sources where ${f}_{\mathrm{obs}}/{{ef}}_{\mathrm{obs}}\gt 2$.

Download table as:  ASCIITypeset image

The observed fluxes are also explicitly compared against these many estimators in Figure 14. A quantification of the mean, median, and 1σ scatter in the logarithmic ratio of the predicted and measured 1.2mm flux densities is presented in Table 6, and it is clear there is substantial scatter between the observed and predicted flux densities. The scatter ranges from 0.65–0.81 dex, with the smallest dispersion found for the $\mathrm{IRX}(\beta ,M)$ and ${({\mathrm{IRX}}_{\beta }{\mathrm{IRX}}_{M})}^{1/2}$ estimators, with only slight increases in the dispersion for the other relations. The $\mathrm{IRX}(\beta ,M)$ and ${({\mathrm{IRX}}_{\beta }{\mathrm{IRX}}_{M})}^{1/2}$ estimators also provide the best predictions of the observed flux densities in the median.

Zoom In Zoom Out Reset image size

As a separate means of evaluating the estimators, we compare the predicted 1.2mm flux densities from these estimators with the measured flux densities using both the detected sources in Table 6 and sources expected to be detected at >2σ averaging the IRX–β and IRX–stellar mass relations derived here (Equations (9) and (8)), i.e., 70 sources in total. For each of these sources, we computed the difference between the measured and predicted flux for each source, i.e., f_obs and f_pred, divided the result by the measurement error ef_obs, and then determined the average as well as the upper and lower quartiles. For almost every estimator, the difference between the upper and lower quartiles is larger than the measurement error by ≳5×.

For each of the estimators, we also computed the differences between the measured and predicted flux densities for the same sources as the previous exercise, divided the result by the rms of the predicted flux densities and flux measurement uncertainties, and finally computed the upper and lower quartiles. This should give an approximate relative uncertainty on the flux density predictions. All of our estimators perform comparably well, with only modest differences between them.

In summary, as with previous work (e.g., Meurer et al. 1999; Reddy et al. 2006), estimators of the IRX tend to be accurate in predicting the obscured SFRs or IR luminosities for the average source and tend to show at least ∼0.65 dex scatter for individual sources. Of those we consider, the different estimators for the IRX all perform comparably, with marginally better performance for the estimators that consider both mass M and β, i.e., $\mathrm{IRX}(\beta ,M)$ and ${({\mathrm{IRX}}_{\beta }{\mathrm{IRX}}_{M})}^{1/2}$, and the ${\mathrm{IRX}}_{M99}$ estimator predicting the IRXs the least well.

It is interesting to compare the present set of ALMA continuum detections to those that were previously reported over the HUDF by Aravena et al. (2016), Bouwens et al. (2016), and Dunlop et al. (2017). The reported detections and tentative detections by Aravena et al. (2016) and Bouwens et al. (2016) made use of the 1 arcmin² pilot for ASPECS, while the Dunlop et al. (2017) results were based on the 1.3mm ALMA continuum observations they obtained over a 4.5 arcmin² region within the HUDF/XDF.

Using the 1 arcmin² pilot observations for ASPECS, Aravena et al. (2016) and Bouwens et al. (2016) detected five $z\gt 1.5$ galaxies and reported tentative detections for three more $z\gt 1.5$ galaxies. Our new observations confirm all of our previously claimed detections at >4σ, making it clear that those detections were real. In addition, one of the tentatively detected sources from our pilot program, i.e., XDFU-2370746171, shows a >4σ detection (40 ± 11 μJy beam⁻¹) in the new data, confirming that the reported tentative detection (34 ± 14 μJy beam⁻¹) from our pilot was real.

The measured flux densities for the two other tentative detections from our pilot, i.e., XDFU-2365446123 and XDFU-2384246384, are −27 ± 17 μJy beam⁻¹ and 8 ± 10 μJy beam⁻¹ versus our measurements of 38 ± 16 μJy beam⁻¹ and 36 ± 14 μJy beam⁻¹, respectively, in the pilot for these sources. Combining the measurements, the flux is 7 ± 12 μJy beam⁻¹ for XDFU-2365446123 and 17 ± 8 μJy beam⁻¹ for XDFU-2384246384. While the new observations do not support the reality of either source, XDFU-2384246384 still shows a tentative 2.1σ detection in the continuum in the combined data set and thus may be real.

In the Dunlop et al. (2017) search, 16 dust-continuum (>3.5σ) detections are identified, 11 of which have an estimated redshift in excess of 1.5 and lie within the ASPECS footprint. 8 of these 11 sources are clearly confirmed with our ASPECS ALMA observations. For the 3 reported continuum detections from the Dunlop et al. (2017) which are not unambiguously confirmed by our ASPECS observations, we measure $-9\pm 21\mu$ Jy (UDF9), $-45\pm 31\mu$ Jy (UDF12), and $-3\pm 9\mu$ Jy (UDF15).

It is interesting to compare the IRX–stellar mass and IRX–β relations we derive with the many previous determinations in the literature. We focus on determinations at z ∼ 1.5–3.5 since this is where our results are the most significant and where most of previous results have been obtained. In Figure 15, we compare the IRX–stellar mass relationship we find at z ∼ 1.5–3.5 with what we obtained in our pilot study (Bouwens et al. 2016) and many other determinations in the literature (Heinis et al. 2014; Pannella et al. 2015; Álvarez-Márquez et al. 2016, 2019; Bourne et al. 2017; Fudamoto et al. 2017, 2020a; Koprowski et al. 2018; McLure et al. 2018; Reddy et al. 2018).

Zoom In Zoom Out Reset image size

Overall, our new IRX–stellar mass results appear to be in agreement with previous results as presented, e.g., by Heinis et al. (2014), Pannella et al. (2015), Bourne et al. (2017), and McLure et al. (2018), or even as given by the consensus relation derived in our pilot study (shown with the gray line). Our best-fit IRX–stellar mass correlation is ∼0.2–0.3 dex higher at 10¹⁰ M_⊙ than found in our earlier study (Bouwens et al. 2016) but consistent within the quoted uncertainties. Thanks to the larger number of dust-continuum detected sources in the current ASPECS study versus our pilot study (18 versus three 4σ detections), we are able to significantly improve our quantification of the IRX–stellar mass relation relative to our previous study.

The slope recovered for our new IRX–stellar mass relation, i.e., ${0.97}_{-0.17}^{+0.17}$, is very close to 1. We had previous adopted a value of unity in Bouwens et al. (2016) for the consensus relation (Equation (3)) based on the IRX–stellar mass results of Reddy et al. (2010), Whitaker et al. (2014), and Álvarez-Márquez et al. (2016). The IRX–stellar mass relation derived by McLure et al. (2018) using the shallower ALMA observations over the HUDF (Dunlop et al. 2017) also find a slope (0.85 ± 0.05), very close to what we find here. At one other extreme, Fudamoto et al. (2020a) recover a much steeper slope (1.64 ± 0.10) for the IRX–stellar mass relation, similar to what we derive using Prospector for our stellar mass estimates (Appendix D). Meanwhile, earlier results obtained from an analysis of Herschel data by Pannella et al. (2015) find a much shallower IRX–stellar mass relation, with a slope of ∼0.64, clearly shallower than what we find here (see also results by Álvarez-Márquez et al. 2019). Given the current strong constraints on the obscured SFR at low masses (<10^9.25 M_⊙) and the challenge that source confusion presents for the lowest mass sources with Herschel, it seems likely that the slope of the IRX is approximately unity or steeper, as essentially all analyses relying on ALMA data have found.

The IRX–stellar mass results we obtain at z ∼ 3.5–10 can be compared with results obtained using a small sample of bright z ∼ 5–6 galaxies from Capak et al. (2015) and Willott et al. (2015) and assuming the dust temperature evolution given in Equation (1). Also included in this comparison are the new ALPINE results from Fudamoto et al. (2020b), both as quoted in the original study (solid colored points) and adopting the fiducial dust temperature evolution adopted here (Equation 1). This comparison is presented in Figure 16. Our own results appear to be most consistent with the consensus IRX–M_* relationship we had derived in our pilot study (Bouwens et al. 2016) and now derived here as z ∼ 1.5–3.5. While this suggests that the IRX–stellar mass relation may extend to z ∼ 5–6 with little or no evolution, the ASPECS field only contains a few bright, massive sources to probe this well. Additionally, this inference depends critically on the dust temperature being relatively high, i.e., ∼50 K, at z ∼ 4–6. If the temperature is instead ∼41 K as Fudamoto et al. (2020b) adopt in their analysis, clearly the IRX–stellar mass relation at z > 3.5 is lower than what is found at z ∼ 1.5–3.5.

Zoom In Zoom Out Reset image size

In Figure 17, we compare the IRX–β relationship we derive for higher-mass, z ∼ 1.5–3.5 galaxies with the results obtained in our pilot study (Bouwens et al. 2016) as well as a wide variety of different determinations in the literature (Heinis et al. 2013; Álvarez-Márquez et al. 2016, 2019; Bourne et al. 2017; Fudamoto et al. 2017, 2020a; Koprowski et al. 2018; McLure et al. 2018; Reddy et al. 2018). Similar to what we found for the IRX–stellar mass relation, the larger number of dust-continuum detections found here (versus from the smaller-area ASPECS pilot) results in the recovery of a steeper IRX–β relation than in our pilot, i.e., ${1.48}_{-0.11}^{+0.09}$ versus ${1.26}_{-0.36}^{+0.26}$ when fixing ${\beta }_{\mathrm{int}}=-2.23$. The only apparently significant difference occurs for our determination at −1.3 where the limit from our pilot program was ${1.31}_{-0.94}^{+0.67}\pm 0.72$ (at $\beta \sim -1.4$) and where our new measurement is ${6.54}_{-4.97}^{+4.88}\pm 0.28$ (at $\beta \sim -1.2$). This difference results both from the larger number of dust-detected sources in the four times larger area probed by ASPECS (versus our pilot) and from our change of the β binning scheme to exploit the larger number of sources to improve our leverage for constraining the IRX–β relation.

Zoom In Zoom Out Reset image size

Relative to various determinations from the literature, the most significant differences occur for the bluest values of β, i.e., $\beta \sim -1.8$, where our own determination of the IRX is some 0.2–1.0 dex lower than the determinations of Reddy et al. (2018), Fudamoto et al. (2017), Bourne et al. (2017), and McLure et al. (2018). It seems likely that the differences here are due to the presence of blue, IR-luminous sources in many previous selections. While blue, IR-luminous galaxies are known to exist (e.g., Reddy et al. 2006; Casey et al. 2014), especially at high IR luminosities (>10¹² L_⊙) where there is less connection between the UV and IR morphologies in galaxies, these sources are not sufficiently common to be well sampled by the ∼2.5 × 10⁴ comoving Mpc³ volume probed by ASPECS at z = 1.5–3.5.

Otherwise, our IRX–β results are broadly in agreement with the results of Reddy et al. (2018), Álvarez-Márquez et al. (2016), and Heinis et al. (2013). For redder values of β, our IRX–β results are lower than the results of McLure et al. (2018); Fudamoto et al. (2017), Bourne et al. (2017), and Fudamoto et al. (2020a) by ∼0.4 dex. We expect that some fraction of these differences, i.e., 0.3 dex, could result from different calibrations used to derive the IR luminosities and obscured SFRs from the measured ALMA fluxes (e.g., Murphy et al. 2011 versus Whitaker et al. 2017).

The purpose of this section is to take advantage of the results of our analyses from the previous subsections to derive dust corrections that we can apply to the general star-forming galaxy population at $z\geqslant 3.5$.

We will focus on deriving these corrections as a function of the UV luminosity of galaxies and derive a distribution of dust corrections that make up each UV luminosity bin. To ensure a significant sampling of each UV luminosity bin, we leverage the large selections of star-forming galaxies Bouwens et al. (2015) identified at z ∼ 4, 5, 6, 7, 8, and 10 over the CANDELS GOODS-North and GOODS-South.

Each of the sources over the CANDELS GOODS-North and GOODS-South fields has sensitive HST optical/ACS and WFC3/IR photometry available to derive UV-continuum slopes for each source in these samples. Another valuable aspect of sources in these fields is the deep Spitzer/IRAC observations that exist from the 200 hr GREATS program (M. Stefanon et al. 2020, in preparation) to provide rest-optical photometry for z ∼ 4–8 galaxies and thus to estimate stellar masses. HST and Spitzer/IRAC photometry is performed on sources in these fields in a similar way to that described in Section 2.2, and UV-continuum slopes β and stellar masses are estimated using the FAST stellar population fitting code as described in Section 2.5.

In deriving dust corrections for each bin in UV luminosity, we make use of the stellar masses and UV-continuum slopes β derived for our large CANDELS samples and utilize the new relation in Equation (13) we derived in Section 3.4.2 for the IRX expressed as a function of both β and stellar mass M_*. To ensure that our extinction estimates are not overly impacted by noise in the photometry scattering lower-mass sources to red β measurements, we force the IRXs of sources with stellar masses less than 10⁹ M_⊙ to be zero, consistent with our derived observational constraints.

For convenience, we present the dust corrections we have derived here in Table 7. If the dust temperatures of z > 3 galaxies are in fact closer to 35 K than given by our fiducial dust temperature model, the dust correction we compute would be approximately half as large. As in Bouwens et al. (2016), we assume that the average dust correction for UV bright (<25.5 mag) galaxies at z ∼ 3 is ∼5 following the findings of Reddy & Steidel (2004).

Table 7. Estimated Dust Corrections to Apply to the UV Luminosity Density Results Integrated to Various Limiting Luminosities

	${{\rm{log}}}_{10}$ Dust Correction
Sample	(>0.05 ${L}_{z=3}^{* }$) b	(>0.03 ${L}_{z=3}^{* }$) b
z ∼ 3	0.37 a	0.34 a
z ∼ 4	0.33	0.31
z ∼ 5	0.30	0.27
z ∼ 6	0.20	0.17
z ∼ 7	0.09	0.07
z ∼ 8	0.07	0.06

Notes.

^a For uniquely the z ∼ 3 sample, we make use of the finding by, e.g., Reddy & Steidel (2004) and Reddy et al. (2010) that the average IRX for galaxies brighter than 25.5 mag at z ∼ 3 is a factor of ∼5. ^b The specified limits 0.05 ${L}_{z=3}^{* }$ and 0.03 ${L}_{z=3}^{* }$ correspond to faint-end limits of −17.7 and −17.0, respectively, which is the limiting luminosity to which z ∼ 7 and z ∼ 10 galaxies can be found in current probes (Schenker et al. 2013; McLure et al. 2013; Ellis et al. 2013; Oesch et al. 2013; Bouwens et al. 2015).

Download table as:  ASCIITypeset image

As in the analysis for our pilot study, we apply the dust corrections we derive in the previous subsection to the UV luminosity densities integrating the UV luminosity function (LF) of Bouwens et al. (2015) to 0.05 ${L}_{z=3}^{* }$ (−17.7 mag) and to 0.03 ${L}_{z=3}^{* }$ (−17.0 mag). As in previous work, the UV luminosity densities are converted into SFR densities using Madau & Dickinson (2014) conversion factor $\kappa \,=1.15\,\times {10}^{-28}{M}_{\odot }\,{{\rm{yr}}}^{-1}{\rm{e}}{\rm{r}}{{\rm{g}}}^{-1}\,{\rm{s}}{\rm{H}}{\rm{z}}$ (see also Madau et al. 1998 and Kennicutt 1998) modified to assume a Chabrier (2003) IMF:

$$\begin{eqnarray}&&{L}_{{\rm{UV}}}=\left(\displaystyle \frac{{\rm{SFR}}}{{M}_{\odot }{{\rm{yr}}}^{-1}}\right)1.4\times {10}^{28}\,{\rm{erg}}\,{{\rm{s}}}^{-1}\,{{\rm{Hz}}}^{-1}\end{eqnarray} \tag{ 15 }$$

This relationship assumes a constant SFR for 100 Myr. We also apply these dust corrections to the Reddy & Steidel (2009) and McLure et al. (2013) LF results. Our quantitative results for the corrected and uncorrected SFR densities at z ∼ 3–10 are presented in Table 8.

Table 8. Star Formation Rate Densities Inferred to −17.0 AB mag (0.03 ${L}_{z=3}^{* }$)

Lyman		${{\rm{log}}}_{10}{ \mathcal L }$	Dust		${{\rm{log}}}_{10}$ SFR Density
Break		(erg s−1	Correction	(M⊙ Mpc−3 yr−1)
Sample	$\langle z\rangle$	Hz−1 Mpc−3) a	(dex) b	Uncorrected	Corrected	Incl. ULIRG b
U	3.0	26.55 ± 0.06	0.44	−1.60 ± 0.03	−1.26 ± 0.09	−1.16 ± 0.09
B	3.8	26.52 ± 0.06	0.39	−1.63 ± 0.06	−1.32 ± 0.06	−1.24 ± 0.06
V	4.9	26.30 ± 0.06	0.32	−1.85 ± 0.06	−1.58 ± 0.06	−1.53 ± 0.06
i	5.9	26.10 ± 0.06	0.20	−2.05 ± 0.06	−1.88 ± 0.06	−1.85 ± 0.06
z	6.8	25.98 ± 0.06	0.07	−2.17 ± 0.06	−2.10 ± 0.06	−2.10 ± 0.06
Y	7.9	25.67 ± 0.06	0.06	−2.48 ± 0.06	−2.42 ± 0.06	−2.42 ± 0.06
J	10.4	${24.62}_{-0.45}^{+0.36}$	0.00	$-{3.28}_{-0.45}^{+0.36}$	$-{3.28}_{-0.45}^{+0.36}$	$-{3.28}_{-0.45}^{+0.36}$

Notes.

^a Integrated down to 0.03 ${L}_{z=3}^{* }$. Based upon LF parameters in Table 2, Section 6, of Bouwens et al. (2015). The SFR density estimates assume $\gtrsim 100\,\mathrm{Myr}$ constant SFR and a Chabrier IMF (e.g., Madau & Dickinson 2014). Conversion to a Salpeter (1955) IMF would result in a factor of ∼1.6 (0.2 dex) increase in the SFR density estimates given here. ^b The contribution indicated here is our fiducial estimate from far-IR bright, ULIRG-like ($\gt {10}^{12}$ L_⊙) galaxies (see Figure 18). The SFR density contribution for the far-IR bright population tends to be either missed completely due to these sources not being selected in Lyman-break galaxy probes (e.g., Simpson et al. 2014) or significantly underestimated due to the IR luminosities underestimated based on their UV properties (e.g., Reddy & Steidel 2009).

Download table as:  ASCIITypeset image

In computing the SFR density, we must account not only for the impact of dust extinction on the UV luminosities themselves but also for the more massive, far-IR bright sources where standard dust corrections are not effective or which are sufficiently faint in the UV to be entirely missed in standard LBG searches (e.g., Reddy et al. 2006, 2008; Swinbank et al. 2014; Casey et al. 2018; Williams et al. 2019; Dudzevičiūtė et al. 2020). Such sources are known to contribute a substantial fraction of the SFR density at z ∼ 0–3 (Hughes et al. 1998; Blain et al. 1999; Lilly et al. 1999; Chapman et al. 2005; Karim et al. 2011; Barger et al. 2012; Magnelli et al. 2013; Madau & Dickinson 2014; Swinbank et al. 2014; Wang et al. 2019; Dudzevičiūtė et al. 2020). Perhaps the best way to account for these galaxies (proposed earlier by Reddy et al. 2008) is to simply include them based on dedicated searches for these sources in the IR.

We consider the results of Magnelli et al. (2013) at z ∼ 0–2 (which build on the results of Caputi et al. 2007 and Magnelli et al. 2009, 2011), the Franco et al. (2020a) results at z ∼ 2–5 from a 69 arcmin² survey area, the Yamaguchi et al. (2019) results at z ∼ 3–5 from the 26 arcmin² ASAGAO survey area, and the Williams et al. (2019) serendipitous discovery of a probable dusty star formation source at z ∼ 5. We compute the SFR density contribution from ultraluminous IR-type galaxies (ULIRGs) at $z\sim 2.5$ from the ASPECS volume by converting the measured ALMA fluxes from detected sources in their survey area to SFRs assuming 100% of the energy comes from star formation, binning the contributions by the derived redshifts for the sources, and then dividing by the cosmic volume within a 4.2 arcmin² survey area, finding 0.036 ± 0.022 M_⊙ yr⁻¹. We use a similar approach to derive the SFR density contribution from ULIRGs at z ∼ 2–5 from the 69 arcmin² Franco et al. (2020a) probe, but given the limited depth of this probe, we treat this derived contribution as a lower limit. The Franco et al. (2020a) probe builds on the earlier Franco et al. (2018) study using deeper search results presented in Franco et al. (2020b).

Additionally, we consider the integrated SFR density derived from MAGPHYS fits to the ∼1 deg² AS2UDS sample by Dudzevičiūtė et al. (2020), who corrected for incompleteness using the number counts from Geach et al. (2017) and extrapolated from the observed 870 μm flux limit of 3.6 mJy for the SCUBA-2 survey to 1 mJy using the slope of the number counts from Hatsukade et al. (2018). See Section 5.4 and Figure 15 from Dudzevičiūtė et al. (2020) for further details. The uncertainties on the values were calculated by resampling the SFR and redshift probability distributions of each source. The approximate survey volume for the AS2UDS results is approximately 7 × 10⁷ Mpc³ and so is likely to be much more representative than smaller volume studies.

Finally, we also estimate the SFR density from ULIRG-type sources by assuming that the star-forming main sequence results of Speagle et al. (2014) apply to the wide-area z ∼ 1.5–3 mass functions of Ilbert et al. (2013) and z ∼ 3–6 mass functions of Davidzon et al. (2017). Encouragingly enough, the estimated SFR contribution provided by ULIRG-type galaxies using the observed mass functions appears to be plausibly consistent with that derived from constraints available from direct searches for ULIRG-type sources at z = 2–4 and ∼0.2 dex higher at z = 4–6.

A summary of the inferred SFR density for all the aforementioned ULIRG probes is presented in Figure 18. As our fiducial estimate of the obscured SFR density from ULIRGs, we adopt the Magnelli et al. (2013) constraints at z ∼ 2, our mass-function-derived estimate at $z\sim 2.75$, the AS2UDS estimates (Dudzevičiūtė et al. 2020) at z = 3.4–6, and the Wang et al. (2019) at $z\gt 6$. We have indicated fiducial obscured SFR densities in Figure 18 with the hatched red area. This fiducial model is most consistent with the dust-poor model from Casey et al. (2018).

Zoom In Zoom Out Reset image size

We combine these SFR densities with those we derived by correcting the UV LFs at z = 3–10 to present our best estimates for the SFR density at z = 3–10 in Table 8 and Figure 19, together with a few previous estimates (Schiminovich et al. 2005; Reddy & Steidel 2009; McLure et al. 2013) in Figure 19. It is interesting to compare the contribution that unobscured and obscured star formation makes to the total SFR density of the universe. Figure 20 shows such a breakdown of the SFR density. The obscured SFR density shown at $z\lt 2$ is from the Magnelli et al. (2009, 2011, 2013), while at z ∼ 2–3, the obscured SFR density shown is the sum of the SFR density from AS2UDS (Dudzevičiūtė et al. 2020) and from Reddy & Steidel (2009). The contribution to the SFR density from ULIRGs is presented for context. From the presented breakdown, we can see that star formation is mostly unobscured at $z\gt 5$, mostly obscured at $z\lt 5$, and z ∼ 5 marks the approximate transition redshift between the two regimes. Previously, Bouwens et al. (2009), Bouwens et al. (2016), and Dunlop et al. (2017) found that the approximate transition point between the two regimes was z ∼ 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, before closing the discussion we provide in this paper on the SFR density, it is interesting to try to estimate the SFR density within the ASPECS HUDF/XDF volume itself. Given the limited volume probed by ASPECS and the impact of large-scale structure, this is an interesting issue to examine to help determine the extent to which conclusions drawn from the HUDF volume are applicable to much larger cosmic (∼10⁶ comoving Mpc³) volumes of the universe (where the impact of large-scale structure is less).

To estimate the approximate SFR density within the ASPECS volume, we rederive the UV LF at z ∼ 2, z ∼ 3, z ∼ 4, z ∼ 5, z ∼ 6, z ∼ 7, z ∼ 8, and z ∼ 10 but only using sources in the ASPECS/HUDF/XDF volume. For simplicity, in deriving this LF, we fix the faint-end slope α and characteristic luminosity M* to that derived from Bouwens et al. (2015) and R. Bouwens et al. (2020, in preparation) at these same redshifts and fit for the normalization ${\phi }^{* }$. The relative normalization we derive for the UV LFs over the ASPECS areas relative to the cosmic average, i.e., ${\phi }_{\mathrm{ASPECS}}^{* }/\langle {\phi }^{* }\rangle$ is 1.17, 0.85, 0.68, 0.74, 1.10, 0.95, 0.98, 1.29, and 0.50 at z ∼ 2, z ∼ 3, z ∼ 4, z ∼ 5, z ∼ 6, z ∼ 7, z ∼ 8, z ∼ 9, and z ∼ 10, respectively. The rms logarithmic scatter in these normalizations are 0.12 dex, i.e., fluctuations of 32% (0.12 dex) in the volume density of galaxies in a given redshift interval of the HUDF relative to the cosmic average. 32% is fairly similar to the expected variations one would expect for sources with a volume density of ∼1 × 10⁻³ Mpc^–3 inside a 2' × 2' × Δz ∼ 1 volume. Figure 21 illustrates how the SFR density we infer from the HUDF might compare with the cosmic average, if we assume that we can apply the LF normalization factors just derived to the SFR density as a whole.

Zoom In Zoom Out Reset image size

To assess the impact of large-scale structure on the present results and other results from ASPECS, it is relevant to compare the observed 0.12 dex scatter with that expected from the relatively small number of dust-detected and CO-detected sources over ASPECS. In cases where the number of sources per unit redshift is in the range 10–15, i.e., similar to the number of dust-detected and CO-detected sources in ASPECS (e.g., Boogaard et al. 2019; González-López et al. 2020), the scatter expected from small number statistics will be comparable to that seen in terms of large-scale structure. This suggests that any conclusions drawn from the HUDF ASPECS volume should be applicable to much larger cosmic (∼10⁶ comoving Mpc³) volumes, with a relatively limited impact from large-scale structure.