The Evolving Interstellar Medium of Star-forming Galaxies, as Traced by Stardust^*

We first consider the "super-deblended" photometric catalog (hereafter SDC1) constructed by Liu et al. (2018) using the FIR and submillimeter images in GOODS-N. These images come from the Herschel Space Observatory (PACS and SPIRE instruments; see Elbaz et al. 2011; Magnelli et al. 2013) and the ground-based facilities SCUBA-2 (850 μm; Geach et al. 2017) and AzTEC+MAMBO (1.16 mm; Penner et al. 2011).

Several novelties introduced in SDC1 are particularly useful for our analysis. First, detections from deep Spitzer IRAC and Very Large Array 20 cm (Owen 2018) are used as a prior for the positions of the blended FIR/submillimeter sources. Second, the SED information from shorter-wavelength photometry is also used as a prior for subtracting lower-redshift sources. This substantially decreases blending degeneracies and allows for a robust photometry extraction of sources at longer wavelengths. Moreover, the authors estimated more realistic photometry uncertainties for each photometric measurement with extensive Monte Carlo simulations. These improvements allow for deeper detection limits and statistically reliable estimates (both measurements and uncertainties) in the FIR+millimeter bands.

The SDC1 contains 3306 priors in total, including over 1000 FIR+millimeter detections. All sources have photometric redshifts and stellar masses inferred by EAZY (Brammer et al. 2008) and FAST (Kriek et al. 2009) respectively, based on the 3D-HST UV-NIR (Skelton et al. 2014) and Pannella et al. 2015 GOODS-N catalogs. Following Liu et al. (2018), we extend the photometric coverage of the published SDC1 catalog to shorter wavelengths by cross-matching with the 3D-HST UV-NIR (Skelton et al. 2014) and Pannella et al. (2015) GOODS-N catalogs. Approximately half of the objects within the catalog are spectroscopically confirmed. However, as mentioned by the authors, the outer perimeter of the GOODS-N area contains objects with high instrumental noise in the 24 μm prior image that may impair the extraction process. We therefore choose to limit our analysis to the central 134 arcmin² with reliable photometry (flag goodArea=1 in SDC1). This reduced our final sample to 2344 objects.

We supplement our study with the "super-deblended" catalog (SDC2 hereafter) presented in Jin et al. (2018). The catalog covers 1.7 deg² in COSMOS, in the same bands as in SDC1, plus additional MAMBO data at 1.2 mm (Bertoldi et al.2007).

In practice, the deblending methodology remains identical to that adopted in SDC1, with one primary difference: an additional step selecting a highly complete sample of priors in K_s band from the UltraVista catalogs (McCracken et al. 2012). The resulting 24 μm detections are then combined with the mass-limited sample of K_s sources in order to fit the remaining bands in the catalog.

The final input data set contains 195,107 priors, with 13% of them having spectroscopic confirmation. The authors highlight that only 11,220 objects are in fact detected over the 100–1200 μm range. Similarly to GOODS-N, we impose the goodArea=1 flag to only include sources that are present in the UltraVista Data Release 4 area (McCracken et al. 2012). We note that for their catalog Jin et al. (2018) used a combination of Laigle et al. (2016) and Muzzin et al. (2013; hereafter M13) catalogs. The M13 catalog has an advantage in that it does not completely remove the sources around saturated stars, which has a positive effect on the number counts; however, the reduced quality of the photometry could lead to unreliable estimates for photometric redshifts (z_phot), as well as any parameters extracted by fitting optical templates. To be on the safe side, our analysis of SDC2 will only focus on the Laigle et al. (2016) sources, which narrows down the number of objects in the input catalog to 186,549 and the total area to 1.38 deg². The COSMOS 2015 catalog (Laigle et al. 2016) also comes with a plethora of UV−optical photometry, spanning an additional ∼20 bands, as well as photometric redshifts and stellar mass estimates by LePhare (Arnouts et al. 1999; Ilbert et al. 2006). We exploit these data by cross-matching the same COSMOS 2015 UV−optical photometry that was used to derive M_* and redshift to SDC2, thus extending our photometric coverage. In total, the merged catalogs consist of ∼40 bands.

For posterity, in Figure 1 we present the UltraVista Data Release 2 map (Laigle et al. 2016; Davidzon et al. 2017) (L16 area), with the regions where the star subtraction took place being clearly identified. In addition to that, on the west border of the survey there exists a number of sources falling outside of the L16 area. These only have UltraVista coverage and lack additional Subaru photometry, which could affect the reliability of the z_phot.

Zoom In Zoom Out Reset image size

The primary parameters that we can derive from observing the rest-frame FIR emission are the total infrared luminosity (L_IR), the dust mass (M_dust), and T_d or, equivalently, the intensity of interstellar radiation field (〈U〉) in the Draine & Li (2007) dust model. To obtain robust estimates for these parameters, an adequate multiwavelength sampling of a galaxy's SED is required. As such, constraining the IR peak is necessary for a robust L_IR estimate, while detections in the long-wavelength regime (R-J) are imperative to capture the emission from cold dust.

With these considerations in mind, after the initial cleaning of the catalogs described above, we perform our sample selection based on the following requirements:

• 1.  
  Detection at a signal-to-noise ratio (S/N) > 3 significance in at least three FIR to submillimeter bands from 100 μm to 1.2 mm. ¹²
• 2.  
  Available z_phot (or z_spec) and M_* estimates inferred by UV to NIR photometry.

After the quality cuts and the selection criteria, we are left with 4331 objects in COSMOS and 585 sources in GOODS-N, which constitute our sample. We have additionally identified 75 objects within SDC2 that fulfill our FIR detection criteria but, despite having a well-sampled UV−optical photometry, do not have either z_phot or M_* estimates. We subsequently fit these sources with EAZY (Brammer et al. 2008), extract their z_phot and other UV−optical properties, and add them to our final catalog. This brings the total amount of COSMOS sources to 4406. The z_phot and z_spec distributions of the final sample in the two fields are presented in Figure 2 (right). As we will discuss later, a third criterion requesting at least one detection at λ_rest > 150 μm will be imposed in order to define a subsample with robust M_dust estimates.

Zoom In Zoom Out Reset image size

Prior to providing the description of our methodology, we believe that it is important to outline and present a brief introduction of the available SED fitting codes that deal with a three-component fitting approach, namely, combining the optical, active galactic nucleus (AGN), and dust emissions. These include the well-established and widely used energy balance routines such as CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019), MAGPHYS (da Cunha et al. 2008; Battisti et al. 2019), and its AGN template extension presented in Chang et al. (2017). These have inspired more novel and sophisticated approaches that optimize the template libraries to achieve significant improvements in computational speeds (SED3FIT; Berta et al. 2013), or adopt Markov Chain Monte Carlo methods when extracting best-fit parameters such as AGNfitter (Calistro Rivera et al. 2016) and Prospector- α (Leja et al. 2018). These efforts are not just limited to published software packages, with many authors implementing their own routines for a panchromatic model analysis (see, e.g., Feltre et al. 2013; Symeonidis et al. 2013).

To model the extensive photometric coverage of the galaxies in our sample, we develop a new, panchromatic SED fitting tool: Stardust. The code performs a multicomponent fit that linearly combines stellar libraries with AGN torus templates and IR models of dust emission arising from star formation (SF-IR). This approach, which is very similar to that presented in Liu et al. (2021), has a number of key differences compared to the currently existing SED fitting codes such as MAGPHYS, CIGALE, and SED3FIT.

First, the three components (stellar, AGN, and SF-IR) are fit simultaneously yet independently from each other, without assuming an energy balance between the absorbed UV/optical radiation and the IR emission. The energy balance approach relies on the assumption that fitted stellar and dust emissions are cospatial, i.e., the process of UV absorption and subsequent reemission at IR wavelength happen in the same environment (da Cunha et al. 2008; see, e.g., their Sections 2.1 and 2.2). However, the detected stellar and dust distributions within a galaxy are not always physically connected. Resolved observations of high-z dusty SFGs (DSFGs; Simpson et al. 2015; Hodge et al. 2016; Franco et al. 2018; Gómez-Guijarro et al. 2018; Hodge et al. 2019; Kaasinen et al. 2020) have revealed spatial offsets between the extent of the dust and stellar emitting sizes of high-z DSFGs, with the former being on average more compact (e.g., Chen et al. 2017; Tadaki et al. 2017; Calistro Rivera et al. 2018; Cochrane et al. 2021). While energy balance is anticipated to apply universally, the aperture and sensitivity limitations of our observations cast a concern on how reliably we can bring these components together in the same panchromatic fit. These observations are also theoretically supported by radiative transfer codes (e.g., SKIRT; Cochrane et al. 2019) and hydrodynamical simulations (e.g., IllustrisTNG; Pillepich et al. 2019; Popping et al. 2021).

Moreover, the detection of "HST-dark" galaxies, i.e., sources that are undetected in the UV/optical bands and thus do not have the photometry to constrain either dust absorption or stellar emission but are bright in the IRAC and FIR/millimeter bands (Wang et al. 2016; Franco et al. 2018), poses another technical challenge to the correct application of the energy balance approach. Finally, the manner in which dust and stellar emission are connected, by assuming a single dust attenuation law and dust composition, can have a significant impact on derived parameters, as these recipes have been shown not to apply universally (Buat et al. 2019). In summary, while the premise of the energy balance routines is undoubtedly elegant and in most cases physically motivated, we choose to use independent stellar, AGN, and dust components to better focus on the dust properties themselves.

Furthering this complex picture, it is important to note the presence of cold diffuse dust, which is being heated by an older stellar population, rather than an ongoing star formation activity (see Boquien et al. 2011; Bendo et al. 2012; Galametz et al. 2014; Hayward et al. 2014). However, when dealing with nonresolved observations, the diffuse dust and R-J tail emissions are highly degenerate, and as such none of the aforementioned codes, nor Stardust, utilize these templates.

