The Physical Drivers of the Luminosity-weighted Dust Temperatures in High-redshift Galaxies

Seven unlensed, spectroscopically confirmed DSFGs were chosen from 400 arcmin² SCUBA-2 450 and 850 μm maps of the inner COSMOS field (Casey et al. 2013). Of the 31 SCUBA-2-detected COSMOS sources with spectroscopic redshifts in Casey et al. (2017), we selected seven of the brightest sources (single-dish S₈₅₀ > 2 mJy) to span a very broad range of dust SEDs from cold (∼18 K) to warm (∼70 K) based on their Herschel SPIRE and SCUBA-2 photometry. The DSFGs were selected at the time for having known spectroscopic redshifts of z = 1 to 3 (Casey et al. 2017), and updated information on this sample gives a total spanned redshift range of 1.4 < z < 4.6. The galaxy initially thought to have the coldest SED, 450.03 (also known as AzTEC2), was identified initially to have a spectroscopic redshift of 1.123 from an optical/infrared (OIR) counterpart (Casey et al. 2017) but is now confirmed to sit at z = 4.61 from a serendipitous detection of [C ii] in our ALMA data and concurrent confirmation via CO (5–4) line emission from NOEMA observations (Jiménez-Andrade et al. 2020). The ALMA centroid is offset from the OIR source by 16 (see Casey et al. 2017 for details). The redshift for source 850.04 was originally estimated to be z = 1.436, but closer inspection reveals a more likely redshift solution of $z={3.31}_{-0.81}^{+0.76}$, which is a photometric redshift from the Laigle et al. (2016) COSMOS catalog. This photometric counterpart is 10 offset from the spectroscopically confirmed source at z = 1.436 (Casey et al. 2017) but significantly closer (03) from the centroid of ALMA emission.

Observations of these seven DSFGs were carried out as part of the ALMA Cycle 3 program 2015.1.00568.S (PI: Casey). With a requested spatial resolution of 01, observations were made in both compact (C36-3) and extended (C36-6) configurations. The compact and extended configurations had maximum baselines of 462.9 m and 1.8 km, respectively, with 36 antennas. The observations targeted dust continuum emission at the nominal band 7 frequency centered at 345 GHz or 870 μm. At this frequency, ALMA has a 173 primary beam (FWHM), and we utilized the "single continuum" spectral mode.

We reduced and calibrated the raw data from the compact configuration to produce u-v data products using the Common Astronomy Software Applications (casa) version 4.5.3. During calibration, data from antennas with irregular amplitudes were flagged. To produce optimum images for total flux recovery, we imaged the compact-configuration data using the clean algorithm of casa with Briggs weighting and robust = 2 (equivalent to natural weighting) and applied a primary beam correction. The resulting images have a resolution of 068 × 049. The compact data set was used in calculating the flux densities, which are given in Table 1.

Table 1. Redshifts, Total Flux Densities, Half-light Radii as Measured by imfit at 870 μm, Rest-frame SED Peak Wavelengths, Infrared Luminosities, Star Formation Rates, and Stellar Masses for Our COSMOS SCUBA-2 and ALESS Sources

Name	z	S870	Re	λpeak	LIR	SFR	M⋆
		(mJy)	( kpc)	(μm)	(L⊙)	(M⊙ yr−1)	(M⊙)
450.03	4.629 a	19.12 ± 0.16	1.35 b ± 0.06	71.4 ± 7.8	(3.28 ± 0.68) × 1013	4900 ± 1000	⋯ c
850.04	3.31	13.67 ± 0.12	1.94 ± 0.12	98 ± 9	(8.9 ± 1.2) × 1012	1330 ± 180	(2.2 ± 0.8) × 1010
450.14	1.523	2.05 ± 0.03	2.8 ± 0.4	130 ± 30	(1.2 ± 0.7) × 1012	190 ± 100	⋯ c
850.60	1.457	2.21 ± 0.04	1.81 ± 0.28	75 ± 9	(1.9 ± 0.4) × 1012	290 ± 60	(2.1${}_{-0.2}^{+0.1}$) × 1011
850.95	1.556	1.21 ± 0.03	6.1 ± 1.4	98 ± 22	(1.2 ± 0.5) × 1012	180 ± 80	(3.8${}_{-1.0}^{+0.4}$) × 1010
850.25	2.438	5.76 ± 0.05	2.76 ± 0.22	117 ± 20	(3.1 ± 1.1) × 1012	460 ± 170	(2.8${}_{-0.4}^{+0.1}$) × 1011
450.09	2.472	6.07 ± 0.07	1.04 ± 0.07	71 ± 10	(8.3 ± 1.5) × 1012	1240 ± 230	(2.2${}_{-0.4}^{+0.1}$) × 1011
ALESS 3.1	4.24	8.3 ± 0.4	1.30 ± 0.05	69 ± 9	(1.64 ± 0.25) × 1013	2460 ± 380	(2.3${}_{-1.6}^{+1.0}$) × 1011
ALESS 5.1	2.86	7.8 ± 0.7	2.00 ± 0.10	106 ± 12	(4.9 ± 1.0) × 1012	730 ± 150	(2.9${}_{-0.8}^{+0.2}$) × 1011
ALESS 9.1	4.50	8.8 ± 0.5	1.50 ± 0.05	71 ± 9	(1.63 ± 0.26) × 1013	2430 ± 390	(5.6${}_{-2.7}^{+1.3}$) × 1011
ALESS 10.1	0.76	5.3 ± 0.5	2.60 ± 0.20	190 ± 9	(2.3 ± 0.4) × 1011	34.5 ± 5.5	(1.2 ± 0.3) × 1010
ALESS 15.1	1.93	9.0 ± 0.4	2.40 ± 0.10	125 ± 5	(2.95 ± 0.24) × 1012	440 ± 40	(3.1${}_{-1.8}^{+1.3}$) × 1011
ALESS 17.1	1.54	8.4 ± 0.5	1.75 ± 0.05	138 ± 5	(1.88 ± 0.18) × 1012	281 ± 27	(2.3${}_{-1.7}^{+0.1}$) × 1011
ALESS 29.1	1.44	5.9 ± 0.4	1.35 ± 0.05	136 ± 12	(1.37 ± 0.23) × 1012	200 ± 30	(3.0${}_{-2.1}^{+1.2}$) × 1011
ALESS 39.1	2.44	4.3 ± 0.3	1.95 ± 0.10	100 ± 7	(3.3 ± 0.4) × 1012	490 ± 70	(4.5${}_{-1.6}^{+1.0}$) × 1010
ALESS 45.1	2.34	6.0 ± 0.5	2.15 ± 0.15	117 ± 10	(2.7 ± 0.5) × 1012	400 ± 70	(5.1${}_{-2.3}^{+1.5}$) × 1011
ALESS 67.1	2.12	4.5 ± 0.4	1.85 ± 0.15	72.2 ± 2.6	(9.4 ± 0.9) × 1012	1400 ± 140	(2.4${}_{-2.1}^{+2.4}$) × 1011
ALESS 112.1	2.32	7.6 ± 0.5	1.90 ± 0.10	109 ± 6	(4.1 ± 0.5) × 1012	610 ± 70	(1.0${}_{-0.5}^{+0.4}$) × 1011

