Light Curves of Type IIP Supernovae from Neutrino-driven Explosions of Red Supergiants Obtained by a Semi-analytic Approach

Core-collapse supernovae (CCSNe) are the spectacular explosions that mark the death of massive stars with zero-age main-sequence masses ($M_{\rm ZAMS}$) greater than $\sim$8–10 $M_{\odot}$ (e.g., Ibeling & Heger, 2013) in the case of single-star progenitors. Understanding the explosion mechanism of CCSNe has become the equivalent of a millennium problem of modern astrophysics. CCSNe have great importance as a source of multi-messenger events¹¹1SN 1987A was the first-ever extragalactic astronomical multi-messenger event. and of compact remnants that they leave behind, and are the origin of most heavy elements in the universe.

Many of the latest three-dimensional (3D) simulations with sophisticated physics inputs, such as accurate modeling for the neutrino transport, have yielded successful CCSN explosions (see, e.g., Lentz et al., 2015; Müller et al., 2017a; Ott et al., 2018; Burrows et al., 2019; Bollig et al., 2021), which supports the view that most CCSNe are powered by the neutrino-driven mechanism aided by hydrodynamical instabilities (see the reviews of Bethe, 1990; Janka, 2012; Burrows & Vartanyan, 2021; Müller, 2020). There is, however, still an ongoing discussion on several key issues. Using phenomenological models, considerable progress has been made in determining how the pre-collapse stellar structure impacts which stars successfully explode and which ones fail and make black holes (O’Connor & Ott, 2011; Ugliano et al., 2012; Ertl et al., 2016; Müller et al., 2016; Sukhbold et al., 2016; Pejcha & Thompson, 2015), but the best structural correlates for “explodability” and the parameter space for neutron star and black hole formation are still debated in supernova theory (Couch et al., 2020; Tsang et al., 2022). Beyond the question of explodability, both detailed multi-dimensional simulations and phenomenological models have shed some light on the relation between progenitors and their explosion properties, such as the remnant mass, explosion energy, and nucleosynthesis yields as critical input for chemogalactic evolution (Nomoto et al., 2013), but the key challenge is now to more rigorously validate the emerging theoretical picture using observational data.

Direct multi-messenger probes of the CCSN mechanism include gravitational waves (GWs) (e.g., Fryer & New, 2011; Evans & Zanolin, 2017; Kalogera et al., 2019; Abdikamalov et al., 2022) and neutrinos (e.g., Janka, 2017; Horiuchi & Kneller, 2018; Müller, 2019a). Current neutrino and GW detectors, however, are only sensitive to events within $\mathord{\sim}100$ kpc (Scholberg, 2012; Abbott et al., 2016). Electromagnetic signals are more readily available, especially in today’s era of large-scale surveys (Bellm, 2014; Chambers et al., 2016; Tonry et al., 2018; Masci et al., 2019; Ivezić et al., 2019).

Type IIP supernovae (SNe IIP) are of particular interest for comparing theoretical CCSN model predictions to observations. They are the most common observed supernova type and originate from hydrogen-rich red supergiants (RSGs; Smartt 2015) that are predominantly, though not exclusively, unaffected by binary mass transfer (Podsiadlowski et al., 1992; Zapartas et al., 2021). They thus represent the supernova sub-population that most closely matches the progenitor models underlying population studies based on phenomenological explosion models²²2Note that phenomenological explosions models for stripped stars in binary systems have also been presented recently by Ertl et al. (2020) and Schneider et al. (2021).. A Type IIP supernova exhibits a $\sim$100-day phase with nearly constant luminosity (“plateau” – P) in its light curve (LC) during the inward propagation of a recombination wave through the shock-heated hydrogen envelope. The plateau luminosity ($L_{\rm pl}$) and duration ($t_{\rm pl}$) are related to the CCSN explosion energy, the progenitor radius, and the mass of the hydrogen envelope (Popov, 1993; Kasen & Woosley, 2009). The plateau phase is followed by an exponential luminosity tail that is at first powered by the radioactive decays of ⁵⁶Ni and ⁵⁶Co and by other radioactive species later on.

