© The Authors 2023

1. Introduction

Galaxies are among the most important structures of the Universe, hosting the majority of its cold baryonic mass. The fundamental process of star formation (SF), which takes place within galaxies, transforms the gas into stars and therefore plays a significant role in all aspects of galaxy evolution and their current state. For example, the current star formation rate (SFR) of galaxies is strongly correlated with their stellar mass (M_⋆; e.g., Elbaz et al. 2007, also known as the main sequence of galaxies), with their galactic winds (e.g., Heckman et al. 1990; Lehnert & Heckman 1996), the number of supernovae (e.g., Mac Low et al. 2005), and their X-ray luminosity originating from high-mass X-ray binaries (e.g., Mineo et al. 2012; Kouroumpatzakis et al. 2020).

Similarly, the M_⋆ of galaxies is strongly correlated with their kinematic and morphological parameters (e.g., Bundy et al. 2005; Genzel et al. 2006), while it is also related to the nature of their stellar content and the properties of their interstellar medium (ISM), for example through the mass–metallicity relation (e.g., Tremonti et al. 2004). Moreover, there is a strong correlation between the M_⋆ of the bulges of galaxies with the mass of their central supermassive black hole (e.g., Kormendy & Ho 2013). Therefore, being able to accurately estimate the SFR and M_⋆ is crucial for most studies involving galaxies, and studies of the evolution of the Universe.

Galaxies come in many forms and shapes, and in different states with respect to their current star-forming activity or ISM properties. This is the main reason for discrepancies in measuring the SFR because the various SFR tracers are based on different emission mechanisms related to SF. These tracers include: the ultraviolet (UV) light coming directly from the photospheres of massive stars, emission lines from nebulae that have been ionized by star-forming activity (e.g., Hα, [O III]), infrared (IR) emission from dust that was heated by the young stars’ UV and optical emission, radio emission from relativistic electrons related to supernovae activity, and several others (for a review see Kennicutt & Evans 2012).

Depending on the availability of observations, different methods are used to estimate the SFR. However, it has been shown that there can be significant discrepancies between these methods, sometimes up to an order of magnitude (e.g., Kennicutt & Evans 2012). In particular, although monochromatic IR-based SFR indicators are good for dusty and high-SFR galaxies, they can lead to a significant underestimation of the true SFR in low-mass or low-metallicity dust-deficient galaxies (e.g., Calzetti et al. 2007; Kouroumpatzakis et al. 2021). Moreover, IR SFR tracers can also be biased by the age of the stellar population (SP; e.g., Cortese et al. 2008; Leroy et al. 2012; Boquien et al. 2014; Cerviño et al. 2016; Nersesian et al. 2019; Kouroumpatzakis et al. 2020) or stochastic heating of the dust (e.g., Camps et al. 2015; Lianou et al. 2019). On the other hand, the use of UV or Hα emission as SFR tracers requires precise measurements of extinction in order to provide reliable estimations.

A possible solution is to use one of the hybrid SFR indicators (e.g., 24μm + Hα, FIR + UV; e.g., Kennicutt et al. 2009; Hao et al. 2011). They combine emission that traces SF in a more direct way but it can be affected by absorption, along with IR emission which is not. Thus, these hybrid tracers account for the energy lost from the UV/optical light that has suffered absorption by measuring the corresponding IR light. However, it is relatively harder to have observations in both wavelength regimes due to the lack of deep UV or spectral observations. Spectral energy distribution (SED) fitting can provide robust estimations for both the M_⋆ and SFR, but they require extended photometry in multiple bands, while the lack of UV or IR photometry can still lead to large biases (e.g., Lanz et al. 2013).

As a result of these limitations, most catalogs providing SFRs and stellar masses for local Universe galaxies were limited in the Sloan Digital Sky Survey (SDSS) footprint where the sources were mainly characterized through SED fitting combining the optical with WISE (Wright et al. 2010) and/or GALEX (Martin et al. 2005) photometry (e.g., Chang et al. 2015; Salim et al. 2016, 2018). Other surveys not limited to the SDSS footprint (e.g., Sheth et al. 2010; Ashby et al. 2011; Bianchi et al. 2018; Leroy et al. 2019; Nersesian et al. 2019) combined IR and UV GALEX observations. However, all of them were limited to relatively nearby galaxies and small samples. Catalogs targeting the whole sky were based solely on IR photometry (e.g., Jarrett et al. 2013; Kovlakas et al. 2021).

A significant problem in measuring SFRs through monochromatic IR emission is that they do not account for the different degrees of extinction between galaxies, especially in low-metallicity, dust-poor galaxies. As a result, an IR SFR indicator would underestimate the SFR for a galaxy with less dust, because in this case, a larger amount of UV photons, produced by the star-forming activity, would have escaped with respect to a higher-extinction galaxy. Moreover, the fact that dust can be excited by other means besides the reprocessing of UV radiation from young stars (e.g., stochastic heating, heating by older SPs, post-AGB stars, etc. Sauvage & Thuan 1994; Draine & Li 2001; Nersesian et al. 2019; Zhang & Ho 2023) adds another complication to the IR SFR indicators. In fact, several studies have shown that in quiescent galaxies the mid/far-IR emission is higher with respect to their SFR (e.g., Davis et al. 2014; Simonian & Martini 2017; Leroy et al. 2019). Therefore, for these cases, an IR SFR indicator would overestimate their SFR.

Similarly, the M_⋆ measurements can be obtained by means of SED fitting, which requires photometry in an extended wavelength range. However, most commonly they are based on near-IR (NIR) photometry (e.g., 2MASS K_S, WISE band-1 or band-2) which traces the thermal/continuum emission of low-mass SPs. Such observations provided M_⋆ measurements for a large amount of relatively nearby galaxies through the all-sky IR surveys (e.g., IRAS, AKARI, 2MASS, WISE). However, it has been shown that the relation between NIR emission and the M_⋆ depends strongly on the age of the SPs. In the past, optical (e.g., Bell et al. 2003; Zhu et al. 2010; Wen et al. 2013) or IR colors (e.g., Jarrett et al. 2013, 2023) were used as tracers of the SP age and thus, as a way to correct for their contribution in the NIR emission. However, these corrections depend on the color of choice and its ability to trace the SP age or the biases induced by interstellar reddening. The existing calibrations are based on the intrinsic optical colors whereas in their application the observed colors are used. Therefore, results from various calibrations show discrepancies and large scatter in the estimations of the M_⋆ (e.g., Bonfini et al. 2021). Of course, not accounting for the effect of extinction further exacerbates these discrepancies.

This work investigates the biases introduced by the SP age and extinction in the widely used calibrations of photometric IR luminosity (L_IR)-to-SFR, and mass-to-light ratio, and proposes new calibrations that try to mitigate these effects. The analysis is based on creating a large grid of galaxy models with different ISM (e.g., metallicity, ionization parameter U) and dust properties, for a large range of star-forming conditions. The obtained SEDs and monochromatic luminosities of these galaxy models are then compared with the model SFRs and stellar masses in order to derive the corresponding calibrations.

In Sect. 2 we discuss the motivation of this work and give an example that reveals the physical reasons behind the discussed biases. The used model suite is presented in Sect. 3. The details about the analysis and the SFR and M_⋆ calibrations are provided separately in Sect. 4. In Sect. 5 the proposed calibrations are compared with previous works and observations, and in Sect. 6 we present the conclusions. This work considers only rest-frame luminosities. Throughout this work, unless stated otherwise, given values correspond to the modes of the distributions and the uncertainties to 68% percentiles.

2. Biases in star formation rate and stellar mass estimations from infrared emission: An example

In order to explore the possible biases in measuring the SFR, and M_⋆ from IR photometry we plot in Fig. 1 the young stellar, old stellar, and dust components from the SEDs of two model galaxies with low (sSFR = 10^−11.75 M_⊙ yr⁻¹/M_⊙) and high (sSFR = 10^−8.5 M_⊙ yr⁻¹/M_⊙) specific SFR (sSFR ≡ SFR/M_⋆) using the code CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019). Plotted in Fig. 1 are also the transmissions of the WISE bands 1 and 3, and JWST NIR-F200W, MIRI-F2100W (Gardner et al. 2006; Gordon et al. 2022) bandpasses showing the parts of the spectrum used to extract information about the M_⋆ and the star-forming activity.

Fig. 1. Intensity as a function of wavelength for the total attenuated intensity (gray), and the components of the young stellar (blue), old stellar (orange), and dust emission (red) for two model galaxies: (a) in the top panel is a galaxy with low star-forming activity (SFR = 0.01 M⊙ yr−1, sSFR = 10−11.75 M⊙ yr−1/M⊙) and low-extinction (E(B − V) < 0.001), and (b) in the bottom panel is a dusty (E(B − V) = 0.55) galaxy with high star-forming activity (SFR = 68 M⊙ yr−1, sSFR = 10−8.5 M⊙ yr−1/M⊙). The dashed-dotted lines with brown and purple colors represent the WISE bands 1 and 3 transmission bandpasses respectively. The dotted fuchsia and red lines represent the JWST NIR-F200W, and MIRI-F2100W bandpasses respectively.

