ASPHERICAL SUPERNOVA SHOCK BREAKOUT AND THE OBSERVATIONS OF SUPERNOVA 2008D

The serendipitous discovery of X-ray outburst (XRO) 080109, associated with SN 2008D, on 2008 January 9 (Berger & Soderberg 2008; Kong & Maccarone 2008; Soderberg et al. 2008) has allowed us to view a stripped-envelope core-collapse supernova (SN) from its earliest stages, at or near the moment of shock breakout from the progenitor star. The radiation burst associated with shock breakout is the first electromagnetic indicator of an SN explosion (Colgate 1968, 1974). Shock breakout occurs when radiation trapped in the vicinity of the SN shock is able to escape ahead of the shock (Klein & Chevalier 1978; Ensman & Burrows 1992; Matzner & McKee 1999). When this happens, the shock transitions from a radiation-mediated shock to a hydrodynamic shock (Katz et al. 2010). SN 2008D was a normal Type Ib SN (Modjaz et al. 2009), indicating a compact progenitor lacking a significant hydrogen envelope. Shock breakout emission from such compact progenitors may retain more information about the nature and shape of the SN driving mechanism than breakouts from larger progenitors with intact envelopes (see, e.g., Couch et al. 2009).

The discovery of XRO 080109 is described by Soderberg et al. (2008). The burst lasted about 500 s and reached a peak Swift XRT count rate of about 7 counts s⁻¹. Based on a Comptonized, non-thermal emission model, Soderberg et al. conclude that the origin of the XRO is shock breakout at a radius of about 7 × 10¹¹ cm. This radius is larger than that of the typical Wolf–Rayet star, and Soderberg et al. argue that this indicates the need for an optically thick wind around the SN progenitor. Radio observations of SN 2008D, however, imply a wind mass-loss rate too low for the wind to be optically thick at a radius of 7 × 10¹¹ cm (Soderberg et al. 2008; Chevalier & Fransson 2008). These conclusions were drawn assuming spherical symmetry and raise questions about the actual radius of shock breakout.

The X-ray spectrum of XRO 080109 can be fit reasonably well with a power law, a blackbody, or a combination of the two (Modjaz et al. 2009). It is thus reasonable to consider both thermal and non-thermal sources of emission in attempting to explain the outburst. Chevalier & Fransson (2008) posit that a thermal source with a blackbody spectrum is plausible within the uncertainties of the observations. Wang et al. (2008) argue, however, that bulk-Comptonization will scatter thermal photons to higher energies creating a power-law spectrum. Although both thermal and non-thermal emission sources may be able to explain the shape of the spectrum, neither can account for the timescale of the X-ray light curve in a spherically symmetric geometry. The characteristic time for both emission types, measured as the FWHM of the observed light curve, is the light-crossing time of the progenitor star, R/c, assuming a spherically symmetric shock breakout. The FWHM of the XRO, about 100 s (Soderberg et al. 2008; Modjaz et al. 2009), indicates a spherical breakout radius of ∼10¹² cm, greater than any plausible progenitor radius and well above where the progenitor wind could be optically thick. This is an additional contradiction that is difficult to explain with a spherical shock breakout model.

In this work, we describe the characteristics of aspherical SN shock breakout and compare these characteristics to the observations of XRO 080109/SN 2008D. Compounding observational evidence indicates that core-collapse SNe, especially those involving envelope-stripped progenitors (Type Ib/c), are not spherical (see, e.g., Wang & Wheeler 2008). In the particular case of SN 2008D, polarization measurements show that the explosion is not spherical, with dramatic asymmetries in the structure of some line-forming regions (Maund et al. 2009; Gorosabel et al. 2010). Theoretically, current models for the explosion mechanism of core-collapse SNe produce inherently aspherical shock waves (Wheeler et al. 2000, 2002; Blondin et al. 2003; Burrows et al. 2006; Obergaulinger et al. 2006; Buras et al. 2006; Burrows et al. 2007). Here we focus on core-collapse SNe driven by bipolar jets (see Khokhlov et al. 1999; Couch et al. 2009), as may arise from a magneto-rotational mechanism (Wheeler et al. 2000, 2002; Burrows et al. 2007). These models have features that may explain many of the observed features of core-collapse SNe that indicate asymmetry (Khokhlov et al. 1999; Wheeler et al. 2000; Wang et al. 2001; Höflich et al. 2001; Wang & Wheeler 2008; Couch et al. 2009). The general features of aspherical shock breakout that we discuss, however, apply to arbitrarily aspherical shocks, not just those produced by bipolar jets.

The absence of spherical symmetry dramatically modifies the observational characteristics of shock breakout and subsequent stages of emission. We assume blackbody emission in our models and we apply a detector response function appropriate for the Swift XRT and account for X-ray absorption due to neutral matter along the line of sight so that we can make a direct comparison to the observations of XRO 080109. We show that the timescale of the light curve is not set by the light crossing time of the progenitor star but by the shock-crossing time. This can account for the length of the XRO associated with SN 2008D with a Wolf–Rayet star progenitor of reasonable parameters. Further, we demonstrate that the spectral shape of aspherical shock breakouts is considerably different from that of a single-temperature, spherically symmetric blackbody even if thermal emission is assumed.

Recently Suzuki & Shigeyama (2010) have reported on their study of aspherical SN shock breakout from blue supergiant progenitors. Suzuki & Shigeyama (2010) present a semi-analytic method for calculating shock breakout light curves based on results of two-dimensional hydrodynamic simulations of aspherical core-collapse SNe. They assume, as we do in this work, that the breakout emission is thermal and calculate bolometric breakout light curves. Suzuki & Shigeyama (2010) find that the asphericity of the explosion and the angle from which the explosion is viewed determine the shapes of the resulting light curves. We find a similar result in this work. Our study, however, is targeted to explaining the observations of XRO 080109/SN 2008D and, as such, we calculate X-ray spectra and light curves that allow a direct comparison to the observations. Also, our simulations are carried out in a more compact Wolf–Rayet progenitor star.

