Accuracy of Stellar Mass-to-light Ratios of Nearby Galaxies in the Near-Infrared

To utilize the SED derived from the observed photometric dataset, we first employ the sample from the DustPedia (Davies et al., 2017; Clark et al., 2018). The multi-wavelength dataset of the DustPedia spanning the UV to far-infrared (FIR) enables robust estimations of both stellar population and dust properties through SED fitting with the Code Investigating GALaxy Emission (CIGALE) code (Boquien et al., 2019). This characteristic makes the DustPedia sample ideal for our study (Nersesian et al., 2019). The DustPedia constructed the SEDs based on the aperture-matched photometry on the 42 bands data covering UV to sub-mm, obtained from the Hershel observation, along with the archival observations from the GALEX, SDSS, 2MASS, WISE, Spitzer, and Planck. The original DustPedia sample contains 875 nearby galaxies with $D_{L}\leq 79$ Mpc. However, this DustPedia sample is slightly biased toward low-mass galaxies ($\langle M_{*}\rangle=10^{9.7\pm 0.95}M_{\odot}$) because of their proximity.

To compensate for this limitation, we additionally adopt a sample of nearby ($0.01\leq z\leq 0.11$) massive galaxies from Li et al. (2023), which employed data from the SDSS Stripe 82 region (S82). Li et al. (2023) conducted the sophisticated photometry on the multi-wavelength imaging data of 2685 massive galaxies ($10^{10}M_{\odot}\lesssim M_{*}\lesssim 10^{11.5}M_{\odot}$). As the initial sample selection was based on stellar mass estimates from the GALEX–SDSS–WISE Legacy Catalog (GSWLC-2; Salim et al., 2018), a small fraction of the S82 sample finally shows $M_{*}<10^{10}M_{\odot}$. These data encompassed UV, optical, NIR, MIR, and FIR data obtained from GALEX, SDSS, 2MASS, WISE, and Herschel-SPIRE, respectively. The matched-aperture and profile-fitting photometry were applied for the shorter-wavelength bands ranging from UV to W2 and the longer-wavelength bands, respectively. In addition, multi-band imaging decomposition was carefully performed to accurately remove the contribution from neighboring galaxies or foreground stars. These procedures resulted in reliable flux measurements and precise estimations of the uncertainties of the multi-wavelength dataset. From the comparison with the photometric data of GALEX, SDSS, and 2MASS extended source catalogs (Jarrett et al., 2000), Li et al. (2023) showed that their measurements agree within $\pm 0.1$ mag in the UV, optical, and NIR bands. However, the photometric measurements of Li et al. (2023) are systematically brighter ($\sim 0.2-1.0$ mag) than those in 2MASS point source catalog (Cutri et al., 2003) and ALLWISE, in which the photometric data were originally estimated from a PSF-fitting method. Notably, these systematic discrepancies are significantly reduced up to $\leq 0.2$ mag when compared with the unWISE catalog (Lang et al., 2016).

For the S82 sample, stellar and dust masses are drawn from Li et al. (2023). These parameters were obtained via SED fitting using the CIGALE. In this SED fitting process, the stellar SED was modeled using the simple stellar population (SSP) from Bruzual & Charlot (2003). The stellar population was modeled with a double-exponential SFH, comprising an old stellar population with an age of 12 Gyr and a young stellar population with an age range from 10 to 5000 Myr. The stellar emission was modeled with metallicities of $0.004-0.02$. The IMF of Chabrier (2003) was adopted. In addition, the starburst attenuation curve from Calzetti et al. (2000) with a modification of the power-law slope was adopted to account for the dust attenuation in the SED fitting.

Dust emission is modeled with $U_{\rm min}=0.1-25$, $q_{\rm PAH}=2.5$ (fixed), $\gamma=0.001-0.2$. Here, $U_{\rm min}$ denotes the minimum radiation field from the stars; $q_{\rm PAH}$ represents the mass fraction of PAHs; $\gamma$ indicates the mass fraction of dust illuminated from the minimum to maximum radiation field. $\alpha$, the power-law slope of the dust continuum over the IR wavelength, is fixed to 2. Nebular emission lines were also included (Inoue, 2011). The AGN component modeled with Fritz et al. (2006) was incorporated. Here, the AGN fraction ($f_{\rm AGN}$) was calculated as the ratio of the AGN luminosity to the total IR luminosity estimated within $1-1000$ $\mu$m.

