© ESO 2021

1. Introduction

The rest-frame ultraviolet (UV) emission of galaxies is one of the more direct probes of star formation in galaxies. However, its interpretation as a star formation rate (SFR) depends on a number of factors, some of which may not be well constrained, including the following: star formation history (SFH), stellar metallicity, stellar atmospheres, stellar tracks, and perhaps most importantly, dust attenuation.

Even though it represents typically less than 1% of the mass of the interstellar medium (ISM; e.g., Rémy-Ruyer et al. 2014), dust can have a considerable effect on the emerging UV radiation of galaxies, both dimming and reddening its emission. The SFR derived from the observed UV luminosity, without any dust corrections, is often referred to as unobscured star formation.

The reliability of the UV emission as a SFR estimator therefore critically depends on the ability to accurately correct for dust attenuation. One approach for achieving this goal is to relate the observed UV color (β) to the infrared (IR) excess (IRX), where the latter is closely related to the intrinsic UV attenuation (e.g., Meurer et al. 1999). Larger UV attenuation translates into a redder UV color so that the difference between the intrinsic and observed UV colors should ideally trace the attenuation. However, a key issue with this approach is that the relation between the observed UV color and the attenuation is strongly dependent on the attenuation curve (Boquien et al. 2009; Salim & Boquien 2019), which is often poorly known and yet spans a wide range of shapes and varies significantly from galaxy to galaxy (Salim & Narayanan 2020). While, the determination of the total SFR using UV and optical spectral energy distribution (SED) modeling (e.g., Conroy 2013) is a more sophisticated approach than the IRX–β technique, it is based on the same fundamental principle (comparison of observed and intrinsic colors), and it is therefore subject to similar systematic and random errors arising from poorly constrained or inappropriately assumed attenuation curves.

Stellar luminosity absorbed by dust grains is reprocessed and re-emitted at longer wavelengths, mainly in the mid-infrared (MIR, λ ≲ 40 μm) and far-infrared (FIR, λ ≳ 40 μm) domains. Dust emission is thus complementary to the UV emission emerging from galaxies as a tracer of star formation (e.g., Daddi et al. 2007; Elbaz et al. 2007). In heavily obscured galaxies, in which only a small fraction of UV photons can escape and therefore a reliable correction of the UV or UV- and optical-based SFR is particularly challenging, the dust emission traces nearly the entirety of star formation (Kennicutt 1998). Overall, the IR emission provides the most reliable “correction” for the unobscured SFR determined from the UV and optical emission. Furthermore, as recognized previously (e.g., Inoue 2002) and further demonstrated in this paper, the IR emission can be used as a reasonably good indicator of the total SFR of normal star-forming galaxies in its own right, even without explicitly accounting for the unobscured star formation.

Both the monochromatic (λ × L_λ(λ)) and the total (bolometric) infrared luminosity (L_TIR) have been employed as SFR estimators (e.g., Calzetti et al. 2007). In this paper, we, among other objectives, tackle the question of which wavelength range is the best SFR indicator. Regardless of the answer to that question, the bolometric dust emission (L_TIR) is a fundamental galaxy property that we wish to be able to determine as reliably as possible, and its measurement represents the principal objective of this work.

Measuring the bolometric dust emission in a galaxy is a difficult and complex affair. The underlying reason for this difficulty is that finely sampled SED (or spectra) that would cover the entire IR range from a few microns to ∼1 mm simply do not exist. Even the relatively well sampled SED, which are by no means common, still require interpolations and extrapolations in order to produce continuous IR spectra. These interpolations and extrapolations require theoretical modeling of dust emission, which is in itself a rather complex task.

The complexity of dust modeling arises in no small part from the intricate emission spectrum of dust, which depends on the nature of the grains (composition and size distribution), as well as the local physical conditions, in particular the shape and intensity of the incident radiation field. The most prominent feature of the dust emission spectrum consists in a gray body component peaking at around 100 μm, which modern models describe as a combination of the emission of large dust grains over a broad range of temperatures, reflecting the temperature distribution throughout a galaxy. At shorter wavelengths, the emission is dominated by Polycyclic Aromatic Hydrocarbons (PAH) that exhibit prominent bands that appear to be strongly dependent on metallicity. PAH emission, as well as the emission from very small grains that dominates at wavelengths between PAH and the peak, comes from a stochastic heating process.

Over the past decades, dust emission models have increased in sophistication and complexity, improving our ability to reproduce and interpret the observations (e.g., Desert et al. 1990; Dale et al. 2001; Draine & Li 2007; da Cunha et al. 2008; Compiègne et al. 2011; Draine et al. 2014; Jones et al. 2017). However, because of their relative complexity and ensuing flexibility, such models are of little help for deriving the total IR luminosity in cases when the observations are confined to only a small wavelength range, or even a single band, that is, when huge extrapolations are required. Such is the situation for hundreds of thousands of galaxies observed by the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010), as will also be the case of the galaxies that will be studied through the eye of the James Webb Space Telescope (JWST, Gardner et al. 2006). In both cases the observations only extend up to ∼25 μm. Fortunately, making inferences about the bolometric emission based on a very limited wavelength coverage is nevertheless possible, but one has to resort to using the templates.

Templates are important because they narrow down the vast multidimensional parameter space of theoretical spectra into a family that is parametrized on a single parameter, which is calibrated based on a sample of galaxies with good SED coverage. Parameterization should ideally be on a parameter that best correlates with the changes in the shape of the SED. The dust emission templates approach was pioneered by Chary & Elbaz (2001), who produced a set of templates dependent on the total infrared luminosity, which in effect provides constraints on the shape of the dust emission spectrum. If the galaxies from which the templates were constructed are representative of the galaxies to which they are applied, unbiased estimates of L_TIR will be possible from limited data. Because in this case the parameterization is on the extensive quantity (i.e., it directly scales with the “extent” of the galaxy), even a single flux point (band) is sufficient to estimate the total dust luminosity.

However, these and other widely used templates (e.g., Dale & Helou 2002; Rieke et al. 2009) were based on galaxies selected from relatively shallow surveys carried out by the Infrared Astronomical Satellite mission (IRAS, Neugebauer et al. 1984) and the Infrared Space Observatory (ISO, Kessler et al. 1996). Such a selection will inevitably favor high-luminosity galaxies that are atypically luminous for their stellar or dust mass, possibly because of the more efficient star formation resulting from mergers that are common among these galaxies. As a result, they will have a warmer SED than normal star-forming galaxies of similar L_TIR, which could lead to offsets in the estimation of dust luminosities of normal galaxies, especially based on observations beyond the peak, even at low redshift (Lin et al. 2016). Furthermore, the relatively small sample sizes (∼100 galaxies) that underlied these efforts precluded the investigation of more complex parameterizations. Such parameterizations may extend the validity of templates over a broader range of galaxies and redshifts.

The upcoming launch of JWST highlights the need for new templates that provide an un-biased view of galaxies at higher redshifts. In the absence of adequate detailed templates, the limited spectral coverage of JWST may negatively affect our ability to estimate the bolometric dust emission to trace star formation. Existing templates from local luminous samples are not necessarily representative of similarly luminous galaxies at higher redshifts (Safarzadeh et al. 2016). Indeed, when applied to higher redshifts, some of the existing templates are suspected to produce systematic offsets in L_TIR derived from MIR monochromatic luminosities (Lin et al. 2016). Efforts to produce IR templates using the actual high-redshift galaxies have made significant progress in recent years, primarily as the result of the stacking of Herschel observations. For now, the focus has been mostly on producing average templates in various redshift bins (Magdis et al. 2012), which cannot capture the diversity of dust emission properties present at a given redshift. The stacking approach was refined in Schreiber et al. (2018), where the templates depend on the redshift but also on the SFR relative to the main sequence. However, lacking a parameterization on an extensive quantity, such templates cannot be utilized when observations exist in only a single band.

To address these concerns with the existing templates based on local galaxies and to provide a more comprehensive physically motivated basis for interpreting IR observations in general, in this work we exploit a significantly larger sample of local galaxies from MIR (WISE) and FIR (Herschel-ATLAS) surveys, and with photometry from the far UV to the submillimetric regime, to gain insight into the drivers of IR SED shape and the ability of different monochromatic luminosities to constrain dust parameters. Based on that analysis, we construct new single and two-parameter templates, where the role of the second parameter is to account for the range of normal star forming galaxies at different redshifts, thus making the new templates specifically well suited for JWST. Furthermore, in this work we highlight the use of explicit relations as a more practical alternative to discrete templates. Our templates and relations are constructed using some new approaches that are physically motivated and homogeneously applied to all the data.

We introduce the multiwavelength data set we build on in Sect. 2. In Sect. 3 we describe the method we have developed for constructing and parameterizing the new templates. We present our results in Sects. 4 and 5 and discuss them in Sect. 6 before concluding in Sect. 7. Throughout this paper we assume a WMAP7 cosmology (H₀ = 70.4 km s⁻¹ Mpc⁻¹, Ω_c = 0.226, Komatsu et al. 2011).

2. Data and sample

2.1. Herschel-ATLAS

The sample for the construction of dust emission templates, or, more broadly, the relations that allow us to estimate the TIR luminosity and other physical properties from IR observations, should have the following characteristics: (1) good wavelength sampling across the full IR range in order to probe the main dust emission features, (2) deep IR observations in order to provide reliable and highly complete photometric measurements, (3) representative sampling of star-forming galaxies, including high-mass galaxies, to take into account variations in dust emission across galaxy types, (4) extensive UV and optical data to constrain the stellar masses and dust-corrected SFR, (5) ability to distinguish and remove AGN, which can contaminate the IR emission, and, finally (6) relatively large size, in order to adequately sample the parameter space and reduce statistical uncertainties. These conditions are fulfilled for a sample selected from the Herschel-ATLAS survey (H-ATLAS, Valiante et al. 2016; Maddox et al. 2018), combined with multiwavelengths data from SDSS (Sloan Digital Sky Survey, York et al. 2000), GALEX (Galaxy Evolution Explorer, Martin et al. 2005), and WISE (Wide-field Infrared Explorer, Wright et al. 2010), which form the basis for the SDSS-GALEX-WISE Legacy Catalog (GSWLC, Salim et al. 2016, 2018).

2.2. Data

H-ATLAS is the largest uniform imaging survey carried out by the Herschel Space Observatory, covering 660 sq. deg. It observed the sky in five FIR and submillimeter bands, using the PACS (100 μm and 160 μm, Poglitsch et al. 2010) and SPIRE (250 μm, 350 μm, and 500 μm, Griffin et al. 2010) cameras. H-ATLAS data release 1 (DR1, Valiante et al. 2016) consists of three equatorial fields with a total area of 160 sq. deg covering the Galaxy and Mass Assembly (GAMA, Driver et al. 2011) survey, whereas two fields covering a total of 500 sq. deg around the north and the south Galactic poles were released as DR2 (Maddox et al. 2018).

The SPIRE photometry of the H-ATLAS catalog is based on detections in the 250 μm image, which has the highest density of sources in any Herschel band (Maddox et al. 2018). This photometry is carried out by successively measuring and subtracting the sources in the order of their brightness, separately in each filter. Most SDSS galaxies are unresolved in SPIRE images, the beam size ranges from 18″ at 250 μm to 35″ at 500 μm, so the default Point Spread Function (PSF) photometry is adequate for them. For extended sources (based on their optical size), which are found among z < 0.1 galaxies, aperture photometry is also performed. The flux we adopt (reported in the catalog as “best”) is whichever of the aperture or PSF flux is the brightest. Because the PACS beam is smaller (7″ at 100 μm and 11″ at 160 μm), only aperture photometry is performed in these bands, centered on the position of the optical counterpart. The H-ATLAS catalogs provide matched counterparts in SDSS, which were identified following the methodology of Bourne et al. (2016). We add in quadrature a flux calibration uncertainty of 7% for PACS and 4% for SPIRE to the flux uncertainties reported in the H-ATLAS catalogs.

For the MIR photometry (12 μm and 22 μm), we use the all-sky survey from WISE, processed by the unWISE project (Lang et al. 2016). They performed forced photometry using SDSS positions and galaxy profiles, thus obtaining the fluxes that are more accurate compared to the standard pipeline photometry, which treats all sources as unresolved. Forced photometry is important in particular for sources that are resolved and/or blended in WISE images. We add in quadrature a 2% flux calibration uncertainty to the formal uncertainties reported in the unWISE catalog.

For deriving the stellar masses and SFR we use near-IR (NIR), optical, and UV photometry. NIR photometry in three bands (J, H, and K_s) comes from the 2MASS (Two Micron All Sky Survey, Skrutskie et al. 2006) Extended Source Catalog, as described in Salim et al. (2016). Optical photometry in five bands (u, g, r, i, and z) comes from SDSS DR10. Far- (FUV) and near-UV (NUV) photometry is obtained from the GALEX data release GR6/7, with the various flux corrections that are discussed in Salim et al. (2016), which also described the matching algorithm to SDSS. We use GALEX data regardless of the depth of UV observations (GSWLC-X2 dataset), covering 90% of SDSS spectroscopic targets.