These two extreme SEDs demonstrate that the emission in these IR bands corresponds to different physical processes when galaxies are in different star formation states. In the NIR (e.g., WISE band-1, JWST NIR-F200W), for the low SF activity galaxy, the stellar emission is orders of magnitude higher compared to the dust emission. Calibrations of mass-to-light ratio are based on emission in this wavelength range by tracing the thermal emission of the SPs. However, we see that in a highly star-forming galaxy, the dust component dominates the emission even in the NIR wavelengths, and therefore if it is not taken into account it will significantly bias the stellar-mass estimations.

Similarly for the high SF galaxy, in the mid-IR (MIR; e.g., WISE band-3, JWST MIRI-F2100W) wavelengths, emission from the dust component dominates the SED and it is orders of magnitude higher compared to the stellar continuum emission. However, we see that for a low-dust/low-sSFR galaxy, the stellar thermal emission can dominate the emission in that bandwidth, enough to significantly bias the SFR estimation. This example shows that both the SP age and extinction have to be taken into account to properly estimate the SFR and M_⋆ through the IR emission across the full range of galactic ISM environments.

It should be noted that Fig. 1 shows that regardless of the level of the star-forming activity the shape of the spectrum attributed to the dust component is not changing. This is because CIGALE is not performing radiative transfer that would allow a more detailed description of the SED based on the amount and chemical composition of the dust, the geometric effects, and the local thermal equilibrium of the dust accounting for the SPs in the region of the dust clouds. However, radiative-transfer analysis is significantly more computationally intensive and could not be performed for this analysis which requires a large number of galaxy models needed to cover the variety of star formation histories (SFHs) and ISM conditions. Therefore, in the following analysis dust temperature variations are not being considered, but we expect them to play a minor role in comparison to the total luminosity due to the fact they affect mainly longer wavelengths than 24 μm (e.g., Nersesian et al. 2021), and the relative scatter in the L_IR–SFR relation being considerably less compared to that caused by extinction.

3. Modeling the SED of galaxies

3.1. Initial library of model-galaxy SEDs

In order to examine the correlation between the IR emission as a tracer of the SFR or the M_⋆, we created a large library of mock-galaxy SEDs that cover a wide range of SFHs, stellar masses, and ISM conditions. Modeling galaxies with a SED-fitting code allows us to know a priori their SFR and M_⋆, which are an outcome of the given SFHs, while the corresponding band-luminosities are given by the convolution of the filter transmissions and the resulting SEDs. For the generation of the model-galaxy SEDs, we adopted the commonly used SED fitting code CIGALE.

CIGALE is based on the principle of energy balance, where the energy corresponding to the UV/optical light that has been absorbed by the dust is re-emitted in the mid and far-IR parts of the spectrum. This makes CIGALE the ideal tool to examine the dependence of IR SFR and M_⋆ tracers on the amount of dust within a galaxy and the corresponding extinction. The generated grid of galaxies covered the extinction range for E(B − V) between 0.001 and 1, with steps of E(B − V) = 0.05 magnitudes. This E(B − V) refers to the nebular, not the stellar-continuum extinction.

Moreover, a wide range of different SFHs, covering continuously the distribution of star-forming activity between passive and highly star-forming galaxies, is required in order to investigate the dependence of the IR tracers of SFR and M_⋆ on the SP age. We adopted the sfhdelayed module, which models the coexistence of an old SP, along with a delayed and recent star formation burst that generates the young SPs. By giving a large range to the modulation parameters (Table A.1) the initial library covers a wider range of SFHs and star-forming conditions. Throughout this analysis, the SFR coming from modeling with the CIGALE code refers to the sfh.sfr10Myrs that corresponds to the average SFR over 10 Myrs. In addition, the generated grid of galaxies covered a wide range of metallicities from extremely low (Z = 0.0001) to high (Z = 0.05). The configuration of the SED creation modules is given in Table A.1. This initial setup generated 1 854 720 SEDs of galaxies. This library was extended by renormalizing the stellar masses, SFRs, and corresponding luminosities creating a grid of model galaxies with a large range of SFRs, M_⋆, and ISM conditions.

Figure 2 shows the relation between the luminosity of WISE band-3 (L_W3) per unit M_⋆ and the sSFR. In this work, in addition to the generally used E(B − V) extinction metric that can be obtained through various ways, we also utilize the ratio of WISE band-4 to the optical u, or g band luminosities as infrared excess tracers where IRXu = log L_W4/L_u and IRXg = log L_W4/L_g. The latter are chosen because of the large coverage for galaxies in these bands through wide-area surveys (e.g., PANSTARSS; Chambers et al. 2016). It should be noted that this IRX is not the same as the one given by CIGALE which refers to the ratio between the dust, and GALEX far-UV (FUV) luminosities.

Fig. 2. WISE band-3 luminosity (LW3) per M⋆ unit, as a function of sSFR for the generated library of galaxy SEDs. All panels show the same sources but in the top left they are color-coded based on their E(B − V), in the top right based on their metallicity (Z), in the bottom left on their infrared excess (IRXu = log LW4/Lu), and in the bottom right based on their infrared excess using the g band instead of u (IRXg = log LW4/Lg).

The generated galaxies show a continuous coverage of sSFRs ranging from extremely low to highly star-forming galaxies. Moreover, Fig. 2 shows that the relation between L_W3/M_⋆ and sSFR: (a) is not continuously linear: the linearity breaks in low sSFRs where beyond a point the L_W3 remains constant while the sSFR drops. Previous works based on observational data (e.g., Salim et al. 2016) showed that this break occurs near sSFR ≃10⁻¹¹ M_⊙ yr⁻¹/M_⊙ which is similar to what this analysis shows; (b) there is a strong dependence on extinction which can result in offsets of the sSFR–L_W3/M_⋆ correlation of up to 2 orders of magnitude; (c) metallicity is an additional reason for scatter. The fact that in low-sSFR galaxies the L_W3 remains almost constant although the sSFR continues to drop is a representation of the contribution of the old SPs in the MIR emission (see also Sect. 2). Omitting the common denominator M_⋆ in the sSFR–L_W3/M_⋆ comparison shows that the above conclusions also stand for the relation between SFR and L_W3.

By binning the set of simulated photometries in different metallicity bins and calculating the scatter with respect to the best SFR calibration in each bin we can disentangle the role of metallicity and extinction in the observed scatter. We note that this metallicity refers to the metallicity of the stars and because in our simulation the SP and the ISM components are independent there is no intrinsic correlation between the gas-phase metallicity (and therefore the extinction) and the stellar metallicity used as a parameter in our analysis. The scatter induced by metallicity in the relation between log SFR and log L_W3 is maximized for extremely low SFRs (SFR ≃ 10⁻⁶ M_⊙ yr⁻¹) and low-extinction (E(B − V)≃0.1) galaxies at about 0.3 dex. Therefore, the metallicity-induced scatter is by far smaller compared to that related to extinction. In addition, because measuring metallicity requires spectroscopic observations which are more costly, we do not account for the metallicity differences in the following analysis. However, we include galaxies that fully cover subsolar to super-solar metallicities as given by CIGALE (0.0001 ≤ Z ≤ 0.05). Thus, the following relations and the evaluation of their scatter encompass the dispersion caused by metallicity differences.

3.2. Fitting sample

The final master sample of model-galaxy SEDs used in the following analysis is part of the initial library, and it was selected to fully and evenly cover the plethora of star-forming conditions between passive and starburst, or dwarf and massive galaxies. It is a result of sampling randomly the initial library in order to produce uniform M_⋆ and SFR distributions whose ratio, the sSFR, is a distribution that can be described as a normal distribution with mean ⟨log sSFR/(M_⊙ yr⁻¹/M_⊙)⟩ = −10.02 and σ = 1.45. This sample of model-galaxy SEDs covers the range 10^6.5 < M_⋆/(M_⊙) < 10¹², 10⁻⁶ < SFR/(M_⊙ yr⁻¹) < 10³, and 10⁻¹³ < sSFR/(M_⊙ yr⁻¹/M_⊙) < 10^−7.5. Thus, this sample would not bias the analysis toward any particular group of galaxies. The SFR-M_⋆ plane for this master sample is shown in Fig. 3. The best linear fit is:

$$\begin{array}{c}\hfill log\frac{\mathrm{SFR}}{(M_{\odot }{\mathrm{yr}}^{-1})}=-8.93(\pm 0.02)+0.89(\pm 0.01)log\frac{M_{\star }}{(M_{\odot })}.\end{array}$$(1)

Fig. 3. SFR, M⋆ plane for the final sample of CIGALE model galaxies used for the rest of the analysis and fitting. The color code indicates their specific SFR [sSFR ≡ SFR/M⋆ (M⊙ yr−1/M⊙)]. The best linear fit is represented with a red line. The green dashed, yellow dashed-dotted, magenta dotted, and blue dashed-dotted lines represent the main sequence for diverse samples of local Universe galaxies from Elbaz et al. (2007), Speagle et al. (2014), Popesso et al. (2019), and Kouroumpatzakis et al. (2021) respectively.