Notes. Measurements of peak wavelengths, infrared luminosities, and SFRs are done by the same method for both the COSMOS sample and the 11 ALESS sources we analyze. Measurements for 450.03 are for the total circularized size and total L_IR of both components. ALESS flux densities are from Swinbank et al. (2014), and sizes are from Hodge et al. (2016), derived in a fully consistent manner as for our COSMOS sample.

^a The precision quoted on the redshift indicates whether or not the source has a spectroscopic redshift (four significant figures) or a photometric redshift (three significant figures). ^b Total circularized size of both components of 450.03. The individual half-light radii measured by imfit are 1.24 kpc for the brightest source and 1.08 kpc for the secondary source. ^c Stellar masses are not provided for 450.03 or 450.14. In the case of 450.03, this is due to significant blending with several foreground galaxies on angular scales that precludes measurement of near-infrared Spitzer flux densities; in the case of 450.14, there is a superposition of two galaxies of different redshifts within a 2'' region, making a proper separation of the background/foreground stellar mass distribution impossible. Both cases are discussed at greater length in Casey et al. (2017).

Download table as:  ASCIITypeset image

We created a higher-resolution set of maps by combining the raw data from the compact and extended configurations and then reducing and calibrating the data using the same method described above. We tested different robust values with the clean algorithm and found the resulting size of each source to be consistent and independent of the robust parameter. Thus, we adopt a robust value of 1 for the full sample, as it maximizes the source signal-to-noise ratio for the most extended source in the sample (850.95) while maintaining a relatively high angular resolution. The resulting resolution is 024 × 027. This does not meet the requested 012 resolution, but even with a robust value of 0, the spatial resolution is still not as good as requested. The combined data were used to calculate the effective radii, which are also given in Table 1. Images of the sample are shown in Figure 1, and overlays with optical/near-infrared data are shown in Figure 2 on the same scale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The data sample used in this analysis does not fully sample the full parent population of galaxy dust SEDs but is representative of the breadth of DSFGs at high z. The COSMOS sources from Casey et al. (2017) were chosen explicitly to have the broadest range of SEDs and are selected at either 450 μm or 850 μm, which together is unbiased with respect to SED dust temperatures in the range of 20–100 K (e.g., Casey et al. 2013; Roseboom et al. 2013). The ALESS sample selected at 870 μm only is more likely to be biased toward colder dust temperatures (Chapman et al. 2004a; Casey et al. 2009), but we see a large range of temperatures well represented in the sample, similar to the COSMOS sample. Furthermore, while 870 μm selection alone is intrinsically biased and favoring sources with colder SEDs, those that have well-sampled FIR SEDs (and high-resolution ALMA imaging) are not biased with respect to λ_peak because detection in the Herschel bands bracketing the SED peak is implicitly required. In other words, many ALESS sources without Herschel SPIRE detections (i.e., those with cold SEDs) would not have well-sampled FIR SEDs, are unlikely to have received high-resolution ALMA follow-up, and thus are excluded from our analysis here.

The sample (both COSMOS and ALESS DSFGs) could be seen as biased toward galaxies with high SFRs, but sources with high SFRs are the only ones for which the high spatial resolution measurements made in this paper are possible with a reasonable investment in ALMA time. This has potential impact on our ability to probe the relationship between λ_peak and L_IR in a fully unbiased way, given the limited dynamic range compared to galaxy populations at large. This is discussed more fully in Section 5.2.2. The combined sample does sample galaxies' relative distance from the main sequence well, in line with the expected dynamic range of much larger populations of studied galaxies. They span two orders of magnitude in sSFR from ∼0.6 to 60 Gyr⁻¹, and about one order of magnitude in orthogonal distance from the main sequence; more methods to quantify galaxies' distance to the main sequence are discussed later in Section 4.4.