To be included in GSWLC, and therefore in our current sample, in addition to being covered by GALEX, SDSS spectroscopic targets need to have an SDSS Petrosian magnitude in the r band brighter than 18 and a redshift within 0.01 < z < 0.3. In other words, they need to be similar to galaxies selected for a statistically complete SDSS Main Galaxy Survey (Strauss et al. 2002).

2.3. Sample

The H-ATLAS catalogs¹ (DR1 v1.2 and DR2 v1.4) contain a total of 10676 GSWLC-X2 galaxies over ∼400 sq. deg. From this initial sample we remove 127 AGN with broad lines (based on their SDSS spectroscopic classification), 3759 galaxies classified as AGN based on the BPT emission line diagram² and a further 29 galaxies not contained in the MPA/JHU catalog. The removal of AGN is essential, as they can significantly contribute to the IR spectrum, in particular in the WISE bands. From the remaining 6761 galaxies we remove 1786 galaxies with negative fluxes in the W3 and/or W4 bands. Furthermore, of the remaining 4975 galaxies we remove 2386 galaxies with a signal-to-noise ratio lower than two at 100 μm and/or 160 μm, and, finally, another five galaxies with sSFR < 10⁻¹² yr⁻¹. Our final sample consists of 2584 galaxies. Its median redshift is 0.08 with a 90 percentile redshift range between 0.02 and 0.17. Our sample is an order of magnitude larger than the samples used in most previous efforts to build dust emission templates³.

2.4. Characteristics

To estimate the physical properties of our sample, we have carried out SED modeling as described in detail in Sect. 3.1. In Fig. 1 we present a comparison of our sample with respect to all GSWLC-X2 galaxies in the sSFR–M_star plane. We see that our sample follows well the main sequence of star-forming galaxies.

Fig. 1. Distribution of the selected sample in the sSFR–Mstar plane (red circles) compared to the full GSWLC-X2 catalog (hexagons, with the density indicated by the color bar).

In Fig. 2 we present the distributions of the physical properties of the sample measured: SFR averaged over the last 100 Myr, stellar mass (M_star), total infrared luminosity (L_TIR, corresponding to the integral of the dust emission between 8 μm and 1 mm), dust mass (M_dust), specific SFR (sSFR), and oxygen abundance. Overall, the sample covers a wide range of non-dwarf, normal star-forming galaxies, including luminous IR galaxies (LIRGs). Furthermore, it includes galaxies with a large range of sSFR at a given mass, up to log sSFR ∼ −9. It does not include more extreme, but relatively rare galaxies such as ultraluminous IR galaxies (ULIRGs, log L_TIR > 12), which can have a markedly different dust emission spectrum compared to more typical star-forming galaxies (e.g., Rieke et al. 2009). In Sect. 6.3 we discuss the cases in which the relations that we derive in Sect. 5 can be used for galaxies that fall outside of the range of properties of our sample.

Fig. 2. From the top-left to the bottom-right, distribution of LTIR, Mdust, SFR, Mstar, sSFR, and the oxygen abundance. For the latter the sample is restricted to 2391 objects rather than 2584 for the other physical properties as the oxygen abundance could not be computed for all the objects.

3. Method

The central aspect of this work is the construction of relations between one or more IR fluxes and general galaxy properties. These relations will serve as estimators of the bolometric dust luminosity and will also allow for the generation of dust emission templates parameterized on the TIR luminosity⁴ and/or other physical properties. To carry out these tasks it is first necessary to estimate both the physical properties (stellar mass, dust mass, SFR, sSFR, etc.) of each galaxy in our sample and their associated dust emission spectrum.

3.1. Estimating the physical properties and spectra of galaxies

One of the most physically motivated techniques to estimate the physical properties of galaxies is via the modeling of their electromagnetic emission using synthesis population models. Since the first models of Tinsley (1972), numerous codes have been developed over the years with increasing sophistication and power. The latest generation of codes can model galaxies from the FUV to the FIR, including dust in absorption and in emission in a consistent way, a key ingredient to break degeneracies and estimate the physical properties more reliably (Burgarella et al. 2005; da Cunha et al. 2008; Noll et al. 2009; Leja et al. 2017; Carnall et al. 2018; Boquien et al. 2019; Bowman et al. 2020; Robotham et al. 2020).

For this work we adopt CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019)⁵ version 2020.0, a highly versatile code that has been successfully used for modeling galaxies and study a broad range of questions at different redshifts (e.g., Buat et al. 2011, 2012, 2018, 2019; Burgarella et al. 2011, 2020; Boquien et al. 2012, 2014, 2016; Ciesla et al. 2015, 2016, 2017, 2018, 2020; Salim et al. 2016, 2018; Salim & Boquien 2019; Hunt et al. 2019; Franco et al. 2020; Dale et al. 2020; Mountrichas et al. 2021). Under its standard operating mode, CIGALE generates and fits a large number of multiwavelength spectro-photometric models to observations, and estimates the physical properties and the related uncertainties from their probability distribution function. Here, however, we aim not just to measure the general physical properties of galaxies, but also to derive the best-fitting dust emission spectrum for each galaxy to be used as a “ground truth” for constructing the aforementioned relations or templates. Both of these aims could in principle be achieved with CIGALE by simultaneously modeling and fitting the UV, optical, and NIR photometry (i.e., the stellar populations) together with the MIR and FIR photometry (dust emission). However, if we wish to use the dust emission models required to fit the IR emission with sufficiently high resolution in model parameters (the model grid), we would end up with an excessively large number of combined (stellar populations and dust) models (in our case one trillion). We have thus decided to adopt and adapt (Sect. 3.1.3) the two-step approach developed in Salim et al. (2018). In a nutshell, we first model only the emission of the dust to estimate the physical properties that are largely independent from the stellar population modeling (e.g., dust mass and luminosity), as well as to obtain the dust emission spectrum. We require a detailed dust emission spectrum as a basis for the construction of relations and templates that are continuous in wavelength, and therefore can be used with any instrument and at any redshift. In the second step, we fit models of the stellar populations to UV, optical, and NIR photometry and estimate different physical properties (SFR, stellar mass, etc.). Importantly, we do so while also including in the fitting the dust luminosity (L_dust) determined in the first step as a critical information to break the age-attenuation degeneracy. Finally, with the knowledge of the stellar spectra, we correct for the stellar contamination in the MIR bands, and repeat these two steps.

3.1.1. Dust modeling

We model the dust emission using the Draine et al. (2014) update of the Draine & Li (2007) models. We base our choice on the fact that these models are physically motivated while providing enough flexibility to accurately fit our 7-band IR photometry. The Draine & Li (2007) models are known to adequately reproduce well-sampled IR SED (e.g., Ciesla et al. 2014; Hunt et al. 2019) and have been successful at reproducing the emission of galaxies (e.g., Aniano et al. 2012, 2020, for the KINGFISH sample that contains diverse galaxies in the nearby universe), while showing great flexibility to simultaneously reproduce the emission of the PAH, and the warm and cold dust components. It is thus unlikely that there is any major issue with this underlying model. The main caveat may be in the 25 μm to 60 μm range as we see in Sect. 6.1, with little data being available in this range to test and constrain the models in the first place. Each model contains four free parameters and is based on a combination of two components. One, illuminated by an incident radiation field of intensity U_min, represents the diffuse dust across the galaxy. The other component models the dust in star-forming regions, with an illumination intensity following a power law dM_dust/dU ∝ U^−α, where M_dust is the dust mass, U the incident radiation field intensity, and α is the adjustable power law index. A coefficient γ sets the mass fraction of the component associated with star-forming regions. We use the entire parameter space provided by the model grid: U_min ranges from 0.1 to 50, α ranges from one to three, and q_PAH, the mass fraction of PAH, ranges from 0.47% to 7.32%. We sample γ from 10^−2.5 to 10^−0.3 in 15 logarithmically-spaced steps, a range over which it has the most pronounced effect on the shape of the spectra. The final grid consists of 124 740 models at each redshift. This modeling is performed with the DL2014 module in CIGALE.

We model the dust emission based on WISE and Herschel data. In the Herschel SPIRE bands, we sometimes only have upper limits. In order to exploit this information we also include these upper limits in the modeling, with the computation of the goodness-of-fit described in detail in Boquien et al. (2019). Since the modeling is performed at the observed redshifts known from SDSS spectroscopy, no K-correction is required. The output of dust modeling includes estimates of L_dust and M_dust, their uncertainties, as well as the best-fitting spectra made of 1001 flux densities at a constant separation in log space of 0.004 dex from 1 μm to 10 mm.

An inspection of the fits revealed that there was a relatively small, but statistically significant systematic offset between the observations and the best-fit model at 160 μm, with the former being 7.3% higher than the latter on average. The cause of this offset is unclear. Fitting the observations using the THEMIS models (Jones et al. 2017), we observe a similar offset, suggesting that it may not have a physical origin. Private communication with the members of the Herschel-ATLAS team did not allow us to eliminate the possibility of an instrumental origin. In any case, deriving the suite of templates and estimators by adjusting the fluxes 7.3% downward, we find that the differences are very small, typically within 0.01 dex for monochromatic estimators. Given the minute differences, we adopt the Herschel-ATLAS fluxes as-is.

3.1.2. Stellar population modeling

With the dust luminosity in hand, we model the stellar populations with CIGALE to estimate other physical properties. The stellar emission is computed adopting the Bruzual & Charlot (2003) single stellar populations (CIGALE module bc03) following a Chabrier (2003) initial mass function and a metallicity ranging from subsolar (Z = 0.004) to supersolar (Z = 0.05). The SFH is described by two decaying exponentials (CIGALE module sfh2exp), one modeling the general stellar populations and the other one the latest episode of star formation. We assume the galaxies to be 10 Gyr old. The exponential describing the general stellar populations has an e-folding time ranging from 850 Myr to 20 Gyr. The onset of the latest episode of star formation, which has a large timescale, occurred between 100 Myr and 5 Gyr before the time of observation and this episode represents between 0% and 50% of the total mass of star formed over the lifetime of the galaxy. Our models also include the nebular emission (CIGALE module nebular) to account for the contamination of broadband fluxes by gas ionized by massive stars. This comprises both emission lines and continuum emission (free-free, free-bound, and 2-photon processes). Finally, we attenuate the stellar and nebular emission using the Noll et al. (2009) generalization of the Calzetti et al. (2000) attenuation law (CIGALE module dustatt_calzleit). In short, the diversity of attenuation curves (Salim & Narayanan 2020) is obtained by multiplying the Calzetti et al. (2000) curve by a power law of index δ ranging from −1.4 (steeper) to 0.4 (shallower) and by adding a bump at 217.5 nm in strengths ranging from zero to twice that of the Milky Way. We allow for a differential reddening of a factor 0.44 between stars younger and older than 10 Myr. The absolute reddening of the younger population goes from 0.0125 mag to 0.7 mag, with a particular emphasis on lower values. Overall this constitutes a grid of 8 156 736 models at each redshift.

For modeling the stellar populations we use the GALEX FUV and NUV bands, the SDSS u, g, r, i, and z bands, and the 2MASS J, H, and K_s bands. In addition, we use L_dust, the bolometric dust luminosity, as estimated in Sect. 3.1.1 as a further constraint.

For some of the analysis, we use oxygen abundances, which have been determined from dust-corrected [N II]6584/[O II]3726,3729 emission lines ratio, using the calibration of Kewley et al. (2002). This particular calibration has an advantage over other common calibrations in that it is less sensitive to the effects of the ionization parameter and that it has a unique mapping between the line ratio and the oxygen abundance.

3.1.3. Correction of stellar contamination in MIR bands

Because stellar emission contributes to MIR bands, we have expanded the two-step approach of Salim et al. (2018) described in Sect. 3.1 into a 4-step strategy to reliably correct for this contamination. First, we fit the observed fluxes from 12 μm to 500 μm, as presented in Sect. 3.1.1, in order to estimate the dust luminosity. With this constraint in hand, we fit the stellar populations as described in Sect. 3.1.2, from which we derived the expected stellar fluxes for each galaxy in the WISE 12 μm and WISE 22 μm bands. Stellar contribution is typically small, around a few percent. The last two steps consisted in repeating steps one and two but with the estimated stellar fluxes subtracted from the WISE bands. The final difference in the derived physical properties is small and neglecting this correction would not have affected the results in any substantial way.

3.2. A new approach to the construction and use of dust emission “templates”