This paper is organized as follows. In Section 2, we describe the hydrodynamic simulations of jet-driven SNe. In Section 3, we describe our method of modeling the shock breakout emission from our simulations. In Section 4, we present the resulting spectra and light curves and compare our spectral and light curve models with the observations of SN 2008D. We discuss our results and give our conclusions in Section 5.

We have carried out high-resolution hydrodynamic simulations of jet-driven core-collapse SNe using the FLASH code, version 2.5 (Fryxell et al. 2000). We use two progenitor model stars in our calculations. The first is model s1c5a from Woosley et al. (1995). This model is a non-rotating, non-magnetic evolved helium star with a pre-SN radius of 1.35 × 10¹¹ cm (1.93 R_☉) and a pre-SN mass of about 2.5 M_☉. The second model used is a stretched version of model s1c5a. For this model, we stretch the radial coordinates at each model grid point of s1c5a according to r_new = r_orig + 6.17 × 10⁻⁸r^1.7_orig. The physical variables, such as density and temperature, at each grid point are left unchanged. This function then leaves the core mass and radius approximately unchanged while extending the envelope of the progenitor. The mass and radius of this model are 6.8 M_☉ and 6.5 × 10¹¹ cm (9.3 R_☉), respectively. Both progenitors are surrounded by a wind with a mass-loss rate of 1.5 × 10⁻⁵ M_☉ yr⁻¹ and a wind velocity of 1000 km s⁻¹. The wind is assumed to be spherically symmetric, which may not be the case for a rotating progenitor. The transition from the progenitor model profile to the wind profile is made linearly over about eight computational zones. Figure 1 shows the density profiles of the two progenitor models.

Zoom In Zoom Out Reset image size

We use an equation of state (EoS) that accounts for contributions to the internal energy and pressure from radiation and gas. In the wind, the gas and radiation are not in thermal equilibrium and so radiation will not contribute to the pressure or internal energy, however our single-temperature code cannot correctly account for this using an EoS that calculates contributions from both radiation and gas. Therefore, we initially set the temperature in the wind to a small value. This is justified because for a strong shock the upstream temperature is unimportant to the downstream thermodynamics. We track seven atomic species in our simulations: ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, and ⁵⁶Ni. We employ the consistent multi-fluid advection scheme implemented in FLASH (Plewa & Müller 1999). We do not include nuclear burning. Gravity is calculated using the multipole Poisson solver with m = 0 and ℓ = 1.

All simulations are carried out in two-dimensional cylindrical geometry. In order to cover the enormous dynamic range from the inner regions of the progenitor star (∼10⁸ cm) to the shock radius several minutes after shock breakout (∼10¹⁵ cm), we have implemented a logarithmically spaced cylindrical mesh. This is achieved through a radially dependent maximum level of refinement limiter. This limiter requires that the grid spacing at radius r, Δx, not fall below some fraction of the radius r. The grid spacing then takes the form Δx>ηN⁻¹_xr, where N_x is the number of zones per block in the x-direction and η is a small number that sets the resolution scale. In effect, ηN⁻¹_x is analogous to the minimum angular resolution in spherical geometry. Additionally, we have set the maximum level of refinement anywhere on the grid to be time dependent; successively higher levels of refinement are dropped from the grid as the simulation proceeds. This has the effect of dramatically increasing the Courant-limited time step at late times, allowing the calculations to be completed in a relatively small amount of computer time. Each simulation described in this paper required approximately 3000 CPU hours to cover 10⁵ s of simulation time. This also negated the need to remap the simulation onto a new grid to continue the simulations to late times (e.g., Couch et al. 2009; Kifonidis et al. 2003, 2006).

The jets that drive the explosions are introduced as time-dependent boundary conditions at the inner boundary of the grid where we inject two identical, oppositely directed energetic flows. In order to facilitate this, an essentially spherical inner hole is excised from the two-dimensional cylindrical grid. Within this hole, the hydrodynamic solution is not calculated. A diode boundary condition was enforced at the edge of the hole (see, e.g., Zingale et al. 2002). This boundary condition is equivalent to an outflow boundary condition when the flux into the hole is positive, but the flux out of the hole is always zero. We include the gravitational effect of the mass initially residing within the hole as a Newtonian point-mass at the center of the grid, and compute the self-gravity of the gas on the grid. The mass that flows into the hole is tracked and included in the calculation of the central point-mass gravitational potential. The radius of the hole expands during the simulation, cutting out the smallest zones where the Courant condition is most limiting and ensuring that the hole radius is always resolved by a large number of zones as the maximum allowed refinement level is reduced. The jet injection velocity, v_jet, varies in time according to

$$\begin{equation} v_{\rm jet}(t) = {\left\lbrace {\begin{array}{@{}l@{\quad }l@{}}v_{\rm max}, & \mbox{$t \le 0.25\ t_{\rm jet} $}, \\ \ms v_{\rm max} \frac{4}{3} (1 - t/t_{\rm jet}), & \mbox{$ 0.25\ t_{\rm jet} < t \le t_{\rm jet}$}, \\ \ms 0, & \mbox{$t > t_{\rm jet}$}, \end{array}}\right.} \end{equation} \tag{ 1 }$$

where v_max is the maximum jet injection velocity and t_jet is the total jet injection time.