The primary goal of this work is to assess the relative strength of correlation between luminosity-weighted dust temperature (as measured observationally through the proxy λ_peak) and other physical characteristics of high-z galaxies, including quantities constrained by the FIR-emitting region as measured by high-resolution ALMA dust continuum imaging. Thus, this set of observations constitutes a well-constructed and relatively unbiased sample for our purposes.

Total integrated flux densities were measured directly from the compact ALMA data using a method similar to the one described in Hodge et al. (2016). We converted the maps from Jy beam⁻¹ to mJy pixel⁻¹ and masked emission less than 2σ. We then summed contiguous emission within an aperture of radius 3 × b_maj, where b_maj is the FWHM (major axis) of the synthesized beam. The resulting total band 7 flux densities are given in Table 1. The median flux density of our seven sources is 5.8 ± 3.7 mJy at 870 μm.

In general, when measuring the flux density, only contiguous emission was summed, with one exception. Source 450.03 is composed of two sources separated by 3'' at the same redshift (Jiménez-Andrade et al. 2020), both of which were included in the total reported flux density.

To check our flux density calculations, we used the imfit tool in casa and found that our calculations agree with the imfit-derived densities for all but one of our sources. Source 850.95 differed with less than 2σ significance of deviation. The average ratio between imfit and aperture flux density measurements was 0.96 ± 0.10.

In addition, we compare the flux densities as measured by ALMA to the single-dish flux densities obtained for the COSMOS sources from SCUBA-2 at 850 μm using the single-dish JCMT and for the ALESS sources from LABOCA at 870 μm. The average ratio of ALMA to single-dish flux densities was 0.94 ± 0.21, as shown in Figure 3. Note that the expected ratio of 850 μm to 870 μm flux is ≈0.9 for RJ tail emission.

Zoom In Zoom Out Reset image size

We quantitatively examine the size of the 870 μm dust emission by fitting each source in the image plane with two-dimensional Sérsic and Gaussian models. The Sérsic model was applied using galfit, and the Gaussian model was applied using casa imfit.

Dust emission half-light effective radii (R_e) were measured using galfit on the ALMA maps in the image plane. galfit is a 2D fitting tool that models the profiles of astronomical objects using parametric functions (Peng et al. 2010). Here we use galfit on ALMA data to test for non-Gaussian source morphologies, which is most easily applied in the image plane. The full width at half maximum (FWHM) was measured for each source using casa imfit, which is also applied in the image plane. Most of our sources are moderately resolved—more than two beams across in size—so size measurements are unlikely to be affected by fitting methods. Thus, making measurements in the image plane rather than the u-v plane is justified.

To address whether the weighting applied with the clean algorithm has an impact on size measurement in the image plane, measurements were also made in the u-v plane using casa uvmodelfit to fit a Gaussian model. The uvmodelfit and imfit measurements agree within uncertainties for four of our sources. There is a discrepancy between u-v and image-plane sizes for 450.03 and 850.25, which we attribute to significant excess emission from the core of the systems, which is evident in the imfit residual images. There was also a discrepancy in the measurements for source 850.95, which we attribute to the imfit model extending beyond the region of significant flux associated with the source. In general, the measured Gaussian FWHMs in the u-v plane are systematically offset from the FWHMs in the image plane measured by imfit such that the median ratio of u-v-plane to image-plane FWHMs is 0.81 ± 0.08. Our scientific results rely on relative size comparisons of objects within the sample (and extended to the ALESS sample) rather than absolute sizes. The results of galfit are shown in Figure 4, along with the residuals and models. It is a known problem that errors from galfit are underestimated; therefore, we used the imfit results in our analysis and report imfit sizes in Table 1.

Zoom In Zoom Out Reset image size

The sizes measured here are the mass-weighted sizes—measured on the RJ tail of the SED—rather than the luminosity-weighted sizes, which we would measure closer to the peak of the SED. The only exceptions to this are those galaxies that sit at z > 3 in our sample (5 of 18); at these redshifts, observed 870 μm probes rest-frame wavelengths shortward of 200 μm, which may be more prone to tracing an optically thick regime of the SED in some cases (dependent on the actual dust column). We explore the impact of these higher-redshift sources on our final results later in Section 5. A dearth of data in the literature makes it difficult to know whether measured size varies substantially along the RJ tail, but the limited one-off measurements that do exist suggest that it is safe to assume that any measurement on the tail should be consistent within measurement uncertainties. In addition, these sizes are distinct from the sizes of the galaxies traced by stellar emission; the two are often spatially distinct in DSFGs, such that the rest-frame UV/optical emission is heavily obscured by dust and thus does not serve as an adequate probe of the characteristic scale of the ISM (Biggs & Ivison 2008; Chen et al. 2015; Hodge et al. 2015, 2016; Lang et al. 2019). It should also be noted that dust continuum sizes are systematically smaller and not entirely consistent with radio continuum sizes (which probe synchrotron emission scattering off of supernova remnants), the measurements of which exist for much larger samples than those with dust continuum size measurements (Miettinen et al. 2017).