Empirical dust emission templates allow the determination of dust luminosities and other dust properties in cases when the wavelength coverage of observations is sparse. Because IR spectroscopy is in general limited to small portions of the IR spectrum and is not widely available, all efforts to construct empirical templates ultimately rely on some theoretical modeling (even if it is just a simple gray body) in order to produce dust emission spectra that are continuous in wavelength. Previous work in this area utilized various combinations of empirical MIR templates and theoretical FIR modeling to achieve this goal. Sophisticated theoretical dust emission models (or model grids), such as the Draine & Li (2007) models used in this work, have the flexibility required for reproducing the important variations in the spectral shape and the strength of the PAH bands seen in star-forming galaxies. Because the model grids are too unconstrained when the SED is not well sampled, and are impossible to use for estimating the dust luminosity in the case when only one IR flux point is available, the purpose of the empirical dust emission templates is to effectively narrow down that parameter space of spectral shapes to what is actually found in galaxies. Templates are typically a family of spectra, discretely dependent on some parameter. The usual choice for this single parameter is either the FIR color (i.e., the dust temperature, Dale & Helou 2002), or the dust luminosity (Chary & Elbaz 2001). The two quantities are correlated, but the use of dust luminosity has the advantage that, being an extensive quantity, it allows a template to be fit and the dust luminosity to be estimated even when only one IR flux point is available. When additional fluxes are available, and are reasonably well separated in wavelength, the family of empirical templates can be used to fit the relative fluxes (i.e., the colors) instead of the absolute flux, conceptually mimicking the usual SED fitting process. In such case, no use is made of the dust luminosity attached to each individual template. We refer to these two ways of estimating the dust parameters from the templates as the traditional approach.

A straightforward method for constructing the templates would be to combine different spectra in bins of a given physical property and average them in some way. This approach has several important drawbacks. The discretization of the templates that results from the binning is somewhat arbitrary, the number of objects in a bin can be highly variable causing nonuniform accuracy, and the stochastic nature of averaging small samples could lead to nonphysical, or at least odd, spectral shapes under some circumstances. Templates have been made to avoid the latter issues by forcing monotonic relations between monochromatic and total luminosity (e.g., Chary & Elbaz 2001), but they are still discrete. In light of these disadvantages, we build our non-discrete templates by producing functional relations that connect monochromatic IR luminosities to any physical property of interest. The IR spectrum corresponding to any value of this physical property can then be obtained by the reciprocal relations. This technique has multiple advantages. First, it allows for the IR spectra to be defined over a continuum of one or more parameters rather than in more or less sparse and arbitrary discrete bins. Since the relations are fit linearly (in logarithm space) to the entire sample, the template spectra iron out the noise from the stochastic diversity of galaxies in any given bin.

Second, while one can export the newly derived IR spectra as discrete templates, and use them in the traditional way described above (essentially, fitting them to one or more bands), our approach makes this extra step unnecessary, because one can derive a desired physical property from the relations directly, as a function of one or more monochromatic luminosities, or even as a function of some additional parameter. In the case of a single flux, the process is conceptually equivalent to, and the results are identical to, the template fitting. When additional fluxes are available, we have found that the functional relations actually provide more accurate estimates of the parameter (specifically, TIR luminosity) than the fitting method. The principal disadvantage of the relations method compared to the fitting method is that the estimation of the uncertainties of the derived parameter (TIR luminosity) would require a Monte Carlo simulation with perturbation of the input fluxes. In what follows, we generally discuss relations, reserving the term templates for the discrete set of SEDs.

In practical terms, the construction of relations starts by adopting, for each galaxy, the rest-frame dust emission spectrum and associated physical properties (both for the dust and the stellar populations) corresponding to the minimum χ² (the best-fitting model). Then for each wavelength λ we compute the relation between each physical property p and the luminosity λL_λ(λ). As the flux density of the emission spectrum of galaxies does not necessarily vary linearly with the physical properties, we have elected to fit simple power laws. In effect, considering the logarithms of the respective quantities, the approach reduces to:

$$\begin{array}{c}\hfill log\lambda L_{\lambda }(\lambda )=\alpha (\lambda )\times logp+\beta (\lambda ),\end{array}$$(1)

with α(λ) and β(λ) the wavelength-dependent coefficients determined from the fits. As we see in Sect. 4, α(λ) is a particularly interesting quantity as it indicates whether the luminosity at wavelength λ changes linearly (α(λ) = 1) or nonlinearly (α(λ) ≠ 1) with p. We will eventually adopt a more complex formulation by considering two parameters simultaneously, but for now we focus on a single-parameter dependence. Naturally, as we see in Sect. 5.2, the reciprocal of Eq. (1) allows us to estimate a physical property p from λL_λ(λ). We use the former formulation to study parameterized IR spectra (Sect. 4) and the latter to provide estimations of the physical properties from one or more IR photometric bands (Sect. 5).

4. Results: physical drivers of the diversity of dust emission spectra

We divide the presentation of our results into two sections. In the current section we explore how IR spectra depend on the physical properties, in particular, dust luminosity, SFR, sSFR, stellar mass, dust mass, and oxygen abundance (gas-phase metallicity). In addition to providing us with a physical insight into what drives the diversity of dust emission spectra, we also investigate how well can each of these properties be determined from different parts of the IR spectrum. Informed by this analysis, in the next section we discuss the estimation of two of these properties (dust luminosity and SFR) using one or more IR bands.

4.1. IR spectra as a function of extensive and intensive properties

In this section we explore relationships between the dust emission spectra and various physical properties, which can be divided into extensive and intensive properties. Extensive properties (SFR, M_star, M_dust, and L_TIR) involve absolute rather than relative fluxes, that is, they directly scale with the extent of the galaxy. Intensive properties (sSFR and oxygen abundance), on the other hand involve normalized or relative fluxes. Even though it is inherently not possible to parameterize dust emission spectra on an intensive property, this becomes possible if they are first normalized by an extensive property. We therefore explore relationships between IR SED normalized by L_TIR and two intensive physical properties.

We present in Fig. 3 the dust emission spectra as a function of the six aforementioned physical properties. The spectra parameterized on extensive physical properties show large absolute variations. First and foremost, their monochromatic luminosities (in log) scale very well with the physical properties, often in a nearly linear way, as we see later. There are more subtle but nevertheless clear relative variations as well, in particular regarding the location of the peak of the modified black body emission, such that it moves to shorter (respectively longer) wavelengths with increasing L_TIR and SFR (resp. M_dust and M_star). The shift is intuitively explained in that higher levels of L_TIR and SFR probably trace more intense radiation fields that heat the dust to warmer temperatures. Conversely, in dustier galaxies the energy is distributed over a larger number of dust grains (everything else being the same) and higher mass galaxies tend to have a softer, less intense radiation field than low mass galaxies, which are optically bluer.

Fig. 3. From the top-left to the bottom-right, dust emission spectra parameterized on LTIR, Mdust, SFR, Mstar, sSFR, and the oxygen abundance. The spectra are calibrated in absolute luminosity for the extensive properties, but they are normalized to LTIR for the intensive properties. The color of each spectrum follows the corresponding physical property as indicated by the color bar above each panel.

Even though the above description provides a satisfactory overall qualitative explanation, the physical properties are not independent from one another and there can be important variations due to secondary parameters. Such variations can actually be seen in the trends of the dust emission spectrum with intensive properties. For instance, a higher sSFR corresponds to a warmer dust temperature with a displacement of the peak toward shorter wavelengths, which is what we would intuitively expect. This strong trend with sSFR shows that neither SFR nor M_star are sufficient by themselves to account for the full range of the dust emission spectral diversity, and indeed the trend with sSFR may be more fundamental (da Cunha et al. 2008; Nordon et al. 2012; Magnelli et al. 2014). Finally, the oxygen abundance has a more moderate impact on the variation of the spectrum. As expected, but still interesting that we see that in these “monotonized” spectra, the galaxies with the higher oxygen abundance have a stronger PAH emission and a slightly colder peak compared to lower metallicity galaxies. However, even at lower oxygen abundance PAH emission remains fairly important. The reason is that we have a dearth of very low metallicity objects in this sample. Indeed the MIR emission of low-metallicity galaxies tends to be more difficult to detect at the shallow depth of large surveys. With over 99.9% of the sample galaxies having an oxygen abundance larger than 8.4, this is consistent with the results obtained for instance by Engelbracht et al. (2005, 2008) with Spitzer.

4.2. How monochromatic IR luminosities scale with physical properties

To gain additional insight into the physical drivers of variations in the dust emission spectra, we show in Fig. 4 how α, the scaling between a monochromatic luminosity and a physical property (see Eq. (1)), varies as a function of wavelength. For extensive quantities, a value of one discriminates between the superlinear (the monochromatic luminosity increases faster than the physical property) and the sublinear (the monochromatic luminosity increases slower than the physical property) regimes.

Fig. 4. Scaling coefficient between the monochromatic luminosity and a physical property versus the wavelength. The scaling coefficient α for the extensive properties LTIR (blue), Mdust (orange), SFR (green), and Mstar (red) is shown in the left panel. The scaling coefficient α′ for the intensive properties sSFR (blue) and the oxygen abundance (orange) is shown in the right panel. For the latter, λLλ has been normalized to LTIR, that is, it shows the residuals at fixed luminosity. The black horizontal lines indicate respectively the threshold between superlinearity and sublinearity (left), and between correlation and anticorrelation (right).

Both L_TIR and the SFR present a remarkably similar behavior over the full wavelength range, showing the strength of the dust emission to estimate the SFR of galaxies. We note that the SFR discussed in this paper is not what is sometimes referred to as the “obscured” SFR, but rather the true, total SFR. Interestingly, the scaling only varies slightly depending on whether the wavelength lies within a PAH line or is situated in the continuum. This constancy helps in cases when the redshift of the source and therefore where the filter is placed in rest frame is not precisely known. At longer wavelengths, the relations become clearly superlinear for both properties, peaking close to 50 μm, before decreasing, becoming briefly linear again around 100 μm. This superlinear range is likely due to the progressively warmer modified-blackbody at higher L_TIR and SFR, which causes a rapid increase of the emission as the peak shifts to shorter wavelengths. Finally, the scaling becomes increasingly sublinear beyond 100 μm. At these wavelengths the emission progressively becomes dominated by the cold dust, which contributes less to L_TIR and is only weakly related to the SFR, as these physical properties more tightly relate to the emission of warmer dust emitting at shorter wavelengths.

Our results regarding the scaling between the monochromatic luminosity and L_TIR are qualitatively similar to those of Rieke et al. (2009), derived for nearby, FIR selected galaxies. They find an α ≈ 0.85 sublinearity at 8 and 12 μm and an α ≈ 1.1 superlinearity at 24 and 60 μm.

The relation based on M_star shows a systematic sublinearity, with α varying typically between 0.75 and 0.95. A sublinearity is expected as higher mass galaxies tend to have progressively redder stellar populations. Therefore an increase in M_star translates into a smaller increase of the energy absorbed and re-emitted by dust. Scaling is closer to linear in PAH lines, possibly due to the role of the interstellar radiation field, and therefore, older stars, in the excitation of these lines (e.g., Haas et al. 2002).

The case of M_dust is more interesting in that α shows a stronger dependence on wavelength than M_star. Up to around 200 μm, α is lower for M_dust than it is for M_star. Dust dominating the emission at these wavelengths represents only a minor fraction of M_dust. There is a change in the regime at longer wavelengths, with α progressively pushing above 0.9 in the submillimeter, reflecting the fact that the bulk of the dust is cold and dominates the emission in the Rayleigh-Jeans regime. This forms the basis for using a single submillimeter band as a measure of the dust mass, and by extension of the gas mass too (e.g., Dunne & Eales 2001; Groves et al. 2015; Scoville et al. 2016; Millard et al. 2020).

In the case of intensive properties, we have normalized the spectra to L_TIR, which means that the scaling coefficient, which we denote α′ to differentiate it from α that is used in the case for extensive properties, probes variations of the shape of the emission spectra, and not their absolute normalization. Relations with α′(λ) = 0, would indicate the normalized emission spectra are invariant at wavelength λ with the considered physical property, that is, no dependence between the normalized monochromatic luminosity and the intensive property. Figure 4 shows that the sSFR and the oxygen abundance are almost perfectly anticorrelated, with their α′(λ) of opposite signs at almost all wavelengths. This behavior is probably the consequence of the anticorrelation between the oxygen abundance and sSFR at a fixed stellar mass, the extension of the mass–metallicity relation (Ellison et al. 2008; Salim et al. 2014). There are however two regions where α′(λ) ≃ 0, slightly below ∼20 μm and around ∼90 μm. With the dust emission being largely independent from these physical properties at these wavelengths, the monochromatic dust luminosity and L_TIR should be in linear relation with one another. This is consistent with what we found previously. Between these two wavelengths, α′(λ) > 0 (α′(λ) < 0) for the sSFR (oxygen abundance), which translates into the increase (decrease) of dust temperature with respect to the sSFR (oxygen abundance). Conversely, at shorter wavelengths we see the strengthening of PAH features with increasing metallicity and their slight decrease with the sSFR. At long wavelengths, the submillimeter emission is strongly dependent on metallicity, probably as a consequence of an increase of the gas-to-dust mass ratio, whereas it weakens at high sSFR as a larger fraction of the luminosity is emitted at shorter wavelengths.

4.3. How monochromatic IR luminosities correlate with physical properties

Besides the linearity, another important quantity to consider is the scatter around the relation, that is, the degree of the correlation between monochromatic luminosities and some physical property. Indeed, a relation that is linear but presents a very large scatter may ultimately not be sufficiently reliable for individual objects and a nonlinear relation but with a smaller scatter will be more desirable for parameter estimation. To investigate this aspect, we present in Fig. 5 the scatter of the residuals around the relations versus the wavelength.

