NUMERICAL MODELING OF THE EARLY LIGHT CURVES OF TYPE IIP SUPERNOVAE

One of the outstanding problems in astrophysics is connecting the variety of core-collapse supernovae (SNe) we observe with the massive progenitors that give rise to them. Ideally, we would use pre-explosion imaging to directly identify these progenitors (e.g., Li et al. 2006; Smartt et al. 2009; Van Dyk et al. 2012). Unfortunately, in most cases such information is not available because the progenitor is too dim or deep pre-explosion imaging is not available at the needed location. For this reason, features of the early light curves can be especially helpful in constraining progenitor properties. Emission dominated by shock cooling should reflect the radius and the density profile of the exploding star (see, e.g., Grassberg et al. 1971; Falk & Arnett 1977; Soderberg et al. 2008; Nakar & Sari 2010, hereafter NS10; Rabinak & Waxman 2011; Bersten et al. 2013). The presence of extended material around the progenitor influences both the photometry and spectroscopy of the early SN light curve (see, e.g., Chevalier & Irwin 2011; Moriya et al. 2011; Ginzburg & Balberg 2014; Nakar & Piro 2014; Svirski & Nakar 2014; Piro 2015 for the theoretical view and Ofek et al. 2010; Garnavich et al. 2016; Khazov et al. 2016 for observational evidence).

The early light curve evolution has been explored in a number of theoretical works semi-analytically (e.g., NS10; Rabinak & Waxman 2011). These studies generally took the approach of assuming an idealized (i.e., polytropic) density profile for the star to make connections between early light curve properties and the properties of the star in a general way. The question, though, is whether in nature such idealized profiles actually occur at relevant depths within the star. On the other hand, numerical models represent a powerful complementary approach to analytic studies of SN light curves, including the early phases (see Eastman et al. 1994; Young 2004; Kasen & Woosley 2009; Bersten et al. 2011; Dessart et al. 2013). For example, Bersten et al. (2012) constrained the radius of the progenitor star for SN 2011dh based on the numerical modeling of its early light curve. In this case, many of these calculations are tailored for specific events. This makes it difficult to generalize these results more broadly, which would be especially useful as current and future transient surveys (e.g., iPTF/ZTF, Rau et al. 2009; Bellm 2014; ASAS-SN, Shappee et al. 2015; BlackGEM, Bloemen et al. 2015; Pan-STARRS, Huber et al. 2015; Schlafly et al. 2012; LSST, Ivezic et al. 2008) find larger samples of SNe.

Motivated by these issues, we undertake a numerical exploration of how rise times of SNe IIP vary with the parameters of their progenitors. We utilize two sets of nonrotating solar-metallicity progenitor models from the stellar evolution codes MESA (Paxton et al. 2011, 2013, 2015) and KEPLER (Weaver et al. 1978; Woosley & Heger 2007, 2015; Sukhbold & Woosley 2014), and explode these models and generate light curves using the SuperNova Explosion Code (SNEC; Morozova et al. 2015). We describe important differences between realistic stellar models and more idealized treatments, and how these impact observable features of the early light curves. In addition, based on the results of our numerical models, we derive a relation for the rise time ${t}_{{\rm{rise}}}$ of SN IIP light curves in the g band as a function of the radius R of the progenitor. This relation may be useful in future analyses of observed early rises, especially when comparing large samples of events (e.g., Gall et al. 2015; González-Gaitán et al. 2015; Rubin et al. 2016; Valenti et al. 2016). Our results show that the observed rise times are in agreement with the observed red supergiant (RSG) radii from Levesque et al. (2005, 2006), and there is no need to restrict the progenitor models to small radii $\lesssim 500\,{R}_{\odot }$, as was suggested in some previous studies (see, e.g., Dessart et al. 2013; Shussman et al. 2016).

The paper is organized as follows. In Section 2, we describe the progenitor models used in this study and the numerical setup of our simulations. In Section 3, we discuss two important aspects of our calculations related to the stellar models and radiative diffusion. This highlights differences with more idealized treatments. In Section 4, we present our full set of explosion models and summarize how the properties of the light curve rise relate to the radius and ejecta mass. In Section 5, we summarize our main results and discuss future work.