The other advantage of Stardust is related to the χ² minimization approach to select the best-fit models. Instead of finding the solution from a vast library of precompiled templates, we devise an optimization method to find the best linear combination of a much smaller set of "basic" templates (similar to eigenvectors in principal component analysis). This is the same approach adopted in the photometric redshift fitting code EAZY (Brammer et al. 2008). In our case, the basic templates are divided into three classes and the linear combination includes a sum of templates from these classes. The models used to create these templates are the following:

• 1.  
  Stellar library. We incorporate an updated version of the stellar population synthesis (SPS) models described in Brammer et al. (2008). Although UV and optical photometry is not always available, the inclusion of the stellar component in the code is important in the NIR regime. In particular, it allows us to better constrain the AGN contribution (see Figure 18 and Appendix A). This stellar library represents an optimized basis set, where the nonnegative linear combinations of models can be considered to be the "principal components" of a much larger parent template catalog (see Brammer et al. 2008; Blanton & Roweis 2007).
• 2.  
  AGN library. We adopt empirically derived templates from Mullaney et al. (2011) describing AGN intrinsic emission from 6 to 100 μm. ¹³ We include both high- and low-luminosity templates (total of 2). Since these can be linearly combined, we do not include the median luminosity AGN template.
• 3.  
  Infrared library. It consists of 4862 Draine & Li (2007, hereafter DL07) templates, with the additional updates from Draine et al. (2014; see also Aniano et al. 2020). These models describe the contribution from warm dust and polycyclic aromatic hydrocarbon (PAH) features in the photodissociation regions (PDRs), together with cold dust in the diffuse part of the ISM. We consider ¹⁴ a wide array of values for the minimum radiation field (${U}_{\min }$) in the 0.1 < ${U}_{\min }$ < 50 range, as well as the fraction of the total dust mass locked in PAHs (q_PAH) between 0% and 10%. We have fixed ${U}_{\max }={10}^{6}$ and α = 2, as described in Magdis et al. (2012a). These templates are not linearly combined within their class; the algorithm instead chooses a single best-fit DL07 template.
• 4.  
  Radio continuum. Data points in radio are not considered by our fitting routine; however, they can be used a posteriori to quantify possible radio excess and further confirm the presence of an AGN, if needed. Our radio model is based on the radio−FIR correlation, described in Delvecchio et al. (2021), with a spectral index of −0.75.

More details on the characteristics of the templates and the motivation for selecting them are provided in Appendix A. With such a configuration, fitting a single object (including the computation of the uncertainties) with Stardust takes less than 10 s, ¹⁵ i.e., a factor of 8–10 faster than software like CIGALE (see Appendix C), based on large precompiled template sets. If we choose to precompile all of our templates, instead of linearly combining them, the resulting model library would contain millions of possible SEDs, with computation time increasing by a significant amount. The code is also highly parallelized, which allows it to run on multiple threads simultaneously, thus achieving significant computation speeds on modern CPUs.

For each object, the input consists of measured flux densities, their corresponding uncertainties, and a redshift estimate. ¹⁶ The user must then choose the corresponding filter curves from the precompiled set, or upload their own. The individual template components can be switched off and on as an additional user input. The algorithm then outputs the best-fit FIR, as well as AGN properties of the source. If the photometry is available, the UV−optical parameters are also produced. These can be summarized as follows:

• 1.  
  The total infrared luminosity integrated over the SF-IR+AGN templates (L_IR,total), ¹⁷ the total infrared luminosity associated with star formation (L_IR,DL07), and the relative contribution of the PDR component (f_PDR) to L_IR,DL07.
• 2.  
  The bolometric IR luminosity of the AGN (L_AGN) and its fractional contribution (f_AGN) to the total IR energy budget. ¹⁸
• 3.  
  The total dust mass (M_dust), the warm dust mass component heated by PDRs (${M}_{\mathrm{dust}}^{\mathrm{warm}}$), the fraction of the total dust mass heated by PDRs (γ), the cold dust mass component (${M}_{\mathrm{dust}}^{\mathrm{cold}}$) in the diffuse ISM, and the fraction of the total dust mass locked in PAHs q_PAH.
• 4.  
  The intensity of the radiation field that the diffuse ISM is exposed to (${U}_{\min }$), which is a proxy of the mass-weighted T_d of the ISM, and the mean radiation field intensity (〈U〉), which is a proxy of the luminosity-weighted T_d.
• 5.  
  The stellar mass (M_*), star formation history (SFH), E(B − V), and the unobscured SFR, if there is available optical photometry.

Figure 3 presents an example fit to one of the COSMOS galaxies, chosen for its prominent AGN contribution. The top panel of the figure shows the data points and the best-fit model, with different colors for the four components listed above; the bottom panel displays the χ² distributions of all relevant IR quantities.

Zoom In Zoom Out Reset image size

In order to estimate the errors associated with the derived quantities during the fitting procedure of Stardust, we consider two main sources of uncertainty: one concerning the linear combination coefficients of the best-fit optimization, and one linked to the broadband photometric data.

To quantify the linear combination uncertainty, we resample the best solution coefficients. A covariance matrix is first created by considering all of the templates that went into the best-fit solution. We draw the coefficients from a multivariate normal distribution whose median is given by the coefficient of the best solution vector and the standard deviation is computed from the diagonalized covariance matrix. This is done 10⁴ times to provide a good balance between robustness of the error estimates and computational speed. We recompute all the relevant FIR properties for each realization of our routine. From resultant distributions we then define our lower and upper uncertainty as the 16th and 84th percentile confidence interval, respectively. However, given the fact that only a single solution with a single coefficient is considered in the IR, the final uncertainty on the FIR properties is underestimated.

We therefore also consider the observational uncertainty, which is primarily driven by the quality of the photometric data. We compute it by considering all possible solutions from our template library that fall within the 68% confidence interval range of the best-fit. This would correspond to a region where the solutions fall within ${\rm{\Delta }}{\chi }_{\nu }^{2}=1$, where ${\rm{\Delta }}{\chi }_{\nu }^{2}={\chi }_{\nu ,i}^{2}-{\chi }_{\nu ,\mathrm{best}}^{2}$, since the nondiagonal terms of the template covariance matrix are not zero. ¹⁹ We show these as red shaded areas on the bottom panel of Figure 3. The observational uncertainty is then derived as simply the width of the shaded region.

The final errors are computed as a quadrature sum of the linear combination uncertainty and the observational uncertainty.

To explore and quantify how the uncertainty in z_phot propagates into the estimates of L_IR and M_dust, we built mock IR SEDs of 1200 galaxies utilizing a suite of $0.1\lt {U}_{\min }\lt 50$ DL07 models and place them uniformly in the 0.03 < z_true < 5 redshift range. Thus, each mock SED is characterized by a set of predefined L_IR,in, M_dust,in, and z_true values. We then infer synthetic broadband photometry in all IR bands available in SDC2 (24–1100 μm) for each simulated galaxy. Since we are interested in the effect of the z_phot uncertainty on the derived FIR properties, to minimize any other possible sources of error (e.g., photometric uncertainty, poor photometric sampling of the SED), we adopt S/N = 5 in all bands and add to our photometric data set the monochromatic flux density of the template at 2.2 mm. We then fit the synthetic photometry of each galaxy by fixing the redshift first to z₋ = z_true − Δz and then to z₊ = z_true + Δz, where Δz = (1 + z_true). For our purposes, and based on the z_phot accuracy of the COSMOS field ( = 0.005–0.03 as quantified in L16), we first adopt = 0.02 and then repeat the analysis for an even more conservative case with = 0.05. The comparison between the extracted L_IR,out and M_dust,out to the input values for each simulation is presented in Figure 4.

Zoom In Zoom Out Reset image size

Our analysis suggests that the effect Δz has in the derivation of the FIR properties is not negligible, even for the idealized case of detailed (24 μm–2.2 mm) and high-quality (S/N = 5) photometric coverage. Indeed, we find that a typical value of = 0.02 ( = 0.05) introduces an extra scatter of ∼12% (25%) and ∼17% (35%) in the derived L_IR and M_dust, which remains rather constant with redshift (at least out to z = 4). At the same time, we also find that a symmetric Δz, such as the one adopted in our simulations, does not inflict a noticeable systematic offset in the extracted FIR quantities (<0.05 dex).

Based on these results, we update the uncertainties of the inferred FIR properties of our z_phot sample (ztype = 0; see Appendix A) by adding in quadrature the extra error arising from a symmetric Δz (assuming = 0.02 for all sources) to the error budget inferred by the SED fitting procedure (photometry+model). While our correction is based on an average value of = 0.02, we note that for catastrophic z_phot failures ( ∼ 0.15) we find a systematic offset of ≤ 30% between the input versus output M_dust, while the L_IR ratios remain uniformly scattered around unity.

It has been well established that for a robust modeling of the M_dust at least one photometric data point longward of the FIR peak, i.e., into the R-J tail of the SED, is necessary (e.g., Draine et al. 2007; Magdis et al. 2012a; Berta et al. 2016; Scoville et al. 2017). Here we attempt to quantify how the accuracy of the derived M_dust estimates varies as a function of the rest-frame wavelength of the last available detection in conjunction with the selection criteria of our sample (i.e., at least three detections at 24 μm < λ ≤ 1200 μm). For this, we perform the following simulations.

We start by building mock IR SEDs of fixed M_dust and L_IR, with 0.1 $\lt {U}_{\min }\leqslant 50$, 0 ≤ γ ≤ 0.5, and q_PAH = 2.8%, using the DL07 library and place them at 0.01 ≤ z ≤ 4.5 with a step of Δz = 0.05. After all of the models are created, synthetic photometry is performed by convolving each mock redshifted SED with a filter transmission curve. We consider all filters redder than MIPS 24 μm available in SDC2, for a total of nine bands, and set all recovered fluxes to an S/N = 3 significance level. At each redshift the algorithm calculates the rest-frame wavelength for each available band, producing a grid of possible rest-frame wavelengths of the last detection (λ_last) after accounting for our selection criterion that requests the availability (and the detection) in at least two bluer bands. For each λ_last it then randomly selects two additional bands blueward of λ_last, producing a final set of three photometric data points. We then fit each set, as well as 50 permutations by varying the original fluxes by 10%, with our code to extract M_dust and L_IR estimates. This procedure is then repeated for all mock SEDs and all acceptable values of λ_last in each redshift.

