HEAD-ON COLLISIONS OF WHITE DWARFS IN TRIPLE SYSTEMS COULD EXPLAIN TYPE Ia SUPERNOVAE

There is strong evidence that Type Ia supernovae (SNe Ia) are thermonuclear explosions of carbon–oxygen white dwarfs (CO-WDs; Hoyle & Fowler 1960). It is widely thought that the explosion is caused by accretion of matter onto the CO-WD, and as it approaches the unstable Chandrasekhar mass limit of ≈1.4 M_☉, a poorly understood mechanism causes the required explosion (see Hillebrandt & Niemeyer 2000, for a review). In these models, the uncertainties in the burning waves allow the introduction of free parameters (such as the deflagration velocity or transition to detonation criteria), tuned to account for the observations. It was recently shown that the rate of head-on collisions of typical field WDs in triple systems may be as high as that of SNe Ia (Katz & Dong 2012; see, however, Hamers et al. 2013), and it was suggested that some or all SNe Ia are due to such collisions. Here we provide evidence supporting a scenario in which the majority of SNe Ia are the result of such head-on collisions of CO-WDs. In this case the nuclear detonation is due to a well understood shock ignition, devoid of the commonly introduced free parameters, and unrelated to the Chandrasekhar limit.

Collisions of WDs were previously suggested to occur mainly in dense stellar systems (e.g., Hut & Inagaki 1985; Sigurdsson & Phinney 1993; Thompson 2011).⁴ While such collisions are believed to have rates which are orders of magnitude smaller than the rate of SNe Ia, they motivated three-dimensional (3D) hydrodynamic simulations of the thermonuclear explosion of colliding WDs (Benz et al. 1989; Raskin et al. 2009, 2010; Rosswog et al. 2009; Lorén-Aguilar et al. 2010; Hawley et al. 2012). While the amount of ⁵⁶Ni (the decay of which powers the observed light; Colgate & McKee 1969) synthesized in most of these simulations was non-negligible, some results were contradictory, with inconsistent amounts of ⁵⁶Ni and different ignition sites for the same initial conditions.

Section 2 describes our simulations of zero-impact-parameter CO-WD collisions, which fully resolve the ignition process. In Section 3, we compare our results to previous works, resolving previous discrepancies. Section 4 discusses several observational tests, which avoid radiation transfer uncertainties, and which demonstrate the consistency of our models with the majority of SNe Ia. Section 5 presents avenues for future research.

We consider the collisions of (equal and non-equal mass) CO-WDs with masses 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 M_☉, covering the range of CO-WD masses. This problem is axisymmetric, allowing the use of two-dimensional (2D) numerical simulations with high resolution (∼few km cell size), which is higher by at least an order of magnitude than those in previous Eulerian approaches (limited to >100 km cell sizes; Hawley et al. 2012) and smooth particle hydrodynamic (SPH) models (Raskin et al. 2010, in which the resolution is estimated from the "cell volume," which corresponds to particle mass divided by the density). We employ two different hydrodynamic codes: FLASH4.0 (Dubey et al. 2009; Eulerian, adaptive mesh refinement) and VULCAN2D (Livne 1993; arbitrary Eulerian–Lagrangian, ALE). FLASH4.0 solves the equations of reactive hydrodynamics by combining the dimensionally-split Piecewise Parabolic Method (Colella & Woodward 1984; Fryxell et al. 2000) with a 19 isotope alpha-chain reaction network (Timmes 1999). The system of equations is closed with the Helmholtz equation of state (Timmes & Swesty 2000) and a multipole gravity solver. VULCAN2D solves the reactive hydrodynamic equations on an ALE grids using a Lagrangian step followed by a remapping step on a new grid with a 13 isotope alpha-chain reaction network. We exploit the flexibility of the grid to refine it near the collision interface, so that at the onset of a detonation a maximal, pre-determined resolution is reached there. The equation of state, gravity solver and boundary conditions are similar to those used in FLASH4.0. Unless otherwise stated, we present the FLASH4.0 results.

Initially the CO-WDs are at contact with free fall velocities. The structure of each CO-WD is obtained from an isothermal stellar model⁵ at T = 10⁷ K and with a uniform composition of 50% carbon and 50% oxygen by mass.

