ZERO IMPACT PARAMETER WHITE DWARF COLLISIONS IN FLASH

W. P. Hawley; T. Athanassiadou; F. X. Timmes

doi:10.1088/0004-637X/759/1/39

1. INTRODUCTION

Supernovae Type Ia (SN Ia) have continued to be foremost probes of the universe's accelerating expansion (Riess et al. 1998; Perlmutter et al. 1999; Riess et al. 2011; Sullivan et al. 2011; Suzuki et al. 2012). While light curves between different SNe Ia vary, the variations generally correlate with distance independent light curve properties, such as the decline from B-band maximum after 15 days (Phillips 1993). Calibration of the light curves onto a standard template yields distance indicators accurate to ∼10% (e.g., Silverman et al. 2012) and are primarily applied to SNe Ia not showing peculiarities (Branch et al. 1993). These "normal" SNe Ia presumably emerge from a homogeneous population of white dwarf (WD) progenitors. While the favored population is thought to be a carbon–oxygen WD accreting matter from a non-degenerate companion star (e.g., Whelan & Iben 1973; Thielemann et al. 1986), recent observations suggest that a fraction of SNe Ia may derive from double-degenerate progenitors (Howell et al. 2006; Hicken et al. 2007; Gilfanov & Bogdán 2010; Bianco et al. 2011).

In view of these and other observations of SNe Ia progenitor systems, recent theoretical studies have explored double-degenerate mergers and collisions of WDs as potential progenitors of SNe Ia (Guerrero et al. 2004; Yoon et al. 2007; Maoz 2008; Lorén-Aguilar et al. 2009, 2010; Raskin et al. 2009, 2010; Rosswog et al. 2009; Pakmor et al. 2010, 2012; Shen et al. 2012). Almost all of these efforts use smooth particle hydrodynamic (SPH) codes to model most of the collision or merger process. SPH and Eulerian grid codes, such as FLASH (Fryxell et al. 2000), have well-known complementary strengths and weaknesses—particle codes are inherently better at angular momentum conservation, whereas grid codes have a superior treatment of shocks. Only Rosswog et al. (2009) included a zero impact parameter WD collision model with FLASH. They used a mirror gravitational potential for one WD at one spatial resolution. They found their FLASH calculations yielded about half as much ⁵⁶Ni as the equivalent SPH calculation (0.32 M_☉ for SPH, 0.16 M_☉ for FLASH).

In this paper, we use the Eulerian adaptive-mesh refinement code FLASH to model the zero impact collisions between 0.64 + 0.64 M_☉ and 0.81 + 0.81 M_☉ carbon–oxygen WD mass pairings. Like the single case studied by Rosswog et al. (2009), our configurations are highly idealized cases of head-on collisions between identical, initially spherical WDs. One aim of our paper is to determine whether or not, given presently available computing resources and numerical algorithms, simulations of collisions can be used to reliably predict the fraction of WD material that is converted by explosive nucleosynthesis into heavier elements such as silicon and nickel. Other efforts have focused on the realism of the initial conditions and subsequent evolution, including, but limited to, in-spiraling from a binary orbit (Rasio & Shapiro 1995; Pakmor et al. 2010; Dan et al. 2011; Raskin et al. 2012), unequal mass collisions (Benz et al. 1989, 1990; Rosswog et al. 2009; Lorén-Aguilar et al. 2010; Raskin et al. 2010; Pakmor et al. 2012), and the final long-term fate of merged systems (van Kerkwijk et al. 2010; Yoon et al. 2007; Shen et al. 2012). Our simulations, through their idealized nature, highlight the essential physics and numerical convergence properties of the simplest possible configuration. In addition, our idealized configurations form a baseline for further studies that incorporate more realistic initial conditions.

Our paper is organized as follows. In Section 2, we describe the input physics, initial conditions, and boundary conditions of our FLASH simulations. In Section 3, we discuss the results of our studies over a range of spatial resolutions and time step choices, and in Section 4 we explore the implications of our results and describe directions for future studies.

2. INPUT PHYSICS, INITIAL CONDITIONS, AND BOUNDARY CONDITIONS

Our three-dimensional simulations are carried out with FLASH 3.2, a three-dimensional Eulerian hydrodynamics code with adaptive-mesh refinement (Fryxell et al. 2000; Calder et al. 2002). We use the included Helmholtz equation of state (Timmes & Swesty 2000), the 13 isotope alpha-chain reaction network that includes isotopes from ⁴He to ⁵⁶Ni to model energy generation from nuclear burning (Timmes 1999), and the multigrid Poisson gravity solver with Dirichlet boundaries (Ricker 2008). All the simulation domain boundaries use a diode boundary condition, which is a zero-gradient boundary condition where fluid velocities are not allowed to point back into the domain. We follow both WDs in three-dimensional rectilinear coordinates throughout calculation, rather than using a mirror gravitational potential and evolving one WD (as used in Rosswog et al. 2009).