This work is primarily focused on SNe IIP, and thus we only consider RSG progenitor models since these are the observationally confirmed SN IIP progenitors⁵ (see, e.g., Smartt et al. 2009; Fraser et al. 2012; Maund et al. 2013). We consider two sets of nonrotating solar-metallicity pre-collapse RSG models. The first set of models comes from the stellar evolution code KEPLER (Weaver et al. 1978; Woosley & Heger 2007, 2015; Sukhbold & Woosley 2014; Sukhbold et al. 2016) and has zero-age main sequence (ZAMS) masses in the range between $9\,{M}_{\odot }$ and $20\,{M}_{\odot }$ in steps of $0.5\,{M}_{\odot }$. We refer the reader to Sukhbold et al. (2016) for a detailed description of these models. The second set of the models we generate⁶ using the stellar evolution code MESA (Paxton et al. 2011, 2013, 2015), revision 7624. These models have ZAMS masses in the range between $11\,{M}_{\odot }$ and $20\,{M}_{\odot }$ in steps of $0.5\,{M}_{\odot }$. Our stellar evolution calculations use the same algorithms and input parameter set used for the unstripped model in Morozova et al. (2015; see their Section 3.2). We do not artificially strip mass from the models presented in this paper because it is not our goal to investigate the impulsive mass-loss events that might occur in the evolution of massive stars. We summarize for completeness the parameters that more directly influence the radius determination. We use the 21 isotope nuclear reaction network approx21.net, and use the Ledoux criterion for convection following Sukhbold & Woosley (2014) for the choice of the free parameters. This corresponds to a mixing-length parameter ${\alpha }_{\mathrm{mlt}}=2.0$, exponentially decreasing overshooting (both for the core and the convective shells) with ${f}_{\mathrm{ov}}=0.025$ and f₀ = 0.05 (see Paxton et al. 2011), and semiconvection efficiency ${\alpha }_{\mathrm{sc}}=0.1$. Wind mass loss is included as in Morozova et al. (2015): we use the Vink et al. (2000, 2001) rate, modified like in Morozova et al. (2015), for the early (blue) phase of the evolution, and the de Jager et al. (1988) rate for RSG mass loss ("Dutch" wind scheme in MESA). We do not apply any modifying efficiency factor to either rate. We set mesh_delta_coeff = mesh_delta_coeff_for_highT = 1.0 for the spatial resolution, and varcontrol_target = 10⁻⁴ for the temporal resolution (see also Paxton et al. 2011, 2013). We note that experiments with MESA show that the pre-collapse RSG radius depends sensitively on wind efficiency (M. Renzo et al. 2016, in preparation), overshooting, and mixing-length parameters. Generally speaking, the higher the wind efficiency, the more mass is removed and the smaller the radius. The larger f_ov, the more massive and luminous the He core and the larger the radius. The larger ${\alpha }_{\mathrm{mlt}}$, the more efficient the energy transport through the envelope and the smaller the radius (see also Dessart et al. 2013). Ultimately, three-dimensional (3D) radiation-hydrodynamic models will be needed to robustly predict RSG radii; e.g., Chiavassa et al. (2009).

Figure 1 summarizes the main features of the considered progenitor models, the radii, and masses at the onset of core collapse. The KEPLER models show a strong correlation between radius and mass. We emphasize that whether or not one set of models is more "correct" (in the sense that it is a closer representation of what actually occurs in nature) is irrelevant to our work. The reason is that we are looking for general trends that are satisfied by the explosion of any stellar model (as we will summarize in Section 4). The fact that so much diversity is seen in Figure 1 is actually a strength and not a weakness for our study.

Zoom In Zoom Out Reset image size

All the stellar models are exploded with SNEC, which is described in detail in Morozova et al. (2015). We excise the inner $1.4\,{M}_{\odot }$ of the models, assuming that this part collapses and forms a neutron star. Later in the paper, the ejecta mass, ${M}_{{\rm{ej}}}$, is defined as the total stellar mass at core collapse minus $1.4\,{M}_{\odot }$. We do not model the fallback of material onto the remnant. For the current study, we use a thermal bomb mechanism for the explosion. We define the final energy, ${E}_{{\rm{fin}}}$, as the final (asymptotic) explosion energy of the model. Note that E_fin is not equal to the energy of the thermal bomb that we use to initiate the explosion. The latter is equal to the difference between ${E}_{{\rm{fin}}}$ and the total (mostly gravitational) energy of the progenitor before explosion. We point out that the absolute value of the initial gravitational energy of our models can be of the same order or even larger than ${E}_{{\rm{fin}}}$. We inject the energy of the thermal bomb in the inner $0.02\,{M}_{\odot }$ of the model (after excision of the inner $1.4\,{M}_{\odot }$) for a duration of $0.001\,{\rm{s}}$ (the choice of bomb duration is discussed in Section 4). "Boxcar" smoothing of the compositional profiles is performed as in Morozova et al. (2015) and the same values for the opacity floor are used. We do not include radioactive ${}^{56}{\rm{Ni}}$ in our models since it does not impact these early phases. We use the equation of state by Paczyński (1983) and solve for the ionization fractions of hydrogen and helium. The numerical gridding of each model is identical to that used in Morozova et al. (2015) and consists of 1000 cells in mass coordinate. However, for a special case described in Section 3, we use a grid consisting of 2000 cells, with increasing resolution toward the surface and the center of the models. All light curves are generated for the first 40 days after explosion.⁷

Before presenting our full set of explosion calculations, it is helpful to discuss some of the key aspects of the stellar models and the properties of the radiative diffusion that determine the rising light curves we find. The two main factors we focus on are the shallow density profile of RSGs and the time-dependent position of the so-called luminosity shell. The location of the luminosity shell is the depth from which photons diffuse to reach the photosphere at a given time after shock breakout. These two issues are very connected, since the density profile determines how the shock accelerates and will eventually set the depth of the luminosity shell as shock cooling sets in.