Fig. 5. Standard deviation of the residuals around the relation between the monochromatic luminosity and a physical property (Eq. (1)) versus the wavelength, for the extensive properties LTIR (blue), Mdust (orange), SFR (green), and Mstar (red).

First, it appears that in addition to being strongly nonlinear, the M_star relation presents a high standard deviation that is typically larger than 0.3 dex. Interestingly, the scatter is smallest in the MIR and then at λ > 100 μm, lending additional support to the notion that the two regions are governed by interstellar radiation field heating (Haas et al. 2002). Up to ∼100 μm the M_dust relation has a behavior that is very similar to that of M_star. However, for M_dust the standard deviation keeps decreasing with increasing wavelength, reaching less than 0.1 dex at 1 mm, confirming once again the reliability of a single long wavelength band to estimate the dust mass in star-forming galaxies, and showing that the longer the wavelength the better the estimator of M_dust it is.

By far, the L_TIR shows the smallest standard deviation of all extensive quantities, with a minimum of slightly more than ∼0.05 dex (corresponding to a relative scatter of ∼13% in linear scale) at ∼95 μm. In the 20−40 μm range however the standard deviation is higher, from 0.15 dex to 0.16 dex. Even though on average the monochromatic luminosity and L_TIR are in linear relation with one another in the MIR, the bulk of L_TIR is emitted around the peak of the modified black body. Thus we can expect that variations in the dust emission emerging in the MIR would have a relatively smaller influence on L_TIR, leading to an increase of the standard deviation of the residuals in this spectral region. It is interesting to note that the region between 20 μm and 50 μm actually has a higher standard deviation than the region between 10 μm and 20 μm, invalidating a simplified notion that going to longer wavelengths toward the peak is always preferred. The reason for this may lie in the prominent role of stochastically heated very small grains, the emission of which peaks in this region (Desert et al. 1990). Only beyond 50 μm does the standard deviation diminish as one moves toward the peak of the emission spectrum, dominated by large grain emission. This consideration may be relevant when designing IR detectors. Also, it is worth pointing out that the region beyond ∼200 μm is actually inferior to MIR, potentially informing the use of ALMA vs. JWST. Altogether, simultaneously taking into account the linearity, the lack of dependence on sSFR, and the small standard deviation, suggests that, even though no wavelength is flawless, the emission of the dust around 90−100 μm and L_TIR trace each other best.

Compared to the L_TIR relation, the SFR relation shows a scatter higher by 0.03 dex to 0.05 dex in the MIR, a region for which we showed the relations are close to being perfectly linear. Overall, the best monochromatic luminosities for estimating SFR are the same that best determine L_TIR (∼18 μm and ∼100 μm). There has been some debate in the literature whether a monochromatic luminosity is a better tracer of SFR than L_TIR.

Since SFR is the total of optically obscured and unobscured contributions, the larger scatter in comparison to L_TIR may be due to the diversity of attenuation properties across the sample (e.g., amplitude of the attenuation, shape of the attenuation curve, differential reddening between different stellar populations) as well as the range of relative contributions of young and old stellar populations to dust heating. The limitations in the ability to estimate the total SFR from dust emission alone could be lifted by using hybrid SFR estimators, which combine IR and UV luminosities (e.g., Elbaz et al. 2007; Daddi et al. 2007; Hao et al. 2011; Liu et al. 2011; Kennicutt et al. 2009; Boquien et al. 2016), a topic that we will explore in detail in a future publication, or more generally by performing SED fitting utilizing the full UV through IR SED.

The overall analysis in the current and the preceding sections provides more detail to our current understanding of the physical drivers of the dust emission, confirming the close connection of the MIR and 80−100 μm spectral regions with both L_TIR and the SFR. The former connection is important in particular for JWST, which will only be able to target the MIR for galaxies up to z ∼ 2. The analysis also highlights the sensitivity of dust parameter estimation on sSFR, which we will address shortly. Finally, we confirm that the submillimeter emission is a reassuringly reliable tracer of the dust mass.

5. Results: estimation of TIR luminosity and SFR from IR photometry

5.1. Parameterization of dust emission templates and relations

A parameterization on an extensive quantity has the added benefit that it makes it possible to derive L_TIR using a single flux point or band. We point out that this extensive quantity does not need to be L_TIR. In particular, it may be useful to parameterize the templates on the total SFR, which, like L_TIR displays a close relation to the shape of the IR spectrum (Fig. 3), and which in our modeling is known from the UV-to-IR SED fitting. While the unobscured star formation is not directly observable in the IR, it is in some instances possible to estimate it reasonably well without the UV data. Furthermore, parameterizing directly on SFR has the advantage of not having to use a fixed conversion factor to translate L_TIR into SFR.

Previous studies have emphasized either the total luminosity (e.g., Chary & Elbaz 2001) or the sSFR (e.g., da Cunha et al. 2008) as the principal drivers of the shape of the dust emission spectrum. What Fig. 4 shows is that, except for particular wavelengths, both are important. From the right panel of Fig. 4 we see that the change in the sSFR of one dex at fixed L_TIR results in the difference in the luminosity at 8 μm of nearly 0.2 dex. If the templates (or relations) only depended on L_TIR, the estimates of the total IR luminosity produced by them would be biased for galaxies with atypically high or atypically low levels of star formation, potentially leading to the systematics when our estimators that are constructed based on low redshift galaxies are applied to objects at higher redshifts where the average sSFR is higher. This motivates us to present our templates and relations parameterized on sSFR as well. In Fig. 3 we see that the peak of the emission shifts from 120 to 70 μm as sSFR increases. This range encompasses the peaks of stacked SEDs going from z ∼ 0 to z ∼ 2 (Béthermin et al. 2015), giving us confidence that by incorporating the sSFR dependence we are effectively producing templates that are applicable to a wide range of redshifts. On the other hand, in the presence of multiple flux points, the color term implicitly accounts for the sSFR dependence, and including the sSFR explicitly in the relations is superfluous, as verified by the tests that we carried out.

Two-parameter templates covering FIR (30−1000 μm) and based on spectra of simulated galaxies have been previously introduced by Safarzadeh et al. (2016). However, their second parameter is M_dust, which cannot be well constrained in the absence of submillimeter data, which are often not available. Magdis et al. (2012) and Safarzadeh et al. (2016) stress the importance of L_TIR/M_dust as a driver of the IR SED shape. Being qualitatively similar to L_TIR/M_dust, our sSFR dependence confers similar benefits, but is more accessible, requiring essentially only the knowledge of M_star (Sect. 5.4).

To conclude, in this work we produce the following four types of average dust emission spectra (templates): templates parameterized on L_TIR, templates parameterized on total (obscured plus unobscured) SFR, templates parameterized on L_TIR and sSFR simultaneously, and templates parameterized on total SFR and sSFR simultaneously.

As pointed out in Sect. 3.2, we argue that the estimation of the physical properties can be more easily carried out by direct, continuous relations rather than the templates. Therefore, we also construct the following relations for estimating L_TIR and total SFR: relations for single-band flux measurements (Sects. 5.2 and 5.3), sSFR-dependent relations for single-band flux measurements (Sect. 5.4), and relations for multiple-band flux measurements (Sect. 5.6).

Finally, we also provide a tool within the software package to generate any set of discrete templates that the user may find useful, specified by L_TIR or SFR, with or without the additional dependence on sSFR (Sect. 5.7).

5.2. Estimation of L_TIR and the SFR from a single IR band

In this section we focus specifically on estimating L_TIR and SFR from the emission in a single band and not taking into account the sSFR dependence. This is equivalent to using templates parameterized on L_TIR or SFR. In our approach, the relations used for estimating L_TIR and SFR are precomputed for particular bands. Specifically, we have selected a set of representative bands that are or will likely be extensively used for constraining the dust emission of galaxies: WISE 12 μm and 22 μm, Spitzer 8.0 μm and 24 μm, Herschel 70 μm, 100 μm, and 160 μm, and JWST 7.7 μm, 10.0 μm, 11.3 μm, 12.8 μm, 15.0 μm, 18.0 μm, 21.0 μm, and 25.5 μm. The fluxes in these bands are computed by integrating each best-fit template through the filter bandpasses. All the analysis in this subsection is performed assuming z ≈ 0, but in the accompanying software package we provide estimators up to z = 4 so that the estimation can be performed without an explicit K-correction. We present how to take variations of the sSFR into account in Sect. 5.4.

We determine the relation by inverting the dependent and independent variables of Eq. (1) and repeating the fitting procedure:

$$\begin{array}{c}\hfill logp=m(b)\times log\lambda L_{\lambda }(b)+n(b),\end{array}$$(2)

with p the physical property to be estimated (L_TIR or the SFR), λL_λ(b) the luminosity in a photometric band b of pivot wavelength λ, and m(b) and n(b) the coefficients obtained from the fit for that band. The luminosities are computed by integrating each spectrum through the corresponding bandpasses.

As an illustration, we show in Fig. 6 the estimates of L_TIR and the SFR based on Eq. (2) and three representative bands, Spitzer 8 μm, JWST 18 μm, and Herschel 100 μm. The data come from our sample used to derive the relations. As expected, all three bands correlate well with both L_TIR and the SFR, however with visibly less scatter for the former, in agreement with analysis presented in Fig. 5. Some reduction in scatter, depending on the quality of the available sSFR estimates, is possible when using sSFR-dependent relations (Sect. 5.4). In Table 1 we provide the coefficients m and n to estimate L_TIR and the SFR for all the bands mentioned previously. These coefficients are valid only for redshifts close to zero. For use at other redshifts we direct the user to the provided software tool (Sect. 5.7). This provides in effect estimators that are continuous in redshift, similar to the MIR SFR estimators of Battisti et al. (2015).

Fig. 6. Estimation of LTIR (left) and the SFR (right) from the Spitzer 8 μm (blue), JWST 18 μm (orange), and Herschel 100 μm (green) bands. The dots represent our H-ATLAS sample and the solid lines the best-fit relation computed from Eq. (2).

We see that in line with our findings presented in Sect. 4.2, the emission around 100 μm is one of the best tracers for both L_TIR and the SFR. Surprisingly, most estimators are sublinear (m < 1), even though the reverse relation is also sublinear (α < 1). Upon closer inspection, the reason for this counterintuitive behavior is because the variance for each band is larger than the covariance with L_TIR or the SFR. We note that both in Sect. 4.2 and here we are using an ordinary least-square fit. In essence, in each case we know the dependent variables perfectly as we rely entirely on the best-fitting models for the emission spectra and the associated physical properties, and depending on the case, the objective is indeed to minimize the scatter either for the emission spectrum (Sect. 4.2) or the physical properties (current section). It is important to keep this aspect in mind when considering the physical interpretations of our results.

5.3. Estimation of M_dust from a single IR band

Another important physical property to estimate is M_dust, which can serve as a proxy for the total gas mass, while also providing constraints on the dust production and destruction processes. We saw in Sect. 4 that the estimation of M_dust at longer wavelengths has an excellent potential, with α approaching 1, while the standard deviation of the residuals are progressively dropping with wavelength. We see in Fig. 7 that there is indeed a reasonably tight relation between M_dust and the luminosities in the Herschel SPIRE bands. We do not consider here the Herschel PACS bands because the shorter wavelengths are also sensitive to warmer dust and do not provide sufficiently precise results (Fig. 5). We give the coefficients for M_dust estimation using Eq. (2) and Herschel SPIRE bands at z ∼ 0 in Table 2.

Fig. 7. Estimation of Mdust from the Herschel 250 μm (blue), 350 μm (orange), and 500 μm (green) bands. The dots represent our H-ATLAS sample and the solid lines the best-fit relation computed from Eq. (2).

All SPIRE bands are almost perfectly linear estimators of M_dust. The standard deviation of the residuals goes from 0.19 dex at 250 μm down to 0.12 dex at 500 μm, an improvement of slightly more than 30%. One caveat to keep in mind using this estimator is that contrary to L_TIR or the SFR, M_dust is sensitive to the details of the underlying dust models, such as the emissivity index. Assumptions differing from that of the Draine & Li (2007) and Draine et al. (2014) models, including, for example, the choice to model carbonaceous grains as amorphous rather than graphite (Schreiber et al. 2018), may yield systematic differences in the resulting M_dust.

5.4. Estimation of L_TIR and the SFR from a single IR band and the sSFR

Even though single-band estimators can provide us with good estimates for low-redshift samples on average, they may suffer from biases in some regions of the parameter space, in particular for galaxies that have atypically low or high sSFR at fixed L_TIR (Sect. 5.1), or for general population of galaxies at higher redshifts. For those cases it is recommended to use sSFR-dependent templates and relations.

Following the approach we used in Sect. 5.2, we have computed estimators as a function of one IR band and the sSFR:

$$\begin{array}{c}\hfill logp=m(b)\times log\lambda L_{\lambda }(b)+n(b)+s\times log\mathrm{sSFR}.\end{array}$$(3)

We give the corresponding coefficients in Table 3.