Historically, supernova explosion and progenitor properties have most often been inferred by fitting individual SN LCs with semi-analytic solutions (e.g., Arnett, 1980, 1982) or radiation hydrodynamic simulations (e.g., Blinnikov et al., 2000; Kasen & Woosley, 2009; Bersten et al., 2011; Dessart & Hillier, 2010). This approach, however, can suffer from a degeneracy of the explosion energy and progenitor mass as key parameters that determine the LC (Dessart & Hillier, 2019). The problem of degeneracies can be reduced by considering larger samples of observed transients from surveys or compilations (Li et al., 2011; Faran et al., 2014; Pejcha & Prieto, 2015a; Martinez et al., 2022; Gutiérrez et al., 2017; Müller et al., 2017b; Martinez et al., 2020, 2022). Most work on inferring progenitor and explosion parameters for larger supernova samples to date has relied on LC fitting, i.e., on reverse modeling (Morozova et al., 2018; Martinez et al., 2020). The complementary approach is to use forward modeling of entire supernova populations for validating or constraining CCSN explosion models. Several recent studies produced a considerable number of LCs derived from phenomenological CCSN explosion models (Sukhbold et al., 2016; Barker et al., 2022a; Curtis et al., 2021). What is still missing, however, is a global comparison between such a suite of theoretical models and a representative, volume-limited supernova sample.

Such a comparison also needs to explore the sensitivity and robustness of explosion parameter and LC predictions to variations in model assumptions. This is particularly important to ascertain the potential for determining physical parameters of individual supernovae or entire populations, e.g., the recent idea to exploit proposed correlations between iron core mass and plateau luminosity (Barker et al., 2022a, b) for use in parameter inference.

In this work, we use the radiation hydrodynamics code SNEC (Morozova et al., 2015b) to calculate LCs of SNe IIP based on two sets of progenitor and explosion models from Müller et al. (2016, hereafter M16) and Sukhbold et al. (2016, hereafter S16). We obtain explosion properties using the efficient semi-analytic model for neutrino-driven explosions from M16. Though both evolved with the stellar evolution code KEPLER (Weaver et al., 1978; Heger & Woosley, 2010), M16 and S16 progenitors have been evolved with slightly different physics assumptions, and illustrate that LC predictions are especially sensitive to model variations that affect the hydrogen envelope. We then quantitatively compare the models to the volume-limited SNe IIP sample of Pejcha & Prieto (2015a, hereafter PP15) to highlight salient points of agreement and disagreement between the predictions and the observations. In particular, we highlight that even though the models reproduce the well-known correlation between plateau luminosity $L_{\rm pl}$ and nickel mass $M_{\rm Ni}$, there are still tensions between models and observations in distribution of plateau luminosity, plateau duration and nickel mass. Similar to Dessart et al. (2013), our results underscore the sensitivity of the LCs to the envelope structure of the progenitor.

Our study has not taken into account the following uncertainties encountered in the theoretical modelling of SNe IIP LCs: the mixing-length parameter and convective overshooting (Maeder & Meynet, 1987; Dessart et al., 2013), mixing and composition (Couch et al., 2015; Dessart & Hillier, 2020), clumping (Kifonidis et al., 2000; Dessart et al., 2018; Dessart & Audit, 2019), large scale asymmetry (Wongwathanarat et al., 2015; Dessart et al., 2021), line blanketing (Kasen & Woosley, 2009; Dessart & Hillier, 2011), and metallicity (Dessart et al., 2014). Also our models do not consider the impact of circumstellar materials on SNe IIP LCs (Chugai et al., 2007; Dessart & Audit, 2019; Dessart & Hillier, 2022). How these complexities affect the SNe IIP population warrants further exploration.

This paper is organized as follows. In §2 we review the semi-analytic approach for obtaining neutrino-driven CCSN explosions and the resulting explosion landscape for supernova progenitors. We compare our approach to two other phenomenological explosion models for the S16 progenitor set in §3. In §4 we present the theoretical SN IIP LCs from radiation hydrodynamic simulations, and they are compared to the observational sample in §5. Our conclusions are given in §6.