A goal of this work is to provide an extinction-dependent calibration that can correct for the variance introduced by extinction in the relation between L_IR and SFR. However, because estimations of extinction are not always available, this work also provides calibrations that do not necessarily require the input of extinction. From the master sample, two subsets were created in order to account for the differences between the extinction-independent and extinction-dependent calibrations. For the extinction-dependent relations the sample followed a uniform E(B − V) distribution in order for the MCMC fits to be evenly affected by galaxies with low and high extinction. Thus, the extinction-dependent calibration can compensate and correct for the dependence of the L_IR–SFR relation on extinction (Fig. 2) using the E(B − V) or the IRX as an input. This E(B − V)-uniform sample used for the extinction-dependent calibration included 397 256 galaxy models.

The extinction-independent calibration does not have as an input an extinction indicator to correct for the variance it introduces in the L_IR–SFR relation. But still, this calibration is being applied to estimate the SFR of galaxies whatever their extinction may be. Thus, the distribution of extinction of the sample used for the fitting of the calibration plays a significant role in determining where the normalization between the L_IR and the SFR will be. The uniform E(B − V) sample has on average more galaxies with higher extinction compared to galaxies found in nature (Fig. 4). Therefore, if it was used for the fitting, it would have shifted the calibration toward dustier galaxies for whom the SFR corresponds to higher L_IR compared to average-extinction galaxies. This would have mistakenly led to lower normalization, and thus, lower SFRs for the same L_IR that would be inappropriate for the majority of the galaxies.

Fig. 4. Distribution of E(B − V) for the uniform (green histogram) and the SDSS-matched samples (red histogram) of model galaxies. With blue color is the nebular E(B − V) distribution based on the Balmer decrement for sources in the SDSS spectroscopic MPA-JHU catalog considering only star-forming galaxies with E(B − V)/E(B − V)err > 3. The ordinate corresponds to density which is equal to the ratio of the count of each bin to the total number of counts multiplied by the bin width with the integral of the area under the histogram being equal to one.

Therefore, for the extinction-independent relations, the distribution of the nebular E(B − V) of the model galaxies was matched to follow the distribution of E(B − V) of the SDSS spectroscopic MPA-JHU catalog (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004). The E(B − V) for the SDSS spectroscopic MPA-JHU catalog was calculated using the flux ratio of the Hα and Hβ emission lines based on the Balmer decrement adopting the conversion of Domínguez et al. (2013):

$$\begin{array}{c}\hfill E(B-V)=1.97\log \left[\frac{(f_{\mathrm{H}\alpha }/f_{\mathrm{H}\beta })}{2.86}\right],\end{array}$$(2)

using the reddening law of Calzetti et al. (2000). Sources with E(B − V)/E(B − V)_err < 3 were omitted in order to keep only galaxies with reliable extinction estimations. The distribution of the E(B − V) color excess for the models and the SDSS galaxies (used as a reference) are shown in Fig. 4. This SDSS-matched sample included 123 908 galaxy models.

Finally, for each of these parent samples, we created by sampling randomly two separate subsets for the fitting and testing of the best-fit models. The fitting and testing samples included 70% and 30% of the parent samples respectively ensuring that there is no overlap between them.

4. Analysis and results

4.1. Fitting methods

In all cases in the following analyses, the calibration models were fitted with the Markov chain Monte Carlo (MCMC) technique using the Python emcee package (Foreman-Mackey et al. 2013). The MCMC fitting used 64 walkers over 10 000 iterations, while the first 500 were for the burn-in phase. The initial model parameters were selected by running a maximum likelihood fit on the parameter space based on the scipy minimization algorithm (Virtanen et al. 2020). A uniform prior was adopted covering a wide range for each parameter. The models and best-fit results are presented separately for the L_IR-to-SFR and the mass-to-light ratios in the following sections.

Additionally, in order to independently test the calibrations given by the MCMC fitting method, we performed machine-learning analysis on the same datasets. We adopted the widely-used supervised machine learning algorithm Random Forest (RF) in a regression mode (for a detailed review see Louppe 2014). The building blocks of RF are the decision trees. Each decision tree is a nonparametric model which can be trained to learn the relation between the target variable (e.g. SFR, M_⋆) and the feature values (e.g observables as W1, W3, E(B − V), etc.) by using a set of continuous nodes in a tree-like structure. During the training, the algorithm searches for the feature and its corresponding value that leads to the best condition of separation at each tree node. This process is repeated recursively until a decision tree reaches its final nodes, namely very well-separated and homogenized subsets of the initial training data set. We adopted the Gini impurity (Baron 2019) for the calculation of the best separation, which is also the default criterion of the sklearn RF routine.

One of the advantages of the RF regression algorithm is that it can handle problems where the parameter space is highly complex and nonlinear, as in our case. This makes the RF method the best independent test of the MCMC calibration results which is the standard fitting method followed in this work.

We trained exactly the same calibration datasets that were fitted by MCMC, for consistency between the comparisons of the two methods. As with the MCMC method, in all RF models, we considered as training dataset the 70% of the initial sample and as the test dataset the 30%. We also tuned the basic hyper-parameters of the RF algorithm (such as the number of trees in the forest, the maximum depth of the tree, etc.) in order to avoid overfitting or underfitting. We used the implementation of the RF classifier sklearn.ensemble.RandomForestClassifier() provided by the scikit-learn¹ version 1.1.2 (Pedregosa et al. 2011) package for Python 3.

4.2. Star formation rate calibrations

In order to account for the contribution of old SPs in the MIR, we fit the relation between the MIR luminosity (L_MIR) and the SFR including a tracer of the M_⋆. The NIR (e.g., K_S, WISE 1 bands) traces part of the stellar continuum emission of low-mass stars and it has been traditionally used to measure the M_⋆ in galaxies (e.g., Cole et al. 2001; Bell et al. 2003). In this work, the L_MIR and SFR are calibrated including an observable component that traces the M_⋆. This allows to partially disentangle the contribution of old SPs in the MIR emission which can be significant, especially in low-SFR galaxies (see also Sect. 2, Fig. 1). The form of this relation is described in the following Eq. (3):

$$\begin{array}{c}\hfill \log L_{\mathrm{MIR}}=\log (\alpha M_{\star }+\beta \mathrm{SFR})\end{array}$$(3)

where we assume that the observed luminosity arises from two components: a young component that scales linearly with SFR and an older component that scales linearly with stellar mass. The older component could be either the tail of the continuum emission of the stars or associated with stochastically heated dust.

As discussed in Sect. 3 the goal of this work is to define a SFR indicator based on NIR and MIR photometry which (a) is applicable to dust-rich and dust-poor galaxies, and (b) galaxies with intense and low-level star-forming activity (where the contribution of an older stellar component in the MIR emission might be not negligible). Based on the photometric bands we adopt, we have that M_⋆ = kL_NIR. We include the effect of dust by parametrizing the parameter β in Eq. (3) as a logistic function: β = δ′/(1 + ϵ′ 10^−ζ′e) where δ′, ϵ′, and ζ′ are fitted parameters and e is an extinction metric (Eq. (4)). We adopt as an extinction metric the color excess E(B − V) or the IRX index (as defined in Sect. 3). However, because stochastically heated dust may also have a contribution in the NIR bands we also include an extinction term in the parameter α. We find that this is best parametrized with a second-order polynomial of the extinction metric: α = α′+β′e + γ′e².

In order to account for cases where there are no available extinction measurements we also fit the model data with an extinction-independent parametrization, after dropping all the extinction-sensitive terms. This extinction-independent analysis is performed on the SED models that include the effect of extinction but follow the E(B − V) distribution of SDSS galaxies (Sect. 3.2).

For the WISE photometric system, as a SFR tracer, we adopted the WISE band-3 (12.1 μm), and as a M_⋆ tracer the WISE band-1 (3.4 μm). The 12.1 μm band covers a polycyclic aromatic hydrocarbons (PAHs) band which has been used as a SFR indicator (e.g., Jarrett et al. 2013), while the 3.4 μm band is mostly dominated by stellar continuum emission with a small contribution by a PAH component (e.g., Fig. 1). For the JWST photometric system, as a SFR tracer, we adopted the MIRI-F2100W (21 μm) which probes thermal emission from dust, and as a M_⋆ tracer the NIR-F200W (2 μm) which is very similar to the K_S band. Both JWST filters have the best effective responses in the wavelength area of interest, and were designed for general purposes thus, are not significantly dominated by emission from PAHs which could lead to additional biases.

Solving Eq. (3) for the SFR term, the calibration of SFR as a function of the measured photometric quantities is:

$$\begin{array}{cc}\hfill \frac{\mathrm{SFR}}{(M_{\odot }{\mathrm{yr}}^{-1})}& =\frac{L_{\mathrm{MIR}}-\alpha L_{\mathrm{NIR}}}{{10}^{21}\beta }\hfill \\ \hfill L_{\mathrm{MIR}}& =\frac{L_{\mathrm{W}3}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})},\mathrm{or}{L}_{\mathrm{MIR}}=\frac{L_{F2100W}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})}\hfill \\ \hfill L_{\mathrm{NIR}}& =\frac{L_{\mathrm{W}1}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})},\mathrm{or}{L}_{\mathrm{NIR}}=\frac{L_{F200W}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})}\hfill \\ \hfill \alpha & ={\alpha }^{\prime }+{\beta }^{\prime }e+{\beta }^{\prime }{e}^2\hfill \\ \hfill \beta & =\frac{{\delta }^{\prime }}{1+{ϵ}^{\prime }{10}^{-{\zeta }^{\prime }e}}\hfill \\ \hfill e& =E(B-V),e=\mathrm{IRXu},e=\mathrm{IRXg}.\hfill \end{array}$$(4)