Figure 1 shows an example of the ignition process and the ejecta formation for a 0.64–0.64 M_☉ collision. For this case we verified that placing the CO-WDs with a separation of four stellar radii between their centers has a minor effect (a few percent) on our results. We used VULCAN2D to evolve several cases, and found the same ⁵⁶Ni mass and the same location of the detonation ignition (see, e.g., the right panel of Figure 2). In both codes, collisions involving CO-WDs with mass ≲ 0.7 M_☉ ignite at the shock region (as in Figure 1), while higher mass CO-WDs ignite at the contact region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The shocked region at the vicinity of the symmetry axis has approximate planar symmetry (see, e.g., panel (c) of Figure 1). This allowed us to develop a one-dimensional (1D) planar model (evolved with Lagrangian and Eulerian schemes) to ensure that 2D simulations correctly resolve the ignition process (panel (a) of Figure 2). The initial CO-WD density of the 1D model equals the density on the axis of symmetry in the 2D model, and the initial velocity is the free fall velocity. The gravitational field is mimicked by an adjustable acceleration, which is constant in time and space. For example, the acceleration g₀ ≈ 1.1 × 10⁸ cm s⁻², approximately reproduces the 2D velocity profiles of the 0.64–0.64 M_☉ collision in Figure 1, and results in a similar location of the ignition (∼1.5 × 10⁸ cm from the contact surface). In particular, this suggests that the curvature of the shock has little effect on the ignition process. By artificially changing the acceleration we verify the robustness of the correct ignition location. Only for a substantially higher value of the acceleration, g₀ ≈ 1.6 × 10⁸ cm s⁻², does the detonation occur at the contact region.

Numerically unstable burning occurs if the energy in a cell is significantly increased in a time shorter than the sound crossing time, t_s ≡ Δx/c_s, where Δx is the length scale of the cell and c_s is the sound speed. This is a severe problem for Eulerian calculations with cell sizes larger than ∼km. Stability is achieved by limiting the energy injection rate from burning to fε/t_s, where ε is the internal energy of the cell and f ≪ 1. This is implemented in all of our simulations by appropriate renormalization of all burning rates within a cell if the limit is exceeded. As the numerical resolution increases, the renormalization becomes less severe, guaranteeing correct convergence while avoiding premature ignition. The limiter does not modify the ignition of a detonation wave process, which is set at lower temperatures, where the stability criterion is automatically satisfied. By comparing the 1D Eulerian model with resolution Δx ∼ few km to a converged 1D Lagrangian model, we verify that the limiter does not affect the ignition and burning for f ≲ 0.4. Therefore, we adopted f = 0.1 in our simulations. In several cases we find that if the limiter is not used, a premature numerical ignition may occur (see, e.g., left panel of Figure 2).

The amount of ⁵⁶Ni synthesized in the various collisions is shown in Figure 3,⁶ which is the main result of this Letter. We find that CO-WD collisions lead to the synthesis of ∼0.1–1 M_☉ of ⁵⁶Ni, covering the range of yields observed in the vast majority of SNe Ia (see, e.g., Figure 5). Furthermore, there is a tight correlation between the ⁵⁶Ni yields and the total mass of the colliding CO-WD masses, which is insensitive to their mass ratio.

Zoom In Zoom Out Reset image size

For equal mass collisions, similar ⁵⁶Ni yields were obtained by most of the SPH simulations of Rosswog et al. (2009) and of Raskin et al. (2010) (except for the 0.5–0.5 M_☉ collision, for which both groups obtained a negligible amount of ⁵⁶Ni), and significantly lower yields (by factor of ∼2) were obtained by FLASH calculations of the same groups (Rosswog et al. 2009; Hawley et al. 2012). However, most non-equal mass SPH results contradict ours. Four such simulations exist, and three of them (the 0.9–0.6 M_☉ collision of Rosswog et al. 2009 and the 0.81–0.64 M_☉ and 1.06–0.64 M_☉ collisions of Raskin et al. 2010) have lower ⁵⁶Ni yields than ours by factors of a few up to orders of magnitude. As a result, the tight correlation between the ⁵⁶Ni yields and the total mass of the colliding CO-WD masses was not seen in any of the previous works.

