Resolved Stellar Mass Estimation of Nearby Late-type Galaxies for the SPHEREx Era: Dependence on Stellar Population Synthesis Models

For the spectral analysis in this study, we obtained the E-MILES templates from the MILES cloud storage³³3https://cloud.iac.es/index.php/s/aYECNyEQfqgYwt4. The E-MILES model provides stellar spectral templates covering wavelengths ranging from $0.16~{}{\rm\mu m}$ to $5~{}{\rm\mu m}$, which is an extended version of the MILES stellar library (Sánchez-Blázquez et al., 2006). These templates offer a high spectral resolution of $\Delta\lambda_{\rm FWHM}\sim 2.5~{}{\rm\AA}$ in the optical wavelength range up to $0.9~{}{\rm\mu m}$, but the resolution decreases to $\Delta\lambda_{\rm FWHM}\sim 15-20~{}{\rm\AA}$ in the Spitzer/IRAC $3.6~{}{\rm\mu m}$ and $4.5~{}{\rm\mu m}$ bands (see Figure 8 in Vazdekis et al., 2016). We selected the E-MILES templates based on the Padova 2000 isochrones (Girardi et al., 2000) and a scaled-solar model for metallicity, assuming ${\rm[Fe/H]=[Z/Z_{\odot}]}$. For our stellar population analysis, we utilized spectral templates for 12 age bins of [0.0708, 0.1122, 0.1778, 0.2818, 0.5012, 0.7943, 1.2589, 1.9953, 3.1623, 5.6234, 8.9125, 14.1254] in Gyr and four metallicity bins of ${\rm[Z/Z_{\odot}]}=[-0.71,-0.40,0.00,0.22]$.

However, the E-MILES models lack spectral templates for stellar ages younger than $63~{}{\rm Myr}$, which can introduce systematic uncertainties when analyzing young stellar populations, as discussed in Section 6.3.2 of Emsellem et al. (2022). To address this limitation, we supplemented our analysis with “superyoung” models, as done in the previous PHANGS-MUSE study on stellar population properties by Pessa et al. (2023). These superyoung E-MILES models differ from their older versions in the sense that they employ the Padova 1994 isochrones (Girardi et al., 1996) and provide a higher metallicity set with ${\rm[Z/Z_{\odot}]=0.41}$ instead of ${\rm[Z/Z_{\odot}]=0.22}$. In addition, these templates do not include IR spectra, covering only wavelengths up to $0.9~{}{\rm\mu m}$. This limitation could potentially affect the accuracy of mass-to-light ratios at IR wavelengths for regions with stellar ages younger than $10^{8}~{}{\rm yr}$. Following Pessa et al. (2023), we included the additional superyoung models with five age bins of [0.0063, 0.0100, 0.0158, 0.0251, 0.0398] in Gyr and four metallicity bins of ${\rm[Z/Z_{\odot}]}=[-0.71,-0.40,0.00,0.41]$.

BC03 is one of the most widely used SPS models for predicting the distribution of stellar populations. In this study, we utilized the updated 2016 version of the BC03 templates, computed using the GALAXEV code⁴⁴4https://www.bruzual.org/bc03. The BC03 models cover a wide wavelength range from $91~{}{\rm\AA}$ to $160~{}{\rm\mu m}$. They combine a high-resolution spectral library ($\Delta\lambda_{\rm FWHM}\sim 3~{}{\rm\AA}$) from STELIB (Le Borgne et al., 2003) for the $3200-9500~{}{\rm\AA}$ range, alongside lower-resolution libraries ($\Delta\lambda_{\rm FWHM}\sim 20~{}{\rm\AA}$) from BaSeL 3.1 (Westera et al., 2002) and Pickles (Pickles, 1998) for the remaining wavelength range. These models are based on the Padova 1994 isochrones, rather than the Padova 2000 isochrones, as detailed in Bruzual & Charlot (2003). BC03 provides spectral templates for a wide range of stellar ages from $0.1~{}{\rm Myr}$ to $20~{}{\rm Gyr}$. For this study, we adopted templates for 17 age bins of [0.0063, 0.01, 0.0158, 0.0251, 0.04, 0.0719, 0.1139, 0.1805, 0.2861, 0.5088, 0.8064, 1.2781, 2.0, 3.25, 5.5, 9.0, 14.25] in Gyr and four metallicity bins of ${\rm[Z/Z_{\odot}]}=[-0.70,-0.40,0.00,0.40]$.