The effective radii measured by galfit range from 013 to 069 with a median of R_e = 032 ± 009. Using the redshift of each source, we convert these angular sizes to physical sizes, which range from 1.1 to 5.9 kpc (median R_e = 2.7 ± 0.8 kpc). Source 850.95 was an outlier at 5.9 ± 3.3 kpc. It is known to be somewhat abnormally large as also traced by Hα kinematics (Drew et al. 2018). Source 850.95 also has the lowest signal-to-noise ratio among our sources. The remaining sources ranged from 1.1 to 3.2 kpc (median R_e = 2.3 ± 0.5 kpc).

The FWHMs (deconvolved from the restoring beam) measured by imfit are shown in Table 1. They range from 026 to 145, which convert to 2.1−12.3 kpc on physical scales. The half-light radii quoted in Table 1 are equivalent to these FWHM measurements divided by two. The median Gaussian profile has an FWHM = 3.9 ± 1.6 kpc. Excluding 850.95, the median is 3.8 ± 1.4 kpc.

The Gaussian profiles of the ALESS sources included in our analysis were measured via the same technique used here (Hodge et al. 2016), where the median FWHM = 3.8 ± 0.5 kpc. The median size of our sources agrees with that of the ALESS sources, despite the fact that the ALESS sample is brighter than all but two of our sources. Overall, the DSFG populations from which the ALESS and COSMOS SCUBA-2 DSFGs are drawn are very similar, with the exception of the ALESS median total 870 μm flux densities being 1.3× brighter. The median of our sources' half-light radii is also consistent within uncertainty with that found by Spilker et al. (2016) and Gullberg et al. (2019).

The Sérsic indices of our sources range from n = 0.3 to 2.1, with a median of n = 0.8 ± 0.4. Our median n is consistent with that of the ALESS sources in Hodge et al. (2016), whose median Sérsic index n = 0.9 ± 0.2 implies a non-Gaussian, exponential disk morphology. Gullberg et al. (2019) also found an exponential disk emission profile for a sample of DSFGs.

Wien's law states that λ_peak = b/T_d (where b = 2.898 × 10⁻³ m K) for an idealized blackbody, but this relationship deviates depending on the adopted opacity model of the blackbody as shown in Figure 20 of Casey et al. (2014). Thus, we parameterize the luminosity-weighted dust temperature (T_d) instead as the observable λ_peak, the rest-frame wavelength where the dust SED peaks. The observable λ_peak is not sensitive to opacity assumptions and is a better trace of the observed peak in the SED in any observational sample without high-quality constraints on the underlying dust opacity.

SEDs in this work were fitted using the technique described in Casey (2012), based on a single dust temperature modified blackbody fit with a mid-infrared power law shortward of ∼70 μm in the rest frame. The mid-infrared power law represents a superposition of modified blackbodies at warm to hot temperatures, where the least amount of dust mass is heated to the hottest temperatures, and the vast majority of the mass of dust in the ISM is at much colder temperatures, represented by the temperature of the fit at rest-frame wavelengths ≳200 μm. Because our sample has somewhat limited rest-frame mid-infrared photometric constraints, we fix the mid-infrared power-law slope to α_MIR = 2, representative of the average measured power-law slope for local LIRGs and ULIRGs that have well-sampled mid-IR SEDs (e.g., U et al. 2012).

The photometric points used in the fitting include ALMA (870 μm), Herschel SPIRE (250, 350, 500 μm), SCUBA-2 (450, 850 μm), and Spitzer (24 μm). Points are weighted relative to their signal-to-noise ratio, with the exception of the Spitzer 24 μm point, whose signal-to-noise ratio is capped at 10 for the purposes of these bulk SED fits. The 24 μm photometry often is a very high signal-to-noise ratio measurement but is not always in agreement with the underlying SED model for dust continuum emission alone, due to the possible contribution of PAH emission features and silicate absorption in the rest-frame mid-infrared but not directly modeled here. We also note that use of the response function of the various FIR/millimeter passbands is not necessary using the Casey (2012) fitting method given the lack of ≈1% precision in the photometry in the FIR regime (see also Casey 2020).

While the bulk SED, including the measured IR luminosity or rest-frame peak wavelength, is not impacted by the adopted opacity model, here we adopt a fixed opacity model for all SEDs such that τ = 1 at 200 μm following Conley et al. (2011). Given that our observations include spatially resolved dust maps, we can test whether this assumption is appropriate for this sample by inferring the wavelengths at which τ = 1 from the dust column density. Using the total 870 μm flux densities from Section 4.1, sizes from Section 4.2, a fixed mass-weighted dust temperature of 25 K, and dust mass absorption coefficients from Li & Draine (2001), we infer that the wavelength at which the SEDs become optically thick varies from 10 to 300 μm (with an average λ(τ = 1) = 80 ± 50 μm). Given the uncertainties in the derivation of dust mass itself, we choose to apply the uniform assumption to our SEDs rather than tuning each fit according to its measured dust column. The SEDs are shown in Figure 5, and the measured peak wavelengths and IR luminosities are also given in Table 1. We also refit SEDs for the ALESS sources described in Hodge et al. (2016), whose multiwavelength photometry is given in Swinbank et al. (2014). Fitting SEDs for both the COSMOS SCUBA-2 DSFGs and the ALESS sources was done in a fully consistent manner.

Zoom In Zoom Out Reset image size