Two quantities are of particular interest, the coefficient s, which scales log sSFR, and σ, the standard deviation of the residuals. As for α′, a value of s close to 0 indicates that the estimator is largely independent from sSFR. For L_TIR this is in particular the case of the MIR bands that do not cover prominent PAH features, or the FIR around 100 μm, close to the peak of the emission. Conversely, bands with strong PAH features and long wavelength emission beyond the peak show the strongest dependence. The reason for the stronger dependence for bands overlapping prominent PAH features is not entirely clear as we would expect the PAH and L_TIR emission to scale linearly with each other over a fairly large range of radiation field intensities (see for instance Fig. 15 of Draine & Li 2007). However this could be an indirect effect of the PAH abundance, since in our sample galaxies with a higher oxygen abundance tend to have a lower sSFR, or the consequence of the role of older populations in heating both the PAH features and large grains (see Sect. 4.2).

The standard deviations of the residuals in Table 3 show systematically reduced values with respect to Table 1. Unsurprisingly the largest reductions correspond to the largest values of s. The relative improvement is greater for the estimation of SFR than L_TIR because the proportionality between these two quantities is itself dependent on sSFR.

The sSFR is a quantity that can be challenging to evaluate, ideally requiring SED modeling. The availability of the sSFR would also in many cases eliminate the need to separately estimate L_TIR or the SFR. However, uncertainties on the sSFR are generally small compared to the full dynamical range of the sSFR, and the shape of the templates varies slowly and monotonically with the sSFR. This means that even an imprecise estimate of the sSFR (or one obtained indirectly using some population age estimate, e.g., the D4000 index or Hα equivalent width) would still prove highly useful to reduce or even eliminate possible biases. As a matter of fact, as long as the stellar mass is available, it is possible to use the SFR estimators iteratively. For instance, the initial SFR could be estimated from the relations given in Table 1. Combined with M_star it would yield an estimate of the sSFR. Then it would be possible to apply the relations provided in Table 3. Applying these simple steps appears to be sufficient to eliminate the dependency with the sSFR. As an example, for the Spitzer 8 μm band, the R² coefficient of the difference between the exact L_TIR from the model and the estimated L_TIR from the estimator with respect to the sSFR goes from 0.1787 to 0.0197. Even though the determination of the SFR is not ideal as the relation itself shows one of the strongest dependencies on the sSFR among all the bands, the approximation remains sufficient to obtain satisfactory results to eliminate systematic biases.

5.5. Explicit estimation of SFR and L_TIR from JWST 21 μm observed at various redshifts

JWST’s Mid-Infrared Instrument (MIRI, Rieke et al. 2015) will be a workforce to estimate the SFR of galaxies from the local universe up to the cosmic noon (e.g., the SFR estimators of Senarath et al. 2018). The software tool that accompanies this paper allows the calculation of SFR and L_TIR from any JWST flux for arbitrary redshift. Nevertheless, to facilitate quick calculation, in this section we provide simple, redshift-dependent formulas to derive SFR and L_TIR from F2100W fluxes. We estimate these quantities at different redshifts for a range of fluxes using sSFR-dependent relations, where we assume the sSFR that corresponds to the main sequence at that redshift at log M_⋆ = 10.5, according to Speagle et al. (2014). In order to simulate the effect of the intrinsic width of the main sequence and the fact that the stellar mass may be different from log M_⋆ = 10.5, we perturb sSFR by 0.5 dex 1-σ Gaussian. Figure 8 shows the dependence of SFR on the observed flux in five redshift bins. The scatter is the result of the assumed scatter in sSFR which makes it possible for the galaxies with the same flux and the same redshift to have different SFR or IR luminosities.

Fig. 8. Estimated total (obscured plus unobscured) SFR as a function of the flux measured in JWST MIRI 21 μm filter (F2100W) and redshift, based on the relations that take into account the evolution of sSFR. The scatter results from assuming a 1-σ scatter of 0.5 dex in sSFR of actively star-forming galaxies at each redshift.

If stellar masses (and therefore sSFR) of galaxies are unknown, one can use the relations to estimate SFR and IR luminosity that assume the sSFR evolution from Speagle et al. (2014):

$$\begin{array}{cc}& log\mathrm{SFR}=0.83log{F}_{F2100W}+2.23logz+2.19,\sigma =0.19,\hfill \end{array}$$(4)

$$\begin{array}{cc}& log{L}_{\mathrm{TIR}}=0.92log{F}_{F2100W}+2.19logz+12.23,\sigma =0.14.\hfill \end{array}$$(5)

If the stellar masses are known, SFR can be estimated using the above relation and the sSFR-dependent relations provided below can then be used to refine the result.

$$\begin{array}{cc}& log\mathrm{SFR}=0.84log{F}_{F2100W}+1.83logz\hfill \\ \hfill & +0.31log\mathrm{sSFR}+5.00,\sigma =0.10,\hfill \end{array}$$(6)

$$\begin{array}{cc}& log{L}_{\mathrm{TIR}}=0.92log{F}_{F2100W}+2.00logz\hfill \\ \hfill & +0.15log\mathrm{sSFR}+13.57,\sigma =0.12.\hfill \end{array}$$(7)

Estimates of the SFR especially benefit from the inclusion of the sSFR term because the relationship between SFR and L_TIR is to first order dependent on the sSFR.

To get the full error on estimated quantities, the relative flux errors should be added in quadrature to the above uncertainties for the relationships. The relations should not be used above z = 2.2 where the rest-frame wavelength covered by F2100W becomes rather short. Our estimates for the relation between F2100W and L_LIR agree well with those of Schreiber et al. (2018): at log L_TIR = 11 the difference is less than 0.1 dex at z < 0.5 and less than 0.2 dex at 1 < z < 2. The principal difference is that we do not assume a linear relation between monochromatic and total dust luminosities. Furthermore, we provide specific relations for SFR which are even more sublinear than the ones for L_TIR.

5.6. Estimation of L_TIR and the SFR from multiple IR bands

Whenever possible, it is desirable to include information provided by several bands to determine L_TIR or the SFR. The fundamental reason being that having more bands provides additional constraints on the shape of the emission spectrum. Following this idea, multiple derivations have been provided in the literature with excellent results using a diverse combination of bands from IRAS (Sanders & Mirabel 1996) to Spitzer (Dale & Helou 2002; Boquien et al. 2010) and Herschel (Boquien et al. 2011; Galametz et al. 2013). We limit our relations to a maximum of four bands, yielding a total of 4047 band combinations. The fitting procedure is similar to that presented in Sect. 5.2 but extended to multiple bands:

$$\begin{array}{c}\hfill logp={\displaystyle \sum _i{m}_i\times log{\lambda }_i{L}_{\lambda }(b_i)+n,}\end{array}$$(8)

with each scaling coefficient m_i corresponding to band b_i, i being the index of the band. Given the particularly large size of the resulting table, we provide these coefficients in an electronic form only.

We must mention, however, that this derivation is made with a small but important modification with respect to the single-band case. Previously it was not necessary to impose any bound on m or n. When considering the case of Eq. (8), a priori nothing would restrict m_i to be negative, which could intuitively be understood as color terms. However, bands providing very similar information are degenerate. This is often the case of bands with close wavelengths, for instance JWST 21 μm and WISE 22 μm. Without imposed bounds, the two m_i will be of opposite sign and have similar and very high absolute values (for instance −50 and +50). This means that in practical cases even a small amount of noise will be amplified to a considerable degree, strongly perturbing the estimate. To address this issue, we have imposed that all m_i are bounded between 0 and 2. Analysis of these estimators with synthetic catalogs injected with noise shows that even though the residuals around the calibration sample are slightly higher when the m_i coefficients are bound, they provide us with much more reliable estimates as they are more resilient to the photometric noise.

Throughout this section we primarily focused on relations for galaxies at z ∼ 0. The behavior of the coefficients in the scaling relations at different redshifts in the case of a single band is presented in Appendix A.

5.7. Data products and software

Given the wealth of data contained in the relations we have presented, their practical use could be a challenge. To address this, we provide a number of products to the community on the web site that hosts GSWLC⁶, the sample upon which this study is based.

First, following the traditional approach, we provide a grid of templates parametrized on the physical properties as described in Sect. 5.1, both as ASCII and FITS files. Both for flexibility and ease-of-use, we also provide a software tool to generate templates for any value of the physical properties and also for any combination of L_TIR and sSFR as a second parameter. An interactive visualization tool of the templates based on these two parameters is also provided.

Furthermore, we provide a tool to estimate L_TIR and SFR from an input table of fluxes and redshifts for any combination of up to four bands. For estimates based on single bands, a value that takes into account sSFR dependence is also provided. The fitting coefficients given in Table 1 along with the corresponding covariance matrices for any combination of up to four bands from z = 0 to z = 4 with steps of 0.01 are provided with this tool.

6. Discussion

6.1. Comparison with previously published templates

In this section we compare our templates to four previously published sets of templates: Chary & Elbaz (2001), Dale & Helou (2002) (in particular the Dale et al. 2014 update), Rieke et al. (2009), and Smith et al. (2012) templates updated to include Herschel PACS bands of a larger sample than in the original paper⁷. Except for Smith et al. (2012), the other templates were derived or constrained using samples, often heterogeneous, selected from shallow FIR surveys. Details of the construction of these four template sets are given in Appendix B. We limit our comparison to template sets that have been parameterized on an extensive quantity (L_TIR).

A straightforward direct comparison between these templates is not always possible as L_TIR may not have been computed over the exact same wavelength range, the templates are generally defined on a different discrete grid of L_TIR, or the templates may not have been parametrized on L_TIR to begin with, as is the case of the Dale & Helou (2002) templates. In order to provide a fair comparison we have therefore taken a number of steps to homogenize these template sets. First, we have parametrized the Dale & Helou (2002) templates against L_TIR by using the relation of Marcillac et al. (2006), which links L_TIR to the 60-to-100 μm ratio, in combination with Table 2 from Dale & Helou (2002), which links this ratio to α_SF, the intrinsic parameter of the templates (not to be confused with the α(λ) coefficient used in the present article). Then we have recomputed the value of L_TIR of each set of templates to correspond to the integral of the emission in the 8 μm to 1 mm wavelength range. Finally, we built an interpolator for each set of templates so that a spectrum could be computed for any L_TIR within their range of validity, enabling different templates to be compared for the same values of L_TIR.

We present in Fig. 9 the comparison of our L_TIR dependent spectral templates to the previously published ones. Differences with respect to other templates, even if small at some wavelengths, are not random, but rather systematic and therefore require some discussion. First of all, the PAH emission in our templates appears to be consistent with what was determined by Chary & Elbaz (2001), Dale & Helou (2002), Dale et al. (2014), and Rieke et al. (2009). Moving to longer wavelengths, the different templates have clear discrepancies starting around 20 μm and extending up to about 70 μm, the region associated with the emission from very small grains (Desert et al. 1990). More specifically, templates other than Smith et al. (2012) tend to show higher emission levels and this excess is luminosity dependent, with the largest discrepancies found for the largest L_TIR. The origin of these differences is not certain. In the case of the Dale & Helou (2002) templates, Fig. 10 from Dale et al. (2001) and Fig. 4 from Dale & Helou (2002) suggest that there are systematic residuals between their observations and their models below 88 μm. Such residuals appear consistent with the difference we see here between their set of templates and our own work.

Fig. 9. Comparison of the templates derived in this work with those from Chary & Elbaz (2001) (top left), Dale & Helou (2002), Dale et al. (2014) (top right), Rieke et al. (2009) (bottom left), and Smith et al. (2012) (bottom right). The solid lines correspond to the new templates and the dashed lines the literature templates. Each pair of lines of the same color corresponds to log LTIR/L⊙ of 9.0, 9.5, 10.0, 10.5, 11.0, and 11.5, following the bar above each panel. Spectra outside of the definition range of a given set of templates have been omitted. The Chary & Elbaz (2001) and Smith et al. (2012) templates include stellar populations, precluding any comparison at short wavelengths where such populations dominate.