The factor 10²¹ was introduced in order to homogenize the values of SFR (M_⊙ yr⁻¹) and the IR luminosities (Watt Hz⁻¹). Following the method outlined in Sect. 4.1 we fit the above relation using the emcee implementation of the MCMC fitting scheme. Because the photometry is based on SED models given the different SFH scenarios, they come without uncertainties. Therefore the fitted likelihood functions do not include terms accounting for uncertainties in photometry or in the SFR. The likelihood function for the extinction-independent model is:

$$\begin{array}{c}\hfill p(y|x,z,\alpha ,\beta )=-\frac{1}{2}{\displaystyle \sum _n}{[y_n-\frac{x_n-\alpha z_n}{{10}^{21}\beta }]}^2,\end{array}$$

and for the extinction–dependent model is:

$$\begin{array}{cc}& p(y|x,z,e,{\alpha }^{\prime },{\beta }^{\prime },{\beta }^{\prime },{\delta }^{\prime },{ϵ}^{\prime },{\zeta }^{\prime })\hfill \\ \hfill & =-\frac{1}{2}{\displaystyle \sum _n}{[y_n-\frac{x_n-({\alpha }^{\prime }+{\beta }^{\prime }e+{\beta }^{\prime }{e}^2)z_n}{{10}^{21}(\frac{{\delta }^{\prime }}{1+{ϵ}^{\prime }{10}^{-{\zeta }^{\prime }e}})}]}^2,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \mathrm{where}:y& =\mathrm{SFR}/(M_{\odot }{\mathrm{yr}}^{-1})\hfill \\ \hfill x& =\frac{L_{\mathrm{W}3}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})},\mathrm{or}x=\frac{L_{\mathrm{F}2100\mathrm{W}}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})}\hfill \\ \hfill z& =\frac{L_{\mathrm{W}1}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})},\mathrm{or}z=\frac{L_{\mathrm{F}200\mathrm{W}}}{(\mathrm{Watt}{\mathrm{Hz}}^{-1})}\hfill \\ \hfill e& =E(B-V),e=\mathrm{IRXu},e=\mathrm{IRXg}.\hfill \end{array}$$

The best-fit results for Eq. (4) are given in Table 1.

4.3. Comparisons with the star formation rates of the model galaxies

In order to quantify the scatter in the derived SFR calibrations and any biases with respect to the input values or SFRs measured with other methods, we calculate the difference between the results of Eq. (4) applied to the model SED photometry, and their true SFRs which are an outcome of the assumed SFHs. Table 2 summarizes the comparisons with the true SFRs for low and high SFRs, and for the uniform, and SDSS-matched extinction samples. Comparisons with SFRs calculated using the calibrations of Jarrett et al. (2013), Chang et al. (2015), and Cluver et al. (2017) are also given as reference. It should be noted that the 68% percentiles provided in Table 2 are dominated by the intrinsic scatter of the data due to the various dependencies that affect the relation between the SFR and the IR luminosities (e.g., extinction distribution of the applied sample, metallicity, ionization parameter, dust content etc.). On the other hand, 68% uncertainties given in Table 1 refer only to the statistical uncertainties of the MCMC fitting procedure and the corresponding parametrization.

Figure 5 shows comparisons between the true SFR based on the assumed SFHs (SFR_true; referring to the average SFR over 10 Myr from CIGALE), and the calculated SFR based on the model photometry and the application of Eq. (4) (with the best-fit results; Table 1), and the RF models for different extinction values. As we see in Fig. 5 and Table 2 the extinction-dependent relation of Eq. (4) results in excellent agreement with the true SFR of the galaxies regardless their extinction over the full extinction range. Our analysis using the extinction-dependent parametrization provides a remarkable improvement in the reliability of the SFR measurements in low extinction galaxies in comparison to the extinction-independent parametrization or other metrics that do not account for extinction. In the case of higher extinction galaxies, the agreement is also excellent (Fig. 5, lower-right panel). In very low SFRs (≲10⁻¹ M_⊙ yr⁻¹) we see worse performance with a larger scatter.

Fig. 5. Comparison between the true SFR as given by the CIGALE output, and the best-fit models of Eq. (4) (MCMC; Table 1), and the random forest machine learning method (RF) using the WISE bands 1 and 3. The y-axis represents the logarithm of the ratio between the SFR of each model and the true SFR, thus the black dashed line at zero represents equality. Blue, red, pink, and purple colors show the extinction-independent, and extinction-dependent (E(B − V), IRXu, IRXg) relations of Eq. (4) respectively. With light-blue and orange colors are the extinction-independent, and extinction-dependent models based on the RF method respectively. The linear log LW3–log SFR calibration by Cluver et al. (2017) is also displayed in green color. The lines represent the distribution modes and the shaded areas show the 68% percentiles of the log (SFRmodel/SFRtrue) as a function of the log SFRtrue. Each panel shows the comparison of galaxy models with one of the discrete E(B − V) values, shown in the bottom left corner of each panel. This comparison is performed with a uniform-E(B − V) sample, therefore, there is an equal number of sources for each extinction value. At the bottom of each panel are the modes and 68% percentiles of the log (SFRmodel/SFRTrue) distribution for all galaxy models with the specific E(B − V), shown with corresponding colors to the relations used to derive the SFR. These comparisons are also summarized in Table 2.

The extinction-independent relation of Eq. (4) results in an overall good agreement, that is excellent for galaxies with average, and high extinction. In very low SFRs the calibrations tend to overestimate the SFRs which, however, cannot be reliably measured also due to SP sampling effects. This behavior affects all SFR calibrations due to stochasticity driven by relatively small numbers of massive stars in the probed SPs (e.g., Kennicutt & Evans 2012). The relations based on the JWST photometric system are not shown in Fig. 5 in order to avoid congestion but they yield similar results as those based on WISE (Table 2).

The RF method agrees well with the MCMC fitting method with respect to its behavior regarding extinction. However, for the extinction-dependent calibrations, this method shows slightly improved behavior for model galaxies with extremely low extinction and slightly worse results for model galaxies with extremely high extinction. Conversely, for extinction-independent calibrations, the RF performs better in high-extinction galaxies and worse in extremely low-extinction galaxies. This is probably because of the large flexibility of the RF method with respect to MCMC because it does not fit an explicit functional form to the data. Due to the behavior of the RF models, in the following part of this work we refer to and provide only analytical calibrations which are a result of the MCMC fitting. However, the fact that the RF method yields similar results to the MCMC indicates that the SFRs from the proposed calibrations of Eq. (4) do not strongly depend on the fitting method.

Equation (4) relations offer excellent agreement while minimizing the scatter. It is worth noting that the results between the true and inferred SFRs for the extinction dependent and independent relations are almost the same for the SDSS-matched-E(B − V) sample. This is mainly because the SDSS-matched-E(B − V) distribution mainly consists of average extinction galaxies (E(B  −  V)  ≃  0.3), while it completely lacks galaxies with E(B − V) < 0.05 (Fig. 4). Thus, the use of the extinction-independent relation will not significantly bias the SFR estimations for medium extinction galaxies, or samples of galaxies similar to that of SDSS (e.g., local Universe). However, for a sample including galaxies with extreme (low or high) values of extinction, the use of the extinction-dependent relations yields better results with less scatter. As also shown in Fig. 5 this can be particularly important for dwarf, and dust-free galaxies.

Table 2 shows the modes and 68% C.I.s of the distributions of the log(SFR_model/SFR_true) for the calibrations of Eq. (4), and that of Jarrett et al. (2013), Chang et al. (2015), and Cluver et al. (2017) for both the SDSS-matched-E(B  −  V), and uniform-E(B  −  V) samples considered here, separated to low, and high-SFR galaxies. The relation of Cluver et al. (2017) shows good agreement for average/high-extinction and high-SFR galaxies but the discrepancy increases as the extinction decreases. These comparisons show that the Cluver et al. (2017) calibration can lead to an underestimation of the SFR for dust-free galaxies by more than an order of magnitude. A probable explanation for these discrepancies is the selection of the sample. The calibration of Cluver et al. (2017) used the SINGS/KINGFISH sample of galaxies supplemented with some dwarf and some ultra-luminous IR galaxies (ULIRGs). On average the agreement is better and the scatter is reduced using the calibrations of Eq. (4) as shown from the modes and 68% percentiles of the distributions of SFR_model/SFR_true (in the bottom of the individual panels; Table 2).