CB19⁵⁵5https://www.bruzual.org/CB19/CB19_chabrier/Mu100 is the latest major revision of BC03, incorporating the PARSEC stellar evolutionary tracks (Bressan et al., 2012) instead of the older Padova isochrones. These updated PARSEC tracks incorporate newly computed models for the evolutionary physics of TP-AGB stars (Marigo et al., 2013) and stellar winds from hot massive stars, such as OB-type and Wolf–Rayet stars (Chen et al., 2015). The CB19 models provide a broad wavelength range from $5.6~{}{\rm\AA}$ to $36,000~{}{\rm\mu m}$ by integrating various spectral libraries: MILES for wavelengths of $\lambda\sim 3540-7350~{}{\rm\AA}$, STELIB for wavelengths of $\lambda=7351-8750~{}{\rm\AA}$, and BaSeL 3.1 combined with several TP-AGB models for wavelengths longer than $8750~{}{\rm\AA}$ (see Appendix A in Sánchez et al., 2022, and references therein). In this study, we used the CB19 templates including 17 age bins of [0.0063, 0.01, 0.016, 0.025, 0.040, 0.070, 0.110, 0.180, 0.275, 0.5, 0.8, 1.3, 2.0, 3.25, 5.5, 9.0, 14.0] in Gyr and four metallicity bins of ${\rm[Z/Z_{\odot}]}=[-0.70,-0.40,0.00,0.30]$.

The FSPS models, developed by Conroy et al. (2009), offer great flexibility in terms of IMFs, isochrones, spectral libraries, and physical treatments for stellar phenomena such as the horizontal branch (HB), blue stragglers, and TP-AGB stars. For this study, we compiled FSPS v3.2⁶⁶6https://github.com/cconroy20/fsps and generated the default FSPS models using the python-fsps⁷⁷7https://github.com/dfm/python-fsps package. These templates are based on the spectral libraries of MILES ($\lambda\sim 3600-7400~{}{\rm\AA}$) and BaSeL 3.1 ($\lambda=91~{}{\rm\AA}-160~{}{\rm\mu m}$, excluding the MILES coverage), the MIST isochrones (Choi et al., 2016), and the interstellar dust emission model from Draine & Li (2007). We used the FSPS models for 14 age bins of [0.0063, 0.01, 0.0158, 0.0251, 0.0398, 0.0708, 0.1122, 0.1778, 0.2818, 0.5012, 0.7943, 1.2589, 1.9953, 3.1623, 5.6234, 8.9125, 14.1254] in Gyr and four metallicity bins of ${\rm[Z/Z_{\odot}]}=[-0.70,-0.40,0.00,0.40]$.

Full-spectrum fitting is a powerful technique for quantifying stellar population parameters, such as age, metallicity, and mass, and for reconstructing the detailed SFH of stellar systems. In this study, we employed the Penalized PiXel-Fitting (pPXF; Cappellari & Emsellem, 2004; Cappellari, 2017, 2023) program⁸⁸8https://pypi.org/project/ppxf/9.1.1 for full-spectrum fitting, which is known for its superior performance in recovering stellar population parameters and computational efficiency (Ge et al., 2018; Woo et al., 2024). For our implementation of pPXF, we applied the four different SPS models described in Section 3 to the optical IFU datacubes rebinned to the SPHEREx pixel size, as described in Section 2.