In this section, we first review the semi-analytic approach of Müller et al. (2016) for neutrino-driven CCSN explosions. Then we apply this semi-analytic model to obtain the properties of CCSN explosions for two sets of progenitor models.

It is currently not feasible to perform 3D simulations with neutrino transport to determine the properties of CCSN explosions for a sufficiently large number of progenitors required for population studies. Our semi-analytic model is among several efficient phenomenological approaches to predict the outcome of collapse (explosion or non-explosion) as well as explosion and remnant properties (O’Connor & Ott, 2011; Ugliano et al., 2012; Pejcha & Thompson, 2015; Perego et al., 2015; Sukhbold et al., 2016; Couch et al., 2020; Ertl et al., 2020; Barker et al., 2022a; Ghosh et al., 2022). Most other studies rely on 1D simulations that mimic the supportive role of multi-dimensional flow instabilities in enabling shock revival either by increasing the neutrino emission, the neutrino energy deposition, or by means of 1D turbulence models (but see Müller 2019b for a critical discussion of this approach). Qualitative and quantitative differences and similarities between the various phenomenological models have been discussed in the literature, and Pejcha (2020) also provides a side-by-side comparison of important outcomes such as the relation between explosion energy and nickel mass or the predicted neutron star mass distribution. Such comparisons can be somewhat skewed by differences in the size, mass range, and input physics of underlying stellar evolution model sets.

For this reason, it is useful to compare our results to those obtained by different 1D simulation studies for the S16 progenitor set, namely from the study of Sukhbold et al. (Sukhbold et al., 2016) and Barker et al. (Barker et al., 2022a). Sukhbold et al. used the P-HOTB code (Ugliano et al., 2012; Ertl et al., 2016) with a gray neutrino-transport scheme and a proto-neutron star core model, and is calibrated by two well-observed CCSNe. Their models are calibrated to inferred explosion properties for SN 1054 and SN 1987A at the respective progenitor masses. The SN 1987A calibration is used for all progenitors with $M_{\rm ZAMS}>12\,M_{\odot}$, and for $M_{\rm ZAMS}<12\,M_{\odot}$ interpolation between the relevant model parameters for the two calibration cases is applied. Barker et al. used the FLASH code with a multi-group two-moment neutrino-transport scheme (O’Connor & Couch, 2018) plus the STIR method for simulating turbulence in 1D (Couch et al., 2020). Their STIR method is calibrated to fit full 3D simulations run in the same code (O’Connor & Couch, 2018).

The comparison is shown in Fig. 4 for $E_{\rm exp}$, $M_{\rm Ni}$ and $M_{\rm NS,by}$. Although with quite different implementations and degrees of approximations, we find considerable agreements among the results from Sukhbold et al., Barker et al. and this work. The agreement is especially remarkable for the baryonic neutron star mass $M_{\rm NS,by}$, which once again confirms the important role of the Si-O shell interface as a natural point for the onset of the explosion and a strong predictor for the final neutron star mass.

Discrepancies are noteworthy mainly in the mass ranges with near-critical explodability (gray shaded bands in Fig. 4), with $M_{\rm ZAMS}\simeq$ $12$–15 $M_{\odot}$ and $22$–25 $M_{\odot}$. For $M_{\rm ZAMS}\simeq$ $12$–15 $M_{\odot}$, Barker et al. predicts no explosion while both Sukhbold et al. and the semi-analytic model obtain explosions. $E_{\rm exp}$ and $M_{\rm Ni}$ in Sukhbold et al. are, however, larger by about a factor of 2.5 than those in this work, which may be related to the change in calibration case of P-HOTB from SN 1054 to SN 1987A at $12\,M_{\odot}$.