We ran a total of six simulations. Two of these simulations are spherical, non-jet-driven, explosions for comparison to the jet-driven cases. The spherical explosions are initiated in an identical manner to the jet-driven cases: injection of energetic material, except that the "jet" opening half-angle is π/2. For the four jet simulations, the opening half-angle of the jets is about π/12. The parameters of the jets are listed in Table 1. The model name labeling scheme is mMrR[cold, hot], where M is the progenitor mass to the nearest solar mass, R is the progenitor radius in units of 10¹¹ cm, and the cold or hot designates the jet parameters used, given in Table 1. For all simulations, the maximum extent of the grid is 10¹⁵ cm and the initial radius of the inner hole is 2 × 10⁸ cm, roughly the radius of the iron core of the progenitor models used. The ambient density and temperature at this inner radius for both progenitors are about 5.2 × 10⁶ g cm⁻³ and 3.3 × 10⁹ K. The jets are assumed to consist entirely of ⁵⁶Ni to facilitate the tracking of the injected jet material. We note, however, the jet parameters in some of our models, e.g., m7r6cold and m7r6hot, would predominantly freeze out into lighter nuclei (e.g., ⁴He) and not into iron group elements (see, e.g., Pruet et al. 2004). Our slower, denser jets would freeze out into the iron group. The true resulting ⁵⁶Ni fraction will then be a strong function of the proton fraction Y_e in the jet that we do not attempt to model. The parameters of the simulations were chosen so that in every case, the injected jet mass is about 0.1 M_☉, a value similar to the ⁵⁶Ni mass estimated from observations of SN 2008D (Soderberg et al. 2008; Mazzali et al. 2008; Modjaz et al. 2009). Also, for each model, except m2r1hot and m2r1sph, the ratio of explosion energy to ejecta mass is about 0.8 (10⁵¹ erg/M_☉), similar to the ratio estimated from measurements of the photospheric velocity of SN 2008D at maximum light (Soderberg et al. 2008; Mazzali et al. 2008). Model m2r1hot and m2r1sph have slightly higher explosion energy to ejecta mass ratios of about 1.4 (10⁵¹ erg/M_☉). The jet parameters for models m2r1cold and m2r1hot approximately correspond to the jet parameters used in Couch et al. (2009) for their models v3m12 and v1m12, respectively. There are 25 levels of refinement at the start of each simulation and the effective angular resolution, ηN⁻¹_x, is π/1024. The simulations are run until 10⁵ s, long after shock breakout in all cases.

Figures 2–5 show density plots of the four jet-driven explosion simulations at three epochs: when jet injection stops, initial shock breakout, and the end of the simulation. In each jet explosion simulation, the jets drive bipolar shocks that expand out from the jet injection sites along the cylindrical axis. The shocks cross in the equatorial plane and establish a dense, hot pancake of unbound material. The shocks in all cases erupt from the surface of the progenitor stars first at the poles. The shocks accelerate into the low-density wind region and sweep around the surface of the progenitor and cross again on the equatorial plane. This happens just before the original equatorial shock structure erupts from the progenitor surface. The prolate shock structure evolves toward sphericity in the wind region as the reverse shock, established by the outgoing shock colliding with the wind, sweeps up an unstable shell of ejecta.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The explosions in the smaller progenitor, models m2r1hot and m2r1cold, reach the surface of the progenitor approximately 50 s after the start of the simulations. For explosion m2r1cold, the shocks take about 30 s to cross the surface of the progenitor and collide along the equatorial plane. The shock surface-crossing time is only 20 s in model m2r1hot because the shock structure is more spherical than in m2r1cold. The shocks reach peak speeds of about 1.4 × 10¹⁰ cm s^-1 immediately following eruption from the progenitor surface and then begin to slow in the wind. At the end of the simulations, around one day after shock eruption, the pole to equator axis ratio for m2r1cold is 1.2 and for m2r1hot is 1.0.

The explosions in the larger progenitor take significantly different amounts of time to reach the surface. Model m7r6cold takes about 225 s to erupt from the progenitor poles, while model m7r6hot takes about 390 s to do the same. The time it takes the shocks to sweep across the progenitor surface is also different, taking 225 s in model m7r6cold and only 125 s in model m7r6hot. As is the case for the smaller progenitor simulations, this is because the hot-jet model is less prolate than the cold. In the larger progenitor, the shocks reach peak speeds after breakout of about 1 × 10¹⁰ cm s⁻¹.

We have carried out a resolution study using jet explosion model m2r1cold. We have run two additional simulations at resolutions of ηN⁻¹_x = π/788 and ηN⁻¹_x = π/1331. The results, compared with those of the fiducial simulation with ηN⁻¹_x = π/1024, are shown in Figure 6. The higher resolution simulation shows large-scale of north–south asymmetry. This is due to small-scale north–south asymmetries near the cylindrical axis early-on that then propagate to large scales as the simulation proceeds to later times. Artificially accelerated growth of instabilities near the axis is a well-known problem in Eulerian calculations carried out in curvilinear coordinate systems and higher-resolution simulations are more susceptible to these artificial instabilities (see, e.g., Fryxell et al. 1991). Additionally, it has been documented that the directionally split piecewise parabolic method for Eulerian hydrodynamics does not conserve symmetries in small-scale structures (see, e.g., Liska & Wendroff 2003; Almgren et al. 2006). The higher resolution simulation produces slightly higher shock velocities momentarily during shock breakout, reaching a maximum speed of 1.5 × 10¹⁰ cm s⁻¹. Within the accuracy of velocity measurements, the low-resolution simulation attains breakout shock speeds equivalent to those of the medium-resolution simulations, about 1.4 × 10¹⁰ cm s⁻¹. The shock velocities in the different resolution simulations quickly become equivalent, as is evident by the similar extent of the shock structures shown in Figure 6. Because of the influence of the numerical artifacts that appear in the high-resolution simulation, the results found in the low- and intermediate-resolution simulations are more reliable. In fact, in terms of calculating the X-ray emission from the simulations, there is negligible difference between the low and intermediate simulations, as shown in Figure 19, indicating that the gross dynamics have effectively converged at the intermediate resolution. The small-scale differences between the low- and intermediate-resolution simulations have little impact on the resulting X-ray emission.

Zoom In Zoom Out Reset image size

In order to understand the observational effects that aspherical shock breakout would produce, we simulate the X-ray emission from our simulations. We calculate the X-ray spectra at each output time from our simulations and produce X-ray light curves and integrated spectra. We also compute the estimated X-ray counts as would be detected by the Swift XRT, accounting for the XRT detector response function and photon absorption along the line of sight. In this section we discuss the details of our approach. In Section 4, we present the simulated spectra and light curves for the explosion simulations and discuss the general observable characteristics of aspherical shock breakout.

