Bursting SN 1996cr’s Bubble: Hydrodynamic and X-ray Modeling of its Circumstellar Medium
SN1996cr is one of the five closest SNe to explode in the past 30 years. Due to its fortuitous location in the Circinus Galaxy at 3.7 Mpc, there is a wealth of recently acquired and serendipitous archival data available to piece together its evolution over the past decade, including a recent 485 ks Chandra HETG spectrum. In order to interpret this data, we have explored hydrodynamic simulations, followed by computations of simulated spectra and light curves under non-equilibrium ionization conditions, and directly compared them to the observations. Our simulated spectra manage to fit both the X-ray continuum and lines at 4 epochs satisfactorily, while our computed light curves are in good agreement with additional flux-monitoring data sets. These calculations allow us to infer the nature and structure of the circumstellar medium, the evolution of the SN shock wave, and the abundances of the ejecta and surrounding medium. The data imply that SN 1996cr exploded in a low-density medium before interacting with a dense shell of material about 0.03pc away from the progenitor star. We speculate that the shell could be due to the interaction of a blue supergiant or Wolf-Rayet wind with a previously existing red supergiant (RSG) wind. The shock wave has now exited the shell and is expanding in the medium exterior to it, possibly the undisturbed continuation of the dense RSG wind. The narrow lines that earned SN 1996cr its IIn designation possibly arise from dense, shocked clumps in the CSM. Although the possibility for an LBV progenitor for this Type IIn SN cannot be completely excluded, it is inconsistent with much of the data. These calculations allow us to probe the stellar mass loss in the very last phases ( years) of a massive star’s life ( years) , and provide another means to deducing the progenitor of the SN.
keywords:circumstellar matter; methods: numerical; techniques: spectroscopic; supernovae: individual: SN 1996cr; stars: winds, outflows; X-rays: individual: SN 1996cr
Core-collapse supernovae (SNe) arise from stars with zero-age main-sequence masses . However, decades of research has failed to determine a direct relationship between the various SNe types (as classified based on their optical spectra and light curves) and the progenitor stars that gave rise to the SNe. It is not clear how the type IIP, IIL, IIb, IIn, Ib and Ic types relate to the properties of their respective progenitor stars. The progenitors themselves, as it turns out, are not very well known, with the handful of identified ones appearing to be associated almost exclusively with Type IIP SNe (Smartt, 2009). It had been surmised that the main progenitors of core-collapse SNe were red supergiants (RSGs) and Wolf-Rayet (W-R) stars (Falk & Arnett, 1977; Podsiadlowski, 1992). The explosion of SN 1987A revealed that blue supergiants (BSGs) could also be SN progenitors (Sonneborn et al., 1987), perhaps in a binary system (Morris & Podsiadlowski, 2007; Podsiadlowski et al., 2007). In the last decade there has been discussion of LBV stars being the progenitors of Type IIn SNe (Chu et al., 1999; Salamanca, 2000; Kotak & Vink, 2006; Vink, 2008; Smith, 2008; Trundle et al., 2008, 2009; Gal-Yam & Leonard, 2009), which is problematic because stellar theorists have mainly placed the LBV stage as an intermediate post-main sequence stage, not as a pre-explosion phase (Schaller et al., 1992; Langer, 1993; Langer et al., 1994; Stothers & Chin, 1996; Garcia-Segura et al., 1996; Maeder et al., 2005; Maeder & Meynet, 2008).
The problem is clear - we rarely know the progenitor star that led to a SN explosion, because it has to be typed from pre-explosion images, leading to significant ambiguity and potential bias. The expansion of the SN shock wave and the resulting emission due to circumstellar interaction (Chevalier & Fransson, 1994) opens up another window into the exploration of the pre-SN star. The thermal emission from this interaction, including the X-ray and optical emission, and to some extent the non-thermal radio emission (Chevalier, 1982b), depends directly on the external density. Thus an accurate analysis and interpretation of this emission acts as a probe of the density profile. In the case of core-collapse SNe, which lose a significant amount of mass prior to collapse, the surrounding medium is formed by material from the pre-explosion star. Decoding the structure of this circumstellar medium therefore will allow us to probe the mass-loss parameters of the pre-SN star, which can then be linked to the stellar parameters. Thus this provides a way to explore the stellar parameters even after the star ceases to exist, and allows us to probe the pre-SN properties of classes of SN that have heretofore lacked progenitor counterparts.
In this paper we pursue this method for the unusual Type IIn SN 1996cr. This SN, which exploded around 1996 but was only discovered about 11 years later, is only the second SN after SN 1987A which shows increasing radio and X-ray emission over a sustained period of a few years (Bauer et al., 2008). Bauer et al. (2008) suggested that the increasing X-ray and radio emission may be associated with a dense shell of material which the expanding SN shock wave interacts with a couple of years after explosion. Herein we carry out hydrodynamical simulations to evaluate this hypothesis, and to decipher the detailed structure of the circumstellar medium into which the shock wave from SN 1996cr is expanding. We compute the X-ray lightcurves and X-ray spectra using non-equilibrium ionization conditions, which we compare with high-resolution Chandra observations. We achieve a detailed agreement which, given that we are using a single model across multiple epochs, affirms the validity of our hydrodynamical model, and allows us to constrain a range of abundances for both the material ejected in the explosion and the surrounding circumstellar medium.
The plan of this paper is as follows: In §2 we review the observational details of SN 1996cr, including our recent High Energy Transmission Grating (HETG) spectra which motivated this analysis. §3 describes the reasoning behind our model of the circumstellar medium (CSM), and the results of the hydrodynamical modelling of the SN ejecta interacting with the CSM. Our techniques for computation of the X-ray light curves and spectra from the hydrodynamical models, and the resultant X-ray emission, are outlined in §4. §5 puts the narrow lines that earned SN 1996cr its Type IIn designation into the context of the overall model. In §6, the implications for the progenitor, and the abundance determinations, are discussed in depth. Finally, §7 summarizes our results and outlines future work in this area.
2 SN 1996cr
This luminous type IIn SN was only identified in the nearby Circinus Galaxy (3.70.3 Mpc; Koribalski et al., 2004) 11 years after it was believed to have exploded (Bauer, 2007; Bauer et al., 2008), but its remarkable evolution was fortuitously captured in archival data from HST, Chandra, XMM-Newton, Spitzer, and several ground-based optical and radio observatories taken for other purposes. Initially undetected at radio and X-ray wavelengths, the SN emission jumped by a factor of at least 300 at radio wavelengths about 2 yrs after explosion. It continued to increase modestly for 8 years before finally exhibiting signs of a potential X-ray and radio turnover in late 2008.
Fig. 1 shows the radio and X-ray evolution of SN 1996cr. Unfortunately, only upper limits exist for the first few years after the explosion. But even these are sufficient to illustrate the low level of emission in the initial years, followed by a steep rise in the radio, and a less steep but no less significant rise in the X-ray emission, which lasted for roughly a decade. This behavior is reminiscent of the more famous SN 1987A, whose emission follows a similar pattern (McCray, 2003, 2007), albeit with a factor of 1000 lower luminosity (Bauer et al., 2008). Emission from SN 1987A was explained by Chevalier & Dwarkadas (1995) as arising from the interaction of the SN shock wave with a dense HII region interior to the circumstellar ring seen in the equatorial region. Further interaction with the denser ring should result in a continuously increasing light curve in the case of SN 1987A (Dwarkadas, 2007b).
Early (2004) HETG X-ray data of SN 1996cr showed hints of the complex velocity structure, similar to that seen in the optical. This, coupled with the similarity with SN 1987A and the need to investigate such an unusual SN further, motivated a deep 485ks Chandra HETG observation of SN 1996cr in early 2009, with regular monitoring thereafter. A companion paper (Bauer et al. 2010, in preparation) describes the analysis of the HETG spectrum. In this paper we concentrate on the interpretation of the X-ray emission, and the computation of the structure of the medium into which the SN shock wave is expanding. The steep rise in the X-ray and radio, and the more recent turnover in the emission (see Fig. 1) point to an interaction of the SN shock wave with high density material, from which the shock wave has now emerged. This possibility was already indicated by an analysis of the radio spectrum earlier in Bauer et al. (2008). Herein we provide a more quantitative assessment of the various properties of the circumstellar medium based on the X-ray data.
The 2009 HETG spectrum (shown later in Fig. 6) is probably one of the most detailed X-ray spectra of a young supernova, next only to SN 1987A. The well-resolved emission lines allow us to estimate the relative abundances of the ejecta to reasonable accuracy. Combined with lower signal-to-noise HETG spectra from Chandra in 2000 and 2004, and a high signal-to-noise XMM-Newton p-n spectrum in 2001, as well as other CCD-quality X-ray data, it provides a strong basis for accurate comparison of our simulated spectra with the observations, allowing us to fine tune our models.
Any model of the medium into which SN 1996cr is expanding must be able to explain the increase in the radio and X-ray lightcurves, and the spectral details. The model that we put forth in this paper is constrained by the following observational parameters: (1) The range of explosion dates, from March 1995 to March 1996 (see Bauer et al., 2008). (2) Observed X-ray upper limits prior to year 2000 and detected fluxes thereafter (Figure 1). (3) A radio VLBI measurement that provides a size measurement, which we take to be the radius of the outer shock. This gives the radius of the outer shock as about 2.8 cm at 12 years (with a statistical error of 20%). (4) Emission-line velocities seen in the X-ray and optical spectra, ranging from 2000-5000 km s. (5) Fits to the temperature (kT) and the absorption column density at various times, obtained from the X-ray spectra at several epochs. (6) Fits to prominent lines seen in the emission spectra, which constrain CSM and ejecta abundances. (7) The radio light curve (Figure 1).
3 Ejecta and CSM Models and Hydrodynamical Calculations
In order to study the interaction of the SN ejecta with the surrounding circumstellar medium, we need to delineate the density profile of the ejecta as well as the surrounding medium. In this first paper we make the assumption that the density structure, of both the ejected material and the surrounding medium, is spherically symmetric, and therefore that a one-dimensional (1D) model will suffice. Although the complexity of line emission seen, especially in the optical spectrum, may indicate otherwise, a spherically symmetric model is a reasonable starting point, given the lack of any imaging constraints. The validity of such an assumption can only be determined by computation of the light curves and spectra resulting from these calculations. As we shall show later, the spherically symmetric assumption appears to be adequate in this regard. The companion paper by Bauer et al (2010) will explore the line shapes and emission in detail for deviations from spherical symmetry. This work will then serve as a launching pad for more complicated explorations of the density structure.
Following the work of Chevalier (1982a) and Chevalier & Fransson (1994), we denote the ejecta density dropping off as a power-law with radius (or velocity). Matzner & McKee (1999) narrows this power-law to a value of around 9 or 10; in the current work we adopt a value of 9. In order to conserve mass and energy, the power-law density profile cannot extend all the way to the origin, therefore below a certain velocity the density is assumed to remain constant. This ejecta profile model therefore has two free parameters - the normalization of the outer density profile and the transition-to-plateau velocity; alternatively, these two parameters can be expressed in terms of the ejecta energy and mass.
Note that for early times, until the reverse shock reaches the plateau region, only the single parameter of the density normalization is involved in the hydrodynamics and there is a degeneracy between ejecta mass and energy. In this case one can set the energy to a canonical value and use the ejecta mass as the free parameter. Changing the explosion energy scales this mass accordingly (Chevalier & Fransson, 1994). Our late time modeling of SN 1996cr shows that in order to match the light curve, the reverse shock must have passed into the plateau region. Fitting the X-ray spectra constrains the ejecta energy and mass to be 10 ergs and 4.45 . The latter mass value is very reasonable for the total ejected material in a core-collapse explosion (Woosley et al., 2002). We note that this gives a characteristic velocity km s which is suggestive of a stripped envelope SN (Maurer et al., 2009). In §6.1 we will discuss the connection with stripped envelope SNe in more detail.
3.1 CSM Density Profile -Power-Law Model
Elucidation of the density structure of the CSM required considerable trial and error. The simplest assumption of a CSM is one with a power-law density profile r, where for a constant density medium and for a wind. A wind whose mass-loss parameters are changing with time can result in other values of . Therefore we first explored the validity of such monotonic-with-radius models.
The light curve indicates that the X-ray emission began increasing linearly about two years after explosion. In theory it is possible to have the X-ray emission increasing with time even if the ambient density is decreasing with radius. This can be shown in the context of the self-similar model (Chevalier, 1982a) for SN evolution. Detailed results for the X-ray emission using the self-similar model, for various values of the parameters mentioned below, and investigating several different assumptions for the forward and reverse shocks, are given in Fransson et al. (1996). Herein we present a simplified analysis that encapsulates the important ideas relevant to this work.
The X-ray luminosity from a source depends on the electron density , the emitting volume and the cooling function
We assume that the density of the medium into which the SN shock wave is expanding goes as . The emission arises from a thin shell of radius at a mean radius , whose volume can be expressed as . Note that for self-similar evolution, , and therefore . For a SN in the early stages, the temperature is going to be much larger than 10K, and the cooling function is assumed to vary as (Chevalier & Fransson, 1994). For a strong shock, , and therefore the cooling function in the self-similar case.
Therefore we get
For s=2, a SN expanding in a wind with constant mass-loss parameters, this gives the well-known result that the emission decreases inversely with time, .
For s=1, we get . In the self-similar case, , and therefore we get that . The parameter for s=1 with (Chevalier, 1982a), which is greater than 0.5 for all . Therefore , and the power-law exponent is always positive, leading to an X-ray evolution that increases with time, even though the density is decreasing as r.
For s=0, . With , (0.4 being the value it would have in the Sedov-Taylor stage), the X-ray luminosity increases with time for all values of .
Theoretically, it is therefore possible to envision an increasing X-ray luminosity with a CSM whose density is constant or even decreasing with radius. The radius can then be expressed in terms of the self-similar solution. We use the self-similar solution to show that this approximation is not a good one in the current situation, for the following reasons:
The radio spectra indicate that the absorption was high even two years after explosion (Bauer et al., 2008), and that as the SN shock expanded within this region the free-free absorption continually decreased. This suggests a region of high absorption into which the shock was expanding two years after the explosion. The absorption decreased as the shock made its way through this region. It would be difficult to explain this with a constant or continually decreasing density profile.
In the discussion above the post-shock temperature is related to the shock velocity, which varies as . Thus the temperature will vary as . For example, in order to get an X-ray luminosity , roughly intermediate between the hard and soft luminosity slopes seen in Figure 1 , we can take and . Then , and from equation 3. This would mean that the temperature would vary as t, which equates to an expected temperature drop of 1.7 between 2000 and 2009. Such a temperature variation is not seen. Instead the X-ray data (Bauer et al. 2010) above 2 keV show little spectral variation, and hence little apparent temperature evolution, over 9 years. The spectrum below 2 keV does vary somewhat, although it is not clear whether this is because of decreasing absorption in the soft X-ray regime or whether it is because the spectrum is actually becoming softer.
A radius of 2.8e17 at a time of 12 years implies an average velocity of 7500 km s. The optical and X-ray spectra do not show velocities exceeding 5000 km s, and in most cases somewhat lower. If we assume spherical symmetry and an average velocity of 4000 km s from 2000 to 2008 (when spectra are available) then it suggests that the average velocity in the first 4-5 years (depending on explosion date) must have exceeded an average velocity of 12,000 km s. A power law model has shock velocity varying as , which would make it extremely difficult, if not impossible, to fit such a velocity profile for reasonable values of . Whereas a model with a lower CSM density early on and a higher density later could fit it easily. Using our earlier example, we can construct a model with and using the self-similar solution, which satisfies the constraint that the shock radius be 2.8e17 at 12 years. This model gives a velocity of 5290 km s at 12 years. While this is somewhat high, it could be considered acceptable within the error bounds. However such a model would also predict a velocity of almost 6000 km s at 8 years and 6450 km s at 6 years, which is not supported by the observations. Note that a constant density model with n=9 would have , which would make this discrepancy worse.
In a model with monotonically decreasing density, the density of the CSM would need to be significantly high throughout the evolution for the shock velocity to decelerate to a few thousand km s in a few years, as suggested by the X-ray and optical spectra. For instance, in the above model with and , the density at 5 cm would be 7.71e-21 g cm. This is about 3800 times greater than what is in our simulations presented in the next section, which would result in the X-ray emission, especially the hard X-rays, being about 10 times greater on average, significantly exceeding the upper limits derived from archival observations.
The radio emission is increasing significantly with time after about 2 years (Figure 1), which is hard to accomplish in a model with decreasing density. In particular free-free absorption by itself cannot account for this increase (Bauer et al., 2008), there still needs to be a flux increase . This can best be implemented with a density increase, as the radio emission, although non-thermal, indirectly depends on the external density (Chevalier, 1996).
The X-ray emission increases significantly after two years. A simple power-law increase in the X-ray flux, as predicted by the self-similar model, cannot match the data, which requires at least two distinct power-laws, with the upward slope becoming much steeper after 2 years. This would then require at least two different density slopes, with the density decreasing slower in the second case. When combined with the VLBI radius constraint, and the velocity and temperature constraints outlined above, it would be almost impossible to make the power-law model consistent with all the available data.
The above arguments suggest that a CSM density profile that shows a simple power-law decrease would not work. Although a self-similar solution is not applicable, many of the same arguments can be used to show that an increasing power-law density profile from the origin (the stellar surface) also would not work. Instead, a density profile where the density changes after 2 years is suggested.
3.2 CSM Density Profile - High Density Shell
The above discussion clearly argues for these basic requirements of any CSM model for SN 1996cr: a CSM whose density is low close to the star, then increases sharply with radius. This would fulfill both the velocity and density criteria outlined above, as well as lead naturally to the observed radio spectra and light curves. One way to arrive at such a distribution is with a low density wind medium followed by a thin, dense shell, as is often seen around massive stars (Weaver et al., 1977; Chu, 2003; Cappa et al., 2003; Chu, 2008), and in the nebula around SN 1987A (Blondin & Lundqvist, 1993). In this paper we consider such a distribution. The inner and outer radius of this shell, the shell density, and refinement of the CSM density structure, are obtained via an iterative procedure that consisted of adopting a value for the shell parameters, computing the hydrodynamic interaction, calculating the X-ray light curves and spectra, and re-iterating until a reasonable fit was obtained between the simulated and observed spectra and light-curves. This took on the order of a dozen iterations before convergence was obtained.
The model that we finally converged on has a shell with density around 1.28 g cm extending from 1 to 1.5 cm. The X-ray light curves and spectra are highly sensitive to the shell parameters, i.e. the inner and outer shell radii and shell density, and so manage to constrain them quite tightly. Even a 20% deviation of the shell radii from these values destroys the agreement with the lightcurve, while a factor of 1.3 in the density has about the same effect. Therefore the shell itself is strongly constrained by the X-ray data. A rearrangement of the mass in the shell, i.e. a modest change in the density profile that does not alter the total mass (Figure 2), can result in modest changes to the early-time X-ray flux, as shown in Figure 5. For the CSM inside of the shell, we have assumed a structure resemblant of a wind-blown bubble, with an initial freely expanding wind whose density drops as r ending in a wind termination shock, which is assumed to be a strong shock. Beyond this is a constant density shocked wind region. The wind density, and location of the wind termination shock, are constrained only by the density of the shocked wind region, and are obtained as described in the Appendix. For the standard case described here and shown in Figure 2, the value of the wind parameter is 6.76 , the wind termination shock is located at 1.63 cm, and the shocked wind region has a constant density of 8.15 .
The structure of the medium beyond the dense shell is difficult to estimate properly, since the shock front has not had much time to expand within it. One assumption is that the dense shell consists of fully swept-up wind material. In this case the CSM density structure resembles a wind-blown bubble formed by the interaction of a supersonic wind with the slower wind from a previous epoch; the dense shell is composed entirely of the slower wind material that has been swept-up. The mass of the shell must equal the total wind mass up to the outer radius of the shell, or , where the subscript refers to the outer wind. Since the shell outer radius is known, and the mass is defined once we have the shell boundaries and shell density, this provides us with the wind parameters . Assuming constant wind parameters, the density of the material outside the shell at any given radius can then be expressed as that of a freely expanding wind . This was one of the two assumptions that we have tried, and since it gives satisfactory results over the short period of time that the shock is interacting with this region, we have continued to use it. (The other assumption involved a constant density profile exterior to the shell. However in that case the shell mass was always found to be larger than the mass of the swept-up material given that density, suggesting a complicated origin for the shell. A constant density profile either gives too small a density close in or too large a density further out, depending on its magnitude. Furthermore, the total absorption column from any constant density profile that provides adequate X-ray emission beyond the dense shell increases too fast, and thus such a profile cannot be sustained beyond 0.4pc in any reasonable model. Since the wind profile anyway fits better, we did not explore the constant density assumption further).
The final density profile for the SN ejecta and CSM that we have used, which forms the initial conditions for the hydrodynamical simulations, is shown in Fig. 2. The thin solid lines superimposed on the constant density shell are intended to show how rearranging the mass to perturb the density in this fashion can lead to a modest change in the X-ray lightcurve (the dashed gray line in Figure 5).
3.3 Hydrodynamic Simulations
We carried out hydrodynamical simulations to compute the interaction of the SN ejecta with the CSM. The simulations were carried out using VH-1, a three-dimensional finite-difference hydrodynamic code based on the Piecewise Parabolic Method of Colella & Woodward (1984). Although cooling is included in these simulations via the inclusion of a cooling function, it was not found to play an important role in the 1D calculations. However the cooling time of the shock within the shell was of the same order of magnitude as the evolution time, and it is possible that in multi-dimensional calculations, where the growth of instabilities can give rise to dense clumps, cooling could prove to be more effective111Our X-ray spectral calculation, as well as the formation of narrow optical lines, does in fact require the presence of dense clumps, but they comprise at most 1% of the mass. In these scenarios cooling is important, and in fact the shock within the clumps is assumed to be radiative in order to account for the narrow-line emission. However this would affect only a small portion of the material. In any case, these clumps are not included in the spherically symmetric numerical hydrodynamics simulations, and they would not affect the general hydrodynamical evolution described herein.. The simulation was started with the initial conditions described above, and the SN ejecta allowed to expand into the surrounding medium.
The evolution of the various dynamical and kinematic quantities over time are shown in Fig. 3. The SN ejecta expand quickly in the low density medium (see Fig. 4) until they reach the higher density shell in about 20 months. The shock, which was expanding very rapidly up to this time, is instantaneously decelerated by the high density, and its velocity drops to about 2000 km s due to the high pressure. As seen in Fig. 3, the high pressure region between the forward and reverse shocks gets highly compressed. The collision of the ejecta with the high density shell sends a lower velocity transmitted shock into the dense shell, and a reflected shock is formed that travels back into the ejecta. The reflected shock quickly overtakes the original reverse shock, thus thermalizing the ejecta faster than the SN shock itself would have done. The transmitted shock gradually increases in velocity as the highly over-pressured region slowly depressurizes, and then attains a more or less uniform velocity of around 5500 km s.
The interaction of the shock wave with the high density shell lasts for about 4.6 years (from t 1.5 to t 6.1 years.). The structure of the interaction region resembles the interaction of a power-law profile with a constant density medium, as expected, and the density decreases from the reflected shock to the contact discontinuity. Since the pressure does not vary much, the temperature is almost inversely proportional to the density, and increases outwards from the reflected shock, until it reaches a large value at the contact discontinuity. We note that the structure is comparable to that described by the self-similar solution for a steep ejecta density profile colliding with a wall (Chevalier & Liang, 1989). The comparison is close but not exact, since the shock in this case imparts some motion to the shell, which is not included in Chevalier & Liang (1989). It is also interesting to note that, given the presence of the dense shell so close to the star, the density into which the forward shock is propagating, while in the shell, is much higher than the density into which the reverse shock is propagating. This is contrary to a general expansion in a power-law or constant density medium, and results in the temperature behind the forward shock () being initially lower than the temperature behind the reverse shock (). Therefore we expect that the major contribution to the X-ray emission in the 0.5-10 KeV wavebands at this early time arises from the shocked circumstellar material. As the shock continues to expand, the temperature behind the forward shock increases, whereas that behind the reverse shock decreases as the shock expands into the increasing steep power-law ejecta density profile. When the transmitted shock exits the shell and begins to reform in the freely expanding wind around it, the reflected shock is then interacting with the constant density ejecta, whose density is higher than that of the wind medium. The density behind the reflected shock is higher than that behind the forward shock, its temperature lower, and the emission becomes dominated by the reverse shocked ejecta. It was found that a good fit to the late-time emission required that the reverse shock enter the plateau ejecta region. The location of the plateau, the velocity below which the ejecta density remains constant, was adjusted to match the observations, fixing the total ejected mass to the value given in the previous section.
Once the transmitted shock exits the shell, it is expanding in the lower density wind, whose density decreases as r. The radius of the forward shock at 12 years is consistent with that obtained from VLBI imaging to within 5%. Thus this approximation has proved suitable so far, but only time and future data can confirm whether this approximation is valid for the further expansion of the shock wave. While most of the X-ray emission at this point is expected to arrive from the reverse shocked ejecta, the radius and velocity of the forward shock, and therefore the size of the remnant, are determined by the interaction with this wind.
It is interesting to note that although the shock is interacting with an r wind, the shock structure is resemblant of the shock interacting with a constant density medium. This is to be expected as it will take a few dynamical times for the shock to reach the self-similar solution in the wind (Dwarkadas, 2005). Over time therefore the entire density structure of the shocked interaction region is predicted to change, with the density at the contact discontinuity going from relatively small to relatively large values. Such a change should be reflected in the X-ray emission over the next few decades.
4 Computation of the X-ray Emission and Comparison with Observations
4.1 Non-equilibrium Ionization and Lagrangian Shells
Our goal is to compute the X-ray emission from the output of the hydrodynamic simulation efficiently and with reasonable accuracy. Plasma that has been recently shocked and whose density is low is presumably not in ionization equilibrium (Hamilton et al., 1983; Borkowski, 2000; Borkowski et al., 2001); as a rule of thumb this is the case if s cm where is the electron density and is the time since the plasma was shocked. Since the plasma impacted by the SN shock wave has a density typically less than 10 cm, the time to reach ionization equilibrium is measured in years and the X-ray emission must be computed using non-equilibrium ionization (NEI) conditions. In general this is done by calculating the evolution of ion fractions using ionization and recombination rates under appropriate simplifying assumptions (Hamilton et al., 1983; Borkowski et al., 1994). Here, however, our approach is to make use of existing NEI codes and appropriately interface them to the hydrodynamics output.
We use NEI emission models available in the XSPEC package, vgnei and vpshock, along with “NEI version 2.0” atomic data from ATOMDB (Smith et al., 2001) which has been updated to include relevant inner-shell processes222Files provided by K. Borkowski and online at http://space.mit.edu/home/dd/Borkowski/. These codes encode the NEI state with only a few parameters; we have chosen to sacrifice some accuracy for speed and simplicity. To calculate the X-ray emission from an individual fluid element these models need information that is based on the element’s history, in particular the ionization age, , the current electron temperature, , and the ionization-age-averaged temperature, (Borkowski et al., 2001). These quantities can be calculated from the set of simulation time steps provided the location of a given fluid element can be determined and followed in time. To this end we regroup the Eulerian output at each time step in a Lagrangian manner, using cumulative mass as the Lagrangian parameter to define and track a set of mass intervals, representing “shells” in 3-D. These mass shells, of order 50 of them, are not as finely spaced as the hydro radial bins (thousands), so we calculate shell-averaged hydrodynamic fluid parameters: pressure , density , and radial velocity , for each shell at the time steps, .
4.2 Abundances, , and
The hydrodynamic variables are converted to values relevant for the X-ray emission, , , through the abundances and the related mean mass per particle, . We assign an abundance to each shell of material: shells interior to the contact discontinuity are assigned ejecta abundances and those exterior are CSM. It is possible (and as we argue in the discussion section, most likely) that the CSM can be divided into two phases, the portion interior to the shell and the portion including the shell and exterior to it. However the evolution of the SN ejecta into the region interior to the shell is not constrained by the data. Therefore a single set of abundances is assumed for the CSM for simplicity. It is likely that the ejecta may be layered, as is expected for a massive star, and different abundance distributions should be expected for different layers. However, we have no a priori indication of any of the abundance distribution, and 4 spectra of different resolution spread over 9 years, only the last two of which are actually tracking reverse-shocked emission according to our calculations, would be unable to distinguish between abundances of various layers. Therefore as a first approximation we have considered a single set of ejecta abundances, although our code is easily capable of assigning a unique abundance to each shell of material if required.
Define as the relative abundance values with respect to a reference number-ratio abundance set, A(Z); here we use the Anders & Grevesse (1989) abundance set. The mean atomic weight can be calculated as:
where is the mass of an ion of element . Similarly the mean charge state averaged over all ions is given by:
where gives the average charge state of ion , i.e., the average number of free electrons per -ion, and T refers to the temperature in collisional ionization equilibrium. This value changes with progressive ionization, however it is dominated by the low- elements which are quickly ionized, so we approximate its value based on ionization equilibrium at a fixed temperature, K.
The mean mass per particle is then:
where here has units of mass. 333Note that it is common to informally express in implied atomic mass units, , e.g., “ = 0.5 for a fully ionized H plasma”. The densities of electrons and ions can now be obtained from the fluid density as: , and for a specific ion: .
We can now compute the average particle temperature in an element of plasma from the hydrodynamic pressure and density: . Since the shocks are collisionless, most of the post-shock energy is transferred to ions, and the electron temperature is generally a fraction of the ion temperature (Ghavamian et al., 2007). This is accounted for by setting with the value set to be as given in Ghavamian et al. (2007); the shock velocity is calculated self-consistently from . Other choices for can be easily implemented if desired. For example, we clipped to have a minimum value of 0.05 as an approximation to the results of van Adelsberg et al. (2008); however, this produced no significant changes in our calculated flux or spectra as the values of the mass-cuts are generally well above 0.1.
Finally, the value is modified to include its evolution due to Coulomb interactions: approaches 1.0 as the ionization age approaches . The dependence on shock velocity is given roughly by: , with = 4 cm s; this follows from the variation of the Coulomb equilibration time constant. For smaller values of , changes such that , e.g., as seen in Figure 4 of Michael et al. (2002). With these relations the value of is increased from its initial post-shock value when the of the mass-cut increases sufficiently.
4.3 Implementation and Embellishments
Using the above equations, for each shell at each time step we can calculate the input values needed for the vgnei model (Borkowski et al., 2001): kT, Tau, <kT>, and the norm , with the shell volume and the source distance. The vgnei abundance values are the appropriate set (ejecta, CSM) scaled so that the H abundance is 1. The calculations, from reading the hydro output files to evaluating the XSPEC spectral model, are all done with custom routines written in S-Lang, the interactive scripting language used in the ISIS package (Houck & Denicola, 2000).
The flexibility of the code also allows us to implement additional modifications to improve the fidelity of the emission calculation as needed. For most mass shells we use the vgnei model with an appropriate value based on the shell history as described above. For the case where a mass shell is partially shocked, the shock front is within the shell and a constant value is not appropriate since there is a range of starting at up to a maximum value. In this case we can calculate the shell’s emission using the vpshock model, which allows us to set a range of from a low (zero) to a high value. Inclusion of these partially-shocked shells adds some line emission of lower ion states and smoothes the light curves.
There are indications that the X-ray emission from SN 1996cr is not simply the sum of optically thin components: the as-fit decreases with observation epoch and the X-ray line shapes at high-resolution are asymmetric. These data are examined in more detail in the companion paper (Bauer et al. 2010, in preparation) but the methods to include them in the emission model are briefly summarized here. To self-consistently model the effect, we can add “internal absorption” proportional to the amount of unshocked CSM material external to the forward shock. In the case of the asymmetric line shapes it seems likely, and our hydro model predicts, that there is substantial opacity through the SN ejecta core. Hence some of the red-side emission is obscured producing asymmetric lines and reducing the observed flux. Based on simple 3D modeling (Dewey & Noble, 2009), a value of 0.75 has been used to account for the amount of flux typically obscured by the core, and a blue-shift and Doppler Gaussian broadening, proportional to the hydro shell velocity, are applied as a first approximation to the line shape from the shell.
4.4 Comparison with the Data
As the previous sections imply, the X-ray emission from a given hydro simulation depends on some additional parameters and assumptions. Perhaps the most influential of these are the elemental abundances assigned to the hydro plasmas. These abundances show up in two main ways: i) they determine the value of the plasma and hence affect the electron temperature and the number of electrons per ion (section 4.2), and ii) the abundances of the higher Z elements are responsible for the X-ray emission lines in the spectra. Because we use the XSPEC vgnei (vpshock) routine, the abundances of the 13 elements currently treated by that model are the ones relevant here: H, He, C, N, O, Ne, Mg, Si, S, Ar, Ca, Fe, Ni. These abundances are specified relative to the Anders & Grevesse (1989) values.
The amount of the lowest-Z elements, especially the ratio of H to He, is most responsible for determining the mean mass per particle and thus changing the electron temperature proportionally. Temperature changes lead to changes in the continuum shape and in the ratios of He-like to H-like line emission from higher-Z elements. Because their line emission is not clearly visible here, yet they contribute most of the X-ray continuum, we refer to the elements H, He, C, N, and O as “continuum elements”. Using only X-ray data it is hard to decide among different possibilities for these “continuum” abundances: H and He produce no lines in the X-ray range, and depending on the temperatures and amount of absorption involved, it can also be difficult to see lines to constrain C, N and even O abundances with X-ray spectra. For this reason we have little handle on the abundances of C, N, and O in SN 1996cr.
With the continuum abundances set, the abundances of the higher-Z elements can generally be varied quite a bit without effecting the value and in this sense the high-Z abundances are relatively uncoupled to the hydrodynamics. However, because line emission from these elements, Ne through Fe, can be clearly seen in our X-ray spectra these abundances are well constrained by the data. Two other “degrees of freedom” in generating the observed X-ray emission from the hydrodynamics are absorption, both foreground and internal, and clumping in the ejecta and/or wind (Chugai & Danziger, 1994). Both of these have little effect on the higher energy flux, 2–8 keV, but can produce 20% to 50% changes in the low-energy flux. In this low energy range there is some degeneracy between the CNO abundances, the , the clumping parameters, and the Fe abundance through Fe-L lines. For our multi-epoch model we have set the at the HETG-00 and HETG-09 epochs to 9.6 cm and 7.8 cm, respectively, up to 6 cm of which is likely from absorption in the Milky Way and Circinus Galaxy disks. We then interpolate to other epochs based on the amount of unshocked CSM. For the swept-up shell clumping, which likely extends over a range of values, we have chosen a density factor of 7.0 and a filling fraction of 0.0020; these are sufficient to generate the Si He-like lines seen at early times.
The output of our calculations provides the emitted X-ray flux versus time in the 0.5-2 keV band and the 2-8 keV band that we have used, as well as spectra and spectral model parameter files at each epoch. Furthermore, we have performed a detailed comparison of the line-shapes between the data and model Si and Fe lines at the HETG-09 epoch, which are described in the companion paper by Bauer et al. (2010).
Figure 5 shows the comparison between the flux calculated from our hydrodynamic simulations and the observations. The abundance list used to generate these plots is given in Table 1. The separate contribution of the shocked ejecta and the shocked CSM to the flux are also shown. The modelled flux is within 20% except for the year 2000 data where it is low by 30-50 %. As Figure 3 (panel 3, 5.02 years) shows, at this time the flux is determined by the inner-half of the dense shell. We have verified (Figure 5) that shifting 10% of the dense shell mass from the outer half to the inner half can produce the needed flux increase with little change at later times. While it is tempting to expand the shape of the dense shell into a set of basis functions and adjust them for the best light curve we feel this would be “over tuning” our simple model and limited data. Certainly in future when IXO will obtain a deep, well-sampled dataset of a similar SN, such an approach could be very fruitful.
It is clear that for about the first 7 years the shocked CSM material dominates the X-ray lines and the continuum, but after that the contribution of the ejecta increases as the forward shock exits the shell. By about 8 years the shocked ejecta becomes the primary contributor. If our approximation of a wind medium exterior to the dense shell is correct, then the shocked ejecta will remain the primary contributor to the total flux for more than 25 years. It is interesting to note that this is exactly the opposite scenario as deduced for another SN with a high X-ray flux, SN 1993J (Nymark et al., 2009; Chandra et al., 2009), where the reverse-shocked ejecta was found to dominate the X-ray emission for the first several decades. The dominant early contribution of the CSM is clearly due to the presence of the bubble and dense shell in this case.
Our method enables us to compute simulated spectra using non-equilibrium ionization, which we can compare to the observations. Observed Chandra HETG spectra are available for 2000, 2004, and 2009, and observed XMM spectra for 2001. One important point of our method is that the same hydrodynamical simulation can produce spectra that match the observed ones at 4 different epochs, without any unnecessary free parameters. Given the transition from CSM to ejecta dominance in the emission, the CSM higher-Z abundances were obtained by using fits to both the HETG 2000 and XMM 2001 spectra, while the ejecta abundances were set by fitting to our 2009 spectrum.
Figure 6 shows the comparison of the simulated and observed spectra at all the 4 epochs at which observed spectra are available. Our simulated spectra are able to fit the continuum emission very well, as well as most of the major lines (H and He-like transitions of the various lines are shown on the plot to aid the reader). We also manage to fit the line transitions in the hard X-ray band adequately. Initially, our simulated spectra did not adequately reproduce the He-like transitions of many elements, especially those that occur in the soft X-ray band below 2 keV. In order to fit these better, we assumed that there existed a component, about 1% by mass, with a density that was about an order of magnitude higher than that of the surrounding medium, and whose temperature was consequently an order of magnitude lower in order to preserve pressure equilibrium. High-density clumps may be responsible for this component. The rich and complicated line structure, at both the X-ray and optical wavelengths, is a good indicator of clumping, and the small percentage required to fit the various lines is well within the realm of possibility.
5 SN1996cr as a Type IIn - Clumps and narrow-line formation
SN 1996cr was classified as a Type IIn supernova, based on the presence of narrow lines in the VLT spectrum from 2006 (Bauer et al., 2008). In the model presented herein, by 2006 the shock had exited the shell and was expanding in the external wind region. Thus the narrow optical lines could not arise due to the shock within the shell, whose maximum velocity in any case was much higher than the FWHM of the narrow component (Bauer et al., 2008). It is possible that the actual FWHM could be even smaller and the lines even narrower, beyond the limit of resolution of that particular observation.
Could the narrow lines arise from recombination in the unshocked, ionized CSM, which we suggest (below) is a RSG wind? The velocity of the putative RSG wind would be low, , and this line would not be resolved. In order to get the observed flux in the narrow lines, the RSG wind must have a certain minimum mass. Using equation 11 in Salamanca et al. (1998), we find that in order to get the required H luminosity, a H number density in excess of 10 is required at the position of the shock, at the time the VLT spectrum was taken (Jan 2006). This is more than 3 orders of magnitude higher than the density in our simulations, suggesting that in our scenario this line could not have come from the unshocked CSM.
It is possible that the lines arise from shocked dense clumps suggested in §4.4. In order to see if this explanation is feasible, and to estimate the density, number and filling fraction of clumps, we follow the method outlined by Chugai & Danziger (1994). From the simulations, the average pressure in the shocked region at 10.5 years is about 3.9 dynes cm. If the shock transmitted into the clump has a maximum velocity of 670 km s (the approximate velocity of the narrow H line; Bauer et al., 2008) and is in pressure equilibrium, then the clump density is 8.69 g cm. Using the average value of the mean molecular weight for the CSM (see Table 1), this would mean a number density . This density is consistent with the “clumps” assumed in the previous section.
The cooling time of the shock in the clumps must be small compared to the destruction time for the clumps by Rayleigh-Taylor and Kelvin-Helmholtz instabilities. The temperature of the shock driven into the clumps, following Chugai & Danziger (1994), is K. The cooling function at that temperature, from Chevalier & Fransson (1994), has a value of 5.26 ergs cm s. Assuming that the cooling time is less than the time-scale on which the instabilities develop () gives a minimum size for the clumps of cm.
As suggested by Chugai & Danziger (1994), the narrow component arises in the cool material behind the clump shock, which radiates due to the reprocessing of X-ray emission from the hot shocked gas. As the shock cools, it will progressively radiate X-ray, then UV and then optical emission from the cooling region behind the cloud shock. An upper limit to the luminosity of the ensemble of shocked clumps is given by:
The luminosity of the narrow component of H is 2.79 ergs s. Using the values of and above gives
where is the clump size in units of 10 cm. Keeping in mind the lower limit on , we adopt cm, which gives . The total mass in clumps is about 0.025 . Note that the (3D) filling fraction of the clumps, , is quite low, and consistent with that obtained from the X-rays given the large uncertainties. The (2D) covering factor 0.027, thus indicating that they will not obscure the SN. However these clumps could add additional column density along certain sightlines.
A consistency check on the above calculations can be done by comparing the ratio of the velocities. Since the clumps are in pressure equilibrium with the interclump medium, the ratio of the interclump to clump shock velocity must be proportional to the square-root of the clump to interclump density. The density of the unshocked gas at this epoch in the simulation is 1.34 g cm. Given the density of the clumps mentioned above, the jump in density between the clumped and interclump medium is a factor of 65. If these are in pressure equilibrium, then the velocity ratios must be equal to the square-root of this, or about 8.06. Therefore, if the velocity of the shock within the clumps is 670 km s, this implies that the shock velocity must be about 5400 km s. This is close to the velocity of the CSM shock in our simulations (see figure 4), as would be expected. However, this is not the velocity of the broad H line mentioned in Bauer et al. (2008), which would be assumed to reflect the shock velocity. The inference is that either the broad H velocity has been underestimated, or that it does not arise from the shocked CSM, but perhaps from the shocked ejecta. Computing the FWHM of H depends on fitting the continuum precisely, which depends on the degree and nature of the reddening towards the SN, which is not accurately known. It is possible that the H has been underestimated, and that a more accurate fit to the H line will yield a FWHM compatible with this value.
The computed values for the number, size and density of clumps are comparable to those found by Chugai & Danziger (1994), and are consistent with the “clumps” referred to in the previous section. The notion that the narrow optical lines arise in dense shocked CSM clumps, which is the same population that is also responsible for the lower-temperature (clumped) X-ray lines, is quite appealing. However the H lines are formed at the systemic velocity, while the X-ray lines are generally shifted by 3000-5000 km s from the systemic velocity. Therefore, it appears that two separate populations may be needed. This will be investigated further in the companion paper by Bauer et al. (2010).
6.1 Implications for the SN progenitor
What does our model suggest for the hydrodynamic evolution of the progenitor star? The derived density profile, with a low-density interior surrounded by a dense shell, and a wind region exterior to it, resembles a wind-blown bubble (Weaver et al., 1977) formed by the interaction of two winds. The formation of wind-bubbles, and the evolution of a SN within such wind-bubbles, has been extensively studied in the past (Ciotti & D’Ercole, 1989; Chevalier & Liang, 1989; Tenorio-Tagle et al., 1990, 1991; Franco et al., 1991; Dwarkadas, 2005, 2007a, 2007b). These calculations have suggested that the interaction of the SN shock wave with the dense shell delineating the boundary of the wind bubble leads to increasing X-ray and radio emission. This is evidently the scenario that best resembles the emission from SN 1996cr. This is only the second case after SN 1987A where a sustained increase in the X-ray and radio light curves over a multi-year period is seen. However, the shell in the case of SN 1996cr is 6 times closer to the star than the dense equatorial ring surrounding SN 1987A. Its density in our model is a few times higher than the ring density deduced for SN 1987A. The most recent data suggest that the shock has exited the shell and the emission is turning over. Thus SN 1996cr provides an early glimpse into what the further evolution of the lightcurves in SN 1987A will look like, as the SN shock in that case is just beginning to interact with the dense equatorial ring of material (Racusin et al., 2009). In this picture, SN 1996cr represents a highly compressed, and much brighter version of SN 1987A.
If we assume that the dense shell consists only of swept-up wind material, then the mass of the shell is equal to the mass of swept-up wind material. We can apply the wind bubble scenario to this model. In this scenario a fast supersonic wind sweeps up the wind from an earlier phase into a thin, dense shell. We denote the earlier phase by subscript 1 and the later phase by subscript 2. The mass of the swept-up shell is about 0.64 . If this is made up of swept-up wind material, then M. Given that the shell outer radius is 1.5 cm, this gives for the external wind a value of = 8.5e15 g cm. If we assume that the wind velocity is 10 km s, as is appropriate for a red supergiant wind, then the mass-loss rate of the star in this phase would be , which is again reasonable for a RSG wind. Thus one possibility is that the star in this phase may have been a RSG. If we assume a larger wind velocity, then the mass-loss rate must scale accordingly. Mass-loss rates on the order of 10 have been deduced for hypergiant stars with a wind velocity of about 40 km s (Humphreys et al., 1997), and it is possible that the star in this phase could have been a hypergiant, consistent with the possibility that the progenitor was a very massive star. Even higher velocities would make it unreasonably large for any phase except perhaps the outburst stage of a Luminous Blue Variable (LBV) phase (Humphreys & Davidson, 1994). But such large mass-loss rates are unsustainable for a large period of time, and we would expect such a phase to be necessarily short-lived (Humphreys & Davidson, 1994). S Doradus type instabilities generally have a timescale of years to decades (van Genderen, 2001). Therefore the existence of an LBV phase would be apparent if the density profile were to change rapidly as the shock expanded outwards, which would show up as a strong deviation in the X-ray light curve from that suggested by us in this paper. We look forward to such measurements in future.
What about the final phase of the wind before the star exploded? Only upper limits to the X-ray emission are available for this period, thus only providing an upper limit to the density profile. In the simulations, we assume a wind with a value of = 6.76 g cm, which we refer to as the standard value of . This means that for a wind velocity of a 1000 km s, the mass-loss rate is on the order of 10 (Figure 7). The wind density is most constrained by the upper limits on the X-ray emission between days 600-900. A factor of 18 increase in the value of B delays the high-flux turn-on and increases the discrepancy for the year 2000 measurements by 20 to 30 % . Thus the maximum mass-loss rate that would be able to produce a reasonable light curve (within 30% of the observed values) and just about satisfy other constraints (for a velocity of 1000 km s) would be about 1.8 .
At the lower end we find that the density could be significantly smaller before the shock collides with the dense shell too early and exceeds the upper limits on days 700-900 (although this can be somewhat mitigated by changing the explosion date, which is uncertain by about a year). For the same 1000 km s velocity, we find a lower-limit to the mass-loss rate of about 3.26 . Thus, although we cannot exactly constrain the inner density, we find that the best fit light curve leads to mass-loss rates between 3. and 2. , for wind velocities of 1000 km s.
At the higher end, this mass-loss rate would be compatible with a Wolf-Rayet star progenitor for the SN. This would make the SN a Type Ib/c. At the lower end, and especially for lower velocities, the mass-loss rate and wind velocity combination approaches that deduced for SN 1987A (Chevalier & Dwarkadas, 1995; Lundqvist, 1999; Dwarkadas, 2007c, b), and may further strengthen the analogy between SN 1996cr and SN 1987A. While we are not suggesting that the progenitor of SN 1996cr was a blue supergiant, it lies within the realm of possibility.
One implication of these low mass-loss rates is that it is unlikely that the progenitor star was an LBV star. LBVs are highly luminous, unstable stars, which can undergo episodes of dramatic mass-loss (Humphreys & Davidson, 1994), leading to the formation of circumstellar nebulae around the star. LBV nebulae are generally of order 1pc (Weis, 2001), much larger than the circumstellar nebula in our simulations. LBV stars have often been suggested as progenitors for Type IIn SNe (Gal-Yam et al., 2007; Smith, 2008; Vink, 2008; Gal-Yam & Leonard, 2009; Miller et al., 2009), based mainly on the fact that they undergo mass-loss episodes with a high mass-loss rate . LBV winds are also slower than W-R winds, with velocities on the order of 200 km s (Vink & Kotak, 2007). Assuming such a velocity would lower the mass-loss rates derived above by a further factor of 5. This would lead to mass-loss rates lying between about 6 and 4 (see Figure 7) for a presumed LBV progenitor velocity. These are significantly lower than mass-loss rates proposed for LBV progenitor stars. Mass-loss from an LBV in the quiescent stage may be consistent with the higher end of the rates deduced here (Humphreys & Davidson, 1994), but then it is impossible to distinguish this from other massive supergiant stars (Stothers & Chin, 1996) based purely on hydrodynamic considerations. In any case by assumption this would indicate that the wind bubble was not the result of an LBV outburst.
We can further constrain the wind parameters by including the shell formation time, t. As above, we assume that the shell was formed by the interaction of two winds with constant mass-loss parameters, the first wind being the earlier RSG wind, followed by the second wind phase. The shell radius is given by R cm (Dwarkadas, 1997) where is the mechanical luminosity of the second wind. With R = 1.5 cm, and g cm, we still have one more unknown, the age of the shell. For an age of 10,000 years, we obtain = 1.37 ergs s. This is on the weak side for a W-R star, but may be comparable to that for the BSG progenitor of SN 1987A. Using the standard value of above, this gives a wind velocity of 160 km s, and mass-loss rate of the second wind about 1.7 . For the lowest value of we get a velocity of 1098 km s and mass-loss rate of 3.6 , while for the highest value we get a velocity of 61 km s and a mass-loss rate of 1.17 . Values of slightly lower than the standard value are comparable to the BSG progenitor of SN 1987A, while other combinations do not seem reasonable for known stars. Thus an age of about 10 years is possible only for low values of , leading to a low velocity and low mass-loss rate progenitor. The expansion velocity of the shell R is about 5 km s (which technically requires that the RSG wind velocity be also reduced, but factors of a few are easily accommodated).
It can be seen that , and therefore if we reduce the age t to 1000 years, we get = 1.37 ergs s. This implies a mass-loss rate of 1.7 and a wind velocity of 1594 km s for the standard value. The expansion velocity of the shell now is around 48 km s. These parameters are reasonable for a W-R star. It can be seen that lowering the values of by more than a factor of 30 increases the wind velocity to 5000 km s, which is incompatible with known wind velocities. The highest value of gives a velocity of 608 km s for a mass-loss rate of 1.16 . In this case it seems that values ranging from about 3% of the standard value to 4 times the standard value give mass-loss parameters compatible with known stars.
If we reduce the age still further, thereby increasing , the wind velocity increases proportionally. Reduction in age by a factor of 2 gives a mass-loss rate of 3.4 for a wind velocity of 3188 km s for the standard value of . Values much lower than the standard value are pretty much ruled out because they give too high wind velocities. At the higher end this gives a wind velocity of 1216 km s with a mass-loss rate of 2.33 , probably reasonable for a W-R star. It is unlikely therefore that the age can be reduced by even a factor of 2 for the lowest values of . For high values of the age could possibly be reduced by a factor of 3-5. An age of 250 years with the highest possible value of gives a wind velocity of 2431 km s with a mass-loss rate of 4.67 . These may just about be compatible with W-R stars. However we note that ages less than 1000 years gives velocities for the swept-up shell that are much larger than observed shell velocities around massive stars. Observed shell velocities of nebulae around W-R stars (Cappa et al., 2003), as well as around BSG stars (Crotts & Heathcote, 1991), are all in the range of 10 km s or less. Although much larger shell velocities are seen in LBV nebulae (Weis, 2001), the mass-loss rate and velocity combinations that are possible here never resemble those seen in LBV outbursts. The conclusion is that the age of the shell lies between 10 and 10 years, with a preference for a smaller age. There is a large range over which a W-R progenitor seems viable, and a smaller range where a SN 1987A-like BSG progenitor may work. The range of values for the mass-loss and wind velocity, given the parameters used and the acceptable age of the wind-blown shell, are shown in Figure 7.
The above analysis then suggests the following scenario for the progenitor: It was a massive O star that evolved off the main sequence (MS) to become a RSG or hypergiant with a mass loss rate of . This star then entered a second, brief post-MS phase wherein it became either a BSG star, or more likely a W-R star. This second post-MS phase lasted between 10 and 10 years, with a preference towards the smaller age for a W-R progenitor, and a larger age for the 87A-like BSG progenitor. The BSG case suggests a star with a zero age main-sequence (ZAMS) of 15-30 (Woosley et al., 1988). The W-R case points to a massive ZAMS star , which went from a RSG to a brief W-R phase before exploding. The abundances derived from the X-ray fitting (§6.2) suggest a preference for the W-R case.
A final phase lasting for upto 10 years is not unlike the case of SN 1987A, where the final BSG phase is assumed to have lasted for about 20,000 years before the star exploded (McCray, 2003, 2007). The small timescale derived for the final pre-SN phase puts SN 1996cr in the category of several other recent SNe, including SN 1987A, SN 2006jc (Foley et al., 2007; Tominaga et al., 2008), SN 2006gy (Smith & McCray, 2007), SN 2008iy (Miller et al., 2009) and SN 1994W (Chugai et al., 2004), which lost about a solar mass of material a short while before the stellar explosion. In either case the current observations and our simulations are helping to shed light on the last years , or of the star’s multi-million year history.
Is it possible that the wind bubble scenario does not apply at all, that the shell is expanding ballistically, and that it was ejected in an LBV eruption? Given the speeds of the winds seen in LBV nebulae (Weis, 2001), and the fact that the shell is expanding ballistically, one would expect it to have speeds of at least 500 km s. The above analysis shows that such an eruption would have occurred less than a 100 years before the SN explosion. However, if the shell is purely wind driven, it is not clear why the interior is such low density, and why the shell is deduced to be so thick. Density contours from radiatively driven wind models of LBV nebulae can be seen in Figures 3 and 4 of Dwarkadas & Owocki (2002), and do not resemble the density structure derived herein, although it is possible that other mass-loss mechanisms may result in different density distribution. In particular the thickness of the shell (50% of the radius) implies some kind of wind interaction, and would be hard to explain in a model where the wind was expanding ballistically. Furthermore, the fact that the mass in the shell is found to be equal to the mass in the external wind of the assumed density upto the outer radius of the shell, would be an extreme coincidence if the shell was not composed of swept-up external wind material. Therefore we consider this to be a less likely possibility.
Table 1 describes the abundances that we have derived in order to compute spectra that match the observations. The statistical 1- range of values for each abundance is also given. Figure 8 graphically displays the range of abundances deduced from our X-ray fitting. The H, He, N, O values are shown with very small error bars to indicate that they are fixed (H, He) or undetermined (N,O). In the figure we also compare the abundances to those derived from LETG and HETG observations of SN 1987A (Zhekov et al., 2009), as well as those given for WCE stars by Mellema & Lundqvist (2002). These comparisons lead to the following observations: (1) The X-ray data provide very little constraint on the ratios of He to H, or , of the plasmas. Agreement with the light curve and spectra was improved with our choice of a H-rich CSM and a more He-rich ejecta. However it is likely that other possibilities, e.g., both plasmas with a 1.0, could be accommodated with some adjustments to the density profile. (2) Our values for the relative abundances of the heavier elements in the CSM are consistent with those for the CSM around SN 1987A derived from fitting the LETG and HETG spectra (Zhekov et al., 2009). Those abundances were found to be consistent with those of the LMC, and with SNRs in the LMC. In the case of the Circinus galaxy, Oliva et al. (1999) arrive at an average metallicity of -0.7 0.3 from considerations of the equivalent width of the stellar CO lines in a square of 100pc on a side. This indicates a metallicity lower than that of the LMC, which would enhance the potential for a BSG to explode. (3) The main difference between SN 1996cr and SN 1987A abundances is in the H and He that make up the “continuum”. Our fits predict H and He abundances that are an order of magnitude lower than those derived for SN 1987A. (4) The low H and He in the ejecta are suggestive of a W-R star which has expelled its H and He envelope, consistent with the analysis from purely hydrodynamical considerations above. The abundance measurements are closest to those for a WCe star as given in Mellema & Lundqvist (2002) (5) The CSM in this case, although mainly composed of H and He, is overabundant in heavy elements like S and Si. The existence of strong Si and S H-like lines, and the very strong Fe H-like line at 7 KeV in the 2001 XMM spectrum, requires a metal-rich plasma. Although RSGs are not expected to have such metal-rich winds, it is interesting to note that the RSG WOH G64 shows the presence of [SII] 6717/6731 in its spectrum (Levesque et al., 2009), while the M supergiant VY CMa shows the presence of Si I lines, and SiO maser emission (Wallerstein & Gonzalez, 2001), along with various signatures of other metal species. Therefore one possibility is that the winds are enhanced in Si and S. An alternate possibility is that the shocked CSM could actually be a mixture of CSM and ejecta material, due to mixing at the interface between the shocked CSM and shocked ejecta, which can lead to Rayleigh-Taylor “fingers” of ejecta propagating into the CSM. In order to test this possibility, we have assumed a modified CSM composition, that is made up of equal parts (by number) of the SN 1987A CSM (Zhekov et al., 2009) and SN 1996cr ejecta. Such a 50-50 mixture has a resulting abundance similar to that of our “CSM”, and, if used in place of our CSM abundance values, gives a similar light curve and spectra. Thus what we are labeling as “CSM” could be contaminated by ejecta mixing, and the derived abundance reflects this mixture. (6) The H abundance is somewhat higher in the CSM while the He abundance is higher in the ejecta. This may indicate that the star gave off most of its H layer but retained somewhat more of its He layer, pointing perhaps to a Type Ib SN. The characteristic velocity computed in §3 is also closest to that of a Type Ib SN (Maurer et al., 2009). The overall abundance of metals is still high in the CSM, suggesting perhaps that most of the H and He was lost in a much earlier stage, and is much further out in radius. This is consistent with the fact that the circumstellar bubble and dense shell with which the SN shock wave was interacting in its first decade was emitted in the very last stages of the stellar lifetime, and that much more CSM material that was released earlier is to be found further out at a larger radius. If so the SN should probe this material over the next few decades. (7) One aspect that we have not considered, and which may play a role in this scenario, is grain chemistry. RSG stars are known to have dusty winds, and stars with higher mass-loss rates are more likely to show the presence of dust (Jones, 2001). It is possible that the dense shell, if composed of RSG wind material, has a significant dust component. Under equilibrium conditions, metals can condense in dust grains, thus depleting certain species and affecting the X-ray plasma abundances. Alternatively, given the presence of the high density shell and the hot gas, which can provide ions with an energy of several keV, it is also possible that sputtering of the dust by hot gas as the shock travels through the dense shell is important. This may alter the X-ray plasma composition.
|Abund.s w.r.t. Si (AG89)||Number fraction||Percent mass|
|Ne||2.09 [1.74 – 2.48]||0.07 [0 – 0.68]||4.30e-03||8.81e-05||2.427||0.104|
|Mg||0.82 [0.70 – 0.96]||0.46 [0.13 – 0.85]||5.23e-04||1.70e-04||0.354||0.242|
|Si||1.00 [0.91 – 1.10]||1.00 [0.81 – 1.26]||5.93e-04||3.48e-04||0.468||0.576|
|S||0.92 [0.77 – 1.08]||1.43 [1.13 – 1.78]||2.49e-04||2.28e-04||0.225||0.431|
|Ar||0.82 [0.51 – 1.15]||1.22 [0 – 2.03]||5.00e-05||4.34e-05||0.056||0.102|
|Ca||0.13 [0 – 0.48]||2.82 [0 – 3.94]||5.04e-06||6.35e-05||0.006||0.150|
|Fe||0.49 [0.44 – 0.55]||0.46 [0.41 – 0.52]||3.83e-04||2.10e-04||0.605||0.693|
We have computed a hydrodynamic model that is able to explain the X-ray light curves and spectra of the young supernova 1996cr. Our intent was not to completely explore the parameter space - indeed, this is practically impossible - but to evaluate the best possible model consistent with all of the available data. Our model incorporates a dense shell of material that lies about 0.03pc from the star and extends radially for about 0.015pc, with a density of about 8 amu/cc. The density interior to the shell is quite low, although it is not well constrained by the observations. The density external to this is well approximated by a wind profile at present. The interaction of the SN ejecta with this shell and wind is able to reproduce the increasing X-ray light curves that have distinguished this supernova, and are also consistent with the behavior of the radio light curves.
We have computed the X-ray emission from our simulations, using non-equilibrium ionization conditions, and taking into account the inequality between the electron and ion temperatures. Our method is both quick and accurate444When referring to the accuracy of the flux obtained from spherically symmetric simulations of a SN which is still unresolved at all wavelengths, we mean to better than about 30%. and allows us to compute simulated spectra at every epoch, that can be directly compared with the observations. We have compared both the full spectrum as well as specific lines. This is better described in the companion paper (Bauer et al. 2010, in preparation). The flux at various epochs, and the computed X-ray spectra, agree well with the observed data, thereby validating our model.
Our hydrodynamic model suggests that the progenitor of SN 1996cr was more likely a BSG or W-R star, with a fast, low density wind. The interaction of this wind with a prior wind, possibly from a RSG phase, gave rise to the wind-blown shell with which the SN shock wave was colliding in its first decade, leading to the increasing radio and X-ray emission. The similarity to SN 1987A in many different aspects is quite striking. An LBV progenitor is inconsistent with much of the data.
Given the similarity with the medium around SN 1987A, which has a clear ring structure surrounding it, it is interesting to ask whether a ring structure, rather than a fully developed spherical shell, would be a possible scenario for SN 1996cr. In our study, we have found that a somewhat denser shell (about 50% higher density), coupled with the assumption that the emission was arising from a region of 2 steradians, would give approximately the same level of flux. Thus it is possible that a donut-shaped structure rather than a spherical shell could be accommodated, although we have not investigated the detailed spectra. It seems unlikely that a denser shell would fit because then the VLBI radius constraint would not be satisfied. Further investigations in this direction are left to a future paper as they would require far more data for comparison than is currently available.
It is possible that the further expansion of the shock may not agree with our supposition that the medium outside the shell is a RSG wind, or indeed any kind of wind. Only further data on the expansion can shed light on this, and we eagerly look forward to continually exploring the X-ray emission of this interesting and enigmatic object with the Chandra and XMM-Newton telescopes e.g. through yearly flux monitoring and and future HETG spectroscopy if feasible. Similarly, a second VLBI measurement, which will provide a second radius measurement and a direct indication of the current shock velocity, is eagerly awaited. Finally, we also have applied for further optical and IR observations that will help to create a multi-wavelength history of this SN, and enable us to understand the various properties of this fascinating object.
The relationship of SNe to their progenitor stars is quite uncertain, as mentioned in the introduction. Surveys that try to pin down the progenitors are doing a fantastic job (see Smartt, 2009, and references within), but the problem is difficult. Theoretical considerations like those used in this paper provide another completely different technique to explore the CSM around the supernova, and link it to its progenitor star. No doubt there is uncertainty in this, as in any other method that is trying to ascertain the properties of a star that no longer exists. But by using the same hydrodynamic simulation to fit multi-epoch X-ray spectra as well as light curves without adding too many free parameters, we have tried to mitigate uncertainties as far as we can. The fact that the simulations reproduce the observed shape of the continuum at each epoch without any additional parameter speaks for the validity of our model. It is only by amalgamating the results from various methods that we can hope to understand the true progenitors of the various types of SNe. In the case of a SN such as 1996cr, where the possibility of looking for a progenitor star at the same location, which faded away after the SN explosion, no longer exists, a theoretical analysis may be the only method of excavating the SN progenitor.
VVD would like to thank Poonam Chandra, Nikolai Chuagi and Don York for suggestions. We are extremely grateful to the referee, Roger Chevalier, for a thorough reading of the manuscript, and for several useful comments that have helped to improve this paper. Support for this work was provided by the National Aeronautics and Space Administration through Chandra Awards Number SAO GO8-9074X, GO9-0086A/B/D and GO0-11095A/B issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. DD was supported by NASA through SAO contract SV3-73016 to MIT for support of the Chandra X-ray Center (CXC) and Science Instruments.
8 Appendix - Calculation of the Progenitor Wind Parameters
A detailed description of the wind parameters in an interacting winds scenario, for various shock strengths and other variables, is given in Chevalier & Imamura (1983). A simplified description, as outlined below, is sufficient for our purposes.
We start with the equation for the radius of the shell in §6.1
Noting that and we get
Note that is the total mass of the emitted wind over time . This includes the mass of the freely expanding wind , which extends upto the wind termination shock () and the mass of the shocked wind, which stretches from to the inner boundary of the shell . If we assume that the density of the shocked wind is constant throughout, a reasonably good approximation, and that the wind termination shock is strong, i.e. with a shock jump of 4, then the density of the shocked wind is 4 times that of the free wind at , , and the mass of shocked wind
Therefore, with we get an equation for
This equation can be solved using standard techniques for all between 0 and to get appropriate values of .
Our calculations constrain the density of the shocked wind , so we compute the values of and that provide us with the constrained values of , and thus set limits on the progenitor wind density.
- Anders & Grevesse (1989) Anders E., Grevesse N., 1989, GeCoA, 53, 197
- Bauer (2007) Bauer F., 2007, Central Bureau Electronic Telegrams, 879, 1
- Bauer et al. (2008) Bauer F. E., Dwarkadas V. V., Brandt W. N., Immler S., Smartt S., Bartel N., Bietenholz M. F., 2008, ApJ, 688, 1210
- Blondin & Lundqvist (1993) Blondin J. M., Lundqvist P., 1993, ApJ, 405, 337
- Borkowski (2000) Borkowski K. J., 2000, in S. J. Arthur, N. S. Brickhouse, & J. Franco ed., Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 9 of Revista Mexicana de Astronomia y Astrofisica Conference Series, Nonequilibrium Ionization Models for Analysis of X-Ray Spectra of Supernova Remnants. pp 288–289
- Borkowski et al. (2001) Borkowski K. J., Lyerly W. J., Reynolds S. P., 2001, ApJ, 548, 820
- Borkowski et al. (1994) Borkowski K. J., Sarazin C. L., Blondin J. M., 1994, ApJ, 429, 710
- Cappa et al. (2003) Cappa C. E., Arnal E. M., Cichowolski S., Goss W. M., Pineault S., 2003, in K. van der Hucht, A. Herrero, & C. Esteban ed., A Massive Star Odyssey: From Main Sequence to Supernova Vol. 212 of IAU Symposium, Radio observations of interstellar bubbles surrounding massive stars. pp 596–+
- Chandra et al. (2009) Chandra P., Dwarkadas V. V., Ray A., Immler S., Pooley D., 2009, ApJ, 699, 388
- Chevalier (1982a) Chevalier R. A., 1982a, ApJ, 258, 790
- Chevalier (1982b) Chevalier R. A., 1982b, ApJ, 259, 302
- Chevalier (1996) Chevalier R. A., 1996, in Taylor A. R., Paredes J. M., eds, Radio Emission from the Stars and the Sun Vol. 93 of Astronomical Society of the Pacific Conference Series, The Circumstellar Interaction Model for Radio Supernovae. pp 125–+
- Chevalier & Dwarkadas (1995) Chevalier R. A., Dwarkadas V. V., 1995, ApJL, 452, L45
- Chevalier & Fransson (1994) Chevalier R. A., Fransson C., 1994, ApJ, 420, 268
- Chevalier & Imamura (1983) Chevalier R. A., Imamura J. N., 1983, ApJ, 270, 554
- Chevalier & Liang (1989) Chevalier R. A., Liang E. P., 1989, ApJ, 344, 332
- Chu (2003) Chu Y., 2003, in K. van der Hucht, A. Herrero, & C. Esteban ed., A Massive Star Odyssey: From Main Sequence to Supernova Vol. 212 of IAU Symposium, Ring nebulae around massive stars throughout the Hertzsprung-Russell diagram. pp 585–+
- Chu (2008) Chu Y., 2008, in F. Bresolin, P. A. Crowther, & J. Puls ed., IAU Symposium Vol. 250 of IAU Symposium, Bubbles and Superbubbles: Observations and Theory. pp 341–354
- Chu et al. (1999) Chu Y., Caulet A., Montes M. J., Panagia N., van Dyk S. D., Weiler K. W., 1999, ApJL, 512, L51
- Chugai et al. (2004) Chugai N. N., Blinnikov S. I., Cumming R. J., Lundqvist P., Bragaglia A., Filippenko A. V., Leonard D. C., Matheson T., Sollerman J., 2004, MNRAS, 352, 1213
- Chugai & Danziger (1994) Chugai N. N., Danziger I. J., 1994, MNRAS, 268, 173
- Ciotti & D’Ercole (1989) Ciotti L., D’Ercole A., 1989, AA, 215, 347
- Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
- Crotts & Heathcote (1991) Crotts A. P., Heathcote S. R., 1991, Nature, 350, 683
- Dewey & Noble (2009) Dewey D., Noble M. S., 2009, in D. A. Bohlender, D. Durand, & P. Dowler ed., Astronomical Society of the Pacific Conference Series Vol. 411 of Astronomical Society of the Pacific Conference Series, Simulation and Fitting of Multi-Dimensional X-ray Data. pp 234–+
- Dwarkadas (1997) Dwarkadas V. V., 1997, PhD thesis, UNIVERSITY OF VIRGINIA
- Dwarkadas (2005) Dwarkadas V. V., 2005, ApJ, 630, 892
- Dwarkadas (2007a) Dwarkadas V. V., 2007a, ApSS, 307, 153
- Dwarkadas (2007b) Dwarkadas V. V., 2007b, in Immler S., McCray R., eds, American Institute of Physics Conference Series Vol. 937 of American Institute of Physics Conference Series, SN Shock Evolution in the Circumstellar Medium surrounding SN 1987A. pp 120–124
- Dwarkadas (2007c) Dwarkadas V. V., 2007c, in Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 30 of Revista Mexicana de Astronomia y Astrofisica Conference Series, Supernova Explosions in Winds and Bubbles, with Applications to SN 1987A. pp 49–56
- Dwarkadas & Owocki (2002) Dwarkadas V. V., Owocki S. P., 2002, ApJ, 581, 1337
- Falk & Arnett (1977) Falk S. W., Arnett W. D., 1977, ApJS, 33, 515
- Field (1974) Field G. B., 1974, ApJ, 187, 453
- Foley et al. (2007) Foley R. J., Smith N., Ganeshalingam M., Li W., Chornock R., Filippenko A. V., 2007, ApJL, 657, L105
- Franco et al. (1991) Franco J., Tenorio-Tagle G., Bodenheimer P., Rozyczka M., 1991, PASP, 103, 803
- Fransson et al. (1996) Fransson C., Lundqvist P., Chevalier R. A., 1996, ApJ, 461, 993
- Gal-Yam & Leonard (2009) Gal-Yam A., Leonard D. C., 2009, Nature, 458, 865
- Gal-Yam et al. (2007) Gal-Yam A., Leonard D. C., Fox D. B., Cenko S. B., Soderberg A. M., Moon D.-S., Sand D. J., Li W., Filippenko A. V., Aldering G., Copin Y., 2007, ApJ, 656, 372
- Garcia-Segura et al. (1996) Garcia-Segura G., Mac Low M.-M., Langer N., 1996, AA, 305, 229
- Ghavamian et al. (2007) Ghavamian P., Laming J. M., Rakowski C. E., 2007, ApJL, 654, L69
- Hamilton et al. (1983) Hamilton A. J. S., Sarazin C. L., Chevalier R. A., 1983, ApJS, 51, 115
- Houck & Denicola (2000) Houck J. C., Denicola L. A., 2000, in N. Manset, C. Veillet, & D. Crabtree ed., Astronomical Data Analysis Software and Systems IX Vol. 216 of Astronomical Society of the Pacific Conference Series, ISIS: An Interactive Spectral Interpretation System for High Resolution X-Ray Spectroscopy. pp 591–+
- Humphreys & Davidson (1994) Humphreys R. M., Davidson K., 1994, PASP, 106, 1025
- Humphreys et al. (1997) Humphreys R. M., Smith N., Davidson K., Jones T. J., Gehrz R. T., Mason C. G., Hayward T. L., Houck J. R., Krautter J., 1997, AJ, 114, 2778
- Jones (2001) Jones A. P., 2001, Royal Society of London Philosophical Transactions Series A, 359, 1961
- Koribalski et al. (2004) Koribalski B. S., Staveley-Smith L., Kilborn V. A., Ryder S. D., Kraan-Korteweg R. C., Ryan-Weber E. V., Ekers R. D., Jerjen H., Henning P. A., Putman M. E. e. a., 2004, AJ, 128, 16
- Kotak & Vink (2006) Kotak R., Vink J. S., 2006, AA, 460, L5
- Langer (1993) Langer N., 1993, Space Science Reviews, 66, 365
- Langer et al. (1994) Langer N., Hamann W.-R., Lennon M., Najarro F., Pauldrach A. W. A., Puls J., 1994, AA, 290, 819
- Levesque et al. (2009) Levesque E. M., Massey P., Plez B., Olsen K. A. G., 2009, AJ, 137, 4744
- Lundqvist (1999) Lundqvist P., 1999, ApJ, 511, 389
- Maeder & Meynet (2008) Maeder A., Meynet G., 2008, in A. de Koter, L. J. Smith, & L. B. F. M. Waters ed., Mass Loss from Stars and the Evolution of Stellar Clusters Vol. 388 of Astronomical Society of the Pacific Conference Series, Mass Loss and the Evolution of Massive Stars. pp 3–+
- Maeder et al. (2005) Maeder A., Meynet G., Hirschi R., 2005, in T. G. Barnes III & F. N. Bash ed., Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis Vol. 336 of Astronomical Society of the Pacific Conference Series, Chemical Abundances and Yields from Massive Stars. pp 79–+
- Matzner & McKee (1999) Matzner C. D., McKee C. F., 1999, ApJ, 510, 379
- Maurer et al. (2009) Maurer I., Mazzali P. A., Deng J., Filippenko A. V., Hamuy M., Kirshner R. P., Matheson T., Modjaz M., Pian E., Stritzinger M., Taubenberger S., Valenti S., 2009, ArXiv e-prints
- McCray (2003) McCray R., 2003, in K. Weiler ed., Supernovae and Gamma-Ray Bursters Vol. 598 of Lecture Notes in Physics, Berlin Springer Verlag, Supernova 1987A. pp 219–240
- McCray (2007) McCray R., 2007, in S. Immler, K. Weiler, & R. McCray ed., Supernova 1987A: 20 Years After: Supernovae and Gamma-Ray Bursters Vol. 937 of American Institute of Physics Conference Series, Supernova 1987A at Age 20. pp 3–14
- Mellema & Lundqvist (2002) Mellema G., Lundqvist P., 2002, AA, 394, 901
- Michael et al. (2002) Michael E., Zhekov S., McCray R., Hwang U., Burrows D. N., Park S., Garmire G. P., Holt S. S., Hasinger G., 2002, ApJ, 574, 166
- Miller et al. (2009) Miller A. A., Silverman J. M., Butler N. R., Bloom J. S., Chornock R., Filippenko A. V., Ganeshalingam M., Klein C. R., Li W., Nugent P. E., Smith N., Steele T. N., 2009, ArXiv e-prints
- Morris & Podsiadlowski (2007) Morris T., Podsiadlowski P., 2007, Science, 315, 1103
- Nymark et al. (2009) Nymark T. K., Chandra P., Fransson C., 2009, AA, 494, 179
- Oliva et al. (1999) Oliva E., Marconi A., Moorwood A. F. M., 1999, AA, 342, 87
- Podsiadlowski (1992) Podsiadlowski P., 1992, PASP, 104, 717
- Podsiadlowski et al. (2007) Podsiadlowski P., Morris T. S., Ivanova N., 2007, in S. Immler, K. Weiler, & R. McCray ed., Supernova 1987A: 20 Years After: Supernovae and Gamma-Ray Bursters Vol. 937 of American Institute of Physics Conference Series, The progenitor of SN 1987A. pp 125–133
- Racusin et al. (2009) Racusin J. L., Park S., Zhekov S., Burrows D. N., Garmire G. P., McCray R., 2009, ApJ, 703, 1752
- Salamanca (2000) Salamanca I., 2000, Memorie della Societa Astronomica Italiana, 71, 317
- Salamanca et al. (1998) Salamanca I., Cid-Fernandes R., Tenorio-Tagle G., Telles E., Terlevich R. J., Munoz-Tunon C., 1998, MNRAS, 300, L17
- Schaller et al. (1992) Schaller G., Schaerer D., Meynet G., Maeder A., 1992, AApS, 96, 269
- Smartt (2009) Smartt S. J., 2009, ARAA, 47, 63
- Smith (2008) Smith N., 2008, in F. Bresolin, P. A. Crowther, & J. Puls ed., Massive Stars as Cosmic Engines Vol. 250 of IAU Symposium, Episodic Mass Loss and Pre-SN Circumstellar Envelopes. pp 193–200
- Smith & McCray (2007) Smith N., McCray R., 2007, ApJL, 671, L17
- Smith et al. (2001) Smith R. K., Brickhouse N. S., Liedahl D. A., Raymond J. C., 2001, ApJL, 556, L91
- Sonneborn et al. (1987) Sonneborn G., Altner B., Kirshner R. P., 1987, ApJL, 323, L35
- Stothers & Chin (1996) Stothers R. B., Chin C.-W., 1996, ApJ, 468, 842
- Tenorio-Tagle et al. (1990) Tenorio-Tagle G., Bodenheimer P., Franco J., Rozyczka M., 1990, MNRAS, 244, 563
- Tenorio-Tagle et al. (1991) Tenorio-Tagle G., Rozyczka M., Franco J., Bodenheimer P., 1991, MNRAS, 251, 318
- Tominaga et al. (2008) Tominaga N., Limongi M., Suzuki T., Tanaka M., Nomoto K., Maeda K., Chieffi A., Tornambe A., Minezaki T., Yoshii Y., Sakon I., Wada T., Ohyama Y. e. a., 2008, ApJ, 687, 1208
- Trundle et al. (2008) Trundle C., Kotak R., Vink J. S., Meikle W. P. S., 2008, AA, 483, L47
- Trundle et al. (2009) Trundle C., Kotak R., Vink J. S., Meikle W. P. S., 2009, in G. Giobbi, A. Tornambe, G. Raimondo, M. Limongi, L. A. Antonelli, N. Menci, & E. Brocato ed., American Institute of Physics Conference Series Vol. 1111 of American Institute of Physics Conference Series, Can LBV’s Be The Direct Progenitors of Core Collapse Supernovae?. pp 311–314
- van Adelsberg et al. (2008) van Adelsberg M., Heng K., McCray R., Raymond J. C., 2008, ApJ, 689, 1089
- van Genderen (2001) van Genderen A. M., 2001, AA, 366, 508
- Vink (2008) Vink J. S., 2008, New Astronomy Review, 52, 419
- Vink & Kotak (2007) Vink J. S., Kotak R., 2007, in Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 30 of Revista Mexicana de Astronomia y Astrofisica Conference Series, Mass loss from Luminous Blue Variables and Quasi-periodic Modulations of Radio Supernovae. pp 17–22
- Wallerstein & Gonzalez (2001) Wallerstein G., Gonzalez G., 2001, PASP, 113, 954
- Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
- Weis (2001) Weis K., 2001, in R. E. Schielicke ed., Reviews in Modern Astronomy Vol. 14 of Reviews in Modern Astronomy, LBV Nebulae: The Mass Lost from the Most Massive Stars. pp 261–+
- Woosley et al. (2002) Woosley S. E., Heger A., Weaver T. A., 2002, Reviews of Modern Physics, 74, 1015
- Woosley et al. (1988) Woosley S. E., Pinto P. A., Ensman L., 1988, ApJ, 324, 466
- Zhekov et al. (2009) Zhekov S. A., McCray R., Dewey D., Canizares C. R., Borkowski K. J., Burrows D. N., Park S., 2009, ApJ, 692, 1190