Moreover, regardless of the extinction, the commonly used linear L_W3–SFR calibrations fail to recover the true SFR for passive galaxies. Figure 5 shows that the Cluver et al. (2017) calibration overestimates the SFRs in the low-SFR regime regardless of their extinction. This is because in low-SFR galaxies (e.g., ≲0.01 M_⊙ yr⁻¹) the increasing contribution of the emission from old SPs becomes dominant in the MIR photometric bands like WISE band-3 (e.g., Figs. 1, and 2). On the other hand, the relation of Cluver et al. (2017) shows better agreement for the high-SFR sample that follows the SDSS extinction distribution although it still has a tendency to overestimate the SFRs as shown from the upper 68% percentile of the distribution. Overall, the SFRs based on the linear L_W3–SFR calibrations show a nonlinear trend between the estimated SFR and SFR_true.

As shown in Table 2 the calibrations of Jarrett et al. (2013), and Chang et al. (2015) tend to underestimate the SFR on average for high SFR galaxies. The 68% percentiles of the ratio between these indicators and the true SFR when they are applied to the model SEDs shows that they tend to overpredict the SFR, especially in the case of low-SFR galaxies (Table 1). This bias is significantly reduced (but not absent) in higher SFR galaxies. This behavior is interpreted as the result of older SPs contributing to the emission of the MIR bands we consider. The calibrations of Jarrett et al. (2013), and Chang et al. (2015) are not shown in Fig. 5 in order to avoid congestion.

4.4. Mid-IR emission dominated by old stellar populations

Although the data used to derive the parameters of Eq. (4) are based on SED fitting tools built on the principle of energy balance, the large scatter in the W3 luminosity at low SFR (SFR < 0.01 M_⊙ yr⁻¹) combined with the negative term in Eq. (4) (introduced to account for the contribution of old SPs) may lead to negative SFR for galaxies with low SFR or low sSFR. This is because Eq. (4) gives an average relation between SFR, MIR, and NIR luminosity while the significant scatter resulting from the wide variety of SFHs and ISM parameters may result in stronger NIR emission with respect to the MIR emission. This effect may be exacerbated by the presence of stochastically heated dust (which is more likely in low sSFR galaxies). The 2–4 μm bands include a strong dust emission feature that can be significantly excited by the general interstellar radiation field (e.g., Camps et al. 2015). For these cases, the SFR can not be calculated through Eq. (4).

The comparisons with the test dataset (Sect. 4.1) show that sources showing negative SFRs using Eq. (4) have on average L_old ≃ 65 L_young. The results are similar for both the extinction-dependent and extinction-independent relations. For the extinction-independent calibration the cutoff value for the W3–W1 color, below which galaxies give negative SFRs, is log(L_W3/L_W1) = − 0.371, while this is not constant for the extinction-dependent relations. The galaxies whose SFR can not be measured are relatively few (∼6%, and ∼12% for the extinction-independent and extinction-dependent relations respectively) considering that this analysis included galaxies with extremely low sSFR down to sSFR ≃10⁻¹³ M_⊙ yr⁻¹/M_⊙. Therefore, galaxies that yield negative SFRs based on Eq. (4) have extremely small star-forming activity and their IR emission is dominated by old SPs. Overall, the Eq. (4) calibrations offer reliable estimations even in low SFRs.

Additionally, we examine the capability of the IR–SFR calibrations to recover the SFR throughout the range of SF activity. Figure 6 shows the ratio of the estimated SFR, from the proposed and existing IR–SFR calibrations, and the true SFR as a function of the model-galaxies sSFR. Previous works already showed that in passive galaxies with sSFR ≲ 10⁻¹¹ M_⊙ yr⁻¹/M_⊙ the dust is not heated by young star and thus, SFR estimations based on the MIR emission for these galaxies are not reliable (e.g., Kennicutt 1998; Salim et al. 2009, 2016). Figure 6 shows that all the IR–SFR calibrations compared here fail to recover and overestimate the SFR in galaxies with sSFR ≲ 10^−10.7 M_⊙ yr⁻¹/M_⊙. However, the proposed calibrations from this work, and especially the extinction dependent, tend to reduce the biases in the SFR estimation. This is demonstrated in Fig. 6 by the slightly better agreement between the modes of the proposed calibrations (red and blue lines) with the equality (black dashed) line in relatively low sSFRs (sSFR ≃ 10^−10.5 M_⊙ yr⁻¹/M_⊙) and in high sSFRs (sSFR ≃ 10⁻⁸ M_⊙ yr⁻¹/M_⊙).

Fig. 6. Logarithm of the ratio of the estimated SFR based on the IR emission over the true SFR as given by the CIGALE models, as a function of their true specific SFR (sSFR). The comparisons show the extinction-dependent and extinction-independent calibrations of Eq. (4), and the LIR–SFR calibrations of Jarrett et al. (2013), Chang et al. (2015), and Cluver et al. (2017) with blue, red, purple, yellow, and green lines respectively for model galaxies separated in groups of low (left), moderate (center), and high (right) extinction.

4.5. Stellar mass calibrations

Similarly to the SFR-L_IR relation, we explore the dependence of the mass-to-light ratio on the SP age and dust content/extinction. For this analysis, we use the same data as with the SFR analysis, as described in Sect. 3.2. Previous studies have identified the SP age as the main cause of uncertainty in the calculation of the M_⋆ (e.g., Rix & Rieke 1993). Here, as an observational tracer of the SP age, we adopt the u − r, and g − r colors, which have also been used in other studies (e.g., Bell et al. 2003; Wen et al. 2013). However, this work also explores the possibility of increased scatter in the inferred mass-to-light ratio induced by the use of extinction-corrected colors as is usually the case. Dust emission, heated by the star-forming activity is mainly in the mid, and far-IR parts of the spectrum.

However, emission from PAH molecules and stochastically heated dust grains may have a significant contribution to the NIR emission. The relation between the mass-to-light ratio with the optical colors shows a large scatter (e.g., Fig. 7). This may be the result of extinction affecting the u − r, and g − r colors of the SPs creating multiple mass-to-light ratio tracks for different values of extinction. In fact, if we account for the extinction (i.e., correct the u − r color for the extinction based on the E(B − V)) the scatter is slightly reduced. The advantage of our analysis is that it accounts for the effect of dust in two ways: (a) by its effect on the colors of the SPs and (b) by its effect on the NIR emission through stochastically heated dust. There is a tertiary effect originating from the contribution of the Hα emission in the r band photometry. This is also indirectly accounted for through the (loose) correlation between extinction in the ISM and star-forming activity.

Fig. 7. Mass-to-light ratio as a function of u − r (left column) and g − r (right column) color for the SDSS-matched sample of the modeled galaxies. In the bottom row, the u − r and g − r colors are corrected for extinction while in the top row, they are not. The galaxies are color-coded based on their E(B − V). The dashed magenta line represents the extinction-independent relation of Eq. (6), and the dashed-dotted lines the extinction-dependent, drawn for some discrete E(B − V) values as shown in the legend using the parameters presented in Table 3. The dotted black line shows the relations of Wen et al. (2013).

Figure 7 shows the mass-to-light ratio of the CIGALE model galaxies, as a function of the u − r color. Each point represents a single model SED and it is color-coded based on its extinction. Figure 7 reveals the scatter induced by the differences in extinction, where especially for blue galaxies it can be more than an order of magnitude. The scatter increases for decreasing u − r color. However, in highly star-forming galaxies (u − r ≲ 1.5 or g − r ≲ 0.75) the scatter is by far larger. In these color regimes, galaxies have higher SFRs, and therefore, IR emission from the dust is brighter, and its respective dependence on extinction is stronger. Moreover, the g − r color has an overall narrower range and shows a larger scatter as a function of the mass-to-light ratio compared to u − r, and therefore is not as good as a SP age indicator for the purposes of this work.

We describe the correlation between the mass-to-light ratio and the optical color as a power law with negative power values, thus as M/L = α(u − r)^β. Moreover, in this analysis, we fit the mass-to-light ratio relation with the optical color taking also into account the extinction. We include the effect of dust by parametrizing the parameter β as a linear function of extinction with the form β = δ′+ϵ′ e where δ′ and ϵ′ are fitted parameters and e is an extinction metric. We adopt as an extinction metric (e) the color excess E(B − V) or the IRX index (as defined in Sect. 3). However, because extinction also affects the u − r color we also include an extinction term in the parameter α which is best parametrized with a second-order polynomial of the extinction metric: α = α′+β′e + γ′e². In order to account for cases where there are no available extinction measurements we also fit the data with an extinction-independent parametrization. These fits involve the SED models that include the effect of extinction but follow the E(B − V) distribution of SDSS galaxies (Sect. 3.2).

The final mass-to-light calibrations are presented in Eq. (6):

$$\begin{array}{cc}\hfill \log \frac{M_{\star }/({\mathrm{M}}_{\odot })}{L_{\mathrm{NIR}}}& =\alpha {(x+1)}^{\beta }+0.5\hfill \\ \hfill x& =u-r,\mathrm{or}x=g-r\hfill \\ \hfill L_{\mathrm{NIR}}& =\frac{L_{\mathrm{W}1}}{(L_{\odot ,\mathrm{W}1})},\mathrm{or}{L}_{\mathrm{NIR}}=\frac{L_{F200W}}{(L_{\odot })}\hfill \\ \hfill \alpha & ={\alpha }^{\prime }+{\beta }^{\prime }e+{\beta }^{\prime }{e}^2\hfill \\ \hfill \beta & ={\delta }^{\prime }+{ϵ}^{\prime }e\hfill \\ \hfill e& =E(B-V),\mathrm{or}e=\mathrm{IRXu},\mathrm{or}e=\mathrm{IRXg}.\hfill \end{array}$$(6)

