Red Supergiant problem viewed from the nebular phase spectroscopy of type II supernovae

When massive stars, with zero-age-main-sequence (ZAMS) mass $M_{\rm ZAMS}$ $>$ 8 $M_{\rm\odot}$ exhaust the nuclear fuel in their core, they will experience iron-core infall, and explode as a core-collapse (CC) supernova (SN). The CCSNe are classified based on the absorption features that emerge in the early phase spectroscopy; those with hydrogen lines, therefore probably possessing a massive hydrogen-rich envelope, are classified as type II supernovae (SNe II), while those without hydrogen lines are classified as stripped-envelope supernovae (SESNe). The readers may refer to Filippenko (1997); Gal-Yam (2017); Modjaz et al. (2019) for the classification of CCSNe.

The first SN II that have directly confirmed red supergiant (RSG) progenitor from pre-supernova imaging was SN 2003gd (Hendry et al., 2005), whose spectral energy distribution (SED) was consistent with those of field RSGs. As the sample of SNe II with directly identified progenitors from pre-SN imaging has grown, it has become increasingly evident that most SNe II originate from RSG progenitors, as predicted by stellar evolution theory. However, a discrepancy has emerged: while the most luminous field RSGs can have bolometric luminosity reaching log $L/L_{\rm\odot}$ = 5.5 dex (Davies et al. 2018)¹¹1Throughout this work, log $L$ refers to bolometric luminosities and are expressed in solar units unless otherwise specified., no evidence exists for such bright RSGs as progenitors of SNe II. Indeed, the most luminous RSG progenitor detected to date is that of SN 2009hd, with log $L$ = 5.24 dex, or $M_{\rm ZAMS}$ $\sim$ 20 $M_{\rm\odot}$, although the conversion from luminosity to $M_{\rm ZAMS}$ is model dependent. This apparent absence of bright and massive RSG progenitors for SNe II is referred to as the RSG problem (Smartt 2009; Walmswell & Eldridge 2012; Eldridge et al. 2013; Meynet et al. 2015; Smartt 2015; Davies & Beasor 2018; Strotjohann et al. 2024).

This observational discrepancy challenges stellar evolution theory. According to single-star models, only stars with $M_{\rm ZAMS}$ $>$ 30 $M_{\rm\odot}$ are predicted to experience sufficiently strong stellar winds to completely strip their hydrogen-rich envelopes (Meynet & Maeder 2000; Sukhbold et al. 2016), and one would expect to observe SNe II with $M_{\rm ZAMS}$ between 20 and 30 $M_{\rm\odot}$. Several theories have been proposed to explain the absence of RSG progenitors within this range. One possibility is the failed SN scenario, which suggests that RSGs within this mass range collapse to form black holes and disappear quietly without producing a bright explosion (O’Connor & Ott 2011; Horiuchi et al. 2014; Pejcha & Thompson 2015; Ertl et al. 2016; Müller et al. 2016; Sukhbold et al. 2016, 2018; Ebinger et al. 2019; Sukhbold & Adams 2020; Fryer et al. 2022; Temaj et al. 2024). Another explanation involves eruptive mass loss, where instabilities in RSGs of this mass range lead to significant mass ejections, stripping their hydrogen-rich envelopes. Such stars may instead explode as SESNe or interacting SNe, rather than SNe II (Smith et al. 2009; Yoon & Cantiello 2010; Smith & Arnett 2014; Meynet et al. 2015; Temaj et al. 2024).

Despite these theoretical investigations, efforts have been made to assess the significance of the RSG problem. Converting pre-SN magnitudes to bolometric luminosities depends on several assumptions, such as the spectral type of the progenitor, circumstellar dust properties (Walmswell & Eldridge 2012; Van Dyk et al. 2024), and bolometric corrections (Davies & Beasor 2018; Healy et al. 2024; Van Dyk et al. 2024; Beasor et al. 2025), all of which introduce substantial uncertainties. Recently, Healy et al. (2024) and Beasor et al. (2025) demonstrated that using a single bandpass for progenitor identification, as was done for many RSGs, can lead to systematic underestimations of luminosities, and found no statistical significant evidence of missing high luminosity RSGs in pre-SN images. Statistical limitations also play a role; Davies & Beasor (2018) argued that the RSG problem might be partly due to the small sample size of observed RSG progenitors. Additionally, Strotjohann et al. (2024) raised concerns about the impact of telescope sensitivity on RSG progenitor detection statistics.

Given these considerations, it is important to investigate other methodologies to infer $M_{\rm ZAMS}$ independently to cross-check the significance level of the RSG problem. One of the most frequently adopted techniques is modeling the light curve at plateau phase (Morozova et al. 2018; Martinez et al. 2022; Moriya et al. 2023). This method involves evolving models with different $M_{\rm ZAMS}$ until the onset of core-collapse and then injecting different amounts of energy and ⁵⁶Ni into the central region to trigger the explosion. The resultant light curve is compared with observation to determine these quantities. By employing KEPLER models as RSG progenitors, Morozova et al. (2018) found an upper $M_{\rm ZAMS}$ cutoff at 22.9 $M_{\rm\odot}$ for a sample of 20 SNe II. In a similar investigation on a larger sample, Martinez et al. (2022) found the upper cutoff at 21.3 $M_{\rm\odot}$. This approach has the advantage of allowing for relatively large samples, however, it also has limitations: the properties of the plateau light curve are mainly determined by the mass of the hydrogen-rich envelope $M_{\rm Henv}$ (when other properties, such as the explosion energy and the radius of the RSG, are fixed) rather than $M_{\rm ZAMS}$ itself (Kasen & Woosley, 2009; Dessart & Hillier, 2019; Goldberg et al., 2019; Hiramatsu et al., 2021; Fang et al., 2025a). The validity of this approach depends on the assumption of a unique relationship between $M_{\rm ZAMS}$ and the envelope mass, which may hold for single stars but can break down in the presence of a binary companion (Heger et al. 2003; Eldridge et al. 2008; Yoon et al. 2010; Smith et al. 2011; Sana et al. 2012; Smith 2014; Yoon 2015; Yoon et al. 2017; Ouchi & Maeda 2017; Eldridge et al. 2018; Fang et al. 2019; Zapartas et al. 2019, 2021; Chen et al. 2023; Ercolino et al. 2023; Fragos et al. 2023; Hirai 2023; Matsuoka & Sawada 2023; Sun et al. 2023, among many others) or uncertainties in stellar winds (Eldridge & Vink 2006; Mauron & Josselin 2011; Meynet et al. 2015; Davies & Beasor 2018, 2020; Wang et al. 2021; Massey et al. 2023; Vink & Sabhahit 2023; Yang et al. 2023; Zapartas et al. 2024, among many others).