Whereas the SED fitting with the CIGALE was also applied to the DustPedia sample (Jones et al., 2013; Nersesian et al., 2019), parameters for the stellar population and dust modeling significantly differ from those of Li et al. (2023). In addition, contrary to the S82 sample, the SED fitting of the DustPedia sample did not consider any AGN components. To avoid possible systematics due to these factors, we perform the SED fitting on the DustPedia dataset. To maintain consistency, we adopt the same fitting parameters and bandpasses used for the S82 dataset. We exclude the photometric data contaminated by artifacts or nearby sources, as well as data obtained from the imaging data that partially covers the target source (Clark et al., 2018). To secure accurate stellar mass measurements from the SED fitting, samples lacking flux measurements in the optical band ($gri$) are further discarded. The systematic effects due to the discrepancy in the fitting parameters between the original DustPedia study (Nersesian et al., 2019) and this study are further discussed in Appendix A.

We discard the targets with a reduced $\chi^{2}$ value exceeding 2 to ensure reliable assessments of the stellar and dust properties. Furthermore, following the recipe from Li et al. (2023), the AGNs defined as $f_{\rm AGN}>0.1$ were excluded, as the flux contribution from the AGN can introduce additional uncertainties. This results in the final samples containing 435 and 2418 galaxies from the DustPedia and S82, respectively. Throughout this study, we utilize the synthetic SEDs obtained from the best fit for magnitude estimations across various filters and wavelengths. The goodness-of-fit, particularly in NIR regions, is examined in Appendix B. The stellar masses derived from the SED fitting agree with those of GALEX-SDSS-WISE Legacy Catalog 2 (GSWLC-2; Salim et al., 2018) with $\Delta M_{*}=0.01\pm 0.11$ dex (Li et al., 2023). However, a moderate difference is observed in the star formation rate when compared to GSWLC-2 ($\Delta\log\ {\rm SFR}=-0.12\pm 0.56$ dex). As described in Li et al. (2023), this discrepancy arises from the differences in the methods of photometric measurements, particularly in the mid- and far-infrared bands.

Figure 1 illustrates the distributions of stellar and dust masses determined from the SED fitting process, demonstrating the complementary characteristics of the two subsamples in the two-parameter space. Furthermore, our sample covers a wide range of SFR at a given stellar mass, from the star-forming main sequence to relatively quiescent galaxies, which is crucial to investigate the variations in the $M_{*}/L$ ratio with various properties of galaxies (Fig. 2).

Conventionally, stellar-mass estimations in the NIR band ($3-5\ \mu$m) have been based on broadband photometry obtained from the Spitzer Space Telescope and the WISE mission. As the first step of our analysis, we compute $M_{*}/L$ at commonly used broadband filters in the NIR (e.g., IRAC1 and W1). Figure 3 shows $M_{*}/L$ at IRAC1 and W1 using the stellar and total luminosities. Overall, we can estimate the stellar mass with an accuracy of $0.14-0.15$ dex solely with the NIR data, which is consistent with the previous studies (e.g. Meidt et al., 2014; Kettlety et al., 2018; Jarrett et al., 2023).¹¹1It is worthwhile to note that this accuracy is calculated relative to the stellar mass derived from the SED fitting performed using the parameterized SFH, potentially introducing a zero-point offset of up to $0.3-0.4$ dex. This bias, however, is not considered in this study (e.g. Conroy, 2013; Lower et al., 2020). However, as the dust continuum and 3.3 $\mu$m PAH emission can contribute to the total luminosity in those filters, the $M_{*}/L$ ratios calculated using the total luminosity exhibit a substantially larger scatter compared to those calculated using the stellar luminosity.