Our initial one-dimensional WD profiles are calculated assuming hydrostatic equilibrium and mass conservation.³ Our initial WD models use the same equation of state as in FLASH, namely, the Helmholtz equation of state. We assume a uniform composition of 50% ¹²C and 50% ¹⁶O. We use WDs with masses 0.64 and 0.81 M_☉, to match the masses used in Raskin et al. (2010), with an isothermal temperature of 10⁷ K. We map the one-dimensional WD profiles for the density, temperature, and composition onto a three-dimensional rectilinear Cartesian grid. Our computational domain is a cubic box chosen to be eight times the WD radius (see Figure 1). The WDs are initially placed four WD radii apart from center to center, which is large enough to allow the subsequent evolution to produce tidal distortions while allowing sufficient numerical resolution in the central regions.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Two-dimensional slice of density through the x–y mid-plane at t = 0.0 s for the 0.64 + 0.64 M_☉ collision. Each tick mark has a value of one white dwarf radius, which is 8.3 × 10⁸ cm. The size of the domain is equal to eight white dwarf radii, and the white dwarfs are positioned four white dwarf radii apart from center to center.
Download figure:
Standard image High-resolution image

The symmetry of head-on collisions between identical, initially spherical WDs suggests two-dimensional axisymmetric simulations may have been sufficient. Our rationale for deploying three-dimensional rectilinear coordinates is threefold. First, we want to explicitly show that FLASH maintains symmetry throughout the collision and subsequent explosion processes. Second, we want to compare our results on this important numerical test case with other existing three-dimensional calculations (both grid and particle). Imposing axisymmetric conditions would have complicated these comparisons because we would not know if differences from existing three-dimensional models were driven by different physics, different numerics, or the imposition of axisymmetry itself. Third, we anticipate exploring unequal mass and non-zero impact parameter collision models with FLASH, both of which violate axisymmetry. To better assess the impact of these effects requires an equal mass, zero impact parameter, three-dimensional benchmark calculation.

We use the free-fall expression for the initial, relative speed of the two WDs, v = [2G(M₁ + M₂)/Δr]^1/2, where M_i are the masses of the constituent WDs and Δr is the initial separation of their centers of mass, which for our initial conditions is 4R_WD. Each WD moves toward the other WD, one in the positive x-direction and the other in the negative x-direction, with half of the relative speed. The centers of both stars lie on the x-axis, and thus the initial velocities are purely in the x-direction.

The surrounding ambient medium is set to the same temperature as the isothermal WDs with a density that is small (10⁻⁴ g cm⁻³) compared to the density of the outermost regions of the WD (∼1–10 g cm⁻³). Table 1 lists the initial conditions for each of our six simulations.

Table 1. Initial Conditions for the Three-dimensional FLASH Models

No.	M₁, M₂	l	R	D	v₁, v₂	f	R_WD	ρ_WD
	(M_☉)		(10⁷ cm)	(10⁹ cm)	(10⁸ cm s⁻¹)		(10⁸ cm)	(10⁶ g cm⁻³)
1	0.64	5	5.19	6.64	±1.59	0.2	8.30	4.51
2	0.64	6	2.59	6.64	±1.59	0.2	8.30	4.51
3	0.64	7	1.30	6.64	±1.59	0.2	8.30	4.51
4	0.81	5	4.32	5.51	±1.97	0.2	6.88	11.2
5	0.81	6	2.16	5.51	±1.97	0.3	6.88	11.2
6	0.81	7	1.08	5.51	±1.97	0.3	6.88	11.2

Note. Columns are the run number, white dwarf masses (M₁, M₂), maximum level of refinement (l), maximum spatial resolution (R), domain size (D), white dwarf initial velocities (v₁, v₂), the value of the time step limiter (f), white dwarf radii (R_WD), and central white dwarf densities (ρ_WD).

Download table as: ASCII Typeset image

Our FLASH models begin with one top-level initial block, where each block contains eight cells in each direction (x, y, z). The blocks are refined or derefined at each time step based on changes in density and pressure. For each successive level of refinement, the block size decreases by a factor of two, creating a nested block structure. At maximum refinement, the smallest block size is determined by R = D/(8 × 2^{l − 1}), where D is the domain size in one dimension and l is the maximum level of refinement as seen in Table 1. At first contact between the WDs, shock waves are sent into the ambient medium, causing the grid in the ambient medium to rapidly become maximally refined. To avoid concentrating resources on these less interesting regions of the models, we use a derefine procedure at first contact that sets a radius equal to 1.2 WD radii beyond which the blocks in the ambient material are forced to be less refined than the blocks in the collision region.

