1. Introduction
Stellar death-throes in supernovae are some of the most spectacular events that stars produce in the Universe. Attempting to understand these events is both a mature and constantly developing field. It is primarily driven by observations. New discoveries occur primarily when new observational techniques are employed, or when the peculiarities of an individual event attract attention and follow-up, leading to discovery of a new transient type. In comparison, theory may lag behind. Some rare transient types have been theoretically predicted before being observed, such as the kilonovae from the merger of two neutron stars, which appeared as theoretical predictions (Li & Paczyński Reference Li and Paczy´nski1998) and a tentative detection associated with a short gamma ray burst (Tanvir et al. Reference Tanvir, Levan, Fruchter, Hjorth, Hounsell, Wiersema and Tunnicliffe2013 before the electromagnetic radiation associated with GW 170817 unambiguously confirmed their existence (Abbott et al. 2017. However these have been deduced primarily from gravitational theory (Kilonovae, Tidal Disruption Events, etc). There has not yet been a clear case of stellar structure and evolution theory leading the way and predicting new supernovae types that were then observed.
The success of such theoretical predictions, verified by observation, in other areas indicates that we should redouble our theoretical work in understanding supernovae. In Eldridge et al. (2018, hereafter paper I), we introduced the supernova lightcurve population synthesis (CURVEPOPS) project with this goal. This comprises lightcurves derived from a large number of supernova progenitor models identified within the Binary Population and Spectral Synthesis project (BPASS, Eldridge et al. Reference Eldridge, Stanway, Xiao, McClelland, Taylor, Ng, Greis and Bray2017) and exploded with the Supernova Explosion Code (SNEC, Morozova et al. Reference Morozova, Piro, Renzo, Ott, Clausen, Couch, Ellis and Roberts2015). By studying the synthetic lightcurves from a population of realistic progenitors drawn from populations including binary stars we were able to demonstrate that binary interactions are the main source of diversity of type II supernova lightcurves.
Here we take the next step which is to validate the CURVEPOPS models against observational data, and hence gain insight into the importance of explosion parameters such as explosion energy, nickel mass and amount of nickel mixing. The most effective way to perform this test is to consider supernovae for which observations of both the progenitor stars and the observed supernovae exist in the archive. Such studies have been carried out before for individual supernovae and large samples (e.g. Bersten et al. Reference Bersten2012, Reference Bersten2014; Morozova et al. Reference Morozova, Piro and Valenti2017, Reference Morozova, Piro and Valenti2018). Here we have focused our study on the sample of type IIP supernova progenitors that have detected progenitor stars as described in Smartt (Reference Smartt2015). We have collated lightcurve data from the Open Supernova Catalogue Footnote a and compared it to a large suite of SNEC explosion models derived from single-star progenitor stellar structures calculated as part of the BPASS project. We then constrain the nature of the progenitor star from fitting the supernova lightcurve and compare these inferences with those in the literature based on fitting the colour and magnitude of the progenitor star in pre-explosion images. This allows us to gain insight into the accuracy and usefulness of our CURVEPOPS models. The SNEC inputs and outputs from this project have been made available at the BPASS website (http://bpass.auckland.ac.nz) and in the PASA datastore.
We note that there are many studies that have investigated the lightcurves of core-collapse supernovae (e.g. Utrobin Reference Utrobin1994, Reference Utrobin2005, Reference Utrobin2007; Utrobin & Chugai Reference Utrobin and Chugai2008, Reference Utrobin and Chugai2009; Dessart et al. Reference Dessart, Livne and Waldman2010, Reference Dessart, Hillier, Livne, Yoon, Woosley, Waldman and Langer2011; Bersten et al. Reference Bersten, Benvenuto and Hamuy2011, Reference Bersten2012; Dessart et al. Reference Dessart2014; Bersten et al. Reference Bersten2014; Morozova et al. Reference Morozova, Piro, Renzo, Ott, Clausen, Couch, Ellis and Roberts2015; Dessart et al. Reference Dessart, Hillier, Woosley, Livne, Waldman, Yoon and Langer2016; Morozova et al. Reference Morozova, Piro, Renzo and Ott2016, Reference Morozova, Piro and Valenti2017, Reference Morozova, Piro and Valenti2018; Martínez & Bersten Reference Martínez and Bersten2018, Reference Martínez and Bersten2019). However what sets this study apart is the large number of SN models that we have created, aiming to model a number of supernovae simultaneously rather than attempting to produce a perfect fit for one SN alone. Furthermore we are attempting to determine how strong a link there is between the time evolution of the observed explosion and the properties of the progenitor star in a lightcurve-derived model. Our hope is that from this work we will be able to improve our CURVEPOPS lightcurves and use relations derived from this work in terms of final mass to explosion energy, nickel mass and nickel mixing to produce more realistic supernova lightcurve populations.
The structure of this paper is as follows, first we describe the creation of our grid of supernova models with BPASS and SNEC. Next we describe the observational sample of supernovae we employ in this project. We then outline the fitting method used to compare the models to observations. We present and discuss our results, before summarizing our conclusions.
2. Creation of supernova simulations
We use the v2 single-star models from the Binary Population and Spectral Synthesis code BPASS (Eldridge et al. Reference Eldridge, Stanway, Xiao, McClelland, Taylor, Ng, Greis and Bray2017). The models are calculated using a version of the Cambridge STARS code that has been adapted to follow binary evolution (see Eldridge et al. Reference Eldridge, Stanway, Xiao, McClelland, Taylor, Ng, Greis and Bray2017, for full details). The results of the code are available from the website, bpass.auckland.ac.nz. We summarize the most important details here.
The models are calculated from the zero-age main-sequence up to the end of core carbon burning. At this point the models stop since the STARS code is unable to compute further due to the increasing complexity and time resolution of late time evolution. We assume that our models are close enough to the time of core-collapse that the parameters we use for our lightcurve models will not vary significantly. We have used a model calculated from the Modules for Experiments in Stellar Astrophysics (MESA-r10398, Paxton et al. Reference Paxton, Bildsten, Dotter, Herwig, Lesaffre and Timmes2011, Reference Paxton2013, Reference Paxton2015, Reference Paxton2018) stellar evolution code for a typical SN IIP progenitor – a single-star with M = 15.6 Mʘ – to explore how different later models might be if evolution is taken closer to core-collapse. The results of this analysis are presented in Appendix C. In brief, the outer core structure changes little while the density of the central core increases; the greatest changes are within the region that is assumed to form the remnant. This may change the late time evolution of our supernova models and so our late-time lightcurves should be assumed to be subject to increased uncertainty. If we compare our models to other supernova models (e.g. Dessart et al. Reference Dessart, Hillier, Waldman and Livne2013; Dessart & Hillier Reference Dessart and Hillier2019; Goldberg et al. Reference Goldberg, Bildsten and Paxton2019) we find that they are generally similar in shape. Differences are greatest at the end of the plateau phase where the extra structure and mixing within the model will have the greatest impact on the resultant lightcurve.
We use models with a metallicity mass fraction of Z = 0.014 which is close to estimates of massive stars in the Solar neighborhood (Nieva & Przybilla Reference Nieva and Przybilla2012). We use every integer initial mass model from 5 to 26 Mʘ as these all have sufficiently massive hydrogen envelopes to give rise to a long lasting plateau in their lightcurves. We do not use the binary models because, as shown in paper I, type IIP supernovae mostly arise from progenitors that have not experienced significant binary interactions. While some stars in binaries may experience interactions and provide type IIP progenitors with interior structures different to single stars, if we were to explode all the binary progenitors from paper I with the range of explosion parameters below this would require calculation of 154 791 supernova simulations. Doing so is beyond the scope of the current pilot study.
We next format the stellar structure models, taken at their last time step, for input into the SuperNova Explosion Code, SNEC. This is an open-source code that is available online from https://stellarcollapse.org/SNEC and which has been applied to modelling various aspects of supernovae (e.g. Morozova et al. Reference Morozova, Piro, Renzo, Ott, Clausen, Couch, Ellis and Roberts2015, Reference Morozova, Piro, Renzo and Ott2016, Reference Morozova, Piro and Valenti2017, Reference Morozova, Piro and Valenti2018). We only include the composition variables that we have within the BPASS stellar evolution code which are hydrogen, helium, carbon, nitrogen, oxygen, neon, magnesium, silicon, and iron. We also include a nickel-56 variable but leave this blank, inputting the nickel within the explosion parameters of the SNEC model. We have made all these input files available on the PASA datastore as well as on the BPASS website (http://bpass.auckland.ac.nz).
We extend the surface layers of the progenitor models to include the stellar progenitor’s wind. Recent studies have shown how including the red supergiant’s wind can have important impact on the early lightcurve of the progenitor star (Morozova et al. Reference Morozova, Piro, Renzo and Ott2016, Reference Morozova, Piro and Valenti2017, Reference Morozova, Piro and Valenti2018; Moriya et al. Reference Moriya, Förster, Yoon, Gräfener and Blinnikov2018). To incorporate this into our progenitor models we determine the mass-loss rate from the stellar model and calculate the wind velocity using the method outlined in Eldridge et al. (Reference Eldridge, Genet, Daigne and Mochkovitch2006). The mass-loss rates used are those of de Jager et al. (Reference de Jager, Nieuwenhuijzen and van der Hucht1988). We do not vary the wind parameters as in Moriya et al. (Reference Moriya, Förster, Yoon, Gräfener and Blinnikov2018) but leave them fixed so that the nature of the circumstellar medium is linked to the initial mass of the progenitor. We attach the outer most stellar mesh point to the wind assuming a beta wind velocity law as used by Moriya et al. (Reference Moriya, Förster, Yoon, Gräfener and Blinnikov2018) with β = 5 which is typical of cool red supergiants (Schroeder Reference Schroeder1985; Gräfener & Vink Reference Gräfener and Vink2016) although we note this value might not be correct for all red supergiants (Ohnaka et al. Reference Ohnaka, Weigelt and Hofmann2017). Varying this parameter will make small changes to the very early lightcurve and future observations, with better sampling of the early time evolution, will be required to determine the best value to use. We have run supernova models with values of β between 2 and 6, and these are shown in Appendix D. We find that the effects of varying β are somewhat degenerate with the total density of the circumstellar material density assumed. The effect of the material however is restricted to the early lightcurve for most stellar masses, with an effect on the late time plateau duration for only the most massive stars. Given that most of the progenitors are expected to have initial masses below about 20 Mʘ, at which the full range of β explored changes the plateau length by <10 d, we do not expect a significant impact on our results. Nonetheless it is clear that this area should be considered in future studies.
We extend the stellar wind out to approximately 10 000 times the radius of the star which is typically shorter than a parsec. We show examples of these lightcurves compared to those without the stellar wind included in Figure 1. We find inclusion of this causes brighter phases in our early lightcurves. This is most significant for the most massive stars above 20 Mʘ. There are also some spurious jumps in the nickel tail of our model as varying shells of material collide at late times. We suggest that this is an artifact of our one-dimensional simulations and that this behaviour would be smoothed out in a 3-dimensions.
Figure 1. Sample synthetic lightcurves, demonstrating how including the circumstellar material of the red supergiant’s wind changes the lightcurve. The solid lines are for lightcurve models with the circumstellar material included as discussed in the text and the dotted lines assume no material surrounding the progenitor star. The dashed line is where the circumstellar material density has been reduced by a factor of 2. While the dash-dotted line is where the circumstellar material density has been increased by a factor of two. The explosion energy, nickel mass and mixing are kept constant (1050.5 erg s−1, 10−1.5 Mʘ and mid-strength mixing). Figures for each initial mass on its own are included in Appendix D.
The mass-loss rates of red supergiants are uncertain as there is a question as to whether the rates of de Jager et al. (Reference de Jager, Nieuwenhuijzen and van der Hucht1988) are correct or not. Two studies, Mauron & Josselin (Reference Mauron and Josselin2011) and Beasor & Davies (Reference Beasor and Davies2018), found that while de Jager et al. (Reference de Jager, Nieuwenhuijzen and van der Hucht1988) have the correct scale of mass-loss rates there are significant departures away from the predicted rates for some stars. Mauron & Josselin (Reference Mauron and Josselin2011) found that red supergiants varied by up to a factor of four away from the de Jager et al. (Reference de Jager, Nieuwenhuijzen and van der Hucht1988) predictions. While Beasor & Davies (Reference Beasor and Davies2018) found more significant deviations, especially that de Jager et al. (Reference de Jager, Nieuwenhuijzen and van der Hucht1988) rates may be overestimating the mass lost during the red supergiant phase significantly (we note in their study that the STARS models, effectively identical to BPASS models, are closest to their estimated mass losses of all the stellar evolution models they consider). In addition to these studies there is evidence from radio observations of supernovae that the mass-loss rates we assume are similar to those of typical red supergiants (Chevalier et al. Reference Chevalier, Fransson and Nymark2006). In light of this uncertainty, in addition to varying the β parameter, we explore this further by testing how varying the circumstellar material density affects our model lightcurves. We show models with double or half the wind density of our fiducial models in Figure 1 and Appendix D. Again we find that, for the majority of the lightcurves, varying the wind density by a factor of two causes minimal changes to the resulting light curves. As noted before, the effects are also in degenerate with the assumed value of the wind acceleration parameter β.
The largest impact of varying the circumstellar density is again observed in the early lightcurves, and particularly for the most massive stars. For many of the observed lightcurves available in the literature, the early evolution is not well sampled or is simply unobserved, leaving models largely unconstrained. The uncertainties in the winds of RSGs, while important to consider, should thus not alter our interpretations to any significant degree. We also note that the most massive stars leading to core collapse supernovae, for which the wind effects are greatest, are also more likely to have massive binary partners and interact, leading to stripped-envelope supernovae. By not accounting for binary star progenitors in these models we are also ignoring variations in the circumstellar medium that would be caused by binary interactions which may increase or decrease apparent mass-loss rates for a progenitor. Therefore rather than varying the circumstellar medium in our fitting grid based on mass, we link the assumed circumstellar medium density to the progenitor at the point of explosion in line with our current best understanding.
Within SNEC, as well as the progenitor structure, other parameters must be specified. These are the explosion energy, nickel mass, amount of nickel mixing and excised mass. In paper I we assumed that the first three parameters were constant while the excised mass was determined by the progenitor structure. Here we compute multiple models with varying explosion energy, nickel mass and nickel mixing. We still derive the excised mass to be the remnant mass computed as described in Eldridge et al. (Reference Eldridge, Stanway, Xiao, McClelland, Taylor, Ng, Greis and Bray2017). Our grid of values is as follows,
1. The explosion energy varies from log(ESN/erg s−1 = 50 to 52 in steps of 0.25 dex.
2. The nickel mass varies from log (MNi/Mʘ) = −3 to −1 in steps of 0.25 dex.
3. The nickel mixing is determined by setting the nickel boundary mass, out to which the nickel is mixed from the excised mass. We set this bound to one of three values: low, mid or max. These are determined as follows,
a. M Ni boundary low = M excised + 0.1M ejecta
b. M Ni boundary mid = M excised + 0.5M ejecta
c. M Ni boundary max = M excised + 0.9M ejecta
Covering all of the parameters of initial mass, explosion energy, nickel mass and nickel mixing requires 5346 SNEC models to be run. For each observed supernova we then use this grid of these models to find the parameters that best match its lightcurve.
We show in Figure 2 a sample of some of explosion models. Here we have shown how varying each of the explosion parameters varies certain aspects of the light curves. In these plots the base-line model is an initially 15 Mʘ progenitor, with an explosion energy of 1050.5 erg s−1, a nickel mass of 10−1.5 Mʘ and the mid strength nickel mixing.
Figure 2. Sample synthetic lightcurves, demonstrating how the explosion parameters change the lightcurve. The upper left panel shows how changing the initial mass of star varies the lightcurve when the stellar structure, explosion energy, nickel mass and mixing are kept constant (1050.5 erg s−1, 10−1.5 Mʘ and mid mixing). In the other panels one of the explosion parameters are varied while the stellar structure and other parameters are kept constant: upper right – explosion energy, lower left – nickel mass and lower right – nickel mixing.
In the panel with varying initial masses we see that most models have strong plateaus in their light curves. However, above 20 Mʘ we see early strong rises in the light curve due to denser winds around our progenitor stars. The plateau phase for these lightcurves is also extended due to the circumstellar medium. The exact shape of this rise will vary on the density of our assumed wind and also the value of β assumed for the wind acceleration as shown by Moriya et al. (Reference Moriya, Förster, Yoon, Gräfener and Blinnikov2018). For the lowest mass progenitors we see a late time rising of the light curves due to growing interaction between the supernova ejecta and our model circumstellar medium. In both case these features rely on our assumption that the surrounding wind is directly linked to the nature of the progenitor star as prescribed by de Jager et al. (Reference de Jager, Nieuwenhuijzen and van der Hucht1988) rather than allowing the mass-loss rates to vary arbitrarily. We stress that this initial work is largely exploratory and a full future study will involve models of varying metallicities, which will have different circumstellar environments, and also include progenitors that are the result of binary interactions. These will be very different to those from single-star evolution alone.
In the other panels we see first that varying the explosion energy changes the length and brightness of the plateau. Second, varying the nickel mass changes the late time evolution of the plateau as well as the brightness of the nickel tail. Finally changing the strength of nickel mixing has a complex impact: if the mixing is low then more nickel remains in the ejecta to power the late time light curve, while stronger mixing allows more of the nickel decay energy to escape without contributing to the V-band magnitude in the nickel tail.
We note that some light-curve tracks end abruptly. This is because of the known problem with SNEC that once a significant fraction of nickel is outside the photosphere the code stops predicting broad-band magnitudes for the supernova. There are also apparent bumps and short term features in some of the model light curves. This is the result of our circumstellar medium as well as the structure of our stellar progenitors. How realistic these features are is difficult to determine, but since most only occur late in the nickel tail we expect they will have only a small impact on our fitting here.
3. Observational sample
This work only considers Type II-P supernovae with observed progenitors as listed in Table 1 of Smartt (Reference Smartt2015). These are supernovae with direct progenitor detection and are all within a distance of approximately 20 Mpc. We list these in Table 1 along their host galaxy identification, assumed distance, foreground Galactic extinction for the host galaxy and the source of the photometric data we have used to reconstruct their lightcuves.
Table 1. List of Supernovae used in project
a The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
b Hendry et al. (Reference Hendry2005b)
c Maund et al. (Reference Maund, Reilly and Mattila2014)
d Fraser et al. (Reference Fraser2011)
e Maguire et al. (Reference Maguire2010)
f Tomasella et al. (Reference Tomasella2013)
g Fraser et al. (Reference Fraser2014)
h Crockett et al. (Reference Crockett, Smartt, Pastorello, Eldridge, Stephens, Maund and Mattila2011)
i Van Dyk et al. (Reference Van Dyk2012)
j Maund et al. (Reference Maund2013)
k Through Open Supernova Catalog: Guillochon et al. (Reference Guillochon, Parrent, Kelley and Margutti2017)
Distances were taken from cross-referencing literature, NEDFootnote b and SIMBAD (Wenger et al. Reference Wenger2000) values for the host galaxies. Where discrepancies exist values from the paper containing progenitor measurements were prioritized. Error calculations were performed using the uncertainties package of Python.
While there are some supernovae with observations from multiple filters in the source data sets, we have focused on the V band as there was plentiful data for each supernovae. In addition, U and B band are more strongly affected by extinction, and emission lines such as Hα may effect the R band magnitudes.
We show our sample lightcurves in Figure 3, after correction for distance and extinction. Each supernova in our sample has a clear plateau phase that lasts from approximately 75 to 125 d. There is a range of plateau luminosities of about two magnitudes and a broader range in the nickel tail of four magnitudes.
Figure 3. Light curves of observed supernova in V-band absolute magnitude. SN host details and sources of data are given in Table 1.
4. Lightcurve fitting
4.1. Preparing the synthetic lightcurves
All the synthetic lightcurves were visually inspected, since for some of the models the input stellar models have a small number of meshpoints that introduced oscillatory behaviour in the simulated lightcurves, where the luminosity of the lightcurve would vary by about a magnitude timestep by timestep, which was set to 6 h. Such oscillations would begin/disappear and the mean behaviour however was still that of a plateau lightcuve. We expect these were numerical in nature due to the sparse number of mesh points in some of our progenitor models, which causes numerical problems with the integration when the timestep was too high. To remove these spurious sequences we smoothed these parts of the light curve with the following simple algorithm:
1. Mask all model points where the neighbouring value differs from the previous by more than 0.1 magnitude.
2. Take the resulting array and mask again if the i + 2 value differs from i by more than 0.1 magnitude.
3. Compress the resulting array and repeat step 1.
The algorithm works as the time domain is very densely sampled and plotting of model light curves show that normally, consecutive model points should not differ from previous by more than 0.1 magnitude. The result removes most oscillations with a few that remained which were manually removed from our sample. These were mainly those with low explosion energies. Oscillatory behaviour also exists above 20 Mʘ, however the light curves at higher masses exhibit different shapes in the early plateau phase so it is unclear whether the behaviour is purely from model computation. Thus, smoothing was only applied to models up to 20 Mʘ.
4.2. Fitting procedure
As the first magnitude measurement of the SNe is not necessarily the explosion date, a search through literature was performed for all SNe to find the minimum and maximum possible explosion dates. The minimum possible date is taken as the date of SNe discovery, and the maximum is the last reported non-discovery. If the difference between minimum and maximum is larger than 10 d, then the explosion epochs tried are found by linearly splitting the range into 10 intervals and trying those 10 values. For SNe with no maximum possible explosion date, the range is set to 100 d. The best-fit date is then found and a second run testing best-fit date ±20 d in integer days is performed.
The model is linearly extended at the tail to the maximum time of the observed data points by using numpy.polyfit() on the last 15 model magnitude values to allow better fit of the radioactive-decay tail.
A minimum χ 2 method was employed to find the best fitting model. The photometric uncertainty associated with each observed data point is assumed to be Gaussian distributed around the measured value y obs, and the probability P for a match between the data and model is calculated as
\begin{equation} \chi^2 = \Sigma_{\rm mod} \left(\frac{y_{\mathrm{mod}}-y_{\mathrm{obs}}}{\sigma_{\mathrm{tot}}}\right)^2. \end{equation}
where σ tot is the quadrature-combined error from magnitude measurement and the error of the model, estimated at 0.25 mag as half of the magnitude difference between models of successive explosion energies. y mod is the model value at the same time t as the data measurement. Because the model is time-wise densely sampled, the model value at t is found by linear interpolation between model values using Python’s numpy.interp() function.
We then identify the model that has the minimum χ2 and estimate the uncertainty in each of the parameters as the range of values of each fit parameter for which Δχ2 < 5.89. Best fit and optimal χ2 values for each supernova are given in the appendix. The parameters we derive for each supernovae are the initial mass, nickel mass, explosion energy and amount of mixing. For the nickel mixing we introduced a numerical value, X, for the mixing of 0.1 for low mixing, 0.5 for intermediate mixing and 0.9 for max mixing. This thus represents the amount of nickel mixing into the ejecta mass as a fraction of the ejecta mass. We report these values in Tables 2 and 3 and include figures showing the distribution of χ2 over these parameters in the Appendix. When the uncertainties on our fits were less the spacing of our rather course grids in the parameters we assumed a minimum uncertainty for all our fit values taken from half the grid spacing of our models. For example, 0.13 (rounded up from 0.125) for our uncertainty in the explosion energy and 0.5 for the error in the initial mass.
Table 2. Reference and Free Fitted Parameters
a From progenitor observations and modelling, as in Smartt (Reference Smartt2015)
b Hendry et al. (Reference Hendry2005a)
c Hendry et al. (Reference Hendry2006)
d Smartt et al. (Reference Smartt, Eldridge, Crockett and Maund2009)
e Anderson et al. (Reference Anderson2014)
f Fraser et al. (Reference Fraser2011)
g Tomasella et al. (Reference Tomasella2013)
h Dall’Ora et al. (Reference Dall’Ora2014)
i Jerkstrand et al. (Reference Jerkstrand2015b)
j Bose et al. (Reference Bose2015)
Table 3. Reference and Progenitor Constrained Fitting Parameters
a From progenitor observations and modelling, as in Smartt (Reference Smartt2015)
b 0: low 0.9M r + 0.1M *, 1 : mid 0.5M r + 0.5M *, 2: max 0.1M r + 0.9M* where M r is mass of the remnant and M * is the mass of the ejecta.
c Hendry et al. (Reference Hendry2005a)
d Hendry et al. (Reference Hendry2006)
e Smartt et al. (Reference Smartt, Eldridge, Crockett and Maund2009)
f Anderson et al. (Reference Anderson2014)
g Fraser et al. (Reference Fraser2011)
h Tomasella et al. (Reference Tomasella2013)
i Dall’Ora et al. (Reference Dall’Ora2014)
j Jerkstrand et al. (Reference Jerkstrand2015b)
k Bose et al. (Reference Bose2015)
5. Results
Two sets of fitting were undertaken, one where the fitting code was free to run over the entire grid of models with masses from 6 to 26 Mʘ, and the other where the initial masses were constrained to vary over the ±1 σ mass range deduced from progenitor observations in Smartt (Reference Smartt2015). These two cases evaluate first what parameters can be derived from the lightcurve alone, and second whether a quantitative difference would be found in these fits in the rare cases where progenitor imaging provides additional constraints. We perform a quality assessment of the best fitting lightcurves (see Appendix) into broad categories, assigning a flag value: A – robust fit with 
$\chi^2 <\frac{1}{2}N_{\rm obs}$, B – good fit but with some features not captured and 
$\chi^2_{\rm min}\sim N_{\rm obs}$, or C - poor fit with 
$\chi^2_{\rm min}>N_{\rm obs}$. In the quantitative analysis that follows, poor quality class C fits are omitted from analysis, although they are shown in the figures for information. Only supernovae 2005cs and 2006my fall into this category due to poor end of plateau matching and a limited number of observed points for 2006my.
The derived parameters of the best-fitting lightcurve model for each supernova are shown in Tables 2 for fitting across the full, unconstrained mass range, and Table 3 for mass-constrained fits. The observational data used, together with the best fit model and the model closest to the fit parameters in each case are shown in Figures A11 and B1 in the appendices.
We again point out that our derived parameters are model dependent and have an inherent uncertainty that results from the assumptions put into SNEC. Most notably is the use of simple bolometric corrections to obtain broad-band magnitudes rather than the more complex method such as those in Dessart et al. (Reference Dessart, Hillier, Waldman and Livne2013). However checking against the results of other studies enables us to have confidence in the derived parameters even if the models have limitations. As an extra test we have compared the photosphere velocity given by SNEC in our best fitting models to the velocities of some spectral lines observed. We find that we predict values of the correct order although underpredict the velocities at early times. This is again complex as we do not calculate spectra for our models and so precise comparisons are difficult.
5.1. Validation against progenitor fitting
We test the validity and precision of our approach by comparing the progenitor parameters we derive from lightcurve fitting alone against those determined by Smartt (Reference Smartt2015) based on direct analysis of progenitors directed in pre-explosion imaging.
In Figures 4 and 5 we demonstrate the results of our two approaches – fits in which the progenitor parameters are unconstrained across our model grid, and fits in which the progenitor properties are permitted only to vary within the 1 σ uncertainty range associated with the Smartt (Reference Smartt2015) progenitor mass values. The former demonstrates the power of CURVEPOPS lightcurve fitting to yield an independent estimate of the parameters. The latter indicates the added value that lightcurve constraints can yield to reduce the range of uncertainties on parameters already estimated from progenitor observations.
Figure 4. A comparison of the progenitor initial mass derived from lightcurve fitting to that derived from analysis of progenitor observations. We show cases where the range of permitted values for lightcurve fitting is constrained by the range found by the progenitor fitting of Smartt (Reference Smartt2015) and where it unconstrained (i.e. where the two model fits are independent). Symbols and colours indicate reliable fits (classified A, blue squares), reasonable fits (B, green diamonds) and poor (C, yellow points). Small offsets are applied to models with the same integer masses for clarity.
Figure 5. A comparison of the progenitor Nickel-56 mass derived from lightcurve fitting to that derived from analysis of progenitor observations, as described in Figure 4.
As the figures demonstrate, the ability of CURVEPOPS to recover progenitor constraints without any dependence on pre-explosion imaging is impressive. Our unconstrained fits and those by Smartt (Reference Smartt2015) are entirely consistent, given the still rather large uncertainties on each parameter. There is a tendency for our more robust fits (classified A) to yield slightly higher progenitor initial masses (by an average of 2.8 Mʘ) although this is reduced (to 1.3 Mʘ) if class A and B fits are considered.
In Tables 2 and 3 we also compare our progenitors to the estimates from analysis of the surrounding stellar population from Maund (Reference Maund2017), the similar study using SNEC of Morozova et al. (Reference Morozova, Piro and Valenti2018) and updated progenitor masses from Davies & Beasor (Reference Davies and Beasor2018). We see that generally there is agreement consistent with the quoted uncertainties, however for certain supernovae there is some disagreement. For supernova 2004A the stellar population age is less than that which we derive and that inferred from the progenitor detection which could indicate that the progenitor was a runaway star and is not associated with the surrounding stellar population. The reverse is true for supernova 2009md where the stellar population inferred mass is too high but the light curve mass and pre-explosion mass match well.
Supernova 2012aw has the greatest mis-match as the Morozova et al. (Reference Morozova, Piro and Valenti2018) mass is significantly higher than that we estimate. Our fit does have regions that do not match the light curve despite the generally good (class A) fit and the agreement of our mass with the progenitor detections. For supernova 2012ec the reverse is true and we predict a higher initial mass than Morozova et al. (Reference Morozova, Piro and Valenti2018). We are uncertain why we achieve such a different fit. As the χ 2 parameter space indicates, there are a number of degeneracies in fitting to the light curves and the results can be dependent on the progenitor models.
A final verification of our fitting can be gained by comparing our masses to those estimated from late-time nebular spectra of some of our supernova sample by Jerkstrand et al. (Reference Jerkstrand2015a). The initial mass constraints they suggest are that SN 2006my was less than 12 Mʘ, SN 2012A was 12 Mʘ, while supernovae 2004et, 2012aw and 2012ec were all in the range between 12 to 15 Mʘ. Again these broadly agree with our derived masses and suggests an accurate future method to estimate progenitor masses.
It is also notable that when we constrain the permitted values to lie within the 68% confidence interval of the progenitor detection model, we are able to substantially reduce the uncertainty on both the initial and nickel masses, with our inferred uncertainties typically dominated by our sampling of these parameters rather than by the lightcurve data.
5.2. Additional explosion parameters
In addition to the progenitor parameters, CURVEPOPS analysis involves fitting over explosion parameters which are not explored in Smartt (Reference Smartt2015). These represent added value from this approach, which cannot be derived from pre-explosion progenitor imaging.
5.2.1. Explosion energy
In Figure 6 we investigate the dependence of assumed supernova explosion energy on the progenitor initial mass for the best fitting CURVEPOPS model.
Figure 6. The explosion energy derived from lightcurve fitting and its dependence on initial mass, with symbols as described in Figure 4.
We find no strong dependence of energy on mass, with a mean log(E exp/ergs) = 50.52 ± 0.10 for our most robust, class A lightcurve fits. When our slightly less-good, class B fits are included, we see some evidence for an increased spread of explosion energies particularly for low progenitor initial masses. It will be interesting to see if this trend is robust in larger samples of lightcurves.
5.2.2. Nickel mixing
In Figure 7 we investigate the dependence of nickel mixing length on the progenitor initial mass for the best fitting CURVEPOPS model. As explained above, we have calculated models with a nickel mixing mass coordinate defined as M excised + X M ejecta where the parameter X = [0.1, 0.5, 0.9] with larger values indicating more mixing. We interpolate between these to find the best fitting value.
Figure 7. The nickel mixing length parameter derived from lightcurve fitting and its dependence on initial mass, with symbols as described in Figure 4.
Unfortunately, the extremely large uncertainties on our inferred Nickel mixing lengths makes definite conclusions from this analysis impossible. We note that our very coarse grid (with only three samples) is likely inflating these uncertainties, although there may well be underlying degeneracies in the fit parameters, as demonstrated by the projections of χ 2 parameter space for individual fits illustrated in the Appendices. We do perhaps see a hint that lower mixing lengths might be favoured for high mass progenitors, but further work with a much finer grid of models will be required to verify this and explore its origin further.
5.2.3. Explosion epoch
During CURVEPOPS fitting, an explosion epoch must be derived, particularly for those events in which the pre- and early post-explosion light curve is poorly sampled. Initial constraints on the SNe explosion epochs are determined by searching the literature for the last reported non-detection of each object to act as the upper limit on its age and the supernova discovery date to be the lower limit. These are refined in the literature as described in section 4.2. The details of constraints, along with fitted explosion epochs, are summarized in Table 4.
Table 4. Constraints and Fitted Results of Explosion Epochs in MJD
6. Discussion
A first important question is to evaluate how well our fitting method has recovered the initial and nickel masses of the SNe we have studied. We can see in Figures 4 and 5 that in general the both of these masses are consistent with those from pre-explosion mass constraints from the literature. For the initial mass we find that at first glance there is little difference between the fits that are constrained by the progenitor mass from pre-explosion images and those that are not. A closer look reveals that the scatter is less in the latter case. This is probably due to the degeneracy in SN light curves between mass and energy. As for the nickel mass we again reproduce the literature values, however for some of the greater values the scatter is larger. We have also compared our masses to those from estimating the age of co-eval stars surrounding the progenitor site as estimated by Maund (Reference Maund2017). Again here we see that we are at least consistent with these masses in most cases. Other indirect measurements of a few of the progenitor mass are also possible, for example Xiao et al. (Reference Xiao, Galbany, Eldridge and Stanway2018), and these in general also agree with the masses we derive.
While our derived masses are comparable to those of Smartt (Reference Smartt2015) in general we derive higher masses. This is agreement with others such as Morozova et al. (Reference Morozova, Piro and Valenti2018) with a similar lightcurve fitting method and Davies & Beasor (Reference Davies and Beasor2018) who re-evaluated the pre-supernova luminosity of the detected progenitor stars.
We note that Dessart & Hillier (Reference Dessart and Hillier2019) and Goldberg et al. (Reference Goldberg, Bildsten and Paxton2019) have pointed out the difficulty of deriving initial masses of supernova progenitors from the lightcurves. We do find our results do have large uncertainties in some cases. However our results are dependent on for example our matching of the circumstellar environment to the stellar progenitor model as well as other caveats. At best our estimates give a relative estimate of whether the progenitors are more or less massive, similar to the broad scheme suggested by Chevalier et al. (Reference Chevalier, Fransson and Nymark2006).
For several of the supernovae in our sample, 2004et, 2005cs, 2008bk, 2009md, 2012aw and 2013ej the fits have some part that does not match the observed light curves. We suspect this is due to the limited scope of this study which only uses single-star progenitor models. We expect that if we were to also explode progenitor models that had undergone binary interactions we should be able to have a greater variety of light curves (see paper I for example). The other parameter we have not varied is initial metallicity of the progenitors. Here we have only used one metallicity, but varying this will change in subtle ways the progenitor structure and also the density of the circumstellar medium, by changing the stellar wind mass-loss rates and the wind velocity. The alternative would be to adopt a similar approach to Morozova et al. (Reference Morozova, Piro and Valenti2018) and instead compute models over a range of circumstellar medium parameters rather than those from the progenitor model. Their work does suggest in the cases where we achieve a poor fit we may be underestimating the amount of circumstellar material. Both these approaches could allow us to find a better fit, however the computational demand of calculating synthetic light curves from more progenitor models or more circumstellar environments, while still varying the explosion parameters as here, for each progenitor is extreme. We have demonstrated here that the CURVEPOPS concept is useful and can thus now begin to undertake this mammoth set of numerical calculations.
Morozova et al. (Reference Morozova, Piro and Valenti2018) found good fits to a number of the same supernova lightcurves by varying the amount of circumstellar material around each of their progenitor models. In each case the amount required was significant, of the order of a few times 0.1 Mʘ. Here we have also included the circumstellar medium and achieve a similar early brightening of our theoretical lightcurves. In contrast to Morozova et al. (Reference Morozova, Piro and Valenti2018) we achieved this with a mass-loss rate and wind velocity determined from the progenitor model. In addition, our use of the the wind acceleration model of Moriya et al. (Reference Moriya, Förster, Yoon, Gräfener and Blinnikov2018) means that the density of the circumstellar medium is higher close to the star before decreasing out to the density expected for a freely expanding wind. This produced similar high densities as found by Morozova et al. (Reference Morozova, Piro and Valenti2018). This suggests how the circumstellar medium is modelled does not matter so much as the mass of material that is close to the progenitor star upon explosion. Future investigations of progenitors at different metallicities and thus different wind parameters will enable us to understand the importance of the circumstellar medium to a greater degree.
Smartt et al. (Reference Smartt, Eldridge, Crockett and Maund2009) proposed that the evolution of explosively-synthesised nickel mass with progenitor initial mass was mediated by the amount of oxygen in the star at the supernova epoch, and specifically considered the ratio between oxygen mass and carbon-oxygen (CO) core mass at explosion. In Figure 8 we consider the same parameters based on light curve fitting alone. In common with Smartt et al., we see a trend to lower nickel masses for supernovae with lower initial progenitor masses. The only outliers from this trend are supernova for which the light curve fit is particularly poor and the inferred parameters unreliable. Overplotted on the figure we show a plausible relation between these quantities. The dashed line indicates the size of the carbon-oxygen mantle that surrounds the forming compact remnant, scaled to match the data at 15 Mʘ. This appears to track the datapoints (subject to the substantial uncertainties) suggesting that this may be an important parameter in determining the nickel mass.
Figure 8. The nickel mass derived from lightcurve fitting and its dependence on fitted initial mass, with symbols as described in Figure 4. Lines show a simple linear relation between the logarithms of initial and nickel mass (dotted) and also a model in which the nickel mass is proportional to the amount of mass in the CO core before collapse which is not removed into a compact remnant (dashed line, see text).
7. Conclusion
This article is the second in the CURVEPOPS series and again shows the utility of calculating large grids of synthetic supernova light curves for comparison to observed supernovae. While not every fit to an observed supernova was good, in general the population of results can be used to find various relationships between initial progenitor mass and explosion energy, nickel masses or nickel mixing.
We stress that there are limitations to our study, which are primarily caused by computational limits on the number of models considered. The only solution is to calculate more models to allow for the full diversity of stellar structure and circumstellar environment around each star. We estimate that to calculate all the necessary synthetic light curves requires would be approximately 25 million supernova simulations which would take approximately 120 million CPU hours. This is certainly possible but even with significant computing resources it still takes time and is beyond the scope of this article.
From considering the fits together we find that,
1. The typical explosion energy to be input into SNEC is log(E exp/ergs) = 50.52 ± 0.10.
2. We find a relation between nickel mass and initial mass which may track the size of the carbon-oxygen core at core collapse.
3. We find suggestions of a weak dependence of nickel mixing on initial mass with less mixing when there is a more massive ejecta and initial progenitor.
4. It is possible to achieve strong constraints on the progenitors of type IIP supernovae from the light curves alone.
5. As found by Morozova et al. (Reference Morozova, Piro and Valenti2018) and Moriya et al. (Reference Moriya, Förster, Yoon, Gräfener and Blinnikov2018) it is important to include the circumstellar material surrounding the progenior stars to correctly model type IIP supernova lightcurves. However exactly how to include this in the lightcurve modelling requires further study.
6. Good fits are not possible for every observed supernovae which suggests that there are other factors at play in specifying the shape of type IIP light curves. These include the initial stellar metallicity and binary interactions.
In summary, we can produce progenitor constraints independent of progenitor observations and rivaling them in quality and, where progenitor imaging exists, we add additional data and tighten the uncertainties on key parameters. Finally the synthetic light curves and SNEC input files are freely available from the BPASS website and PASA data store as a resource for the community to use. This data will be continually added to as more simulations are computed.
Acknowledgements
This work would not have been possible without use of the NeSI Pan Cluster, part of the NeSI high-performance computing facilities. New Zealand’s national facilities are provided by the NZ eScience Infrastructure and funded jointly by NeSI’s collaborator institutions and through the Ministry of Business, Innovation & Employment’s Research Infrastructure programme. URL: https://www.nesi.org.nz. JJE acknowledges support from the University of Auckland and thanks the OCIW Distinguished Visitor Program at the Carnegie Observatories, which funded JJE’s visit to gain advice from Anthony L. Piro on using SNEC and including the circumstellar media around the progenitors. ERS acknowledges support from the University of Warwick. NYG and NR both acknowledges support from University of Auckland’s Summer Scholarships. LX acknowledges support from the China Postdoctoral Science Foundation (grant No.2018M642524). LX acknowledge the grant from the National Key R&D Program of China (2016YFA0400702), the National Natural Science Foundation of China (No. 11673020 and No. 11421303), and the National Thousand Young Talents Program of China.
Appendix A. Results from a free fit of the light curves
Here in Figures A1 to A8 we present the best fitting V-band magnitude lightcurves for the supernovae when a free fit across the full range of modelled initial masses is allowed as well as corner plots showing how the χ2 varies with the five parameters we fit. We also include plots showing how χ2 only depends on one parameter at a time as shown by the solid black line while the horizontal dashed lines show the Δχ2 for the 1σ, 2σ, and 3σ. In the bottom figure of the lightcurves we plot the best fitting model (black line) along with the lightcurves at are within the 1σ uncertainty in grey while the observations are shown in red.
Figure A1. SN2003gd Free-Fit, corner plots showing how χ2 varies over the five parameters as well a plot comparing the observed lightcurves to thematching theoretical models.
Figure A2. SN2004A Free-Fit, as in Figure A1
Figure A3. SN2004et Free-Fit, as in Figure A1
Figure A4. SN2005cs Free-Fit, as in Figure A1
Figure A5. SN2006my Free-Fit, as in Figure A1
Figure A6. SN2008bk Free-Fit, as in Figure A1
Figure A7. SN2009md Free-Fit, as in Figure A1
Figure A8. SN2012A Free-Fit, as in Figure A1
Figure A9. SN2012aw Free-Fit, as in Figure A1
Figure A10. SN2012ec Free-Fit, as in Figure A1
Figure A11. SN2013ej Free-Fit, as in Figure A1
Appendix B. Results from a progenitor constrained fit of the light curves
Here in Figures B1 to B8 we present the best fitting V-band magnitude lightcurves for the supernovae when a constrained fit across the 68% uncertainty range of modelled initialmasses allowed by the masses derived by Smartt et al. (2009). We also include corner plots showing how the χ 2 varies with the five parameters we fit. We also include plots showing how χ 2 only depends on one parameter at a time as shown by the solid black line while the horizontal dashed lines show the Δχ 2 for the 1σ, 2σ and 3σ. In the bottom figure of the lightcurves we plot the best fitting model (black line) along with the lightcurves at are within the 1σ uncertainty in grey while the observations are shown in red.
Figure B1. SN2003gd Constrained, corner plots showing how χ 2 varies over the 5 parameters as well a plot comparing the observed lightcurves to the matching theoretical.
Figure B2. SN2004A Constrained, as in Figure B1
Figure B3. SN2004et Constrained, as in Figure B1
Figure B4. SN2005cs Constrained, as in Figure B1
Figure B5. SN2006my Constrained, as in Figure B1
Figure B6. SN2008bk Constrained, as in Figure B1.
Figure B7. SN2009md Constrained, as in Figure B1
Figure B8. SN2012A Constrained, as in Figure B1
Figure B9. SN2012aw Constrained, as in Figure B1
Figure B10. SN2012ec Constrained, as in Figure B1
Figure B11. SN2013ej Constrained, as in Figure B1
Appendix C. Test of changes at end of stellarmodels
Our input grid of stellar evolution and structure models, based on the BPASS modification of the STARS code, was selected for its utility in modelling binary interactions and its compatibility with a large population synthesis modelling project, and thus will allow future generalisations of this study. However an important drawback in the BPASSmodels is that they only follow stellar evolution to the end of core carbon burning. It is entirely possible that important structural changes may occur to massive stars in the very rapid and complex late stages of core nuclear burning, immediately before supernova. This is impossible to explore within BPASS, but is accessible to other stellar evolution codes.
To evaluate the principal effects of late core-burning stages, we have considered the evolution of a 15.6 Mʘ single-star model (i.e. a typical SN IIP progenitor) with the Modules for Experiments in Stellar Astrophysics (MESAr10398, Paxton et al. Reference Paxton, Bildsten, Dotter, Herwig, Lesaffre and Timmes2011, Reference Paxton2013, Reference Paxton2015, Reference Paxton2018) stellar evolution code, which follows the stars through to core collapse. In Figure C1 we illustrate the difference in stellar structure observed as a function of initial mass for the models captured at three stages: at core carbon ignition, the end of core carbon burning and at core collapse.
Figure C1. An initially 15.6 Mʘ stellar model evolved with MESA v10398. The models are taken at the beginning of carbon burning, the end of core carbon burning and the final model output before core-collapse. The BPASSmodels we use are taken after the point of carbon burning completing in the stellar core.
As the figure illustrates, structural differences are observed in the very inner regions of the stellar core (i.e. within R<0.01 R *). In our formalism, these regions are extremely likely to be subsumed within the stellar remnant and so would contribute little if anything to the evolution of the supernova lightcurve. Nonetheless, we note this as a limitation of our modelling which may be addressed in future work.
Appendix D. Effect of circumstellarmedium on supernova lightcurves
Here we show the evolution of absolute V band magnitude for our lightcurve models, as in Figure 1, but to aid clarity each panel only shows one initial mass progenitor star, with either five lines indicating the effect of changing the β acceleration parameter of the stellar wind (Figure D1), or four lines showing the effect of different assumed circumstellar medium densities (Figures D2 and D3).
Figure D1. Model lightcurves in the V-band for stars with initial masses of 21 Mʘ and above. Each panel now only shows one stellar initial mass, while the different lightcurves have varying values of β the wind acceleration parameter.
Figure D2. As in Figure D1, but here we vary the density of the circumstellar medium around the star.
Figure D3. As in Figure D2 but now showing stars in the initial mass range 6<M/Mʘ <20, for which changing the circumstellar medium has little affect.