We reproduce the low ⁵⁶Ni yields of the discrepant simulations by running simulations with similarly low resolutions. We are able to identify the main cause for the discrepancy in most cases. In a 0.64–0.64 M_☉ simulation with a similar setup as in Hawley et al. (2012), a premature ignition occurs at the contact region that prevents the correct shock ignition at later times. This occurs mainly due to the fact that Hawley et al. (2012) did not implement the required burning limiter to avoid numerical unstable burning (see Section 2). We also performed 0.5–0.5 M_☉, 0.81–0.64 M_☉, and 0.9–0.6 M_☉ simulations with similar resolutions to the cell volumes of the SPH simulations. We find that the propagation of the detonation wave is significantly altered due to the low resolution. For example, in our high-resolution 0.81–0.64 M_☉ simulation, following ignition in the 0.64 M_☉ star, the detonation wave propagates inside the 0.64 M_☉ star in the r-direction up to a height of ∼3000 km from the symmetry axis, and then crosses into the 0.81 M_☉ star. However, the detonation crossing, which is responsible for the large amount of ⁵⁶Ni (∼0.4 M_☉), is not resolved by Raskin et al. (2010), in which the resolution is lower than our high-resolution runs by a factor of ∼100 in the crossing location (density of few × 10⁶ g cm⁻³). As a consequence, the low-resolution runs obtain small amounts of ⁵⁶Ni. Given that similar dynamics are observed in our 1.0–0.6 M_☉ collision, we suspect that the low ⁵⁶Ni yield obtained for the 1.06–0.64 M_☉ collision of Raskin et al. (2010) is due to the same reason.

The amount of ⁵⁶Ni obtained in our simulations (Figure 3), implies that collisions of typical CO-WDs in the mass range 0.5–1 M_☉ produce ⁵⁶Ni masses in the range $0.1-1\,M_{\odot }$ in agreement with the inferred range of SNe Ia (see, e.g., Figure 5). The tight correlation between the ⁵⁶Ni yields and the total mass of the colliding CO-WD masses, which is insensitive to their mass ratio, may explain the fact that many of the SNe Ia properties are tightly correlated.

One constraint for explosion models is the observed substantial amount of high-velocity intermediate mass elements in the outer layers of the ejecta. Such elements are naturally produced in the collisions (e.g., Rosswog et al. 2009; Raskin et al. 2010), as can be seen in Figure 1.

Our models pass two additional quantitative, nontrivial, and robust observational tests, which are independent of the complicated optical radiation transfer.

The main channel by which the radioactive decay chain ⁵⁶Ni → ⁵⁶Co → ⁵⁶Fe converts its energy to observed light, is through the energy deposition of γ-ray photons and positrons into the ejecta. The heated ejecta reprocesses this energy into observed optical light. The late-time (≳ 60 days since explosion) bolometric luminosity equals the instantaneous γ-ray energy deposition rate, which is calculated by a Monte Carlo code for the Compton-scattering-dominated transport of γ-rays (with a small contribution from positrons which are assumed to deposit their kinetic energy locally). The ejecta is taken from the calculations at a sufficiently late time where the expansion is homologous. We find that the deposited fraction of γ-rays as a function of time is the same in all of our models (within 30%, where our results converge to better than the numerical statistical noise of ∼3% ⁷) and in excellent agreement with the sample of 24 low extinction (Galactic+host E(B − V) < 0.3) bolometric light curves of Stritzinger et al. (2005), when appropriately normalized by the (time-weighted) integrated luminosity (see Figure 4). Given that the ejecta is not spherically symmetric, the emission may depend on viewing angle, possibly widening the scatter of expected late-time luminosities. While this possible correction is beyond the scope of this Letter, the angle-average value should be accurately represented by our results.

Zoom In Zoom Out Reset image size

To illustrate that this agreement is nontrivial, we calculated the light curves of simple Chandrasekhar toy models (see, e.g., Jeffery et al. 1992; Nugent et al. 1995; Woosley et al. 2007), with an exponential density distribution as function of velocity ($\rho \propto e^{-v/v_{e}}$, where v_e is set by the total release of energy in the explosion) and varying ⁵⁶Ni masses (0.15–0.8 M_☉). While the late-time light curve of the 0.8 M_☉ ⁵⁶Ni Chandrasekhar-model lies at the high end of observed light curves, the 0.15 M_☉ ⁵⁶Ni Chandrasekhar-model has a much higher normalized light curve and is inconsistent with observations (see Figure 4). The strong emission at late time is due to the efficient deposition of γ-rays in the slowly moving, massive 1.25 M_☉ ejecta surrounding the ⁵⁶Ni. To confirm that the main problem with the Chandrasekhar toy model is the high ejecta mass, we apply a "collision-like" spherical model. This model has the same exponential density distribution, but with a smaller ejecta mass (1.0 M_☉), roughly corresponding to the relevant head-on collision for the 0.15 M_☉ ⁵⁶Ni yield (the kinetic energy is set by the same collision). The "collision-like" model has a lower normalized light curve, and shows a much better agreement with the observations (Figure 4). We note that Chandrasekhar models in which the ⁵⁶Ni mass is located far from the center may alleviate this particular problem. More details will be provided in a subsequent paper (S. Dong et al. in preparation).