The nuclear reaction network in FLASH uses constant thermodynamic conditions over the course of a time step. However, the Courant-limited hydrodynamic time step may be so large compared to the burning timescale that the nuclear energy released in a cell may exceed the existing specific internal energy. To ensure the hydrodynamics and burning remain coupled, as well as to capture the strong temperature dependence of the nuclear reaction rates, we limit the time step as a result of nuclear burning by a factor f, which constrains the maximum allowable change in specific internal energy. The overall time step is dt_{n + 1} = min [dt_hydro, dt_burn], where

$\begin{equation} dt_{\mathrm{burn}} = dt_{n} \times f \times \frac{u^{i} _{n-1}}{u^{i} _{n}-u^{i} _{n-1}} {,} \end{equation} \tag{ 1 }$

where the subscript n refers to the time step number, dt_hydro is the hydrodynamic time step, dt_burn is the burning time step, and uⁱ is the specific internal energy of the ith cell. Table 1 lists the nominal values of f used for our six simulations, and the effects of using different values of f are discussed in Section 3.2.

3. RESULTS

3.1. General Features of the Collision Models

Zero impact parameter, or head-on, WD collisions undergo four distinct phases of evolution. First, the WDs become tidally distorted as they approach each other. For the 0.64 + 0.64 M_☉ case (hereafter 2 × 0.64), the velocity gradient across the WD at first contact ranges from about 3500 km s⁻¹ to 5000 km s⁻¹. Second, a shock wave is produced normal to the x-axis at first contact. The shock stalls because the speed of infalling material and the sound speed are comparable. Third, nuclear burning is initiated within the stalled shock region. Finally, the nuclear energy released unbinds the system, leading eventually to homologous expansion.

An overview of the evolution of the 2 × 0.64 collision is shown in Figure 2. The three-dimensional calculation has been sliced through the x–z mid-plane to show detail at the center of the collision. Due to the symmetry of a head-on collision, a cut through the x–z mid-plane will look identical to a cut through the x–y mid-plane. The top panel of the figure represents four times in the collision from the start of the simulation (t = 0.0 s), to first contact (t = 4.0 s), to the formation of the stalled shock region (t = 5.0 s), and finally to the jettisoning of material orthogonal to the x-axis (t = 5.5 s).

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Three-dimensional images cut through the center of the y-axis of the 2 × 0.64 collision density from first contact at 0.0 s to after ignition at 6.9 s. For scale, each white dwarf has a radius of 8.3 × 10⁸ cm. The density color bar is logarithmic and extends from 10⁻⁴ to 10¹⁰ g cm⁻³.
Download figure:
Standard image High-resolution image

Given the WD radius and initial velocity shown in Table 1 for the 2 × 0.64 collision, the time to first contact would be 2R/v = 5.2 s if the initial speed was constant and the WDs remained spherical. However, the initial speed increases due to gravitational acceleration and tidal distortion causes the WDs to become elongated along the x-axis. As a result the two WDs experience first contact sooner, at about 4.0 s.

The bottom panel represents the further progression of the collision from the continued jettison of material (t = 6.0 s), to just before ignition (6.5 s), to just after ignition (t = 6.8 s), and finally to the spread of nuclear burning through the WDs (t = 6.9 s). These steps are discussed in further detail below.

Figure 3 shows the thermodynamic, mechanical, and morphological properties of the 2 × 0.64 head-on collision model. At 6.00 s after the beginning of the model, the WDs are past first contact but have not yet begun runaway nuclear burning. The right panel shows the mass density profiles of a slice through the simulation in the x–y plane. In addition to the ambient medium (white in the figure), there are two distinct regions of density: the uncollided WD material and the stalled shock region. The density and temperature are not yet high enough to fuel runaway burning. The lower left panel of Figure 3 shows these quantities as well as the sound speed and velocity in the x-direction along a line connecting the centers of the two WDs and parallel to the x-axis. The sound speed is lower than the infall velocity speed, causing the stalling of the shock region. The temperature profile peaks at ≈10⁹ K, which is not hot enough to reach the carbon burning threshold. The upper left panel of Figure 3 shows the state of the collision in the density–temperature plane. The color of the points represents the primary composition of the corresponding cell; green for ¹²C, blue for ²⁸Si, and red for ⁵⁶Ni. Material with T ≈ 10⁷ K represents the cold and dense parts of the two stars that have not yet collided. The region with 10⁷ < T < 10⁹ K and 10⁴ < ρ < 10^6.5 g cm⁻³, represents the shocked material. At this point in the collision, "tracks" run from the lower left to the upper right, representing tori of material orthogonal to the x-axis at the center of the collision. In this case, ²⁸Si and ⁵⁶Ni have not yet been produced, thus all the cells are primarily composed of ¹²C.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Analysis images of the 2 × 0.64 collision at t = 6.00 s, after first contact, but before ignition. Top left: locations of all cells in the density–temperature plane. The color of the points represents the primary composition of the corresponding cell: green for ¹²C, blue for ²⁸Si, and red for ⁵⁶Ni. The data are binned into 100 equally spaced bins in logarithmic density and temperature. Bottom left: temperature, x-velocity, density, and sound speed along the x-axis. Right: a two-dimensional slice of density through the x–y mid-plane.
Download figure:
Standard image High-resolution image

