HYDRODYNAMICAL MODELS OF TYPE II PLATEAU SUPERNOVAE

Our code follows a one-dimensional, Lagrangian prescription that solves for the variables: radius (r), velocity (u), density (ρ), and temperature (T) as a function of the Lagrangian mass coordinate m, using the following equations which simply express conservation laws:

$$\begin{equation} \frac{\partial r}{\partial t}= u, \end{equation} \tag{ 1a }$$

$$\begin{equation} V= \frac{1}{\rho }= \frac{4 \pi }{3} \frac{\partial r^3}{\partial m}, \end{equation} \tag{ 1b }$$

$$\begin{equation} \frac{\partial u}{\partial t}= - 4 \pi r^2 \,\frac{\partial }{\partial m}(P+q) - \frac{G m}{ r^2}, \end{equation} \tag{ 1c }$$

$$\begin{equation} \frac{\partial E}{\partial t}= \epsilon _{\rm Ni} - \frac{\partial L}{\partial m} - (P+q) \, \frac{\partial V}{\partial t}, \end{equation} \tag{ 1d }$$

$$\begin{equation} L = -(4 \pi r^2) ^2 \, \frac{\lambda ac}{3 \kappa } \, \frac{ \partial T^4}{\partial m}. \end{equation} \tag{ 1e }$$

Here, V is the specific volume, P is the total pressure (of gas and radiation), and q is the artificial viscosity, which is included in the equations to spread the pressure and energy over several mass zones at the shock front. There are many expressions for the artificial viscosity, all dependent on the velocity gradient, which aim at providing a convenient interpolation scheme between unshocked and shocked fluid. We adopt the expression given by Von Neumann & Richtmyer (1950). E is the internal energy per unit of mass (of gas and radiation) and _Ni is the energy deposited by the radioactive decay of nickel as we describe in Section 2.3. We do not consider other sources of cooling or heating in Equation (1d), such as losses due to neutrino processes or energy released by thermonuclear reactions. Even if neutrinos are very important in the formation of the SW, as the explosion depends noticeably on the efficiency of their energy transfer, most of them are emitted before the SW reaches the stellar photosphere, so they have no effect on later epochs of the SN evolution (see Hillebrandt 1994; Burrows 1991; Janka et al. 2007). The energy released by explosive nucleosynthesis is much less than the energy of the SW, as has been previously established (see Imshenik & Nadezhin 1965; Woosley & Weaver 1994; Arnett 1996). T is the temperature of both matter and radiation, κ is the Rosseland mean opacity, and L is the luminosity. Finally, λ is the so-called flux-limiter, included in the equation of radiative transfer in the diffusion approximation to ensure a smooth transition between diffusion and free-streaming regimes to assure causality. The expression adopted for λ is

$$\begin{equation} \lambda = \frac{ 6 + 3 R}{6 + 3 R + R^2}, \end{equation} \tag{ 2 }$$

where

$$\begin{equation} R= \frac{ \mid \nabla T^4 \mid }{\kappa \rho \ T^4}= \frac{4 \pi r^2}{\kappa T^4} \bigg| \frac{ \partial T^4}{\partial m} \bigg|. \end{equation} \tag{ 3 }$$

The quantities E, P, q, and κ are functions of ρ, T, and chemical composition. Further details on these relations are given in Section 2.2.

As boundary conditions we use u = 0 near the center (at m = M_core, where typically we adopt M_core = 1.4 M_☉), and P_gas = ρ = 0 at the surface (m = M).

We discretize the previous equations using a space-centering discretization with the extensive quantities evaluated at the interfaces and the intensive quantities in the midpoints of the grid zones. Two time steps are adopted in each cycle: one to advance the velocity, and the other to advance the material state variables. The time step is chosen depending on limitations of stability and accuracy (see the discussion below). Special care is taken in the centering of the opacity used in the discretization of Equation (1e) in order to prevent numerical noise from appearing due to the propagation of the radiation flow at the steep front where the opacity changes significantly (Christy 1967).

We use an explicit scheme for the integration of the hydrodynamic equations, but a semi-implicit scheme for the temperature, similar to the one used by Falk & Arnett (1977). The equations are linearized in δT and solved iteratively for each time step using the tridiagonal method. The discretization typically uses 300 mesh points, with a finer sampling for the outer layers typically smaller than 10⁻⁶ M_☉. This value was chosen based on tests performed with our model and following previous works (Woosley 1988; Ensman & Burrows 1992) that show that the early LC is sensitive to the mass zoning in the outer layers when a coarse grid is used.

Given that we are using an explicit hydrodynamic scheme, the time step should be chosen as a fraction of the minimum of the Courant condition for all zones in order to achieve stability. We note that we have also imposed additional conditions on the time step: we have required that changes in temperature, density, and flux over one time step be less than 5%.

The formation of the SW following core collapse is simulated by artificially adding internal energy ("thermal bomb") almost instantaneously in the central region of the core. We usually set a mass cut of 1.4 M_☉ which is the material assumed to collapse and form a neutron star or black hole. Specifically, we employ an exponential function both in mass and in time to distribute the injected energy across several layers and some time steps, which helps to improve the numerical treatment. With the scale factors used in this work for the exponential functions almost all of the explosion energy was completely injected in a shell of 0.1 M_☉ and within less than 1 minute, i.e., a short time compared with the hydrodynamical timescale. Although there may be some differences in the velocity and temperature profiles during the shock propagation on the scale factors used, our values were chosen such that the differences in the velocity profiles using this value or lower values were not significant. We have also tested explosions generated by injecting kinetic energy and we have obtained similar results. The latter method, however, leads to very short numerical time steps, and therefore to slower calculations. We have checked the accuracy of our calculation method by testing the conservation of energy. The total amount of energy is conserved within 0.6% but if we consider the conservation of energy between two consecutive time steps, this is even better, within 4 × 10⁻⁶. We consider it very acceptable for our purposes.