Prior to running pPXF, we manually masked all emission lines, including Balmer lines and strong forbidden line doublets (e.g., [O III]$\lambda\lambda 4959,5007$, [N II]$\lambda\lambda 6548,6584$, and [S II]$\lambda\lambda 6717,6731$), to minimize contamination from nebular emission in the stellar light components. Following the methods used in the previous PHANGS-MUSE works (Emsellem et al., 2022; Pessa et al., 2023), we implemented the pPXF tasks in two steps, i.e., deriving the stellar kinematics and estimating the stellar populations with the fixed kinematic parameters. For both steps, we used a spectral range of $4850-7000~{}{\rm\AA}$ to avoid unstable sky subtraction at longer wavelengths, as noted by Emsellem et al. (2022). These procedures ensured consistency with previous PHANGS-MUSE results.

To examine the stellar kinematics, we measured four parameters: the radial velocity along the line of sight ($v_{\rm rad}$), stellar velocity dispersion ($\sigma_{v}$), and third- and fourth-order Gauss–Hermite coefficients ($h_{3}$ and $h_{4}$, respectively).

The kinematics was determined through three iterations of the pPXF process to obtain stable solutions. The first iteration provided initial guesses for the kinematics and constructed noise models by comparing the best-fit templates with the observed data. The second iteration refined these initial guesses and noise models, and the third iteration determined the final kinematic parameters. Throughout these processes, we used additive Legendre polynomials with a degree of 12 (degree=12), no multiplicative Legendre polynomials (mdegree=0), and a regularization error of 0.2 (regul=5) , consistent with earlier PHANGS-MUSE works. Following Emsellem et al. (2022), we used only the E-MILES library with the BaSTI isochrones (Pietrinferni et al., 2004) to explore the kinematics, which differs from the template sets used for stellar population analysis.

To investigate the stellar populations, we applied the four different SPS models in the pPXF tasks, using the fixed kinematic parameters derived from the previous step. Similar to the kinematic analysis, we conducted four iterations of the pPXF tasks. The first and second iterations were focused on generating and improving the noise models for the input spectra. The third iteration was aimed at deriving internal extinction magnitudes ($A_{V}$) using the reddening keyword, applying the Calzetti et al. (2000) extinction law. The spectra were then corrected for internal extinction using the derived $A_{V}$ values. In the final iteration, we used the extinction-corrected spectra to derive the weighting coefficients of the input spectral templates, which were then used to determine the stellar population parameters. Throughout these processes, we used no additive polynomials (degree=-1), multiplicative polynomials with a degree of 12 (mdegree=12), and no regularization (regul=0), except during the third iteration for internal extinction. These configurations were consistently set according to Pessa et al. (2023). In the third iteration, we used no multiplicative polynomials (mdegree=0) and activated the reddening option (similar to Choice-1 in Lee et al., 2024).

Figure 2 illustrates example results of the full-spectrum fitting process for several rebinned spaxels in the central regions of the PHANGS-MUSE galaxies. The zoomed-in panels show the best-fit model spectra (blue lines) compared with the observed spectra (orange lines) within the given wavelength ranges. As all SPS models used in this study cover NIR wavelengths up to $5~{}{\rm\mu m}$, the best-fit model spectra can be extended to these wavelengths by applying the output stellar population parameters and the weighting coefficients of the input spectral templates. The zoomed-out panels display the best-fit results across the full wavelength ranges of the SPS templates, from UV to NIR. For these panels, the photometric data of 2MASS and Spitzer/IRAC were obtained from the pixel values in the images reprojected to the SPHEREx pixel size. These wavelength extensions of the best-fit results were used to estimate the resolved stellar mass and mass-to-light ratios at NIR wavelengths, as will be mentioned in Section 4.2. For our analysis, we excluded the SPHEREx-binned spaxels with a signal-to-noise ratio below 155 (equivalent to a mean signal-to-noise ratio of 5 for the original unbinned pixel scale) and unreliable kinematics measurements with $v_{\rm rad}<500~{}{\rm km~{}s^{-1}}$ or $\sigma_{v}>200~{}{\rm km~{}s^{-1}}$.

