No Redshift Evolution of Galaxies' Dust Temperatures Seen from 0 < z < 2

We build our lowest-redshift sample, which we will refer to as the IRAS sample, by performing a spatial cross match between two separate catalogs. The first catalog contains IRAS (Neugebauer et al. 1984) sources from the Revised IRAS Faint Source Redshift Catalog of galaxies (RIFSCz; Wang et al. 2014). This sample is selected at 60 μm and contains observations at 12, 25, 60, and 100 μm covering galactic latitudes ∣b∣ > 20° in unconfused regions of the sky at 60 μm (where "unconfused" corresponds to fewer than eight sources per beam with 6' radius) with order of 60,000 sources. The second catalog is data releases 1 and 2 of the Herschel Astrophysical Terahertz Large Area Survey catalog, which covers ∼120,000 sources (H-ATLAS; Bourne et al. 2016; Valiante et al. 2016; Smith et al. 2017; Furlanetto et al. 2018; Maddox et al. 2018). This sample is selected at 250 μm and contains observations at 100, 160, 250, 350, and 500 μm from the SPIRE and PACS instruments covering 660 deg².

To create the RIFSCz, Wang et al. (2014) selected galaxies detected with IRAS from the IRAS Faint Source Catalog (FSC; Moshir et al. 1993) that have S/N > 5 in the 60 μm band. The FSC contains ∼1.7 × 10⁵ sources with only about 40% of sources detected at 60 μm. In order to obtain accurate positional estimates, Wang et al. (2014) used the likelihood ratio technique (e.g., Sutherland & Saunders 1992) to match FSC sources to 3.4 μm counterparts in the Wide-field Infrared Survey Explorer (WISE) All-Sky Data Release (Wright et al. 2010). This wavelength is used because it is the deepest band in WISE and it samples a galaxy's SED at a wavelength where it is likely to be detected from either warm dust or stellar emission mechanisms. The authors then cross match the IRAS sources with WISE astrometry (the latter is more precise) and then to several other wide-field surveys: the Sloan Digital Sky Survey Data Release 10 (SDSS DR10; Ahn et al. 2014) using a matching radius of 3'', the Galaxy Evolution Explorer All-Sky Survey Source Catalog (Seibert et al. 2012) using a matching radius of 5'', the Two Micron All Sky Survey (2MASS Huchra et al. 2012), which is already included with the WISE catalog, and the first public release of the Planck catalog of compact sources (Planck Collaboration 2013; Planck Collaboration et al.2014) using a matching radius of 3'. Each matching radius used by Wang et al. (2014) is chosen based on the spatial resolution of each match. Sources are retained in the RIFSCz catalog even if they are not matched to detections in these ancillary catalogs. See Wang et al. (2014) for more details. In the present work, focusing on IR SEDs, we use IRAS and WISE flux densities and redshifts for our analysis.

We further cross match RIFSCz sources with the H-ATLAS catalog (Bourne et al. 2016; Valiante et al. 2016; Smith et al. 2017; Furlanetto et al. 2018; Maddox et al. 2018) to fill out the SEDs at longer wavelengths. H-ATLAS presents sources detected at S/N > 4 at 250 μm (corresponding to a 4σ depth of 29 mJy) in an area of the sky spanning 660 deg². The catalog contains two data releases: DR1, covering the North Galactic Pole (NGP; 500 deg²), cross matched with SDSS DR10 (Ahn et al. 2014) for redshifts and optical counterparts and positions, and DR2, covering the South Galactic Pole (SGP; 162 deg²).

Sources in the RISFCz catalog have astrometric positions based on the best available precision where reliable counterparts are identified in Wang et al. (2014): either from SDSS DR10, 2MASS, or WISE. From the Wang et al. (2014) astrometry, we further cross match RIFSCz sources with both DR1 and DR2 of the H-ATLAS catalog using search radii that correspond to the radius of the average point-spread function (PSF) of the best available ancillary catalog provided by Wang et al. (2014). H-ATLAS DR1 and the NGP field of DR2 are matched to RIFSCz sources with a search radius of 1'' because the RIFSCz catalog in this region is based on SDSS DR10 astrometry, precise to 1''. However, there is no SDSS DR10 coverage of the SGP field in H-ATLAS DR2, so sources in the SGP field are matched using a less restrictive search radius of 16 for counterparts in the 2MASS catalog whose astrometry is slightly less well constrained (154 sources are matched). In the absence of 2MASS counterpart matches, sources from the SGP are matched to the WISE astrometry given in the RIFSCz using a matching radius of 3'' (14 sources are matched this way). Note that contamination from false positives is negligible within this search radius given the low source density of both RIFSCz sources (∼1 deg⁻²) and possible 2MASS counterparts (∼50 deg⁻²). Cross matching sources from the RIFSCz to H-ATLAS data releases 1 and 2 results in 780 matches or 2% of RIFSCz sources, which is the expected proportion of the full RIFSCz catalog because H-ATLAS only covers 2% of the solid angle spanned by the RIFSCz.

We limit the IRAS sample to only include galaxies with spectroscopic redshifts (spec-zs). This requirement reduces the sample from 780 to 750 galaxies (96% of the original) spanning the redshift range z < 0.48 with median $\left\langle z\right\rangle =0.050\pm 0.002$. This reduction does not show a bias in redshift within the sample. Requiring spec-zs reduces the sample size by only 4% and thus is a reasonable choice because uncertainties in adopted redshift will translate into uncertainties on our measured rest-frame peak wavelengths and IR luminosities.

Because the aim of this work is to look for redshift evolution in the correlation between the peak wavelength and the IR luminosities of galaxies, we want a clean snapshot of the L_IR–λ_peak correlation without influence from any potential redshift evolution. We therefore further restrict the sample to only include galaxies in a redshift range that does not evolve strongly with cosmological time and yet is beyond the local volume. We limit our sample to galaxies between 0.025 <z < 0.19. The lower-redshift bound is specifically chosen to remove spatially resolved sources. It corresponds to the redshift where a galaxy with a radius of 10 kpc would be resolved across two 250 μm Herschel beams (18'' × 2 = 36''). This removes some of the uncertainty associated with the more extreme aperture corrections required to derive accurate photometric constraints in the local volume. The upper-redshift bound is chosen based on the time frame over which the star formation rate of galaxies on the star formation "main sequence" evolves by less than a factor of two from our lower-redshift bound (Speagle et al. 2014). In other words, the galaxies' average SFRs (at fixed stellar mass) do not evolve by more than two times over this interval. This factor of two is chosen because it corresponds to the intrinsic width of the main sequence itself (∼0.3 dex). Thus, in this snapshot in time, the main sequence does not evolve by more than its intrinsic width. This additional redshift restriction leaves 511 out of 750 galaxies remaining in the sample (68%). We also tested setting the lower-redshift bound to the redshift where a galaxy with radius 10 kpc fits within a single beam (at z = 0.01) and found this does not change the overall results of this paper but does increase the uncertainty in our measurement of the L_IR–λ_peak relation, as it limits the IRAS sample to only 257 out of the 780 original galaxies (33%). We therefore continue with the less restrictive 0.025 < z < 0.19 selection to maximize the dynamic range of luminosities in our sample. We find that this sample is >90% complete above >340 mJy and >81% complete above >179 mJy at 60 μm. Figure 1 shows the distribution of redshifts versus the IR luminosities for the IRAS sample along with the other samples that we will describe below in Sections 2.2 and 2.3.

Observational constraints included in the IRAS sample are WISE 12 μm, IRAS 12 μm, WISE 22 μm, IRAS 25 μm, IRAS 60 μm, IRAS 100 μm, Herschel 100 μm, Herschel 160 μm, Herschel 250 μm, Herschel 350 μm, and Herschel 500 μm, with inclusion in the sample predicated only on a detection at 60 μm. These wavelengths allow us to probe approximately the full range of dust emission. Because the 12 μm band may capture PAH emission at the redshifts covered by this sample, we test that the inclusion of two 12 μm bands—from both IRAS and WISE–does not bias our measurement of IR luminosities and peak wavelengths. We find that the median IR luminosities and peak wavelengths of galaxies in this sample differ by less than a percent when excluding either set of observations at 12 μm. This difference is smaller than the average uncertainty in these parameters, so we conclude that the inclusion of both 12 μm bands does not introduce additional bias.

To summarize, our selection function for this sample is galaxies that are (1) detected at S/N > 5σ at 60 μm in the IRAS FSC, (2) are present in the H-ATLAS catalogs (S/N > 4σ at 250 μm), (3) have a spec-z in either RIFSCz or H-ATLAS, and (4) are in the redshift range 0.025 < z < 0.19, taken to represent a snapshot in time of the local universe. See Table 1 for a list of the wavelengths included and selection limits for this sample as well as the other samples described below. Figure 2 illustrates the detection limits for each sample in L_IR–λ_peak space at their median redshifts with contours showing the distribution of sources. The goal of these selections is to create samples that are minimally biased with respect to SED dust temperature at a given IR luminosity. We will discuss this further in Section 5.

Zoom In Zoom Out Reset image size

Finally, in order to test that the IRAS sample selection and observed filter set is not biased with respect to dust temperature, we generate mock SEDs covering a grid in L_IR–λ_peak space at peak wavelengths and IR luminosities that are broader than those we find in the sample, namely peak wavelengths between 35 and 200 μm, which corresponds to dust temperatures between 103 and 12.5 K, respectively. We find no offset between the inputs to these noise-perturbed mock sources and the measured values for L_IR and λ_peak at the signals to noise representative of the IRAS sample with the wavelength sampling of the IRAS sample. In other words, over a range of mock SEDs that is more diverse than the actual observed IRAS sample, we find no bias due to the wavelength coverage or S/N of the IRAS sample. If there were galaxies that populated parameter space in this grid we would be able to accurately measure their luminosities and peak wavelengths. We return to a discussion of this test, with respect to any bias in this combination of filters, in Section 4.2 after introducing the mock galaxy samples in more detail.

Our second sample, which we will refer to as the H-ATLAS sample, expands the redshift horizon of this study beyond z ∼ 0.2, probing redshift evolution out to z ∼ 0.4. The sample is defined by galaxies detected at 250 μm with S/N > 4 over an area of 161.6 deg² from data release 1 of H-ATLAS (Bourne et al. 2016; Valiante et al. 2016). For this sample, we make exclusive use of H-ATLAS DR1 sources because Bourne et al. (2016) performed a careful cross matching with WISE photometry (Wright et al. 2010), giving us important constraints on the sample's short wavelength dust SED. There are of the order of 70,000 galaxies in H-ATLAS DR1. We further restrict our analysis to galaxies that are detected at 100 μm with S/N > 3, which significantly helps constrain the galaxies' λ_peak and reduces the sample to 3407 galaxies (5% of all H-ATLAS DR1 galaxies). The sky coverage with PACS and SPIRE in the H-ATLAS DR1 is roughly the same, but the PACS instrument is much less sensitive than SPIRE. This accounts for the relative dearth of galaxies with S/N > 3 at 100 μm. We will discuss the impact of this 100 μm detection cut further in Section 4.2, but to summarize here, it effectively reduces the dynamic range in IR luminosity to the same range used for the IRAS sample while also dramatically improving the SED fit quality by adding a precise constraint near the peak. In other words, the requirement of a 100 μm detection removes low-L_IR sources below the IR luminosity range studied in this manuscript.

Although H-ATLAS DR1 includes WISE observations at 12 and 22 μm, we use only the 22 μm band in our SED fits because PAH emission is expected to contribute to the 12 μm band, depending on the redshift, and our modeling does not account for PAH emission. In the IRAS sample, we included WISE 12 μm, IRAS 12 μm, WISE 22 μm, and IRAS 25 μm because we have as many as five to six constraints at wavelengths shorter than the peak wavelength of emission, so possible PAH contamination was downweighted by the other observations. On the short wavelength end, the H-ATLAS sample only has observations at 12 μm, 22 μm, and 100 μm, so PAH contamination in the 12 μm band is more likely to skew our analysis. While we do not require a detection with WISE, we do limit our sample to galaxies observed with WISE. Excluding those without WISE coverage removes 723 galaxies, reducing our sample to 2679 galaxies (by 22%). This sample covers 0 < z ≲ 0.8, but we limit our analysis to galaxies between 0.05 < z < 0.4, due to a significant drop in H-ATLAS survey sensitivity beyond z ∼ 0.4; i.e., H-ATLAS is only sensitive to ULIRGs (>10¹² L_⊙) beyond z ∼ 0.4, but their sky density at those redshifts is significantly lower, making average dust SEDs difficult to constrain. As we did for the IRAS sample, we use a lower-redshift bound to avoid including spatially resolved sources. For this sample, we choose a lower limit of z = 0.05 because this is the redshift where a galaxy with radius 10 kpc is unresolved by Herschel at 250 μm (i.e., within one beam). This lower limit reduces the number of galaxies to 2301, and the upper limit reduces the number of galaxies to 2167 (3% of all H-ATLAS DR1 galaxies).

Valiante et al. (2016) estimated the completeness of DR1 to be greater than 95% for galaxies with flux densities above 42.15 mJy at 250 μm. This corresponds to 94% of galaxies in our sample. To address the possible inclusion of false positives, or galaxies whose FIR emission is incorrectly matched to an optical/near-infrared source within the redshift range of interest at z < 0.4, we evaluate the likelihood of random, nonphysical associations in existing optical/near-infrared imaging in the H-ATLAS survey. We find the likelihood of a random match to an OIR source at the depth of SDSS DR10 data (22nd magnitude) at z < 0.4 to be 0.7% within a 3'' matching radius, which would correspond to 15 sources. We also investigate whether it is likely that there are sources present in H-ATLAS whose optical counterparts are not detectable in the available optical imaging from SDSS DR10. ² The infrared excess, the ratio of the IR to the UV luminosities (L_IR/L_UV ≡ IRX) of sources at z = 0.4 that have U-band flux densities at the detectability limit of SDSS DR10 (AB magnitude 22) and Herschel 250 μm fluxes at the detectability limit of H-ATLAS, corresponds to IRX ∼ 4 ×10⁵, which is much higher than a more typical value of IRX for DSFGs in the literature ∼10³–10⁴ (e.g., Howell et al. 2010; Casey et al. 2014). Adopting an IRX value of 10⁴ implies an AB magnitude limit of 18, implying all OIR counterparts should be detected in SDSS DR10.

The included wavelengths of observation for the H-ATLAS sample are then 22, 100, 160, 250, 350, and 500 μm. See Figure 1 for the redshift distribution versus IR luminosity of the H-ATLAS sample. Figure 2 shows the detection limits at the median redshift of the H-ATLAS sample with contours showing the distribution of sources in L_IR–λ_peak space. To summarize, our selection function for this sample is galaxies that have (1) S/N > 4 at 250 μm and S/N > 3 at 100 μm in H-ATLAS DR1, (2) are observed with WISE (though no detection is necessary), and (3) lie between 0.05 < z < 0.4.

Our third sample, which we will refer to as the COSMOS sample, extends out to z = 2. It contains Herschel-selected galaxies from the COSMOS field (Scoville et al. 2007) originally presented in Lee et al. (2013; hereafter L13). L13 measures Spitzer MIPS flux densities from Sanders et al. (2007), Herschel SPIRE and PACS flux densities from Oliver et al. (2012) and Lutz et al. (2011), 1.1 mm AzTEC (Scott et al.2008; Aretxaga et al. 2011), and 1.2 mm MAMBO (Bertoldi et al. 2007) flux densities. They use positional priors from 24 μm (Le Floc'h et al. 2009) and 1.4 GHz (Schinnerer et al. 2007) following the cross-identification (XID) linear inversion technique (Roseboom et al. 2010, 2012). The authors then match the 24 μm positions to K_s -band counterparts from the COSMOS-WIRCam near-infrared imaging survey (McCracken et al. 2010) using a 2'' search radius. Finally, these K_s -band positions are matched with a 1'' search radius to the photometric redshift catalog of Ilbert et al. (2009) to derive photometric redshift constraints. This results in a total of 39,333 sources. Most of these are not significantly detected in the IR/millimeter and so are excluded from the analysis in this paper upon applying further selection cuts.

In the present paper, we update the L13 sample and methodology using more recent COSMOS data releases. We start with the 4220 galaxies in the L13 catalog with S/N > 3 in at least two PACS or SPIRE bands (11% of all L13 galaxies). L13 uses only galaxies with S/N > 5 in at least two PACS or SPIRE bands (1810 galaxies, using confusion limits as a floor on the errors on flux density). We explored both S/N cuts for our final sample and found that pushing down to S/N > 3 added significant dynamical range to the luminosity limit of our analysis and does not detract from the quality of our SED constraints. We then cross match L13 positions for these 4220 galaxies with Atacama Large Millimeter/submillimeter Array (ALMA) sources in the A³COSMOS catalog (Liu et al. 2019) given the significant uncertainty in the initial 24 μm/1.4 GHz association with Herschel flux densities. Magdis et al. (2011) showed that 98% of Herschel sources at z < 2 will be detected in deep Spitzer 24 μm imaging. Similarly, ALMA detections centered on Herschel-selected sources are physically associated due to the rarity of both populations on the sky. Therefore the likelihood that an ALMA source located within the MIPS 24 μm PSF radius is not associated with the galaxy responsible for 24 μm emission at z < 2 is low (<1%). A³COSMOS contains all available ALMA archival continuum pointings in the COSMOS field with detections of 756 unique galaxies across ∼280 arcmin². We use a search radius of 3'', the beam radius of the Spitzer 24 μm band, to match 24 μm positions to ALMA-detected counterparts. ALMA detections were found for 147 sources in L13 (4% of the original 4220 sources). Of these, we find that 20 sources (14% of the 147 matches) were misidentified in the L13 sample based on the adopted positional prior. In other words, the original OIR counterpart identified within the 24 μm beam by L13 was found to be spatially offset from an ALMA source, also within the 24 μm beam, by larger than 20 kpc for 14% of the subsample with ALMA measurements available. These mismatches did not affect the photometry of these sources, as the XID linear inversion technique simply measures the fluxes in all bands at the positions of 24 μm Spitzer sources. A new matched position within the beam of 24 μm emission only results in a new matched OIR counterpart and hence a new redshift.

After adopting the new positions from the A³COSMOS subset, we cross match all source positions (including those without ALMA data) with the COSMOS2020 photometric redshift catalogs (Weaver et al. 2022). Their photometric redshifts are updated to make use of the near-IR and optical data from the UltraVISTA and Subaru HSC surveys' deeper imaging in the field. The COSMOS2020 photometric redshift (photo-z) catalog consists of two separate photometric catalogs generated using different source extraction techniques as well as two photo-z fitting techniques applied to those two catalogs, resulting in four unique photo-z catalogs. See Weaver et al. (2022) for more details. Following the procedure adopted by L13 and Le Floc'h et al. (2009), we perform our spatial cross match from long-wavelength astrometry (ALMA if available, otherwise Spitzer 24 μm) separately on all COSMOS2020 catalogs using a 2'' search radius, and adopt the nearest neighbor match. These works found that by adopting a 2'' search radius, they were able to reduce the fraction of multiple matches while still taking into account most of the positional uncertainty associated with the 3'' width of the MIPS 24 μm PSF. For sources with ALMA positions, we adopt a search radius of 1'', as this is the size of the average ALMA beam in A³COSMOS. Thus, from COSMOS2020, we have between one and four photometric redshift estimates for each galaxy. There are 4011 galaxies matched out of the original 4220 sources with S/N > 3 in at least two Herschel bands (95%). Next we remove galaxies that are marked as masked in both source catalogs (for proximity to a bright star that obfuscates the OIR source photometry), leaving only 2295 galaxies of the original 4011 with matches (57%). Rather than fitting these masked galaxies with the original L13 photo-zs, we choose to only retain the subset of sources matched to COSMOS2020 to retain consistency in how redshifts are measured within the sample. We do not find that this cut hinders our final analysis in any way, as the remaining sample is still sufficiently large.

Where spectroscopic redshifts are available in the COSMOS2020 catalog, we use them instead of photometric redshifts. Of the 2295 sources, 285 have spec-zs (12%). For the remaining galaxies, we adopt a single photo-z. For those with two photo-z estimates from the COSMOS2020 catalogs, we use the lower z because galaxies in this sample are statistically more likely to be at lower z. For galaxies with three photo-z estimates, we adopt the median redshift. And finally, for galaxies with four photo-z estimates, we adopt the redshift that is the closest to the median. We do this rather than adopting the true median so that we may retain the redshift uncertainties included with the COSMOS2020 catalogs. There are 2005 sources with four photo-z estimates, 83 with three photo-z estimates, 206 with two photo-z estimates, and one source with one photo-z estimate. The median standard deviation of photo-z estimates is 0.02 for sources with four estimates. Note that only a minority of sources (423 of 2295, or 18%) have standard deviations in their two to four photometric redshift estimates greater than 0.1. Additionally, the majority of our galaxies have very similar photometric redshift estimates to those from L13, with σ_Δz/(1+z) = 0.009. See Figure 1 for the distribution of COSMOS sample redshifts and IR luminosities relative to the other samples.

Galaxies in the COSMOS sample span 0.02 < z < 4.8. We limit our analysis to 0.15 < z < 2 due to the poor sensitivity of Herschel for IR-luminous sources beyond z ∼ 2 and potential bias in the dust temperature measured from those SEDs beyond that epoch. Unlike the H-ATLAS sample where a lower-redshift bound of 0.05 is chosen to exclude galaxies resolved with Herschel, for the COSMOS sample, we increase our lower-redshift bound to 0.15 because the COSMOS field is much smaller than that of H-ATLAS, and a small area provides a low number of sources per given bin of Δz below z < 0.15. Whereas there are 72 sources between 0.05 < z < 0.15 in the COSMOS sample, there are 1314 sources in the H-ATLAS sample over this same range. Low number statistics between 0.05 < z < 0.15 are likely to result in a bias to our analysis, particularly given the lower IR luminosities of that subsample, so we require >100 galaxies per Δz = 0.1 bin on the low-redshift end of the COSMOS sample. We achieve this by adopting a lower-redshift bound of 0.15. This reduction takes the COSMOS sample to 1990 sources.

The last step in our sample selection procedure is to cross match the COSMOS sample with sources from the S2COSMOS catalog (Simpson et al. 2019) using a search radius of 65, the beam radius of SCUBA-2 at 850 μm, in order to fold in their 850 μm data to our SED analysis. This results in 124 detections with S/N > 3, and nine with S/N > 5. For galaxies in the COSMOS sample that do not have a matched S2COSMOS source, we adopt the flux in the S2COSMOS 850 μm map at the location of our source. For the uncertainty on this flux density, we adopt the rms from the S2COSMOS 850 μm rms map. The wavelengths of observation for our COSMOS sample are then 24, 70, 100, 160, 250, 350, 500, 850, 1100, and 1200 μm. See Figure 1 for the redshift distribution of sources in the COSMOS sample versus IR luminosity, and see Figure 2 for the detection limits at the median redshift of the COSMOS sample with contours showing the distribution of sources in lirlp space. We find that this sample is >90% complete at flux densities >10 mJy at 250 μm. To summarize, our selection function for this sample is galaxies from L13 that have (1) S/N > 3 in at least two of five Herschel bands (100, 160, 250, 350, and 500 μm), and (2) have redshifts between 0.15 < z < 2.0.

The IRAS sample serves as our low-redshift anchor to which we compare the H-ATLAS and COSMOS samples that extend farther out to high redshift. The left panel of Figure 4 shows the L_IR–λ_peak correlation fit for the IRAS sample. The right panel of Figure 4 shows a histogram of the log of the peak wavelength measurements minus the log of the best-fit model peak wavelengths, illustrating the scatter about the L_IR–λ_peak relation. The distribution is well fit by a Gaussian with a tail of galaxies at shorter peak wavelengths, or hotter dust temperatures. AGNs are a possible extra source of heating for dust in galaxies. They are capable of heating torus dust to temperatures from 50 to 100 K to as high as the dust sublimation temperature, ∼1000–1500 K (e.g., Donley et al. 2012). The galaxy-integrated dust temperature may or may not be shifted to warmer temperatures because of AGNs, depending on the emergent AGN luminosity. To remove bias in our fits to the L_IR–λ_peak correlation due to the presence of the hotter tail of sources in the IRAS sample, we iteratively sigma-clip sources at >3σ, consisting of 27 out of 511 sources (5%). We find that approximately 26% of these outliers (seven out of 27) contain known AGNs by cross referencing with the first bright quasi-stellar object survey (Gregg et al. 1996), the 13th edition of the catalog of quasars and active nuclei (Véron-Cetty & Véron 2010), and the all sky catalog of WISE-selected AGNs from Secrest et al. (2015). Note that the fraction of cold ("inliers") similarly identified as AGNs in the same catalogs is significantly lower (7%, or 35 of 484) than 26%. For the other samples, we do not sigma-clip, as we do not find a substantial population of outliers like those in the IRAS sample. The reason they are present in the IRAS sample is likely due to the 60 μm selection. Figure 2 shows how a 60 μm selection is more sensitive to low-luminosity warm outliers compared to the selection functions that drive the H-ATLAS and COSMOS samples. Though AGNs are present in our samples and have been shown to account for about 10%, on average, of the total IR luminosities of their host galaxies (e.g., Iwasawa et al. 2011), they do not seem to affect the global relationship measured here. Note that if we were not to clip these outliers from our sample, the derivation of the Θ_IRAS parameter set would still be consistent with our findings within measurement uncertainty; however, we find that this clipping provides the most appropriate fit by excluding sources known not to follow the bulk relationship drawn from a normal distribution in $\mathrm{log}({\lambda }_{\mathrm{peak}})$ about a given IR luminosity.

Zoom In Zoom Out Reset image size

The best-fit parameters for the IRAS sample, Θ_IRAS, are given in Table 2. In the following subsections, we will use these values as our basis for comparison with our other samples to investigate redshift evolution.

Beyond the nearby universe, the detection limits of FIR/submillimeter facilities can, depending on the wavelength of observation, be biased with respect to the galaxies' dust temperatures, rendering either very cold or very warm SEDs undetectable. This is already demonstrated in our analysis by the presence of likely AGN-heated hot sources in the IRAS sample (driven by the 60 μm selection) not present in the H-ATLAS and COSMOS samples. If one tries to fit the L_IR–λ_peak correlation for subsamples of galaxies selected at different wavelengths without first accounting for these relative selection biases, the measured fit parameters Θ may themselves be biased. Such a bias is a function of sample selection wavelength, detection limit, and redshift. Figure 2 shows the detection limits at the median redshift of each sample plotted with contours of detected galaxies across all redshifts in each sample. For galaxies to be detected in a given sample, they need to have higher IR luminosities than indicated by the thick colored lines, which denote the flux density selection process outlined in Table 1. For example, to be detected in the COSMOS sample, galaxies must lie to the right of at least two curves in the rightmost panel of Figure 2: each curve represents the detection limit in each of the five bands used to define the catalog. For the H-ATLAS sample, the inclusion of the 100 μm detection requirement acts like a cut in IR luminosity; it is insensitive to dust temperature. We find this additional selection greatly improves our constraints on the measurement of peak wavelengths of galaxies in this sample.

To determine the degree to which measuring Θ at any given redshift is biased by selection functions, we generate mock galaxies to represent each data set in this paper. These mock galaxies are modeled such that they follow the local IRAS sample derived L_IR–λ_peak correlation. In other words, we assert in the mock samples that the L_IR–λ_peak correlation does not evolve with redshift. These mock galaxies are generated using the Θ parameter set fit to the IRAS sample, and are therefore used to test the null hypothesis, i.e., that there is no intrinsic redshift evolution in Θ parameters.

In order to generate representative galaxies for each sample, we use the infrared luminosity function from Zavala et al. (2021) based on the Casey et al. (2018b) model to calculate the expected number of galaxies of each IR luminosity that exist in the given survey solid angle. Note that these IR luminosity functions are consistent with many other IR luminosity functions in the literature at z ≲ 2 (e.g., Magnelli et al. 2009; Casey et al. 2012a, 2012b; Gruppioni et al. 2013; Magnelli et al. 2013). We randomly assign IR SEDs to the mock galaxies by sampling a normal distribution in log(λ_peak) with width σ_λ around the best-fit L_IR–λ_peak correlation and with α and β values drawn from the posterior distributions of the IRAS sample. Photometry is then generated for the mock galaxies at the wavelengths where observations exist for each sample with appropriate instrument and confusion noise added. We then remove mock galaxies that fall below the detection limits of the survey as defined by the thresholds in Table 1 in each redshift bin. Thus, we generate mock versions of each sample in both survey volume and selection biases as they would exist if the null hypothesis were true. Then, if the measured Θ from the mock samples are statistically consistent with the Θ measured from the observed samples in a given redshift bin, then that redshift bin is consistent with the null hypothesis, i.e., that there is no redshift evolution in the L_IR–λ_peak relationship since z ∼ 0.

We test for selection bias in the IRAS sample in this same manner, though by construction we might expect agreement, given that mock galaxies are drawn from the derived parameter set Θ fit to the IRAS sample. This would not be the case if there was a strong selection bias, however. We find that the fits to the mock IRAS sample are consistent within errors with the measurements from the observed sample. This suggests that the detection limits for the IRAS sample are not presenting a strong bias for our measurement of Θ in the nearby universe. Additionally, as is detailed in Section 2.1, we find no offset between the inputs and best-fit values of mock galaxies that we generate in a grid covering a larger range of L_IR–λ_peak space than exists in the IRAS data. Together, these two tests suggest that there is not a strong bias for our measurement of Θ in the nearby universe.

With mock samples that match our observed H-ATLAS and COSMOS samples in hand, we investigate whether or not there is evidence of redshift evolution in the L_IR–λ_peak relation. To do this, we divide both mock and observed samples into discrete redshift bins for which we measure Θ. We choose bin widths that provide an adequately large dynamic range in the subsamples' L_IR (≳1 dex) and at least a few hundred galaxies per bin, both of which are needed to accurately constrain the parameters of Θ after accounting for the added uncertainty due to selection biases. For the H-ATLAS sample, we use redshift bins of width Δz = 0.1 between 0.05 < z < 0.4 with the lowest redshift bin truncated to 0.05 < z < 0.1. For the COSMOS sample, we use bins with width Δz = 0.5 ranging between 0.15 < z < 2, again with a truncated lowest redshift bin with Δz = 0.35.

Figure 5 shows the H-ATLAS and COSMOS samples in the L_IR–λ_peak plane in purple alongside mock galaxies made to mimic them in orange. We find that for the most part, the mock galaxies qualitatively resemble the distribution of real galaxies in both the H-ATLAS and COSMOS samples. The black dashed line shows the best fit to the IRAS sample, and the orange and purple lines show 100 random samples from the MCMC fits to real and mock galaxies, respectively. Because the mock galaxies are drawn from a distribution mirroring the nearby IRAS sample, we test for redshift evolution by comparing the best-fit parameter sets of the real samples to the mock galaxy samples. If the two fits are statistically consistent within a 95% confidence interval (i.e., 2σ), that would constitute direct observational evidence that the L_IR–λ_peak relation does not evolve.

Zoom In Zoom Out Reset image size

Figure 6 compares the best-fit parameter set Θ = (λ_t, η, σ_λ ) from Equation (3) between real and mock galaxies in the H-ATLAS and COSMOS samples. Because we find that variations in Θ associated with generating a volume limited set of mock galaxies are larger than the uncertainties derived when fitting any one set of mock galaxies (∼2× higher than uncertainties), we simulate each redshift subsample for H-ATLAS and COSMOS 100 times using input Θ that are drawn from the posterior distribution of the IRAS sample, Θ_IRAS. We then measure Θ in each trial to generate a distribution of measured Θ for each subsample. Similarly, we also fit 100 bootstrap iterations of the real data for each subsample to generate realistic joint probability distributions of Θ.

Zoom In Zoom Out Reset image size

To statistically assess the agreement or disagreement of our measured parameter set Θ for a given sample and redshift bin, we compute the ratio of slopes and intercepts between mock and real samples, η_mock/η_real and λ_t,mock/λ_t,real. If mock and real galaxies are statistically consistent, these ratios will be approximately equal to one. We find that indeed, η_mock/η_real and λ_t,mock/λ_t,real are consistent with 1 within the 95% confidence interval (2σ) for all subsamples and redshifts. In other words, we rule out the possibility of an evolving L_IR–λ_peak relationship at a confidence of >99.9997% (or 4.6σ significance), which we calculate by integrating the two-dimensional probability density distribution functions in derived parameters η and λ_t at values that would be >3σ offset from those measured in Θ_IRAS. Thus we conclude for all subsamples that the data are consistent with no redshift evolution in the correlation between the IR luminosities of galaxies and the peak wavelength of their dust emission between 0 < z < 2.

The rightmost panels of Figure 6 show the width of the L_IR–λ_peak distribution in each redshift bin for all samples. The best fit to the IRAS is shown as a green dashed line, the scatter of observed galaxies about the L_IR–λ_peak correlation, σ_λ , is shown in purple, and the scatter of the distributions of mock galaxies is shown in orange. The measured σ_λ of the real galaxies are consistent with the measurements both from the IRAS sample and the mock samples for each redshift bin.

We interpret our results to mean that IR-luminous galaxies at high-redshifts are not dramatically different in size or geometry than those in the local universe, assuming the dust opacity does not significantly vary between sources or does not evolve with redshift. The z ∼ 0 relation we measure in L_IR–λ_peak from the IRAS sample exhibits only a subtle IR–luminosity dependency on λ_peak as evidenced by measured slope of η = − 0.09 ± 0.01. For context, over four decades of luminosity, from 10⁹ L_⊙ to 10¹³ L_⊙, the mean λ_peak measured is only expected to change by a factor of 2.3 ± 0.3. The intrinsic scatter about the relation, measured via the parameter σ_λ , is of the order of the same factor (×1.1). The subtlety of the relationship and large intrinsic scatter can be interpreted with some understanding of galaxies' dust-emitting sizes and source geometries.

The size and geometry of a dust-emitting region relate directly to the emergent IR luminosity and peak wavelength as would be expected from the Stefan-Boltzmann Law (e.g., Hodge et al. 2016; Simpson et al. 2017), where luminosity is proportional to dust temperature to the fourth power and to surface area (for an optically thick blackbody). A nice discussion of the applicability of the Stefan-Boltzmann Law to galaxies' dust emission is offered in Burnham et al. (2021), who show that IR luminosity surface density is the most fundamental tracer of a galaxy's integrated dust temperature, or λ_peak (see also Chanial et al. 2007; Lutz et al. 2016). Using the equation relating the best-fit surface density of IR luminosity to the peak wavelength of the dust SED from Burnham et al. (2021), $\mathrm{log}{{\rm{\Sigma }}}_{\mathrm{IR}}\approx 20.18-4.19\mathrm{log}({\lambda }_{\mathrm{peak}}/\mu {\rm{m}})$, we find that for a galaxy with L_IR =10¹² L_⊙, the difference between galaxies at peak wavelengths ±1σ_t corresponds to a difference of only 0.2 dex (<2×) in the effective radius of the emitting dust. A wide array of dust geometries in individual galaxies, from compact and clumpy to extended and diffuse, then accounts for the intrinsic scatter in the L_IR–λ_peak relation at fixed L_IR. This scatter translates to ±5 K temperature differences at λ_peak = 100 μm, and an effective half-light radius difference of 0.4 kpc for galaxies of typical radii of 0.8 kpc, assuming the case of an optically thick blackbody. Because galaxy sizes (as measured by dust continuum; e.g., Hodge et al.2016; Simpson et al. 2017) only vary by a factor of a few and not orders of magnitude, the L_IR–λ_peak relation itself has only shallow dependence on L_IR. In other words, in our best-fit relation for the IRAS sample, η = − 0.09 ± 0.01 implies a size–luminosity relationship such that ${R}_{\mathrm{eff}}\propto {{L}_{\mathrm{IR}}}^{0.31}$ (consistent with ${R}_{\mathrm{eff}}\propto {{L}_{\mathrm{IR}}}^{0.28\pm 0.07}$, measured by Fujimoto et al. 2017).

Though this work does not directly measure the size of the dust-emitting region in our samples, nor does it quantify the dust geometry through morphological analysis, our finding of consistent nonevolving SEDs in the L_IR–λ_peak plane argues for unchanging dust sizes from z ∼ 2 to z ∼ 0. Similarly, it argues that the breadth of sizes in those galaxies also does not evolve. This naturally follows when considering the IR-luminous population (>10¹¹ L_⊙) is, on a whole, representative of fairly massive galaxies that have likely established most of their mass reservoirs—whether in stars or gas supply that will be converted to stars—and will not continue to grow substantially in size, whether they are at z ∼ 2 or z ∼ 0.

Our finding that the L_IR–λ_peak relationship does not evolve from 0 < z < 2 is perhaps seen to be in conflict with other works in the literature that claim galaxies' temperatures evolve with redshift. One set of claims argues that higher-redshift galaxies have hotter dust temperatures compared with galaxies at z ∼ 0. This is based on observations of galaxies between 1 ≲ z ≲ 5 (e.g., Magdis et al. 2012; Magnelli et al. 2014; Béthermin et al. 2015; Faisst et al. 2017; Schreiber et al. 2018; Faisst et al. 2020), and theoretical modeling, primarily at z ≳ 5 (e.g., De Rossi et al. 2018; Arata et al. 2019; Liang et al. 2019; Ma et al. 2019; Sommovigo et al. 2020). These higher temperatures are attributed to higher sSFRs, higher star formation surface densities, or lower dust abundances relative to z ∼ 0 samples.

Another set of claims argues that dust temperatures are colder at high redshift than at z ∼ 0 (e.g., Chapman et al. 2002; Pope et al. 2006; Symeonidis et al. 2009; Hwang et al. 2010; Kirkpatrick et al.2012; Symeonidis et al. 2013; Magnelli et al. 2014;Kirkpatrick et al. 2017). This is attributed to galaxies having higher dust masses (e.g., Kirkpatrick et al. 2017), higher dust opacities or emissivities, or more extended spatial distributions of dust at higher redshift (e.g., Symeonidis et al. 2009; Elbaz et al. 2011; Rujopakarn et al. 2013) as compared with z ∼ 0.

A third set of claims argues that the dust temperatures of galaxies do not evolve with redshift. These works find no evidence of redshift evolution in L_IR–λ_peak space for their samples (e.g., Casey et al. 2018b; Dudzevičiūtė et al. 2020; Reuter et al. 2020), which we also find. How is it possible that there are such seemingly contradictory conclusions in the literature?

Our results are, indeed, consistent with more than just the third set of these claims. Observational works such as those by Magnelli et al. (2014), Béthermin et al. (2015), and Schreiber et al. (2018) found that dust temperatures increase between low and high redshift by studying samples of galaxies at fixed stellar mass on the main sequence. The star formation rates of galaxies on the main sequence increase as redshift increases (e.g., Speagle et al. 2014), which naturally results in an increase in IR luminosity. Figure 7 highlights how a nonevolving L_IR–λ_peak relation then translates to a perceived evolution in dust temperatures when drawing from main-sequence galaxies at fixed stellar mass. At high redshift, the SFR (thus L_IR) at fixed stellar mass is higher (e.g., Speagle et al. 2014), and galaxies with higher IR luminosities have hotter SEDs (e.g., Casey et al. 2018b; Burnham et al. 2021). The tracks of fixed stellar mass in this figure are derived from our L_IR–λ_peak relationship, where L_IR is converted to a star formation rate (Kennicutt & Evans 2012) then to stellar mass (using the Speagle et al. 2014 relation). We then convert λ_peak to dust temperature using a simple optically thin modified blackbody, as in Bouwens et al. (2020). We see remarkable agreement between our projected L_IR–λ_peak curves at fixed stellar masses and measurements from literature works claiming measurement of hotter SEDs at higher redshifts (e.g., Magnelli et al. 2014; Béthermin et al. 2015; Schreiber et al. 2018). While the conversion between λ_peak and dust temperature is reliant on an understanding of the underlying opacity of dust, there is roughly a linear scaling between the dust temperatures estimated using the optically thin assumption and other opacity models. So while the absolute dust temperatures may not be known, the trend toward "hotter" SEDs as perceived in Figure 7 should be relatively robust. An illustration of this can be seen with the Béthermin et al. (2015) sample. Their calculated dust temperatures (black circles) do not agree with those that we calculate for their sample (having refit them using MCIRSED; open circles), though both show monotonic increase in temperature with redshift. This offset is due to differences in the model assumed to translate λ_peak to T_d between our two works. We attribute the increase in temperature with redshift solely to the increasing SFR of galaxies on the main sequence, in line with expectations, given a nonevolving L_IR–λ_peak relation.

Zoom In Zoom Out Reset image size

In contrast to works claiming that higher-redshift galaxies have hotter temperatures, there are also a number of literature works that show higher-redshift galaxies may have cooler SEDs at fixed L_IR (e.g., Chapman et al. 2002; Pope et al. 2006; Symeonidis et al. 2009; Hwang et al. 2010; Kirkpatrick et al. 2012; Symeonidis et al. 2013; Magnelli et al. 2014; Kirkpatrick et al. 2017). Figure 8 shows the L_IR–λ_peak correlation as reported from a handful of different works in the literature, color coded roughly by redshift. Blue points represent samples at low redshift, red points represent samples of high-redshift galaxies, and cream points are representative of samples at intermediate redshifts. We have converted dust temperatures to λ_peak using the fixed values of β and λ₀ from each sample, where provided by the authors; otherwise, we use Wien's Law to convert, as was done in the case of Chapman et al. (2005), Magnelli et al. (2012), and Hwang et al. (2010). We plot these works because we are interested in the general trend with redshift, not the absolute scaling of dust temperature. From what is presented here from the literature, there is a clear offset between the lowest-redshift samples and the high-redshift samples, with some of the intermediate-redshift samples aligning with the low-redshift relation and some with the high-redshift relation.

Zoom In Zoom Out Reset image size

In this paper we find no such evolution. We attribute the difference in conclusions to the manner of constraining these galaxies' SEDs to measure λ_peak. Some of these works, in an effort to compare to the most complete local samples of galaxies, used the IRAS RBGS without requiring submillimeter photometry just redward of the peak (e.g., Casey et al. 2014). Indeed, such observations did not exist for large numbers of galaxies until the launch of Herschel. Figure 9 shows that when fitting SEDs to IRAS data alone from local samples, the measurement of λ_peak is biased toward warmer dust temperatures than when fitting to data covering the full range of the dust SED. Without data in the 160–500 μm range, intrinsically colder SEDs (with λ_peak ≳ 100 μm) become untethered, with very poor λ_peak constraints. On average, we find that IRAS photometry alone will systematically underestimate λ_peak for such systems, i.e., skewing fits toward hotter dust SEDs than they may intrinsically have. Furthermore, we find a systematic offset between IRAS 100 μm measurements and Herschel PACS observations at 100 μm. The median ratio of flux at 100 μm from PACS to that from IRAS for our IRAS sample is 1.13 ± 0.01), shown in Figure 9. This systematic bias for the sample's 100 μm flux densities results in warmer fitted dust temperatures when SEDs are fit to IRAS data exclusively rather than the full suite of IRAS, WISE, and Herschel described in Section 2.1. With the introduction of Herschel's photometric constraints on the Rayleigh–Jeans tail, such biases are eliminated. The stringent requirement we place on our local universe IRAS sample (to have Herschel-SPIRE constraints from H-ATLAS) translates to colder net SEDs in the L_IR–λ_peak relation than had been calculated previously. We have effectively removed the bias whereby colder systems had been fit to SEDs with warmer temperatures for lack of long-wavelength data.

Zoom In Zoom Out Reset image size

One last comparison we offer is a comparable fit to the Arp 220 system, one of the best-known ULIRGs in the nearby universe. Literature works often describe Arp 220 as potentially having a hotter dust SED than the average ULIRG at high redshift (λ_peak,Arp220 ∼ 60 μm versus λ_peak,SMG ∼ 70–80 μm; e.g., Dudzevičiūtė et al. 2020). However, when we refit Arp 220's well-sampled IR/millimeter SED (e.g., Eales et al. 1989; Rigopoulou et al. 1996; Sargsyan et al. 2011; Brown et al. 2014) using MCIRSED, we find that, in fact, it is well within the expected scatter for galaxies of its luminosity at any redshift 0 < z < 2 (see Figure 8).

We do not analyze galaxies beyond z ∼ 2 because Herschel-SPIRE observations beyond that redshift are limited to very bright, rare sources (L_IR ≳ 10¹³ L_⊙). Extending the analysis to higher redshift would require sufficient dynamic range in L_IR with full SED sampling. Sources in our samples at z > 2 do not provide enough dynamic range per redshift bin to perform this analysis, and other literature samples either lack statistics or an unbiased sampling of SEDs. We have demonstrated the importance of using observations that bracket the peak wavelength of the dust SED when trying to accurately measure λ_peak. Unfortunately, other samples do not bracket the peak wavelength of dust emission until z ∼ 7, when the observed-frame peak reaches ∼850 μm, which is measurable from the ground to great depth thanks to ALMA. Though there exist a number of works that fit dust temperatures to galaxies beyond z > 2, including SCUBA-2 and South Pole Telescope samples (e.g., Dudzevičiūtė et al. 2020; Reuter et al. 2020), these samples do not have enough galaxies per redshift bin to allow for the analysis presented in the present work to be performed. Investigating whether there is evolution at z > 2 will require millijansky sensitivity in the 250–500 μm wavelength range, which does not currently exist (1 mJy at, e.g., 500 μm corresponds to a galaxy of L_IR = 8.4 × 10¹⁰ L_⊙ at z = 4). Earth's atmosphere is opaque at those wavelengths whenever the atmosphere is not extremely dry. This constraint limits surveys to small areas and, hence, small samples of galaxies. One such survey, reaching a ∼1 mJy depth at 450 μm, is STUDIES (Wang et al. 2017), which covers a 151' region in the COSMOS field and detects too few sources to perform an analysis such as the one presented in this paper. Future space-based FIR missions like the proposed Origins Space Telescope or the 2020 decadal recommendation of an FIR probe could provide these constraints for the z > 2 universe. In the absence of direct constraints on galaxy SEDs down to ∼1 mJy, stacking analyses can provide helpful insight, but the different stacking techniques used in the literature in highly confused FIR maps can result in a broad range of predictions that are not especially constraining (e.g., Viero et al. 2013; Béthermin et al. 2015).

The analysis presented in this paper has shown that disparate techniques in IR/millimeter SED fitting and galaxy sample selection can manifest in very different interpretations of galaxy dust SEDs. Our finding of broad consistency between SEDs from the local universe out to z ∼ 2 was shown only by accounting for both sample selection wavelengths and by approaching SED fitting using a uniform technique. Our SED fitting technique has been tuned to be broadly applicable to most dusty galaxies in the literature that may lack extensive photometric constraints in the IR/millimeter.

Based on work in this paper, we have assembled a list of "best practices" in implementing IR/millimeter SED fitting for the community, which will facilitate easier cross-reference comparisons and eventual measurement of broad physical trends in galaxy populations out to higher redshifts. These recommendations are:

• 1.  
  In most cases, galaxies have a dearth of IR/millimeter data (three to five photometric constraints) and as a result, should not be fit to models with more free parameters than data. Though this is standard scientific practice, it is sometimes violated given the intrinsically complex nature of SED models and the relative dearth of data, especially in this wavelength regime. Most galaxy IR/millimeter SEDs are superpositions of modified blackbodies that may be approximated as a dominant cold blackbody joined with an MIR power law with no more than five free parameters.
• 2.  
  The number of free parameters in IR/millimeter SED fitting should be adjusted according to existing data constraints. Gaussian priors for parameters should be used where data have low signal-to-noise (i.e., S/N < 5), and parameters should be fixed in the absence of data. For example, the MIR power-law slope (when using MCIRSED or other piecewise modified blackbody with an infrared power law) should be set to 2.3 ± 0.5 if fewer than two photometric constraints exist at rest-frame wavelengths λ ≲ 80 μm. Similarly if fixing λ₀ = 200 μm (as we have done here), β should be set to 1.96 ± 0.43 if there are fewer than two photometric constraints longward of rest frame ≳250 μm. These are the median parameter set values for our fitted IRAS sample.
• 3.  
  λ₀ should be fixed if there are no direct constraints on the dust column density via spatially resolved observations used to constrain dust mass surface density.
• 4.  
  Dust temperatures cannot be directly constrained for the vast majority of galaxies that have poorly sampled SEDs and no spatially resolved dust continuum observations. In lieu of dust temperatures, SED fitting should focus on measurement of λ_peak, the rest-frame peak wavelength of the SED in S_ν . Measurements of λ_peak, which serves as a proxy for dust temperature, are directly comparable between works, even when different SED modeling techniques are employed (assuming the dust opacities are roughly comparable). The same is not true of dust temperatures.

We have made the MCIRSED code publicly available ³ to help facilitate the use of these recommendations. While this list should be regarded as a useful guide, it requires custom attention, given the above recommendations. Similarly, galaxies with especially thorough IR/millimeter coverage (e.g., like some local (U)LIRGs, and high-redshift lensed objects and quasars) may not be fit well by our modeling technique, in which case more detailed radiative transfer models may be appropriate. In all instances, we recommend that community members first reflect on their goals in deriving an IR/millimeter SED for a galaxy or sample of galaxies before approaching the task of SED fitting, and use the appropriate tool for the type of measurement that needs to be made.

We choose to fix λ₀ = 200 μm for all galaxies analyzed in this paper for a few reasons. First, our observations—like the vast majority of all galaxies—do not have sufficient wavelength coverage near the peak of dust emission to constrain λ₀ as a free parameter. Second, we do not have resolved dust continuum measurements for these galaxies. Third, fixing λ₀ = 200 μm allows our peak wavelength constraints to be more easily compared with other works in the literature that do the same. Last, our fitting function is, at face value, phenomenological rather than physical. In reality, each galaxy has a distribution of dust temperatures rather than a single "cold" temperature representative of the bulk luminosity of the ISM. A galaxy with a distribution of luminosity-weighted dust temperatures may emit radiation in such a way that mimics opacity effects from the self-absorption factor, (1 − e^−τ ) or a shallower measured dust emissivity index, β. Typical values of λ₀ measured for similarly IR-luminous galaxies or adopted in the literature range from ≈100–200 μm (e.g., Blain et al. 2003; Draine 2006; Conley et al. 2011; Rangwala et al. 2011; Greve et al. 2012; Conroy 2013; Spilker et al. 2016; Simpson et al. 2017; Zavala et al. 2018; Casey et al. 2019).

MCIRSED employs Markov Chain Monte Carlo (MCMC) to construct posterior distributions for our fit parameters and derive confidence intervals. We use the Python package, PyMC3 (Salvatier et al. 2016), with its Hamiltonian Monte Carlo fitter (HMC) and the No-U-Turn Sampler (NUTS; Hoffman & Gelman 2011) to efficiently explore parameter space without requiring the user to tune the HMC parameters. The NUTS algorithm takes advantage of the gradients of the Bayesian likelihoods in order to converge more quickly than other sampling methods for models with high dimensionality. One of the key strengths of using the PyMC3 package with NUTS is that no prior knowledge or testing of step size or number of steps per sample is required for a given galaxy. We refer the reader to Hoffman & Gelman (2011) for the NUTS algorithm and more detail on the tuning phase with the NUTS sampler.

PyMC3 requires the user to specify the model to fit to the data, the prior distributions for each parameter associated with the model, the number of samples to run, and the number of tuning steps to use. Tuning is a phase of fitting done prior to sampling in which the step size and number of steps per sample are automatically varied to minimize the sampling time while trying to reach the target acceptance rate. We choose 3000 tuning steps for all of our fits and exclude those sampling steps from our final analysis. We find that adding more tuning steps does not change the fits for our galaxies. We sample the posterior distributions of each galaxy 5000 times because with this number we find our sampling always converges with >1000 effective samples.

For the IRAS sample, we use flat bounded priors on α, β, and T_d. The bounds are chosen to restrict sampling to parameter space that corresponds to conditions that are physical though not restrictive. The dust temperature may not be lower than the temperature of the cosmic microwave background (CMB) at the redshift of the galaxy being fitted and not hotter than 1200 K, a temperature that corresponds to some of the hottest quasars' host galaxies' ISMs (e.g., Glikman et al. 2006).

The MIR power-law slope, α, is constrained to be between 0 < α < 15. Higher α values asymptotically approach the case of a pure Planck function with α = 15 being nearly indistinguishable. See the left panel of Figure 10 for a plot of example SEDs with a range of α values from 0.5 to 15. The α = 15 curve and the Planck function are indistinguishable.

Zoom In Zoom Out Reset image size

The dust emissivity index, β, is constrained to be between 0.5 < β < 5.5. Some of the highest observed dust emissivity values are β = 3.7 (Kuan et al. 1996; observed in the giant molecular cloud, Sagittarius B2), so we choose a limit beyond that of 5.5. MCIRSED allows the user to set the prior boundaries if one wishes to use higher or lower bounds. For the galaxies analyzed in the present work, we find none of the bounds affect our sampling of the posteriors, as would be evidenced by truncated posterior distributions. See the right panel of Figure 10 for a plot of example SEDs with 0.5 < β < 5.5.

For the H-ATLAS and COSMOS samples, which do not have as high-quality constraints as the IRAS sample, we use the measured average IRAS sample parameters as a basis, adopting normal priors on α and β with medians and standard deviations from the posterior distributions from the IRAS sample. The adopted values are μ_α = 2.3, σ_α = 0.5, μ_β = 1.96, and σ_β = 0.43. The median emissivity, β = 1.96 is fully consistent with several recent publications that found similarly steep values for the emissivity spectral index in various high-redshift galaxies with high-quality FIR constraints (e.g., Kato et al. 2018; Casey et al. 2021). We choose to adopt these priors because the minority of galaxies that produce unphysical fits when run with flat priors have more physical fits using normal priors. For instance, we find that 7.5% of COSMOS galaxies have β that are higher than the highest measured value in the literature of 3.7 (Kuan et al. 1996) when fitted with flat bounded priors. However, the adoption of normal priors results in much lower β values (median 2.03). The reason for the high β values is that the S/Ns of observations that constrain β are low. In the absence of firm constraints, the Gaussian priors return more reasonable values of β, rather than something closer to an average of all allowed values between 0.05 and 5.5. For sources with high-S/N observations on the RJ tail, the best-fit values of β are unaffected by the choice of a flat bounded prior or that of a Gaussian prior taken from the posterior of the IRAS sample fits. Fitting galaxies with Gaussian priors in the case of low-S/N observations on the Rayleigh–Jeans tail, for instance, is less restrictive than fixing β.

Though it is not necessary to account for CMB heating of dust SEDs in this paper (as all sources are at z < 2), we have build MCIRSED to be used for galaxies at any redshift, including those at z > 5, where CMB heating has a nonnegligible effect on the SED. MCIRSED allows users to choose whether to correct for the effects of observing a galaxy against a nonnegligible CMB temperature. By default, the CMB correction is off, but one may toggle it on when calling the MCMC fitting function. Work by da Cunha et al. (2013) demonstrates that if the CMB at the redshift of the observed galaxy is nonnegligible compared with the dust heating due to the interstellar radiation field, corrections need to be made (1) to the observed dust temperature in order to isolate the dust heating due only to the interstellar radiation field, and (2) to the luminosity of the galaxy observed against the nonnegligible background continuum. When this setting is enabled in MCIRSED, in addition to the parameters mentioned above, the ${T}_{\mathrm{dust}}^{z=0}$ parameter (the temperature the dust would have if bathed in the z = 0 CMB) and ${S}_{\nu }^{\mathrm{intr}}$ (the flux density intrinsic to the source without contribution from the CMB) will also be sampled. These quantities are defined as:

$$\begin{eqnarray}&&{S}_{\nu }^{\mathrm{intr}}=\displaystyle \frac{{S}_{\nu }^{\mathrm{obs}}}{1-\tfrac{{B}_{\nu }({T}_{\mathrm{CMB}}(z))}{{B}_{\nu }(T)}},\end{eqnarray} \tag{ A1 }$$

where ${S}_{\nu }^{\mathrm{obs}}$ is the observed flux density, B_ν is the Planck function, T_CMB(z) is the temperature of the CMB at the redshift of the galaxy, and T is the dust temperature measured from the best fit to the data, as above.

$$\begin{eqnarray}&&{T}_{\mathrm{dust}}^{z=0}={\left[{\left(T\right)}^{4+\beta }-{\left({T}_{\mathrm{CMB}}^{z=0}\right)}^{4+\beta }[{\left(1+z\right)}^{4+\beta }-1]\right]}^{\tfrac{1}{4+\beta }},\end{eqnarray} \tag{ A2 }$$

where ${T}_{\mathrm{CMB}}^{z=0}$ is the temperature of the CMB at z = 0, and β is the dust emissidvity index. See da Cunha et al. (2013) for more detail.