In this work, we investigate the RSG problem using late-phase ($nebular$ phase) spectroscopy of SNe II, taken on $\sim$ 200 days after the explosion. During this phase, the spectroscopy is dominated by emissions lines, of particular importance is the oxygen emission [O I] $\lambda\lambda$6300,6363. The [O I] emission is considered as an important tool for measuring the oxygen content in the ejecta (Fransson & Chevalier 1989; Maguire et al. 2012; Jerkstrand et al. 2012, 2014; Kuncarayakti et al. 2015; Silverman et al. 2017; Dessart & Hillier 2020; Dessart et al. 2021; Fang et al. 2022), which is monotonically dependent on $M_{\rm ZAMS}$ and therefore the luminosity of the RSG progenitor Sukhbold et al., 2018; Takahashi et al., 2023. As a result, the $M_{\rm ZAMS}$ inferred from the strength of the [O I] line can be considered as an independent view point on the RSG problem from pre-SN images.

This paper is organized as follows: In §2,we describe the nebular spectroscopy sample and the methods used to process them. In §3, we introduce the method to determine $M_{\rm ZAMS}$ for individual SNe from [O I] emission, and establish the $M_{\rm ZAMS}$ distribution of the full sample. In §4, we correlate the $M_{\rm ZAMS}$, determined in §3, with the luminosities of the RSG progenitors from pre-SN images for a sub-sample of SNe II. This calibrated mass-luminosity relation is applied to the full sample to establish the luminosity distribution of their RSG progenitors, which is modeled with a power law function in §5 to assess the significance of the RSG problem. We discuss the physical implications in §6. Finally, we summarize our conclusions in §7.

In this work, we compile nebular spectroscopy of SNe II from the literature that meets the following criteria: (1) that the wavelength range must cover 5000 to 8500 $\rm{\AA}$, (2) that the spectra must be obtained more than 200 days after the explosion to ensure the nebular phase is reached, but not later than 450 days to allow for comparison with spectral models; (3) that the spectra are available on the Open Supernova Catalog (Guillochon et al. 2017), the Weizmann Interactive Supernova Data Repository (WISeREP; Yaron & Gal-Yam 2012) or the Supernova Database of UC Berkeley (SNDB; Silverman et al. 2012). For objects with multiple nebular spectra available, we select the one closest to 350 days post-explosion. This phase is chosen because it is late enough to ensure all SNe II are fully nebular, yet not so late that flux contributions from shock-circumstellar material (CSM) interaction become significant (Dessart et al. 2021; Rodríguez 2022; Dessart et al. 2023). The final sample consists of 50 SNe II, which are listed in Table A1 in the Appendix. While this sample does not encompass all SNe II nebular spectra in the literature, a size of $N\,=\,$50 is sufficient for statistical analysis.

The absolute or relative strengths of the [O I] emission line that emerges in the nebular spectroscopy of SNe II are useful indicators of the carbon-oxygen (CO) core mass and, consequently, the ZAMS mass of the progenitor. In this work, we use the fractional flux of the [O I] line within the wavelength range of 5000 to 8500 $\rm{\AA}$, $f_{\rm[O\,I]}$, as a diagnostic for the oxygen mass in the ejecta. We compare these measurements with model spectra from Jerkstrand et al. (2012) and Jerkstrand et al. (2014), following a methodology similar to that of Barmentloo et al. (2024) and Dessart et al. (2021). The wavelength range is chosen to encompasses all the observed spectra in the sample, and cover most of the main emission features in the optical band. This approach has an important advantage: because $f_{\rm[O\,I]}$ measures relative fluxes, it is unaffected by distance and flux calibration, and is insensitive to extinction in the host environment as long as it is not highly extincted, which typically constitutes one of the largest sources of uncertainty. However, to measure $f_{\rm[O\,I]}$ and make a meaningful comparison with the models, the observed spectra must first be standardized, as described below.

The nebular spectra of SNe II consist of multiple prominent emission lines, including [O I] $\lambda\lambda$6300,6363, H$\alpha$, and [Ca II] $\lambda\lambda$7291,7323, superimposed on a so-called pseudo-continuum formed by thousands of weak spectral lines. Figure 1 shows the model spectra of SNe II taken from Jerkstrand et al. (2012), normalized to their integrated flux within the wavelength range 5000–8500 $\rm{\AA}$. Hereafter, we refer to models from $M_{\rm ZAMS}$ = 12 $M_{\rm\odot}$ as M12 models, while those from $M_{\rm ZAMS}$ = 15 $M_{\rm\odot}$ and 19 $M_{\rm\odot}$ are referred to as M15 and M19 models, respectively. In the left panel of Figure 1, four spectral regions are specifically highlighted: 5450–5500, 6020–6070, 6850–6900, and 7950–8000 ${\rm\AA}$. These wavelength ranges do not contain strong emission lines, allowing the fluxes in these regions to be treated as pure pseudo-continuum (Barmentloo et al., 2024). As shown in the right panel of Figure 1, the average fluxes $\overline{f}({\lambda}_{\rm i},t)$ within these regions are independent of $M_{\rm ZAMS}$ but well determined by the spectral phase $t$, which are fitted by quadratic functions (the solid lines in the right panel of Figure 1) to estimate $\overline{f}({\lambda}_{\rm i},t)$ at arbitrary phases. This forms the basis for the approach to addressing contamination from the host environment.