The pPXF routines provide luminosity-weighted coefficients ($w_{i}^{L}$) for the $i$-th SSP template in the input SPS models for all spaxels of the SPHEREx-rebinned PHANGS-MUSE datacubes. To calculate stellar mass-to-light ratios from the best-fit pPXF results, we converted the output luminosity weights to mass weights ($w_{i}^{M}$) through division by the mean flux of each SSP template within a normalizing wavelength range of $\lambda=5070-5950~{}{\rm\AA}$, corresponding to the typical FWHM range of the $V$-band. Subsequently, we derived the mass-to-light ratio for a specific band ($M_{\ast}/L_{\lambda}$) from the best-fit pPXF results of each spaxel using the following formula:

where $m_{\ast,i}$ represents the stellar mass of the $i$-th SSP spectral template in solar mass ($M_{\odot}$), and $L_{\lambda,i}$ is the specific luminosity at a given band for the $i$-th SSP template in specific solar luminosity ($L_{\lambda,\odot}$). The values for $m_{\ast,i}$ are predetermined for all SSP templates by each SPS model, including the masses of living stars and stellar remnants such as white dwarfs, neutron stars, and black holes. The luminosity at a specific band for the best-fit template was computed using the PyPhot package.

Using Equation 1, we computed the stellar mass-to-light ratios for all spaxels, ranging from the SDSS optical bands to the Spitzer/IRAC NIR bands. The stellar mass of each spaxel was determined by multiplying the specific luminosity in the SDSS $r$ band by the corresponding mass-to-light ratio derived from its best-fit SED template. Through these procedures, we created resolved stellar mass maps for the 19 PHANGS-MUSE galaxies, with an example shown in Figure 3. These maps serve as the fiducial stellar masses used throughout this study. We did not account for the inclination effects of disk galaxies in this study.

Since the full-spectrum fitting in this study was performed on SPHEREx-binned IFU data, it is important to assess the potential influence of this binning on the pPXF results. To quantify the effect of SPHEREx binning, we compared the stellar kinematics and stellar population parameters obtained in this study with those from previous PHANGS-MUSE analyses (Emsellem et al., 2022; Pessa et al., 2023), which were based on Voronoi-binned datacubes. Although stellar kinematics cannot be measured with SPHEREx data due to its low spectral resolution, this kinematics comparison is necessary for testing the reliability of the stellar population parameters because stellar kinematics can affect the stellar population analysis during the pPXF fitting process. For the kinematics comparison, we used the value-added maps from PHANGS-MUSE data release v1.0, which provide stellar radial velocity and velocity dispersion maps derived using the pPXF configurations specified by Emsellem et al. (2022). For stellar populations, we digitized the data of luminosity-weighted stellar ages and metallicities from Figures 2 and 7 of Pessa et al. (2023), as these data were not publicly available. To maintain consistency with previous studies, all comparative analyses described in this section were based on the E-MILES models.

Figure 4 presents comparisons of stellar radial velocity (left) and velocity dispersion (right) for 19 PHANGS-MUSE galaxies as a function of galactocentric distance, normalized by effective radius (Table 4 in Leroy et al., 2021). The radial velocities derived from the SPHEREx-binned data in this study exhibit a slight underestimation compared with the Voronoi-binned data across all radial bins, with a median discrepancy of $\sim-5~{}{\rm km~{}s^{-1}}$. Given that the systematic radial velocities for the PHANGS-MUSE sample range from $650~{}{\rm km~{}s^{-1}}$ to $2400~{}{\rm km~{}s^{-1}}$, this discrepancy is considered negligible. In terms of velocity dispersion, the results from this study show a slight overestimation of $\Delta\sigma_{v}\lesssim 5~{}{\rm km~{}s^{-1}}$ within the regions up to $r/R_{\rm e}\lesssim 2$, consistent with those from the previous study. However, velocity dispersions from the SPHEREx-binned data tend to be higher in the outer regions than those from original PHANGS-MUSE data, with a discrepancy up to $\Delta\sigma_{v}\sim 10~{}{\rm km~{}s^{-1}}$ at $r/R_{\rm e}>3$. This overestimation may be associated with the smearing effect due to the SPHEREx-binning process. This discrepancy is non-negligible, considering the median velocity dispersion of $49~{}{\rm km~{}s^{-1}}$ and its uncertainty of $17~{}{\rm km~{}s^{-1}}$ from the original datacubes. Nonetheless, the effect of SPHEREx binning on velocity dispersion appears to be limited to a small number of pixels in the outer regions ($r/R_{\rm e}>3$).