The accuracy with which we can recover M_dust estimates for each λ_last is then quantified by the scatter of the ratio of output to input M_dust presented in Figure 5. As expected, we find a decreasing scatter in M_dust,out/M_dust,in at longer λ_last, which drops from a factor of ∼2 (for 68% of the simulated galaxies) at λ_last = 100 μm to a factor of ∼1.1 at λ_last = 400 μm.

Zoom In Zoom Out Reset image size

Based on these results, we choose to define the subsample of "M_dust–robust" galaxies, for which at least one detection at λ_last ≥ 150 μm is available. This threshold was chosen as an optimal compromise between the number of the rejected sources and the precision of the derived M_dust, which for λ_last ≥ 150 μm is ≥70%. Indeed, while λ_last ≥ 150 μm is evidently not deep into the R-J, it seems that the addition of the two extra data points blueward of λ_rest = 150 μm is adequate to anchor the general shape, and eventually the M_dust, of the templates.

As a sanity check for the effectiveness of our criterion, we cross-match the "M_dust–robust" sample with the A3COSMOS ALMA photometric catalog presented in Liu et al. (2019a) and refit the 233 galaxies that we find in common, adding this time the extra ALMA data point to the input photometry. The comparison of the inferred M_dust,A3 to our M_dust estimates yields a very good agreement between the two estimates with a median log(M_dust,A3/M_dust) ≈ −0.04 ± 0.06, further supporting our analysis. Nevertheless, we do identify a handful of sources for which the addition of the ALMA data point results in lower M_dust estimates by a factor of ≥3. An inspection of the SEDs of these extreme outliers reveals that the discrepancy originates either from possible catastrophic blending of the SPIRE 500 μm photometry or alternatively from overresolved ALMA photometry. ²⁰

The emerging "M_dust–robust" sample consists of 3312 sources drawn from the same M_* and redshift distributions as the 4991 originally selected galaxies. Finally, we note that our simulations operate under the assumption that the DL07 models are a good representation for the FIR emission of the real galaxies. Variations in dust composition, dust emissivity, and dust absorption coefficients that could result in systematic offsets in the inferred M_dust (e.g., Dale et al. 2012; Magdis et al.2012a; Berta et al. 2016; Scoville et al. 2017) cannot be modeled with our approach. As is the case for any other M_dust analysis, the relative rather than the face-value estimates bare more physical significance.

We now attempt to compute the completeness threshold of our sample in terms of M_dust and L_IR as a function of redshift. Again, we build a grid of mock IR SED in the z = 0−5 range using the same approach as described above. However, instead of considering the full range of possible 〈U〉 ∝ L_IR/M_dust values of the DL07 models, the constructed templates this time follow the 〈U〉–z relation of MS galaxies presented in Béthermin et al. (2015). At each redshift we then create a grid of templates normalized to log(M_dust/M_⊙) = 6−10 in steps of 0.1.

As before, the templates are used to produce synthetic photometry for each template in all bands available in the SDC2 catalog. For each band we then adopt an rms based on the depth of the corresponding survey at each wavelength in the COSMOS field (Jin et al. 2018) and impose the same selection criteria on the mock photometry as those applied to the real catalogs. Following Section 4.1, we also request that the simulated sources have λ_last ≥ 150 μm. If a galaxy of given M_dust fulfills these criteria at a given redshift, the algorithm moves to a lower M_dust until the object is rejected by our selection. By following these steps at different redshifts, we thus obtain a limiting M_dust as a function of z, which we coin lim(M_dust)(z). In the process we also consider the scatter of the 〈U〉–z relation of Béthermin et al. (2015), in order to account for the variation of 〈U〉 among MS galaxies at a given redshift. The derived lim(M_dust)(z) can then be converted to lim(L_IR)(z), via

$$\begin{eqnarray}&&\langle U\rangle =\displaystyle \frac{{L}_{\mathrm{IR},\mathrm{DL}07}}{125\,{M}_{\mathrm{dust}}},\end{eqnarray} \tag{ 1 }$$

as described in Draine & Li (2007).

We also repeat our simulation for SBs, by fixing 〈U〉 = 40 (e.g., Magdis et al. 2012a; Tan et al. 2014; Béthermin et al. 2015). We note, however, that the 〈U〉 of SBs can vary substantially to lower or higher values (e.g., Magdis et al. 2012a, 2017; Schreiber et al. 2018; Jin et al. 2019; Cortzen et al. 2020). Therefore, the chosen 〈U〉 = 40 is only representative, but not necessarily unique. Nevertheless, 〈U〉 < 40 templates are represented in the simulation of the MS galaxies, while galaxies with 〈U〉 > 40 are rather rare.

The results of our simulations are presented in Figure 6, where we show the derived M_dust and L_IR as a function of z for the whole sample, along with the evolution of lim(M_dust) and lim(L_IR). We see that the limiting M_dust increases toward high z, peaking at z ∼ 2, and remains flat afterward, signifying that the balance between cosmological dimming and negative K-correction is achieved beyond that point. Following the black line, we could infer the M_dust threshold above which our sample should be 100% complete, assuming an MS-like population of galaxies. However, since our sample is not limited to MS galaxies, we naturally also find sources that fall below our limiting M_dust track. As we will discuss later, these are predominately starbursting galaxies that exhibit elevated 〈U〉 with respect to the MS. Our secondary lim(M_dust) trend for an SB-like population displays that with the SDC2 detection limits it is possible to detect a low-M_dust object only if it is also an SB.

Zoom In Zoom Out Reset image size

We also find a similar trend for MS galaxies when considering the evolution of lim(L_IR). In this case, however, we do not observe a flattening at z ∼ 2, and the trend continues to rise into the early universe. The balance between cosmological dimming and the negative K-correction is not being achieved here, since the wavelengths that are required to reliably constrain the L_IR are positioned to the left of the FIR peak.

Admittedly, depths of FIR surveys are not the only limiting factors of sample selection. A requirement to have a photometric redshift and M_* would mean that the optical photometry has to be sufficiently sampled, to allow such an analysis. Moreover, the deblending procedure itself goes through various selection stages, including both brightness and mass cuts. Various IR studies (Wang et al. 2016; Franco et al. 2018) have revealed substantial populations of optically dark sources at 2 < z < 4, which are otherwise bright in IRAC and FIR bands, which would be unintentionally excluded from our analysis. Moreover, even at low z, objects that are faint in the K band would also be missed. Indeed, a combination of these factors creates significant obstacles in our completeness analysis; we address this in more detail in Section 6.2.

We now focus on the FIR properties of the "M_dust–robust" galaxies described in Section 4.1, which should represent the most reliable sample for the exploration of the dependency of the L_IR, M_dust, and 〈U〉 on redshift, cosmic age, and ΔMS. The emerging results are presented in Figure 9, where for completeness and to facilitate comparisons we also include the inferred properties of the full sample.

Zoom In Zoom Out Reset image size

Both L_IR and M_dust are found to increase smoothly as a function of ΔMS. At the same time, we also find that for MS galaxies 〈U〉 evolves as $(3.2\pm 1.3)\times {\left(1+z\right)}^{1.2\pm 0.3}$, which is in excellent agreement with the stacking analysis of Béthermin et al. (2015). The fact that our individually detected galaxies appear to follow the same 〈U〉–z relation as the much deeper stacked ensembles reinforces the notion that the adopted "M_dust–robust" subsample does not introduce a significant bias toward colder objects.

Since 〈U〉 is proportional to L_IR,DL07/M_dust, and also a proxy for T_d, our analysis provides further evidence that dust in MS galaxies becomes warmer toward higher redshifts, a trend that has already been recovered in previous studies (most of them based on stacking analysis; see, e.g., Magnelli et al. 2014; Schreiber et al. 2015; Davidzon et al. 2018). Similarly, our data also confirm a progressive increase of 〈U〉 (or T_dust) with an increasing elevation above the MS (e.g., Magdis et al. 2017; Jin et al. 2019).

It is worth noticing that the full sample follows the same general trends albeit with a considerably larger scatter (∼×2) in M_dust and 〈U〉. The reduced scatter for the "M_dust–robust" subsample is driven by the imposed λ_last ≥ 150 μm selection criterion that primarily removes the locus of sources with very cold fitting solutions (〈U〉 ≲ 1). We highlight that the rejection of these objects should not introduce a bias in our sample since such low 〈U〉 values are more indicative of poor photometric coverage/quality (lack of available data point in the R-J) rather than of realistic, extremely cold ISM conditions. However, we note that not all of the extremely cold solutions have been removed from the "M_dust–robust" sample by our selection, as ∼200 objects with 〈U〉 < 1 meet our λ_last > 150 μm criterion. These can be easily identified in the 〈U〉 − z plot and as the outliers populating the secondary blue cloud of points in the M_dust−ΔMS plot in Figure 9. As we will discuss later, these could be sources with unreliable z_phot estimates, failures of the deblending in the SPIRE bands, or, more interestingly, gas giants or very compact galaxies with optically thick FIR emission.

The SED decomposition introduced in Section 3 can also provide constraints on the relative contribution of the warm (PDR, ${L}_{\mathrm{IR}}^{\mathrm{warm}}$, ${M}_{\mathrm{dust}}^{\mathrm{warm}}$) and cold (diffuse, ${L}_{\mathrm{IR}}^{\mathrm{cold}}$, ${M}_{\mathrm{dust}}^{\mathrm{cold}}$) ISM components to the total L_IR output and the total M_dust budget of the galaxies in our sample. In particular, it is worth investigating if and how the relative contribution of the components varies as a function of ΔMS. Indeed, if SBs are experiencing elevated star formation activity per surface area (Elbaz et al. 2011, 2018; Valentino et al. 2020), then one would expect to see an increased fraction of L_IR (and M_dust) originating (and being heated) from the "PDR" component, where the radiation intensity ranges from ${U}_{\min }$ to ${U}_{\max }$ (Draine & Li 2007).