More generally, relatively few normal galaxies have been observed in the spectral region between the mid and the far IR. IRAS and ISO had a gap between 25 μm and 60 μm. Likewise Spitzer also had a gap between 24 μm and 70 μm in terms of broadband observations and Herschel had no capability below 70 μm. The IRS spectrograph on-board Spitzer had a coverage up to 40 μm and was used by Rieke et al. (2009). However this was only done on a small number of objects. The LIRG and ULIRG templates may have a markedly different spectral shape compared to the galaxies in our sample and for galaxies with log L_TIR/L_⊙ < 11 their templates are based on four noise-free IRS spectra from Smith et al. (2007) parametrized on the ratio of the 12 μm and 25 μm bands. This may explain why we observe a significant deviation from log L_TIR/L_⊙ = 11. More generally, Dale et al. (2009) published IRS spectra of nuclear and extranuclear regions of SINGS galaxies (Kennicutt et al. 2003), showing that there is a clear diversity in the dust emission spectra in this wavelength range. Overall, it is likely that the relative lack of strong constraints makes models and templates uncertain over this domain. Furthermore, some templates show abrupt changes with L_TIR. This is the case in particular of the templates of Rieke et al. (2009), where the peak of the emission shifts rapidly toward shorter wavelengths for log L_TIR/L_⊙ increasing from 10.5 to 11.5 (see in particular their Fig. 6). Conversely, in our templates the transition occurs at an almost constant rate of ∼ − 11.5 μm dex⁻¹ in L_TIR. This abrupt change may be the result of different methodologies applied by Rieke et al. (2009) for galaxies above and below log L_TIR/L_⊙ = 11. Namely, their high-luminosity templates use a modified black body to model the FIR regime, whereas low-luminosity templates are based on Dale & Helou (2002) templates in that region (see Appendix B for details). More generally, what we see is that at fixed L_TIR, the templates of Chary & Elbaz (2001), Dale & Helou (2002), Dale et al. (2014), and Rieke et al. (2009) have a bluer (warmer) peak than our models, which may be the result of their basing of templates on FIR-selected samples from shallow surveys which would have preferentially detected galaxies with warmer dust. These elements point towards a strong influence of a small number of extreme objects at high sSFR for the most luminous objects in the templates of Chary & Elbaz (2001), Dale & Helou (2002), Dale et al. (2014), and Rieke et al. (2009). Conversely, our large sample containing more moderate galaxies should not be affected by a small number of extreme objects. In order to test the influence of the sSFR at fixed L_TIR, we show in Fig. 10 our templates parameterized on both L_TIR and the sSFR. We see that higher sSFR templates tend to have a stronger emission below 90 μm. This is qualitatively in agreement with the aforementioned literature templates and confirms the effect of a few extreme objects. The availability of our templates parametrized on both L_TIR and sSFR makes them highly flexible and applicable beyond just normal star-forming galaxies, as we discuss in more detail in Sect. 6.3.

Fig. 10. Comparison of the templates parametrized on LTIR at log LTIR/L⊙ = 9.0, 9.5, 10.0, 10.5, 11.0, and 11.5 (black lines), and the templates parametrized on both LTIR and sSFR for sSFR/yr−1 = −11, −10.5, −10, −9.5, and −9. Templates corresponding to a high sSFR have a stronger emission below 90 μm.

Smith et al. (2012) templates, which are based on H-ATLAS, as for our sample, but without MIR photometry from WISE, show up to a ∼0.5 dex lower emission in the 20−40 μm range compared to our new templates, whereas they have a stronger emission in the Rayleigh-Jeans regime. The templates of Smith et al. (2012) are in stark contrast with the other sets of templates, which present a reasonably good agreement with each other at longer wavelengths. The large difference in the warm dust emission with respect to previous templates was already noted by Smith et al. (2012). A possibility is that the lack of data below 100 μm for most galaxies in their sample coupled with modeling using MAGPHYS, did not allow for the spectra to be meaningfully constrained during the SED modeling process. MAGPHYS (da Cunha et al. 2008) models the IR emission through a number of independent isothermal modified black bodies, whereas each Draine & Li (2007) model includes a specific temperature distribution (see Appendix B). It follows that even in the absence of observational constraints over a certain wavelength range, indirect constraints are provided through observations at other IR wavelengths. The discrepancy of Smith et al. (2012) with respect to other templates at longer wavelength was not readily apparent previously, but because here we compare spectra at identical L_TIR, a discrepancy in a given wavelength range has to be compensated with an opposite discrepancy in another range to maintain the respective luminosities equal. Ultimately, our modeling with CIGALE that combined the Draine & Li (2007) models and WISE bands at 12 μm and 22 μm, provides stronger constraints over this wavelength range.

Overall, we are confident that this new set of templates provides an excellent characterization of normal star-forming galaxies, and when sSFR dependent templates are considered, this characterization encompasses a full range of star-forming galaxies, from relatively quiescent to intensely star-forming, including those that are considered “normal” at higher redshift.

6.2. Comparison of L_TIR with previously published relations

In Sect. 6.1 we compared different sets of templates to one another, without applying them to any observed data set. Here we focus on the KINGFISH sample of nearby star-forming galaxies, for which Hunt et al. (2019) derived L_TIR with CIGALE and data from 32 bands from the FUV to 850 μm. We compare this reference L_TIR to the estimates based on our own single-band relations and on single-band relations from Galametz et al. (2013). For our estimates we do not take sSFR dependence into account because the KINGFISH sample and our H-ATLAS sample have a similar range of sSFR (Kennicutt et al. 2011). We adopt the photometric fluxes and distances of Dale et al. (2012, 2017).

The agreement of the 70 μm and 100 μm estimates with the actual L_TIR values is nearly perfect. The agreement for the 70 μm estimate is particularly important since we did not have observations probing this wavelength. The difference is small at 160 μm with perhaps a small offset of ∼0.1 dex for our estimator with respect to actual values at lower L_TIR. The main difference occurs for the Spitzer 24 μm band, in particular at lower luminosities where L_TIR from our estimators are in line with observations, whereas using the estimator of Galametz et al. (2013) gives results ∼0.2 dex lower. An important point to keep in mind is that the definition ranges for these relations are different. Being based on a nearby sample, the Galametz et al. (2013) relations extend to fainter objects that are out of reach of large area surveys. This difference at 24 μm may suggest that there is an important change of regime in the MIR emission in fainter objects that is not apparent at longer wavelengths, which may generate a difference between the estimators at this wavelength.

If we now compare L_TIR estimated using the relations given in Table 3 to estimates using the previously published templates, we find the mean differences to lie typically within 0.1 dex. However, in the case of the templates of Smith et al. (2012) the mean difference reaches ∼0.25 dex at MIR wavelengths and goes even higher in some luminosity regimes. For instance, our estimates are consistent with those obtained with the templates of Rieke et al. (2009) in the JWST 10 μm filter up to log L_TIR/L_⊙ = 11, but they abruptly deviate beyond this point. Multiple examples of such discrepancies in certain luminosity regimes can be found when comparing to different templates, translating some of their intrinsic characteristics probably originating from the techniques and the smaller samples used to build them. This emphasizes the importance of having estimators applicable over a broad range in L_TIR, as we see in the next section.

6.3. Applicability of relations and templates to diverse galaxy populations

By construction, when applied to our full sample, the relations should recover the ground truth L_TIR without a bias. We confirm that for various single-band and multiband estimates the systematic offset is always smaller than 0.01 dex. This is, of course, the most “favorable” dataset to test for systematic errors. Analysis in Sect. 6.2 demonstrated that the luminosities of KINGFISH galaxies are recovered reasonably well using single-band relations. A more stringent evaluation of the applicability of our templates and relations would be for the data sets that contain galaxies significantly different from the typical galaxies in our sample, either in terms of having L_TIR outside of the range of our sample, or by having different M_dust, even if L_TIR (and even sSFR) are the same, as seems to be the case with local starbursts. For the results obtained based on single-band fluxes, we focus on sSFR-dependent relations as they produce more accurate results in all cases.

6.3.1. Dusty galaxies in the green valley

While the sample from which we derive the relations consists overwhelmingly of actively star-forming galaxies with logsSFR/yr⁻¹ > −11 (Fig. 1), we have confirmed that the relations are also mostly applicable for non-AGN galaxies in our H-ATLAS sample that lie in the green valley region (−12 < logsSFR/yr⁻¹ < −11). Their L_TIR estimated from single-band fluxes at 12, 24, and 100 μm, and utilizing sSFR-dependent relations, shows no bias with respect to the reference L_TIR from full seven IR fluxes CIGALE fits using the Draine et al. (2014) grids.

6.3.2. Low-luminosity galaxies

We next analyze the performance of our relations using a sample of very nearby galaxies, which therefore includes galaxies well below the IR luminosity range spanned by our sample from which the relations were constructed (9 < log L_TIR/L_⊙ < 12). Specifically, we utilize Spitzer observations in 8, 24, 70 and 160 μm bands from the Local Volume Legacy (LVL) survey (Dale et al. 2009). LVL contains 258 galaxies, of which 195 have detections in all four bands. Whereas the IR luminosities of these LVL galaxies span a very wide range (6 < log L_TIR/L_⊙ < 11), 3/4 are below our sample range (log L_TIR/L_⊙ < 9). We derive four single-band estimates of L_TIR based on different bands, various estimates including two bands, and one estimate based on all four bands. For the single bands, we estimate L_TIR using sSFR-dependent relations, where sSFR is based on the SFR estimated from the relation itself (value after the first iteration), and the stellar mass, taken from Cook et al. (2014), is based on the 3.6 μm luminosity. At 70 μm and 160 μm we used PACS instead of MIPS bandpasses, which has no appreciable effect on the results. We compare our L_TIR estimates to the reference values that we obtain by fitting all four bandpasses with full Draine et al. (2014) grid, with identical range of grid parameters as used to fit our H-ATLAS sample in Sect. 3.1.1. Our reference IR luminosities agree (0.01 dex offset, 0.04 dex standard deviation) with the ones given by Dale et al. (2009), which they derive using a relation between MIPS bands and L_TIR from Dale & Helou (2002).

The performance of our relations for LVL is presented in Fig. 11. For the luminosity range for which the relations are constructed (L_TIR > 9) all our single-band estimates of L_TIR are in excellent agreement with the reference values, with offsets no greater than 0.04 dex, and a scatter smaller than 0.13 dex. If we were to use relations without the sSFR dependence, the offset would in some bands exceed 0.1 dex. In the low-luminosity regime, the single band estimates are reasonably good, with almost no offset at 8 μm and 160 μm. Moving on to two-band estimates, and in particular the estimate based on 24 μm plus 160 μm, we see a remarkable improvement in the scatter at all luminosities. Finally, the four-band estimate follows the reference values almost perfectly. Apparently, the fact that our multiple-band relations are not sSFR dependent is accounted for by the presence of color terms. From this analysis we conclude that our relations are generally applicable even down to log L_TIR/L_⊙ = 6, though some caution in the low-luminosity regime may be warranted for single-band estimates based on bands in the intermediate wavelength range.

Fig. 11. Accuracy and precision of the recovery of the total IR luminosities for galaxies from the Local Volume Legacy (LVL) survey using our relations. Gray points represent galaxies that fall below the LTIR range of our H-ATLAS sample from which the relations were derived. Values in panels indicate the mean offset and standard deviation of the residuals for galaxies with log LTIR/L⊙ > 9. Estimates based on single bands take into account the sSFR dependence of the relations. The red line represents the 1:1 relation. Reference LTIR values come from using all IR bands and Draine et al. (2014) models.

6.3.3. Dusty local LIRGs

We evaluate the applicability of the new relations with respect to the galaxies from the Great Observatories All-Sky LIRG Survey (GOALS, Armus et al. 2009), which contains 181 LIRGs (11 < log L_TIR/L_⊙ < 12) and 21 ULIRGs (12 < log L_TIR/L_⊙ < 12.4). Since the majority of low-redshift U/LIRGs are merger driven and are much dustier than normal galaxies with similar L_TIR, we are testing potentially different IR SED shapes from any galaxy in our sample.

First, we remove 64 galaxies (32% of the sample) that are associated with an AGN in the SIMBAD database. The removal of AGN also removes essentially all ULIRGs in this sample, leaving galaxies that, at least in terms of L_TIR fall within the range covered by our H-ATLAS sample. Furthermore, in order to use sSFR-dependent relations for single-band estimates, we require stellar masses, which we take from U et al. (2012). This leaves a final sample of 39 galaxies. For our analysis we use IRAS fluxes in the 12, 25, and 100 μm bands taken from the IRAS Revised Bright Galaxy Sample (RBGS, Sanders et al. 2003) and PACS/SPIRE Herschel photometry from Chu et al. (2017). For galaxies where Herschel resolves the IRAS source into two components, we sum the component fluxes. Distances are taken from Armus et al. (2009), which also provides luminosities derived using the Sanders & Mirabel (1996) formula, which was constructed from single-temperature dust emission models and utilizes all four IRAS bands (12, 25, 60, and 100 μm). We use CIGALE with the Draine et al. (2014) models, with the parameters chosen the same way as in Sect. 3.1.1, to derive new L_TIR based on all of IRAS and Herschel photometry (9 bands in total). Our own estimates of L_TIR and those from Armus et al. (2009) have a mutual scatter of only 0.02 dex and no systematic bias. To evaluate the relations derived in this work we continue by using our L_TIR values.

Since we have not calculated our relations for IRAS bands, we use WISE W3 and Spitzer MIPS 24 μm bandpasses instead. The performance of our relations for GOALS LIRGs is presented in Fig. 12. We produce single-band estimates based on 12, 25, 70, 100, and 160 μm (100 μm is from Herschel) and using sSFR-dependent relations. Systematics exist for PACS bands, but are relatively small (up to 0.1 dex). We believe these systematics emerge because nearby LIRGs, selected from shallow IRAS survey, would be biased towards higher IR luminosities for their sSFR and dust mass compared to a more typical population of galaxies with such dust mass and sSFR. In terms of the precision, the 70 μm has the smallest scatter, 0.05 dex, but there is a slight offset. For the estimate based on the 8 μm emission, the scatter is reasonably small (0.12 dex) and there is no bias.

Fig. 12. Accuracy and precision of the recovery of the total IR luminosities for galaxies from the Great Observatories All-Sky LIRG Survey (GOALS) using our relations. Values in panels indicate the mean offset and standard deviation of the residuals. Estimates based on single bands take into account the sSFR dependence of the relations. The red line represents 1:1 relation. Reference LTIR values come from using all IR bands and Draine et al. (2014) models.