3.1. Density Profiles of RSGs

RSGs have convective envelopes, which means that nominally their density ρ obeys the power-law $\rho \propto {(R-r)}^{n}$, where r is the distance from the center of the star, R is the radius of the star, and n = 1.5. For this reason, previous analytic studies paid special attention to density profiles with n = 1.5 (see Matzner & McKee 1999; Rabinak & Waxman 2011; NS10). To test this assumption, we plot in Figure 2 $\,{\mathrm{log}}_{10}\,\rho$ as a function of $\,{\mathrm{log}}_{10}(R-r)$ for a representative subset of the models described in Section 2. The rest of the models have a similar structure and they are not plotted for clarity. Black dashed lines show fits of the power-laws $\rho \propto {(R-r)}^{n}$ to different regions of each profile. The indices n₁ and n₂ give the power-law exponents obtained for the outer and inner parts of the profile, respectively. The innermost and outermost boundary of the fits for each model are chosen in such a way that the inner power-law exponent is 1.5 for all models since the bulk of the envelope is clearly convective. Note that the innermost boundaries of our fits do not reach the transition to the helium core, which corresponds to densities $\rho \sim {10}^{-6}\mbox{--}{10}^{-7}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ for different models from our sets. We do not attempt to fit a power law to the innermost parts of the envelopes and the cores, where n is significantly larger than 1.5.

Zoom In Zoom Out Reset image size

From the comparison shown in Figure 2, there are multiple key points to take away. First, indeed at sufficiently large depths within the envelope, the density profile obeys an n = 1.5 polytrope. Second, in both MESA and KEPLER models, the outer density profile is different from n = 1.5 and in general shallower (although at sufficiently short distances to the surface, the KEPLER models are steeper; these regions do not impact the light curves we study). This region can cover a few hundred solar radii. As for the stellar radii mentioned in Section 2, the difference in the outer density profiles of MESA and KEPLER models is controlled by such parameters of stellar evolution as wind efficiency, mixing, and overshooting. Whether or not the shallow profile region has an important impact on the rise depends on the depth where photons are diffusing from during the shock-cooling phase. We address this next.

3.2. Position of the Luminosity Shell

The shock-cooling phase (from a few to ∼20 days) and the plateau phase (from ∼20 to ∼100 days) of SNe IIP light curves are powered by the energy of the post-explosion shock wave. This energy diffuses out of the expanding envelope and at each moment of time we see the photons coming from the luminosity shell. The position of the luminosity shell at time t is defined by the condition ${\hat{t}}_{{\rm{diff}}}=t-{t}_{0}$, where t₀ is the time of shock breakout, ${t}_{{\rm{diff}}}={t}_{{\rm{diff}}}(r,t)$ is the diffusion time at each time at each depth, and the hat indicates the value of this quantity taken specifically at the luminosity shell. Once the position of the luminosity shell is found, the observed bolometric luminosity can be estimated as $\hat{E}/(t-{t}_{0})$, where $\hat{E}$ is the internal energy of the luminosity shell.

Given the relatively shallow outer envelope density profiles we find in realistic stellar models, we explore where the luminosity depth is at each time. In integral form, the diffusion time at the luminosity shell is (for the details of the derivation, see the Appendix)

$$\begin{eqnarray}&&{\hat{t}}_{{\rm{diff}}}={\int }_{\hat{r}}^{R}\displaystyle \frac{3\tau {dr}}{c},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&\tau (r)={\int }_{r}^{R}\kappa \rho {dr}.\end{eqnarray} \tag{ 2 }$$

The top panel of Figure 3 shows the position (in mass coordinates) of the luminosity shell found in our calculations for the MESA model with ${M}_{{\rm{ZAMS}}}=14\,{M}_{\odot }$ and the final energy ${E}_{{\rm{fin}}}={10}^{51}\,{\rm{erg}}$. The radius of this model is $R=848\,{R}_{\odot }$ and the ejecta mass is ${M}_{{\rm{ej}}}=10.93\,{M}_{\odot }$. The plotted quantity ${\hat{m}}_{{\rm{sh}}}$ is the difference between the total mass of the model and the mass coordinate of the luminosity shell as a function of time. The green curve finds ${\hat{m}}_{{\rm{sh}}}$ using our integral definition for ${\hat{t}}_{{\rm{diff}}}$ given in Equation (1). From this, one can see where photons are diffusing from at a given time. The blue curve finds ${\hat{m}}_{{\rm{sh}}}$ using the condition $\hat{\tau }=c/\hat{v}$ for defining the luminosity shell, where v is the velocity of matter. This alternative relation (in comparison with the integral we use in Equation (1)) arises from simplifying the integral definition for the diffusion time by taking

$$\begin{eqnarray}&&{t}_{{\rm{diff}}}\approx \displaystyle \frac{\hat{\tau }\hat{d}}{c},\end{eqnarray} \tag{ 3 }$$