The observed spectra of SNe 2014cx and 2004dj, normalized to their integrated flux within the wavelength range 5000–8500 $\rm{\AA}$, are illustrated in the left panels of Figure 2. The spectrum of SN 2014cx shows significant contamination from its host environment, as indicated by its unusual slopes, similar to SN 2012ec (Jerkstrand et al., 2015). Although the case of SN 2004dj is less extreme, the average fluxes in the aforementioned wavelength ranges consistently exceed those predicted by the spectral models at the same phase. This discrepancy suggests that the background emission might not have been completely removed during the processing of the raw observational data. Before measuring $f_{\rm[O\,I]}$, these residual fluxes are removed as follows:

1. 1.

   the observed flux $f_{\rm obs}$ in the rest frame are transformed to the standardized (or normalized) flux $f_{\rm norm}$ assuming:

   $$f_{\rm norm}=A\,\times\,(f_{\rm obs}-f_{\rm con}),$$

   (1)

   here $f_{\rm con}$ is the fluxes of the residual from the host, and $A$ is a normalized constant. The destination function, normalized flux $f_{\rm norm}$, should meet the requirement that when normalized to unity, the fluxes from the aforementioned 4 regions should be close to the model spectra at the same phase, which are estimated from the quadratic fits in the right panel of Figure 1. For most SNe in the sample, $f_{\rm con}$ is assumed to be a constant for simplicity. However, for 3 objects (SNe 2012ch, 2012ec, and 2014cx) that exhibit unusual spectral slopes due to the contamination of their bright host environments, a quadratic form of $f_{\rm con}$ is applied:

   $\displaystyle f_{\rm con}=b_{\rm 2}\,\lambda^{2}+b_{\rm 1}\,\lambda+b_{\rm 0},$

   where $\lambda$ is the wavelength.

2. 2.

   The Python package scipy.optimize is imported to find the pair {$A$, $f_{\rm con}$} (or {$A$, $b_{\rm 0}$, $b_{\rm 1}$, $b_{\rm 2}$} for $f_{\rm con}$ in quadratic form) that minimize the following quantities:

   	$\displaystyle r_{\rm 0}$	$\displaystyle=\int_{\rm 5000\,\AA}^{{\rm 8500\,\AA}}f_{\rm norm}~{}d\lambda-1$
   	$\displaystyle r_{\rm i}$	$\displaystyle=\overline{f}_{\rm norm}({\lambda}_{\rm i})-\overline{f}_{\rm model}({\lambda}_{\rm i},t),~{}i=1,2,3,4.$

   Here $\overline{f}_{\rm norm}({\lambda}_{\rm i})$ is the average fluxes of the observed spectra after normalization within the 4 selected wavelength ranges, and $\overline{f}_{\rm model}({\lambda}_{\rm i},t)$ is the average pseudo-fluxes of the spectral models at the same phase $t$, estimated from the quadratic fit (the solid lines in the right panel of Figure 1). These procedures ensure that, after the transformation described by Equation 1, the integral of the normalized flux within 5000 to 8500 ${\rm\AA}$ equals unity. Moreover, the fluxes of the pseudo-continuum within the specific regions are aligned with those of the models, allowing for fair comparison.

The above procedures are applied to SN 2014cx and 2004dj, where quadratic and constant forms of $f_{\rm con}$ are assumed respectively, and the resultant $f_{\rm con}$ are shown as the black dashed lines in the left panels of Figure 2. The normalized fluxes, after $f_{\rm con}$ is removed, are shown in the right panels and compared with the spectral models at similar phases, showing that the pseudo-continuum fluxes of the normalized observed spectra are consistent with those of the models. After the spectra is normalized to $f_{\rm norm}$, we fit the lines in the range of 6100 to 6800 ${\rm\AA}$ with multi-Gaussian functions: (1) two Gaussians with the same standard deviation and peaks separated by 63 ${\rm\AA}$ to represent the [O I] doublet, (2) one Gaussian centered near 6563 ${\rm\AA}$ (with a small allowed shift) to represent H$\alpha$, and (3) an additional Gaussian with an arbitrary center to account for a spectral feature commonly observed between [O I] and H$\alpha$ in many SNe II nebular spectra. The fluxes of [O I] and H$\alpha$ are measured by integrating the fitted profiles. Figure 1 illustrates this fitting procedure using SNe 2014G and 2023ixf as examples. Although SN 2023ixf exhibits a complex [O I] line profile (see, e.g. Ferrari et al. 2024; Fang et al. 2025b), likely reflecting intricate ejecta geometry (Fang et al. 2024), this study focuses solely on the integrated flux, and these complexities are not considered.

The uncertainties in the fractional flux of [O I] (H$\alpha$), $f_{\rm[O\,I]}$ ($f_{{\rm H}\alpha}$), come mainly from uncertainties in subtracting the contamination flux $f_{\rm con}$. This is quantified using a Monte Carlo method: the original observed spectra are first smoothed, and the fluxes in the four fitting regions are replaced with the smoothed fluxes, augmented by random noises. The noise level in each region is estimated as the standard deviation of the difference between the original flux and the smoothed flux within it. The determination of $f_{\rm con}$ and the measurement of the [O I] (or H$\alpha$) fluxes are then repeated 1000 times, following the same procedure. In each trial, the pseudo-continuum fluxes are allowed to randomly vary by at most 20% (0.08 dex; see the scatter lever in the right panel of Figure 3). The standard deviations of these measurements are taken as the uncertainties in the fractional line fluxes.

In addition to $f_{\rm[O\,I]}$, we further define its regulated form as

i.e., the fractional flux of [O I] after the H$\alpha$ line is subtracted from the spectrum. The motivation for introducing $f_{\rm[O\,I],reg}$ is to address the uncertainties associated with H$\alpha$ arising from various factors, which will be discussed in §3.

