Connecting the Light Curves of Type IIP Supernovae to the Properties of their Progenitors

Observations of core-collapse supernovae (CCSNe) reveal a wealth of information about the dynamics of the supernova ejecta and its composition but very little direct information about the progenitor. Constraining properties of the progenitor and the explosion requires coupling the observations with a theoretical model of the explosion. Here, we begin with the CCSN simulations of Couch et al 2020 ApJ 890 127, which use a non-parametric treatment of the neutrino transport while also accounting for turbulence and convection. In this work we use the SuperNova Explosion Code to evolve the CCSN hydrodynamics to later times and compute bolometric light curves. Focusing on SNe IIP, we then (1) directly compare the theoretical STIR explosions to observations and (2) assess how properties of the progenitor's core can be estimated from optical photometry in the plateau phase alone. First, the distribution of plateau luminosities (L$_{50}$) and ejecta velocities achieved by our simulations is similar to the observed distributions. Second, we fit our models to the light curves and velocity evolution of some well-observed SNe. Third, we recover well-known correlations, as well as the difficulty of connecting any one SN property to zero-age main sequence mass. Finally, we show that there is a usable, linear correlation between iron core mass and L$_{50}$ such that optical photometry alone of SNe IIP can give us insights into the cores of massive stars. Illustrating this by application to a few SNe, we find iron core masses of 1.3-1.5 solar masses with typical errors of ~0.05 solar masses. Data are publicly available online (\url{https://doi.org/10.5281/zenodo.6631964}).


INTRODUCTION
Core-collapse supernovae (CCSNe) are the explosive deaths that result from the ends of stellar evolution for massive stars with zero-age main sequence (ZAMS) masses M ZAMS 8M . The current understanding suggests that some fraction of possible progenitors will successfully produce CCSNe while others will fail and produce a black hole  Ertl et al. 2016;Sukhbold et al. 2016;Adams et al. 2017;Sukhbold et al. 2018;Couch et al. 2020). The details of the explosion mechanism have been the subject of decades of work with current work favoring, for most progenitors, the delayed neutrino-driven mechanism (Bethe & Wilson 1985). For an in-depth review of the CCSN explosion mechanism and related problems, see recent reviews (e.g., Bethe 1990;Janka et al. 2007Janka et al. , 2012Janka et al. , 2016Burrows 2013;Hix et al. 2014;Müller et al. 2016;Couch 2017;Pejcha 2020).
CCSNe are detectable by three primary messengers -EM waves, neutrinos, and gravitational waves (GWs). Neutrino and GW signals have the very desirable property that they are arXiv:2102.01118v2 [astro-ph.HE] 10 Jun 2022 emitted directly from the core of the star at the time of collapse and may reveal information about the structure there (e.g., Pajkos et al. 2019Pajkos et al. , 2020Warren et al. 2020;Sotani & Takiwaki 2020), unlike photons which are emitted from the photosphere in the outer layers of the supernova ejecta until the remnant phase. However, to date there has been only one detection of neutrinos from a supernova (Arnett et al. 1989. With modern neutrino detectors, only CC-SNe occurring within our galaxy may be detectable (Scholberg 2012). Similarly, there have been no confirmed detections of GW emission from a CCSNe. The current suite of detectors (aLIGO, Virgo, and KAGRA) can only detect GWs from a CCSNe if it occurs within a distance of ≤ 100 kpc (Abbott et al. 2016). It is the case, however unfortunate, that the overwhelming majority CCSNe will only be observed in EM signals.
The focus of this paper is connecting EM signals to progenitor properties for SNe IIP. These events have been shown to originate from red supergiant progenitors (Van Dyk et al. 2003;Smartt 2009;Van Dyk et al. 2019). Despite being the most common type of CCSNe, their diversity of observable features -such as light curve morphologies -is still not fully understood (e.g., Anderson et al. 2014;Valenti et al. 2016). The connection between SNe IIP and IIL supernovae, for example, still remains an open question -whether IIL's are the limit of IIP's as the H envelope is depleted or a separate class (Barbon et al. 1979;Blinnikov & Bartunov 1993;Faran et al. 2014; Morozova et al. 2015).
Understanding the connection between SNe IIP light curves and stellar progenitors has a new urgency. Coming next-generation telescopes such as the Vera C. Rubin Observatory and its primary optical photometry survey, The Rubin Observatory Legacy Survey of Space and Time (LSST) (Ivezić et al. 2019), will allow for extremely deep imaging of the entire sky every couple of nights. The LSST will allow for statistical studies of populations of CCSNe of an unprecedented scale (for recent statistical studies see, e.g., Anderson et al. 2014;Sanders et al. 2015;Gutiérrez et al. 2017a,b).
Ultimately, properly characterizing the diversity in SNe II supernova light curve morphology will require the union of observation and theory. On the theory side, this comprises realistic stellar evolution models including the core collapse, following the resulting explosion with robust physics, and calculating EM light curves (as well as neutrino and GW signals). The gold standard is full three-dimensional (3D), selfconsistent simulations. Core-collapse supernovae and their progenitors are truly 3D in nature and the key to understanding the diversity of light curve morphology lies in faithfully modeling these asphericities (Wongwathanarat et al. 2013Dessart & Audit 2019;Stockinger et al. 2020;Sandoval et al. 2021). 3D simulations are, however, computationally expensive to perform and, as such, are limited in number and the range of parameter space that they cover. Spherically symmetric (1D) simulations remain necessary for understanding the CCSNe explosion mechanism and their observables by surveying landscapes of possible CCSNe. Great progress has been made in the last few years regarding 1D CCSN simulations (Ebinger et al. 2017;Sukhbold et al. 2016;Couch et al. 2020), allowing for successful explosions in 1D using neutrino-driven explosions across wide ranges of progenitor masses. These 1D simulations allow for large parameter studies performing potentially thousands of simulations spanning ranges of progenitor masses, equations of state, and metallicities, for example.
Light curve calculations are the final, crucial piece of the theoretical process of understanding these explosions. Commonly, calculations of synthetic bolometric light curves of CCSNe invoke a thermal bomb or piston-driven model, where energy is artificially injected into a thin region above a user-specified mass cut within the progenitor (see, e.g., Bersten et al. 2011;Morozova et al. 2015;Ricks & Dwarkadas 2019). In these models, the explosion energy is a user-set parameter instead of being determined by the structure of the progenitor and explosion physics. The calculations cannot determine whether a given progenitor will result in a successful supernova or fail to revive its stalled shock and collapse to a black hole. The explodability has been shown to have non-trivial behavior across a large range of ZAMS mass progenitors Ebinger et al. 2017;Sukhbold et al. 2018;Couch et al. 2020) and cannot be captured with more simplified models. The clear next step is the coupling of high fidelity CCSN simulations with bolometric light curve calculations.
Light curves contain information about their progenitor and the explosion -properties such as the composition of the ejecta, mass or radius of the progenitor, or explosion energy may be inferred (Litvinova & Nadezhin 1985;Popov 1993;Kasen & Woosley 2009;Sukhbold et al. 2016). The process of inferring progenitor and explosion properties from light curves has been shown to be highly degenerate (Goldberg et al. 2019;Dessart & Hillier 2019) with many combinations of properties being capable of producing a given light curve. Of particular interest, however, is the early time light curve dominated by radiation streaming form the shock heated outer envelope. This early time behavior can be compared to shock cooling models to put constraints on the stellar radius (Nakar & Sari 2010;Tolstov et al. 2013;Shussman et al. 2016;Kozyreva et al. 2020). Recently Morozova et al. (2016); Rubin & Gal-Yam (2017) explored the effectiveness and temporal limitations of these models and these methods have been widely used for constraining the progenitor preexplosion radius (e.g., Rabinak & Waxman 2011;Gall et al. 2015;González-Gaitán et al. 2015;Sapir & Waxman 2017;Soumagnac et al. 2020;Vallely et al. 2021). These early time observations may help to break the degeneracies between progenitor and explosion properties (Goldberg & Bildsten 2020).
In this work, we calculate the bolometric light curves of the recent 1D simulations done with the FLASH 1 (Fryxell et al. 2000;Dubey et al. 2009) code using the new Supernova Turbulence In Reduced-dimensionality (STIR) model ). This 1D convection scheme has the benefit of being more consistent with some properties of full physics 3D CCSN simulations -such as explosion energies and landscapes -while leaving the neutrino physics unaltered. Like any 1D method, it remains a simplification of the full picture and is not without its shortcomings (e.g., Müller 2019). Similar 1D schemes have also been used to study Rayleigh-Taylor instabilities in supernova remnants (Duffell 2016). The initial conditions of these models are set by the 1D stellar evolution models of Sukhbold et al. (2016), which make up a suite of 200 solar metallicity, non-rotating massive stars between 9 and 120 M . We couple the final state of the STIR simulations with the SuperNova Explosion Code (SNEC) (Morozova et al. 2015), which follows the explosion through the rest of the star and through the plateau and nebular phases of the light curves. We will demonstrate that using a more sophisticated 1D explosion model to determine a distribution of explosion energies consistent with 3D simulations imparts non-trivial features to observables and thus properties inferred from them, highlighting the importance of the explosion model used. With this set of light curves, we make available a new set of theoretical predictions to compare directly with observations. Furthermore, we investigate direct correlations between progenitor properties and light curve properties. We recover known correlations, and we quantify the dependence of SNe IIP luminosity on the progenitor iron core mass at time of collapse -thus providing a way of obtaining core properties from EM signals without the need for the much rarer neutrino and GW signals.
This paper is laid out as follows: in Section 2 we discuss the various progenitors, codes, and statistical methods that are used in this study. Section 3 presents our results: 3.1 presents observable properties of our light curves and their trends across ZAMS mass, 3.2 presents preliminary comparisons to observations of SNe IIP, 3.3 shows correlations found between light curve and progenitor properties. In Section 4, we summarize our results and briefly discuss comparison to other theoretical light-curve calculations and prospects for future work.