where $\hat{d}$ is the width of the luminosity shell, and then setting $\hat{d}=\hat{v}(t-{t}_{0})$. This simpler description is often used in analytic and semi-analytic works (NS10; Piro 2012). In general, the integrated diffusion depth is larger at any given time. This is because the relatively shallow density profile of realistic models makes the full integral crucial for deriving the optical depth. In contrast, if the density profile was steeper, only the conditions at the luminosity depth would really be important and the integral would not be as important.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 3 compares different ways for calculating the bolometric luminosity of the same model. The magenta solid curve shows the bolometric luminosity at the photosphere ($\tau =2/3$), as returned by SNEC. Blue and green dashed curves show the light curves computed as $\hat{E}/(t-{t}_{0})$ using the two different conditions discussed for determining the location of the luminosity shell. To compute $\hat{E}$, we use the SNEC output for the internal energy at each grid point and find the total amount of this energy between the luminosity shell and the photosphere. Figure 3 shows that the bolometric luminosity at the photosphere agrees well with the luminosity calculated as $\hat{E}/(t-{t}_{0})$ until $\sim 10\mbox{--}20$ days after shock breakout, provided the condition ${\hat{t}}_{{\rm{diff}}}=t-{t}_{0}$ using Equation (1) is employed to find the location of the luminosity shell. At the same time, the luminosity $\hat{E}/(t-{t}_{0})$ calculated from the condition $\hat{\tau }=c/\hat{v}$ considerably underestimates the photospheric luminosity and has a different slope. This difference is because at larger luminosity depth (see the top panel of Figure 3) there is more energy available and thus the shock cooling is more luminous when the depth is calculated correctly.

We have already seen that at early times the light curves are determined by a polytrope with $n\ne 1.5$, but what about at later times when n = 1.5? Whether or not this is ever satisfied in the shock-cooling phase depends on when recombination starts being important. Figure 4 shows the light curves of the models from Figure 2 for final energy ${E}_{{\rm{fin}}}={10}^{51}\,{\rm{erg}}$. Colored circles indicate the time ${t}_{{\rm{conv}}}$ when the luminosity shell of each model computed as ${\hat{t}}_{{\rm{diff}}}=t-{t}_{0}$ enters the convective part of the model's envelope. Colored triangles indicate times when the effective temperature of the radiation goes down to $7500\,{\rm{K}}$, which we take as a rough criterion for the onset of recombination. From Figure 4 it is clear that for most models the time it takes for the luminosity shell to reach the convective part of the envelope is comparable to the time when recombination sets in. We therefore conclude that there is rarely much time during which the emission from shock cooling would be consistent with coming from a stellar structure described by an n = 1.5 polytrope. Comparing the two sets of models, one can see that the time it takes for the luminosity shell to reach the convective parts of the envelopes is shorter for the MESA models than for the KEPLER models. As was shown in Figure 2, the non-convective parts of the envelopes in the MESA models are on average smaller in size and have lower density than the non-convective parts of the envelopes in the corresponding KEPLER models (both factors reducing the diffusion time). As a consequence, for a few of the MESA models, the n = 1.5 parts of their envelopes control their cooling emission for $\sim 5\mbox{--}10$ days.

Zoom In Zoom Out Reset image size

During the last several decades, numerical modeling of SN IIP bolometric light curves was focused on reproducing their gross properties, such as the length and the luminosity of the plateau (see Popov 1993; Kasen & Woosley 2009; Bersten et al. 2011). However, as was shown in the analytical works of NS10 and Goldfriend et al. (2014), the slopes of the light curves during the shock-cooling and plateau phases encode important information on the structure of the density profiles of the progenitor stars. Increasing abundance and quality of observational data has made it possible to conduct systematic studies of the early slopes for large sets of SN IIP light curves (Anderson et al. 2014; González-Gaitán et al. 2015; Garnavich et al. 2016; Valenti et al. 2016). Modeling the slopes of the observed light curves and deducing the characteristics of the progenitor stars based on these slopes will be a very important task for future research.

In the current work, we instead focus on simply modeling the rise times of SNe IIP. Being more robust and easier to measure than the power-law exponent of the bolometric light curve, the rise time can still give us important information on the progenitor characteristics. Our work is also motivated by a number of recent works considering large sets of early SN IIP light curves and focusing on their rise times (e.g., Gall et al. 2015; González-Gaitán et al. 2015; Rubin et al. 2016).

