Supernova lightCURVE POPulation Synthesis II: Validation against supernovae with an observed progenitor

Abstract

We use the results of a supernova light-curve population synthesis to predict the range of possible supernova light curves arising from a population of single-star progenitors that lead to type IIP supernovae. We calculate multiple models varying the initial mass, explosion energy, nickel mass and nickel mixing and then compare these to type IIP supernovae with detailed light curve data and pre-explosion imaging progenitor constraints. Where a good fit is obtained to observations, we are able to achieve initial progenitor and nickel mass estimates from the supernova lightcurve that are comparable in precision to those obtained from progenitor imaging. For 2 of the 11 IIP supernovae considered our fits are poor, indicating that more progenitor models should be included in our synthesis or that our assumptions, regarding factors such as stellar mass loss rates or the rapid final stages of stellar evolution, may need to be revisited in certain cases. Using the results of our analysis we are able to show that most of the type IIP supernovae have an explosion energy of the order of log(E_exp/ergs) = 50.52 ± 0.10 and that both the amount of nickel in the supernovae and the amount of mixing may have a dependence on initial progenitor mass.

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 1998) and a tentative detection associated with a short gamma ray burst (Tanvir et al. 2013 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. 2017) and exploded with the Supernova Explosion Code (SNEC, Morozova et al. 2015). 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. 2012, 2014; Morozova et al. 2017, 2018). Here we have focused our study on the sample of type IIP supernova progenitors that have detected progenitor stars as described in Smartt (2015). We have collated lightcurve data from the Open Supernova Catalogue ^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 1994, 2005, 2007; Utrobin & Chugai 2008, 2009; Dessart et al. 2010, 2011; Bersten et al. 2011, 2012; Dessart et al. 2014; Bersten et al. 2014; Morozova et al. 2015; Dessart et al. 2016; Morozova et al. 2016, 2017, 2018; Martínez & Bersten 2018, 2019). 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. 2017). The models are calculated using a version of the Cambridge STARS code that has been adapted to follow binary evolution (see Eldridge et al. 2017, 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. 2011, 2013, 2015, 2018) 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. 2013; Dessart & Hillier 2019; Goldberg et al. 2019) 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 2012). 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. 2015, 2016, 2017, 2018). 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. 2016, 2017, 2018; Moriya et al. 2018). 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. (2006). The mass-loss rates used are those of de Jager et al. (1988). We do not vary the wind parameters as in Moriya et al. (2018) 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. (2018) with β = 5 which is typical of cool red supergiants (Schroeder 1985; Gräfener & Vink 2016) although we note this value might not be correct for all red supergiants (Ohnaka et al. 2017). 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 (10^50.5 erg s⁻¹, 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. (1988) are correct or not. Two studies, Mauron & Josselin (2011) and Beasor & Davies (2018), found that while de Jager et al. (1988) have the correct scale of mass-loss rates there are significant departures away from the predicted rates for some stars. Mauron & Josselin (2011) found that red supergiants varied by up to a factor of four away from the de Jager et al. (1988) predictions. While Beasor & Davies (2018) found more significant deviations, especially that de Jager et al. (1988) 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. 2006). 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. (2017). Our grid of values is as follows,

1. 1. The explosion energy varies from log(E_SN/erg s⁻¹ = 50 to 52 in steps of 0.25 dex.

2. 2. The nickel mass varies from log (M_Ni/M_ʘ) = −3 to −1 in steps of 0.25 dex.

3. 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

   Where M _excised is the mass which collapses into the compact remnant, and M _ejecta is the mass of material ejected by the supernova explosion. This therefore assumes significant mixing of nickel into the envelope in all cases, but more mixing in some cases than others.

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 10^50.5 erg s⁻¹, 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 (10^50.5 erg s⁻¹, 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. (2018). 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. (1988) 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 (2015). 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. (2005b)

^c Maund et al. (2014)

^d Fraser et al. (2011)

^e Maguire et al. (2010)

^f Tomasella et al. (2013)

^g Fraser et al. (2014)

^h Crockett et al. (2011)

ⁱ Van Dyk et al. (2012)

^j Maund et al. (2013)

^k Through Open Supernova Catalog: Guillochon et al. (2017)

Distances were taken from cross-referencing literature, NED^b and SIMBAD (Wenger et al. 2000) 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. 1. Mask all model points where the neighbouring value differs from the previous by more than 0.1 magnitude.

2. 2. Take the resulting array and mask again if the i + 2 value differs from i by more than 0.1 magnitude.

3. 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 χ ² 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

(1) \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 χ² and estimate the uncertainty in each of the parameters as the range of values of each fit parameter for which Δχ² < 5.89. Best fit and optimal χ² 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 χ² 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 (2015)

^b Hendry et al. (2005a)

^c Hendry et al. (2006)

^d Smartt et al. (2009)

^e Anderson et al. (2014)

^f Fraser et al. (2011)

^g Tomasella et al. (2013)

^h Dall’Ora et al. (2014)

ⁱ Jerkstrand et al. (2015b)

^j Bose et al. (2015)

Table 3. Reference and Progenitor Constrained Fitting Parameters

^a From progenitor observations and modelling, as in Smartt (2015)

^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. (2005a)

^d Hendry et al. (2006)

^e Smartt et al. (2009)

^f Anderson et al. (2014)

^g Fraser et al. (2011)

^h Tomasella et al. (2013)

ⁱ Dall’Ora et al. (2014)

^j Jerkstrand et al. (2015b)

^k Bose et al. (2015)

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 (2015). 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. (2013). 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 (2015) 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 (2015) 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 (2015) 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 (2015) 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 (2017), the similar study using SNEC of Morozova et al. (2018) and updated progenitor masses from Davies & Beasor (2018). 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. (2018) 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. (2018). We are uncertain why we achieve such a different fit. As the χ ² 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. (2015a). 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 (2015). 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 χ ² 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 (2017). 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. (2018), and these in general also agree with the masses we derive.

While our derived masses are comparable to those of Smartt (2015) in general we derive higher masses. This is agreement with others such as Morozova et al. (2018) with a similar lightcurve fitting method and Davies & Beasor (2018) who re-evaluated the pre-supernova luminosity of the detected progenitor stars.

We note that Dessart & Hillier (2019) and Goldberg et al. (2019) 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. (2006).

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. (2018) 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. (2018) 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. (2018) 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. (2018) 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. (2018). 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. (2009) 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. 1. The typical explosion energy to be input into SNEC is log(E _exp/ergs) = 50.52 ± 0.10.

2. 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. 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. 4. It is possible to achieve strong constraints on the progenitors of type IIP supernovae from the light curves alone.

5. 5. As found by Morozova et al. (2018) and Moriya et al. (2018) 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. 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.