We calculate the spectra of our simulation results in the Swift XRT bandpass (0.1–10 keV) as a post-processing step. The AMR data are first merged onto a uniformly spaced two-dimensional cylindrical grid using volume-weighted averaging. In order to calculate the emission as seen from multiple viewing angles, the data are rotated prior to calculating the optical depths (see Figure 7). The electron scattering and absorption optical depths, τ_es and τ_abs, are then calculated via integration along rays directed from the assumed location of the observer. We assume the observer's line of sight to be along the cylindrical R-coordinate in the post-rotated data. This line of sight is then not normal to the simulation data symmetry axis for non-zero viewing angles, as shown in Figure 7. We define the thermalization depth, where the radiation and matter temperatures equilibrate, to be where the effective optical depth $\tau _{*} = \sqrt{3 \tau _{\rm abs} \tau _{\rm tot}} = 2/3$, and τ_tot = τ_es + τ_abs (see, e.g., Rybicki & Lightman 1986; Ensman & Burrows 1992). We assume a blackbody emission spectrum is formed with temperature equal to the matter temperature at the thermalization depth. The emission directed toward the observer is then the blackbody flux from the thermalization depth times the surface area of the thermalization depth projected toward the observer.

Zoom In Zoom Out Reset image size

The thermalization depth, and hence the emissivity, is strongly dependent on the scattering and absorptive opacities. We use multi-group opacities obtained from the TOPS database maintained by Los Alamos National Lab (Magee et al. 1995). These opacities are temperature and density dependent and use 32 photon energy groups spaced logarithmically in the XRT bandpass. We assume a chemical mixture comprised mostly of helium, as would be relevant to the Type Ib SN 2008D, with a solar mix of metals at half of the solar metal abundance (i.e., X = 0, Y = 0.992, and Z = 0.008). The inclusion of metals is critically important to obtaining accurate values for the absorptive opacities, as we discuss in Section 4.1. The TOPS database gives values for the number of free electrons per atom, N_e, the Rosseland mean opacity, and the Planck mean opacity at each temperature, density, and photon energy. We assume the total absorptive opacity, κ_abs(E_γ), to be the Planck opacity. The electron scattering opacity is κ_es = 0.1N_e cm² g⁻¹, for a predominately helium gas. The optical depths are then

$$\begin{equation} \tau _{\rm [es,abs]}(E_\gamma,z) = \int ^{R}_\infty \kappa _{\rm [es,abs]}(E_\gamma,R^{\prime },z) \rho (R^{\prime },z)\ dR^{\prime }, \end{equation} \tag{ 2 }$$

where R and z are the cylindrical radius and height, and E_γ is photon energy. A schematic diagram of the optical depth integration is shown in Figure 7. The radius of the thermalization depth as a function of photon energy and z is then

$$\begin{equation} R_{\rm *}(E_\gamma, z) = R({\rm where}\ \tau _*(E_\gamma,z) = 2/3), \end{equation} \tag{ 3 }$$

and the temperature at the thermalization depth is T_*(E_γ, z) = T(R_*).

The intensity of the X-ray emission is assumed to have a blackbody spectral energy distribution with a color temperature equal to T_*(E_γ, z), i.e.,

$$\begin{equation} I(E_\gamma, z) = \frac{2E_\gamma ^3}{h^2 c^2} \frac{1}{e^{E_\gamma /k T_*}-1}. \end{equation} \tag{ 4 }$$

The emergent flux is then the intensity multiplied by the area projected toward the observer. The projected area in cylindrical geometry for a line of sight perpendicular to the symmetry axis is just A_p(E_γ, z) = 2R_*,0(E_γ, z)dz, where R_*,0 is the transformation of R_* into the frame in which the simulation symmetry axis corresponds with the cylindrical axis. For other viewing angles the projected area becomes

$$\begin{equation} A_p(E_\gamma, z) = \pi \sin ^2{\alpha }\ \vert R_{*,0}\vert dR + 2 f \cos ^2{\alpha }\ R_{*,0} dz, \end{equation} \tag{ 5 }$$

where R_*,0 is always measured from the symmetry axis (see Figure 7), f is 0 for R_*,0 < 0 and 1 otherwise, and α is the angle between the symmetry axis and the z-axis, or equivalently the angle between the line of sight and the normal to the symmetry axis. In this way, a surface area element is treated as a cylindrical ring and we account for shadowing effects. We have assumed that the effective optical depth, τ_*, is axisymmetric. This is only approximately correct as it neglects limb-darkening effects at latitudes away from the symmetry axis. The right panel of Figure 7 shows a graphical representation of the integration for a non-zero viewing angle α.

The total specific luminosity projected toward the observer is

$$\begin{equation} L_p(E_\gamma) = \sum _z I(E_\gamma, z) A_p(E_\gamma,z). \end{equation} \tag{ 6 }$$

The total luminosity directed toward the observer is ∫L_p(E_γ) dE_γ. In the calculation of L_p(E_γ), we correct for light-travel time effects. For the arbitrary geometries we consider, this is accomplished by assuming that the observed time cadence is the same as the simulation output cadence, t^obs_n = t^sim_n = t_n. Each emitting surface area element is assumed to have a constant luminosity in the interval dt_n. The arrival time of the energy emitted by a surface area element in a given time interval, i.e., L_p(E_γ)dt_n, is calculated based on the geometry of that area element, and the amount of energy the observer would see in time dt_n is summed-up. The observer's measured luminosity at time t_n is then this energy divided by the time interval dt_n. This approach to correcting for light travel time is applicable to arbitrary geometries of emitting surfaces.

The emitting regions of our simulations are typically sampled by about 1000 lines of sight along which the specific luminosities are calculated. We restrict our analysis to the Swift XRT bandpass, 0.1–10 keV. Since the absorptive opacity is photon energy dependent, the thermalization depth is different for each photon energy group. The total spectrum is thus not a single-temperature blackbody but the superposition of many blackbodies at different temperatures and with different emission areas. In order to make a more direct comparison with the XRT observations, we convolve our model X-ray spectra with the XRT detector response function and account for X-ray absorption due to neutral matter along the line of sight. We assume a distance to SN 2008D of 31 Mpc.