Figure 5 shows g-band light curves for a representative set of models from Figure 2 for a final energy ${E}_{{\rm{fin}}}={10}^{51}\,{\rm{erg}}$. In our calculations, the SN emits as a blackbody at all times. Although the Rayleigh–Jeans part of the spectrum should not be strongly affected by this assumption during the shock-cooling phase (Tominaga et al. 2011), this will be more directly addressed in future work (T. Shussman et al. 2016, in preparation). The light curves generated by SNEC are plotted with solid lines. The colored circles indicate the local maxima of the light curves. Here we define the rise time, ${t}_{{\rm{rise}}}$, as the time between shock breakout and the local maximum of the color light curve in a given band. The error bars show the intervals of time during which the magnitude differs from the maximum by less than $0.01\,{\rm{mag}}$. Later in this paper (Figures 6 and 7) we take this criterion as a definition of the uncertainty with which the rise time can be measured. We choose the g band for the current study because in the u band the light curves are more sensitive to non-thermal effects like iron-group line blanketing (see Figure 8 of Kasen & Woosley 2009), which are not properly taken into account in SNEC. Furthermore, the g band is a good match to many current and future optical surveys. In the r band, the light curves become very flat, which makes it difficult to robustly identify the maximum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that in some of our calculations we find that the rise time depends on the duration of the thermal bomb. To investigate this dependence, we modeled the light curves of the KEPLER models for bomb durations of 0.1, ${10}^{-2}$, ${10}^{-3}$, and $5\times {10}^{-4}\,{\rm{s}}$. The vast majority of models show rise times independent of the bomb duration. The largest variations in ${t}_{{\rm{rise}}}$ ($\sim 4\,{\rm{days}}$) are seen for models with ZAMS masses in the range between $\sim 13\,{M}_{\odot }$ and $\sim 16\,{M}_{\odot }$ and at ${E}_{{\rm{fin}}}=2\times {10}^{51}\,{\rm{erg}}$. The rise times for these model converge (do not change with further decrease of the bomb duration) only at ${t}_{\mathrm{bomb}}\leqslant {10}^{-3}\,{\rm{s}}$. We believe that the observed dependence is due to the way the outer core structure in these models reacts to abrupt versus to more gradual energy input. In first principles core-collapse supernova simulations (e.g., Bruenn et al. 2016), the initial energy input by the expanding shock is abrupt, but most of the explosion energy builds up only gradually over ∼$1\,{\rm{s}}$. Since our thermal bomb approach does not allow us to fully realistically model the energy injection, we chose a thermal bomb duration of ${10}^{-3}\,{\rm{s}}$, which gives converged rise times for all our models.

For comparison with SNEC results, the dashed lines in Figure 5 show the light curves obtained from Equations (29) and (31) of NS10 for RSGs. For the values of mass and energy in these equations we use the ejecta mass ${M}_{{\rm{ej}}}$ and the final energy ${E}_{{\rm{fin}}}$ of each model. We choose the radius in these equations, ${R}_{{\rm{adj}}}$, in such a way that the maximum of the analytical light curve coincides in time with the maximum of the corresponding numerical SNEC light curve (the colored triangles indicate the maxima of the analytical light curves). By doing so, we mimic the way in which the progenitor radii were derived from the rise times in Section 5.6 of González-Gaitán et al. (2015). One can see that the values of ${R}_{{\rm{adj}}}$ are on average a factor of $\sim 1.5\mbox{--}2.5$ smaller than the actual radii of the RSG models. This result shows the risk of using existing analytical equations for estimating the radii of RSGs based on the value of the rise time only. Note that the same conclusion was demonstrated in the work of Rubin et al. (2016).

As an aside, we note that from Figure 5 it may seem that the numerical and analytical light curves have different breakout times (note the offset of the light curves at early times). This is explained by the fact that during the planar phase of post-explosion expansion (first one to two days), the g-band flux of the NS10 model goes down and starts to rise once the spherical phase of the expansion begins (this is also seen in the right hand panel of Figure 1 of González-Gaitán et al. 2015).

Figure 5 also gives an idea about the correlation between the maximum brightness and the rise time of our light curves. There is no definite answer in the literature on whether or not this correlation is present in the observational data. Valenti et al. (2014) suggested that slowly rising SNe IIP have higher maximum luminosities (see also Gal-Yam et al. 2011). The same trend was found in Gall et al. (2015). At the same time, Faran et al. (2014b) and Rubin et al. (2016) do not find this correlation in their samples, and only a weak correlation is found in the sample of Valenti et al. (2016). Closely related to this question is the possible correlation between the rise time and the post-maximum decline of the light curve, which could help clarifying the nature of the Type IIP and Type IIL SNe. From Figure 5, both MESA and KEPLER models suggest a positive dependence between the rise time and the maximum brightness of the light curves in the g band for a given value of final energy. At the same time, we do not consider the post-maximum slopes and therefore cannot distinguish between the IIP-like and IIL-like light curves in our models.

We next run a full collection of explosions using all stellar models and final energies of $0.5\times {10}^{51}$, ${10}^{51}$, and $2\times {10}^{51}\,{\rm{erg}}$. The main results of this set of calculations are summarized in Figures 6 and 7, where we plot g-band ${t}_{{\rm{rise}}}$ for all considered models as a function of R and ${M}_{{\rm{ej}}}$, respectively. Figure 6 demonstrates that there is a strong correlation between the rise time and the radius of the progenitor. Furthermore, both the MESA and KEPLER models show similar dependencies on R. On average, larger final energies have longer rise times, but this dependence is less strong. On the other hand, when we compare the rise time to the ejecta mass in Figure 7, the results are more mixed. The KEPLER models (bottom panel) show some correlation, but the spread is much larger than in the radius comparison. In addition, some of this dependence clearly comes from the strong correlation between mass and radius in the KEPLER models (Figure 1). The MESA models (top panel) show even less correlation with different masses sometimes having the same rise time.