The distribution of widths of the (⁵⁶Ni mass-weighted) line-of-sight (LOS) velocities for the collision models is compared to those inferred from nebular-phase observations of SNe Ia in Figure 5. For each viewing angle, the model widths (v_mod) are obtained by fitting the distribution of ⁵⁶Ni mass per unit LOS velocity in that angle using a quadratic function, $dM_{^{56}{\rm Ni}}/dv_{\rm LOS}\propto 1-v_{{\rm LOS}}^{2}/v_{{\rm mod}}^{2}$ (e.g., Mazzali et al. 1998). The observational widths are obtained by fitting the observed late-time (>150 days) spectra⁸ in the range 4800–5700 Å to the convolution of the low-width spectra templates of 1991bg and 1999by with the same quadratic velocity distribution. We note that while several properties of 1991bg and 1999by are uncommon, their nebular spectra in the above range are very similar to other SNe (e.g., Mazzali et al. 1998, except for the line widths) and should suffice as templates for measuring the line widths of other SNe (see B. Katz et al. in preparation, for more details). The amount of ⁵⁶Ni is obtained by fitting the bolometric light curves at t ∼ 60 days to the universal injection function presented in Figure 4, which is well described by L_deposit = (1 + (t_d/40)³)^−2/3L_decay, where t_d is the time since explosion in days and L_decay is the energy injection in γ-rays by ⁵⁶Ni and ⁵⁶Co. There is a clear correlation between the observed ⁵⁶Ni masses and the nebular-phase velocity widths. Both the correlation and the scatter of this Mazzali relation (Mazzali et al. 1998) are well reproduced by the collision model (10 viewing angles for each calculation, equally spaced in cos (i); these widths are converged to a level of ∼20%). As can be seen in Figure 5, the simple Chandrasekhar-model predicts a similar correlation, somewhat in offset. The amount of SNe Ia with nebular spectra and well described bolometric light curves is limited. In the top panels larger samples⁹ are used to show the continuous correlation of observational features of SNe Ia with ⁵⁶Ni yields. The strong correlation between ⁵⁶Ni yield and the total mass of the colliding CO-WDs, Figure 3, is suggestive as the source for these correlations, possibly including the Philips relation (Phillips 1993).

Zoom In Zoom Out Reset image size

A future robust and detailed test is the time dependent MeV scale γ-ray spectrum, which can be computed in a straightforward manner. Other tests, closely related to the nebular velocity widths, include the distribution of non-zero nebular line shifts, which are expected in the collision of non-equal mass CO-WDs, as well as the unique line shapes resulting from the nontrivial velocity distribution in many collision scenarios. Finally, the early-time spectra and light curves can be compared to observations using suitable radiation transfer models (see encouraging attempts in Rosswog et al. 2009).

We stress that our calculations were restricted to the zero-impact-parameter. Only if the effect of the non-zero-impact-parameter is small for a large range of impact parameters, do our calculations represent a significant fraction of collisions (the distribution of impact parameters is uniform for collisions in triple systems; Katz & Dong 2012, e.g., 50% of collisions have impact parameters <0.5(R₁ + R₂), where R_{1, 2} are the radii of the CO-WDs). Raskin et al. (2010) showed that ignition of a detonation can be triggered for several non-zero-impact-parameter collisions in SPH simulations. However, their simulations may suffer from numerical problems which are similar to the ones we identified for their zero-impact-parameter runs in Section 3, making these results quantitatively uncertain. Non-zero-impact-parameter collisions should be studied using 3D hydrodynamical codes of similar accuracy to the 2D calculation presented here. High-resolution 3D calculations are also desirable for constraining 3D effects for the zero-impact-parameter collisions.

We thank A. Gould, R. Kirshner, J. Prieto, M. Zaldarriaga, E. Waxman, and F. Dyson for discussions. D.K., S.D., and R.F. are supported by NSF grant AST-0807444. B.K. is a John N. Bahcall Fellow, and supported by NASA through the Einstein Postdoctoral Fellowship awarded by Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. S.D. was supported through a Ralph E. and Doris M. Hansmann Membership at the IAS. FLASH was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. Computations were performed at PICSciE and IAS clusters.