We tested various two-band estimates, and, as expected, they provide certain improvements. In particular, estimates based on MIR alone (12 and 25 μm) and on the combination of the 25 and 160 μm fluxes (both shown in Fig. 12) are both unbiased. Finally, the four-band estimate comprising of MIR IRAS bands and two two shorter PACS bands achieves a scatter of just 0.03 dex with respect to the reference luminosities.

In conclusion, the analysis in this section provides reassurance that the application of single-band relations to JWST observations of higher-redshift galaxies, even if submillimeter selected (i.e., analogous to nearby LIRGs), should yield unbiased estimates. In general, however, high-redshift luminous galaxies (assuming no AGN is present) are more likely to resemble the high-luminosity galaxies in H-ATLAS, so our relations should be even more appropriate for them.

7. Summary and conclusion

We have measured the physical properties of 2584 low-redshift star-forming galaxies by modeling their emission from the FUV to 500 μm with the spectral energy distribution modeling code CIGALE (Boquien et al. 2019), using in particular the stellar populations from Bruzual & Charlot (2003), flexible dust attenuation curves, and the dust emission models from Draine & Li (2007), Draine et al. (2014). Our sample spans a wide range of stellar masses and sSFR, and excludes AGN. The modeling facilitated a detailed study of the dependence of dust emission spectra on different physical properties (M_star, total SFR, sSFR, L_TIR, M_dust, and the oxygen abundance), from which we conclude the following:

• The overall shape of dust emission spectra varies strongly with L_TIR (or, alternatively, with SFR), but also on sSFR, confirming the findings of other studies.

• The dependence of the shape of the dust emission spectra on sSFR is independent from the L_TIR dependence, that is, sSFR is a second parameter.

• Monochromatic IR luminosities best constrain L_TIR at 90 μm rest frame, with the minimum dispersion of only ∼0.05 dex.

• The current SFR is better constrained by the total IR luminosity (0.11 dex) than by any monochromatic luminosity, though the luminosity circa 90 μm rest-frame comes close in precision (0.12 dex).

• Precise estimates of L_TIR are possible with only a single band. In particular, single-band photometry in two IR regions provides constraints on L_TIR that are better than 0.1 dex: 12−17 μm and 55−130 μm. Performance is significantly worse between 20 and 40 μm and beyond 200 μm (all wavelengths are rest-frame). Thus the prospects for JWST MIRI are excellent.

• Monochromatic IR luminosities constrain the total (obscured plus unobscured) SFR less well (minimum dispersion of 0.12 dex) than L_TIR because of the diversity of stellar populations and attenuation curves at a given luminosity. The knowledge of the population age (sSFR) reduces this uncertainty to 0.09 dex. These precisions are in either case remarkable.

• Dust mass can be constrained reasonably well (dispersion < 0.2 dex) only at rest-frame λ > 200 μm, with longer wavelengths giving progressively more accurate estimates (∼0.1 dex at 1 mm).

Informed by the above results, we constructed single-parameter IR templates parameterized on L_TIR and (separately) on SFR, as well as the functional relations that allow the determination of these two properties using between one and up to four commonly used bandpasses from Spitzer, WISE, Herschel, and, in the future, JWST. These templates and relations can be used in a large number of circumstances, in particular when there is only a limited coverage of the dust emission that precludes the fitting of theoretical dust models that contain many free parameters, such as those of Draine & Li (2007) or da Cunha et al. (2008). Furthermore, the dependence of the dust emission spectra on sSFR at fixed L_TIR motivated the development of two-parameter dust templates and relations. Two-parameter templates increase the reliability of estimation of dust parameters (including L_TIR) for galaxies with a wide range of star formation activity, from relatively quiescent local galaxies, to galaxies that, with their high sSFR, resemble typical star-forming populations at higher redshifts. Based on this we additionally conclude the following:

• Previously published dust emission templates based on samples selected from shallow FIR surveys generally have a bluer peak than ours, consistent with our sample being on average more representative of normal star forming galaxies. Our two-parameter templates do, however, become progressively bluer for higher sSFR, as expected.

• Our relations provide improvements over some literature relations at MIR wavelengths, in particular accounting for the effect of variations of the dust temperature.

• Our relations used with the individual single bands in the MIR and FIR yield fairly unbiased estimates of L_TIR (within 0.1 dex) for a diverse population of galaxies, including local LIRGs and dusty galaxies in the green valley.

• Utilizing multiple bands, in particular from different wavelength ranges, considerably reduces the scatter in estimating L_TIR. Relations involving multiple bands were tested to be unbiased even down to log L_TIR/L_⊙ = 6.

Finally, we make available software tools to generate L_TIR-dependent and L_TIR+sSFR-dependent spectra at any value of these parameters, eliminating the need for the traditional discrete templates, as well as the tools to derive L_TIR and SFR from a wide selection of bands at 0 < z < 4. We also provide explicit formulae to calculate L_TIR and the total SFR directly from JWST 21 μm (F2100W) fluxes and redshifts.

Acknowledgments

We thank the anonymous referee for helpful comments that have contributed to improving and clarifying the manuscript. We also thank Dan Smith for kindly providing us his updated galaxy templates, Steve Maddox and Steve Eales for their suggestions regarding PACS photometry, and the entire H-ATLAS team for their wonderful survey. M.B. was partially supported by FONDECYT regular grants 1170618 and 1211000, and the work of S.S. was partially supported by NASA award 80NSSC20K0440. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

Appendix A: Redshift-dependent relations for estimating L_TIR and SFR from a single IR band

In Sect. 5.2, we examined the case of galaxies in the local universe. However with increasing redshifts, the K-correction becomes a major factor and a significant source of uncertainty, in particular given the complex shape of the dust emission. At the same time, more distant galaxies tend to have fewer bands available, increasing the need to have efficient single-band estimators of L_TIR and the SFR. To fill this need, we have computed the coefficients for redshift-dependent relations, adding the redshift as a variable to Eq. 2:

$$\begin{array}{c}\hfill logp=m(b,z)\times log\lambda L_{\lambda }(b)+n(b,z),\end{array}$$(A.1)

with m(b, z) and n(b, z) the coefficients obtained from the fit. The other terms are defined similarly as in Eq. 2. We have built a grid of estimators up to z = 4 by fitting the redshifted best-fit models with this equation. We show the evolution of m(b, z) in Fig. A.1 and A.2 for the main Spitzer, WISE, Herschel, and JWST bands. Naturally, the values at z = 0 correspond to what was presented in Table 1. We observe that for the most part the curves appear to be qualitatively similar, but offset from one another with respect to the redshift, which is expected as different filters probe the same rest-frame wavelength at different redshifts and the derived estimators depend on the rest-frame wavelength. The exact filter shape only has a limited effect, as can be seen from the evolution of m(b, z) as well as for the standard deviation of the residuals. Overall we find that at all redshifts the FIR bands sampling the peak of the emission show a better performance, with always a band simultaneously with m ∼ 1 and a low scatter in the residuals. The performance of MIR bands tends to be significantly worse. The bands in the 20 μm to 25 μm range appear to work best at z ∼ 0.5, with a progressive degradation with increasing redshift as they become heavily dominated by PAH bands. All MIR bands present a sublinearity at all redshifts. The SFR estimators present a globally similar behavior.

Fig. A.1. Top row. Evolution of m(b, z) from Eq. A.1 as a function of the redshift to estimate the TIR luminosity from the main WISE, Spitzer, Herschel, and JWST bands. Bottom row. Standard deviation of the residual of the TIR luminosity estimated from Eq. A.1. Fewer bands are available at higher redshifts as the computation is only carried out for λ/(1 + z) > 6 μm, with λ the pivot wavelength of the filter.

Fig. A.2. Same as Fig. A.1, but for the SFR.

Appendix B: Construction of select previously published templates

B.1. Chary & Elbaz (2001) templates

Chary & Elbaz (2001) used the following approach to construct their templates. They started with a collection of data from local galaxies observed at 6.7, 12, 15, 60, and 850 μm with ISOCAM, IRAS and SCUBA, as well as estimates of their total and FIR luminosities from IRAS bands (12, 25, 60, and 100 μm). From this they construct linear (or piecewise linear, in log) relations between various monochromatic and integrated IR luminosities. The samples are not necessarily the same ones in different relations, and vary in number from ∼50 to a few hundred. The galaxies generally have total IR luminosities that lie in the same range (9 < log L_TIR/L_⊙ < 12) as in our study, but probably with a greater share of starburst galaxies (i.e., galaxies with higher L_TIR than typical for their stellar and dust masses). In order to get templates that are continuous in wavelength, Chary & Elbaz (2001) used Silva et al. (1998) theoretical SED (but updated with real MIR spectra) that describe four nearby galaxies representative of different IR SEDs. These “proto-templates” were constructed by Silva et al. (1998) by fitting a model that includes three grain species and three dust environments. Chary & Elbaz (2001) then split the four proto-templates into two wavelength regimes at 20 μm and interpolated them over many luminosities. The final 105 luminosity-dependent templates were chosen among the interpolated proto-templates that best fit the observed relations described above. The two wavelength regimes were then joined, and the L_TIR of each template was determined by integration. The templates span 0.1–1000 μm, that is, they also include the (highly uncertain) stellar emission. The templates provided to the community span 8.4 < log L_TIR/L_⊙ < 13.5, that is, they have been extrapolated into the ULIRG range.

B.2. Dale & Helou (2002) and Dale et al. (2014) templates

The templates of Dale & Helou (2002) represent a refinement of the templates produced by Dale et al. (2001), which in turn rely on models of Desert et al. (1990). Specifically, Dale et al. (2001) start by producing their local (as opposed to galaxy-wide) model templates by combining theoretically predicted emission from large grains, emission from very small grains and an average PAH spectrum from actual observations. This PAH spectrum includes an empirically motivated damping factor that increases as the heating intensity (U) increases. A range of global model SED were obtained by summing up local SED according to the power-law exponent α_SF, which represents the varying mix of contributions of regions with different levels of activity (i.e., heating intensities, U). Dale & Helou (2002) empirically refine the FIR portion of the Dale et al. (2001) templates (λ > 100 μm) by allowing the emissivity to vary as a function of U in a way that minimizes residuals with respect to 850 μm observations. Unlike the templates of Chary & Elbaz (2001), the templates of Dale et al. (2001) are constrained by, but not directly fit to the data. Nevertheless, they show that their models well reproduce the average SED of a sample of 69 normal galaxies with 8.2 < log L_TIR < 12.0 and spanning 6.7–100 μm. Relating the power-law exponent α_SF via FIR color to the total IR luminosity allows Dale & Helou (2002) templates to be used with single-band measurements. Such a calibration by Marcillac et al. (2006), based on IRAS Bright Galaxy Sample, associates 64 templates of Dale & Helou (2002) with 8.3 < log L_TIR/L_⊙ < 14.3. More recently, Dale et al. (2014) updated the PAH spectrum of the Dale & Helou (2002) templates and allowed the addition of an AGN-heated component.

B.3. Rieke et al. (2009) templates

Rieke et al. (2009) construct what they refer to as average templates separately for galaxies above and below log L_TIR/L_⊙ = 11, using somewhat different (but mostly empirical) methods for each group. For high-luminosity (log L_TIR/L_⊙ > 11) templates, they first produce continuous SED for 11 individual galaxies (local LIRGs and ULIRGs), which we again refer to as proto-templates. These proto-templates are constructed from a combination of synthetic SED in the optical and near-IR, empirical spectra in the MIR (6-35 μm), and modified blackbodies fit to the photometry in the FIR (> 60 μm). The proto-templates are then combined with different weights in order to match the empirical colors obtained from the relations between luminosities in different bands (8, 12, 24 and 60 μm) and L_TIR of ∼70 galaxies spanning 10 < log L_TIR/L_⊙ < 12.3. To produce intermediate luminosity (log L_TIR/L_⊙ < 11) templates, Rieke et al. (2009) rely on the combination of Dale & Helou (2002) models for the FIR and empirical IRS templates for the MIR. Dale & Helou (2002) models are implemented at > 70 μm through the correlation between α_SF and L_TIR, whereas the MIR templates (< 37 μm) are associated to L_TIR via the 12-to-25 μm color. The two were combined by matching to the 70 and 160 μm photometry and to the 25-to-60 μm color. The final set consists of 14 luminosity-dependent templates with 9.75 < log L_TIR/L_⊙ < 13.0 and covering a wavelength range from 4 μm to 30 cm.

B.4. Smith et al. (2012) templates