The equation of state (EOS) is calculated using simple expressions for T, ρ, composition, and for the ionization degrees of hydrogen and helium corresponding to local thermodynamical equilibrium (LTE). Therefore, the degree of ionization is determined by solving the corresponding set of Saha equations for ionization of hydrogen and the first and second ionization of helium. Ionization of heavy elements is neglected in our EOS (but, of course, it is taken into account in the calculation of the opacity). The degeneracy pressure is also included in the EOS. We tested our results using a more sophisticated EOS, such as the one available from the Los Alamos Tables (Rogers et al. 1996), and we did not find any significant differences with respect to the results obtained using our simple EOS. This is expected because of the low densities attained in most of the layers of the models.

The Rosseland mean opacities, κ, used in our calculations are derived using the OPAL opacity tables (Iglesias & Rogers 1996, and references therein). As these tables are given for T > 6 × 10³ K, we complement them with the opacity table provided by Alexander & Ferguson (1994) for lower temperatures, which includes molecular opacities. The tables are interpolated to each other in order to guarantee a smooth transition at T = 10⁴ K.

These tables allow us to calculate opacities for several metallicities. Also, for a fixed metallicity, different mixtures of H, He, C, and O can be used. Although in this work we adopt a fixed value of Z = 0.02, departures from this value (as expected in the inner regions of the object) are taken into account as excesses of C and O with respect to the values of the adopted metallicity, at the expense of He.

Figure 1 shows the opacity given by the tables as a function of temperature for Z = 0.02 and different densities. The range of temperature shown corresponds to the values reached during the evolution of an SN II-P while the values of density are restricted to those achieved during the plateau phase. Also shown in the plot is the contribution to the opacity from electron scattering. Note that electron scattering is the dominant source of opacity for T > 10⁴ K and ρ < 10⁻¹⁰ g cm⁻³.

Zoom In Zoom Out Reset image size

The Rosseland mean opacity includes scattering and absorption processes. Scattering is the dominant process in the SN ejecta during the plateau phase when the densities are ρ < 10⁻¹⁰ (see Figure 1), and this will be so until most of the electrons recombine with ions and absorption processes become an important source of opacity. On the other hand, in rapidly expanding envelopes where large velocity gradients are present, the Rosseland mean opacity underestimates the true line opacity (Karp et al. 1977), which hinders the estimation of the actual opacity in the outermost (recombined) layers. Another effect that is not included in the calculation of κ is the non-thermal excitation or ionization of electrons that are created by Compton scattering of γ-rays emitted by radioactive decay of ⁵⁶Ni and ⁵⁶Co. The LTE ionization used in the calculation of κ considerably underestimates the true ionization. The correct way to treat these effects is to calculate the actual contribution of non-thermal ionization to the opacity and to include the expansion opacity of lines. However, such treatment is beyond the scope of this paper. We adopt an alternative approach to the problem that has been extensively used in the literature which consists in using a minimum value of the opacity (or "opacity floor") to partially solve the shortcoming in the Rosseland mean opacity. Given the dependence of the opacity on composition, usually two values of the opacity floor are used: one for the H-rich envelope material and another for the metal-rich core. The opacity at each time and mesh point is chosen as the maximum value between the tabulated Rosseland mean opacity and the opacity floor for the corresponding composition.

Values of the opacity floor should be based on contributions to the opacity that are not included in the Rosseland mean opacity tables. There are differences in the values adopted in the literature. For example, Herzig et al. (1990) adopted a value of 0.01 cm² g⁻¹ for the whole structure. Shigeyama & Nomoto (1990) used minimum values for the bound–free and bound–bound opacities of k_bf = 9 × 10⁻³ Y cm² g⁻¹ and k_bf = 1 × 10⁻² Z cm² g⁻¹, where Y and Z denote the mass concentration of helium and heavy elements. Note that it is needed to add scattering and free–free absorption effects to the latter in order to obtain the Rosseland mean opacity. Swartz et al. (1991) employed 0.05 cm² g⁻¹ for the envelope material, and 0.1 cm² g⁻¹ for the metal-rich core material. And Young (2004) used a value of 0.25 cm² g⁻¹ for the helium-rich core, and 0.01 cm² g⁻¹ for the hydrogen-rich envelope. We do not intend to be exhaustive but the examples above provide a summary of the values adopted in the field of SNe II.

In this work, the minimum opacity values adopted are 0.01 cm² g⁻¹ for the envelope material, and 0.24 cm² g⁻¹ for the metal-rich core material. These values were chosen based on a comparison performed between our code and the STELLA code (Blinnikov & Bartunov 1993; Blinnikov et al. 1998) using the same initial model provided by Umeda & Nomoto (2005). Note that STELLA is an implicit hydrodynamic code that incorporates multi-group radiative transfer, and which additionally uses different opacity tables and includes the effect of line opacities (Sorokina & Blinnikov 2002). In spite of our approximations, it is important to remark that we obtain an excellent overall agreement of the LC given by both codes (both in terms of the duration of the plateau and of the morphology of the bolometric LC) when the opacity floor is set the values given above.

The radioactive decay of ⁵⁶Ni →⁵⁶Co → ⁵⁶Fe constitutes a relevant source of energy that heats the envelope of the SN. Several studies of SNe II-P have usually ignored the contribution of such a component during the early SN evolution and have only considered it beyond the end of the plateau phase. Such an approach would be correct only if ⁵⁶Ni were very deeply concentrated in the ejecta, in which case the local deposition of gamma rays is a good approximation because gamma photons can hardly diffuse out from the regions where they are formed. Nevertheless, as was shown in studies of SN 1987A, there is no reason to assume this type of ⁵⁶Ni distribution (Shigeyama et al. 1988; Woosley et al. 1988; Arnett 1988; Blinnikov et al. 2000, among others) and therefore we need to calculate the diffusion of gamma rays from the location where they are emitted to the outer regions. To calculate this, we solve the gamma-ray transfer in the gray approximation for any distribution of ⁵⁶Ni, assuming that gamma rays interact with matter only through absorption. It has been shown by comparison with Monte Carlo simulations of Type Ia SNe that the complex scattering process between gamma rays and electrons can be satisfactorily approximated as an absorptive process (Sutherland & Wheeler 1984; Swartz et al. 1995). The value for the gamma-ray opacity that we adopt here is κ_γ = 0.06 y_e cm² g⁻¹, where y_e is the number of electrons per baryon.