SFRs were calculated using the total infrared luminosities derived from the SED fits in Section 4.3 and the logarithmic constant of conversion, $\mathrm{log}({C}_{\mathrm{IR}})=43.41$ from Kennicutt & Evans (2012). We do not have reliable stellar mass estimates for 450.03 and 450.14 owing to foreground−background source contamination (see Casey et al. 2017), so they are omitted from all analysis requiring measurements of M_⋆. The stellar masses of the COSMOS sample are presented in Casey et al. (2017), and for the ALESS sources they are presented in da Cunha et al. (2015). All stellar masses used in this study are derived with MAGPHYS (da Cunha et al. 2008, 2015) in a uniform fashion using stellar population models from Bruzual & Charlot (2003). For the two sources for which the redshifts were corrected, we reran the MAGPHYS modeling accordingly. Note that the stellar masses of DSFGs are often highly uncertain, though not necessarily systematically so, as they can also be overestimated by applying an inappropriate star formation history (see discussions in Hainline et al. 2011; Michałowski et al. 2012). The uncertainties on the stellar masses are consistent with the broad uncertainties on DSFG stellar masses in general. We note that our MAGPHYS fits are used exclusively to derive stellar masses, and the MAGPHYS fits are not used in this work to infer FIR SED characteristics like dust temperature or SFR; this is because of the potential for the rest-frame optical/near-infrared SED to be highly decoupled from FIR emission in DSFGs, as well as our desire to directly test the FIR fits alone and independently of the rest-frame optical/near-infrared emission used to measure the stellar masses. We also calculated the sSFR, where sSFR = SFR/M_⋆.

Figure 6 shows the COSMOS and ALESS sources in relation to the SFR–M_⋆ relation, also called the galaxy main sequence. In red, we show the main sequence of star-forming galaxies with a quadratic function of the form used in Whitaker et al. (2014) applied to COSMOS data (1.5 ≤ z ≤ 2.5) from Laigle et al. (2016). We applied a Markov Chain Monte Carlo (MCMC) approach to derive the parameters of the fit for the broader COSMOS data set. Also shown are the cosmic time-dependent best fits from Speagle et al. (2014) for redshifts z = 1, 2, 3, 4, and 5.

Zoom In Zoom Out Reset image size

For each galaxy in our sample, we also compute a parameterized "distance" to the main sequence, or D_MS. D_MS represents the orthogonal distance between the main sequence at the source's redshift and its measured value in the $\mathrm{log}{M}_{\star }\mbox{--}\mathrm{log}\,\mathrm{SFR}$ plane. We compute orthogonal distances to the main sequence rather than projected distances in SFR or sSFR to minimize the impact of measurement uncertainties in both M_⋆ and SFR for each galaxy. The main sequence used to calculate D_MS for each source were the cosmic time-dependent curves given in Speagle et al. (2014); however, we note that computing D_MS using a quadratic form of the main sequence as in Whitaker et al. (2014) produces similar results, ultimately not changing the conclusions of our work no matter which definition of D_MS is adopted. Similarly, as we discuss later in Section 5, we test for correlation between λ_peak and ΔsSFR, finding results that are consistent with our adopted D_MS.

Most galaxies for which dust temperatures have been constrained lack spatially resolved measurements needed to investigate the effects of dust geometry on dust temperature. Thus, the underlying physical driver of a galaxy's globally averaged dust temperature has not been directly constrained. Before ALMA, there were small samples of high-z galaxies that had resolved radio or millimeter sizes (Chapman et al. 2004b; Biggs & Ivison 2008; Casey et al. 2009; Younger et al. 2009), but the observations needed to determine accurate sizes for larger samples of DSFGs require the order-of-magnitude improvements in sensitivity brought by ALMA (see Hodge et al. 2016; Rujopakarn et al. 2016; Simpson et al. 2017; Elbaz et al. 2018; Gullberg et al. 2019). In this paper, we have presented a sample of DSFGs with both high-resolution dust size constraints and good dust SED constraints, allowing a more thorough investigation into the physical drivers of galaxies' dust SEDs.

This study is theoretically motivated by the application of the Stefan–Boltzmann law on galaxy scales. Specifically, it relates the emergent IR luminosity to the effective size of the blackbody (R_e) and temperature (T_d) of an optically thick spherical blackbody:

$$\begin{eqnarray}&&{L}_{\mathrm{IR}}=4\pi {R}_{{\rm{e}}}^{2}\sigma {T}^{4}.\end{eqnarray} \tag{ 1 }$$

The Stefan–Boltzmann constant is σ = 5.670 × 10⁻⁵ erg s⁻¹ cm⁻² K⁻⁴. Galaxies are certainly not as simple as stars in the application of such a relation. Nevertheless, approximating galaxies' IR emission as an optically thick shell of dust surrounding an obscured point source (i.e., deeply embedded OB associations in star-forming regions) would result in the following expected relation between the dust's temperature, IR luminosity, and dust shell radius:

$$\begin{eqnarray}&&\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{IR}}=\mathrm{log}(4\sigma )+4\mathrm{log}({T}_{\mathrm{dust}}),\end{eqnarray} \tag{ 2 }$$

where Σ_IR is the IR luminosity surface density (=${L}_{\mathrm{IR}}/\pi {R}_{{\rm{e}}}^{2}$). Of course, galaxies' geometry is unlikely to be well represented as a spherical shell, and the surface area of the emergent IR flux does indeed play a role in the expected relation. If instead of a spherical shell we assume some uniform planar distribution of UV emission underneath an optically thick screen (as in a galaxy disk, with half-light radius R_e ), we would expect the following relationship instead:

$$\begin{eqnarray}&&\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{IR}}=\mathrm{log}(\sigma )+4\mathrm{log}({T}_{\mathrm{dust}}),\end{eqnarray} \tag{ 3 }$$

which is a factor of four different than the spherical shell: one factor of two given the modified surface area of a disk as $2\pi {R}_{{\rm{e}}}^{2}$, and one factor of two accounting for R_e ² being a measured half-light radius. Then, given our use of optically thick SEDs (described in Section 4.3), the relationship between dust temperature and λ_peak can be approximated as

$$\begin{eqnarray}&&\mathrm{log}({T}_{\mathrm{dust}}/{\rm{K}})\approx 3.756-1.048\,\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}})\end{eqnarray} \tag{ 4 }$$

in the 15–70 K range. ⁹ This gives a predicted relationship between λ_peak and Σ_IR of

$$\begin{eqnarray}&&\mathrm{log}({{\rm{\Sigma }}}_{\mathrm{IR}}/[{L}_{\odot }\,{\mathrm{kpc}}^{-2}])\approx 20.18-4.19\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}}).\end{eqnarray} \tag{ 5 }$$

Note that the expectation using Wien's law directly instead of the approximation given in Equation (4) would yield

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({{\rm{\Sigma }}}_{\mathrm{IR}}/[{L}_{\odot }\,{\mathrm{kpc}}^{-2}]) & = & \mathrm{log}(\sigma {b}^{4})-4\,\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}})\\ & & \approx 19.00-4\,\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}}),\end{array}\end{eqnarray} \tag{ 6 }$$

where b = 2.898 × 10³ μm K.

One key caveat here is that the Stefan–Boltzmann law should only be applicable to optically thick blackbodies, and the dust in galaxies' ISM is often far from optically thick at all wavelengths (of course, it is a huge benefit that ISM dust is often optically thin on the RJ tail of blackbody emission, allowing for a direct measurement of galaxies' dust mass from measured flux density). We return to the conundrum of dust opacity again in Section 5.2.1, where we derive the relationship between Σ_IR and λ_peak empirically.

Here we investigate the relative strengths of correlations between λ_peak (our observational proxy for dust temperature) and each of the following physical properties of galaxies: sSFR, D_MS, SFR, and Σ_IR. All of these quantities have been argued to correlate tightly with galaxies' dust temperatures, for example, sSFR and D_MS (e.g., Magnelli et al. 2014), SFR or L_IR (e.g., Chapman et al. 2003b; Symeonidis et al. 2013), and Σ_IR (e.g., Chanial et al. 2007; Lutz et al. 2016). For example, Magnelli et al. (2014) report that T_d correlates more strongly with sSFR and parameterized distance to the main sequence (D_MS) than L_IR (they lacked direct measures of R_e and so could not investigate L_IR R_e ⁻²). In addition, they deduced that cold galaxies (≲25 K) sit on the main sequence and warm galaxies (≳30–80 K) lie above the main sequence (Magnelli et al. 2014).

Here we use this high-z sample with well-measured sizes and SEDs to directly constrain which of these physical correlations is tightest and may therefore be the most fundamental to the internal ISM. This may be of particular use to the community when only some of these quantities can be measured (e.g., galaxies lacking high-resolution dust continuum imaging). We calculate the significance of each correlation as the deviation of the model from no correlation, or a horizontal line in the given parameter space. The results are shown in Figure 7.

Zoom In Zoom Out Reset image size

We find that the λ_peak versus D_MS and λ_peak versus sSFR relationships correlate with 3σ and 7σ significance, respectively. In contrast, the λ_peak versus SFR and λ_peak versus Σ_IR relationships carry an 18σ and 24σ correlation significance, respectively. Below, we interpret the underlying physics of the correlations.

To investigate whether the extended redshift range of our sources affects our results, we recalculated the correlations based on a reduced sample with 1.44 ≤ z ≤ 2.86, removing the high-z and low-z extrema from the sample. At z < 3, the reduced sample's continuum observations are on the optically thin portion of the RJ tail of dust emission (i.e., rest-frame wavelengths longer than ∼200 μm), probing cold dust in the ISM. With the redshift-restricted sample, the λ_peak–Σ_IR and λ_peak–SFR correlation significances are reduced somewhat to 16σ and 18σ, respectively, while the correlations between λ_peak and D_MS and between λ_peak and sSFR increase to 4σ and 8σ, respectively. However, the end results are the same: the λ_peak correlations with Σ_IR and SFR are significantly stronger than those with sSFR and D_MS.

5.2.1.  λ_peak versus Σ_IR

We empirically measure the following relationship between λ_peak and Σ_IR, or the star formation surface density, using our sample:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({{\rm{\Sigma }}}_{\mathrm{IR}}/[{L}_{\odot }\,{\mathrm{kpc}}^{-2}]) & = & (18.78\pm 0.33)\\ & & -(3.80\pm 0.16)\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}}).\end{array}\end{eqnarray} \tag{ 7 }$$

The overall strength of the correlation is 24σ. The slope of this relation is a 2.4σ deviation from expectation in Equation (5), and the intercept deviates by 4.3σ. Comparing to the expectation using Wien's law directly (Equation (6)), our empirical relation only deviates by 1.3σ in the slope and 0.7σ in the intercept.