Figure 5 presents comparisons of luminosity-weighted mean stellar ages (left) and metallicities (right) for all galaxies. For this figure, digitized data from Pessa et al. (2023) were interpolated to match our radial bins. The stellar age estimates derived in this work roughly agree with those presented in Pessa et al. (2023) within a scatter of $\sim 0.2~{}{\rm dex}$. This large scatter is expected because the data from Pessa et al. (2023) represent luminosity-weighted mean values in their radial bins, which are not ideal for pixel-to-pixel comparisons. Notably, in the regions within $r/R_{\rm e}<2$, the stellar ages from this study tend to be systematically overestimated by $0.1-0.2~{}{\rm dex}$ compared with those from Pessa et al. (2023), while this systematic offset appears to diminish in the outer regions ($r/R_{\rm e}>2$) despite a larger scatter. These findings suggest that the binning process marginally affects the accuracy of stellar ages in the inner part of galaxies, likely due to the mixing of flux from older spaxels with higher surface brightness and younger spaxels with relatively lower surface brightness, which can bias the age estimates toward older values. In contrast, the radial distribution of metallicity exhibits a good match with that from Pessa et al. (2023) across all radial bins, with a minor scatter of $\sim 0.1~{}{\rm dex}$, revealing that binning has little impact on metallicity measurements. This may be associated with the flatter gradient of stellar metallicity compared to that of stellar ages in most PHANGS-MUSE galaxies, as also noted in Pessa et al. (2023).

In conclusion, the SPHEREx binning procedure in this study has negligible or minimal effects on stellar kinematics and metallicities. This implies that binning does not significantly affect the overall results from full-spectrum fitting, in particular for the resolved stellar mass estimation. However, in the case of stellar ages, we find a tendency for overestimation by $\sim 0.1-0.2~{}{\rm dex}$ with a large scatter in the regions within $r/R_{\rm e}<2$, suggesting a marginal effect of binning on the age distribution.

The choice of SPS models can influence the derived stellar populations obtained through full-spectrum fitting (Ge et al., 2019). In this section, we compare the stellar properties derived from four SPS models to clarify the systematic effects of model selection. Figure 6 displays the distributions of (luminosity-weighted) stellar ages, metallicities (${\rm[Z/Z_{\odot}]}$), and internal extinction magnitudes ($A_{V}$) for all SPHEREx-binned pixels from 19 PHANGS-MUSE galaxies. In this figure, different horizontal panels correspond to the four SPS models. The associated statistics are listed in Table 1.

Across all SPS models, the mean stellar ages exhibit broad distributions, ranging from ${\rm log~{}Age~{}(yr)}\sim 7.9-9.7$, indicating the presence of a considerable amount of young stellar populations in the late-type galaxy sample of PHANGS-MUSE. Stellar metallicity and internal extinction are correlated with stellar age: metallicity increases while dust extinction decreases as stellar populations become old. These trends align with conventional concepts of the relationships between stellar populations and galaxy evolution. Stellar ages and metallicities are known to positively correlate with stellar masses, as the SFH and chemical enrichment are closely related to the stellar mass growth, resulting in mass–age and mass–metallicity relations (Tremonti et al., 2004; Gallazzi et al., 2005; González Delgado et al., 2015). The negative correlation between internal dust extinction and age is also expected, as dust is more abundant in younger star-forming regions than in older stellar regions.