In the upper panel of Figure 4, the measured fractional fluxes of [O I], $f_{\rm[O\,I]}$, are compared with the spectral models from Jerkstrand et al. (2012). A significant proportion of the SNe in the sample cluster between the M12 and M15 tracks. Compared to the results of Valenti et al. (2016), the size of the nebular spectroscopy sample is substantially larger, and more objects are found to lie above the M15 tracks. However, no single object provides evidence for a progenitor more massive than the M19 model. Notably, with the exception of SNe 2015bs (Anderson et al. 2018) and 2017ivv (Gutiérrez et al. 2020), no individual SN in the sample has a derived $M_{\rm ZAMS}$ exceeding 17 $M_{\rm\odot}$.

The fractional fluxes of the M9 models from Jerkstrand et al. (2018), computed from the progenitor model with $M_{\rm ZAMS}$ = 9 $M_{\rm\odot}$ using KEPLER in Sukhbold et al. (2016), are plotted for reference. These models include two sets: one representing the total SN spectrum (M9, SN) and another where the core material is replaced with H-zone material (M9, H-zone). For the former case, the fractional [O I] flux exceeds the M12 track, whereas in the latter, it falls below the M12 track. Despite having a lower oxygen mass, $f_{\rm[O\,I]}$ of the M9 models are not consistently lower than the M12 models.

This non-monotonic behavior may stem from several factors. First, the M9 progenitor undergoes a strong silicon flash, which could alter its core properties compared to the M12, M15, and M19 models, where nuclear burning proceeds stably throughout evolution (Sukhbold et al. 2016). Furthermore, these models assume lower explosion energies, leading to relatively narrow emission lines with a full width at half maximum (FWHM) of $\sim$ 1000 km s^-1, a feature not commonly observed in most SNe II.

Given these unique characteristics, the non-monotonic [O I] flux behavior, and uncertainties regarding whether the full SN models or pure H-zone models should be applied, the M9 models are not used to refine $M_{\rm ZAMS}$ estimates below the M12 track. Notably, SNe 2005cs, 2008bk, and 2018is, which were compared with M9 models in Jerkstrand et al. (2018)²²2For SN 1997D, the wavelength range and the spectral phase do not meet the criteria in this work and Dastidar et al. (2025), exhibit fractional [O I] fluxes below the M12 track, suggesting that the M12 models already capture the low-mass nature of their progenitors.

The comparison of $f_{\rm[O\,I],reg}$ and the spectral models are shown in the lower panel of Figure 4. The motivation for the introduction of $f_{\rm[O\,I],reg}$, i.e., removing H$\alpha$ from $M_{\rm ZAMS}$ estimates, stems from the fact that this line can be affected by many factors as will be discussed below.

Growing evidences suggest that some SNe II experienced partially stripping of the hydrogen-rich envelope prior to their explosions. A reduced envelope mass $M_{\rm Henv}$ is necessary to explain the short duration light curves observed in several individual SNe, such as SNe 2006Y, 2006ai, 2016egz (Anderson et al. 2014; Hiramatsu et al. 2021), 2018gj (Teja et al. 2023), 2020jfo (Teja et al. 2022), 2021wvw (Teja et al. 2024), 2023ixf (Fang et al. 2025b; Hsu et al. 2024) and 2023ufx (Tucker et al. 2024; Ravi et al. 2024). Fang et al. (2025a) further suggest that, to account for the observed diversity in SNe II light curves (e.g., Anderson et al. 2014; Valenti et al. 2016; Anderson et al. 2024), approximately half of SNe II must have stripped their envelopes to about $\sim$ 4.0 $M_{\rm\odot}$ (see also Hiramatsu et al. 2021).

The spectral models from Jerkstrand et al. (2012) assume a massive hydrogen-rich envelope, and therefore may not accurately reproduce the H$\alpha$ flux if the hydrogen-rich envelope is partially removed. Although the exact relationship between H$\alpha$ flux and $M_{\rm Henv}$ has not yet been fully established, a pioneering study by Dessart & Hillier (2020) shows that for SNe II with $M_{\rm Henv}$ $<$ 3 $M_{\rm\odot}$, the H$\alpha$ flux decreases dramatically, while other parts of the spectrum are much less affected (see also Ravi et al. 2024). From observation, SNe II with low $M_{\rm Henv}$, inferred from plateau light curve modeling, also show relatively weak H$\alpha$ in nebular spectroscopy (Teja et al. 2022, 2023; Ravi et al. 2024; Fang et al. 2025b), and faster declines in their radiative tails during the late phases, compared to the cases with full $\gamma$-ray trapping (Anderson et al. 2014; Gutiérrez et al. 2017).

Furthermore, H$\alpha$ can be illuminated by shock-CSM interaction, a process not included in the models of Jerkstrand et al. (2012), introducing another source of uncertainty in the H$\alpha$ flux. As demonstrated in Dessart et al. (2023), the contribution from shock-CSM interaction can increase the integrated flux by enhancing the H$\alpha$ line in nebular phase, and it can explain the H$\alpha$ features for several objects, including SNe 2014G and 2013by. However, since currently we are still lacking consistent interacting models of SNe II in the nebular phase, it remains challenging to quantify the size of this effect.

Given these considerations, we introduce $f_{\rm[O,I],reg}$ as a simple yet effective approach to reduce the uncertainties in H$\alpha$ by removing it from the measurement. Indeed, $f_{\rm[O,I],reg}$ and $f_{\rm[O,I]}$ represent two limiting scenarios: (1) the first scenario, based on $f_{\rm[O\,I],reg}$, assumes that the change in the integrated flux, whether due to the increase in $\gamma$-ray leakage from reduced $M_{\rm Henv}$ or the contribution from interaction power, only affects the H$\alpha$ flux, leaving other spectral features unaffected; (2) the second scenario, based on $f_{\rm[O\,I]}$, assumes that the [O I] flux scales directly with the integrated flux. In this case, if the integrated flux decreases (due to additional $\gamma$-ray leakage when $M_{\rm Henv}$ is low) or increases (due to interaction power) by 50%, the [O I] flux is also varied by the same fraction. Taking the low $M_{\rm Henv}$ models in Dessart et al. (2021) and the interacting SNe II models in Dessart et al. (2023) as references, the actual situation likely falls between these two extremes.