In Figure 10 we plot the inferred properties of the warm and cold ISM components as a function of ΔMS. We find that for a fixed L_IR (or equally SFR), SBs tend to have lower amounts of ${M}_{\mathrm{dust}}^{\mathrm{cold}}$, with ${M}_{\mathrm{dust}}^{\mathrm{cold}}$/L_IR showing a tight anticorrelation with ΔMS. The same, however, does not apply to ${M}_{\mathrm{dust}}^{\mathrm{warm}}$/L_IR, which exhibits a significantly larger scatter and only a very weak dependence on ΔMS. This is a consequence of an increasing fraction of warm to cold M_dust and L_IR between MS and SB galaxies (Figure 10).

Zoom In Zoom Out Reset image size

The observed trends suggest that, compared to MS galaxies, SBs have a larger fraction of their total M_dust exposed to the intense stellar radiation fields of the PDRs, in agreement with expectations discussed above. Our result could indeed reflect an increase in the compactness of the star formation activity for increasing distance from the MS as suggested in recent high-resolution studies. Finally, under the assumption that ${M}_{\mathrm{dust}}^{\mathrm{cold}}$ is proportional to M_gas and L_IR is proportional to SFR, our results point to enhanced star formation efficiencies and shorter gas depletion timescales for sources residing above the MS, as already reported in the literature (e.g., Daddi et al. 2010; Tacconi et al. 2010, 2020; Magdis et al. 2012a, 2017; Sargent et al. 2014; Silverman et al. 2018).

To infer the evolution of f_dust, we adopted the formula described in Liu et al. (2019b), which parametrizes f_dust in terms of z, M_*, and ΔMS. Compared to more simple log-space linear fitting models (e.g., Scoville et al. 2017), this formulation recovers trends that are more physically meaningful and also explores how these parameters are covariant and degenerate with each other in a multidimensional fitting space. As an initial check, we performed a Spearman correlation test and found M_dust to be mildly correlated with logΔMS (ρ = 0.40) and strongly correlated with M_* and t_age (ρ = 0.63 and −0.80, respectively). We consider the following functional form:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}({f}_{\mathrm{dust}})=({a}_{0}+{a}_{1}\,\mathrm{log}({M}_{* }/{10}^{10}{M}_{\odot }))\mathrm{log}{\rm{\Delta }}\mathrm{MS}\\ \,+\,b\,\mathrm{log}({M}_{* }/{10}^{10}{M}_{\odot }))\\ \,+\,({c}_{0}+{c}_{1}\,\mathrm{log}({M}_{* }/{10}^{10}{M}_{\odot })){t}_{\mathrm{age}}(z)\\ \,+\,d,\end{array}\end{eqnarray} \tag{ 2 }$$

where t_age is the cosmic age at a given redshift in Gyr and M_* is the stellar mass in M_⊙. For fitting we use the Python package scipy.optimize.curve_fit, which finds the solution based on the least-squares method. We also consider how the extreme outliers can affect our results and thus only fit the medians in a given redshift bin. The best-fit values are as follows:

$$\begin{eqnarray*}\begin{array}{cc}{a}_{0}=-0.90 & \quad {a}_{1}=1.31\\ b=-0.98 & \quad {c}_{0}=-0.23\\ {c}_{1}=0.11 & \quad d=-0.64,\end{array}\end{eqnarray*}$$

with the uncertainties being computed from the covariance matrix. We then used the functional form given by Equation (2) to renormalize all galaxies to lie on the MS (ΔMS = 1) and M_* = 5 × 10¹⁰ M_⊙, in order to directly compare with our best-fit function in two dimensions.

The normalized data and the best-fit relation presented in Figure 11 are in very good agreement with the collection of similar trends drawn from the literature (e.g., Scoville et al. 2017; Tacconi et al. 2018; Liu et al. 2019b; Magnelli et al. 2020). We note that any apparent discrepancies between the slope and the normalization of our recovered relation to that of Scoville et al. (2017) and Magnelli et al. (2020) are driven by the model parameterization, as for the latter the multivariable functional forms do not consider the covariance between the fitted values.

Zoom In Zoom Out Reset image size

In Figure 12 we also show how the one-dimensional relation between M_dust and M_* compares to the multidimensional fit within six redshift bins. Since our analysis can be affected by the completeness of our sample in terms of M_* and M_dust, we also consider the underlying selection effects based on our lim(M_dust)(z) derivation and assuming that our catalog is complete at M_* > 10¹⁰ M_⊙. We find that for a fixed redshift range the inferred M_dust–M_* relation shadows the multiparameter fit up to z ∼ 1.1, in moderate agreement with the trends reported in Liu et al. (2019b) and Magnelli et al. (2020). However, the increasing incompleteness and the low number statistics do not allow us to extend our analysis at higher redshifts.

Zoom In Zoom Out Reset image size

Based on these results, we will limit the subsequent DMF analysis to the 0.2 < z < 1.1 range.

After obtaining the functional form for f_dust(t_age, M_*, ΔMS), we are now in a position to examine the shape and the evolution of the DMF with redshift. For this we will restrict our analysis to the COSMOS sample (SDC2), as we do not have enough statistics or coverage of the GOODS-N field to reliably constrain the DMF.

The "super-deblending" procedure that went into producing the SDC2 catalog creates a significant obstacle when attempting to consider all of the incompleteness effects of the sample. The objects that end up in our sample go through several selection stages, both before and after the deblending procedure. These include both the brightness and mass cuts of the parent catalog (Laigle et al. 2016), the availability of infrared coverage (Jin et al. 2018), and the selection criteria imposed in our study. As such, it is not possible to robustly assess the properties and the number of objects that end up being "lost." Therefore, we select an alternative approach in computing the DMF for our objects, namely, utilizing the derived f_dust parameterization along with the available SMFs in the literature. Later in this section, we also attempt to account for the incompleteness effects and compare the SMF-derived DMF to the observed number density of galaxies per M_dust bin.

6.2.1. DMF from SMF

For our analysis we adopt the SMF computed by Davidzon et al. (2017), which covers the entire COSMOS field and the z = 0.2−4 range. Their mass function already accounts for the Eddington bias, so we do not need to consider any additional corrections. Since the vast majority of galaxies in our sample are star-forming, we adopt the derived parameters for the "active" SMF only.

The galaxy mass function is normally expressed as a Schechter function (Schechter 1976), which in logarithmic form can be written as

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}(\mathrm{log}M)d(\mathrm{log}M)=\mathrm{ln}(10)\,\exp (-M/{M}^{* })\\ \quad \times \,{{\rm{\Phi }}}^{* }{\left(-M/{M}^{* }\right)}^{(\alpha +1)}d(\mathrm{log}{M}_{* }),\end{array}\end{eqnarray} \tag{ 3 }$$

where α is the slope of the faint end, Φ^* is the normalization, and M^* is the characteristic mass, indicating the position of the "knee" of the Schechter function. To convert the SMF to DMF, we first postulate that the number of galaxies in a given redshift bin is the same, regardless of whether we integrate over M_* or M_dust, namely,

$$\begin{eqnarray*}&&{\int }_{0}^{\infty }{\rm{\Phi }}(\mathrm{log}{M}_{* })d(\mathrm{log}{M}_{* })={\int }_{0}^{\infty }{\rm{\Phi }}(\mathrm{log}{M}_{\mathrm{dust}})d(\mathrm{log}{M}_{\mathrm{dust}}),\end{eqnarray*}$$

and the integrands can be rearranged to obtain

$$\begin{eqnarray*}&&{\rm{\Phi }}(\mathrm{log}{M}_{\mathrm{dust}})={\rm{\Phi }}(\mathrm{log}{M}_{* })d(\mathrm{log}{M}_{* })/d(\mathrm{log}{M}_{\mathrm{dust}}).\end{eqnarray*}$$

We then differentiate Equation (2), to obtain $d(\mathrm{log}{M}_{* })/d(\mathrm{log}{M}_{\mathrm{dust}})={\left(b+{c}_{1}\times {t}_{\mathrm{age}}+1\right)}^{-1}$. In order to perform the final conversion, we also transform all the M_* bins into M_dust bins, by inverting our formulation of f_dust, taken at ΔMS = 1.

6.2.2. Accounting for the Eddington Bias

6.2.3. Comparison to the Observed Number Density

Now we would like to compare the SMF-derived DMF to the observed number density of galaxies in the "M_dust–robust" sample. To this end, we first apply the widely used nonparametric 1/${V}_{\max }$ method to correct for the Malmquist bias of our sample (Schmidt 1968). This method relies on assigning the ${V}_{\max }$ to each redshift bin, based on the detection limits of the survey. Effectively, this correction accounts for the fact that in a given volume a brightness-limited survey is more likely to pick up the bright sources, while missing faint galaxies, which would populate the low-M_dust end. We explicitly highlight that the ${V}_{\max }$ correction only accounts for the FIR flux rms cuts and not the selection criteria outlined in Section 2.3.

To calculate the ${V}_{\max }$, we use the prescription from Weigel et al. (2016), which provides a volume correction for each individual source. As a first step, and for a given redshift bin, we split our sources into 0.4 dex bins in the log(M_dust/M_⊙) = 6–11 interval. Given that the median uncertainty on the M_dust is ∼0.3 dex, the following bin spacing will ensure that there is very little to no cross-contamination between mass bins. The ${V}_{\max },i$, where i denotes an individual galaxy, can be then calculated as

$$\begin{eqnarray}&&{V}_{\max }=\displaystyle \frac{A}{3}\,({d}_{c}{\left({z}_{\max ,{\rm{i}}}\right)}^{3}-{d}_{c}{\left({z}_{\min ,{\rm{i}}}\right)}^{3}),\end{eqnarray} \tag{ 4 }$$

where d_c is the comoving distance and A is the area, which in our case is equal to 1.38 deg². Following Weigel et al. (2016), the ${z}_{\min ,{\rm{i}}}$ is given simply by the lower boundary of the bin. The ${z}_{\max ,{\rm{i}}}$, on the other hand, can be calculated empirically, either through detection limits of individual bands or by considering a limiting mass of the survey. It, however, cannot exceed the maximum redshift of the bin.