Meanwhile, considering only the stellar light, the scatter in $M_{*}/L$ can be significantly reduced down to $\sim 0.10$ dex. While the average $M_{*}/L$ ratio at IRAC 1 [$M_{*}/L\sim 0.58$ or $\log\ (M_{*}/L)\sim-0.24$] is in good agreement with the prediction for the old stellar population ($M_{*}/L\sim 0.60$; Meidt et al. 2014; Querejeta et al. 2015), that at the W1 band [$M_{*}/L\sim 0.54$ ($\log\ (M_{*}/L)\sim-0.27$) and $\sim 0.59$ ($\log\ (M_{*}/L)\sim-0.23$] from the total and stellar luminosity, respectively] is significantly larger than the empirically derived value for the nearby galaxies ($M_{*}/L\sim 0.35$; Jarrett et al., 2023). Note that $M_{*}/L$ in Jarrett et al. (2023) is derived from the W1 total luminosity. With the EAGLE simulations, Norris et al. (2016) examined the M/L distribution of simulated early-type/quiescent galaxies and found $M_{*}/L\sim 0.85$, which is in good agreement with a peak in the higher $M_{*}/L$ found in the distribution of $M_{*}/L$ (Fig. 3). While the origin of this discrepancy remains rather unclear, it can be partially attributed to differences in flux measurements across the chosen datasets. Specifically, the dedicated flux measurements of the WISE dataset obtained from Li et al. (2023) are systematically brighter by $0.2-0.5$ mag compared to the ALLWISE measurements based on profile-fitting results, which was used to estimate the $M_{*}/L$ values in Jarrett et al. (2023). Interestingly, in our study, the $M_{*}/L$ values calculated for W1 and IRAC1 are consistent with each other, which supports the reliability of our estimates. It is also worthwhile noting that, unlike other studies, Jarrett et al. (2023) used the observed total luminosity, instead of the stellar luminosity, to calculate the $M_{*}/L$ values, naturally resulting in a low $M_{*}/L$.

As the NIR spectral data will be obtained through future space missions, such as SPHEREx, investigating the $M_{*}/L$ in these spectral elements is worthwhile. To this end, we employ 68 spectral components with a spectral resolution of 40 (i.e., $\lambda/\Delta\lambda\approx 40$, where $\Delta\lambda$ is the bandwidth of each spectral component), covering a wavelength range of 0.75$-$5.0 $\mu$m. For simplicity, the transmission curve of each spectral component is modeled using a Gaussian profile. Although this assumption does not exactly match the spectral channels of the SPHEREx, it is adequate for examining the overall trend of the $M_{*}/L$ ratio as a function of the wavelength.

As expected, the $M_{*}/L$ ratio is inversely proportional to the wavelength because the $M_{*}/L$ ratios in the shorter wavelengths are more sensitive to the young stellar population compared to those at the longer wavelength. However, the scatter of the $M_{*}/L$ ratio exhibits a U-shaped pattern. The $M_{*}/L$ values calculated using the total luminosity are significantly affected at longer wavelengths due to the existence of the dust continuum and PAH emission (Fig. 4). In particular, the excess in the scatter is remarkable around 3.3 $\mu$m, highlighting the need for the careful removal of the PAH emission to facilitate accurate stellar mass estimations. This effect can be rigorously quantified using the spectral dataset provided by the SPHEREx mission, particularly for star-forming galaxies (e.g. Xie et al., 2018a, b; Zhang et al., 2021).

In addition, the $M_{*}/L$ values in the DustPedia sample are systematically smaller than those in the S82 sample. The $M_{*}/L$ values are known to be strongly dependent on the specific SFR (e.g., Portinari et al., 2004). The DustPedia sample tends to have a larger specific SFR than the S82 sample (see §4.2), which may be the main cause of smaller $M_{*}/L$ values for the DustPedia sample.

The scatter in the $M_{*}/L$ ratios calculated from the stellar luminosity is significantly reduced (Fig. 5). At wavelengths longer than 1$\mu$m, the overall scatter ranges between 0.1 and 0.11 dex, demonstrating a weak dependence on the wavelength. This indicates that the NIR spectrum is valuable for constraining the stellar mass almost regardless of the wavelength if the dust component can be adequately removed. The $M_{*}/L$ values and their scatter for both subsamples are summarized in Table 1.