Figure 4 has the same format as Figure 3 at 6.60 s when runaway nuclear burning has begun. On the right panel, there are three distinct regions of the collision at this point in time: the WD material which has not yet experienced the collision; the lenticular, nearly isobaric, stalled shock region; and the central region where a detonation has begun to propagate. The detonation front is outlined by the darker colored (higher density) oval region in Figures 4. Our FLASH simulations do not resolve the initiation of the detonation. Instead, at all spatial resolutions investigated, the center-most cell in the 2 × 0.64 head-on collision model undergoes runaway carbon burning which begins to propagate a detonation.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Same format as Figure 3, when the model is at t = 6.60 s, right after ignition.
Download figure:
Standard image High-resolution image

In the lower left panel of Figure 4, again, three distinct regions are visible—the unshocked WDs the stalled shock, and the center-most detonation region. The temperature in the unshocked WD material rises smoothly from the initially imposed background temperature of 1 × 10⁷ K to ≈3 × 10⁷ K at the centers of both WDs because of low-amplitude velocity waves sloshing around the WD interiors. However, 3 × 10⁷ K is well below the carbon burning threshold, does not lift the electron degeneracy of the material, and does not impact our results. In the unshocked region, the infall speed of material is greater than the local sound speed. The material behind the stalled shock reaches temperatures that are sufficient to lift electron degeneracy and are just below the carbon burning threshold of ≈2 × 10⁹ K. The density in the stalled shock region reaches a peak of ≈2 × 10⁷ g cm⁻³. In the innermost region where a detonation front has traveled ∼5 × 10⁷ cm from the center, the temperature is ≈6 × 10⁹ K and the density dips to ≈1 × 10⁷ g cm⁻³. In the upper left panel of Figure 4, hot, dense material with T > 10⁹ K and ρ > 10⁷ g cm⁻³ from the central regions of the collision are in the upper right corner where the original carbon material has burned to ²⁸Si and ⁵⁶Ni.

Figure 5 has the same format as Figure 3 and the right panel shows the density profile when the detonation front has traveled outward from the center and the densest parts of the WDs are about to enter the stalled shock region. The upper left panel indicates that more ¹²C material is present in the high density regime with ρ > 10⁷ g cm⁻³, and being burned to ²⁸Si and ⁵⁶Ni. The lower left panel shows that the sound speed in the burned region is comparable with the speed of the infalling material, and the width of the detonation has expanded.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Same format as Figure 3, when the model is at t = 6.79 s, as the stalled shock region slightly expands and the densest parts of the white dwarfs begin to enter the stalled shock region.
Download figure:
Standard image High-resolution image

As additional energy from nuclear burning is added, the double WD system eventually becomes gravitationally unbound. Figure 6 has the same format as Figure 4. The right panel shows the density distribution of the system slightly before the explosion reaches homologous expansion. The innermost 10⁹ cm reaches a nearly constant temperature of ≈3 × 10⁹ K with a slowly varying density distribution that peaks at ≈5 × 10⁶ g cm⁻³. The density–temperature plot in the upper left panel indicates larger amounts of high density, high temperature material with ρ > 10⁶ g cm⁻³ and T > 10⁹ K. More material has achieved the conditions necessary to synthesize ²⁸Si (blue) and ⁵⁶Ni (red). The lower left panel shows the sound speed is always greater than the infall speed of the remaining material.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Same format as Figure 3 at t = 7.38 s, just before the system becomes gravitationally unbound.
Download figure:
Standard image High-resolution image

The 0.81–0.81 M_☉ (hereafter 2 × 0.81) collision models evolve through a similar set of stages as the 2 × 0.64 collision models, except the larger kinetic energy at impact is sufficient for the initial shock to raise the temperature well above the ¹²C+¹²C threshold. Figure 7 shows that the entire stalled shocked region burns rapidly to a state of nuclear statistical equilibrium and achieves a nearly isothermal state. Central ignition does not occur because the ¹²C+¹⁶O material has already been burned to nuclear statistical equilibrium. We discuss this difference in additional detail in Section 3.3.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Same format as Figure 3 at t = 4.02 s for the 2 × 0.81 collision model.
Download figure:
Standard image High-resolution image