Zoom In Zoom Out Reset image size

In principle, one should be able to derive the dependence of the rise time on the radius and ejecta mass of the model in the form $\mathrm{log}{t}_{{\rm{rise}}}={a}_{1}\mathrm{log}{M}_{{\rm{ej}}}+{a}_{2}\mathrm{log}R+{a}_{3}$ with constant coefficients a₁, a₂, and a₃. Unfortunately, correlations between radius and mass as shown in Figure 1 make it difficult to isolate these dependencies. As a consequence, the values of ${t}_{{\rm{rise}}}$ for our models form nearly a line in the 3D space $(\mathrm{log}{M}_{{\rm{ej}}},\mathrm{log}R,\mathrm{log}{t}_{{\rm{rise}}})$, and we find that it is not conclusive to fit a surface of the form $\mathrm{log}{t}_{{\rm{rise}}}={a}_{1}\mathrm{log}{M}_{{\rm{ej}}}+{a}_{2}\mathrm{log}R+{a}_{3}$. In the future, this dependence could possibly be inferred by fitting a much larger range of models with different masses and radii.

On the other hand, comparing the MESA models in the top panels of Figures 6 and 7 reveals that there is a much clearer mapping between the radius and rise time than between ejecta mass and rise time. This is especially apparent in the mass dependence of the rise time in the MESA models (top panel of Figure 7), which shows that the same ejecta masses can have different rise times when the radius is different. This was generally expected from previous analytic work (NS10; Piro & Nakar 2013), but it is helpful how clearly it is seen here. Furthermore, Figure 6 shows that both the MESA and KEPLER models have similar relationships (both slope and normalization) between radius and rise time. Motivated by this, we find a linear fit for $\mathrm{log}{t}_{{\rm{rise}}}$ as a function of $\mathrm{log}R$. Colored lines in Figure 6 show the fits taken separately for each set of models and each value of the final energy. The black dashed line is the common fit in the g band for all models and all final (explosion) energies, which has the following numerical coefficients,

$$\begin{eqnarray}&&\mathrm{log}R\,[{R}_{\odot }]\,\\ &&\qquad =(1.225\pm 0.178)\mathrm{log}{t}_{{\rm{rise}}}\,[{\rm{day}}]+1.692\pm 0.490.\end{eqnarray} \tag{ 4 }$$

The errors quoted are standard errors computed from the error bars given in Figure 6. From Figure 6, one can see that the slope of the radius dependence is fairly insensitive to the model type (KEPLER or MESA) and final energy. Equation (4) can be used as a tool to infer the radii of SN IIP progenitors from current and future transient surveys. This will be especially useful for analyzing large collections of events where it is impractical to do individual modeling.

An analogous relation to Equation (4) between the rise time in the optical part of the spectrum and the progenitor characteristics was obtained in Gall et al. (2015). From the analytical model of Arnett (1980, 1982), they found that the optical rise timescales as ${t}_{\mathrm{opt} \mbox{-} \mathrm{rise}}\propto {R}^{1/2}{T}_{{\rm{peak}}}^{-2}$, while for the Rabinak & Waxman (2011) model, ${t}_{\mathrm{opt} \mbox{-} \mathrm{rise}}\propto {R}^{0.55}{T}_{{\rm{peak}}}^{-2.2}{E}_{51}^{0.06}{M}_{{\rm{ej}}}^{-0.12}$, where ${T}_{{\rm{peak}}}$ is the effective temperature at the peak and E₅₁ is the final explosion energy in units of ${10}^{51}\,{\rm{erg}}$. In comparison, we find a somewhat stronger relation of ${t}_{{\rm{rise}}}\propto {R}^{0.816}$, although we do not include the temperature dependence. In the future, by running a large set of models with more diversity in their ejecta-mass–radius relations, we will be able to better study how ${t}_{{\rm{rise}}}$ depends on other factors besides radius.

Assuming that the rise times we find are representative of what is found in nature, we can also ask what should be expected for the rise times in a larger sample of events. Since the distribution of massive stars is understood in terms of an initial mass function (and not a radius function), we have to assume a mass–radius relation to explore this. Here, we choose to exploit the mass–radius relation of the KEPLER models because it is stronger than the mass–radius relation of the MESA models⁸ (see Figure 1). For the purpose of this calculation, we use a linear fit of the rise times as a function of ZAMS mass for KEPLER models at final energy ${E}_{{\rm{fin}}}={10}^{51}\,{\rm{erg}}$. Using Monte Carlo techniques, we generate a large ($\gt {10}^{5}$ mass values) sample of massive stars with ZAMS masses in the range between $8\,{M}_{\odot }$ and $20\,{M}_{\odot }$ (motivated by the masses approximately inferred by pre-explosion imaging of SNe IIP, Smartt et al. 2009) with a Salpeter initial mass function of

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{{dM}}_{{\rm{ZAMS}}}}\propto {M}_{{\rm{ZAMS}}}^{-2.35},\end{eqnarray} \tag{ 5 }$$