L_W1 is measured in the in-band equivalent solar luminosity (see e.g., Jarrett et al. 2013). It is equal to scaling νL_ν by 22.883. The L_F200W is measured in solar luminosities. The best-fit results based on the MCMC fitting method for the parameters α,  β and α′, β′, γ′, δ′, ϵ′ are given in Table 3. Moreover, we performed a RF model for the same datasets and the same sets of features with the MCMC in order to independently compare the result of the MCMC method.

As mentioned in Sect. 4.2, because the SEDs are an output of modeling, their photometries, and their stellar masses come without uncertainties, and thus they are not included in the likelihoods of the MCMC fitting. The likelihood function for the extinction-independent model is:

$$\begin{array}{c}\hfill \log p(y|x,\alpha ,\beta )=-\frac{1}{2}{\displaystyle \sum _n}{(y_n-\alpha x_n^{\beta })}^2\end{array}$$

and for the extinction–dependent model is:

$$\begin{array}{cc}& \log p(y|x,e,{\alpha }^{\prime },{\beta }^{\prime },{\beta }^{\prime },{\delta }^{\prime },{ϵ}^{\prime })\hfill \\ \hfill & =-\frac{1}{2}{\displaystyle \sum _n}{[y_n-({\alpha }^{\prime }+{\beta }^{\prime }e+{\beta }^{\prime }{e}^2)x_n^{({\delta }^{\prime }+{ϵ}^{\prime }e)}]}^2\hfill \end{array}$$(7)

where:

$$\begin{array}{cc}& y=\log \frac{M_{\star }/L_{\mathrm{W}1}}{(M_{\odot }/L_{\odot ,\mathrm{W}1})}-0.5,\mathrm{or}y=\log \frac{M_{\star }/L_{F200W}}{(M_{\odot }/L_{\odot })}-0.5\hfill \\ \hfill & x=u-r+1(\mathrm{mag}),\mathrm{or}x=g-r+1(\mathrm{mag})\hfill \\ \hfill & e=E(B-V),\mathrm{or}e=\mathrm{IRXu},\mathrm{or}e=\mathrm{IRXg}.\hfill \end{array}$$

4.6. Comparisons with the stellar masses of the model galaxies

Figure 7 shows the extinction-independent relation of Eq. (6) with a dashed magenta line. The dashed-doted lines show the extinction-dependent relation for some representative extinction values. The extinction-independent relation follows the median position of the points. The lines drawn based on the extinction-dependent relation follow well the data points of similar extinction values. As a reference, we also compare with the calibration of Wen et al. (2013) which is empirically calibrated based on the W1 emission and stellar masses estimated from the MPA-JHU analysis. The mass-to-light ratio in this relation is parametrized as a polynomial function of the u − r or g − r colors. The relation of Wen et al. (2013) shows good agreement for average u − r colors (u − r ≃ 1.5), as well as, for moderate extinction [E(B − V)≃0.5]. However, it fails to recover the M_⋆ for galaxies with low or high extinction, while it shows larger discrepancies near the lower and upper limits of u − r color.

Figure 8 shows the logarithm of the ratio between the calculated stellar masses for the MCMC, and RF models (using the WISE band-1) over the true M_⋆, as a function of the true M_⋆. The true M_⋆ is an output of CIGALE based on the given SFHs of the modeled galaxies. We also plot the calculated stellar masses based on the relation of Wen et al. (2013) as a comparison reference. Figure 8 does not show the results based on the JWST-NIR F200W band for clarity because these are very similar to those using the WISE band-1. Table 4 gives the comparisons for the proposed relations for all the models considered in this analysis, separated in low u − r (younger SPs or lower extinction), and high u − r (older SPs or higher extinction) color. These comparisons are also given separately for the sample having a uniform distribution in extinction, and the SDSS-matched sample (Sect. 3.2) in order to reveal possible biases based on the selected sample of galaxies.

Fig. 8. Comparison between the true M⋆ as given by the CIGALE output, and the best-fit models of Eq. (6) (MC; Table 3), and the random forest machine learning method (RF) using the WISE band-1 and the u − r color. The y-axis represents the logarithm of the ratio between the M⋆ of each model and the true M⋆, thus the black dashed line at zero represents equality. Blue and red colors show the extinction-independent, and extinction-dependent relation of Eq. (6) respectively. With light-blue and orange colors are the extinction-independent, and extinction-dependent models based on the RF method respectively. The M⋆ calculated based on the mass-to-light conversion from Wen et al. (2013) is also shown with green color. The lines represent the distribution modes and the shaded areas the 68% percentiles as a function of the log M⋆ true. Each panel shows the comparison of galaxy models with one of the discrete E(B − V) values, shown in the bottom left corner of each panel. This comparison is performed with the uniform-E(B − V) sample, therefore, there is an equal number of sources for each extinction value. At the bottom of each panel are the modes and 68% percentiles of the log M⋆model/M⋆, True distribution for all galaxy models with the specific E(B − V), shown with corresponding colors to the relations used to derive the M⋆. These comparisons are also summarized in Table 4.

The MCMC relations show robust estimations throughout the M_⋆ range. Figure 8 demonstrates that the use of the extinction feature in the M_⋆ calculations is important in low and high-extinction galaxies. The extinction-independent relation of Eq. (6) is in good agreement for mid-range and high-extinction sources. However, in all cases, the scatter is increased when the extinction is not taken into account. Overall, the RF yields similar results to the MCMC indicating the results of the proposed calibrations are independent of the method.

5. Comparison with observations

The comparisons of the SFR and mass-to-light calibrations of Eqs. (4) and (6) show excellent results in retrieving the SFR and M_⋆ of the mock galaxies used in our analysis. However, in order to test these calibrations in a more realistic setting we compare them with observations of galaxies and their SFRs and stellar masses inferred from different methods.

5.1. Comparison of the star formation rate calibrations

Figure 9 shows comparisons between SFRs based on our method of extinction-dependent and extinction-independent calibrations using Eq. (4) and the best-fit relations in Table 1 with various methods including: (a) SED fitting from Chang et al. (2015), Salim et al. (2018), or MPA-JHU (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004), and (b) those based on monochromatic WISE band-3 emission (Jarrett et al. 2013; Chang et al. 2015; Cluver et al. 2017). These comparisons involve star-forming and unclassified but exclude LINER and AGN galaxies as they were classified by MPA-JHU (BPTCLASS ≠ 4 and BPTCLASS ≠ 5; Brinchmann et al. 2004). The adopted IR fluxes for the WISE bands are from the AllWISE source catalog² as given by Chang et al. (2015) before the photometric corrections they applied. Distances used to calculate the luminosities were based on the redshifts provided by the MPA-JHU catalog. Their E(B − V) was estimated through the Balmer extinction based on the Hα and Hβ fluxes from MPA-JHU and the conversion of Domínguez et al. (2013). Galaxies with signal-to-noise ratio (S/N) < 5 in the WISE bands 1 and 3 fluxes and the E(B − V) were omitted ensuring good quality data.

Fig. 9. Comparisons between SFRs estimated through (a) the extinction-dependent and extinction-independent relations of Eq. (4) (Table 1) which utilize WISE bands 1 and 3, b) calibrations based only on WISE band-3 (Jarrett et al. 2013; Chang et al. 2015; Cluver et al. 2017), and (c) based on SED fitting from Chang et al. (2015), Salim et al. (2018) and MPA-JHU (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004). In each subplot, the ordinate indicates the logarithm of the ratio (log (SFRy/SFRx)) between the compared SFR estimations, and the abscissa the logarithm of the SFR as indicated in the axes labels. Red error bars represent the modes and 68% percentiles of the distribution of galaxies within bins of log SFR/M⊙ yr−1 ≃ 0.4 (12 bins within the range −2.5 < log SFR/M⊙ yr−1 < 2.5). The log SFRy/SFRx modes and 68% percentiles for all the galaxies compared in each panel are shown at the top of each subplot. The points are color-coded based on the logarithm of the absolute value of the equivalent width (EW) of the Hα line of the galaxies. The black dashed line represents equality.

The format of Fig. 9 is similar to that of Fig. 5 where the ordinate indicates the logarithm of the ratio of the compared SFRs (log SFR_y/SFR_x), and the abscissa the logarithm of the SFRs as indicated in the axes labels. Therefore, the equality lies at value zero of the y-axis. The points in Fig. 9 are color-coded based on the absolute value of the equivalent width of the Hα line (EW_Hα) of the galaxies. The Hα EW is strongly correlated with the sSFR of the galaxies (e.g., Belfiore et al. 2018) and it is used here as an independent method for tracing their current star-forming activity with respect to their stellar component.