The rate of energy per gram released by the Ni–Co–Fe decay is

$$\begin{eqnarray} \epsilon _{\rm rad} &=& 3.9 \times 10^{10} \, \exp (- t/\tau _{\rm Ni}) + \, 6.78 \times 10^{9} [\exp (- t/\tau _{\rm Co}) \nonumber \\ && - \exp (- t/\tau _{\rm Ni})] \;{\rm erg\;g}^{-1} \;{\rm s}^{-1}, \end{eqnarray} \tag{ 4 }$$

where τ_Ni = 8.8 days and τ_Co = 113.6 days are the mean lifetimes of the radioactive isotopes. The amount of energy deposited at each point is given by the solution of the gamma-ray transfer multiplied by the previous expression. In Figure 2, we show the gamma-ray deposition⁵ for the case of a polytrope with index n = 3, initial mass of 10 M_☉, and different initial radii, assuming a constant distribution of ⁵⁶Ni up to 3 M_☉. Note that the diffusion of gamma rays from the region where they form becomes more noticeable as the object becomes more extended and diluted.

Zoom In Zoom Out Reset image size

The possibility to use an arbitrary distribution of ⁵⁶Ni allows us to study different types of mixing and their effect on the LC. We have found, as shown in Section 3.1.4, that radioactivity becomes an important source of energy even during the plateau phase if we allow extensive ⁵⁶Ni mixing.

There are two different types of initial (or pre-SN) models: those coming from stellar evolution calculations (or "evolutionary" models), and those from non-evolutionary calculations where the initial density and chemical composition are parameterized in a convenient way. In this work, we use double polytropic models in hydrostatic equilibrium as non-evolutionary pre-SN models. Although a single polytrope may represent very well the envelope of the pre-SN object, the inner, dense part that is expected for this type of object is not well reproduced. One way to improve this situation is to consider two polytropes: one representing the inner dense core, and the other accounting for the outer extended envelope. As initial composition we use parametric chemical profiles with mixing between layers of different chemical composition (see Figure 4 in Section 3.1 for SN 1999em). Mixing is expected to occur during the explosion due to hydrodynamical instabilities which can carry the hydrogen very deep into the core and the ⁵⁶Ni out into the hydrogen envelope as previous studies of SN 1987A have shown (Dotani et al. 1987; Shigeyama et al. 1988; Woosley et al. 1988; Arnett 1988; Haas et al. 1990; Blinnikov et al. 2000, among others). Unfortunately, we cannot properly take into account this effect in our one-dimensional prescription. To ameliorate this we implicitly include this effect by imposing mixing in our initial chemical profile. It is important to note that even if mixing is expected only after the shock propagation, the effect on the dynamics of imposing it at the initial time is not noticeable until the shock breakout, when the recombination front recedes into the ejecta.

In order to calculate numerically a double polytropic model, it is necessary to solve the equations of hydrostatic equilibrium and mass conservation assuming a polytropic approximation where the pressure has a dependence on density of the form: P = Kρ^γ = Kρ^n/(n+1), where K is a constant and n is the polytropic index. Two different polytropic relations with different indices (n_i and n_o) and constants (K_i and K_o) are adopted to mimic the characteristic structure of a red supergiant. To calculate the composite polytrope, the previous equations are numerically integrated outward from the inner border (stellar core) to a pre-selected point (fitting point) using the internal polytropic relation. A second inward integration is done from the outer border (stellar surface) to the fitting point using the external polytropic relation. The problem is equivalent to a two-point boundary condition for the case where there are unknown free parameters at both ends of the domain. In our case, the free parameters are n_i, K_i, and K_o, while the external index is fixed to n_o = 3. Starting from initial guesses, the values of the free parameters are iteratively sought so that the solution joins smoothly at the fitting point. This process is performed using a shooting method. Once the values of n_i, K_i, and K_o are found, the density (or pressure) distribution and mass distribution can be calculated. The initial temperature profile is calculated in an iterative fashion using our EOS to ensure hydrostatic equilibrium. This prevents the formation of a spurious SW.

Using this method, a variety of initial models in hydrodynamical equilibrium can be calculated by changing the fitting point and the boundary conditions for a given mass, radius, and central density, which allows us to study how the internal structure of the initial model affects the LC.

Alternatively to the parametric initial model, we have tested our code using initial models from different stellar evolution calculations. We found several features in the resulting LC that are not present in the observations, as was previously noted in the literature (see Utrobin & Chugai 2008). Specifically, several bumps and wiggles appear during the plateau phase. Regardless of the differences among evolutionary models, all of them have the common characteristic of showing a sharp boundary between layers with different chemical composition and a steeper jump in density between the helium core and the hydrogen envelope than in our double polytropic models. These two characteristics are the ones responsible for the unobserved bumps in the LC. We note that a possible reason for the relatively poor results we obtained from evolutionary initial models may be the use of one-dimensional calculations, which cannot take into account effects that produce mixing, such as Rayleigh–Taylor (R-T) instabilities. This point deserves further scrutiny, although it is beyond the scope of this paper.

In summary, we decided to use double polytropic initial models based on the following facts. (1) We have generally found a better agreement with observations using our parameterized profiles, in accordance to previous studies (Utrobin 1993; Baklanov et al. 2005; Utrobin 2007, among others). (2) There is a limited set of pre-SN models from stellar evolutionary calculations available to us. (3) Parameterized initial models allow us to easily vary physical properties and study their effect on the LC, which is critical for our goal of studying a large sample of SNe. However, we emphasize that the choice of polytropic initial models makes the connection with the structure of a progenitor star difficult to assess.