METHODS
For this work, we begin with massive stellar progenitors evolved up to the point of core collapse in Sukhbold et al. 1 https://flash.rochester.edu/site/ (2016) evolved using KEPLER. The core collapse and following explosion or collapse to BH are simulated using the FLASH simulation framework (Fryxell et al. 2000;Dubey et al. 2009) with the Supernova Turbulence In Reduceddimensionality (STIR) model , the details of which are discussed in Section 2.2. The output of the STIR simulations are mapped into the SuperNova Explosion Code (SNEC) (Morozova et al. 2015(Morozova et al. , 2016(Morozova et al. , 2018 to generate bolometric light curves as discussed in Section 2.3. In Section 2.4 we present the methods used to analyze statistical relationships between properties of the progenitor and observables.

Progenitors
We begin with the 200 non-rotating, solar metallicity models of Sukhbold et al. (2016). These models cover a range of ZAMS masses from 9 -120 M and were created with the KEPLER code assuming no magnetic fields or rotation and single star evolution. Progenitors with ZAMS masses above 31M experienced significant mass loss during their lifetimes and did not explode as SNe IIP and this is the upper limit on our mass range (see Sukhbold et al. 2016, for details on their stellar evolution). The more massive Type I SNe progenitors are too few in number to perform a meaningful statistical analysis and we defer their analysis for future work. This leaves 136 progenitors producing SNe IIP supernovae used in this work.
These progenitors span a wide range of possible CCSNe progenitor properties. Figure 1 shows the mass of the Hrich envelope as a function of pre-supernova radius (top) and the stellar pre-supernova mass as a function of ZAMS mass (bottom). Here we show only models that successfully exploded in Couch et al. (2020) and are included in this work. Gaps in this figure -such as that from about 12-15M -represent models which failed to explode and are not included in this work. These progenitors become mass-loss dominated around 23M , as seen in the bottom panel of Figure 1. This complicates correlations between quantities of interest and tends to cause them to deviate from monotonicity. This is key to investigating observable trends in light curves across a wide range of progenitors, as we demonstrate later.
The progenitors in Sukhbold et al. (2016) were further investigated in Sukhbold et al. (2018) using a set of high resolution stellar evolutionary models. They showed that the features of these progenitors -notably the compactness landscape -was not numerical in nature and was present in their high resolution models. Similarly, other, recent works have found similar trends in the presupernova mass and compactness (e.g., Laplace et al. 2021, using MESA (Paxton et al. 2011, 2013, 2015, 2019). We note that while general trends may be reproduced, details such as the apparent 'chaos' seen in Sukhbold et al. (2016Sukhbold et al. ( , 2018 are sensitive to implementation details of stellar evolution and may not appear in other studies (Chieffi & Limongi 2020). Using a different set of progenitors with a different compactness curve would likely affect the explosion landscape but isn't expected to affect the results presented here.
A key result of systematic 1D studies of CCSNe are the socalled "islands of explodabilty" . The final result of stellar core collapse -a successful or failed explosion -is not a monotonic function of ZAMS mass. Instead, the explodability of the progenitor is sensitive to the core structure at the time of collapse. While the placement of these islands is sensitive to the explosion model and the progenitors used, it is a feature that has now been seen amongst many groups (O'Connor & Ott 2011;Perego et al. 2015;Sukhbold et al. 2016;Ertl et al. 2016;Ebinger et al. 2019;Couch et al. 2020). However, studies using exclusively thermal bomb driven explosions uninformed by neutrino driven calculations cannot reproduce the explosion/implosion fate of a progenitor and are insensitive to this feature. Any systematic study of light curves from populations of SNe must capture this complex behavior.

FLASH
The CCSN simulations were conducted in Couch et al. (2020) using the FLASH code framework with the STIR turbulence-aided explosion model. This model is a new method for artificially driving explosions in 1D CCSN simulations.. Turbulence is key in simulating successful, realistic explosions, as turbulence may constitute 50% or more of the total pressure behind the shock (Murphy et al. 2013;Couch & Ott 2015) and turbulent dissipation is important for postshock heating (Mabanta & Murphy 2018). The combined impact of these effects is to aid the explosion. The inclusion of turbulent effects allows for successful explosions in 1D simulations while reproducing the results seen in 3D simulations from various groups ) without the need for parametrized neutrino physics.
STIR models turbulence using the Reynolds-averaged Euler equations with mixing length theory as a closure. This model has one primary scalable parameter, the mixing length parameter α Λ , inherited from mixing length theory which scales the strength of convection. The mixing length parameter has been tuned to fit STIR simulations to full 3D simulations run with FLASH and reproduces 3D results seen in FLASH and other codes, noting particularly good agreement with the 3D works of (Burrows et al. 2020). We use the fiducial value found in Couch et al. (2020) for the mixing length parameter, α Λ = 1.25. STIR also includes four additional diffusion parameters that control the convective mixing of internal energy, turbulent kinetic energy, composition, and neutrinos. As in Couch et al. (2020), all four of these diffusion coefficients are set to 1/6, a value consistent with comparison to fully 3D simulation of convection in massive stars (Müller et al. 2016). We note that the convective dynamics are insensitive to the choice of diffusion coefficients and, thus, impacts on the explosion are negligible (Müller et al. 2016;Boccioli et al. 2021). FLASH with the STIR model has the desireable benefit that there is no need to tune the model to match a specific observation. Instead, its one primary parameter is tuned to be consistent with multi-physics 3D CCSN simulations, reducing the possibility of inserting biases into the results. STIR includes neutrino transport using a state-of-the-art two moment method with an analytic "M1" closure (Shibata et al. 2011;Cardall et al. 2013;O'Connor 2015;O'Connor & Couch 2018). We simulate three neutrino flavors: ν e ,ν e and ν x , where ν x combines the µτ neutrino and antineutrino flavors. M1 transport requires no tuning and has no free parameters (up to the choice of a closure for the highorder radiation moments), allowing for truly physics-driven explosions. The STIR simulations use the now commonly adopted, empirically-motivated "SFHo" equation of state for dense nuclear matter (Steiner et al. 2013) which is able to replicate observed neutron star masses.
At the end of the STIR simulations, the explosion energies for all but the highest-mass progenitors have asymptoted. It is commonplace in CCSNe work to define the explosion energy as the sum total energy, from all sources, of material that has both positive total energy and positive velocity (e.g., Bruenn et al. 2016) at the end of the simulation. This is zero during the stalled shock phase, when all of the material is still gravitationally bound, and becomes positive if/when the shock begins to move outward again due to neutrino heating and other effects. This energy, once it has reached its asymptotic value, represents the energy that is injected into the rest of the star to drive the explosion and unbind the stellar material. When discussing the combined STIR + SNEC simulations, this is the explosion energy that we will reference. It is important to note that this energy is different than the energy that would be used in hydrodynamical modelling (e.g., thermal bomb explosions). In the thermal bomb regime, a user set energy is deposited at t = 0 over defined temporal and spatial extent, and assumes that the energy of the shock comes directly from the core-bounce which is inconsistent with the physical picture of CCSNe. In the case of high fidelity simulations, a large amount of material has already been gravitationally unbound by the shock when the explosion energy is measured. A thermal bomb model with the "same" energy injected into the inner zones would, by the time the same amount of material is unbound, be less energetic by exactly the binding energy of the material. Care should be taken when comparing energetics from these two approaches. While the physics of these two explosion methods are inconsistent with each other, the thermal bomb energetics can be made consistent with neutrino-driven energetics by correcting the bomb energy by the binding energy of the material between the shock and the PNS surface (this material is already unbound when the explosion energy is calculated as above, but in thermal bomb or piston-driven explosions it is not). Without this correction, a thermal bomb model using energetics from neutrino-driven simulations will have less energy available for the explosion, impacting observables. Figure 2 shows the explosion energies obtained with STIR (black) alongside the explosion energy with the progenitor's overburden energy removed (blue). The progenitor's overburden energy is the (negative) total energy above the shock that the explosion must overcome to unbind the star (Bruenn et al. 2016). The total energy, which we compute as the total energy on the computational domain after the explosion has set in is closer to what will characterize the ejecta. Gaps in mass, such as from about 12M to 15M , indicate regions where progenitors failed to successfully launch an explosion in STIR. The bottom panel shows the explosion energy as a function of the iron core mass. These explosion energies are set largely by the structure of the cores of their progenitorseffects which can only be seen by employing neutrino driven explosions (for recent examples of the impacts of core structures on explosions and observable signatures, see Warren et al. 2020;Burrows et al. 2020). The emerging picture from high fidelity simulations is that there is no simple relationship between explosion energy and ZAMS mass, instead requiring multi-physics simulations to determine robustly Ebinger et al. 2017Ebinger et al. , 2019Sukhbold et al. 2018;Couch et al. 2020;Burrows et al. 2020;Ertl et al. 2020). The explosion energy is more closely related to the pre-supernova mass and properties of the core, such as the compactness parameter or the mass of the iron core.

SNEC
We simulate light curves for all of the models that successfully produced explosions in Couch et al. (2020) (see their Figure 6, middle row). This is all but about 50 of the original 200 progenitors. This limits our study to light curves obtained from progenitors that actually explode, allowing us to explore solely relationships that come from physically-driven explosions. At the end of the STIR simulations, the final states are mapped into the SuperNova Explosion Code (SNEC) (Morozova et al. 2015). SNEC is a spherically symmetric, Lagrangian, equilibrium flux-limited diffusion radiation-hydrodynamics code and is publicly available 2 . Unlike STIR, SNEC does not include any form of general relativistic gravity, neutrino transport, or dense matter EoS, which are all important for modeling the explosion but not necessarily for computing the light curve. Instead, it follows the basic physics needed for predicting bolometric supernova light curves. SNEC includes Lagrangian Newtonian hydrodynamics with artificial viscosity following the formulation in Mezzacappa & Bruenn (1993) and a stellar equation of state following Paczynski (1983) that includes contributions from radiation, ions, and electrons with approximate electron degeneracy. This is used in tandem with a Saha ionization solver that can follow ionization of any number of present elements. At high temperatures SNEC uses OPAL Type II opacities (Iglesias & Rogers 1996) suitable for solar metallicity. These opacities are supplemented by those of Ferguson et al. (2005) at low temperatures.
1D modeling cannot properly capture the mixing at compositional interfaces due to Rayleigh-Taylor and Richtmyer-Meshkov instabilities, for example. Without mixing, sharp compositional gradients appear that produce features in light curves that are not observed in nature (Utrobin 2007). In these mixing processes, shock propagation outwards can cause light elements to mix inwards and heavy elements to mix outwards . Of particular importance is the mixing of radioactive 56 Ni, whose mixing extent affects the light curve properties (Morozova et al. 2015). SNEC applies boxcar smoothing that smooths out compositional profiles, simulating mixing and avoiding unphysical light curve bumps. We use the fiducial parameters of Morozova et al. (2015) for our boxcar smoothing.
In the present work we follow the ionization of 1 H, 3 He, and 4 He, similarly to Morozova et al. (2015). H and He make up the majority of the energy contributions from recombination relevant for producing bolometric SNe IIP light curves. Our STIR simulations do not currently track detailed compositional information in their output. When mapping into SNEC, we fill the composition in the STIR part of the domain to be pure 4 He. This has no noticeable effect on the 2 http://stellarcollapse.org/SNEC  light curves in this study (see Appendix A). Figure 3 shows mass fractions of 1 H (light blue), 4 He (dark blue), 12 C (gold), 16 O (red), and 56 Ni (black, dot-dashed line). The solid lines show the unmixed profiles that are input to SNEC. Notably, the gray region shows the STIR domain where the composition, prior to mixing, is set to pure 4 He. The dashed lines show the compsotion after boxcar smoothing is applied. The bottom panel shows the radial density profile in the STIR domain (solid line) and in SNEC after mapping (dashed line). The final, critical ingredient for powering a SNe light curve is radioactive heating from the 56 Ni → 56 Co → 56 Fe decay chain. Radioactive 56 Ni is produced in explosive nuclear burning during the first epochs of the explosion in the inner parts of the star. Hydrodynamic instabilities mix the 56 Ni outward. Gamma-rays and positrons emitted from the decay process diffuse outward and provide an additional source of energy. Capturing this is crucial as, after the end of the plateau phase, the light curve is powered entirely by this radioactive decay. SNEC follows the radiative transfer of gamma-rays from the 56 Ni and 56 Co decays using the gray transfer approximation (Swartz et al. 1995) and the result-ing energy release is coupled to the hydrodynamics independently from the rest of the radiation.
Currently, neither our STIR models used here nor the public version of SNEC include a nuclear reaction network. To alleviate this issue, SNEC allows for a user specified amount of 56 Ni to be input by hand throughout a specified mass coordinate. Sukhbold et al. (2016) simulate the explosions of these progenitors including a large nuclear reaction network, and we use the 56 Ni yields as a function of explosion energy from their work (see their Table 4, Figure 17) to estimate a mass of 56 Ni from that relationship to be distributed by SNEC. For all but the lightest progenitors, they find around 0.07M of 56 Ni. We disperse the 56 Ni up to about 75% of the way through the He shell -avoiding mixing into the H envelope. This provides control amongst the progenitors. As the mixing extent must be set by hand, any further treatment would require a large parameter study. In recent high fidelity models, mixing of radioactive 56 Ni into the H-rich envelope is realized (Utrobin et al. 2015(Utrobin et al. , 2017Stockinger et al. 2020;Utrobin et al. 2021;Sandoval et al. 2021) and is expected to occur in at least some of our models. In Morozova et al. (2015), they showed that variations in these distributions had little effect on the light curve, especially on the plateau (see their Figure 6). Goldberg et al. (2019) show slight variations in the light curve as it falls off the plateau depending on the extent of mixing (see their Figure 10). Kozyreva et al. (2019) explore the effects of mixing prescriptions for 56 Ni, such as uniform or boxcar, on light curves, showing differences on the plateau between these methods. The lack of a reaction network consistently incorporated into the calculations forms a weakness of the current work, despite being based on nucleosynthetic calculations and tuned to our explosion energies. However, the main results of this work (see Section 3.3) use quantities measured on the plateau where they are less sensitive to reasonable variations in 56 Ni mass and distribution. Future work will include nucleosynthesis calculations with the STIR input models to properly seed the SNEC calculations.
Typically, high fidelity CCSN simulations do not simulate the entire star -instead focusing on the inner 15,000 km or so necessary for launching the explosion. We must stitch the STIR simulation data, with the explosions developed on the grid, onto the progenitor pre-explosion profile outside the STIR boundary (15,000 km) in order to simulate the full star. Below the shock, mass profiles are taken from STIR. Above the shock mass coordinate, mass profiles are taken from the progenitor profiles. These smooth, combined STIR -preexplosion progenitor profiles are used as the inputs to SNEC as detailed in (Morozova et al. 2015). One advantage of using the STIR models as the initial conditions to SNEC is that the high fidelity equation of state and neutrino transport yield a physically realistic remnant mass to motivate the mass cut -an amount of material not included in light curve simulations that should be close to the remnant mass. We place a mass cut outside the PNS at the point where the total energy becomes positive -removing both the PNS and a small amount of still gravitationally bound material above it (of order 0.0001M ). For all of the simulations we use 1000 cells in the SNEC domain using a geometric grid, as in Morozova et al. (2015), that places higher resolution in the core around the shock and at the outer domain to resolve the photosphere. Our grid is slightly modified from that of Morozova et al. (2015) to place added resolution in the core over the already existing explosion. Simulations were run until 300 days when possible to adequately sample both the plateau and the tail for all events.
To simulate CCSNe directly from progenitors, SNEC has the ability to artificially drive an explosion with a piston or thermal bomb. One of the primary qualities of our method is to eliminate the need for this and thus eliminate user input explosion energies, which can take any range or distribution, replacing them with physically motivated energetics. However, for some of the more massive progenitors in this study, the explosion energies were still increasing by the time the shock reached the outer boundary. Eventually, energy generation from neutrino heating and other sources will slow as the shock expands and the explosions energies will asymptote. Since our computational domain is limited to 15,000 km, some progenitors do not reach their "true" explosion energies. In order to fully capture the energy of the explosion in STIR, we integrate the neutrino heating in the gain region at the end of the FLASH simulations to estimate the asymptotic explosion energy and add the difference -at most about 0.3×10 51 erg -as a thermal bomb over the shocked region. These additions are most necessary in the region of high energy between about 21M and 25M where the final energies were still readily increasing. This energy is what is displayed in Figure 2.
The light curves presented in this work represent those 136 progenitors (of the suite of 200) that both successfully launch an explosion (Section 2.2) and have light curves that would be identified as a SNe IIP, which we find is simply a mass cut of M ZAMS ≤ 31M .

Correlations
We are interested in uncovering correlations between observable properties of the explosion and properties of the progenitors. The size and fidelity of the sample allows us to address these connections necessary to understand light curve diversity. Our robust treatment of the explosion physics combined with large sample of progenitors makes us uniquely situated to address correlations in a novel way. We proceed similarly to Warren et al. (2020), wherein the correlations between observed neutrino and GW signals with progenitor properties were addressed.
We measure correlations with the Spearman's rank correlation coefficient. The Spearman correlation coefficient measures any monotonic relationships between variables, in contrast to the Pearson coefficient which measures only linear correlations. It is important that we are able to access nonlinear relationships that are seen in the data. The combined effect of a wide range of stellar progenitors with mass loss effects and non-linear, non-monotonic explosion energetics over the range of progenitors produces robust and realisticbut not necessarily linear -relationships.
The Spearman coefficient is obtained by first ranking the data by replacing the values by their indices after sorting. For example, the data (1.5M , 1.4M , 1.6M ) would transform to (2, 1, 3). Then, the Spearman rank correlation coefficient is obtained by computing the Pearson correlation of the transformed data, calculated by for ranked variables x and y withx andȳ being the mean values. This process of first ranking the data is what allows the Spearman process to produce a more robust correlation metric. We note that the above equation is the same as that for the Pearson correlation coefficient and, when used on nonranked data, will produce the Pearson correlation coefficient. A value of +1(-1) represents an exact monotonic correlation (anticorrelation) and a value of 0 indicates no monotonic relationship. We consider values |ρ| 0.5 to indicate strong statistical correlation, values 0.3 |ρ| 0.5 to be moderate correlation, and |ρ| 0.3 to be a weak correlation, as is standard practice. Correlation coefficients were calculated using Python's scipy.stats.spearmanr package.
We strived to limit observables considered to those reasonably detectable with current facilities -mostly photometric and early time (meaning, in this context, on the plateau but not requiring observations within days of explosion) features. Ultimately plateau duration, plateau luminosity, and ejecta velocity -all at early times -proved to be the most useful and accessible. We explored numerous properties of the progenitors for correlations with observables -shell masses, density structures, core compactness, and envelope mass to name a few. Most of these parameters had weak relationships with observable properties. Ultimately, we settled on the mass of the iron core as the most meaningful and useful progenitor property, as we will see in the next section.

Light curve fitting
A common method for estimating CCSN progenitor properties is to construct a grid of models with varying masses, explosion energies, and 56 Ni masses and distributions and select the progenitor from that grid that best fits an observed light curve (see Morozova et al. 2018;Martinez & Bersten 2019;Martinez et al. 2020, for recent examples). We accomplish this by finding the progenitor which minimizes the average relative error ε of a quantity f where f * t * is the observed quantity at time t * , f t * is the synthetic quantity at the same time, and N is the number of observational data points. We compare the synthetic and observed data only at the times where observational data is available, using the closest synthetic data to the observational data, which is always within 0.02 days with the output frequency used with SNEC. That is, we do not interpolate between observational data points. We not not consider uncertainties in the explosion epoch in the current work. We seek models that match both observed bolometric luminosity and velocity evolution, i.e., we seek a model minimizing the combined error metric ε(L 50 ) + ε(v Fe ).
Other approaches have been used, such as (historically) simply fittng by eye, χ 2 minimization (Morozova et al. 2018), and Markov chain Monte Carlo methods (Martinez et al. 2020). We implemented several minimization approaches and found that the above method worked best for the current work. This is discussed more in Section 3.2.2.

RESULTS
We consider the properties of the bolometric light curves followed through the end of the plateaus and into the radioactive tails and the ejecta velocities for models with ZAMS masses 9M ≤ M ZAMS ≤ 31M for a total of 136 progenitors. In an effort to find relationships with observables that are easily detectable, we consider primarily the photometric and spectroscopic properties in the plateau phase. The primary quantities that we consider are the plateau luminosity at day 50 (L 50 ), the plateau duration (t p ), and the ejecta velocity at day 50 (v 50 ). These quantities are commonly used when inferring explosion properties from observations (e.g., Litvinova & Nadezhin 1985;Popov 1993;Pejcha & Prieto 2015) and so their trends from realistic models are of particular interest. These quantities are easily detectable by current and next generation facilities without the need for late time observations or particularly high cadences, acknowledging that the photosperic velocity will not be as easily observable for most sources. This will allow for a relationship to be obtained between these quantities and properties of the core of the progenitor that is both robust and easily detectable with standard measurements.

Landscape properties across ZAMS mass
Here we present global trends in photometric properties to test the impact of our explosion calculation on light curve features. As we will see, these properties exhibit nonmonatonic features as a function of ZAMS mass and thus introduce degeneracy into attempts to infer progenitor properties from direct comparisons to light curves. Figure 4 shows the bolometric luminosity at day 50 (on the plateau for all progenitors) for all masses. The imprint from the distribution of explosion energies is readily seen in the plateau luminosities, with more energetic explosions yielding brighter plateaus. A consequence of this is the highly degenerate mapping between plateau luminosity and ZAMS mass following the explosion energy distribution (Figure 2). Figure 5 shows the plateau duration for the STIR + SNEC models. We follow Valenti et al. (2016) and Goldberg et al. (2019) and compute the plateau duration by fitting part of the light curve near the end of the plateau to a combined Fermi-Dirac -linear function of the form This avoids biases or inconsistencies that may possibly be introduced by determining the plateau durations by eye for a large sample of light curves. The physical significance of the various fitting parameters is described in detail in Valenti et al. (2016) and Goldberg et al. (2019). Importantly, the parameter t p is taken to be the plateau duration and tends to be placed about halfway through the drop off of the plateau. Also of interest are a 0 and w 0 which describe the luminosity drop at the end of the plateau and the width of the drop, respectively. Fitting was done using Python's scipy.optimize.curvefit package starting shortly before the end of the plateau. For a few of the high mass models between 27 and 28M , timestep restrictions made it difficult to simulate the explosions into the radioactive tails. Most made it to the end of the plateau and began drop off, but two progenitors were unable to reach the end of the plateau. For the former case, the fitting is unable to work properly and the plateau duration is set by hand in a way that was consistent with the fitting routine. For the two progenitors that could not reach the end of the plateau -27.4M and 27.5M -we omit them in comparisons involving the plateau duration.
Clearly, the distribution of the explosion energies imparts a resulting morphology on the plateau durations that cannot be reproduced without energetics informed by neutrino-driven explosions. We note that many of the plateaus here are quite long, greater than 150 days or so, which is not very common. These plateaus originate from very massive progenitors, around 20M , which are rare in nature. Moreover, these models retain quite massive H-rich envelopes (see Figure 1) and have reduced explosion energies (see Figure 2). The combination of massive H-rich envelope with reduced explosion energy results in extended plateaus (Popov 1993). Some uncertainty in the plateau duration remains through the prescription for setting the mass and mixing of radioactive 56 Ni, as it lengthens the plateau slightly (Kasen & Woosley 2009;Morozova et al. 2015;Sukhbold et al. 2016;Goldberg et al. 2019;Kozyreva et al. 2019). These uncertainties, however, should be on the order of days (see, e.g., Figure 13 2019)). It is also somewhat difficult to fairly compare plateau durations to observed works, as many authors present the length of the optically thick phase duration (e.g., Gutiérrez et al. 2017b) which may be smaller than our measurement by another 5-10 days or more. For these reasons, we defer further comparisons to observational data of the plateau durations to future work.
All of this directly impacts the ability to reliably extract progenitor features from light curves. Without a distribution of explosion energies that is set by a physically realistic explosion model, any sort of arbitrary distribution of light curve properties may be recovered, even with the same diversity of progenitors used. While STIR is not a perfect or parameter free description of the explosion -no 1D model ever will be -it matches well with 3D results and provides a large set of such physically motivated explosion energies for these studies.
Another quantity of interest -albeit not a directly observable one -is the time to shock breakout. Figure 6 shows the time for the shock to breakout from the stellar surface for the STIR + SNEC models. This is particularly important, as the time to shock breakout sets the on source window for electromagnetic follow-ups of gravitational wave and neutrino events from core-collapse supernovae (Abbott et al. 2020). The time to shock breakout is sensitive to the structure of the  progenitor and the explosion energy and may be significantly over-or under-estimated if an incorrect explosion energy is used.
With the next galactic CCSNe and prospects for detecting their gravitational wave and neutrino signals, the time to shock breakout becomes a measurable quantity through the difference between GW or neutrino detection time and first light from the SNe. The SuperNova Early Warning System (SNEWS) (Adams et al. 2013;Kharusi et al. 2021) will alert observatories to trigger an EM followup after a neutrino detection, and knowing the shock breakout time will be an important factor for the followup study. Combined with constraints from the GW detection (Abbott et al. 2020) and constraints from other EM observations, the time to shock breakout could help to place additional constraints on the SNe progenitor -provided that adequate energetics are used. Similarly, constraints on the shock breakout time after an EM signal may be used to look back at GW and neutrino data, assuming a nearby event.
The previous figures highlight the strong dependence on the distribution of explosion energies used to drive the explosion. This leads to degeneracies when mapping from observables to ZAMS mass with many progenitors of varying masses being capable of producing a given observation.

Comparisons with observations
In this section, we compare our light-curves to observations of SNe IIP both through global properties of many SNe and fits to the light-curves of individual SNe IIP that have M ZAMS determined through pre-explosion imaging data.

Comparison with a large observational sample
Apart from photometric observations, spectroscopic observations may also be used to constrain progenitor properties. While we have not computed full synthetic spectra in this work, we can approximate standard line velocities. Figure 7 shows ejecta velocity at day 50 (v 50 ) versus plateau luminosity at day 50 (L 50 ) for all progenitors that exploded as SNe IIP. Also plotted are data presented in Gutiérrez et al. (2017a,b) 3 . All ejecta velocities here are inferred from the Fe II (5169Å) line. In our models, this velocity is calculated in post-processing as the velocity of the ejecta at the point where the Sobolev optical depth (τ Sob ) is unity, with where q e and m e are the electron charge and mass, n F e is the number density of iron atoms, η i is the ionization fraction relevant for the transition of interest, f = 0.023 is the atomic oscillator strength, t expl is the time since explosion, and λ 0 is the wavelength associated with the transition. For material in homologous expansion, this measures the strength of a particular line (Mihalas et al. 1978;Kasen et al. 2006) and the point where τ Sob = 1 has been shown to match better to observational measurements than the τ = 2/3 electron scattering photosphere ( 2018). To estimate the ionization fraction, we use a table of η i as a function of density and temperature that is now publicly available in MESA (Paxton et al. 2018). We choose to use this metric for the velocity evolution because ultimately we seek to compare with observables. While the standard τ = 2/3 photosphere -and its velocity -are simple to compute, they are not simple to observe. On the other hand, the Fe II 5169Å line is commonly measured. Therefore, we seek to estimate the location in the ejecta where this line is measured, using the Sobolev approximation that has been readily used in recent works (Paxton et al. 2018;Goldberg et al. 2019;Martinez et al. 2020). However, this approach to estimating the iron line velocity is ultimately an approximation and there are physical uncertainties associated with this method. Paxton et al. (2018) investigated the effects of the choice of the Soboloev optical depth used and found relatively small differences when compared to using the traditional photospheric velocity. In lieu of full spectral calculations this method provides an estimate of the desired velocity but more work may be needed to robustly compare to observed ejecta velocities.
The sample of luminosities and velocities from our models matches well with the observational sample, but reach higher in luminosity than the observed set. These high luminosity events are from some of the higher pre-supernova mass stars around the transition to the mass-loss dominated regime (see Figure 1). These high mass stars are less common than their lower mass companions. The highest ZAMS mass stars, those 23M dominated by mass-loss, dip back down and left in luminosity-and velocity-space, obtaining similar luminosities but slightly lower velocities than lower mass pro-genitors. Ultimately, we are able to reproduce observed distributions quite well without having to tune to observations, instead following the explosions from self-consistent simulations.

Determination of progenitor properties for individual events
It is commonplace to estimate supernova progenitor parameters using a grid hydrodynamical models (i.e., codes similar to SNEC using a thermal bomb) with varying initial masses, thermal bomb energies, and other parameters, and determining the best fitting model (see, e.g., Utrobin & Chugai 2008, 2009Pumo et al. 2017;Morozova et al. 2018;Martinez & Bersten 2019;Martinez et al. 2020;Eldridge & Xiao 2019). We attempt to match our set of explosions with 7 observed bolometric light curves from Martinez & Bersten (2019); Martinez et al. (2020) 4 . Bolometric luminosities are calculated using the bolometric correction method of Bersten & Hamuy (2009), which requires only BVI photometry to estimate the bolometric correction. Figures 8 and 9 show observed bolometric light curves (left) and velocity evolution (right) for (top to bottom) SN 2004A, SN 2004et, SN 2005cs, SN 2008bk, SN 2012aw, SN 2012ec, and SN 2017eaw. Dark blue lines show bolometric luminosity and velocity evolution for best fit progenitors from our sample using the STIR + SNEC model using the fitting described in Section 2.5. For the velocity evolution, dashed lines show approximate Fe II λ5169Å line velocities estimated through the methods described in Section 3.2 and solid lines show the proper τ = 2/3 photospheric velocity. Gold lines are for ZAMS mass models corresponding to estimates from pre-explosion imaging. We use the ZAMS mass estimates from Davies & Beasor (2018) for SN 2004A, SN 2004et, SN 2008bk, SN 2012aw, and SN 2012ec. Properties of these SNe are discussed in detail in Martinez & Bersten (2019) and Martinez et al. (2020). For SN 2005cs, Davies & Beasor (2018) estimated an initial mass of about 7 M -well below the minimum mass we consider to produce a CCSNso we use the estimate from Smartt (2015). Finally, we use the mass estimate for SN 2017eaw from Eldridge & Xiao (2019). In all cases we use the optimal value of the initial mass when possible, or the closest value within the reported range that was both on our mass grid and produced an explosion.
We determine the best fit progenitor by minimizing the total relative error of both luminosity and velocity across the entire light curve after day 30 as discussed in Section 2.5. We also tried minimizing χ 2 , as was done in Morozova et al. (2018), but found unsatisfactory performance compared to our method (see Appendix B for an example using SN2017eaw). We did not consider the errors associated with the observations in our fitting. The inverse variance weighting typically used in χ 2 minimization gave stronger significance to the radioactive tail, as this region has much smaller error compared to the plateau. The result was the selection of models that fit the tail nicely, but fit the plateau very badly. We do not consider data before 30 days post shock breakout, as very early time bolometric luminosities may be heavily influenced by interactions with circumstellar material (CSM) for some SNe (Morozova et al. 2018) and we have not included CSM effects in this work.
We do not expect to find close fits for all observed CCSNe. In this work, we have progenitors that cover a wide range of ZAMS masses with explosions driven by turbulence-aided neutrino radiation hydrodynamics simulations, but are limited in scope in other regards, such as rotation, metallicity, 56 Ni mass and distribution, and possible effects of binarity. Moreover, we do not have models with masses lower then 9M , which may contribute to CCSNe. For example, SN 2008bk is very underluminous with low expansion velocities as is very likely a lower mass progenitor than we have in our set (Mattila et al. 2008;Van Dyk et al. 2012;Lisakov et al. 2017;Lisakov et al. 2018;Martinez & Bersten 2019;O'Neill et al. 2021). With these limitations in mind, we still find good fits for two observed CCSNe, notably 16.0 M for SN 2012aw and 20.0 M for SN 2017eaw. Our best fit progenitors tend to have larger ZAMS masses than those estimated from pre-explosion imaging, for example by about 5M for SN2017eaw. This difference of about 5M is not uncommon - Goldberg & Bildsten (2020), for example, find a possible ZAMS mass for SN2017eaw of 10.2M -also about 5M from the value obtained from pre-explosion imaging. We have presented light curves for which we do not find particularly good fits for the sake of completeness and to show the strengths and weakness of the current progenitor set. As previously mentioned, there is no reason for this progenitor set to perfectly fit any specific light curve.
The differences highlighted in Figures 8 and 9 show the inherent degeneracy involved in extracting CCSNe progenitor properties. As shown in Goldberg et al. (2019); Dessart & Hillier (2019), there are familes of progenitor properties that can lead to a given light curve. This further highlights that light curve fitting is extremely degenerate -not only in the ways explored in previous works, but also in the method used to drive the explosion. Thus, we do not claim that these progenitors necessarily reflect the true progenitors, they simply match the observations given a set of neutrino-driven explosions. It has become clear that more work in needed to infer progenitor properties. Matching an observed SN is a necessary, but not sufficient, condition for inferring progenitor and explosion properties.
Finally, we summarize our best fitting models for those light curves for which we see good agreement (SN2012aw and SN2017eaw) in Table 1 alongside various other sources.

Correlations
In this section, we address the primary goal of this study, which is to connect light curve properties to progenitor properties using a statistically significant sample of simulations. Figure 10 shows the Spearman's correlation matrix for the observable quantities and progenitor properties that we consider for the STIR + SNEC models. Our goal is to assess direct correlations between individual quantities, and for this reason we do not consider correlations with ZAMS mass because it does not correlate with any single quantity. In many cases, we are simply recovering well-known correlations, which provide a sanity check on our methods. For example, relationships between ejecta velocity and luminosity have been used in SNe IIP supernova cosmology (Hamuy 2005;Nugent et al. 2006;Poznanski et al. 2009). Relationships between photometric and spectroscopic observables, L 50 , v 50 , and t p , and properties of the progenitor, such as R 500 (the pre-supernova progenitor radius in units of 500R ) in addition to the explosion energy are used in scaling relationships, such as those in Popov (1993) We first consider some typical observables of SNe IIP light curves -the plateau luminosity (L 50 ), plateau duration (t p ), and ejecta velocity measured through the Fe 5169Å line during the plateau phase (v 50 ). These observables correlate with each other and are expected to correlate with properties of the progenitors, such as the presupernova radius (R 500 ) and envelope mass (M env ). We observe significant correlations between t p , L 50 , and R 500 . Correlations with R 500 tend to be non-monotonic (see, e.g., Figure 1), which is why they tend to have weaker values of the correlation coefficient. There is a moderate correlation between the L 50 and v 50 and the presupernova mass (M preSN ).
The explosion energy (E expl , see Section 2.2) is expected to correlate with both progenitor properties and observable properties. Correlations between E expl and observable properties are monotonic relationships (i.e., always increasing or always decreasing, but not necessarily linear), for example with a correlation coefficient of 0.97 for L 50 -E expl . This is because in the self-consistent STIR + SNEC models, the explosion energies are the total positive energies of unbound material as liberated by neutrino heating and is thus correlated with properties of the core (and thus, the rest of the progenitor properties through stellar evolution) of the progenitor.
Finally, we turn our attention to connections between properties of the core of the progenitor and observable quantities. Motivated by connections between explosion energy    Here we consider the following quantities: iron core mass (MFe), progenitor radius (R500), explosion energy (E expl ), ejecta velocity at day 50 (v50) as determined from the Fe II (5169Å) line, log of the plateau luminosity at day 50 (L50), and plateau duration (tp). The lower left half of the matrix shows the Spearman rank correlation coefficient for each pair of quantities. and the compact remnant, we explore correlations with the iron core mass (M Fe ). Progenitors with more massive iron cores tend to liberate more gravitational binding energy, have higher neutrino luminosities, and ultimately are associated with more energetic explosions for progenitors that successfully explode. The origins of this correlation can be seen in the bottom panel of Figure 2 through the connection between iron core mass and explosion energy. This correlation, therefore, once again highlights the need for realistic physics in explosion models even in 1D. Equipped with this correlation, and the previously mentioned relationships between explosion energy and observables, one might expect some imprint of the iron core mass on the observables. Indeed, for the STIR + SNEC models we observe a very strong, linear relationship between iron core mass and plateau luminosity. We note that the compactness parameter ξ 2.5 (O'Connor & Ott 2011) produces a stronger correlation. This, however, is of little practical use, as the 9-12M progenitors have nearly zero values of the compactness parameter (≤ 0.02), breaking the trend for the most common progenitors, and the iron core mass is a more physical quantity (i.e., does not depend on the exact choice of mass coordinate for the measurement). The compactness parameter and iron core mass are very tightly correlated and both provide a measure of the gravitational binding energy available in the explosion. A relationship between iron core mass and supernova observables helps constrain stellar evolution models and characterize the diversity of supernova light curves. Figure 11 shows iron core mass versus plateau luminosity at day 50. Higher luminosity events tend to originate from progenitors with more massive iron cores. Ultimately, more massive stellar cores collapse to form more massive proto-neutron stars, liberating more gravitational binding energy in the process and resulting in higher neutrino luminosities emanating from the PNS surface. All of this results in a more energetic explosion and a brighter supernova. In Table 2 we report the fits coefficient for the M Fe -L 50 relationship and the associated variances and covariances for a linear fit of the form y = ax + b.
This correlation, though simple, has a profound implication that we can constrain core structure from optical photometry alone. While not necessarily providing a precise measure of the iron core mass for individual events due to observational error and uncertainties on the fit parameters from scatter, which we quantify below, it provides a method for comparing the cores of virtually all SNe IIP simultaneously. Furthermore, these parameter estimates can be used to constrain stellar evolution models for CCSN progenitors. We find a similar, although slightly weaker, correlation between the ejecta velocity at day 50, v 50 , as well, but most LSST sources won't have a spectroscopic follow-up so this is of limited use.
For the case of Sukhbold et al. (2016), fewer massive iron cores produced explosions, and the explosions had a tendency to be brighter. Using their data we find slope and intercept parameters of 0.033 and 1.344, respectively.
For any relationship of this type to be useful, error must be taken carefully into account. The optimal fit parameters were obtained with a least squares method. However, it is known that the covariances provided by least squares methods are not appropriate for a wide range of problems, including those with a non-Gaussian intrinsic scatter among other criteria (see, e.g., Clauset et al. 2009, and references therein). For this reason, we resort to a bootstrapping method (Efron 1979) to obtain the errors on the fit parameters. This method has the advantage of making no assumptions about the underlying distribution of the data. Instead, bootstrapping operates by resampling the data M times with replacement. For each resampling, a new fit is made and those fit parameters stored. Then, estimates of the variance and covariance of parameters u and v are given by where u and v are the optimal fit parameters and each of u j , v j are the fit parameters for each of the M resamples. These error estimates tend to be, for this application, somewhat smaller than parameter errors obtained through a simple least squares method. The full set of fit parameters, variances, covariances, and adjust coefficient of determination are supplied in Table 2. We note that the fit presented uses the non-log plateau luminosity as its independent variable, as opposed to the log luminosity presented in other parts of the paper. Then, given errors on the fit parameters it is straightforward to compute the error on an iron core mass estimate. For a linear fit, we propagate the combined observationalfit parameter uncertainty in the following way: where L 50 is the luminosity at day 50 in erg s −1 and where we have included explicitly the covariance of the fit parameters a and b. In order to further account for intrinsic scatter in the relationship, we have included σ res which is the 67% percentile on the residual distribution r i = |M Fe −M Fe |, wherê M Fe is computed from the fit.
As an example, we estimate iron core masses for six wellobserved SNe, shown in Table 3. Data for SN1999em, SN2003hl, and SN2007od are taken from Gutiérrez et al. (2017a,b). Data for SN2004et, SN2012aw, and SN2017eaw are taken from Martinez et al. (2020).

DISCUSSION AND CONCLUSIONS
We present synthetic bolometric light curves for 136 solar metallicity, non rotating CCSNe progenitors and consider statistical relationships for those with ZAMS masses ranging Table 2. Linear fit paramaters for iron core mass (MFe) to plateau luminosity (L50)in units of 10 42 erg s −1 . The first two rows shows the optimal fit parameters. The next two rows shows the error on each parameter. The next row shows the covariance between the parameters and the residual error accounting for intrinsic scatter. The final row shows the adjusted coefficient of determinationR 2 for the fit.  Couch et al. (2020) as the initial condition. This allows for light curves obtained without a user-set explosion energy. Our 56 Ni yields were fit from Sukhbold et al. (2016) who exploded the same progenitors with an expansive reaction network coupled to the evolution. This is sufficient for the current work, and future work with FLASH will include detailed nucleosynthesis calculations. These light curves, as well as the SNEC initial profiles and necessary parameters, are provided online 5 . We also include the necessary binding energy of our progenitors to correct STIR's explosion energy to produce identical results with a thermal bomb explosion. In the online resources, we furthermore provide the light curves for the M ZAMS > 31M models that successfully explode. For progenitors that explode with STIR, we follow the explosions in SNEC to produce bolometric light curves, forming a large, statistically significant set of CCSN light curves followed from high-fidelity explosions allowing us to address relationships between progenitor properties and properties of the explosion in a statistical way. We consider the full shape of these light curves, but also reduce them to characteristic quantities such as the plateau luminosity, plateau duration, and ejecta velocity. Next, we show that global trends in light curve properties -such as plateau duration and plateau luminosity -depend sensitively on the explosion model and require explosion energies set by robust physics. To demonstrate this, we compute bolometric light curves for the same set of progenitors using two different thermal bomb models with SNEC. The distribution of explosion energies plays a leading role in setting the distribution of observables across a large sample of progenitors. Thus, the ability to identify global trends in light curve properties and extract progenitor features from them depends sensitively on the determination of explosion energy, underscoring the need for explosions driven with highfidelity multi-physics models.
We present a simple best-fit procedure to individual, observed CCSN light curves (Martinez et al. 2020). The usual procedure for estimating progenitor properties of observed CCSNe is to construct a large grid of "hydrodynamical models" -usually in ZAMS mass, explosion energy, and perhaps 56 Ni mass and distribution -and find a best fit model. This approach results in known degeneracies, for example, as shown by Goldberg et al. (2019); Dessart & Hillier (2019) wherein there are certain families of progenitor and explosion parameters (such as ejecta mass, explosion energy, and ejecta velocity) that produce a given light curve, though preexplosion radius measurements may help to resolve this degeneracy (Goldberg & Bildsten 2020;Kozyreva et al. 2020). Our approach differs in that we do not control the explosion properties, instead following a dense set of various ZAMS mass progenitors from neutrino driven explosions. While this does not solve the light curve degeneracy problem, it could reduce the size of the family of explosion properties for a given light curve, as some combinations of explosion energy and stellar mass are not realizable. Although the explosions are not calibrated to observed data we still find great agreement both when comparing to large samples of events and for some individual cases. Intriguingly, we find best-fit ZAMS masses that are greater by as much as ≈7M than those estimated from pre-explosion imaging in tandem with stellar evolution modeling. The fact that hydrodynamic models have tended to find ZAMS masses in agreement with preexplosion imaging estimates for these CCSNe (Morozova et al. 2018;Martinez & Bersten 2019;Martinez et al. 2020) may indicate the danger of exploring too large a parameter space instead of knowing which regions are physically realizable, though we note that some hydrodynamic models have also found noticeably higher masses in better agreement with our conclusions (e.g., Utrobin & Chugai 2008, 2009). Ulti-mately, the set of solutions for matching a given observed light curve is degenerate, with many progenitors being capable of producing a given light curve.
Despite the progenitors and explosions in this study not being crafted to reproduce specific events, we find good qualitative agreement with SN2012aw and SN2017eaw. Notably, the luminosity evolution of SN 2012aw is fit by our 16.0M progenitor remarkably well. The best fit progenitors for the observed light curves in this study are not necessarily the progenitors that these explosions originated from -they simply reproduce the observables. We have demonstrated that beyond the now understood light curve degeneracies, there are additional degeneracies inherited from the choice of explosion model. This result is complementary to the recent findings by Farrell et al. (2020) where they showed that a star's final temperature and luminosity cannot be reliably traced back to the star's ZAMS mass -that very different mass stars may end up at the same temperature and luminosity. These results together show that much more work is needed before a SN IIP progenitor's ZAMS mass can be reliably determined -the path from stellar birth to death is not a one-to-one function.
The light curves here present avenues for future work to explore the discussion surrounding explosion energy. There is tension between explosion energies realized in 3D CCSN simulations and energies inferred from fitting hydrodynamical models to observations. The energies from these two methods differ, with those inferred from hydrodynamical modeling being significantly larger (see Murphy et al. 2019, which discusses this tension in detail). On one hand, 3D simulations of very massive progenitors have often simply not asymptoted to their final values within the simulated time. There is also still physics left to include, such as the recently demonstrated affects of magnetic fields on neutrinomatter interactions (Kuroda 2021) and improved neutrino pair-production rates (Betranhandy & O'Connor 2020) on the explosion mechanism, neutrino mixing, among other affects, all of which will likely play a role in setting the final energy. On the other hand, solutions using thermal bomb models have been shown to be degenerate, and these studies access a very large area of this degenerate parameter space and may not necessarily find physically realizable solutions. The methods described here could illuminate or even weaken the tension between these energies by limiting the parameter space spanned by hydrodynamical modeling studies and by using physically-motivated explosions.
The final aim of this study is to leverage the large number of light curves to perform a statistical investigation of relationships between progenitor and explosion properties. Focusing our investigation to SNe II, 136 light curves, a number of correlations between the light curves and their progenitors are found. We find a robust relationship between the iron core mass of the progenitor and the luminosity on the plateau of the SNe. This relationship allows one to, for the first time, constrain properties of the stellar interior from photometry alone. We provide an analytic approximation to the observed correlation, including error, for future use with large survey data such as LSST.
Recently, Curtis et al. (2020) presented synthetic light curves and spectra from a sample of 62 CCSNe with the 1D PUSH model (Ebinger et al. 2017) and SNEC to obtain the light curves. Our results complement one another in several ways -notably, the size and composition of our samples differ. Our sample contains 148 light curves -136 of which are analyzed in this work -from the same metallicity, compared to their 62 light curves from three different metallicity populations ranging from zero to solar. This allows us to more robustly survey global explosion properties of progenitors from similar origin within the nearby universe. These studies, together, survey a vast range of progenitor properties. The CCSNe simulations in our work are perfmormed with FLASH using the STIR model. Notably, STIR requires no tuning to observations, eliminating the potential for biases when simulating progenitors different than the one used for tuning. Importantly, the results from STIR are consistent with 3D simulations. The explosion energies, explodability, and the shape of each as a function of ZAMS mass differ non-trivially for STIR and PUSH (see Couch et al. 2020;Ebinger et al. 2019) and this could impact global trends in explosion properties. On the other hand, Curtis et al. (2020) obtained their 56 Ni distributions using a nuclear reaction network in conjunction with their CCSN simulations. As aforementioned, we estimated 56 Ni mass from the explosion energy, informed by KEPLER yields. Curtis et al. (2020) also have a larger diversity of supernova types through their inclusion of sub-solar and zero metallicity progenitors. To keep the scope of the current work contained, we have not produced synthetic spectra for these explosions, whereas Curtis et al. (2020) calculated spectra for their supernovae.
Similarly, Sukhbold et al. (2016) present a sample of synthetic light curves of the same statistical size and originating from the same progenitors using a different parametrized, neutrino-driven explosion mechanism. Using these simulations they present scaling relations to determine explosion and progenitor properties from observables. The outcomes of these simulations -both the explosions and resulting light curves -differ from STIR and this work, having a tendency to be brighter than those produced in this work. It would be interesting, for future work, to investigate the affect of these differences in explosion mechanism when applied to populations of observed CCSNe and implications for inferred properties such as explosion energy.
This work is part of a larger context to understand and predict full multi-messenger signals from realistic CCSNe.
Understanding how variations in progenitors properties tie into variations of different observables will ultimately help to constrain real populations. This work, in tandem with the work of Couch et al. (2020) and Warren et al. (2020), gives us explosion fates, energies, neutron star mass distributions, neutrino signals, approximate GW signals, and now EM signals for a massive suite of neutrino driven CCSNe. It is only through advanced methods -studying in detail all messengers from first principles simulations -used in tandem with growing observational data that we can truly understand these phenomena.
Software: FLASH (Fryxell et al. 2000;Dubey et al. 2009), SNEC (Morozova et al. 2015(Morozova et al. , 2016 Nuclear Security Administration) that are responsible for the planning and preparation of a capable exascale ecosystem, including software, applications, hardware, advanced system engineering, and early testbed platforms, in support of the nation's exascale computing imperative. This work was supported in part by Michigan State University through computational resources provided by the Institute for Cyber-Enabled Research. We collectively acknowledge that MSU occupies the ancestral, traditional and contemporary lands of the Anishinaabeg -Three Fires Confederacy of Ojibwe, Odawa and Potawatomi peoples. APPENDIX A. COMPOSITIONAL DEPENDENCE For our light curves, we modified the compositional profile in the FLASH part of the domain to be pure 4 He, as full composition is not currently tracked in the output. In this appendix, we provide comparisons of select light curves using thermal bombs with FLASH explosion energies for both the modified compositional profile and the original compositional profile. Figure 12 shows light curves with the unaltered (orange) and modified (blue) compositional profiles for 9, 15.2, 25, and 30M progenitors. For the cases considered here, the difference in luminosity on the plateau is bounded above by 0.1 dex, which has no meaningful affect on the iron core mass estimates and distributions of Section 3.3. 30.0 Figure 12. Light curves using a thermal bomb driven explosion with STIR explosion energies using the modified compositional profile (blue) and unaltered profile (orange). We show light curves for 9, 15.2, 25, and 30M progenitors.