We have calculated shock breakout X-ray spectra and light curves at various observer viewing angles for the four jet-driven explosion models. The light curves are presented in Figures 12–15 and the time-averaged spectra in Figures 17 and 18. All emission models are calculated assuming a neutral matter column depth along the line of sight, N_H, of 1.7 × 10²⁰ cm⁻², corresponding to the Galactic value in the direction of SN 2008D (Dickey & Lockman 1990). X-ray photon absorption due to neutral matter could occur both in the host galaxy (but well beyond the scales we simulate) and locally in the Milky Way. All models are calculated using opacities appropriate for a helium gas with a solar mix of metals with metal abundances that are half of the solar values (i.e., Z = 0.5 Z_☉ = 0.008). The influence of varying both N_H and Z is discussed in Section 4.1. The general light curve characteristics are given in Table 2.

We show the XRT data of XRO 080109 in Figures 8 and 9. To prepare these data, we downloaded Swift observation 00031081002 from the High Energy Astrophysics Science Archive Research Center⁷ and reprocessed the XRT data with the "xrtpipeline" task using the latest calibration files. We extracted events from SN 2008D using xselect and produced a spectrum and light curve. Because SN 2008D was piled up during this observation, we used an annular extraction region of 707 inner radius and 4081 outer radius, corresponding to the 40% encircled energy radius of the point-spread function (PSF) and the 90% encircled energy radius, respectively. The response files were produced using the "xrtmkarf" task, which took into account the fact that we extracted events from only 50% of the PSF. These response files were also used in the post-processing of our simulation results to compare the models directly to the 2008D data. Therefore, the simulation light curves and spectra can be thought of as having been observed with only ∼1/2 of the XRT effective area.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 12–15 show both the total X-ray luminosity and the predicted X-ray count rate, corrected for detector response, X-ray absorption, and distance. As can be seen, the shapes of the count rate curve and the luminosity curve are not the same. This is because the underlying spectrum is varying in time. At times when the spectra are softer, more photons are produced at lower energies where the effects of X-ray absorption and detector response are stronger. Thus it is not appropriate to use a constant count rate-to-luminosity conversion factor at all times.

In each explosion model, the spectrum varies significantly throughout the burst. This is true even in the spherical explosions (see Figure 10), though to a lesser extent, because as the radius of the thermalization depth increases and the temperature there cools adiabatically, the reduction in emissivity is somewhat balanced by the increased emitting surface area. Figure 16 demonstrates the temporal variability of the spectrum for model m2r1cold. This figure shows the instantaneous spectrum of m2r1cold at four times during the first 150 s of breakout emission. The thermalization layer in each explosion lies in the region between the forward and reverse shocks. At later parts of the X-ray bursts, the thermalization layer coincides roughly with the contact discontinuity between the swept-up wind and the reverse-shocked ejecta, where the density increases significantly. Therefore, the nature of the breakout emission can depend significantly on the character of the wind. The absence of a progenitor wind (or the limit of a very low density wind) may result in the thermalization layers being driven deeper into the ejecta, possibly below the contact discontinuity. Since the temperature is highly discontinuous at the contact discontinuity, this could dramatically change the emission.

Zoom In Zoom Out Reset image size

The light curves for these jet-driven explosions may be generally described as having two peaks. The first is associated with the initial shock breakout at the poles. The dip in the light curves following this first peak is due to the decreasing temperature of the thermalization layer as the shock moves out further into the wind. The second peak is attributed to emission from the equator as the bipolar shocks cross there and create a pancake of very hot, twice-shocked gas. As the explosions evolve more toward sphericity at later times, the light curves decay exponentially. The timescales of the light curves are related to the shock crossing times of the progenitors. Figure 11 shows the post-breakout evolution of the temperature and thermalization depth at 0.75 keV for model m2r1cold. Before shock breakout the location of the thermalization layer is the transition region between the wind and the progenitor star, where the density increases rapidly. Emission from these regions is, however, negligible in the XRT band due to the relatively low temperatures there. After shock breakout the thermalization layer lies in between the forward and reverse shocks, the region comprise of shocked, accelerated wind material.

Zoom In Zoom Out Reset image size

Light curves of explosion models m2r1cold and m2r1hot (the smaller progenitor) are shown in Figures 12 and 13 for three viewing angles: along the equator (0°), 45°, and along the axis of symmetry (90°). The widths of the light curves are roughly 100 s at all angles for both models and the total X-ray burst time is around 500–600 s. These values are very close to what was observed for SN 2008D, shown in Figure 8. Examining the shape of the model light curves shows that for 0° viewing angles the light curves rise very quickly. This is due to two effects. The first is simply because the emission from both poles is visible. The second is due to the nature of an aspherical breakout and the effects of light travel corrections. As the bipolar shocks continue to erupt from the progenitor surface, the brightest emission is coming from where the shocks are just reaching the surface (essentially two rings moving across the stellar surface from the poles to the equator). Thus, the brightest emitting regions are moving rapidly toward the observer causing a pile up of emission once light travel time corrections are made. This effect is reduced, or eliminated, at higher viewing angles because the brightly emitting rings are no longer moving so much toward the observer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2 lists the resulting light curve characteristics for our models m2r1cold and m2r1hot, as well as m2r1sph for comparison. As expected, the spherical explosion has a higher peak luminosity, but a much shorter FWHM and overall burst time, Δt. The width of the light curve for m2r1sph is set by the light crossing time of the progenitor. The simulated XRT count rates for models m2r1hot and, especially, m2r1cold are very similar to those we find for XRO 080109/SN 2008D. The light curve rise time and FWHM for model m2r1cold at a viewing angle of 45° (red dashed curve in Figure 12) is a good match to SN 2008D.

The X-ray spectra for models m2r1cold and m2r1hot are shown in Figure 17. These spectra are corrected for the detector response function and for X-ray absorption (assuming N_H = 1.7 × 10²⁰). They are also averaged over the burst time, Δt, as was done for the observations of XRO 080109/SN 2008D (see Figure 9; Soderberg et al. 2008; Modjaz et al. 2009). The shapes of the spectra at the different viewing angles are very similar. Indeed there is not much difference between the two different models. The time-averaging washes away the major differences that are apparent in the light curves. The shape of the spherical explosion spectra is also very similar. They are generally brighter, but this is due to a shorter averaging time. The spectra in each case are softer than the spectrum of XRO 080109 (see Figure 9).