The agreement between Equations (5)–(6) and Equation (7) should be somewhat surprising, given that the geometry of galaxies' ISM is quite complex, and that the ISM is far from optically thick! How might we interpret the remarkably close correlation between λ_peak and Σ_IR?

Such a tight correlation would be expected if the peak of dust luminosity (traced directly by λ_peak) originates from regions of the ISM that are, indeed, optically thick. This may not be in direct conflict with more diffuse, cold dust in galaxies' ISM being optically thin, as that dust may not directly dominate the bulk of emission at the peak of IR emission.

We note that Simpson et al. (2017) provide an analysis of several (other) high-z DSFGs that have resolved ALMA sizes and measured SED dust temperatures; their Figure 5 shows the relationship between T_d and Σ_IR, suggesting tension between predictions from the Stefan–Boltzmann law and galaxies' dust temperatures, when assumed to be optically thick (they find more agreement when using the best-fit temperatures assuming optically thin SEDs). The difference in interpretation between our work and that of Simpson et al. is the presumed normalization of the Stefan–Boltzmann law itself and the geometry of dust: their application of the Stefan–Boltzmann law is drawn from a presumed spherical shell of dust surrounding a single star-forming region (as in Equation (2)), while our normalization accounts for a more disk-like configuration of a dust screen with several spatially distinct deeply embedded star-forming regions (as in Equation (3)). After accounting for the different geometries, we suggest that the Simpson et al. results, like ours, are quite consistent with a galaxy-scale Stefan–Boltzmann approximation for the emergent IR luminosity. Given the strength of the correlation in our data set and the plausible underlying physical interpretation of this relation, we determine that the λ_peak versus Σ_IR relation is the most fundamental of those we explore in this paper.

We note that the correlation strength of 24σ found for λ_peak–Σ_IR is unaffected by calculating Σ_IR with circularized radii or effective radii from Sérsic fitting rather than FWHMs from Gaussian fitting. The strengths of the correlations are also not affected by removing from our analysis the two sources (450.03 and 450.14) for which we do not have reliable stellar mass estimates from our analysis.

It is important to note that our measured sizes trace the dust mass as probed on the RJ tail of blackbody emission in what is expected to be an optically thin regime. This may be a different size than would be probed at the peak SED, which traces SFR density more directly. In other words, one might not expect the mass-weighted size to map directly to the SFR-weighted size. The relative agreement of the two is suggestive of a linear relationship between the effective radii over which the dust mass and star formation exist (i.e., R_e for SFR is similar to R_e for M_dust), which in turn may imply a linear slope for the Kennicutt–Schmidt law of SFR and gas surface density in these galaxies (Kennicutt 1998; de los Reyes & Kennicutt 2019). Further high-resolution follow-up work is needed to explore this in greater detail.

5.2.2.  λ_peak versus L_IR or SFR

We measure the following empirical relationship between L_IR and λ_peak in this sample:

$$\begin{eqnarray}&&\mathrm{log}({L}_{\mathrm{IR}}/{L}_{\odot })=(18.2\pm 0.3)-(2.81\pm 0.16)\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}}).\end{eqnarray} \tag{ 8 }$$

This correlation carries an overall significance of 18σ. The relationship with SFR can then be derived using the Kennicutt & Evans (2012) scaling. Figure 7 demonstrates that the trend broadly follows that of the star formation or IR luminosity surface density, Σ_IR.

We note that this relation is quite a bit steeper than found elsewhere in the literature for much larger samples of galaxies lacking high-resolution dust observations (e.g., Lee et al. 2013; Symeonidis et al. 2013; Kirkpatrick et al. 2017); for example, Casey et al. (2018b) modeled the L_IR– λ_peak relationship by a power law such that log₁₀( λ_peak) $=\,\mathrm{log}({\lambda }_{0})+\eta (\mathrm{log}({L}_{\mathrm{IR}})-12)$. Through empirical measurement of large samples of IR-luminous galaxies from 0 < z < 5, they found η = −0.068 ± 0.001 and $\mathrm{log}({\lambda }_{0}/\mu {\rm{m}})=2.012$. Though our findings here and the Casey et al. findings both show the same general trend—whereby higher-luminosity sources are hotter—here we find that the higher-luminosity sources are substantially hotter than those with fainter luminosities. As touched on earlier in Section 3.3, this is likely due to our selection of sources well-suited for high-resolution dust continuum follow-up (and therefore relatively high L_IR for a wide swath of λ_peak); therefore, it is likely that the more extensive samples from Casey et al. (2018a, and references therein) provide a more well-calibrated measurement of the λ_peak–L_IR relationship over a much larger dynamic range in L_IR.

Because the Casey et al. (2018a) model does not account for galaxy size as our analysis in Equation (7) does, we can use the two in conjunction to infer an overall size dependency of L_IR. Combining our Equation (7) with their λ_peak–L_IR relationship, we find that the size dependence can be modeled as a power-law function of L_IR such that R_e ∝ L_IR ^0.37±0.03. Note that this is roughly consistent with the size–luminosity relationship derived in Fujimoto et al. (2017), who found R_e ∝ L_IR ^0.28±0.07.

5.2.3.  λ_peak versus sSFR or Distance to the Main Sequence

Because both sSFR and the distance of a galaxy from the SFR–M_⋆ relation (the galaxy main sequence) both fundamentally trace physics dependent on both SFR and stellar mass, we discuss them here together. While we do find a subtle relationship between λ_peak versus sSFR and λ_peak versus D_MS, they are substantially weaker than direct correlations between λ_peak and Σ_IR or L_IR.