and calculate the rise times for each value of ${M}_{{\rm{ZAMS}}}$ using the derived fit. The median value of the resulting set of rise times is 7.27 days with a median absolute deviation of 1.3 days, which is similar to the observed median value of 7.5 ± 0.3 days found by González-Gaitán et al. (2015). This implies a median radius of $559\,{R}_{\odot }$ for SNe IIP progenitors. We want to emphasize again that these values were derived assuming the mass–radius relation of the KEPLER models, while Equation (4) does not imply any mass–radius relation and was obtained from both sets of models.

It is interesting to compare our results to the observational samples of massive stars (e.g., Levesque et al. 2005, 2006; Arroyo-Torres et al. 2013). For example, the observed galactic RSGs have radii in the range $100\mbox{--}1520\,{R}_{\odot }$ (Levesque et al. 2005), while Magellanic Cloud RSGs span radii of $470\mbox{--}1310\,{R}_{\odot }$ (Levesque et al. 2006). For the smallest observed radius of $100\,{R}_{\odot }$, Equation (4) predicts ${t}_{{\rm{rise}}}=1.78\,{\rm{days}}$, while for the largest observed radius of $1520\,{R}_{\odot }$ it gives ${t}_{{\rm{rise}}}=16.45\,{\rm{days}}$. From Figures 6 and 7 of González-Gaitán et al. (2015), one can see that the bulk of the observed rise times in the g band lies between these two values with only few outliers. Therefore, we conclude that Equation (4) predicts the correct range of the rise times for the observed RSG radii. At the same time, our study does not let us make conclusions concerning the advantage of one set of RSG models with respect to another. For the considered ZAMS masses and MESA parameters, both the MESA and KEPLER models have radii within the observed range. The most important result for the current work is that the two sets of models yield a consistent dependence between the radii and the rise times of their light curves, as given by Equation (4).

It should be mentioned that González-Gaitán et al. (2015) also present numerical simulations of the light curves from RSGs and obtain results different from ours. In particular, they use explosion models by Tominaga et al. (2009, 2011) generated with the radiation hydrodynamics code STELLA (Blinnikov & Bartunov 1993; Blinnikov et al. 1998). They show that the light curves of RSG progenitor models with radii in the range $500\mbox{--}1400\,{R}_{\odot }$ have rise times too long compared to observations. Instead, they found better agreement with observational data when using models with a modified H/He envelope structure and radii in the range $80\mbox{--}400\,{R}_{\odot }$. The difference between our results and the results of González-Gaitán et al. (2015) may be partially explained by the fact that we consider lower ZAMS masses of $10\mbox{--}20\,{M}_{\odot }$, while their models span the range $13\mbox{--}30\,{M}_{\odot }$. Another important factor is the difference between the SNEC and STELLA codes. Although STELLA uses a multi-frequency approach to radiation transport, SNEC, uses very complete OPAL opacity tables. In order to fully understand the differences between SNEC and STELLA results, a detailed comparison study will be required. We defer such a study to future work.

Using SNEC (Morozova et al. 2015), we have exploded a set of massive star models from both the KEPLER and MESA stellar evolution codes to study how the rise of the SN light curve depends on progenitor and explosion characteristics. We find that the strongest correlation is between the g-band rise time and the radius at the time of explosion, and provide a formula that relates these two properties in Equation (4). This can be used in future analyses of SNe IIP observations.

To better understand what properties of the progenitor are controlling the early light curve, we examined the envelopes of RSGs obtained with both stellar evolution codes. We find that their convective envelopes do not extend all the way to the surface. In fact, all the considered models have regions close to their surface where the power-law exponent n is smaller than 1.5 and its value varies between different ZAMS masses. These regions are important for the early light curve, since the luminosity shell passes through them during the first 10–25 days after shock breakout. Due to the shallow density profiles in these regions, the simple estimate $\hat{\tau }=c/\hat{v}$ for the optical depth at the luminosity shell is inadequate. This explains the differences we find with previous semi-analytic treatments of the early light curve.

The results obtained in this paper add to a long-standing discussion about the progenitor radii of SNe IIP. Based on the results of their numerical models, Dessart et al. (2013) showed that the color properties of SNe IIP may be explained by small progenitor models ($\sim 500\,{R}_{\odot }$), while larger progenitors would produce light curves that remain too blue for too long. González-Gaitán et al. (2015) came to a similar conclusion based on a comparison of the observed rise times to analytical models of the shock-cooling phase and to hydrodynamical models. The radii they obtain are a factor of ∼2 smaller than the observed radii of RSGs (see Levesque et al. 2005, 2006), and have the median value $\sim 350\,{R}_{\odot }$. Progenitor radii $\lesssim 500\,{R}_{\odot }$ are also suggested by the results of Shussman et al. (2016) and Garnavich et al. (2016). At the same time, the works of Valenti et al. (2014) and Bose et al. (2015) for two SNe IIP deduce large radii (up to $\sim 800\,{R}_{\odot }$) for their progenitor stars. Rubin et al. (2016) find that the progenitor radii are weakly constrained by comparison to analytical shock-cooling models (at least based on R-band photometry alone). In this paper, we have shown that deriving the progenitor radius by fitting the rise times of the observed light curves with analytical shock-cooling models may underestimate the radii. On the other hand, exploding the MESA and KEPLER stellar evolution models without reducing their radii, we find a reasonable agreement with the median value of the observed g-band rise times from González-Gaitán et al. (2015).