The NIR colors (IRAC1-IRAC2 or W1-W2) have been utilized to approximate the light contribution from the dust, which needs to be adequately subtracted to estimate the stellar mass precisely. The basic assumption of this procedure is that the dust continuum is intrinsically redder than the stellar continuum. According to this assumption, Querejeta et al. (2015) and Jarrett et al. (2023) demonstrated that $M_{*}/L$ is inversely proportional to both the IRAC1-IRAC2 and W1-W2 colors, respectively. Here, we examine the validity of this finding for our dataset. We find two sequences in the plane of $M_{*}/L$ and W1-W2 color for our sample (Fig. 6). In one sequence, consistent with the findings of previous studies, $M_{*}/L$ is mildly anti-correlated with the W1-W2 and IRAC1-IRAC2 color. Although this correlation broadly agrees with the general trend in earlier studies (e.g., Cluver et al., 2014; Jarrett et al., 2023), the observed relation between the two quantities deviates from those of previous studies. We attribute this to the fact that the photometric measurements of W1 and W2 adopted in this study are systematically brighter than the ALLWISE measurements (Li et al. 2023).

However, in another sequence, $M_{*}/L$ is almost independent of the W1$-$W2 color, with a substantially large scatter. To comprehend this finding, we compare the NIR color with the ratio of the flux from the dust components to the total flux ($f_{\rm dust}$) at the W1 band. Figure 7 shows that the dust fraction at the W1 band is only weakly correlated with the W1$-$W2 color. The hot and warm dust ($T\sim 300-500$ K) from the star-forming regions can naturally increase the W1$-$W2 color (e.g., Xie et al., 2018a). Conversely, the prominent 3.3 $\mu$m PAH emission can boost the flux only at the W1 band, making the W1$-$W2 color bluer. Because the strength of the 3.3 $\mu$m PAH emission is proportional to the SFR but with a significant amount of the scatter, the correlation between the W1$-$W2 color and the dust fraction can be excessively intricate due to the competition between the dust continuum and PAH emission. This complexity cannot be quantified without a spectroscopic dataset (e.g., Yamada et al., 2013; Shim et al., 2023). A similar trend (i.e., a non-linear correlation between the dust fraction and the NIR color) is also shown in the IRAC1$-$IRAC2 color.

As SPHEREx will provide the spectral information, the spectral color can serve as a better tracer for the dust contribution. Here, we adopt a color index between 3 $\mu$m and 4 $\mu$m ($m_{\rm 3\mu m}-m_{\rm 4\mu m}$) to avoid the strong 3.3 $\mu$m emission and compare it with the flux ratio of the dust at 3 $\mu$m ($f_{\rm dust,3\mu m}$). We find that $f_{\rm dust,3\mu m}$ is strongly correlated with $m_{\rm 3\mu m}-m_{\rm 4\mu m}$ expressed as $f_{\rm dust,3\mu m}=-0.05+0.41(m_{\rm 3\mu m}-m_{\rm 4\mu m})$ (Fig. 8), indicating that the dust contribution can be approximately estimated with the spectral color index in SPHEREx dataset.

Given that the $M_{*}/L$ ratios can be highly affected by the stellar properties of the galaxies, adopting a second parameter is often essential for robustly estimating the stellar mass. Previous studies have extensively adopted the color term as the second parameter (e.g., Bell et al., 2003; Into & Portinari, 2013). To identify an alternative in the spectral analysis, we investigate the parameter predominantly deriving the scatter in $M_{*}/L$. Initially, we find that $M_{*}/L$ is weakly correlated with the stellar mass itself, likely attributed to the dependence of the SFH on stellar mass. To account for this dependence, we perform a linear regression between the stellar mass and $M_{*}/L$ at each wavelength as follows: $\log\ (M_{*}/L)=a+b\log\ (M_{*}/M_{\odot})$ (see Figure 9 for an example). The fitting results as a function of the wavelength are summarized in Table 2. With this treatment, the scatter is mildly reduced by $\sim 0.01-0.02$ dex. Thus, overall, with the NIR spectrum, the stellar mass can be estimated with a typical uncertainty of $\sim 0.09-0.1$ dex (Fig. 10).