The light curves for the simulations in the larger progenitor (Figures 14 and 15) are characterized by much longer timescales and overall less bright emission. The FWHM for models m7r6cold and m7r6hot range from around 100–200 s at various angles while the total burst times are from 300 to 1000 s (see Table 2). The peak luminosities and total radiated energies are less than in the analogous explosions in the smaller progenitor. While the ratio of explosion energies to ejecta masses is roughly equivalent across all simulations, the reduced luminosities in the larger progenitor can be explained by a slightly lower shock velocity during the burst and a lower wind density at the radius of shock breakout (i.e., the radius of the progenitor star). As we discuss in more detail in Section 4.2, the density of the wind plays an important role in the strength of the X-ray emission because the thermalization depths during the bursts lie in the shocked wind. The light curve for m7r6cold is similar in shape to m2r1cold and m2r1hot but with longer timescales. Model m7r6hot, however, exhibits a dramatically double-peaked light curve.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As discussed in Section 2, we ran simulations of model m2r1cold at three different resolutions. Figure 19 shows the X-ray light curves for these three simulations plotted together. The light curves of the low and fiducial resolution simulations are extremely similar. The high-resolution case varies only slightly from the lower resolution light curves. The early peak is brighter, due to the slightly greater shock speeds just after breakout (see 4.2). The second peak is slightly dimmer, due to the lesser degree of extension of the southern shock structure in the high-resolution simulation, leading to a smaller emitting area as seen from a viewing angle of 0°.

Zoom In Zoom Out Reset image size

4.1. Dependence on Metallicity and X-ray Absorption

The absorptive opacities depend strongly on the metallicity of the absorbing gas. To illustrate this we have calculated simulated spectra and light curves of our explosion models using metallicity values of 0.1 Z_☉ and 0.0 Z_☉ (i.e., metal-free). Figures 20 and 21 show the light curves for model m2r1cold for metallicities of Z = 0.1 Z_☉ and Z = 0. The spectra for these cases are shown in Figure 22. The lower metallicity drives the thermalization depth deeper into the explosion where the temperatures are higher, resulting in significantly increased emission (note the difference in scale in Figure 21). For the case of metal-free gas, this increase is dramatic. Due to the significantly reduced absorptive opacities, the thermalization depths are pushed down below the reverse shock and into the deep, very hot regions of the ejecta. The calculated emission for this metallicity is orders of magnitude greater than the other models and the observations of SN 2008D. The difference in the emission characteristics between our fiducial models with Z = 0.5 Z_☉ and the models with Z = 0.1 Z_☉ is less drastic. This is because for Z = 0.1 Z_☉, the thermalization depths are still above the reverse shock and the temperature in between the forward and reverse shocks does not vary greatly. The metal-free case spectra also noticeably lack the deep "absorption" features at ∼0.4 keV and ∼0.9 keV. These "absorption" features are caused by a significant increase in the absorptive opacities at these energies (due to the presence of metals), pushing the thermalization depths to larger, cooler radii. We note that our fiducial metallicity of Z = 0.5 Z_☉ is consistent, within the accuracies, with three independent measurements of the metallicity of the region around SN 2008D: Soderberg et al. (2008), Thöne et al. (2009), and Modjaz et al. (2010).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We plot XRT count rates and spectra for model m2r1cold using different values of N_H in Figures 23 and 24. If the location of SN 2008D in its host galaxy were particularly dense, X-ray absorption by neutral matter atoms in the host galaxy may be significant, warranting an N_H beyond the Galactic value of 1.7 × 10²⁰ cm⁻². Increased hydrogen column depth increases the absorption of lower energy X-ray photons. Since our simulated spectra for the non-zero metallicity cases are dominated by emission from below about 1 keV, increasing N_H dramatically reduces the resultant count rates. Due to a significantly greater amount of hard emission in the Z = 0 model, the light curve and spectrum are effected very little by an increase in N_H. It is possible to increase the neutral matter column depth so that the peak count rates for the metal-free model are similar to those observed for XRO 080109 (<10 counts s⁻¹). This requires, however, an unrealistically high value of N_H, greater than 10²³ cm⁻² and results in a very narrow light curve as the softer emission after the first 20 s is almost entirely absorbed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Color Temperature Enhancement

In our simulations, the thermalization depths, where the X-ray spectrum is formed, lie in the region of shocked wind, behind the forward shock and ahead of the reverse shock. Since we assume local emission with a blackbody spectrum, the resultant luminosity is a strong function of the temperature used in calculating the spectrum. In the region in between the forward and the reverse shock, this temperature is dependent on the forward shock velocity and the wind density. To see this, assume that enthalpy flux is conserved at the shock front such that, in the frame of the shock,

$$\begin{equation} \rho _1 v_1\big(\mbox{$\frac{1}{2}$}v_1^2 + \epsilon _1+P_1/\rho _1\big) = \rho _2 v_2\big(\mbox{$\frac{1}{2}$}v_2^2 + \epsilon _2+P_2/\rho _2\big), \end{equation} \tag{ 7 }$$

where ρ is gas density, v is gas velocity, is the specific internal energy of the gas, P is the gas pressure, and subscript 1 denotes pre-shock values and subscript 2, post-shock values. We can assume that the post-shock internal energy is dominated by contributions from radiation, ₂ = 3P₂/ρ₂ = a_radT⁴₂/ρ₂, and the density jump at the shock is ρ₂ = 7ρ₁, for a strong shock. So long as the shock is strong, the upstream internal energy is negligible, ₁ ≈ P₁ ≈ 0. Conservation of mass flux at the shock also gives v₂ = v₁ρ₁/ρ₂ = v₁/7, in the shock frame. In the frame of the progenitor star, v₁ is the shock speed, v_sh, as long as the shock is moving fast relative to the pre-shock gas. Substituting these relations into Equation (7) and solving for the post-shock temperature yields

$$\begin{equation} T_2 \approx 7.6\times 10^{5}\ \rho _{-11}^{1/4}\ v_{\rm sh, 10}^{1/2}\ {\rm K}, \end{equation} \tag{ 8 }$$