On the other hand, for $M_{\rm ZAMS}\simeq$ $22$–25 $M_{\odot}$, Sukhbold et al. and our semi-analytic model predict no explosion, whereas Barker et al. yields relatively large explosion energies $E_{\rm exp}$ ($\geq 2\times 10^{51}$ erg). The explodability of these critical models is still under debate with state-of-the-art 3D simulations (e.g., Ott et al., 2018; Melson et al., 2020; Burrows et al., 2020). The mass distribution of observed SN IIP progenitors (Smartt, 2015) and first observational evidence for the quiet disappearance of a RSG (Adams et al., 2017), presumably by stellar collapse favor a lower probability of explosion in this mass range⁴⁴4A potential exception is SN 2015bs whose progenitor was inferred to have a metallicity of $\leq 0.1Z_{\odot}$ and a $M_{\rm ZAMS}$ of 17-25 $M_{\odot}$ (Anderson et al., 2018)..

The overall trends and patterns in explosion energy are qualitatively compatible between the three phenomenological models outside the gray-shaded areas. They all predict low explosion energies at the low-mass end, a general trend towards higher explosion energies in the range of $15$–22 $M_{\odot}$ with considerable scatter at higher masses. Above 25 $M_{\odot}$, the agreement is less convincing. It is noteworthy, however, that even in the region of $22$–25 $M_{\odot}$, where Barker et al. disagrees qualitatively with the other two models, the high explosion energies reflect a similar pattern in Müller et al. (2016) with parameter choices that increase explodability (i.e., higher turbulent pressure in the gain region or a higher accretion efficiency for neutrino emission).

The situation for the nickel masses, $M_{\mathrm{Ni}}$, which are only available for Sukhbold et al. and our semi-analytic model, is similar to the explosion energies. There is rather good agreement between Sukhbold et al. and our work below $22\,M_{\odot}$, which is rather striking considering the relatively simple model for nickel production used in our approach.

These results demonstrate that predictions of explosion and remnant properties from the three phenomenological models are quite robust to differences in the methodology, once some form of calibration (e.g., for one or two specific supernovae or for the typical energy range of observed explosions) is applied.

With the explosion properties ($E_{\rm exp}$, $M_{\rm Ni}$ and $M_{\rm NS,by}$) obtained in §2, we utilize SNEC (Morozova et al., 2015b) to generate LCs of SNe IIP from M16 and S16 progenitors. SNEC is an open-source spherically-symmetrical radiation hydrodynamics code with the capability to follow the shock propagation through the stellar envelope. It solves the Lagrangian hydrodynamics equations supplemented with a radiation diffusion term. Note that SNEC assumes local thermal equilibrium between matter and radiation, which fails during the shock breakout and nebular phase, but is reasonably reliable for LCs during the plateau phase (Blinnikov & Bartunov, 1993) that is of interest here. We refer to the code paper (Morozova et al., 2015b) and documentation (Morozova et al., 2015a) for details on the numerical implementation.

We employ the default settings of SNEC, such as the equation of state, ionization treatment and opacities. The newborn NS with $M_{\rm NS,by}$ is excised from the numerical grid and a thermal bomb is used to initialize the shock. The sum of $E_{\rm exp}$ and binding energy of the mass content above the excised NS is spread into the $0.1\,M_{\odot}$ above the excised boundary so that the final explosion energy equals the desired value $E_{\rm exp}$ (Morozova et al., 2015a). For the mixing of nickel, we simply spread $M_{\rm Ni}$ homogeneously up to $3\,M_{\odot}$ as our semi-analytic approach cannot treat the mixing. The mixing of nickel is beyond the scope of this paper but its impact on SNe IIP LCs may be worth further investigation (see, e.g., Utrobin et al., 2017). We evolve all the models to $\sim 200$ days, by which time all models have reached the radioactively-powered tail phase. For comparison to observations, we are particularly interested in two LC parameters: the plateau luminosity $L_{\rm pl}$ and the plateau duration $t_{\rm pl}$. We take the bolometric luminosity at 50 days after the shock break out as $L_{\rm pl}$. The determination of $t_{\rm pl}$ is more tricky; we tentatively pick the time of the steepest gradient of the $B$-band magnitude as the end of plateau phase. We also present the photospheric velocity at 50 days $v_{50}$ after the shock breakout, which is a proxy for the mean expansion rate of the ejecta and $E_{\rm exp}/M_{\rm ejecta}$, with $M_{\rm ejecta}$ being the total mass of the ejecta. Here we use the SNEC definition for the location of photosphere, i.e. by the optical depth $\tau=2/3$. The key LC and explosion parameters for all models are publicly available at Zenodo: doi:10.5281/zenodo.7354733 in the same form as listed in Table 1 .