Furthermore, the $M_{*}/L$ ratio is tightly correlated with the SFR divided by the stellar mass, namely the specific star formation rate (sSFR). Figure 11 shows that this relation at $\sim 1.6\ \mu$m is well fitted with a smoothly broken power law, which is defined as:

where $a$ is the normalization factor [i.e., $\log\ (M_{*}/L)$ at ${\rm sSFR}\ll\delta$], $\delta$ is the pivot sSFR, $\alpha$ is the power-law index at ${\rm sSFR}>\delta$, and $\beta$ is the smoothness parameter, which determines the smoothness of the power-law slope change around the pivot sSFR. Note that the unit of sSFR is Gyr^-1. Fitting results for various wavelengths are summarized in Table 3. This trend is rather expected from the adopted SFH with two components (old and young stellar populations) because the sSFR approximately traces the mass ratio of the old stellar population to the young stellar population. Considering this dependency, the scatter in $M_{*}/L$ dramatically decreases to $\sim 0.02$ dex (Fig. 12). This result clearly demonstrates that more precise stellar mass estimations can be achieved using the known SFR. Note that the SFR used in this analysis is an instantaneous SFR provided by the CIGALE based on the star formation history. As the commonly used SFR indicators (e.g., hydrogen recombination lines) trace the average SFR over the last $\sim 10$ Myr, the observed SFRs may differ from those provided by CIGALE to some degree, which can introduce additional uncertainty in the $M_{*}/L$ estimation (e.g., Byun et al., 2021).

Based on our calculations of $M_{*}/L$ under various conditions, we attempt to identify an optimal method to estimate the stellar mass based on the NIR spectral data of nearby galaxies. If the contribution of the dust emission can be adequately removed from the observed spectrum, the scatter in the $M_{*}/L$ ratio becomes almost independent of the wavelength. However, it is challenging to fit the PAH emissions and the dust continuum robustly, even using the NIR/MIR spectroscopic data (e.g. Xie et al., 2018a; Zhang & Ho, 2023). Therefore, if available, NIR data around a wavelength of $1.6\ \mu$m are preferable, as will also be presented in Lee, J. H. et al. (submitted). This is because it is relatively free from dust contamination unless the AGN component is prominent, and the resulting scatter in $M_{*}/L$ is smaller or comparable to that achieved using NIR data.

As outlined in the previous section (§4.2), the scatter in $M_{*}/L$ can be slightly reduced if the dependence of $M_{*}/L$ on stellar mass is considered. To this end, the stellar mass can be estimated from $M_{*}/L$ using the monochromatic luminosity and subsequently subjected to iterative fine-tuning. Specifically, the initial stellar mass is estimated from a constant $M_{*}/L$, following which the correlation between the $M_{*}/L$ ratio and stellar mass can be corrected using the initial stellar mass measurement to achieve a more precise stellar-mass estimation. However, this iterative process can reduce the uncertainty by only $0.01-0.02$ dex. Conversely, if the SFR can be estimated based on the NIR spectral data, for example, from the fluxes of hydrogen recombination lines (e.g., Kennicutt & Evans, 2012) or 3.3 $\mu$m PAH emissions (e.g. Kim et al., 2012; Lai et al., 2020; Belfiore et al., 2023), the scatter in the stellar mass estimation can be dramatically reduced based on the strong correlation between $M_{*}/L$ and sSFR. It is worthwhile to note that 3.3 $\mu$m PAH-based SFR can be more uncertain compared to other SFR indicators (e.g., hydrogen lines, UV, and IR) as its brightness relative to IR luminosity is correlated with physical parameters of host galaxies, such as metallicity (e.g. Shim et al., 2023; Whitcomb et al., 2024). However, this correction must again follow an iterative procedure, as a prior stellar mass estimate is essential to compute the sSFR. An alternative argument is that SFR, rather than sSFR, could be used to trace $M_{*}/L$. However, we find that the scatter in the correlation between SFR and $M_{*}/L$ is significantly larger than that for sSFR. Therefore, sSFR provides a more reliable estimate of $M_{*}/L$. In conclusion, the NIR spectrophotometry provided by the SPHEREx can improve the accuracy of stellar mass estimations.