Several approximations are made in the equations of radiation hydrodynamics. We assume that the fluid motion can be described by one-dimensional, radially symmetric flow. The explosion mechanism of CCSNe is not well known but it may be a very asymmetric process. However, for this particular subtype of SNe with very extended hydrogen envelopes, the asymmetries expected from the explosion mechanism itself appear to be smoothed. This is supported by recent spectropolarimetric studies (Leonard & Filippenko 2005).

We use the equilibrium diffusion approximation to describe the radiative transfer. This approximation assumes that radiation and matter are strongly coupled with a single characteristic temperature and a spectral energy distribution described by a blackbody function (BB hereafter). The approximation breaks down at shock breakout and at late phases when the ejecta is completely recombined and the object becomes transparent. Fortunately, during the plateau phase—which is our main interest—this approximation is adequate. Also note that non-LTE calculations of SN II spectra show that deviations from LTE have significant effects on the lines but not on the overall continuum (Baron et al. 1996; Dessart & Hillier 2008). At late phases, for SNe that experience little interaction with the interstellar medium, as is the case of SNe II-P, the bolometric luminosity can be simply approximated as the luminosity deposited by the radioactive decay of ⁵⁶Co, at least during the early part of the radioactive tail. This is supported by the fact that the luminosity declines obeying an exponential law with a very similar rate to that of the decay of ⁵⁶Co. There was no attempt to apply our model beyond ∼180 days.

During the transition between plateau and radioactive tail, the envelope is fully recombined and the notion of a photosphere loses meaning. Thus, this transition regime is poorly described with our radioactive transfer prescription and therefore a detailed description of this phase cannot be assessed here. As stated above, we have performed a comparison of our calculations with those of the STELLA code, obtaining an excellent agreement between the bolometric luminosities derived by both methods. This is very satisfactory considering that the STELLA code involves a more sophisticated treatment of the radiative transfer by solving these equations using a multi-group prescription and including the effect of line opacities (Blinnikov & Bartunov 1993; Sorokina & Blinnikov 2002).

One of the main motivations for this work is the unprecedented database of BVI LCs of ∼30 SNe which offers the opportunity to significantly improve our understanding of SNe II-P and their progenitors. Given that our code produces bolometric LCs, we need to compute bolometric luminosities for our data set. In Bersten & Hamuy (2009), we found a tight correlation between bolometric correction and colors, which allowed us to calculate bolometric luminosities from BVI colors with an uncertainty of only ∼0.05 dex.⁶

Another critical parameter to compare with observations is the photospheric velocity yielded by our models. Expansion velocity estimates from the minimum of several spectral lines (H_α, Fe ii λ5169, H_β, and H_γ) of our data set were given by Jones et al. (2009). Since each line forms at a different shell, it is not straightforward to compare the observed velocities with photospheric velocities. Jones et al. (2009) got around this problem and derived congruent calibrations between the velocity derived from the absorption minimum of such lines and the photospheric velocity using two independent atmosphere models (Eastman et al. 1996 and Dessart & Hillier 2005; E96 and D05, hereafter, respectively). As noted by Jones et al. (2009), H_β provides a very good proxy to the photosphere velocity as it is not highly saturated as H_α, and is present over most of the SN evolution. Therefore, in the comparisons with our models, we chose to use the calibration given by Jones et al. (2009) for this particular line. However, in spite of the satisfactory behavior of this calibration in the high-velocity regime (v ≳ 5000 km s⁻¹), for lower velocities the atmosphere models fail to reproduce the behavior shown by the data (see Figure 9 of Jones et al. 2009). Thus, we decided to include also in our comparison the velocities estimated from the Fe ii λ5169 line, which is present at later phases of the SN evolution when the velocities are below the limit imposed by the calibration of hydrogen lines.

For consistency with the atmosphere models, we define the photospheric position as the layer where the total continuum optical depth is τ = 2/3. Note that to calculate τ we do not include the opacity floor because such a minimum is only included to take into account line effects, and therefore is a bound–bound opacity that does not contribute to form the continuum. With this definition, the photosphere follows the recombination wave (RW) as we show in Section 3.1.3. Note that this behavior of the photosphere is due to the fact that electron scattering is the dominant source of opacity (see Figure 1). Had we included the opacity floor in the definition of τ, the τ = 2/3 surface would be pushed well above the recombination front.

We remark that with this definition, the photosphere is essentially the surface of last scattering. However, the surface where the continuum is actually formed is located in a deeper layer called the "thermalization depth" and this is due to the dominance of the electron scattering over the absorption processes (Sobolev 1980; Hoflich 1991; Montes & Wagoner 1995). With our simple prescription of radiative transfer we cannot accurately determine the location of the thermalization depth because there radiation decouples from matter. Note that the color temperature, determined by a blackbody fit to the broadband photometry, is nearly coincident with the temperature of the thermalization depth but is greater than the effective temperature defined by T_eff = L/4πσR²_ph, where R_ph is the photospheric radius. In Bersten & Hamuy (2009), we derived a calibration between T_eff and color based on the E96 and D05 models. Using this calibration we study the evolution of T_eff for our sample of SNe II-P and we compare it with T_eff derived from our code.

In this section, we use a pre-SN model with initial parameters consistent with the optimal hydrodynamical model of Utrobin (2007). Specifically, we adopt an initial mass of 19 M_☉, radius of 800 R_☉, and explosion energy of E = 1.25 foe (1 foe = 1 × 10⁵¹ erg). This energy was released as thermal energy near the core of the object in a very short timescale as compared with the hydrodynamic timescale of our model. We also assume a nickel mass of 0.056 M_☉ which was determined by the luminosity of the radioactive tail. This parameter is quite different from the nickel mass of 0.036 M_☉ used by Utrobin (2007), which was determined using the quasi-bolometric luminosity given by Elmhamdi et al. (2003). Their luminosity was based on the integration of only UBVRI photometry with a constant value of 0.19 dex added to take into account the infrared luminosity, while our bolometric luminosity is based on a quantitative study of the bolometric correction for SNe II-P (see Section 2.5). Note also that our value for the nickel mass is closer to the estimate of 0.06 M_☉ given by Baklanov et al. (2005). In the following analysis, we denote this model as m19r8e1.25ni56.