Remarkable discrepancies in stellar population parameters appear depending on the adopted SPS models. First, the E-MILES model with the superyoung templates tend to predict older stellar ages and lower metallicities than the other SPS models by $0.1-0.2~{}{\rm dex}$. This indicates that results from previous PHANGS-MUSE studies (Emsellem et al., 2022; Pessa et al., 2023) might vary slightly if SPS models other than the E-MILES templates are used. Second, the metallicities derived using the BC03 model exhibit a much wider distribution than those from the other SPS models, with the 16th and 84th percentile interval being approximately $0.42~{}{\rm dex}$ for BC03, which is $0.1-0.2~{}{\rm dex}$ broader than that for the other models. This wide distribution could be attributed to the tendency of the BC03 model to predict a higher fraction of pixels with super-solar metallicity. The remaining two SPS models, CB19 and FSPS, yield similar results for age and metallicity.

The $V$-band extinction magnitudes exhibit consistent distributions across the SPS models. For the BC03 model, however, the extinction magnitudes were fixed to match those from the CB19 model, as represented in Figure 6 and Table 1. Before this adjustment, we observed that BC03 significantly underestimated the $A_{V}$ values by $\sim 0.3~{}{\rm mag}$ compared with other SPS models, with around $30\%$ of all spaxels showing $A_{V}=0.0$, which appeared unphysical. As this discrepancy of $A_{V}$ can introduce a potential systematic bias in the other stellar population parameters, we manually adjusted the $A_{V}$ for BC03 throughout this study.

In this section, we explore the distributions of the NIR color index, specifically Spitzer/IRAC $[3.6]-[4.5]$, predicted by four different SPS models. This analysis is helpful in evaluating the effect of the choice of SPS models on the prediction of stellar properties in the NIR. Since the $[3.6]-[4.5]$ color index is sensitive to the treatment of TP-AGB stars and molecular absorption features in SPS models, this color serves as a useful diagnostic tool for assessing the NIR SEDs of stellar populations. Moreover, the $[3.6]-[4.5]$ color has been widely used to trace old stellar populations and disentangle stellar light from non-stellar contributions, such as dust and polycyclic aromatic hydrocarbon (PAH) emissions, in nearby galaxies. For instance, studies from the ${\rm S^{4}G}$ survey demonstrated that the $[3.6]-[4.5]$ color is effective to derive “uncontaminated” stellar masses by distinguishing stellar light from total IRAC $3.6~{}{\rm\mu m}$ fluxes, with a typical range of $[3.6]-[4.5]$ for old stars between $-0.2$ and $0.0$ (Meidt et al., 2014; Querejeta et al., 2015). In this context, we examine the systematic differences in $[3.6]-[4.5]$ color distributions across SPS models to evaluate their NIR SED predictions.

To highlight these variations, Figure 7 shows the distribution of $[3.6]-[4.5]$ colors as a function of stellar age for each SPS model. These colors were derived from pPXF fitting by applying IRAC bandpasses to the best-fit stellar templates, considering only stellar emissions. Among the SPS models, E-MILES exhibits a particularly narrow color distribution between $-0.15$ and $0.0$. This color range is consistent with the stellar IRAC colors provided in the ${\rm S^{4}G}$ survey results (see Table 1 in Querejeta et al. (2015)). In contrast, the other SPS models show substantially redder and broader color distributions, ranging from $-0.05$ to $\sim 0.2$. The age dependence of the $[3.6]-[4.5]$ colors also varies with the SPS models. While E-MILES and FSPS show nearly no correlation between color and stellar age, BC03 and CB19 exhibit slightly age-dependent trends, with colors becoming the reddest around ${\rm log~{}Age}=8.5-9.0$. These discrepancies may be attributed to the variations in the treatment of intermediate-age stellar emission at NIR wavelengths. Although E-MILES systematically overestimates stellar ages by $0.1-0.2~{}{\rm dex}$ as presented in Section 5.1, this age difference cannot fully explain the substantial color discrepancies between E-MILES and the other SPS models, as the $[3.6]-[4.5]$ color in E-MILES consistently falls between $-0.15$ and $0.0$ across all age ranges. Thus, these NIR color distributions indicate intrinsic differences in the NIR SEDs among the SPS models.