In contrast to some literature claims, we do not find clear evidence to support the theory that galaxies' dust temperatures correlate strongly with a galaxy's position on the galaxy main sequence. We note that the dynamic range of sSFR and D_MS in our sample is comparable to similar data sets at matched redshifts (∼2 dex in sSFR and ∼1 dex in D_MS; Barger et al. 2014). Specifically, we find the following correlations between λ_peak and sSFR or D_MS (where we have defined the latter as the orthogonal shortest distance to the Speagle et al. 2014 fits at the source's redshift):

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}(\mathrm{sSFR}/[{\mathrm{Gyr}}^{-1}]) & = & (3.4\pm 0.4)\\ & & -(1.43\pm 0.20)\,\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}})\end{array}\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&\mathrm{log}({D}_{\mathrm{MS}})=(1.5\pm 0.4)-(0.54\pm 0.21)\,\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}}),\end{eqnarray} \tag{ 10 }$$

which have overall correlation significance of 7σ and 3σ, respectively. Note that alternate definitions of D_MS (Δlog(SSFR)_MS and/or D_MS using the Whitaker et al. 2014 quadratic form of the main sequence) result in similarly low correlation strengths ranging from 1σ to 4σ.

The fact that both correlations between dust temperature and sSFR or D_MS are relatively weak suggests that the inclusion of stellar mass effectively dilutes the underlying physics driving the observed value of λ_peak. In fact, if we explicitly test for correlation between stellar mass and λ_peak, indeed, we find no significant correlation.

With integrated dust SED temperature lacking strong direct correlation with the main sequence, if the assertion is made that the main sequence informs on the evolutionary stage of a galaxy, then we infer that T_d (or λ_peak) is not a viable indicator of evolutionary class (starburst or secular disk) and that higher temperature does not necessarily indicate the presence of a starburst (Bothwell et al. 2010; Hodge et al. 2012; Drew et al. 2020). Rather, our results suggest that the underlying dust geometry within the ISM plays a much more substantial role. In addition, ISM dust can be heated by active galactic nuclei whether or not the galaxy hosts an ongoing starburst (Kirkpatrick et al. 2012). Additional recent cosmological simulation work supports the production of intense luminosities through secularly evolving disk systems from either gas infall, minor mergers, or disk instabilities (Davé et al. 2010; Hodge et al. 2012; Narayanan et al. 2015; Hayward et al. 2018; Tadaki et al. 2018).

Magnelli et al. (2014) arrived at different results than we have, suggesting significant correlation between D_MS and T_d. Instead of using our definition of D_MS, they evaluated D_MS to be Δlog(SSFR)_MS, where Δlog(SSFR)_MS = log[SSFR/SSFR_MS(M_⋆, z)]. We also investigated the correlation between Δlog(SSFR)_MS and λ_peak for our sample and found it to have the same correlation strength as the correlation between D_MS, as defined above, and λ_peak. Without uncertainties on individual measurements for Δlog(SSFR)_MS in Magnelli et al. (2014), we are unable to directly compare those results to our analysis on the relative significance of the correlation. We also note that the stacking techniques used in their analysis on Herschel data may be prone to temperature-dependent biases. For instance, the stacking results of Viero et al. (2013) result in fundamentally different (and overall colder) SEDs than the stacking results of Béthermin et al. (2015) despite the use of very similar data sets (P. M. Drew et al. 2021, in preparation). Further analysis of SEDs drawn from broad literature sources is beyond the scope of this work.

Note that our findings—that Σ_IR correlates most closely with λ_peak of a dust SED—is in agreement with previous work at lower redshifts, for which a relationship between the star formation surface density, Σ_IR, and dust temperature exists (Lehnert & Heckman 1996; Chanial et al. 2007; Lutz et al. 2016; see also Díaz-Santos et al. 2010, which shows that nearby galaxies with higher L_IR have distinctly more compact sizes). While the results of Chanial et al. (2007) find a relatively shallow relationship between Σ_IR and T_d, our results are more consistent with the earlier (steeper) dependence suggested by Lehnert & Heckman (1996); the relative discrepancies could be due to measurement of sizes conducted in radio continuum versus millimeter continuum (of which the literature suggests a discrepancy; e.g., Miettinen et al. 2017) or differences in the sample selection.

In particular, we find it useful to analyze the results of Lutz et al. (2016) in context, which is based on FIR size measurements from Herschel/PACS of local galaxies; they derived relationships between Σ_IR, L_IR, and sSFR/sSFR_MS and a proxy measure of dust temperature, the FIR color $\mathrm{log}({S}_{70}/{S}_{160})$, which we will call C. Their Table 2 provides a summary of the derived relationships, allowing us to directly infer the relative correlations between the FIR color C and Σ_IR, L_IR, and Δlog(SSFR)_MS (even though FIR color C is not equivalent to λ_peak). What the Lutz et al. scaling relations suggest is that C and Σ_IR are correlated with a measurement significance of 30σ, C and L_IR are correlated with a measurement significance of 20σ, and C and Δlog(SSFR)_MS are correlated with a measurement significance of 6σ. While the exact significance of the correlations is likely dependent on sample details, the overall trend is similar to our findings: that Σ_IR is the most closely linked physical quantity to dust temperature (or λ_peak or C) and that correlation to the galaxies' main sequence (here probed by Δlog(SSFR)_MS) is the weakest of the correlations.