With $f_{\rm[O\,I]}$ and $f_{\rm[O\,I],reg}$, this work measures $M_{\rm ZAMS}$ for individual SNe II in the sample following the below procedures:

• •

  For $f_{\rm[O\,I]}$ above the M12 track: $M_{\rm ZAMS}$ is determined through linear interpolation with the model tracks. For example, if an SNe is observed at phase $t$ with $f_{\rm[O\,I]}$ between the M12 and M15 tracks, we first estimate $f_{\rm[O\,I]}$ of the models at $t$ through interpolation, and then interpolate between the models again to estimate $M_{\rm ZAMS}$ based on $f_{\rm[O\,I]}$;

• •

  For $f_{\rm[O\,I]}$ below the M12 track: $M_{\rm ZAMS}$ is assumed to follow a uniform distribution between 10 and 12 $M_{\rm\odot}$, with a Gaussian tail (standard deviation of 1 $M_{\rm\odot}$) extending to higher masses.

• •

  The same methodology is applied to measurements using $f_{\rm[O\,I],reg}$;

• •

  The final adopted $M_{\rm ZAMS}$ combines the measurements with $f_{\rm[O\,I]}$ and $f_{\rm[O\,I],reg}$: all $M_{\rm ZAMS}$ values within the range defined by the $M_{\rm ZAMS}$ measured from $f_{\rm[O\,I]}$ and $f_{\rm[O\,I],reg}$ are assumed to be uniformly distributed. Figure 5 shows 3 examples of the final adopted $M_{\rm ZAMS}$ distributions for individual SNe. For individual SNe II, its $M_{\rm ZAMS}$ is not a single value with Gaussian uncertainty but follows an irregular distribution that cannot be described analytically. Throughout this work, $M_{\rm ZAMS}$ and its uncertainty represent the median value and the 68% confidence interval (CI; determined by the 16th and 84th percentiles) of this distribution.

In Figure 6, we compare the $M_{\rm ZAMS}$ estimated from $f_{\rm[O\,I]}$ and $f_{\rm[O\,I],reg}$. The result shows that $f_{\rm[O\,I],reg}$ generally predicts lower $M_{\rm ZAMS}$, with differences reaching up to $\sim\,1\,M_{\rm\odot}$ in some cases. This suggests that spectral models may overestimate the H$\alpha$ flux.

The distribution of $M_{\rm ZAMS}$ of the full sample is shown in Figure 7, calculated using a Monte Carlo method similar to that described in Davies & Beasor (2018): in each trial, for each individual SN, a mass is randomly sampled from its $M_{\rm ZAMS}$ distribution. For those with upper limit, a mass is randomly sampled from a distribution that is uniform between 10 to 12 $M_{\rm\odot}$ (as shown in Figure 5). The simulated sample is then sorted. This process is repeated 10,000 times. The median $M_{\rm ZAMS}$ for each rank, from the SN with the lowest $M_{\rm ZAMS}$ to the one with the highest $M_{\rm ZAMS}$, is calculated and represented by the black line in Figure 7. The shaded regions indicate 95 and 99.7% CI, while the 68% CI is not filled for illustration purposes.

In the previous section, we developed a method to estimate $M_{\rm ZAMS}$ from nebular spectroscopy (hereafter denoted as $M_{\rm ZAMS,neb}$ for clarity). The aim of this study is to assess the significance of the RSG problem, i.e., the lack of luminous RSG progenitor for SNe II. For this purpose, our next step is to convert $M_{\rm ZAMS,neb}$ into a luminosity scale, which can, in principle, be done using the $M_{\rm ZAMS}$–luminosity relation (hereafter referred to as MLR) derived from stellar evolution models. However, $M_{\rm ZAMS,neb}$ mainly reflects the oxygen content in the ejecta. Inferring $M_{\rm ZAMS}$ from nebular spectroscopy relies on the underlying relation between the synthesized oxygen mass $M_{\rm O}$ and $M_{\rm ZAMS}$. While $M_{\rm O}$ is a monotonic function of $M_{\rm ZAMS}$, with more massive stars generally synthesizing more oxygen, the transformation from $M_{\rm O}$ (nebular spectroscopy) to $M_{\rm ZAMS}$, and subsequently to log $L$ strongly depends on the microphysics of the stellar evolution code, particularly the assumptions about internal mixing. Converting $M_{\rm ZAMS,neb}$ to log $L$ essentially reflects a $M_{\rm O}$-log $L$ relation, which, as will be demonstrated in the discussion in §4.1, is subject to significant uncertainties. To address this, an empirical MLR based on observations is established in §4.2, and its robustness is tested in §4.3.

In this section, we discuss the global properties of the log $L$ distribution of the RSG progenitor, estimated from $M_{\rm ZAMS,neb}$ using the calibrated MLR. There are mainly two methods to assess the significance of the RSG problem, (1) Comparison of the LDF of RSG progenitors with that of RSGs in other galaxies (field RSGs). See Rodríguez (2022) as an example; (2) Modeling the LDF of RSG progenitors using analytical functions, typically power-law forms with upper and lower cutoffs. These methods have distinct focuses as well as pros and cons. While method (1) has the advantage of being model independent, it also has several limitations:

• •

  Key focus: Using method (1), Rodríguez (2022) concludes that the $N$(log $L$ $>$ 5.1)/$N$(log $L$ $>$ 4.6) ratio of SNe II progenitors is smaller the that of the RSGs in Large Magellanic Cloud (LMC) and other galaxies, implying fewer bright RSG than expected. Here $N$(log $L$ $>$ 4.6) refers to the number of RSGs with log $L$ $>$ 4.6 dex, and similarly for $N$(log $L$ $>$ 5.1). However, this approach does not directly address the existence of an upper luminosity cutoff in the progenitor population. The RSG problem fundamentally concerns the presence of such a cutoff, which could result from factors like explodability or pre-supernova mass-loss mechanisms (discussed later in §6).