To obtain the ${z}_{\max }$, we consider the best-fit SEDs for our sources and the rms of the parent catalogs in order to redshift the sources to the point where they no longer fulfill our selection criteria as outlined in Section 2. Using this method, we, however, found that, for an overwhelming majority of sources, the computed ${z}_{\max ,{\rm{i}}}$ exceeds the upper boundary of the bin they belong to. Therefore, the ${V}_{\max }$ correction that we apply becomes effectively bound between the lower and the upper redshift of the bin. We find that this method works best in the lowest (0.2 < z < 0.5) redshift bin, with the $\langle V/{V}_{\max }\rangle$ test returning a value of 0.47. The remaining two redshift bins are significantly incomplete, with the ratio returning 0.83 and 1.38, respectively.

Among the other sources of incompleteness, as discussed in Section 6.2, here we can attempt to account for lost sources due to the sensitivity limits of the survey and failures in the deblending procedure. We therefore multiply our points by the loss fraction in each redshift bin that is computed as the ratio between sources in our catalog over the sources that have SN-IR > 1 in the parent catalog. ²¹ The SN-IR parameter, described in greater detail in Liu et al. (2018) and Jin et al. (2018) and references therein), considers a combination of FIR bands starting with 100 μm. We thus expect that in this context our SN-IR threshold can indicate whether a galaxy is intrinsically dusty.

The comparison of the SMF-derived DMF with and without the Eddington bias taken into account, along with the observed volume-weighted number density of galaxies, derived as described above, for three redshift bins, is presented in Figure 13. We find that the DMF without the Eddington bias is insufficient to reproduce the observed dynamic range of M_dust, particularly in the higher-redshift bins, exactly as we have predicted in Section 6.2.2. On the other hand, the SMF-derived DMF with the artificially induced bias, through the f_dust scatter, is in good agreement with the data in the high-mass end, further highlighting the necessity of accounting for the observational biases when inferring relations (i.e., M_*–M_dust in our case) through the observed mass (or luminosity) functions. Although our model underpredicts the high-mass end data in the higher-redshift bin, both still agree within the error bars. At the same time, though, in the low-mass regime our data significantly underestimate the number density of galaxies mirroring the incompleteness of our sample in this M_dust regime. It is worth noticing, though, that the turnover of the observed data perfectly coincides with the independent estimates of lim(M_dust), offering an indirect validation of our simulations presented in Section 4.2. For the analysis in the next sections, we adopt the SMF-derived DMF, which has the Eddington bias corrected out, as the final result against which we will compare previous observationally driven DMFs and theoretical predictions. ²²

Zoom In Zoom Out Reset image size

Similarly to f_dust, we also explore the dependence of f_gas on cosmic age, ΔMS, and M_*. We utilize the same multiparameter fitting function as before (see Equation (2)) and focus on the "M_dust–robust" sample. We calculate the Spearman rank correlation between our variables and find M_gas to be mildly correlated with logΔMS and M_* (ρ = 0.47 and 0.53, respectively) and strongly negatively correlated with t_age (ρ = −0.82). The fitting procedure yields the following best-fit parameters:

$$\begin{eqnarray*}\begin{array}{cc}{a}_{0}=-0.73 & {a}_{1}=1.16\\ b=-1.02 & {c}_{0}=-0.20\\ {c}_{1}=0.09 & d=1.39.\end{array}\end{eqnarray*}$$

The best-fit f_dust–t_age (or redshift) relation and our data, both normalized to ΔMS = 1 and M_* = 5 ×10¹⁰ M_⊙, are presented in Figure 14 (top). Similar to previous studies (Scoville et al. 2017; Tacconi et al. 2018; Liu et al. 2019b; Magnelli et al. 2020), we find a sharp increase in the gas fraction up to z = 2, followed by a milder evolution at higher redshifts, a change that is noticeable only in the f_dust–z parameter space. We note, however, that due to poor statistics and lack of spectroscopic redshifts, our data cannot reliably constrain the high-z evolution of f_gas at z > 2. We also detect a population of sources that display significantly elevated gas reservoirs for their redshift (log(f_gas) > 0.5). Some of those objects either have only a z_phot estimate available or appear to be blended in the SPIRE bands and therefore could have an erroneously large M_dust and subsequently M_gas estimate assigned to them. However, among these, we do identify some sources with spectroscopic redshifts, which are also "clean"/isolated in the IR maps. In particular, we find ∼40 such sources with log(f_gas) > 0.5) and six with log(f_gas) > 1 that, as we discuss later, we coin as "gas giants."

Zoom In Zoom Out Reset image size

Finally, we focus our attention on the depletion time τ_depl = M_gas/SFR = 1/SFE. We employ the same fitting technique as before and explore the evolution of τ_depl with cosmic age, ΔMS, and M_* in a multidimensional parameter space, for the "M_dust–robust" sample. The τ_depl correlates mildly with age and ΔMS (Spearman ρ =−0.46 and −0.58, respectively) and weakly with M_* (ρ = −0.23). The best-fit parameters are as follows:

$$\begin{eqnarray*}\begin{array}{cc}{a}_{0}=-1.65 & {a}_{1}=1.42\\ b=-0.95 & {c}_{0}=-0.02\\ {c}_{1}=0.12 & d=0.37.\end{array}\end{eqnarray*}$$

We show the best-fit τ_depl−t_age relation and our data, both normalized to ΔMS = 1 and M_* = 5 × 10¹⁰ M_⊙, in Figure 14 (bottom). In line with the previous studies of τ_depl by Scoville et al. (2017), Tacconi et al. (2018), and Liu et al. (2019b), we recover a relatively weak decrease of depletion time (or increase in SFE) with redshift.

The recovered trends between M_dust, M_*, and M_gas and their evolution with redshift offer a test bed against which theoretical and previous observationally driven studies can be compared to. As shown in Figures 11 and 14, our analysis yields f_dust and f_gas evolutionary tracks consistent with those presented in Tacconi et al. (2018), Liu et al. (2019b), and to a smaller degree those reported in Scoville et al. (2017) and Magnelli et al. (2020). As discussed earlier, the mild tension between the latter works and our results could be primarily attributed to the choice of the fitting function.

The f_dust and f_gas in our sample of SFGs increase rapidly from z = 0 to z = 1, peak around z ∼ 2−3, and then remain roughly constant. It is, however, a point of contention whether the latter is driven by actual physical processes or is a consequence of the scarcity of data at z > 2. It it also worth mentioning that our analysis points toward a milder evolution of f_dust (−0.8 dex) compared to f_gas (−1.3 dex) from z = 2 to z = 0, with the latter dropping ∼ 3 × faster. This is aligned with the evolution of ρ_dust and ρ_gas derived by the ALMA stacking analysis of Magnelli et al. (2020) and could in fact reflect the evolution of metallicity, and thus of δ_GDR, for fixed M_* toward lower redshifts.

At the same time, the decrease of M_dust with decreasing redshift, for fixed M_*, can be attributed to either the destruction of dust grains by interstellar radiation fields or their incorporation into the stellar population. This is discussed in more detail in Donevski et al. (2020), who also report a decreasing f_dust from earlier cosmic age to the present epoch. We note that the observed trend could also mirror the overall decline in the SFR density in the universe from z = 2 to the present day, which points toward lower star formation activity and thus lower dust production at later cosmic epochs. Finally, at a fixed redshift, both f_dust and f_gas decrease as a function of M_* (as indicated by a negative value of the fitting parameter b; see Equation (2)), in line with previous studies (e.g., Magdis et al. 2012a; Magnelli et al. 2020 and references therein).

In addition to observational studies, we can also compare our results to theoretical predictions. To this end, we consider the HAGN hydrodynamical simulations in the COSMOS field (Dubois et al. 2014; Laigle et al. 2019) and draw a sample of SFGs (ΔMS > 0.3) in the z = 0.2–0.5 range, selected to meet the M_* completeness of the COSMOS 2015 survey, and which fall within a simulation box of 143 Mpc per side (Dubois et al. 2014). To measure ΔMS for each galaxy, we considered the M_* and the 100 Myr averaged SFR from the simulations. Also, since M_dust is not an explicit parameter of HAGN galaxies, we used a constant δ_GDR = 100 to convert the M_gas values, as derived from the simulations, to M_dust. We then use the M_gas, M_*, and M_dust of the simulated galaxies to infer f_gas and f_dust. The median values and their scatter, renormalized to MS (ΔMS = 1) and M_* = 5 × 10¹⁰ M_⊙ in four redshift bins, are presented and compared to the real data in Figures 11 and 14. We find a good agreement between the theoretical predictions and our observationally driven trends (in the 0.2 < z < 0.5 range at least), indicating that the HAGN simulation can successfully reproduce the baryonic components of the galaxies and its evolution with redshift. Conversely, the agreement of our results with both theoretical and observational studies provides an extra indirect validation for the performance of our new SED fitting code.

As with f_dust and f_gas, our recovered trends, which connect τ_depl to redshift, ΔMS, and M_*, show similar behavior to the ones presented in Scoville et al. (2017), and to a lesser extent those in Tacconi et al. (2018) and Liu et al. (2019b). The dependence of τ_depl on M_* is relatively weak across all studies; however, similarly to Liu et al., we find that the depletion time for high-mass galaxies increases from early cosmic ages toward present times, while low-mass galaxies display an opposite trend of decreasing τ_depl with cosmic age. As discussed in Liu et al. (2019b) and Hodge & da Cunha (2020), this could be a signature of downsizing, meaning that more massive galaxies evolve at earlier times.