As representative examples, we plot in Fig. 5 the bolometric LCs of SNe IIP from the pre-SN models shown in Fig. 1, with their respective $L_{\mathrm{pl}}$ and $t_{\mathrm{pl}}$ given in Table 1. It is clear at a first glance that the M16 models are brighter than S16 models during the plateau phase for the same $M_{\rm ZAMS}$, despite the similar explosion properties (also listed in Table 1). This feature is further exemplified in Fig. 6, which compares $L_{\rm pl}$, $t_{\rm pl}$ and $v_{50}$ as a function of $M_{\rm ZAMS}$ between all M16 and S16 models that successfully explode. Whereas $t_{\rm pl}$ is quite similar for the two sets of models, $L_{\rm pl}$ of M16 models is in general larger by a factor of $\mathord{\sim}2$ than that of S16 models. This difference cannot be accounted for even by appealing to large uncertainties in the explosion energy. Similar values of $L_{\rm pl}$ as in S16 can only be realized for M16 models by artificially dividing $E_{\rm exp}$ by three, which is unrealistic and would affect $t_{\rm pl}$ considerably. Indeed, the difference in $L_{\rm pl}$ reflects the systematically different envelope structure between M16 and S16 progenitors (see the density profiles in Fig. 1 and the pre-SN masses and radii in Fig. 2). The slightly larger $v_{50}$ for M16 models reflects the smaller $M_{\rm prog}$ of the M16 progenitors which leads to a smaller $M_{\rm ejecta}$ (cf. panel (a) in Fig. 2). As we shall see in §5, the comparison with observations suggests a preference for the S16 models as realistic progenitors as they match the observed plateau luminosities better.

Lastly, we compare our results to analytic scaling relations often used by observers to infer the properties of progenitor and explosion from LC parameters, both to guide the interpretation of our results and to check the validity of the analytic relations. For $L_{\rm pl}$, we use the relation derived in Popov (1993)

where $E_{51}$ is the explosion energy in units of $10^{51}$ erg, $M_{10}$ is the mass of the hydrogen envelope (the progenitor mass minus the helium core mass) in units of 10 $M_{\odot}$, and $R_{\rm 0,500}$ is the pre-SN stellar radius $R_{\mathrm{prog}}$ in units of $500\,R_{\odot}$. Our preferred values of $L_{0}$ are $1.69\times 10^{42}\,$erg s^-1 and $1.51\times 10^{42}\,$erg s^-1 for M16 and S16 models, respectively. The left panel of Fig. 7 shows that Eq. (3) predicts $L_{\rm pl}$ well overall, with a relative error $\lesssim 10\%$ for most models. The discrepancy for models with a large $L_{\rm pl}$ with a relative error up to $40\%$ is due to their short plateau for which $L_{\rm bol}$ at 50 days may not well represent $L_{\rm pl}$.

The scaling relation for the plateau duration from Popov (1993) assumes no energy input from radioactive decay of nickel and cobalt and reads

Following Sukhbold et al. (2016), we use a modified relation for $t_{\rm pl}$ that takes into account that energy input from radioactive decay can prolong the plateau,

where we set the constant $C_{f}=21$ as suggested in Sukhbold et al. (2016). Comparing the LCs from SNEC to Eq. (5) is more appropriate, as SNEC includes the energy release from radioactive decay. The fitted $t_{0}$ are $93.0\,\mathrm{d}$ and $89.7\,\mathrm{d}$ for M16 and S16 models, respectively. The right panel of Fig. 7 shows that Eq. (5) predicts $t_{\rm pl}$ well at $t_{\rm pl}\gtrsim 100$ days, with a relative error $\lesssim 15\%$. For $t_{\rm pl}\lesssim 100$ days, the relative error can be up to $\sim 25\%$.