To further examine these differences, we illustrate example SEDs predicted by the four different SPS models with a specific age and metallicity in Figure 8. The SEDs are modeled with ${\rm log~{}Age~{}(yr)}=8.9$ and solar metallicity using each SPS model, which approximately corresponds to the median values of the derived stellar population parameters (see Table 1). This figure clearly reveals that the SED deviations between SPS models are larger in the NIR compared to optical wavelengths ($<0.75~{}{\rm\mu m}$). In particular, the CO band at $4.2-4.5~{}{\rm\mu m}$, originating from late-type ($T_{\rm eff}\lesssim 6000~{}{\rm K}$) giant stars (see Figure 7 in Röck et al., 2015), is more prominent in the E-MILES template compared with the other SPS templates. This results in a bluer $[3.6]-[4.5]$ color owing to the reduction of the IRAC $4.5~{}{\rm\mu m}$ fluxes (Peletier et al., 2012). In contrast, in the other SPS models, these CO absorption lines are either absent (BC03 and FSPS) or significantly weaker (CB19). The $[3.6]-[4.5]$ colors from E-MILES align best with observed global $[3.6]-[4.5]$ colors of $-0.2\lesssim[3.6]-[4.5]\lesssim 0$ (Pahre et al., 2004; Peletier et al., 2012), reflecting the detailed treatment of the CO absorption features.

Variations in the treatment of TP-AGB stars, which dominate NIR emission, also contribute to the differences in NIR SEDs among the SPS models. BC03, in particular, is known to exhibit considerable differences in TP-AGB emissions with other SPS models (Maraston et al., 2006; Conroy et al., 2009; Röck et al., 2016). In addition, CB19 shows a more gradual SED slope at wavelengths longer than $4~{}{\rm\mu m}$, compared to other SPS models. These differences in the SED features might contribute to the redder $[3.6]-[4.5]$ colors of BC03 and CB19 shown in Figure 7. Meanwhile, FSPS shows similar SED slope to E-MILES across the NIR wavelengths, but they still differ in the CO absorption features, leading to systematic color variations. Overall, these findings suggest that different SPS models have systematic differences in modeling stellar NIR emissions, which may impact the estimation of stellar mass-to-light ratios in the NIR. Future spectrophotometric data from SPHEREx, covering a spectral range of $0.75-5~{}{\rm\mu m}$ will enable the assessment of the reliability of SPS models and allow for updates to their prescriptions in the NIR range.

In this section, we investigate the dependence of mass-to-light ratios in the NIR on stellar population properties and underlying SPS models. Figure 9 displays the distributions of mass-to-light ratios at IRAC $3.6~{}{\rm\mu m}$ for all spaxels in 19 PHANGS-MUSE galaxies, derived using four SPS models. This figure indicates notable differences in mass-to-light ratios among the SPS models. The mass-to-light ratios derived from the E-MILES and FSPS models are broadly consistent, with FSPS showing a slight underestimation by a median offset of $\sim 0.04~{}{\rm dex}$ ($\sim 9\%$). In contrast, BC03 and CB19 significantly underestimate mass-to-light ratios, with median offsets of $\sim 0.36~{}{\rm dex}$ ($\sim 56\%$) and $\sim 0.26~{}{\rm dex}$ ($\sim 45\%$), respectively, compared to E-MILES. These discrepancies are larger than the typical standard deviations of mass-to-light ratios within individual galaxies (median $\sim 0.15~{}{\rm dex}$), highlighting the critical impact of the underlying SPS models on the derived mass-to-light ratios.