During our analysis, we find ${\tau }_{\mathrm{depl}}\sim {\left(1+z\right)}^{-1.07}$, which is more reflective of the scaling relations derived in Scoville et al. (2017) (${\tau }_{\mathrm{depl}}\sim {\left(1+z\right)}^{-1.04}$), rather than weaker dependencies (${\tau }_{\mathrm{depl}}\sim {\left(1+z\right)}^{-0.62}$ and ${\left(1+z\right)}^{-0.58}$) found by Tacconi et al. (2018) and Liu et al. (2019b), respectively. As expected, and in line with the literature results, we also find that galaxies above the MS (at a fixed M_* and z) form stars with a much higher efficiency (lower τ_depl) than their MS counterparts, with τ_depl ∼ ΔMS^−1.68. We would also like to caution the reader and highlight the fact that Tacconi et al. and Scoville et al. use functional forms that are different from ours, when fitting for evolution of τ_depl. For example, Tacconi et al. consider additional dependence on the effective radius R_e, which might inadvertently carry some redshift dependence. As such, the fitted exponents are not necessarily directly comparable. The differences between evolutionary trends could also be attributed to the different samples used (see, e.g., Hodge & da Cunha 2020).

Presumably, the existence of these outliers can be explained by an increased SFE, which results from major-merger events (see, e.g., Scoville et al. 2017; Cibinel et al. 2019). In fact, galaxies that lie above the MS are also found to have increased gas fractions (Dekel et al. 2009; Tacconi et al. 2020), which is attributed to a more efficient gas accretion from the cosmic web, but the enhanced gas reservoirs are still not large enough to explain significantly enhanced sSFR. The debates regarding the exact reason, which results in an onset of an SB-like mode of star formation, are still ongoing; however, it seems very likely that it is a combination of both increased gas fractions and enhanced SFE. We find that our sample supports this notion, with galaxies above the MS having both large gas reservoirs with median log(f_gas) = 0.15, meaning that gas mass reservoirs take up ∼59% of the total baryonic matter, and relatively short depletion times of ∼400 Myr. Our M_gas values were, however, derived with a general FMR, assuming solar-like metallicities. This, however, might not be applicable for SBs, which can display elevated metallicities owing to the increased sSFR. In fact, it has been shown (see, e.g., Silverman et al. 2015) that if SBs had supersolar metallicities, it would drive down δ_GDR, together with f_gas, and in turn result in increased SFE, thus implying that only the SFE is responsible for galaxies being elevated above the MS.

With the derived DMF in hand, we are also in position to bring together our findings with those presented in previous observationally driven studies and provide a direct comparison to the theoretical predictions as inferred by recent simulations. For our purposes, we focus on the 0.2 < z < 0.5 redshift interval that contains the majority of our objects and offers the most robust statistical analysis. These results are shown in Figure 15.

Zoom In Zoom Out Reset image size

We first compare our DMF to that presented in Pozzi et al. (2020), based on a PACS 160 μm selected sample of SFGs. In Figure 15, the two DMFs appear to be in tension in both the high- and low-mass ends, with our compilation overpredicting the number density of galaxies with high M_dust and underpredicting that of less dusty sources. The discrepancy between the two DMFs can be attributed to the choice of fitting methods/templates, the adopted κ_ν to infer M_dust, and selection effects. For example, the DL07 templates adopt a κ_{250 μm} of 0.51 m² kg⁻¹, while the analysis in Pozzi et al. uses κ_{250 μm} = 0.4 m² kg⁻¹, which would result in 0.1 dex smaller M_dust estimates. It is also important to point out that Pozzi et al. compared their modified blackbody (MBB) SED fitting to MAGPHYS, finding that their MBB method recovers systematically lower M_dust. Indeed, the choice of the fitting methodology can induce up to a factor of two difference in derived M_dust (see, e.g., Magdis et al. 2013 and an in-depth comparison in Berta et al. 2016). Moreover, for a flux-limited survey, the mere selection at λ_obs = 160 μm could introduce a bias toward warmer sources that for fixed L_IR have lower M_dust and that could explain the small number density of sources with log(M_dust/M_⊙) > 8. While it is not possible to correct for the effects of selection and broader SED fitting methodology, we have rescaled the Pozzi et al. DMF to have the same κ_{250 μm} as was adopted in our analysis.

We also compare our results to Dunne et al. (2011), who computed a DMF based on a sample of 250 μm selected galaxies. For our comparison, we have rescaled their DMF by −0.24 dex, to account for the difference in κ_ν . Contrary to Pozzi et al. (2020), we now find that Dunne et al. (2011) overpredict the number density of dusty galaxies at high dust masses. This again can be understood in terms of selection effects since the 250 μm selection could bias the sample toward cold sources and thus to higher M_dust values (again, for a flux-limited survey). While our criterion for at least one detection at λ_rest > 150 μm could be perceived as similar to a 250 μm selection at z ∼ 0.3, we note that the requirement for two extra detections at λ_rest < 150 μm and the super-deblended catalogs, which allow for the detection of fainter than the nominal confusion noise in the SPIRE bands, ease any bias toward either intrinsically cold or warm objects. This is further supported by the fact that our SMF-derived DMF, where the Eddington bias has already been corrected, falls directly between the calculations from Pozzi et al. and Dunne et al. (Figure 15). In conjunction with the derivation of a 〈U〉–z relation that is in excellent agreement with the stacking analysis of Béthermin et al. (2015), this suggests that the careful treatment of selection criteria and of the detection limits of our parent sample has allowed us to gain a unique and unbiased perspective on the evolution of dust properties of COSMOS galaxies.

Finally, we compliment our analysis by comparing our DMF to the theoretical predictions of the HAGN and IllustrisTNG simulations (Millard et al. 2020). In order to produce an HAGN DMF, we define a simulated sample following the procedure described above, bin the galaxies in 0.4 dex intervals of M_dust, and normalize by the volume of the simulation ($4\times {\left(142\right)}^{3}$ Mpc³). For the IllustrisTNG simulation, Millard et al. (2020) consider multiple TNG100 snapshots in a box size of 106 Mpc per side, comparable to the HAGN simulated subset presented earlier. The TNG-DMF is constructed through the M_dust values of the simulated galaxies, derived in post-processing through a fixed dust-to-metals ratio of 0.5.

Unlike real data, simulations do not suffer from observational bias and as such should be compared to DMF derived from the SMF, without adding the Eddington bias. As shown in Figure 15, both HAGN and IllustrisTNG are in excellent agreement with the high-mass end of our SMF-derived DMF. Notably, the HAGN DMF is also consistent with our results at the low-mass end down to M_dust ∼ 10⁷ M_⊙. We recall that, for simplicity, when converting M_gas to M_dust for the HAGN sample, we considered a universal δ_GDR = 100. However, for sources with lower M_* (< 10⁸ M_⊙) and thus with subsolar metallicities, a larger δ_DGR (∼150) is probably more applicable (Rémy-Ruyer et al. 2014). This would translate into a ×1.5 downward correction for the low-mass HAGN bins, bringing them into exact agreement with our DMF down to M_dust ≈ 10⁷ M_⊙. We note that this M_dust, assuming an average M_gas/M_dust ≈ 100, corresponds to the M_* completeness limit of the simulation (M_* ≈ 10⁹ M_⊙). Therefore, the observed decline of the number density of the HAGN galaxies at M_dust ≤ 10⁷ M_⊙ is fully consistent with the expectations.

In comparison to the TNG-DMF, though, we predict a factor of 2.5 fewer objects at the low-mass end. This tension could arise from the incompleteness of our sample at the low-mass end, which leaves the slope of the M_dust−M_* relation at M_dust < 5 ×10⁷ M_⊙ largely unconstrained. We are thus unable to ascertain whether this discrepancy is caused by the limitations of our sample, or whether the TNG simulations overpredict the number density of the galaxies in the low-mass end.

Put together, these comparisons indicate that, at least down to M_dust ≈ 5 ×10⁷ M_⊙, our M_dust−M_*−z relation and the resulting DMFs are robust and fully consistent with the theoretical expectations.

As briefly discussed in Section 7, during our analysis we identified some extreme outliers from the average f_dust and f_gas evolutionary trends (Figures 11 and 14), which typically have log(f_gas) > 0.5, i.e., their gas mass reservoir takes ∼75% of their baryonic matter. Since z_phot could be a major source of uncertainty in both M_* and M_gas, before looking further into this population of "gas giants," we first narrow down our sample to spectroscopically confirmed sources. We then examined the individual SEDs and the cutout images of the remaining sources in order to identify either poor coverage of the FIR peak or blending issues that could result in erroneously large M_gas estimates. With the above considerations, we are left with 41 objects whose extreme f_gas can only be explained by gigantic M_gas reservoirs. This population spans a wide range in redshift (0.21 < z < 4.05, 〈z〉 = 1.34), with $9.0\lt \mathrm{log}({M}_{* }/{M}_{\odot })\lt 11.3$, $\langle \mathrm{log}({M}_{* }/{M}_{\odot })\rangle =10.3$ and 0.11 < ΔMS < 14.2, 〈ΔMS〉 = 1.8. The best-fit SEDs of two such objects are presented in Figure 16. We also note that these two sources are otherwise unremarkable and have what can be considered "typical" values for the log(M_*) ∼ 10.7, and they also do not appear to be strong SBs (ΔMS = 3.8 and 2.2, respectively, for 10041706 and 10100707). Furthermore, the cutouts presented in Figure 16 indicate that these sources do not appear to be blended; therefore, the only unusual characteristic that they possess seems to be an elevated M_gas.

Zoom In Zoom Out Reset image size

A possible explanation for the very high M_dust and subsequently M_gas estimates for these galaxies, other than an extremely low δ_GDR, could be an optically thick FIR emission. In this scenario, the attenuation of the emission in the Wien part of the spectrum makes the galaxy appear cold, leading to an overestimate of its true M_dust (e.g., Jin et al. 2019; Cortzen et al. 2020). Since the DL07 models assume that the galaxy is optically thin at λ_rest > 1 μm, to test this scenario we employed MBB models of general opacity, leaving the effective wavelength (λ_eff) at which τ = 1 as a free parameter (see, e.g., Casey et al. 2012 for the functional form). We fixed the R-J slope to β = 1.8 and only fit the available photometry of each source at λ_rest > 40 μm. Due to the large number of free parameters in this model, we further limit our sample to sources with five or more IR detections. Out of the 41 "gas giants," we are thus able to constrain λ_eff for only 19. The distribution of the inferred λ_eff values is presented in Figure 17.

Zoom In Zoom Out Reset image size