where ρ₋₁₁ is the wind density in units of 10⁻¹¹ g cm⁻³ and v_sh,10 is the shock velocity in units of 10¹⁰ cm s⁻¹, appropriate for shock breakout into a Wolf–Rayet wind. Thus, the post-shock temperature may be increased by enhancing the wind density or the energy of the explosion (which, in turn, increases v_sh).

In order to demonstrate the influence that an increased post-shock temperature has on our simulated spectra and light curves, we present the spectrum and light curve for model m2r1cold calculated by using a temperature that had been enhanced by a factor of 1.8 above the temperature found in our hydrodynamic simulations. Such an increase in temperature would result from a factor of 10 increase in ρ₁v²_sh. Tanaka et al. (2009) find that one-dimensional explosion models with kinetic energy to ejecta mass ratios of 1.4–1.7 (10⁵¹ erg/M_☉) fit the late-time spectrum and light curve of SN 2008D best. This ratio is about a factor of 2 greater than what we have used in our most simulations, including m2r1cold. This alone could account for a ten-fold increase in the shock ram pressure ρ₁v²_sh because the mass averaged-velocity will scale roughly as $\bar{v}_{\rm ej}^2 \propto E_{\rm K} / M_{\rm ej}$ and the shock speed will far exceed $\bar{v}_{\rm ej}$.

The results of this enhanced temperature calculation for model m2r1cold are shown in Figure 25. The simulated light curves and spectra are all for a metallicity of 0.5 Z_☉ and a viewing angle of 0°. Figure 25 shows the behavior of the simulated emission with increased N_H. The enhanced temperature calculation yields a peak X-ray luminosity of 4.9 × 10⁴³ erg s⁻¹ and radiates a total of 2.0 × 10⁴⁶ erg. For the case of N_H = 1.7 × 10²⁰ cm^-2, the X-ray count rate peaks at 65 s⁻¹, much higher than was observed for XRO 080109. Increasing the column depth brings down the peak count rate, narrows the width of the light curve, and hardens the spectrum. For N_H = 2 × 10²¹ cm^-2, the peak count rate is 10 s⁻¹, slightly greater than the observed value. The light curve at this column depth is somewhat narrower and rises more quickly than that of XRO 080109. The spectrum is also still too soft and devoid of significant count rates above 2 keV. In Figure 26, we show the angular dependence of the simulated light curve and spectra for the enhanced temperature version of m2r1cold calculated with a metallicity of 0.5 Z_☉ and N_H = 1.7 × 10²⁰ cm^-2. The spectrum is harder than the fiducial models, computed without temperature enhancement, but is still softer than that of XRO 080109.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The color temperature of the emission may also be enhanced by non-LTE effects. The radiation escaping from the shock during breakout from a Wolf–Rayet star is not in thermal equilibrium with the gas at τ_* = 2/3 (Katz et al. 2010; Nakar & Sari 2010). In fact, during the early parts of the shock breakout the radiation temperature may be orders of magnitude greater than the gas temperature at τ_* = 2/3 (Nakar & Sari 2010). Over the course of the first tens of seconds of the X-ray burst, the radiation temperature will quickly drop down closer to the gas temperature. This time-dependent enhancement of the radiation temperature, if accounted for in our calculations, could produce the hard X-rays that are missing in our spectra while not dramatically increasing the luminosity at later times in the light curve.

4.3. SED Shape

The shapes of the spectral energy distributions of our simulated shock breakout models are not well fit by single-temperature and radius spherical blackbodies. This is the case for both the jet-driven explosion simulations and the spherical explosion models. Figures 27 and 28 show the X-ray spectra of models m2r1cold and m2r1sph along with spherical single-temperature and radius blackbody SEDs corrected for XRT detector response and X-ray absorption. As is shown, no single-temperature and radius blackbody can fit both the soft and hard parts of the spectrum in either the jet-driven or spherical explosions. This demonstrates that, regardless of shock break out geometry, such simple blackbody models are not applicable. Shock breakout is too dynamic a process to be modeled by a single average radius and a single average temperature blackbody, even for spherical breakouts. Simple dynamical models with radii expanding at constant velocities and temperatures cooling adiabatically could provide vast improvements over simple static models and could provide a better relation to the physical process of shock breakout.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have modeled the emission arising from the breakout of an aspherical SN shock assuming that the emission is thermal in nature. Our hydrodynamic simulations do not account for radiative effects such as disparate ion and radiation temperatures, radiative pre-acceleration of gas ahead of the shock, escape of photons from the shock at optical depths greater than 1, etc. Despite the limitations of our simulations, we have demonstrated that in the more general case of a non-spherical explosion, the breakout emission is dramatically different from that expected from a spherically symmetric shock breakout. One of the most important results of this work is that the timescales of the light curves for aspherical shock breakout are not set by the light travel time across the progenitor star, but are instead related to the much longer shock crossing time of the progenitor. Thus, shock breakout light curves contain information about the geometry of the shock structure as well as the radius of breakout. We also show that for aspherical shock breakouts, the observer's viewing angle can play an important role in determining the shape of the observed light curve. These general results apply to any, arbitrarily aspherical SN, not just a jet-driven, bipolar SN.

The fiducial emission models are generally under-luminous and have soft spectra compared to the actual observations. Tanaka et al. (2009) advocate explosions with higher energies than what we have simulated and Soderberg et al. (2008) posit that the wind around the progenitor of SN 2008D was optically thick, thus very dense. Either an increased explosion energy or increased wind density could increase the X-ray luminosity of shock breakout; however, even if the wind and explosion energies were increased to give luminosities commensurate with those inferred for XRO 080109, the spectrum would still remain lacking in sufficient luminosity above 2 keV to match the observed spectrum (see Figure 25). This seems to indicate the need for a scattering of the thermal X-ray photons to higher energies, as prescribed by Soderberg et al. (2008) and Wang et al. (2008), or other, non-thermal processes (Katz et al. 2010; Nakar & Sari 2010). This could serve to harden our simulated spectra and also lengthen the burst time since, in our models, the later parts of the light curve are made up of a larger fraction of soft, highly absorbed photons.