Otherwise, the stages of the 2 × 0.81 collision are very similar to the evolution of the 2 × 0.64 collision seen above, with the 2 × 0.81 collision producing a greater amount of ⁵⁶Ni.

3.2. Numerical Convergence

To assess the numerical convergence, we performed the 2 × 0.64 and 2 × 0.81 simulations at three different spatial resolutions. Each increase in spatial resolution is a factor of two more refined in one dimension, a factor of eight in volume (see Table 1), and takes at least twice as many time steps depending on the burning time step. As the spatial resolution increases, the cells that are burning carbon to heavier elements become smaller in volume and the time step decreases, leading to improved coupling between the hydrodynamics and nuclear burning.

Table 2 lists the ejected masses for each of the six convergence simulations and Figure 8 shows the convergence behavior of ¹²C+¹⁶O, ²⁸Si, ⁵⁶Ni yields, as well as the internal energy, kinetic energy, and the total energy at the end of the simulation. The upper plot in Figure 8 shows that for the 2 × 0.64 collision the ⁵⁶Ni mass (dashed red) increases, the ²⁸Si mass (dashed blue) decreases, and the ¹²C+¹⁶O (dashed green) decreases as the spatial resolution increases. The percent change in ⁵⁶Ni production is 138% between the R = 5.19 × 10⁷ cm and R = 2.59 × 10⁷ cm models, and 3.3% between the R = 2.59 × 10⁷ cm and R = 1.30 × 10⁷ cm models. Although higher resolution models are needed to reach numerical convergence, the ⁵⁶Ni mass is approaching convergence at 0.32 M_☉ (see Table 2). The internal energy (solid green) at the end of the 2 × 0.64 collision simulation decreases with increasing spatial resolution, but the kinetic energy (solid blue) when the model terminates increases with increasing spatial resolution. The net result is that the total energy (solid red) is nearly constant over the range of resolutions explored.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Convergence plot for the 2 × 0.64 (top) and 2 × 0.81 (bottom) head-on collisions. Points from left to right correspond to 5-, 6-, and 7-level runs for each collision. The dashed line colors represent different isotopes, where blue corresponds to ²⁸Si, green to ¹²C+¹⁶O, and red to ⁵⁶Ni. The solid line colors represent internal energy (green), kinetic energy (blue), and total energy (red).
Download figure:
Standard image High-resolution image

Table 2. Ejected Masses

M₁, M₂	R	¹²C + ¹⁶O	²⁸Si	⁵⁶Ni
(M_☉)	(10⁷ cm)	(M_☉)	(M_☉)	(M_☉)
0.64	5.19	0.29	0.45	0.13
0.64	2.59	0.21	0.37	0.31
0.64	1.30	0.19	0.37	0.32
0.81	4.32	0.19	0.41	0.62
0.81	2.16	0.19	0.50	0.45
0.81	1.08	0.18	0.53	0.39

Download table as: ASCII Typeset image

The lower panel in Figure 8 shows that for the 2 × 0.81 collision the ⁵⁶Ni mass decreases, the ²⁸Si mass increases as the spatial resolution increases, and the ¹²C+¹⁶O slightly decreases. Although convergence has not been achieved, the ⁵⁶Ni mass is approaching convergence at 0.39 M_☉ (see Table 2). We discuss the reason for the different convergence trends between the 2 × 0.64 and 2 × 0.81 cases below. The internal energy and kinetic energy at the end of the 2 × 0.81 collision simulations appears to be oscillating toward convergence as the spatial resolution is increased. As a consequence of the internal energy and kinetic energy being out of phase, the total energy is nearly constant over the range of resolutions explored.

Although strict numerical convergence has not been achieved with these six simulations, some trends can be seen. As the total mass of the binary system increases in zero impact parameter collisions, the ⁵⁶Ni mass increases, indicating that larger mass collisions will produce more ⁵⁶Ni. For both mass pairs at highest resolution, ≈0.2 M_☉ of unburned ¹²C+¹⁶O was ejected.

Higher numerical resolutions are desirable, but prohibitively expensive for this study, as our most resolved three-dimensional models required at least 200,000 CPU hours per run. Simulations with higher spatial resolution are not possible in the context of the current study because doubling the grid resolution in a three-dimensional simulation effectively increases the number of cells by a factor of ≈2³ and the number of time steps by a factor of two, meaning over an order-of-magnitude increase in computational time. Although these effects can be ameliorated somewhat by adopting more aggressive derefinement criteria, further restricting the computational domain size, or relaxing the time step controller f, we expect that increasing the maximum resolution another factor of two (6.5 × 10⁶ cm for the 2 × 0.64 models and 5.04 × 10⁶ cm for the 2 × 0.81 models) would require ≈2 million CPU hours per run, which is beyond our capabilities here.