In this paper, we presented $\sim 1100$ light curves of SNe Type IIP generated by SNEC from two sets of single-star solar-metallicity progenitor models in M16 (Müller et al., 2016) and S16 (Sukhbold et al., 2016), with very high resolution in ZAMS mass grid as fine as $0.01\,M_{\odot}$ in the former set. We assume that SNe IIP are driven by neutrinos and calculate the key explosion parameters $E_{\rm exp}$, $M_{\rm NS,by}$ and $M_{\rm Ni}$ using a semi-analytical approach derived in Müller et al. (2016).

The explosion parameters agree well globally between the M16 and S16 model sets and between the semi-analytic model and alternative phenomenological explosion models from previous studies of exploding S16 models (Sukhbold et al., 2016; Barker et al., 2022a). In particular, the agreement between the prediction of the semi-analytic model and the 1D simulations of Sukhbold et al. (2016) for the same progenitor set is striking. The plateaus of SNe Type IIP are systematically fainter by a factor of $\mathord{\sim}2$ in bolometric luminosity for the S16 set due to denser hydrogen envelopes of S16 progenitors. The more extended envelope structure of the M16 models lead to brighter plateaus and is likely artificial because of simplification of the surface boundary condition in the stellar evolution calculations. This reinforces previous findings on the sensitivity of Type IIP explosions to the envelope structure (Dessart et al., 2013) and implies that difference in theoretical light curves may rather reflect assumptions about stellar structure and evolution, in particular those that affect the structure of the convective RSG envelope, than the modeling of the explosion engine. As already pointed out by Dessart & Hillier (2019), this may cause problems in inferring progenitor properties from observables, e.g., inferring the ZAMS mass from the plateau luminosity (Barker et al., 2022b). It is important to highlight that even among available stellar evolution models computed with the same code, there may be subtle different in the treatment of the convective envelope and outer boundary due to code improvements and model parameter choices that may have significant repercussions for supernova light curve modeling. To fully exploit the diagnostic potential of SNe Type IIP light curves, more theoretical and observational work on RSG envelopes and environments is critical.

We compare the parameters of the predicted light curves to the volume-limited PP15 sample of well-observed SNe IIP (Pejcha & Prieto, 2015a). We construct a mock supernova population from the two progenitor sets by weighting the models with the Salpeter IMF. Based on the mean value of the plateau luminosity $L_{\rm pl}$ and a K-S test, the S16 models fit the observed brightness distribution in the PP15 SNe IIP sample better. We find a similar correlation between $L_{\rm pl}$ and $M_{\rm Ni}$ in both model sets and in the PP15 sample. However, we find tensions with the observational data for both model sets, which may either indicate an incomplete understanding of the progenitor-explosion connection or of the pre-supernova progenitor structure. Both progenitor sets lack models with explosions that produce the very small nickel masses $M_{\rm Ni}<0.01\,M_{\odot}$ that are observed in some IIP explosions. This discrepancy may be related to the uncertainties in models at the low-mass end (Woosley & Heger, 2015b) or to models close to black-hole formation. The comparison of the plateau duration $t_{\rm pl}$ remains beset with ambiguities in the definition of the plateau length, but we tentatively find a significantly larger spread in plateau duration in the models compared to the observed sample. Furthermore, the models predict an anti-correlation between $L_{\rm pl}$ and $t_{\rm pl}$, which is not found in the PP15 sample. Indications of an anti-correlation have, however, been found in other samples (Faran et al., 2014), and future studies need to assess whether bigger volume-limited samples confirm the tension between theory and observations.

These results provide an interesting lead for further comparisons of the theoretical models and observational data for SNe IIP LCs. In particular, the predicted correlation between plateau luminosity and plateau duration would present a challenge to current stellar evolution models of massive stars, if the tension to observations can be corroborated using larger volume-limited transient samples. Future studies should explore variations of single-star and binary evolution models and pit them against bigger volume-limited transient samples from recent and upcoming surveys. By obviating the need for time-critical 1D (let alone 3D) supernova simulations, our semi-analytic model may be useful for conducting such large-scale comparisons more efficiently with little loss of accuracy, given the remarkable agreement with the explosion properties obtained by Sukhbold et al. (2016).