One of the main strengths of our technique is that it does not require that the stellar models we use exactly replicate the progenitors that exist in nature. Instead, having a more diverse sample of models with different mass–radius relations gives us an increasingly better handle on how the early rise depends on the progenitor properties. Therefore, for future work it will be useful to explode an even larger set of stellar models. This would help better test the relation we find between rise time and radius, but also help us to better understand how sensitive the rise time is to the ejecta mass. In addition, one could explore whether ${T}_{{\rm{peak}}}$, the temperature at peak, is useful for tightening this relation, as has been found in previous semi-analytic work (see Gall et al. 2015). In this way, a fairly easy additional observable, namely the colors at peak, could be used to make these radius measurements more robust. This will be useful as a tool for future transient surveys, as well as for comparison with pre-explosion imaging of SNe and studies of massive stars that hope to connect these progenitors to the SNe they will eventually make.

We would like to thank our referee, S. González-Gaitán, for valuable comments,which helped to improve the quality of our paper. We acknowledge helpful discussions with and feedback from D. Clausen, L. Dessart, R.J. Foley, S.R. Kulkarni, K. Maeda, O. Pejcha, R. Sari, D. Radice, B.J. Shappee, and T. Sukhbold. We thank T. Sukhbold for providing us with the pre-collapse stellar evolution models from KEPLER. We thank Ehud Nakar for helpful feedback on our paper and also for sharing a nearly complete draft of his own work (Shussman et al. 2016). SNEC and the light curves presented here are available from https://stellarcollapse.org/Morozova2016. This work is supported in part by the National Science Foundation under award Nos. AST-1205732 and AST-1212170, by Caltech, by the Sherman Fairchild Foundation, and by the International Research Unit of Advanced Future Studies, Kyoto University. The computations were performed on the Caltech compute cluster Zwicky (NSF MRI-R2 award No. PHY-0960291), on the NSF XSEDE network under allocation TG-PHY100033, and on NSF/NCSA Blue Waters under NSF PRAC award no. ACI-1440083. This paper has been assigned Yukawa Institute for Theoretical Physics report number YITP-16-32.

We estimate the diffusion timescale for a static, radiation-pressure-dominated medium described by the energy balance equation

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial U}{\partial t}=-\displaystyle \frac{\partial L}{\partial m},\end{eqnarray} \tag{ 6 }$$

where the work by expansion is omitted and $\rho =\rho (r,{t}_{0})$ is the density of matter at a given time t₀. Here, $U={{aT}}^{4}$ is the energy density of the radiation at temperature T, a is the radiation constant, m is the mass coordinate related to the radial distance as

$$\begin{eqnarray}&&\displaystyle \frac{\partial r}{\partial m}=\displaystyle \frac{1}{4\pi {r}^{2}\rho }.\end{eqnarray} \tag{ 7 }$$

The radiative luminosity L in the diffusion approximation is

$$\begin{eqnarray}&&L=-{(4\pi {r}^{2})}^{2}\displaystyle \frac{{ac}}{3\kappa }\displaystyle \frac{\partial {T}^{4}}{\partial m},\end{eqnarray} \tag{ 8 }$$

where κ is the Rosseland mean opacity and c is the speed of light.

Using Equations (7)–(8), Equation (6) can be rewritten in the form of a diffusion equation for the radiation energy density

$$\begin{eqnarray}&&\displaystyle \frac{\partial U}{\partial t}=\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{c}{3\kappa \rho }\displaystyle \frac{\partial U}{\partial r}\right).\end{eqnarray} \tag{ 9 }$$

From the form of Equation (9), the time ${\rm{\Delta }}t$ it takes radiation to diffuse through a distance ${\rm{\Delta }}r$ can be estimated as

$$\begin{eqnarray}&&{\rm{\Delta }}t\sim \displaystyle \frac{3\kappa \rho }{c}{({\rm{\Delta }}r)}^{2}.\end{eqnarray} \tag{ 10 }$$

Therefore, the integral form for the characteristic diffusion time from the radial coordinate r to the surface r = R is

$$\begin{eqnarray}&&{t}_{{\rm{diff}}}(r,{t}_{0})={\int }_{r}^{R}\displaystyle \frac{3\tau {dr}}{c},\end{eqnarray} \tag{ 11 }$$

where $\tau (r,{t}_{0})={\int }_{r}^{R}\kappa \rho (r,{t}_{0}){dr}$ is the optical depth.