Reducing the time step limiting factor, f, and thereby reducing the time step during nuclear burning changes the amount of ⁵⁶Ni produced. For example, changing from f = 0.5 to f = 0.1 in the 2 × 0.64 simulation with a spatial resolution of R = 2.59 × 10⁷ cm causes the ⁵⁶Ni production to increase by approximately 0.1 M_☉, a 30% change. Figure 9 shows the evolution of the hydrodynamic time step (solid lines), burning time step (dotted lines), and ⁵⁶Ni mass (dashed lines) for the 5-level (red), 6-level (green), and 7-level (blue) 2 × 0.81 collisions. We use f = 0.2 for the 5-level run and f = 0.3 for the 6- and 7-level runs to force the burning time step to fall below the hydrodynamic time step during the ⁵⁶Ni production phase. In all our simulations, we set f such that dt_burn ≈ 0.01 dt_hydro during the phase of evolution when nuclear burning is significant. Setting f to smaller values greatly increases the computing time without having a significant effect on the nucleosynthesis yields.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Evolution of the hydrodynamic time step (solid line), burning time step (dotted line), and ⁵⁶Ni mass (dashed line) for the 5-level (red), 6-level (green), and 7-level (blue) models of the 2 × 0.81 collision.
Download figure:
Standard image High-resolution image

Figure 8 shows that the 2 × 0.64 collision produces more ⁵⁶Ni as spatial resolution increases. To understand this behavior we examine profiles along the x-axis of the density and temperature for 5-, 6-, and 7-levels of refinement. The upper panel of Figure 10 shows the three models with different spatial resolutions for the 2 × 0.64 collision at 5.6 s. The shocked region is widest in the 5-level model, and narrower in the 6- and 7-level models. The density is smaller (just below 3 × 10⁶ g cm⁻³) and nearly constant for the 5-level model, larger for the 6-level model than the 5-level model (3.5 × 10⁶ g cm⁻³), and slightly larger for the 7-level model than the 6-level model (3.6 × 10⁶ g cm⁻³). Both the 6- and 7-level models show a small valley in the central density. The peak temperature is smaller for the 5-level (≈10⁹ K), and slightly higher for the 6-level and 7-level models (both above 10⁹ K).

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Density and temperature profiles along the x-axis for the 2 × 0.64 collision at 5.6 s, 6.4 s, and 6.9 s for different levels of refinement.
Download figure:
Standard image High-resolution image

At 6.4 s (middle panel in Figure 10), the density and temperature profile patterns as described for 5.6 s generally still hold, but the peak temperature is now the same for all three levels of refinement (≈1.5 × 10⁹ K). Detonation occurs just after this time frame (as seen below). Thus, we expect to see more ⁵⁶Ni produced for the 6-level and 7-level models than for the 5-level model because there is more material in the shocked region with high density (>10⁷ g cm⁻³) at the same ignition temperature. We also expect only a small difference in ⁵⁶Ni production between the 6- and 7-level models because the peak density is only slightly larger for the 7-level model and the width of the density profile is approximately the same. This explains the pattern in the abundance yields with spatial resolution in Figure 8.

At 6.9 s (lower panel in Figure 10) when the detonation is underway, the Mach number is larger in the 7-level model than the 5- and 6-level models because the pre-detonation density is larger. This causes the 7-level temperature profile to be wider then the 6-level or the 5-level. That is, the burning front travels farther for the same amount of time.

Unlike the 2 × 0.64 collision, the 2 × 0.81 collision produces less ⁵⁶Ni as level of refinement increases (see Figure 8). The top panel of Figure 11 shows the three models with different spatial resolutions for the 2 × 0.81 collision at 4.0 s. The temperatures for all three resolutions are high (>3 × 10⁹ K), indicating the energy generated by burning is large. The shocked region is widest for the 7-level model and narrowest for the 5-level model with the 6-level model in-between. The 7-level model has the lowest, and nearly constant, density (≈8 × 10⁶ g cm⁻³) in the stalled shock region, and has the largest magnitude spikes in the density (just below 2 × 10⁷ g cm⁻³) at the edges of the stalled shock. The spikes occur because the density of material is highest just behind the shock front. The 6-level model has a slightly larger (≈8 × 10⁶ g cm⁻³), nearly constant, density in the middle, and slightly smaller spikes in the density (just below 2 × 10⁷ g cm⁻³) at its edges. The 5-level model has its density spikes (just above 10⁷ g cm⁻³) close enough together that a nearly constant density in the middle is barely reached.

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Density and temperature profiles along the x-axis for the 2 × 0.81 collision at 4.0 s, 4.4 s, 4.5 s, 4.7 s, and 4.9 s for different levels of refinement.
Download figure:
Standard image High-resolution image