• •

  Evolutionary presentation: The field RSGs are typically less evolved than the progenitors of SNe II, which raises questions about whether they really represent the final evolutionary states of massive stars. For example, RSGs with luminosities as high as log $L\,\sim$ 5.5 dex may remain near the Hayashi line during helium burning and could be observed as a field RSG, but it could transition away from this state in later evolutionary phases due to processes like mass stripping from eruptive activities. As a result, they may not explode as SNe II, potentially relaxing the observed difference in the $N$(log $L$ $>$ 5.1)/$N$(log $L$ $>$ 4.6) ratio between progenitor RSGs and field RSGs.

Given these considerations, we employ method (2) for the statistical analysis of the LDF, similar to the one applied in Davies & Beasor (2020) (sorting method). For further details, readers are encouraged to refer to that paper. Here, we briefly describe the workflow:

(1) In the previous section, we have established the LDF using a Monte Carlo method. Especially, we derive the distribution of the luminosity of the $i$-th SN, denoted as $P_{i,\rm obs}$ (log $L$);

(2) Next, we construct the model LDF, which is in the power-law form

bounded by the lower limit $L_{\rm low}$ and the upper limit $L_{\rm up}$. For each set of {$\Gamma_{L},{\rm log}\,L_{\rm low},{\rm log}\,L_{\rm up}$}, we again estimate the LDF using the same method described in the previous section: we perform 10,000 Monte Carlo simulation, and in each trial, we draw 50 (the observed sample size in this work) log $L$ values from the bounded power-law distribution (representing the scatter points in Figure 10). The uncertainties are assigned according to the ranks. For example, for the faintest progenitor in the random sample, it is assumed to follow the same log $L$ distribution as SN 2013am established through Monte Carlo sampling in §4.2, but shifted by a constant to align the median value. For each simulated SN, a value is randomly drawn from its associated log $L$ distribution, and the full sample is sorted again. This second sorting step mimics the ranking method used in Davies & Beasor (2020). The luminosity distribution of the $i$-th SN from the 10,000 trials is denoted as $P_{i,\rm model}$ (log $L$);

(3) For a fixed set of {$\Gamma_{L},{\rm log}\,L_{\rm low},{\rm log}\,L_{\rm up}$}, the probability that the model produces the observed log $L$ of the $i$-th progenitor is

The likelihood function is written as

After the likelihood function is established, we use the Python package emcee to infer the optimized {$\Gamma_{L},{\rm log}\,L_{\rm low},{\rm log}\,L_{\rm up}$} and the associated uncertainties (Foreman-Mackey et al. 2013). The setup is as follow: we use 60 walkers (20 walkers per parameter), and their initial positions {$\Gamma_{L},{\rm log}\,L_{\rm low},{\rm log}\,L_{\rm up}$}₀ are randomly sampled from {$U{\rm(-5,2)},U{\rm(4.0,4.8)},U{\rm(4.8,6.0)}$}, the prior distributions of these parameters. Here $U$($a,b$) represents uniform distribution between $a$ and $b$. We then allow the walkers to explore the parameter space freely using emcee.run_mcmc to run for 1000 chains, after which the result converges. The corner plot of the parameters and their uncertainties are shown in Figure 12. The optimized parameters are:

It is interesting to see that, although a different SNe sample and different method are employed to derive the LDF of the RSG progenitor, the results match quite well with Davies & Beasor (2020), where they find

Although the optimize ${\rm log}\,L_{\rm up}$ is lower than 5.5 dex, consistent with the RSG problem, the 97.5 percentile (+2$\sigma$) of the log $L_{\rm up}$ distribution is 5.44 dex, and the 99.8 percentile (+3$\sigma$) is 5.63 dex. This suggests that the significance of this issue is at the 2$\sigma$ to 3$\sigma$ level.

As discussed in Davies & Beasor (2020), part of the uncertainties in ${\rm log}\,L_{\rm up}$ and ${\rm log}\,L_{\rm low}$ arises from their degeneracy with $\Gamma_{L}$. Specifically, for a steeper power-law distribution, observations are more likely to detect faint objects, reducing the probability of identifying bright progenitors. This reduced detection probability allows a higher ${\rm log}\,L_{\rm up}$ to remain consistent with the data. Similarly, a sharp power-law implies a rapid increase in probability density as ${\rm log}\,L$ approaches ${\rm log}\,L_{\rm low}$. To prevent divergence at the lower end, the cutoff ${\rm log}\,L_{\rm low}$ shifts upward to maintain normalization. If the power-law index $\Gamma_{L}$ is fixed to -1.675, as would be expected if the progenitors in the sample follow the Salpeter form of initial mass function (IMF; characterized by $dN/dM_{\rm ZAMS}\,\propto\,M_{\rm ZAMS}^{-2.35}$; Salpeter 1955; Davies & Beasor 2020), the similar emcee routine returns

as shown in Figure 13. This result shows that fixing $\Gamma_{L}$ leads to higher optimized values for {${\rm log}\,L_{\rm low},{\rm log}\,L_{\rm up}$} compared to cases where $\Gamma_{L}$ is allowed to vary freely. However, the corresponding uncertainties in these parameters decrease because the degeneracy between $\Gamma_{L}$ and the cutoffs is removed. As a result, the RSG problem persists at a significance level of 2$\sigma$ (5.49 dex at 97.5 percentile) to 3$\sigma$ (5.57 dex at 99.8 percentile).

One major concern is the completeness of the sample. Indeed, if the progenitor sample adheres to the Salpeter IMF, the expected $\Gamma_{L}$ is -1.675, a value within the $1\sigma$ uncertainty reported by Davies & Beasor (2020). In this work, the sample size is increased to $N$ = 50. While the uncertainties of the luminosity cutoffs are comparable to Davies & Beasor (2020), $\Gamma_{L}$ is now constrained to a narrower range, and the expected value of -1.675 falls outside the 68% CI. The relatively flat LDF suggests that the sample is probably incomplete, likely missing some low-luminosity progenitors. To address this, we account for sample incompleteness in two ways:

(1) Excluding Low-Luminosity Progenitors. For the first attempt, we only consider bright progenitors. Previous studies suggest that for progenitors with log $L\gtrsim 4.6$–4.7, the LDF is consistent with those of field RSGs (Rodríguez, 2022; Strotjohann et al., 2024; Healy et al., 2024)⁴⁴4It should be cautious that, field RSGs are less evolved than SN progenitors, raising questions about whether they truly reflect the properties of RSGs at the onset of core collapse.. Further, for SN progenitors with log $L$ $<$ 4.6, the measured luminosity may be unreliable, otherwise some SNe would have $M_{\rm ZAMS}$ $<$ 8 $M_{\rm\odot}$, below the minimum mass required for a single star to undergo core collapse (Heger et al., 2003). Based on these considerations, we adopt a cutoff at log $L=4.6$, retaining only bright SNe progenitors. This reduces the sample size to $N$ = 33, comparable to the sample size of Davies & Beasor (2020). Repeating the emcee routine gives

The smaller sample size relaxes parameter constraints, increasing uncertainties to levels comparable to Davies & Beasor (2020). Notably, the expected value of -1.675 now falls within the uncertainty of $\Gamma_{L}$, while the optimized value of ${\rm log}\,L_{\rm up}$ increases by 0.04 dex (10% in linear scale). Consequently, the significance of the RSG problem is reduced to below 1$\sigma$;

(2) Including pseudo SNe for low-luminosity progenitors. In the second approach, we address the faint side of the LDF. While these low-luminosity progenitors do not directly affect the estimate of log $L_{\rm up}$, they play a crucial role in constraining $\Gamma_{L}$ and log $L_{\rm low}$. By changing the overall shape of the LDF, they can indirectly affect the estimate of log $L_{\rm up}$. For instance, decreasing $\Gamma_{L}$ can result in higher values of ${\rm log}\,L_{\rm up}$ as discussed above.

To investigate this effect, we introduce pseudo SNe, artificial data points representing the missing low-luminosity progenitors. These pseudo SNe are assigned the same $M_{\rm ZAMS,neb}$ distribution as objects below the M12 track (middle panel of Figure 5). The exact luminosity distribution of these missing progenitors is unknown, and our assignment is somewhat arbitrary. However, it is reasonable: as shown in Figure 9, the faintest progenitors correspond to SNe located below the M12 track, most of which are low-luminosity SNe II. These SNe are characterized by low ⁵⁶Ni production, making them faint and difficult to detect during the nebular phase (Pastorello et al., 2004; Spiro et al., 2014; Murai et al., 2024), and thus our nebular spectra sample in this range is probably incomplete. By varying the number of pseudo SNe ($N_{\rm pseudo}$) added to the observed sample, we repeat the previous analysis to infer the optimized parameters {$\Gamma_{L},{\rm log}\,L_{\rm low},{\rm log}\,L_{\rm up}$} with the emcee routine. The results, along with their 68% and 95% CIs, are shown in Figure 14.

As expected, $\Gamma_{L}$ and log $L_{\rm low}$ decrease with the increase of $N_{\rm pseudo}$, while ${\rm log}\,L_{\rm up}$ remains largely unaffected. For $N_{\rm pseudo}$ = 30, when the number of the pseudo SNe becomes comparable to the observed sample, we obtain

At first glance, the result might seem contradictory to the earlier discussion, where fixing $\Gamma_{L}$ to -1.675 resulted in larger log $L_{\rm up}$. However, the key distinction lies in sample size. Adding pseudo SNe effectively enlarges the sample. While lower $\Gamma_{L}$ (a sharper power-law distribution) biases detections toward low-luminosity progenitors, allowing for larger ${\rm log}\,L_{\rm up}$, a larger sample size reduces the likelihood of missing high-luminosity progenitors. The competition between these factors stabilizes ${\rm log}\,L_{\rm up}$ at an approximately constant value. Thus, the RSG problem persists at a significance level of $2\sigma$ to $3\sigma$, as shown in the right panel of Figure 14.

In Morozova et al. (2018) and Martinez et al. (2022), the cutoff masses of the $M_{\rm ZAMS}$ mass distribution for large samples of SNe II were investigated through plateau-phase light curve modeling. Converting the luminosity cutoffs derived in this work into $M_{\rm ZAMS}$ values is essential for a direct comparison with these studies and for assessing the robustness of the results. We perform this conversion using the MLR from the KEPLER models, i.e., first convert the luminosity scales (say, log $L_{\rm up}$) to $M_{\rm He\,core}$ via Equation 2, which is subsequently converted to $M_{\rm ZAMS}$ using the $M_{\rm ZAMS}$-$M_{\rm He\,core}$ relation in Sukhbold et al. (2016). Although this introduces uncertainties associated with different stellar evolution codes (as discussed in §4), the progenitor models of Morozova et al. (2018) and Martinez et al. (2022) follow a similar MLR, allowing for a fair comparison.

The upper mass cutoff ($M_{\rm up}$) converted from log $L_{\rm up}$ derived in this work is 20.63${}^{+2.42}_{-1.64}$ $M_{\rm\odot}$, where the quoted uncertainties define the 68% CI. Similarly, we find a lower mass cutoff at $M_{\rm low}$ = 8.95${}^{+0.27}_{-0.32}$ $M_{\rm\odot}$.⁵⁵5For log $L\,<\,$4.3 (corresponding to $M_{\rm ZAMS}$ = 9 $M_{\rm\odot}$), the estimation is based on extrapolation. In Figure 15, we compare the $M_{\rm up}$ and $M_{\rm low}$ with the measurements from plateau phase light curve modeling (Morozova et al. 2018; Martinez et al. 2022). The purple strip is the $M_{\rm ZAMS}$ converted from the maximum log $L_{\rm prog}$($L_{\rm[OI]}$) in Rodríguez (2022) (5.063 $\pm$ 0.077; see Table 9 in that paper).