Binned SED (templates) of Smith et al. (2012) are constructed by median averaging individual best-fit spectra obtained with UV, optical, and FIR SED fitting with MAGPHYS. Modeling of the IR SED emission in da Cunha et al. (2008), and implemented in MAGPHYS SED fitting code, is mostly phenomenological. First, they model IR SED as a sum of the emissions from two sources: birth clouds and the diffuse ISM. Birth cloud emission consists of three components: PAH emission lines, hot MIR continuum of very small grains (peaking around 20 μm), and warm grain continuum (peaking around 50 μm). Emission of the diffuse ISM has these three components plus the cold grain continuum (peaking around 100 μm). The PAH spectrum is empirical and of fixed shape, whereas continuum components are modified black bodies of fixed temperature, except for warm birth cloud grains and cold ISM grains, for which the temperature is adjustable. Each component contributes to some degree to the L_TIR budget of birth clouds or the diffuse ISM. For the diffuse ISM the contributions of PAH, hot and warm dust are mutually fixed. In the end, the model of da Cunha et al. (2008) has six adjustable parameters: one controlling the relative contribution of the diffuse ISM to the total dust luminosity (f_μ), two describing the diffuse ISM emission (the relative contribution and temperature of cold grains) and three describing the birth clouds emission (two relative contributions plus the temperature of warm grains). The f_μ parameter is correlated with the sSFR and plays a strong role in the overall shape of the SED in the FIR since it de facto determines the relative strengths of the warm and cold components.

Smith et al. (2012) SED are based on a 250 μm-selected sample from an initial portion of H-ATLAS. Some 20% of their initial sample of ∼1000 was detected at 160 μm and 10% at 100 μm. WISE data were not available at the time. Given that da Cunha et al. (2008) models contain several flexible components in the wavelength range not covered by the data, the shape of those parts of the SED is not directly constrained. In this work, we compare our results with an updated version of Smith et al. (2012) templates, based on significantly larger number of galaxies and with greater fraction of PACS detections, but still without the constraints in the MIR.

Appendix C: Impact on the accuracy of L_TIR and SFR when using multiple IR bands

In order to quantify the effect of using an increasing number of bands, we have computed the mean standard deviation of the residuals for each band and when one, two, or three other bands are also used in addition to this reference band. For completeness, we use all possible combinations for each of the 18 reference bands when adding one (total of 153 combinations), two (816 combinations), and three (3060 combinations) bands. The results for L_TIR and the SFR are presented in Fig. C.1. As expected, the standard deviation of the residuals diminishes with the increasing number of bands at all wavelengths. The effect is especially important for the determination of L_TIR with sustained gains even with a larger number of bands. The gain is more moderate for the SFR beyond the addition of a second band, and the scatter always remains above 0.10 dex. This floor is set by the intrinsic scatter between L_TIR and the SFR that is due to the variations in the SFH (stellar populations) and of the attenuation.

Fig. C.1. Mean standard deviation of the residuals of LTIR (left) or the SFR (right) using one (blue, 18 bands), two (orange, 153 band combinations), three (green, 816 band combinations), or four bands (red, 3060 band combinations).

This general overview, however, provides limited information on the added value of a given band because it averages over all possible combinations. To investigate this, in Fig. C.2 we plot the reduction in the standard deviation of the residuals when adding one band to a combination with one, two, or three reference bands at other wavelengths. Unsurprisingly some bands have a strong impact on the reduction of the standard deviation. In the case of L_TIR, the Herschel 70 μm and 100 μm bands generate the strongest improvement, reducing the standard deviation of the residuals by 0.1 dex, on average, when used as a second band. Even as a third or fourth band, they still bring on average an improvement better than 0.05 dex. This is a stronger added value than for the MIR bands that only bring a more limited improvement that is on average always under 0.05 dex, even as a second band. With the increasing number of MIR bands their value strongly decreases. However their weaker performance should be interpreted with caution. Because of the large number of MIR filters, there are numerous combinations of only MIR bands. As mentioned earlier, bands that are close in wavelength tend to bring very similar information. For instance, if we take the Herschel 250 μm band as a reference, the standard deviation of the residuals is 0.17 dex. The addition of the JWST 18 μm band reduces it by 0.08 dex compared to a reduction of 0.12 dex when we rather include the Herschel 70 μm band. For the SFR, the improvement brought by adding more bands is limited. Even in the most favorable of cases, the improvement is less than 0.05 dex and the marginal improvement readily drops below 0.02 dex when adding a third band.

Fig. C.2. Reduction of the mean standard deviation of the residuals of LTIR (left) or the SFR (right) when adding a given band to a combination of one (blue, 153 band combinations), two (orange, 816 band combinations), or three (green, 3060 band combinations) bands at other wavelengths.

Overall for evaluating L_TIR at z ∼ 0 the best band combinations are:

• One band: Herschel 100 μm (σ = 0.0504 dex).

• Two bands: JWST 25.5 μm and Herschel 100 μm (σ = 0.0231 dex).

• Three bands: JWST 18 μm, and Herschel 70 μm and 500 μm (σ = 0.0151 dex).

• Four bands: JWST 12.8 μm and 25.5 μm, and Herschel 70 μm and 500 μm (σ = 0.0142 dex).

For the SFR the best band combinations are:

• One band: Herschel 100 μm (σ = 0.1250 dex).

• Two bands: JWST 25.5 μm and Herschel 100 μm (σ = 0.1148 dex).

• three bands: JWST 18 μm, and Herschel 70 μm and 500 μm (σ = 0.1126 dex).

• Four bands: JWST 18 μm and 25.5 μm and Herschel 70 μm and 500 μm (σ = 0.1124 dex).

This shows that the estimation of L_TIR strongly benefits from the combination of various bands that sample different regions from the spectrum, in order of decreasing priority, the peak of the emission by the warm dust around 100 μm, then the hot dust emission in the MIR, and finally the cold dust in the submillimeter range. In effect these relations embed the non linear variations of the shape of the dust emission spectrum with L_TIR. For the SFR a single band already provides results close to the best performance, probably because of the intrinsic scatter between L_TIR and the SFR as mentioned earlier. We note that the assessments in this section do not take into account different levels of quality of data from different facilities and in different bands, and assumes all data having zero error.

Appendix D: Degeneracies and information contributed by individual bands

The analysis in Appendix C suggests that when different bands bring similar information, it is common that all the weight is given to a single band, probably the one that is found to be slightly more informative. This means that the fitting process to determine m_i and n effectively eliminates some bands that do not bring additional information by assigning them m_i ≃ 0. In order to quantify this, we show in Fig. D.1 a matrix giving the relative difference between the m_i coefficients of estimators using two bands for estimating L_TIR (the figure is essentially identical for the SFR estimators). It is readily apparent that when two bands probe a similar physical component (e.g., the PAH or the warm dust in the MIR, or the Rayleigh-Jeans tail in the FIR), one of these bands dominates, which effectively reduces these estimators to single-band estimators. For instance, if we consider the JWST 7.7 μm and Spitzer 8.0 μm bands, m_7.7 < 10⁻¹² and m_8.0 = 0.9018, which corresponds exactly to the coefficient for the monochromatic estimator for the Spitzer 8.0 μm band. By itself, this provides useful knowledge on which bands to preferentially acquire for estimating L_TIR or the SFR when a broad sampling of the dust emission spectrum is not possible.

Fig. D.1. Relative difference between the mi coefficients for estimators using two bands. The lighter shades indicate that they are fairly commensurate with one another whereas the darker shades indicate that one is much larger than the other. This happens in particular when one of the two bands does not bring additional information. This is especially common for bands that are close in wavelengths while probing similar physical components, such as the PAH or the warm dust in the MIR or the Rayleigh-Jeans tail in the FIR.

All Tables

All Figures

	Fig. 1. Distribution of the selected sample in the sSFR–Mstar plane (red circles) compared to the full GSWLC-X2 catalog (hexagons, with the density indicated by the color bar).
In the text

	Fig. 2. From the top-left to the bottom-right, distribution of LTIR, Mdust, SFR, Mstar, sSFR, and the oxygen abundance. For the latter the sample is restricted to 2391 objects rather than 2584 for the other physical properties as the oxygen abundance could not be computed for all the objects.
In the text

Fig. 3. From the top-left to the bottom-right, dust emission spectra parameterized on LTIR, Mdust, SFR, Mstar, sSFR, and the oxygen abundance. The spectra are calibrated in absolute luminosity for the extensive properties, but they are normalized to LTIR for the intensive properties. The color of each spectrum follows the corresponding physical property as indicated by the color bar above each panel.

In the text

Fig. 4. Scaling coefficient between the monochromatic luminosity and a physical property versus the wavelength. The scaling coefficient α for the extensive properties LTIR (blue), Mdust (orange), SFR (green), and Mstar (red) is shown in the left panel. The scaling coefficient α′ for the intensive properties sSFR (blue) and the oxygen abundance (orange) is shown in the right panel. For the latter, λLλ has been normalized to LTIR, that is, it shows the residuals at fixed luminosity. The black horizontal lines indicate respectively the threshold between superlinearity and sublinearity (left), and between correlation and anticorrelation (right).

In the text

	Fig. 5. Standard deviation of the residuals around the relation between the monochromatic luminosity and a physical property (Eq. (1)) versus the wavelength, for the extensive properties LTIR (blue), Mdust (orange), SFR (green), and Mstar (red).
In the text

	Fig. 6. Estimation of LTIR (left) and the SFR (right) from the Spitzer 8 μm (blue), JWST 18 μm (orange), and Herschel 100 μm (green) bands. The dots represent our H-ATLAS sample and the solid lines the best-fit relation computed from Eq. (2).
In the text

	Fig. 7. Estimation of Mdust from the Herschel 250 μm (blue), 350 μm (orange), and 500 μm (green) bands. The dots represent our H-ATLAS sample and the solid lines the best-fit relation computed from Eq. (2).
In the text

	Fig. 8. Estimated total (obscured plus unobscured) SFR as a function of the flux measured in JWST MIRI 21 μm filter (F2100W) and redshift, based on the relations that take into account the evolution of sSFR. The scatter results from assuming a 1-σ scatter of 0.5 dex in sSFR of actively star-forming galaxies at each redshift.
In the text

Fig. 9. Comparison of the templates derived in this work with those from Chary & Elbaz (2001) (top left), Dale & Helou (2002), Dale et al. (2014) (top right), Rieke et al. (2009) (bottom left), and Smith et al. (2012) (bottom right). The solid lines correspond to the new templates and the dashed lines the literature templates. Each pair of lines of the same color corresponds to log LTIR/L⊙ of 9.0, 9.5, 10.0, 10.5, 11.0, and 11.5, following the bar above each panel. Spectra outside of the definition range of a given set of templates have been omitted. The Chary & Elbaz (2001) and Smith et al. (2012) templates include stellar populations, precluding any comparison at short wavelengths where such populations dominate.

In the text

	Fig. 10. Comparison of the templates parametrized on LTIR at log LTIR/L⊙ = 9.0, 9.5, 10.0, 10.5, 11.0, and 11.5 (black lines), and the templates parametrized on both LTIR and sSFR for sSFR/yr−1 = −11, −10.5, −10, −9.5, and −9. Templates corresponding to a high sSFR have a stronger emission below 90 μm.
In the text

Fig. 11. Accuracy and precision of the recovery of the total IR luminosities for galaxies from the Local Volume Legacy (LVL) survey using our relations. Gray points represent galaxies that fall below the LTIR range of our H-ATLAS sample from which the relations were derived. Values in panels indicate the mean offset and standard deviation of the residuals for galaxies with log LTIR/L⊙ > 9. Estimates based on single bands take into account the sSFR dependence of the relations. The red line represents the 1:1 relation. Reference LTIR values come from using all IR bands and Draine et al. (2014) models.

In the text

Fig. 12. Accuracy and precision of the recovery of the total IR luminosities for galaxies from the Great Observatories All-Sky LIRG Survey (GOALS) using our relations. Values in panels indicate the mean offset and standard deviation of the residuals. Estimates based on single bands take into account the sSFR dependence of the relations. The red line represents 1:1 relation. Reference LTIR values come from using all IR bands and Draine et al. (2014) models.

In the text

Fig. A.1. Top row. Evolution of m(b, z) from Eq. A.1 as a function of the redshift to estimate the TIR luminosity from the main WISE, Spitzer, Herschel, and JWST bands. Bottom row. Standard deviation of the residual of the TIR luminosity estimated from Eq. A.1. Fewer bands are available at higher redshifts as the computation is only carried out for λ/(1 + z) > 6 μm, with λ the pivot wavelength of the filter.

In the text

	Fig. A.2. Same as Fig. A.1, but for the SFR.
In the text

	Fig. C.1. Mean standard deviation of the residuals of LTIR (left) or the SFR (right) using one (blue, 18 bands), two (orange, 153 band combinations), three (green, 816 band combinations), or four bands (red, 3060 band combinations).
In the text

	Fig. C.2. Reduction of the mean standard deviation of the residuals of LTIR (left) or the SFR (right) when adding a given band to a combination of one (blue, 153 band combinations), two (orange, 816 band combinations), or three (green, 3060 band combinations) bands at other wavelengths.
In the text

Fig. D.1. Relative difference between the mi coefficients for estimators using two bands. The lighter shades indicate that they are fairly commensurate with one another whereas the darker shades indicate that one is much larger than the other. This happens in particular when one of the two bands does not bring additional information. This is especially common for bands that are close in wavelengths while probing similar physical components, such as the PAH or the warm dust in the MIR or the Rayleigh-Jeans tail in the FIR.

In the text