The second panel of Figure 11 shows at 4.4 s the width of the shocked, burning region is larger for all three resolutions, because the energy generated by burning in the hot shocked region is sufficient to overcome the standing shock formed from material moving inward. That is, the shocked burning region is growing. The temperature is nearly the same and constant (≈5 × 10⁹ K) across all three resolutions, but with small spikes at the edges of each shocked region. The nearly constant density in the central region of the 7-level model is still smaller and wider than in the 6-level model. The 5-level still has its two spikes near the center, thus a nearly constant central density region is not reached.

The third panel of Figure 11 shows at 4.5 s the 7 level model begins to detonate, but the 6-level and 5-level models have not yet detonated. The same patterns in density and temperature described for previous time points still hold. By 4.7 s, the fourth panel of Figure 11 shows the 6-level model begin to detonate, but the 5 level model has not yet detonated. The width of the burning region for the 6-level model is about the same width as in the 7-level model when the 7-level model detonated 0.2 s earlier.

At 4.9 s (bottom panel of Figure 11), the 5-level model is the last to detonate. The 5-level model has finally reached a state of nearly constant density in the central region with spikes at the edges. This nearly constant density of 2 × 10⁷ g cm⁻³ is larger than the nearly constant density reached by either the 6-level or the 7-level models (both about 10⁷ g cm⁻³), but it has reached about the same width. Since the 7-level model detonates at the lowest density (but the same mass since all reach about the same width before detonating) and soonest in time, the 7-level model should produce the least amount of ⁵⁶Ni. The 5-level model detonates at a higher density (and same mass) and latest in time, thus should produce the most ⁵⁶Ni. This explains the pattern in the abundance yields with spatial resolution in Figure 8.

3.3. Similarities and Differences between the Explosion Models

Whether the explosion is initiated along the edge of the stalled shock region (as in the 2 × 0.81 collisions) or in the central regions of the stalled shock (as in the 2 × 0.64 collisions) is controlled by the initial masses of the WDs, as the masses set the infall speed (escape velocity). The infall speed determines the strength of the initial shock, and thus the initial post-shock temperature. In turn, the initial post-shock temperature determines the amount of initial burning behind the shock, and hence the temperature profile of the shocked region. Comparing the temperature profiles in the shocked region between the 2 × 0.64 and 2 × 0.81 collisions, we see that the 2 × 0.64 temperature barely reaches 10⁹ K, while the 2 × 0.81 temperature is a hot ≈5 × 10⁹ K over an extended region. The difference in temperature between the two-model collisions is a direct consequence of the kinetic energy of the collision (larger kinetic energy corresponding to larger temperature).

The initial shock in the 2 × 0.64 collision models barely raises the temperature above the ¹²C+¹²C threshold. As carbon burning proceeds, the central shocked burning region increases in temperature. The top panel of Figure 12 shows the system cannot yet explode because the temperature is not high enough to overcome the ram pressure from the infalling material, which continues to increase due to density profile of the WD. Only when the peak of the WD density profiles enters the shocked region does the central peak undergo thermonuclear runaway, which creates enough pressure to overcome the now decreasing ram pressure (see bottom panel of Figure 12).

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Density, temperature, pressure, and ram pressure profiles along the x-axis for the 2 × 0.64 collision at 5.55 s and 6.55 s.
Download figure:
Standard image High-resolution image

In contrast to the 2 × 0.64 collision model, the 2 × 0.81 collision model is energetic enough that the initial shock raises the temperature well above the ¹²C+¹²C threshold. The entire stalled shocked region burns rapidly to a state of nuclear statistical equilibrium and achieves a nearly isothermal state. Central ignition cannot occur because the ¹²C+¹⁶O material has already lost nearly all of its energy in the burn to nuclear statistical equilibrium. The top panel of Figure 13 shows, similar to the 2 × 0.64 collision models, that the 2 × 0.81 collision model cannot yet explode since the ram pressure from the infalling material is greater than the pressure of the hot burned material pushing outward. Only when the pressure inside the hot burned region is larger than the ram pressure does the system explode, giving the appearance of an edge-lit ignition (see bottom panel of Figure 13).

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Density, temperature, pressure, and ram pressure profiles along the x-axis for the 2 × 0.81 collision at 3.60 s and 4.60 s.
Download figure:
Standard image High-resolution image

4. DISCUSSION

We have performed the first systematic study of zero impact parameter collisions between two WDs with an Eulerian grid code (FLASH). Our simulations spanned a range of effective spatial resolutions for collisions between two 0.64 M_☉ WDs and two 0.81 M_☉ WDs. However, even the highest resolution studies did not achieve strict numerical convergence.