The SED fitting of Chang et al. (2015) included optical (SDSS), and IR (WISE), but not UV photometry. The SED fitting of Salim et al. (2018) was based on a SED + L_IR fitting method utilizing except the SDSS optical and GALEX UV in the SED fitting, corrections based on the 22 μm or 12 μm photometry from unWISE. The SFR estimations from MPA-JHU were based on the galaxies’ emission lines and the method of Kauffmann et al. (2003) for galaxies within the SDSS spectroscopic fiber. For galaxies that did not fit within the fiber, SFRs were estimated using the ugriz photometry based on SED fitting following Salim et al. (2007).

The proposed SFR calibrations of Eq. (4) show the best on-average agreement compared to the SFRs from SED fitting. The extinction-dependent relation shows an overall slightly better agreement with the other indicators compared to the extinction-independent calibration. Moreover, Eq. (4) calibrations show good agreement throughout the SFR range compared to the SFRs from Salim et al. (2018) with on-average no large deviations in the low or high SFRs. As expected, the extinction-dependent relation shows better behavior in the lower and higher SFR regimes. The scatter between Eq. (4) calibrations compared to results from SED fitting is reduced with respect to other IR-based SFRs.

The comparisons between the monochromatic IR–SFR calibrations (Jarrett et al. 2013; Chang et al. 2015; Cluver et al. 2017) reveals some differences between them. While the best agreement is between the calibrations of Jarrett et al. (2013), and Chang et al. (2015), the largest difference is between Jarrett et al. (2013), and Cluver et al. (2017; about 0.5 dex on average). As also shown from the comparisons with the mock photometry used in our analysis (Fig. 5), comparing Eq. (4) with these calibrations show differences ranging from −0.25 dex (compared to Cluver et al. 2017) to +0.25 dex (compared to Jarrett et al. 2013) on average and nonlinear dependence on the value of the SFR.

As revealed by the color-coding of Fig. 9, the galaxies where the Eq. (4) calibrations infer lower SFRs compared to other IR–SFR calibrations have very low EW (EW_Hα ≲ 20 Å) indicating that these are passive galaxies which are mainly dominated by older SPs. Moreover, only the extinction-dependent calibration of Eq. (4) does not show a monotonic EW_Hα gradient in the comparisons with the SED-fitting SFRs. This indicates the ability of the proposed calibration to account for the contribution of older SPs and extinction in the dust heating while all other IR–SFR calibrations tend to overestimate the SFR of passive galaxies and underestimate the SFR of highly star-forming galaxies.

The SED fitting and the use of emission lines have been proven powerful methods for deriving the characteristics of galaxies. However, not using simultaneously UV and IR photometry can lead to biases with respect to the estimation of extinction due to not covering the energy balance between the absorbed UV light, and that which is re-emitted in the IR (e.g., Lanz et al. 2013). This can be more important in the extremes of the SFR, and M_⋆ range, for instance in dust-free galaxies when the UV is not taken into account, or in dusty/high-SFR galaxies when the IR emission is not taken into account. The fact that the calibrations of Eq. (4) show excellent agreement with the SFRs from SED fitting demonstrates the advantage of this method with respect to monochromatic IR SFR indicators and shows that their results can be considered robust over a wide SFR range.

5.2. Comparison of the mass-to-light ratio calibrations

Figure 10 shows comparisons between stellar masses calculated using various methods including: (a) the extinction-dependent, and extinction-independent relations of Eq. (6) which combine the WISE band 1 and the u − r, or g − r color, (b) SED fitting from Chang et al. (2015), Salim et al. (2018), and MPA-JHU (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004), and c) the mass-to-light calibration of Wen et al. (2013) which is using the WISE band-1 and the u − r color. Similarly to Fig. 9, this comparison involves star-forming or unclassified galaxies but excludes LINER and AGN galaxies as they have been classified by MPA-JHU (BPTCLASS ≠ 4 and BPTCLASS ≠ 5). It involves galaxies with S/N > 5 in the L_W1, SDSS u, g, r magnitudes, and E(B − V) extinction metric. WISE band fluxes were based on the AllWISE source catalog and were adopted from Chang et al. (2015). The E(B − V) was estimated through the Balmer extinction based on the Hα and Hβ fluxes from MPA-JHU and the conversion of Domínguez et al. (2013). Distances used to calculate the luminosities were based on the redshifts provided by the MPA-JHU catalog.

Fig. 10. Comparisons between stellar masses estimated through (a) the extinction-dependent, and extinction-independent relations of Eq. (6) based on the MCMC fitting (Table 3) which utilize the WISE band-1 photometry and the u − r or the g − r color, (b) the calibration of Wen et al. (2013) which uses the WISE band 1 and the u − r color and (c) based on SED fitting from Chang et al. (2015), Salim et al. (2018), and MPA-JHU (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004). In each subplot, the ordinate indicates the logarithm of the ratio (log (M⋆, y/M⋆, x)) between the compared M⋆ estimations, and the abscissa the logarithm of the M⋆ as indicated in the axes labels. Red error bars represent the modes and 68% percentiles of the distribution of galaxies within bins of log M⋆/M⊙ ≃ 0.6 (8 bins within the range 7 < log M⋆/M⊙ < 12). The log M⋆y/M⋆x modes and 68% percentiles for all the galaxies compared in each panel are shown at the top right of each subplot. The points are color-coded based on the galaxies’ nebular E(B − V). The black dashed line represents equality.

Overall, these comparisons show that calibrations using monochromatic NIR photometry and an optical color agree on average with the stellar masses based on SED fitting with a small standard deviation of about 0.1 dex. The extinction-independent relation is in excellent agreement with the stellar masses derived using the Wen et al. (2013) calibration showing ∼ − 0.1 dex lower M_⋆ on average without strong deviations in low or high stellar masses. However, it shows a slightly better agreement in low and high M_⋆ galaxies compared to results based on SED fitting, while the u − r relation of Wen et al. (2013) shows a small overestimation for galaxies with large M_⋆ and high extinction.

The extinction-dependent relation of Eq. (6) is in good agreement with stellar masses estimated through SED fitting. These comparisons show slightly better agreement throughout the M_⋆ range, even for very high and low-SFR galaxies where the effect of reddening and SP age becomes gradually more important. Moreover, as indicated by the color code of Fig. 10 the extinction-dependent relation of Eq. (6) tends to correct the overestimation of high extinction galaxies, and vice-versa the underestimation of low extinction galaxies, compared to extinction independent relations. The comparison with the model SEDs, where we know a priori the true M_⋆ of the galaxies (Sect. 4.5), indicates that the inferred stellar masses based on the extinction-dependent relation of Eq. (6) and on SED fitting are closer to the correct M_⋆. These comparisons also demonstrate the importance of taking extinction into account in the estimation of the stellar mass.

6. Conclusions

We used the SED fitting code CIGALE to create a large sample of mock galaxy SEDs that covered a wide range of SFRs, stellar masses, and ISM conditions. We compared the SFR and M_⋆ of the model galaxies with their IR luminosity for photometric bands of WISE and the JWST MIRI and NIRCam instruments. This analysis showed that there is a strong dependence of the IR tracers of galaxies’ SFR and M_⋆ on the age of their SPs, and their extinction. However, the combined use of IR bands that trace SPs of different ages can provide more robust calibrations and minimize the biases and the scatter they introduce in the SFR, and M_⋆ calibrations. This is particularly important for low-SFR galaxies where the contribution of old SPs can be significant when measuring their SFR, and for high-SFR galaxies where the contribution of dust emission can bias the M_⋆ measurements. The addition of extinction-dependent calibrations also offers more reliable SFR, and M_⋆ estimations with less scatter which are particularly better for low-extinction/dust-free (e.g., dwarf) galaxies.

In summary, this work provides extinction-dependent and extinction-independent calibrations and quantification of their scatter for measuring the:

1. SFR based on the WISE band-3 and band-1 utilizing as extinction indicators the E(B − V), the L_W4/L_u, and the L_W4/L_g ratio (Eq. (4), Table 1);

2. SFR based on the JWST NIR-F200W and MIRI-F2100W bands utilizing as extinction indicator the E(B − V) (Eq. (4), Table 1);

3. Mass-to-light ratio based on the WISE band-1 and the u − r or the g − r color utilizing as extinction indicators the E(B − V), the L_W4/L_u, and the L_W4/L_g ratio (Eq. (6), Table 3);

4. Mass-to-light ratio based on the JWST NIR-F200W and the u − r or g − r color utilizing as extinction indicator the E(B − V) (Eq. (6), Table 3).

The comparisons with both modeled and observed samples of galaxies show that the proposed calibrations offer robust SFR, and M_⋆ estimations for a wide range of these values, minimizing the scatter. The random forest analysis yields similar results to the MCMC method showing that the provided calibrations are method independent.

Acknowledgments