Despite the use of different methods, including SN progenitor luminosity analysis (Davies & Beasor 2018, 2020), plateau-phase light curve modeling (Morozova et al. 2018; Martinez et al. 2022), and nebular-phase spectroscopy (Rodríguez 2022; this work), the inferred $M_{\rm up}$ consistently falls within the range of $\sim$ 18 to 23 $M_{\rm\odot}$, well below $M_{\rm ZAMS}$ $\sim$ 29.4 $M_{\rm\odot}$ converted from log $L$ = 5.5. Although the discrepancy between $M_{\rm up}$ and the threshold 29.4 $M_{\rm\odot}$ is at the level of 1 to 3$\sigma$, its persistence across different methodologies suggests that it may reflect a real physical problem. It would be interesting to explore the physical implications of this discrepancy.

Several theories have been proposed to explain the dearth of SNe II with luminous progenitor. The first one involves the ”explodability” of massive stars. Sukhbold et al. (2018) studied the core structures of a large grid of progenitors with varying physical parameters and found that the upper mass limit for SNe II appears to converge at $\sim$ 20 $M_{\rm\odot}$. In their models, progenitors with $M_{\rm ZAMS}$ $\sim$ 20 to 25 $M_{\rm\odot}$ experience collapse into black holes, resulting in failed SNe (see also Table 6 and Figure 19 of Sukhbold et al. 2016). However, multi-dimensional models from Burrows et al. (2024a) suggest successful explosions in this mass range. In a subsequent study, Burrows et al. (2024b) showed that even for a 23 $M_{\rm\odot}$ progenitor that forms a black hole, the model produces an explosion rather than remaining quiescent. Thus, the role of explodability in explaining the absence of massive progenitors for SNe II remains a topic of ongoing debate.

Another hypothesis involves the surface properties of RSG progenitors. Under this scenario, massive stars may still explode but as SNe types other than SNe II. According to the RSG models in Sukhbold et al. (2016), single-star evolution predicts that the hydrogen-rich envelope is fully removed by stellar winds only in stars with $M_{\rm ZAMS}$ $\gtrapprox$ 30 $M_{\rm\odot}$. However, if mass-loss from stellar winds is stronger than assumed in current stellar evolution models, this threshold could be reduced to $\sim$ 20 $M_{\rm\odot}$ (Meynet et al. 2015). There is also observational evidence suggesting that more massive stars are more likely to form in close binary systems (Moe & Di Stefano 2017; Moe et al. 2019). Such systems can efficiently strip the hydrogen-rich envelope, skewing the mass (and luminosity) distribution of SNe II progenitors toward lower values. Further, a luminosity of log $L\,\sim$ 5.2 dex (log $L_{\rm up}$ found in this work and Davies & Beasor 2020) is already sufficient to trigger pulsation (Soraisam et al. 2018; Goldberg et al. 2020). While pulsation-driven mass-loss is not included in stellar evolution codes, Yoon & Cantiello (2010) demonstrated that once initiated, it becomes a runaway process capable of stripping nearly the entire hydrogen envelope. In this scenario, the pulsating RSGs would eventually explode as stripped-envelope SNe or interacting SNe (Smith et al. 2009) rather than SNe II. In a recent work, Cheng et al. (2024) propose an eruptive mass-loss mechanism, under which progenitor models with $M_{\rm ZAMS}$ $\gtrapprox$ 20 $M_{\rm\odot}$ will end their lives as blue supergiant. This could explain the apparent absence of RSG progenitors above this mass range. Further investigations on the instabilities of massive stars are important to fully understand these processes and their potential connections with the RSG problem.

The RSG problem—namely, the observed absence of luminous RSG progenitors for SNe II—raises fundamental challenges about our understanding of massive star evolution and SN explosion mechanisms. In this work, we investigate this topic through a statistical analysis of nebular spectroscopy for a large sample of SNe II. Since nebular spectroscopy provides an independent estimate of the $M_{\rm ZAMS}$, it offers an important cross-check on the RSG problem, complementing results from pre-SN imaging of RSG progenitors.

To achieve this, we first estimate $M_{\rm ZAMS,neb}$ for individual SNe by comparing the fractional flux of the [O I] emission emerging in the nebular spectroscopy with spectral models. The resulting $M_{\rm ZAMS,neb}$ values are then compared with the observed luminosities of RSG progenitors for SNe with detected progenitors, revealing a strong and statistically significant correlation. Using this empirically calibrated mass-luminosity relation, we establish the progenitor luminosity distribution function (LDF) for the full sample. The brightest progenitor in our sample is that of SN 2017ivv, with log $L$ = 5.33${}^{+0.21}_{-0.18}$ dex. Notably, there is no progenitor exceeding log $L$ $>$ 5.5 dex—the empirical upper luminosity limit for field RSGs—at a significance level of approximately 1$\sigma$.

The LDF is modeled using a power-law function with upper and lower luminosity cutoffs, adopting a Monte Carlo method similar to that of Davies & Beasor (2020). Despite differences in sample selection and methods for measuring log $L$, our results are consistent with Davies & Beasor (2020). We find an upper luminosity cutoff of log $L_{\rm up}$ = 5.21${}^{+0.09}_{-0.07}$ dex, with the absence of progenitors above log $L$ = 5.5 dex being statistically significant at the 2$\sigma$ to 3$\sigma$ level.

Finally, we convert the luminosity cutoffs derived in this work back to the $M_{\rm ZAMS}$ scale using the mass-luminosity relation from KEPLER models and compare these with constraints from plateau light curve modeling. Despite methodological differences, both approaches consistently indicate an upper mass cutoff for SNe II progenitors below 29 $M_{\rm\odot}$ at a significance level of 1–3$\sigma$. While each individual method provides only marginal significance, their consistency suggests that the lack of luminous RSG progenitors is likely a real physical problem. This finding highlights the need for further investigations on the explodability of high-mass progenitors and the late-phase instabilities of massive stars to fully understand their potential connections to the RSG problem.