The lack of convergence in the simplest configuration (zero impact parameter, equal masses) suggests that quantitative predictions of the ejected elemental abundances that are generated by binary WD collision and merger simulations should be viewed with caution. However, the convergence trends do allow some patterns to be discerned.

We found the 2 × 0.64 collision model head-on collision model produces 0.32 M_☉ of ⁵⁶Ni, 0.38 M_☉ of ²⁸Si, and 0.2 M_☉ of unburned ¹²C+¹⁶O. Rosswog et al. (2009) included one FLASH based model of a zero impact parameter collision of two 0.60 M_☉ WDs in their study. They reported an ⁵⁶Ni mass of 0.16 M_☉, about one-half of what we find. While Rosswog et al. (2009) used slightly less massive WDs than our study, both sets of FLASH simulations used the same equation of state. The FLASH model of Rosswog et al. (2009) achieved about a factor of 2.6 smaller spatial resolution than our study R = 4.9 × 10⁶ cm versus R = 1.3 × 10⁷ cm), due to their evolving one WD and deploying a mirrored gravitational potential. This difference in the maximum spatial resolution could account for the different ⁵⁶Ni masses between the two calculations, although the convergence trend shown in the upper panel of Figure 8 suggests spatial resolution might not be the only reason for the difference. Another potential reason for the difference in the ⁵⁶Ni masses is the choice of the time step, and thus the coupling between the operator split processes of hydrodynamics and the nuclear burning. In all our simulations, we limited the time step to ≈0.01 of the Courant-limited hydrodynamic time step during the nuclear burning phases (see Figure 9). We found changing the allowed time step can change the ⁵⁶Ni mass produced by 30%–40%.

We find our FLASH-based, zero impact parameter, collision models systematically produce less ⁵⁶Ni and more silicon-group elements than collisions models calculated with SNSPH by Raskin et al. (2010). This difference between particle and grid based codes was first found by Rosswog et al. (2009), who suggested that differences in nuclear reaction networks or advection effects could be responsible for the different yields. While our FLASH models used the same equation of state and nuclear reaction network as Raskin et al. (2010), and we checked the same output values were returned for the same input values, a detailed investigation of the differences between our FLASH model results and the Raskin et al. (2010) results with SNSPH are beyond the scope of this paper.

Red and dim SNe Ia such as SN 1991bg (Leibundgut et al. 1993; Turatto et al. 1996; Hachinger et al. 2009), SN 1992K (Hamuy et al. 1994), SN 1999by (Garnavich et al. 2004), and SN 2005bl (Taubenberger et al. 2008) are characterized by M_V ≈ −17. The light curves of underluminous SNe Ia decline even more rapidly than expected from a linear luminosity to decline-rate relation among normal SNe Ia (Phillips et al. 1999; Taubenberger et al. 2008; Blondin et al. 2012). Spectroscopically, 91bg-like SNe Ia show low-line velocities around B-magnitude maximum (Filippenko et al. 1992) and clear spectral signatures of Ti-II, indicating lower ionization (Mazzali et al. 1997). Taken together, these properties together are consistent with ∼0.1 M_☉ of newly synthesized ⁵⁶Ni. Our FLASH models suggest 2 × 0.64 M_☉ and 2 × 0.81 M_☉ head-on collision models produce ⁵⁶Ni masses below that needed for normal SNe Ia, but are within a range consistent with observations of underluminous SNe Ia. In addition, either mass pairing produces ∼0.2 M_☉ of unburned C+O, which may be a unique signature of mergers and collisions between WDs.

Future studies should include a survey of non-zero impact parameter WD collisions with FLASH, an exploration of unequal mass collisions, and an investigation why Lagrangian particles codes and Eulerian grid codes continue to find about a factor of two difference in the mass of ⁵⁶Ni ejected. The zero impact parameter is insightful as an upper limit on ⁵⁶Ni production, but a non-zero impact parameter study will likely give a range of ⁵⁶Ni yields for different collision configurations. An exploration of unequal mass collisions could provide a broader physical parameter space and allow an improved quantitative description of how SNe Ia luminosity scales with mass pairings.

This work was supported by the National Science Foundation under grant AST 08-06720 and through the Joint Institute for Nuclear Astrophysics (JINA) under grant PHY 02-16783. All simulations were conducted with Arizona State University Advanced Computing Center and Extreme Science and Engineering Discovery Environment (XSEDE) compute resources. FLASH was in part developed by the DOE-supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago. W.P.H. thanks Brandon Mechtley for his invaluable computing assistance.

ZERO IMPACT PARAMETER WHITE DWARF COLLISIONS IN FLASH

Article metrics

Permissions

Share this article

Dates

ABSTRACT