The simulated spectra may also be hardened by relaxing the assumption of thermal equilibrium between the matter and radiation. During shock breakout from a massive, compact Wolf–Rayet progenitor the matter and radiation are not in LTE (Katz et al. 2010; Nakar & Sari 2010). Our method does not account for non-LTE effects. As discussed by Nakar & Sari (2010), if the matter and radiation are not in LTE, the observed color temperature could be much greater than the matter temperature at τ_* ≈ 2/3. Because the internal energy would not be similarly enhanced above the LTE case, the bolometric luminosities for non-LTE breakouts will not depend strongly on the coupling between matter and radiation, as the color temperature does. The color-temperature enhancing effects in non-LTE breakouts discussed by Nakar & Sari (2010) are strongly time dependent with the greatest difference in the temperatures of radiation and matter being at the instant of shock breakout. If non-LTE effects, such as described by Nakar & Sari (2010), were included in our simulations, the spectra would be significantly hardened, especially at early times in the breakout, and the X-ray luminosities modestly enhanced. This could account for the missing hard X-ray emission in our simulated spectra and increase the XRT count rates without the need to invoke a very dense wind or a greatly enhanced explosion energy. We note, however, that Nakar & Sari (2010) do not include the effects of a wind surrounding the progenitor. In our calculations, the presence and character of the wind is an integral factor in the breakout emission formation since the thermalization depths in our simulations lie at, or ahead of, the contact discontinuity between the ejecta and the wind.

The presence of a double-peaked oxygen line in the spectra of SN 2008D at about 109 days may be evidence for asphericity in the explosion (Modjaz et al. 2009). Modjaz et al. (2009) suggest that this may be evidence for a ring-like distribution of oxygen, similar to the conclusions reached by Modjaz et al. (2008) for several other Type Ib/c SNe with double-peaked oxygen lines.⁸ A double-peaked oxygen line may also result from an asymmetric distribution of radioactive cobalt powering the excitation of oxygen (Modjaz et al. 2008; Gerardy et al. 2000). If the cobalt distribution were aspherical, however, other radioactively excited lines would be expected to show a similar double-peaked behavior, which is not the case for SN 2008D or the Type Ib/c SNe discussed in Modjaz et al. (2008). Tanaka et al. (2009) find that a fraction of nickel must be mixed outward to adequately model the spectra of SN 2008D, further indicating asphericity of the explosion. Figure 29 shows the nickel and oxygen mass fractions from our model m2r1cold at 10⁵ s. In our simulations, both the nickel and oxygen distributions are aspherical, as is shown in Figure 29. This may be able to account for the double-peaked oxygen line in the spectra of SN 2008D, however detailed spectral synthesis calculations are needed to be sure. Note that our oxygen distribution is prolate, not ring-like as Modjaz et al. (2008) recommend for explaining double-peaked oxygen lines in SNe.

Zoom In Zoom Out Reset image size

The distribution of intermediate-mass elements such as oxygen we find in our simulations differs from what some other groups find in similar studies of aspherical SNe. Maeda et al. (2002) present hydrodynamic simulations of aspherical CCSNe targeted to explaining the observations of SN 1998bw. They show that the intermediate elements are ejected from the explosion in an equatorial torus. We find that the intermediate elements are ejecta in a bipolar geometry, as shown in Figure 29. The difference in the two results comes from the manner in which the explosions are initiated. In Maeda et al., the explosions are started by depositing kinetic energy asymmetrically in the center of the progenitor. This pushes the intermediate elements outward while simultaneously the more energetic material near the poles compresses them into a toroidal geometry. In our simulation the explosions are driven entirely by the bipolar jets. The intermediate-mass material near the progenitor's equator is accreted into the central engine. Intermediate elements are entrained in the jets and carried out into a configuration that resembles that of the jets themselves. The equatorial torus in our simulations is comprised primarily of helium. The final distribution of intermediate-mass elements thus depends on the mode of asymmetric energy input. Determining the distribution observationally may thus help to constrain models.

Maund et al. (2009) present early spectropolarimetric observations of SN 2008D. They find that the continuum polarization is relatively small, indicating that the SN photosphere may be only slightly aspherical. They also find that there is significant polarization in certain spectral lines, indicating that the line-forming regions of various elements are markedly aspherical. This is in qualitative agreement with the results of our simulations. The photospheres at late times are nearly round, while the detailed composition structure is dramatically asymmetric, as shown in Figure 29. A late-time shock structure that is nearly round, or at least not dramatically aspherical, is also consistent with the radio measurements of Bietenholz et al. (2009).

The number of observed SN shock breakouts has been increasing, and this trend is likely to continue, opening a new window for exploring core-collapse SN through observation. A multi-dimensional interpretation of these observations will be critically important for gaining an accurate understanding of SN shock breakout. We find that aspherical shock breakout can account for the light curve timescales observed for XRO 080109/SN 2008D without requiring an extremely dense wind or an abnormally large Wolf–Rayet progenitor star. Our models with the smaller progenitor, m2r1cold and m2r1hot, have light curves with FWHMs of about 100 s, roughly matching XRO 080109. We find that our models are generally under-luminous and have spectra too soft to match XRO 080109; however, this is likely due to our assumption of LTE. The shock breakout from a Wolf–Rayet progenitor is not in LTE and the radiation temperature in the non-LTE case could be significantly enhanced above the matter temperature (Nakar & Sari 2010). Accounting for non-LTE effects, then, could brighten our simulated light curves and harden the simulated spectra. Inclusion of non-LTE effects will be presented in future work.

The authors thank Peter Höflich and Ehud Nakar for helpful discussion concerning the physics of shock breakout. S.M.C. is very grateful for support provided through the NASA Earth and Space Science Fellowship program. J.C.W. is supported in part by NSF grant AST-0707769 and in part by DOE funding to the Center for Laser Experimental Research (CLEAR) through the University of Michigan. The software used in this work was in part developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high-performance computing, visualization, and data storage resources that have contributed to the research results reported within this paper.