We find that the vast majority of these objects have λ_eff > 100 μm, which implies that these galaxies could be optically thick in the FIR. The unusually high f_dust and f_gas can therefore be incorrect simply as a result of the optical depth effects. Indeed, a comparison between the inferred M_dust estimates for an optically thin and optically thick case yields an average ratio of ∼×1.8 for our sample. However, while this correction would reduce the average f_gas of the "gas giants" from $\langle \mathrm{log}({f}_{\mathrm{gas}})\rangle =0.72$ to $\langle \mathrm{log}({f}_{\mathrm{gas}})\rangle =0.47$, this is still substantially larger with respect to the average population of SFGs. Finally, the real f_gas could in fact be lower if the M_* of the sources is underestimated, a scenario that is indeed in line with a dusty, optically thick ISM.

To understand whether these objects do indeed host unusually high gas mass reservoirs and shed light on their nature, additional observations, with either ALMA or NOEMA, of M_gas tracers (e.g., CO and CI) are necessary.

In our fitting routine we utilize the dust models of Draine & Li (2007), with the updated opacity from Draine et al. (2014). These models aim for a robust and physically motivated approach to SED fitting in both MIR and FIR, as well as allow us to calculate the amount of luminous dust.

The description of the dust locked in the ISM is one of a mixture of carbonaceous and amorphous silicate grains, with their sizes and distributions following the extinction law in the Milky Way, the Large Magellanic Cloud (LMC), and the bar region of the Small Magellanic Cloud (SMC). The carbonaceous grains behave similarly to the PAH molecules, with their properties given by the PAH index q_PAH, which is defined as a fraction of dust mass locked into the PAH grains.

The models provide a bimodal description of the environments containing the interstellar dust: the diffuse ISM and the PDRs. The bulk of the dust mass is thought to be located in the cold and diffuse part of the ISM, which is being heated by a radiation field of a constant intensity ${U}_{\min }$. A smaller proportion of the mass budget described by the γ index is exposed to a gradient of radiation intensities ranging from ${U}_{\min }$ to ${U}_{\max }$ and is supposedly located in the warmer PDRs. Although these warm regions normally contain only a small fraction of the total dust mass, they can make a substantial contribution to the luminosity in the MIR SEDs. As described by DL07, the infinitesimal proportion of dM_dust exposed to radiation fields between U and U + dU can be modeled by a power-law distribution and, in the case of the diffuse ISM where ${U}_{\min }={U}_{\max }$, by a Kronecker δ-function. This leads to the following description:

$$\begin{eqnarray}&&\begin{array}{l}{{dM}}_{\mathrm{dust}}=(1-\gamma )\delta (U-{U}_{\min })\\ \,+\,\gamma \,{M}_{\mathrm{dust}}\displaystyle \frac{\alpha -1}{{U}_{\min }^{1-\alpha }-{U}_{\max }^{1-\alpha }},\end{array}\end{eqnarray} \tag{ A1 }$$

with ${U}_{\min }\leqslant {U}_{\max }$ and α ≠ 1. The parameter γ is the fraction of dust mass locked into the high starlight intensity regions described by the power law, α gives the distribution of radiation intensities in the PDRs, and M_dust is the total dust mass.

The methods described in DL07 allow us to compute a distribution of temperatures for all particles: ones that are small so their size makes them susceptible to the effects of quantized heating, and the larger ones where the steady-state temperatures dictated by the stellar radiation and radiative cooling equilibrium take hold. One can compute an averaged IR emission for a given grain type by first considering their temperature distribution and cross sections and then summing everything up to obtain the specific mass-weighted power that is being radiated by the dust exposed to starlight of intensity U. By integrating these numerical recipes from ${U}_{\min }$ to ${U}_{\max }$, one can obtain the power per unit frequency per unit mass ${p}_{\nu }({q}_{\mathrm{PAH}},{U}_{\min },{U}_{\max },\alpha )$.

In line with DL07, one can then model the galaxy SED as a linear combination of the diffuse ISM and the PDRs. This can be written as follows:

$$\begin{eqnarray}&&{j}_{\nu }=(1-\gamma )\,{j}_{\nu }[{U}_{\min },{U}_{\max }]+{j}_{\nu }[{U}_{\min },{U}_{\max },\alpha ],\end{eqnarray} \tag{ A2 }$$

where j_ν is the emissivity per hydrogen nucleon. If one now considers a galaxy at some distance D_L(z), the received flux density can be written as

$$\begin{eqnarray}&&{f}_{\nu }=\displaystyle \frac{{M}_{{\rm{H}}}}{{m}_{{\rm{H}}}}\displaystyle \frac{{j}_{\nu }}{{D}_{{\rm{L}}}^{2}(z)}.\end{eqnarray} \tag{ A3 }$$

Since j_ν is the quantity contained in the DL07 models, the normalization extracted from the fitting represents the total number of hydrogen nucleons and can be then converted to the luminous dust mass.

The total luminosity contained in both dust components can be written as

$$\begin{eqnarray}&&{L}_{\mathrm{dust}}=\langle U\rangle {P}_{0}{M}_{\mathrm{dust}},\end{eqnarray} \tag{ A4 }$$

with 〈U〉 representing a mean intensity of the radiation field, given by

$$\begin{eqnarray}&&\langle U\rangle =(1-\gamma ){U}_{\min }+\displaystyle \frac{\gamma \,\mathrm{ln}\,({U}_{\max }/{U}_{\min })}{{U}_{\min }^{-1}-{U}_{\max }^{-1}},\end{eqnarray} \tag{ A5 }$$

for α = 2 and P₀ denoting the power absorbed per unit dust mass in a radiation field of intensity U = 1.

In principle, one could think of these models as having six effective free parameters—q_PAH, ${U}_{\min }$, ${U}_{\max }$, α, γ, and M_dust—acting as the normalization as described above. It has been shown in Draine et al. (2007) that the parameter space is insensitive to the adopted dust model (MW, LMC, and SMC) and the values of α and ${U}_{\max }$. It is thus possible to recover a wide range of properties of various SEDs by fixing these to α = 2 and ${U}_{\max }={10}^{6}$. The values of ${U}_{\min }$ below 0.7 correspond to temperatures below 15 K, and while it is expected to find very few systems that exhibit this behavior, we have decided not to limit the range of ${U}_{\min }$ and allow it to vary between 0.1 and 50, to capture even the most extreme cases. It has been shown that, at least in the case of local galaxies in the Spitzer Nearby Galaxy Survey (SINGS), the ${U}_{\min }$ can be limited between $0.7\leqslant {U}_{\min }\leqslant 25$, and an MW-like dust model can be adopted to limit q_PAH between 0.004 and 0.046. Incorporating the use of the optimized set of parameters, which has been done in similar studies (e.g., Magdis et al. 2012a; Magnelli et al. 2012; Santini et al. 2014), might have a positive effect on computational speeds; however, we need to consider that this reduced parameter space has only been robustly verified for nearby solar-like metallicity populations of galaxies and might otherwise risk underestimating (overestimating) the dust masses for extremely cold (warm) systems.

We can thus extract the following physical parameters from the fit: γ, q_PAH, ${U}_{\min }$. The dust mass is simply computed from the normalization, while the L_IR is L_ν integrated over the 8–1000 μm range. As an additional parameter, we can obtain 〈U〉 by utilizing Equation (A5), or alternatively, as prescribed by Magdis et al. (2012a), we can use Equation (A4), where we set L_dust = L_IR and P₀ ≈ 125.

AGNs can have a significant impact on the ISM of galaxies that host them. They possess an ability to halt star formation by heating up the gas and dust or completely quenching the galaxy by stripping away its fuel. Under the common assumption (Antonucci 1993; Urry & Padovani 1995; Tristram et al. 2007), AGNs are surrounded by dusty tori, which, similarly to the ISM dust, can absorb the UV/optical light from the AGN and reradiate it at redder wavelengths, normally peaking in the MIR regime at 20–50 μm. It has also been shown that for select extreme cases (Mullaney et al. 2011) the AGN emission dominates the SED of a galaxy even at 60 μm, which presents a new challenge when calculating an infrared luminosity of a source. Infrared-derived SFR estimates rely on a robust understanding of the L_IR. Therefore, it is imperative to separate the energy contributions from hot dust in the ISM and a possible AGN.

In order to account for the effects of the IR contamination by AGNs when calculating L_IR, as well as to identify all the possible systems that might contain an active nucleus in our sample, we have decided to adopt a set of AGN templates from Mullaney et al. (2011, hereafter M11). These templates have been empirically derived by assuming an MBB function and fitting it to a set of Swift-BAT AGNs and IRAS spectra. The obtained models describe intrinsic AGN emissions in the range spanning from MIR to FIR (6–100 μm). In this case, a typical AGN SED could be thought of as a broken power law at ≤40 μm, which rapidly vanishes when moving above 40 μm. The average intrinsic AGN emission can be described as follows:

$$\begin{eqnarray}{F}_{\nu }=\left\{\begin{array}{ll}{\lambda }^{1.8} & \mathrm{at}\quad 6\,\mu {\rm{m}}\lt \lambda \lt 19\,\mu {\rm{m}}\\ {\lambda }^{\alpha } & \mathrm{at}\quad 19\,\mu {\rm{m}}\lt \lambda \lt 40\,\mu {\rm{m}}\\ {\nu }^{1.5}{F}_{\nu }^{\mathrm{BB}} & \mathrm{at}\quad \lambda \gt 40\,\mu {\rm{m}}\end{array}\right.,\end{eqnarray} \tag{ A6 }$$

where $F{{}_{\nu }}^{\mathrm{BB}}$ is the MBB function and α is the spectral index. In our procedure we utilize the high- and low-luminosity templates, with α = 0.0 and 0.4, respectively. In addition to that, we allow for a linear combination between the two, with varying coefficients, thus expanding the existing template space.

The contribution of stellar emission in the observed NIR bands, such as Spitzer IRAC 1-4, can be quite significant, especially when we move to higher z. Therefore, if no libraries representing the luminosity from dust-attenuated stellar light are available, one might either underestimate the slope of the AGN or overestimate its normalization. This can lead to erroneously assigning more luminosity to the AGN component in the MIR and therefore underestimating the L_IR.

To avoid this, we additionally incorporate a library of stellar emission models, which are an updated version of the templates described by Brammer et al. (2008, hereafter GB08). These are based on the SPS models (Conroy et al. 2009), which were optimized for deep optical−NIR broadband surveys. The models form a basis set of a larger library and were derived by using the "nonnegative matrix factorization" algorithm that was described in Blanton & Roweis (2007). The method attempts to reproduce the full library of templates, by finding a nonnegative linear combination of a smaller number of models. These can be considered as the "principal component" blueprint of the larger catalog. In total, there are 12 optical SEDs, which include both dust-attenuated and nonattenuated starlight. We have incorporated these models into our fitting routine, and if the UV−optical photometry is available, this allows us to extract properties such as M_*, SFR, E(B−V), and SFH. If no UV−optical data are available, the addition of these templates is still useful, as they can account for the excess flux in the rest-frame NIR bands, in conjunction with AGN and IR templates.

To test how our fits behave without the stellar component, we have isolated ∼100 objects with z > 2, so that our bluest available band traces λ_rest ∼ 1 μm, where the contribution from stellar emission becomes nonnegligible. We then exclude GB08 templates and refit our objects. This results in two outcomes that can be seen in Figure 18. We find that, by removing the additional component, we tend to either overestimate the AGN contribution (blue squares), with the median being ∼1.5, or alternatively erroneously assign a galaxy to contain an AGN (red points).

Zoom In Zoom Out Reset image size

In order to model an SED of a galaxy and extract the physical parameters, we first transport all three components—the stellar, the dusty torus AGN, and the infrared dust emission—to a common wavelength grid spanning the range from 10⁻⁴ to 10⁵ μm. Our method relies on a linear combination of these models, thus making it imperative for them to share a common range, so that co-adding them is made possible. In certain cases where this grid falls outside the original range of the template, we extrapolate blueward and redward by using a steeply declining power law. This was done to ensure that the resultant galaxy emission is continuous without any sudden breaks, which could interfere with the fitting, where one of the templates has ran out of range. These added power laws do not introduce any additional emission, as the flux density contribution from them is orders of magnitude lower than that of the original template. We then redshift the wavelength grid on a per-galaxy basis and normalize the templates to ensure that they are not separated by tens of orders of magnitude. In Figure 19 we show all the templates used in Stardust, normalized to the K band, and in Table 2 we list all the relevant parameters of the models.

Zoom In Zoom Out Reset image size

Table 2. Template Parameters Used for Stardust Fit

Parameter	Value
Dust Emission: DL07 (Updated in Draine et al. (2014))
	$\left[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,\right.$
${U}_{\min }$	1.0,1.7,2.0,3.0,4.0,5.0,7.0,8.0,
	$\left.10.0,12.0,15.0,20.0,25.0\right]$, 26 values
qPAH	[0; 0.1], 11 values in steps of 0.1
	[0,0.001, 0.0025, 0.005, 0.0075]
γ	+[0.01; 0.1], 9 values in steps of 0.01
	+[0.2, 0.35, 0.5], 17 values in total
Optical Emission: Brammer et al. (2008) a
AV (Calzetti et al. 2000)	[0.6, 0.12, 0.19, 0.29, 1.05, 2.68,
	0.11, 0.36, 0.98, 1.54, 1.97, 2.96]
M/LV	[0.38, 0.76, 1.68, 4.01, 6.45, 44.48,
	0.12, 0.21, 0.33, 0.64, 1.57, 4.00]
log10(sSFR)	[−10.75, −11.37, −11.90, −12.53, −12.05,
	−12.47, −8.37, −8.60, −8.50, −8.57, −8.93, −8.90]

Note.

^a Please refer to Brammer et al. (2008) for a more detailed description of the creation and selection of these basis-set templates. See Blanton & Roweis (2007) for a methodology regarding the SFH.

Download table as:  ASCIITypeset image

Subsequently, we perform a synthetic photometry on all three components separately in all observed bands where data are available. This is done by convolution of the filter transmission curves with the model SEDs. The resultant synthetic fluxes for each template and available observed band are then all combined into a two-dimensional matrix and passed onto the nonnegative least-squares (nnls) algorithm in the Python scipy.optimize.nnls package, which finds the best solution vector for that object. The nnls is a simple minimization algorithm that in our case takes the form of

$$\begin{eqnarray}&&\min | | {A}_{b,t}{x}^{t}-{y}^{b}| {| }_{2}\quad \mathrm{where}\quad x\geqslant 0,\end{eqnarray} \tag{ A7 }$$

where A is the uncertainty-weighted template matrix, y is the S/N for each band, x is the solution vector, and ∣∣. ∣∣₂ signifies the Eucledian norm. Instead of finding the best-fitting template, the algorithm computes the best-fit coefficients x, with x ≥ 0. The number of simultaneously fit templates is user defined. For our purposes we chose to fit all 12 GB08 and 2 M11 templates at the same time, each time combining them with a single DL07 template, iterating through all the possible combinations, and then finding the best fit. Following Equation (A7), this is done by building the matrix A, for each band b and each template t, where the first 14 templates in the matrix remain fixed, while the last element is being continuously replaced with a new DL07 model and looped over. In the end we obtain 4862 possible best-fit vectors—x^t , one per each DL07 model. The final result is then extracted from the χ² distribution of all best-fit solutions. The advantage of this method is in avoiding progressively looping over all possible template combinations to find the best solution, and instead only choosing to loop through DL07 templates. This approach significantly reduces the amount of required computational resources. The resultant solution vector encodes the individual contributions from each template to the total emitted flux in each band. These are then added together to return the best-fit solution. This three-component split allows us to predict exactly how much each component contributes to the source's L_IR, thus allowing us to differentiate between the AGN and the warm ISM dust emissions. In addition to that, the normalization of DL07 templates leads us directly to the number of hydrogen nucleons, linking it to dust mass via Equation (A4). The radio data points are not being considered by our fitting routine; however, we add a power-law radio slope with a spectral index of −0.75, as described by the FIR−radio correlation in Delvecchio et al. (2021), mainly for visualization purposes, but also to detect the existence of the AGN radio excess.

Table 3. Structure of the Best-fit Catalog ^a

Column Name	Units	Description
Catalog	⋯	Catalog identifier code
ID	⋯	ID of the object from the parent catalog
R.A.	deg	R.A. coordinate, as given in the parent catalog
Decl.	deg	Decl. coordinate, as given in the parent catalog
Area	⋯	Same as goodArea flag in the parent catalogs
z	⋯	Redshift used for fitting
ztype	⋯	Redshift type, 1 for spectroscopic, 0 for photometric, and 2 for zphot from EAZY
LIR_total	L⊙	Total FIR luminosity, obtained as the sum between AGN and DL07 components
eLIR_total	L⊙	Uncertainty on the total FIR luminosity
Lagn_total	L⊙	AGN luminosity
eLagn_total	L⊙	Uncertainty on the AGN luminosity
Lir_draine	L⊙	Luminosity given by the DL07 template
eLir_draine	L⊙	Uncertainty on the luminosity given by the DL07 template
MD	M⊙	Dust mass as predicted by the best-fit DL07 template
eMD	M⊙	Uncertainty on the dust mass
deltaGDR	⋯	Gas-to-dust ratio
MG	M⊙	Gas mass computed from δGDR and Mdust
eMG	M⊙	Uncertainty on gas mass
Mstar	M⊙	Stellar mass, equal to the one in the original catalog
lastdet	μm	Last band that has an S/N > 3 detection (λlast); given in rest frame
chi2	⋯	The χ2 of the best-fit coefficients
f_agn	⋯	AGN fraction, given as LAGN/LIR,total
efagn	⋯	Uncertainty on the AGN fraction
fgas	⋯	Gas fraction computed as Mgas/Mdust
fgas_FMR	⋯	Gas fraction computed assuming δGDR = 100
Umin	⋯	Best-fit ${U}_{\min }$
gamma	⋯	Best-fit γ
qpah	⋯	Index of the best-fit qPAH value
U	⋯	Average radiation field intensity 〈U〉
sU	⋯	Uncertainty on 〈U〉
deltaMS	⋯	SFR/SFRMS, where SFR = LIR,draine ×10−10 and SFRMS is given by Schreiber et al. (2015)
e_deltaMS	⋯	Uncertainty on ΔMS

Notes.

^a Stardust also returns M_*, A_V, and UV−optical SFR; these are not included in the release version of the catalog but are available upon request. ^b We have used the 7970d55 (28/06/21) version of Stardust to produce these catalogs.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

To derive luminosity estimates, we integrate the three-component summed SED in the 8–1000 μm. This gives us the L_IR,tot that contains within itself the energy emitted by the ISM dust and AGN torus, if present. The contribution of stellar emission at λ_rest > 8 μm is negligible; therefore, we do not go through an additional step of subtracting those models. We then integrate the best-fit template with just the AGN contribution in it, to obtain the L_AGN and L_IR = L_IR,tot−L_AGN. We also compute the f_AGN = L_AGN/L_IR,tot to estimate how strongly the infrared SED of a galaxy is contaminated by AGN activity. In addition to that, it allows us to separate our sample into objects that have an active nucleus and those that do not. The conditions of the ISM in these different environments may vary quite significantly and would affect the extracted scaling relations if not treated correctly.

The normalization of the DL07 models returns the number of hydrogen nucleons as one of its free parameters, from which we can obtain the fiducial M_H. We then compute M_dust by converting this quantity assuming a fixed gas-to-dust ratio of δ_GDR = 100, as prescribed in DL07. It is important to note that this ratio is encoded into the models and does not represent an actual physically meaningful conversion factor. We also compute a T_dust proxy in the form of the average radiation field intensity 〈U〉, from Equation (A4), by assuming a P₀ = 125. We note, however, that this quantity represents the luminosity-weighted dust temperature and has little to no bearing on the temperature of the cold dust or gas (see, e.g., Scoville et al. 2016, 2017). The structure of the final output catalog is presented in Table 3.