The authors thank the anonymous referee for providing comments and suggestions that improved the clarity of this work. KK is supported by the project “Support of the international collaboration in astronomy (Asu mobility)” with the number: CZ.02.2.69/0.0/0.0/18_053/0016972. Asu mobility is co-financed by the European Union. KK and AZ acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement n. 617001 (A-BINGOS). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie RISE action, grant agreements No 691164 (ASTROSTAT), and No 873089 (ASTROSTAT-II). EK acknowledges support from the Public Investments Program through a Matching Funds grant to the IA-FORTH. SS is supported through NASA awards NNX12AE06G and 80NSSC20K0440. JS acknowledges financial support from the Czech Science Foundation under Project No. 22-22643S. This research has made use of: (a) data products from the Wide-field Infrared Survey Explorer (WISE), which is a joint project of the University of California, Los Angeles, and JPL, California Institute of Technology, funded by NASA; (b) observations made with the Spitzer Space Telescope, which was operated by JPL, California Institute of Technology under a contract with NASA; (c) the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory (JPL), California Institute of Technology, under contract with NASA; (d) the NASA/IPAC Infrared Science Archive (IRSA), which is funded by NASA and operated by the California Institute of Technology; Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the US Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

References

Appendix A: Appendix A

Table A.1 shows the configuration of the CIGALE SED creation modules. See text of Boquien et al. (2019) for a detailed description of the parameters.

All Tables

All Figures

Fig. 1. Intensity as a function of wavelength for the total attenuated intensity (gray), and the components of the young stellar (blue), old stellar (orange), and dust emission (red) for two model galaxies: (a) in the top panel is a galaxy with low star-forming activity (SFR = 0.01 M⊙ yr−1, sSFR = 10−11.75 M⊙ yr−1/M⊙) and low-extinction (E(B − V) < 0.001), and (b) in the bottom panel is a dusty (E(B − V) = 0.55) galaxy with high star-forming activity (SFR = 68 M⊙ yr−1, sSFR = 10−8.5 M⊙ yr−1/M⊙). The dashed-dotted lines with brown and purple colors represent the WISE bands 1 and 3 transmission bandpasses respectively. The dotted fuchsia and red lines represent the JWST NIR-F200W, and MIRI-F2100W bandpasses respectively.

In the text

Fig. 2. WISE band-3 luminosity (LW3) per M⋆ unit, as a function of sSFR for the generated library of galaxy SEDs. All panels show the same sources but in the top left they are color-coded based on their E(B − V), in the top right based on their metallicity (Z), in the bottom left on their infrared excess (IRXu = log LW4/Lu), and in the bottom right based on their infrared excess using the g band instead of u (IRXg = log LW4/Lg).

In the text

Fig. 3. SFR, M⋆ plane for the final sample of CIGALE model galaxies used for the rest of the analysis and fitting. The color code indicates their specific SFR [sSFR ≡ SFR/M⋆ (M⊙ yr−1/M⊙)]. The best linear fit is represented with a red line. The green dashed, yellow dashed-dotted, magenta dotted, and blue dashed-dotted lines represent the main sequence for diverse samples of local Universe galaxies from Elbaz et al. (2007), Speagle et al. (2014), Popesso et al. (2019), and Kouroumpatzakis et al. (2021) respectively.

In the text

Fig. 4. Distribution of E(B − V) for the uniform (green histogram) and the SDSS-matched samples (red histogram) of model galaxies. With blue color is the nebular E(B − V) distribution based on the Balmer decrement for sources in the SDSS spectroscopic MPA-JHU catalog considering only star-forming galaxies with E(B − V)/E(B − V)err > 3. The ordinate corresponds to density which is equal to the ratio of the count of each bin to the total number of counts multiplied by the bin width with the integral of the area under the histogram being equal to one.

In the text

Fig. 5. Comparison between the true SFR as given by the CIGALE output, and the best-fit models of Eq. (4) (MCMC; Table 1), and the random forest machine learning method (RF) using the WISE bands 1 and 3. The y-axis represents the logarithm of the ratio between the SFR of each model and the true SFR, thus the black dashed line at zero represents equality. Blue, red, pink, and purple colors show the extinction-independent, and extinction-dependent (E(B − V), IRXu, IRXg) relations of Eq. (4) respectively. With light-blue and orange colors are the extinction-independent, and extinction-dependent models based on the RF method respectively. The linear log LW3–log SFR calibration by Cluver et al. (2017) is also displayed in green color. The lines represent the distribution modes and the shaded areas show the 68% percentiles of the log (SFRmodel/SFRtrue) as a function of the log SFRtrue. Each panel shows the comparison of galaxy models with one of the discrete E(B − V) values, shown in the bottom left corner of each panel. This comparison is performed with a uniform-E(B − V) sample, therefore, there is an equal number of sources for each extinction value. At the bottom of each panel are the modes and 68% percentiles of the log (SFRmodel/SFRTrue) distribution for all galaxy models with the specific E(B − V), shown with corresponding colors to the relations used to derive the SFR. These comparisons are also summarized in Table 2.

In the text

Fig. 6. Logarithm of the ratio of the estimated SFR based on the IR emission over the true SFR as given by the CIGALE models, as a function of their true specific SFR (sSFR). The comparisons show the extinction-dependent and extinction-independent calibrations of Eq. (4), and the LIR–SFR calibrations of Jarrett et al. (2013), Chang et al. (2015), and Cluver et al. (2017) with blue, red, purple, yellow, and green lines respectively for model galaxies separated in groups of low (left), moderate (center), and high (right) extinction.

In the text

Fig. 7. Mass-to-light ratio as a function of u − r (left column) and g − r (right column) color for the SDSS-matched sample of the modeled galaxies. In the bottom row, the u − r and g − r colors are corrected for extinction while in the top row, they are not. The galaxies are color-coded based on their E(B − V). The dashed magenta line represents the extinction-independent relation of Eq. (6), and the dashed-dotted lines the extinction-dependent, drawn for some discrete E(B − V) values as shown in the legend using the parameters presented in Table 3. The dotted black line shows the relations of Wen et al. (2013).

In the text

Fig. 8. Comparison between the true M⋆ as given by the CIGALE output, and the best-fit models of Eq. (6) (MC; Table 3), and the random forest machine learning method (RF) using the WISE band-1 and the u − r color. The y-axis represents the logarithm of the ratio between the M⋆ of each model and the true M⋆, thus the black dashed line at zero represents equality. Blue and red colors show the extinction-independent, and extinction-dependent relation of Eq. (6) respectively. With light-blue and orange colors are the extinction-independent, and extinction-dependent models based on the RF method respectively. The M⋆ calculated based on the mass-to-light conversion from Wen et al. (2013) is also shown with green color. The lines represent the distribution modes and the shaded areas the 68% percentiles as a function of the log M⋆ true. Each panel shows the comparison of galaxy models with one of the discrete E(B − V) values, shown in the bottom left corner of each panel. This comparison is performed with the uniform-E(B − V) sample, therefore, there is an equal number of sources for each extinction value. At the bottom of each panel are the modes and 68% percentiles of the log M⋆model/M⋆, True distribution for all galaxy models with the specific E(B − V), shown with corresponding colors to the relations used to derive the M⋆. These comparisons are also summarized in Table 4.

In the text

Fig. 9. Comparisons between SFRs estimated through (a) the extinction-dependent and extinction-independent relations of Eq. (4) (Table 1) which utilize WISE bands 1 and 3, b) calibrations based only on WISE band-3 (Jarrett et al. 2013; Chang et al. 2015; Cluver et al. 2017), and (c) based on SED fitting from Chang et al. (2015), Salim et al. (2018) and MPA-JHU (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004). In each subplot, the ordinate indicates the logarithm of the ratio (log (SFRy/SFRx)) between the compared SFR estimations, and the abscissa the logarithm of the SFR as indicated in the axes labels. Red error bars represent the modes and 68% percentiles of the distribution of galaxies within bins of log SFR/M⊙ yr−1 ≃ 0.4 (12 bins within the range −2.5 < log SFR/M⊙ yr−1 < 2.5). The log SFRy/SFRx modes and 68% percentiles for all the galaxies compared in each panel are shown at the top of each subplot. The points are color-coded based on the logarithm of the absolute value of the equivalent width (EW) of the Hα line of the galaxies. The black dashed line represents equality.

In the text

Fig. 10. Comparisons between stellar masses estimated through (a) the extinction-dependent, and extinction-independent relations of Eq. (6) based on the MCMC fitting (Table 3) which utilize the WISE band-1 photometry and the u − r or the g − r color, (b) the calibration of Wen et al. (2013) which uses the WISE band 1 and the u − r color and (c) based on SED fitting from Chang et al. (2015), Salim et al. (2018), and MPA-JHU (Kauffmann et al. 2003; Brinchmann et al. 2004; Tremonti et al. 2004). In each subplot, the ordinate indicates the logarithm of the ratio (log (M⋆, y/M⋆, x)) between the compared M⋆ estimations, and the abscissa the logarithm of the M⋆ as indicated in the axes labels. Red error bars represent the modes and 68% percentiles of the distribution of galaxies within bins of log M⋆/M⊙ ≃ 0.6 (8 bins within the range 7 < log M⋆/M⊙ < 12). The log M⋆y/M⋆x modes and 68% percentiles for all the galaxies compared in each panel are shown at the top right of each subplot. The points are color-coded based on the galaxies’ nebular E(B − V). The black dashed line represents equality.

In the text