In our calculations, we remove the central 1.4 M_☉ which is assumed to form a neutron star. The initial density profile as a function of mass and radius for model m19r8e1.25ni56 is shown in Figure 3. Note that the initial structure is composed of a dense core and an extended envelope characteristic of a red supergiant. The chemical composition profiles are shown in Figure 4. The outer parts of the envelope, M > 9 M_☉, have a homogeneous composition with mass fractions of X = 0.735, Y = 0.251, and Z = 0.02. From there inward, hydrogen and helium are mixed in order to prevent a sharp boundary between the H-rich and the He-rich layers. Such sharp boundaries are characteristic of stellar evolution models but fail to reproduce the observations. Note that we allow H to mix very deep inside the core and ⁵⁶Ni is mixed out in the envelope until ∼15 M_☉. This type of mixing is presumably due to R-T instabilities that occur behind the shock front, as is obtained in multi-dimensional hydrodynamic calculations (Mueller et al. 1991; Kane et al. 2000) and supported by studies of SN 1987A (Shigeyama et al. 1988; Woosley et al. 1988; Arnett 1988; Blinnikov et al. 2000, among others). The presence of H in the core leads to a smooth transition between the plateau and the radioactive tail. The distribution of ⁵⁶Ni to external layers helps to reproduce a plateau as flat as that observed in SN 1999em (see Section 3.1.3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows a comparison between the bolometric LC (solid line) obtained with our code for model m19r8e1.25ni56, and the observations of SN 1999em (dots). The luminosity due to ⁵⁶Ni → ⁵⁶Co → ⁵⁶Fe (dashed line) is also shown. Note the very good agreement between model and observations. The largest differences appear during the earliest phase and the transition to the radioactive tail. At the earliest epochs, our bolometric corrections have the largest uncertainties because during this time the UV flux—which is the main contribution of the luminosity—is not well constrained as noted in Bersten & Hamuy (2009). During the transition between the plateau and the radioactive tail, the diffusion approximation breaks down because the object is almost completely recombined and the photosphere is not well defined. During the tail, the bolometric luminosity is completely determined by the luminosity of radioactive decay, which makes our calculations more reliable. Although not shown here, our model shows that the shape of the LC at the end of the plateau is very sensitive to the properties of the core, such as mixing and the form of the density transition between the helium-rich core and the hydrogen-rich envelope. On the other hand, the presence of an outer atmosphere can affect the shape of the LC at the earliest stages. During these epochs, and until the hydrogen recombination sets in, the LC samples the outermost layers as the photosphere is located far out in the object (see Section 3.1.3). So, even if it was possible to find a better fit to the observations by modifying the innermost or outermost structures, such efforts are rendered meaningless due to the uncertainties introduced by our approximations.

Zoom In Zoom Out Reset image size

Figure 6 shows the photospheric-velocity evolution of our model compared with observed photospheric velocities as explained in Section 2.5. We have also included in the plot the spectroscopic velocities measured from the absorption minimum of the Fe ii λ5169 line (Jones et al. 2009). Note the very good overall agreement between model and observations. Fe ii velocities match quite well the model photospheric velocities except at the latest times. However, it is important to remark that the photospheric velocity is to be considered a good discriminator between models only at the early phases of the evolution while the object is not completely recombined. This is because, line velocities become poor photospheric velocity indicators with time, and the photosphere begins to lose its meaning in our models as the ejecta becomes nearly completely recombined at the end of the plateau.

Zoom In Zoom Out Reset image size

As shown in Figures 5 and 6 we obtain a remarkably good agreement with observations (bolometric LC and photospheric velocity evolution) for SN 1999em despite the simplifications used in our code. Note that in spite of the mismatch between the model and observations during the early LC, the photospheric velocities, which are better determined than the LC, show a good agreement at such epochs. These results are very encouraging and give us confidence in the ability of our code to infer physical parameters to study their effect on the observed quantities and ultimately to understand the physics of SNe II-P.

In the following subsections, we describe in some detail the evolution of the SN for this particular model. At a glance, the LC is distinguished by three phases: (1) an outburst followed by strong cooling, (2) a plateau, and (3) a radioactive tail. Each phase is essentially determined by the interplay between the main heating and cooling mechanisms. We thus focus our discussion on these processes along the SN evolution. As a reference, Table 1 gives a summary of some properties of the model at characteristic times during the evolution. Even though some uncertainties may arise on the detailed propagation of the SW due to the assumed mixing and energy injection, and also on the shock breakout because of the adopted radiative transport, we think that it is instructive to provide a description of the processes which occur during such phases. We emphasize that the details of the shock propagation do not directly affect the resulting LC and photospheric velocities, which we are mainly interested in modeling.

The initial phase of shock breakout is discussed in Section 3.1.1, the adiabatic-cooling phase where the homologous expansion is reached is addressed in Section 3.1.2, the cooling and recombination phase is described in Section 3.1.3, and finally the radioactive-decay processes that power the tail are explained in Section 3.1.4.

3.1.1. Shock Wave Propagation and the Early Evolution

A powerful SW begins to propagate outward through the envelope when we artificially inject energy near the center of the star (assumed to occur at t = 0). This energy is initially released as internal energy, and part of it is rapidly transformed into kinetic energy (e.g., by t = 0.5 days, the total energy is approximately equally divided between kinetic and internal energy). The velocities acquired by matter are so high that they exceed the local speed of sound, leading to the formation of a SW. The SW heats and accelerates the matter depositing mechanical and thermal energy into successive layers of the envelope until it reaches the surface, where photon diffusion dominates the energy transfer, and energy begins to be radiated away.

Figure 7 shows the effect of the SW propagation on different physical quantities (velocity, density, and temperature) inside the star. The shock front manifests as a sudden change in these quantities. As the shock moves outward, the material behind is accelerated and heated. Note that the outermost layer, with a sharp decline in density, acquires very high velocities. It represents a small fraction of the star. It is also clear from the velocity profiles that the inner parts of the object are decelerated. This deceleration is due to the interaction of the dense core with the extended hydrogen-rich envelope. The same figure also shows the changes in radius for different shells. From this we can deduce the location of the shock front at any time as that of the innermost shell, which has constant radius.

Zoom In Zoom Out Reset image size

At t = 1.36 days, the SW reaches the surface of the object, which produces the first electromagnetic manifestation of the explosion (although neutrinos and gravitational waves escape well before). The effective temperature and bolometric luminosity suddenly rise and reach their maximum values a few hours after breakout, specifically, at t = 1.47 days with values of L_peak = 3 × 10⁴⁴ erg s⁻¹ and T_peak = 1.1 × 10⁵ K. At these temperatures, the peak of the emitted spectrum is in the UV range. Also note that the color temperature is even higher than the effective temperature, as mentioned in Section 2.5.

Hereafter, the star begins to expand and cool very quickly, leading to an increase in photospheric radius and a decrease in temperature in the external layers. The bolometric luminosity abruptly decreases but, according to the decrease in effective temperature and the consequent shift of the emission peak to longer wavelengths, the luminosity in the optical range increases. As a result, a sharp peak in bolometric luminosity and temperature is produced, as shown in Figure 8. For example, in our model we obtain a decrease of 1.5 dex in luminosity only 6.8 hr after peak brightness. During this time the total energy radiated is 2.2 × 10⁴⁸ erg, emitted essentially as a UV flash. The short duration of the breakout explains why so few SNe II-P have been observed during this phase.

Zoom In Zoom Out Reset image size

At temperatures as high as those left by the passage of the SW, the stellar matter is completely ionized, which implies that the breakout is accompanied by a strong increase in opacity. The position of the photosphere during the outburst nearly coincides with the outermost shell. This behavior continues until the onset of recombination. Therefore, at early stages prior to recombination, the velocity of matter in the photospheric position samples the very high velocities of the outermost layers and reaches values close to 1.2 × 10⁴ km s⁻¹.

From the energetics point of view, by t = 0.5 days the total energy is approximately divided in equal proportions between internal, dominated by radiative contribution, and kinetic components. Short after breakout, the kinetic energy completely dominates the energetics (the energy radiated away at any given time is less than 1.5% of the total energy).

It is possible to estimate the average velocity during the SW propagation: given that the SW takes 1.36 days to emerge and considering a radius of R = 5.5 × 10¹³ cm (800 R_☉), we obtain an average speed of v = 4680 km s⁻¹. It is also interesting to calculate the expected time for breakout using the analytic expression given by Shigeyama et al. (1987):

$$\begin{equation} t_{\rm bk} \simeq 1.6 \left(\frac{R_0}{50 \;R_\odot } \right) \times \left[ \left(\frac{M_{\rm ej}}{10 \, M_\odot } \right) \bigg/ \left(\frac{E}{1 \times 10^{51} \,\,{\rm erg}} \right) \right]^{1/2} \,\, {\rm hr}, \end{equation} \tag{ 5 }$$

where R₀ is the initial radius, M_ej is the ejected mass, and E is the explosion energy. Using the values for our initial model we obtain t_bk = 1.26 days, in good agreement with our numerical calculation.

3.1.2. Adiabatic Cooling and Homologous Expansion

The breakout is followed by a violent expansion, resulting in the cooling of the outermost layers. During expansion, only a small fraction of the photon energy can diffuse into the surroundings. Therefore, it is possible to consider the cooling process to be approximately adiabatic and this approximation remains valid while the timescale for radiation diffusion is much longer than the expansion timescale.⁷ Note that there are two mechanisms to cool an SN: (1) loss of photons or diffusion cooling and (2) its own expansion or "adiabatic cooling." If one of these processes dominates we say that the cooling is carried out by that process. In our model, more than 90% of the decrease of internal energy is due to adiabatic cooling up to 18 days after the explosion, when the first layers with neutral hydrogen appear. After that, diffusion cooling begins to be significant and the adiabatic cooling phase comes to an end.

The internal temperature (and internal energy) decreases almost adiabatically, i.e., proportional to r⁻¹ due to the dominance of the radiative term, and quickly reaches a value near the recombination temperature of hydrogen. Quantitatively, the effective temperature goes from values close to 10⁵ K at the time of the burst when the matter is totally ionized to T_eff = 10⁴ K five days later. At these temperatures hydrogen begins to recombine. At the same time, the luminosity reaches L = 1.9 × 10⁴² erg s⁻¹ and the photospheric radius rapidly increases to R_ph = 5.2 × 10¹⁴ cm. After that, the luminosity decreases slightly as a result of the slower decrease in temperature and continuous increase in radius. By day 18, when the first layers of neutral hydrogen appear, we have T_eff = 6960 K, L = 1.2 × 10⁴² erg s⁻¹, and R_ph = 9.1 × 10¹⁴ cm. Note the slight change in luminosity, of only 0.14 dex, between day 7 and 18.

A few days after breakout, the acceleration of the material comes to an end and the expansion becomes homologous. The matter reaches velocities of the order of 1.2 × 10⁴ km s⁻¹ in the outermost layers and the object enters a state of free expansion where forces of pressure and gravitation do not have any dynamical effect on the system. The homologous regime is characterized by a constant velocity in each layer, a linear growth of the radial coordinate with time (r ∝ t), and a density distribution decreasing with time as ρ ∝ t⁻³. This behavior is clearly shown in Figure 7 for t > 3 days and in Figure 9. The nearly constant shape of the density and temperature profiles at late times are a result of the expansion being approximately homologous. This is also evident in the linear behavior of the radial coordinates for different mass shells. The condition of constant velocity is very nearly satisfied for each shell, although the transition between acceleration and homologous expansion happens at slightly different times for different mass shells (see Figure 9).⁸

Zoom In Zoom Out Reset image size

As mentioned in Section 3.1.1, at the time of breakout the internal energy of the envelope is small compared to the kinetic energy. Moreover, the object expands by a factor of ∼17 to a radius of R_ph ∼ 9 × 10¹⁴ cm, before it becomes transparent. Thus, most of the internal energy left by the shock passage is degraded by the adiabatic expansion, and only a small residual will be released as observable radiation in the following phase.

3.1.3. Cooling and Recombination Wave

The appearance of regions with neutral hydrogen determines the onset of an RW. The duration of the RW is related to the total hydrogen mass and how deep hydrogen has been mixed in the initial model. As time goes on, hydrogen recombination occurs at different layers of the object as a wave propagates inward (in Lagrangian coordinate; see Figure 10, top left). Because the opacity is dominated by electron scattering, it strongly decreases outward of the recombination front (Figure 10, top right), increasing the transparency of these layers and allowing the radiation to easily go away. Consequently, the internal energy (mostly of the radiation field) is efficiently radiated away, and the temperature drops sharply from ∼10,000 K to ∼5500 K at the recombination front (Figure 10, bottom left). In other words, a cooling wave associated with the transparency and induced by an RW propagates through the envelope. This is usually called "cooling and recombination wave" (CRW).

Zoom In Zoom Out Reset image size

In our model, the recombination begins approximately at day 18. After that, a CRW develops which moves inward in mass until all the matter is completely recombined by day 120 when the recombination front arrives at the innermost layers of the H-rich envelope (at m = 2 M_☉ for our model), and the luminosity suddenly drops, defining the end of this phase. The propagation of the CRW is clear from a glance at Figure 10, where the evolution of the fraction of ionized hydrogen and temperature profiles as a function of mass for selected times are shown. Note also the behavior of the opacity (see Figure 10, top right), which depends strongly on the ionization state of the matter. In our model, a strong drop in the opacity of the outer layers is seen at day ∼18, which leads to a considerable decrease in the optical depth of these layers and an inward drift in mass of the photosphere—defined at a fixed optical depth of τ = 2/3. The photosphere begins to follow the CRW, as shown in Figure 10 where the points indicate the photospheric position. It is important to note that we define the position of the photosphere using only continuum opacity sources, i.e., excluding the opacity floor (see Section 2.5).

The photosphere, as defined here, is nearly coincident with the outer edge of the CRW. Therefore, it is just the location of the photosphere with respect to the CRW what sets the value of the photospheric temperature close to the temperature of hydrogen recombination (T ∼ 5500 K). However, note that the effective temperature does not necessarily have a constant value. Since the effective temperature is defined as T⁴_eff = L/πσR²_ph and the luminosity remains nearly constant outside the recombination front (see Figure 11), any changes in effective temperature are related to changes in photospheric radius.

Zoom In Zoom Out Reset image size

The CRW divides the object into two distinct regions: (1) an inner zone which is hot, optically thick, and ionized and (2) an outer zone which is relatively cold, optically thin, and completely recombined. Matter in the inner zone is opaque: radiative transfer is too inefficient to produce any appreciable flow of energy (this would be strictly fulfilled if ⁵⁶Ni were confined to the innermost layers; see the discussion) and the matter cools down almost adiabatically. On the other hand, the external layers are transparent and practically do not radiate. Therefore, it is within the CRW where almost the entire radiant flux is released. When matter passes through the CRW, particles begin to be cooled by radiation, emitting more light than they absorb, and the radiant flux increases. This way, the radiative flux emerging from the CRW front carries away internal energy of the matter that is cooled by the wave. More details on the properties of the CRW are given by Grassberg et al. (1971).

Note that the bulk of the radiation does not diffuse to the photosphere; the photosphere instead moves inward, allowing the radiation to escape sooner than it would for a photosphere fixed at the outer boundary of the ejected mass. That is, the recombination process is responsible for the energy release during this phase. It should be noted that most of this energy comes from the energy deposited by the SW and not from the recombination itself. In order to test the previous statement, we ran our code without including the energy released by recombination of ions with electrons. No appreciable change was found; quantitatively, the differences between both calculations are less than 0.04 dex during all the evolution. We have included a figure in the online edition (Figure 17) showing this comparison.

Note, however, that in model m19r8e1.25ni56 where we assumed an extended ⁵⁶Ni mixing, the flux of energy inside the CRW is not negligible (see the left panel of Figure 11). This is due to the fact that the photosphere meets regions with radioactive material earlier than in the case where ⁵⁶Ni is confined to the innermost layers. Therefore, the energy deposited by radioactive decay provides additional power for the LC. A different behavior is seen when ⁵⁶Ni is confined to the innermost layers (inside 2.5M_☉; see the right panel of Figure 11). In this case, the statement that the energy flux inside the CRW is very small is fulfilled at least until day ∼110 where the radiative diffusion of the radioactive decay from the central layers begins to dominate. At earlier times, it is possible to see the outward diffusion of radioactive energy. Note also that for both models the luminosity outside the front is nearly constant.

The assumption of extended mixing was necessary in order to obtain a plateau as flat as the one observed for SN 1999em. Figure 12 shows a comparison of the bolometric LC for three cases: (1) with extended ⁵⁶Ni mixing, (2) no mixing, and (3) without ⁵⁶Ni. For case (1), ⁵⁶Ni begins to affect the LC by day ∼35 while for case (2) the effect of nickel heating is delayed until day ∼75. On the other hand, in the former case there is a less direct energy deposition at the center of the object, and the plateau declines earlier and steeper to the tail than in case (2). Our assumption of mixing of ⁵⁶Ni into the hydrogen envelope is not unreasonable as shown in studies of SN 1987A (Shigeyama et al. 1988; Woosley et al. 1988; Arnett 1988; Blinnikov et al. 2000, among others). However, it is important to mention that Utrobin (2007) found an excellent agreement with observations of SN 1999em by confining ⁵⁶Ni to the innermost layers and good agreement with observations of SN 1987A assuming moderate ⁵⁶Ni mixing (Utrobin 1993, 2004).

Zoom In Zoom Out Reset image size

The evolution of the velocity of matter at the photospheric position v_ph, the photospheric radius (R_ph), and the mass above the photosphere (M_ph) are shown in Figure 13. Initially, the photospheric velocity evolves rapidly, the photospheric position follows a linear behavior, and there is very little mass above the photosphere. Later on, when the recombination sets in, R_ph begins to differ noticeably from the linear behavior, the mass above the photosphere increases, and v_ph decreases because it samples increasingly inner, slower material. Finally, when all the matter is recombined, R_ph and v_ph sharply turn down and the luminosity undergoes a rapid decrease. Note, however, that the photospheric radius does not drop immediately. This is due to the heating caused by radioactive decay which produces some ionization of the gas. On the other hand, the luminosity undergoes a rapid decrease to values close to the luminosity of radioactive decays. If there was no ⁵⁶Ni (case 3), the SN luminosity would abruptly vanish at this point, as shown with a long-dashed line in Figure 12.

Zoom In Zoom Out Reset image size

In conclusion, the CRW has a duration of ∼100 days. The luminosity experiences a small change during the CRW propagation. In order to produce a plateau as flat as that observed for SN 1999em, we need to invoke Ni mixing in the H-rich envelope. Thus, the plateau can be seen as a combination of CRW properties plus some additional energy provided by radioactivity. Defining, from an empirical point of view, the plateau phase as the period of time when the luminosity remains constant within 0.5 mag of the value at day 50, its duration is 103.5 days, and the total energy emitted during this phase is ∼1.12 × 10⁴⁹ erg. Note that with such a definition the plateau phase includes the final stages of the adiabatic cooling phase as well.

3.1.4. Radioactive Tail

The late behavior of the LC (at t > 130 days) is dominated by the energy released from radioactive decay. Without radioactive material, the luminosity would abruptly vanish when hydrogen gets completely recombined as shown in Figure 12. Instead, the observed LC decreases to values close to the instantaneous rate of energy deposition by the radioactive decay:

$$\begin{equation} L\sim 1.43 \times 10^{43} M_{\rm Ni}/M_\odot \,\,\exp (-t/111.3). \end{equation} \tag{ 6 }$$

The decline of the LC at t > 130 days is nearly linear, in concordance with the exponential decline of the radioactive decay law. That is, the diffusion time for optical photons becomes small enough to radiate away the decay energy instantaneously, while the gamma-ray optical depth is still sufficiently large in order to allow a nearly complete local deposition of the decay energy. The bolometric luminosity in this part of the LC is a direct measure of the Ni mass synthesized in the explosion.

In this section, we compare our results with those of previous hydrodynamical studies of SN 1999em. As mentioned in Section 3, there are two previous hydrodynamical studies of SN 1999em: one given by Baklanov et al. (2005, BBP05, hereafter) and another by Utrobin (2007, U07, hereafter). In Table 2, we summarize the physical parameters obtained for SN 1999em in the three studies along with the distance, explosion time, and the chemical composition assumed in each model. It is interesting to note that the three studies employ initial density distributions and chemical composition distributions of non-evolutionary models. In our work, we assume H and He profiles very close to those adopted by U07, yet a very different ⁵⁶Ni distribution. While we assume an extended, uniform mixing of Ni (see Figure 4), U07 confined ⁵⁶Ni to the innermost layers (see their Figure 2). In the case of BBP05, H and He were assumed to be uniformly mixed throughout the envelope and a radial distribution of ⁵⁶Ni up to ∼15 M_☉ was adopted (see their Figure 9). Therefore, our ⁵⁶Ni distribution and the resulting ⁵⁶Ni mass are in better agreement with those of BBP05 than those of U07. We also note that the mass of the compact remnant—which is left out of the calculations—is different in all three works. We assume a compact remnant of 1.4 M_☉, while U07 assumed 1.58 M_☉ and BBP05 removed everything within a radius of R_C = 0.1 R_☉.

The approaches used in each work are very different. U07 used a hydrodynamical code with a one-group approximation for the radiative transport, including non-LTE treatment of opacities and thermal emissivity, non-thermal ionization, and expansion opacity. Their calculations, although involving a more sophisticated radiative transfer treatment, yielded bolometric LCs as in our code. However, the procedure to produce the observed bolometric LC for SN 1999em was quite different in both cases. While U07 based their calculations on the integration of UBVRI photometry with a constant value of 0.19 dex added to take into account the infrared luminosity, we derived UBVRIJKLM bolometric luminosities from color-dependent bolometric corrections. As noted by Bersten & Hamuy (2009), the inclusion of the L and M bands in the calculation of the bolometric correction during the tail phase is very important and it most likely leads to the difference in the nickel mass estimated in this work and in U07. BBP05, in turn, used a multi-group hydrodynamical code which allows the calculation of LC in different photometric bands. This code also included the expansion opacity effect. Therefore, BBP05 were able to compare their model with UBVRI LC separately without being affected by the uncertainties in the calculation of bolometric luminosities.

As shown in Table 2, the parameters yielded by our calculations are intermediate between those estimated by BBP05 and U07. The largest differences are 2.6 M_☉ in mass, 500 R_☉ in radius, 0.3 foes in energy, and 0.045 M_☉ in ⁵⁶Ni mass. It is quite satisfactory that our simple prescription for the radiative transfer yields physical parameters similar to those obtained from more sophisticated codes. However, it should be taken into account that (1) the models and methods used in each study are quite different, (2) there are differences in assumed quantities, such as distance, explosion time, mass cut, photometry, and (3) there may be a degree of degeneracy among explosion energy, initial radius, and initial mass within each model. Therefore, the present comparison does not help to provide a clear assessment of the validity of our physical assumptions and of the predictive power of our model.