Systematic variations in mass-to-light ratios are also observed across different galaxies, mainly driven by stellar age, as previously noted by Emsellem et al. (2022). In Figure 9, galaxies are organized by their median stellar age, ranging from ${\rm log~{}Age~{}(yr)\sim 8.6}$ in NGC 2835 (leftmost) to ${\rm log~{}Age~{}(yr)\sim 9.4}$ in NGC 3351 (rightmost). A clear positive correlation is evident between NIR mass-to-light ratios and stellar ages across all SPS models: galaxies with older stellar populations consistently exhibit higher mass-to-light ratios. Interestingly, the strength of this correlation varies among the SPS models. CB19 shows the strongest age dependence, with maximal differences in $M_{\ast}/L_{\rm 3.6\mu m}$ reaching $\sim 0.4~{}{\rm dex}$ between the oldest and youngest galaxies, while E-MILES and FSPS show milder correlations with maximal differences of $\sim 0.2~{}{\rm dex}$. These findings imply that mass-to-light ratios in the NIR significantly depend on the choice of SPS model and the stellar ages derived from the model.

Figure 10 further examines the relationship between NIR mass-to-light ratios and stellar ages by plotting all spaxel data derived using different SPS models. This figure is divided vertically into panels based on metallicity bins, defined using the 33rd and 67th percentiles of the metallicity derived from each SPS model. As expected, a strong positive correlation between stellar ages and mass-to-light ratios is evident across all panels, with Spearman’s rank correlation coefficients ($r_{S}$) ranging from $0.4-0.7$ for E-MILES and FSPS, and $0.7-0.9$ for BC03 and CB19. For old stellar populations (${\rm log~{}Age~{}(yr)>9.5}$), the pPXF-derived mass-to-light ratios from each SPS models seem to align well with each other. However, in regions with young stellar populations, the NIR mass-to-light ratios rapidly decrease to ${\rm log}~{}M_{\ast}/L_{3.6~{}{\rm\mu m}}\sim-0.5$ to $-1.0$, with different slopes for each SPS model. This variation in the age–$M_{\ast}/L_{\rm 3.6\mu m}$ relation highlights the importance of considering stellar ages and SPS models for stellar mass estimation. E-MILES and FSPS tend to yield higher stellar masses with weaker correlations with stellar ages compared with other models. Conversely, BC03 and CB19 tend to predict substantially lower stellar masses in regions with young stellar populations (${\rm log~{}Age~{}(yr)}<9$). The relations derived from BC03 and CB19 exhibit steeper slopes ($\sim 0.4-0.6$) with slightly smaller scatters compared to those from E-MILES and FSPS ($\sim 0.2-0.3$). This trend naturally leads to the underestimation of stellar mass-to-light ratio for BC03 and CB19 in younger stellar regions relative to E-MIELS and FSPS.

Stellar metallicity seems to have a minor impact on the age–$M_{\ast}/L_{\rm 3.6\mu m}$ relation. Although the slopes may slightly increase in high-metallicity bins across all SPS models, the effect is subtle and less pronounced than the overall scatter of the relation. In low-metallicity bins, the scatters of the relation increase and the Spearman’s correlation coefficients decrease, suggesting that NIR mass-to-light ratios in metal-poor regions are prone to deviate from a simple linear relation with stellar age. Nonetheless, this metallicity dependence is far less significant than the age dependence in determining stellar mass-to-light ratios, as indicated by much lower correlation coefficients ($r_{S}\lesssim 0.1$) for the metallicity–$M_{\ast}/L_{\rm 3.6\mu m}$ relation compared to the age–$M_{\ast}/L_{\rm 3.6\mu m}$ relation. In conclusion, these findings highlight the importance of carefully considering stellar ages to achieve accurate stellar mass estimation using any SPS models, especially in late-type galaxies.