THE IGNITION OF CARBON DETONATIONS VIA CONVERGING SHOCK WAVES IN WHITE DWARFS

The framework used to calculate planar post-detonation conditions is often referred to as the CJ solution, after the original work of Chapman (1899) and Jouguet (1905). From mass, momentum, and energy conservation, and the assumption that the burned ashes move at the speed of sound in the shock's rest frame, the CJ detonation velocity is $v_{\rm CJ} = \sqrt{ 2 (\gamma ^2\hbox{--}1) q }$. Here q is the energy per mass released from burning the fuel to ash, and the equation of state exponent, γ, is typically near 1.4 for our relevant conditions, but must be calculated self-consistently.

For our CJ calculations, we take the end state of burning to be ²⁸Si because our successful detonations quickly burn the C/O to a state of quasi-nuclear statistical equilibrium (quasi-NSE), consisting of isotopes with binding energies near that of ²⁸Si. Given enough time, the material in the propagating detonations will burn all the way to NSE, which involves a mix of IGEs. However, this occurs on lengthscales and timescales that are orders of magnitude larger than required to burn to quasi-NSE. Since we are only concerned with small volumes surrounding the focal point, the successes of the detonations in this paper are determined solely by the binding energy released from converting C/O to quasi-NSE. Table 1 shows the CJ results under this assumption for our different initial conditions.

We now consider the simplified case of spherical imploding shock waves in the absence of chemical reactions. These types of shocks are referred to as CCW shocks due to the pioneering studies of Chester (1954), Chisnell (1955, 1957), and Whitham (1957, 1958). Thorough analytic and numeric work on these imploding shock waves has already been performed (e.g., Guderley 1942; Stanyukovich 1960; Zel'dovich & Raizer 1967; Landau & Lifshitz 1987; Ponchaut et al. 2006; Kushnir et al. 2012), so we only summarize their results.

The strength of a spherically symmetric shock wave in a converging medium increases as the surface area of an imploding shock front decreases, and decreases as an exploding shock expands. The Mach number of the shock wave, M, scales with the shock's distance from the focal point as M∝r^α (Chisnell 1957), where

$$\begin{eqnarray*} &&\alpha = - 2 \frac{M^2\hbox{--}1}{\lambda M^2 }, \nonumber \end{eqnarray*}$$

$$\begin{eqnarray*} &&\lambda = \left(2 \sigma + 1 + \frac{1}{M^2} \right)\left[ 1 + \frac{2 (1-\sigma ^2)}{\sigma \left(\gamma +1 \right)} \right]{\rm,\ and} \nonumber \end{eqnarray*}$$

$$\begin{eqnarray} &&\sigma ^2 = \frac{ \left(\gamma -1 \right)M^2+2}{ 2 \gamma M^2- \left(\gamma -1 \right)}. \end{eqnarray} \tag{ 1 }$$

The equation of state exponents, sometimes referred to as Γ₁ and Γ₃, are presumed to be equal, as appropriate for the cases of ideal gas and radiation, and are denoted as γ. In the limit of a strong shock with γ = 1.4, which will be relevant for our future calculations, the Mach number scales as M∝r^−0.39. Because α depends implicitly on M, it must be calculated numerically.

These results are derived under the assumption that the evolution of the previously shocked material does not significantly affect the shock properties and that γ remains the same before and after the shock. In spite of these assumptions, the analytic relations compare very well with numerical hydrodynamics results, as we show in the next section.

As a first test of the numeric code, we compare the evolution of a purely hydrodynamic (i.e., non-reactive) imploding shock wave to the analytic results from Section 3. Figure 1 shows radial pressure (top panel) and velocity (bottom panel) profiles of an imploding shock wave at 15 snapshots in time. The initial density is 3.2 × 10⁷ g cm⁻³, and the snapshots are separated by 5 × 10⁻⁸ s and begin 7.5 × 10⁻⁷ s before the shock wave has reached the focal point. Figure 2 shows the second half of this calculation, beginning just after the shock wave has reached the focal point and reversed its direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the post-shock pressure versus the position of the shock front for the same calculation shown in Figures 1 and 2. The relative pressure jump across a shock front is

$$\begin{eqnarray} &&\frac{P_1}{P_0} = \frac{2 \gamma M^2 - (\gamma -1)}{\gamma +1}, \end{eqnarray} \tag{ 2 }$$

so the power-law scaling of post-shock pressure with shock radius will be roughly twice as strong as the Mach number's scaling; thus, P∝r^−0.79 for a strong shock with γ = 1.4. This relation is shown as a red dotted line in Figure 3. Both the imploding and exploding shock fronts follow this scaling fairly well. The outgoing shock has a higher normalization than the ingoing because it propagates outward into previously shocked and converging material.

Zoom In Zoom Out Reset image size

When nuclear reactions are included, burning begins as the imploding shock wave approaches the focal point and the post-shock conditions reach burning temperatures and densities. Figures 4 and 5 show a calculation with the same initial conditions as the run in Figures 1 and 2, but with nuclear reactions turned on. As Figure 4 demonstrates, nuclear reactions yield an outwardly propagating detonation before the imploding shock wave has reached the focal point. The outgoing detonation can be seen as a growing spike in the pressure and velocity profiles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Throughout this paper, we use the radius where the imploding shock velocity equals the planar CJ detonation velocity, r(v_shock = v_CJ), as a proxy for the strength of the initial imploding shock. A simulation with a larger value of r(v_shock = v_CJ) has a stronger initial shock. This radius is essentially equivalent to the point at which heating due to nuclear reactions overtakes compressional heating due to the converging shock flow. For the fiducial calculation shown in Figures 4 and 5, this radius is r(v_shock = v_CJ) = 96 cm and is shown as a dotted line.

The spherically symmetric geometry allows us to resolve the very small lengthscales that characterize C detonations at these densities. Figure 6 shows the post-shock structure of a successful detonation at the end of the calculations shown in Figures 4 and 5. While the initial density of the material was 3.2 × 10⁷ g cm⁻³, the upstream material at this stage of the calculation has a density of 5.9 × 10⁷ g cm⁻³ due to the converging shock flow. Solid lines denote the mass fractions of the most abundant isotopes as labeled, while the dashed line shows the energy generation rate, , normalized to its maximum value of 4 × 10²⁶ erg g⁻¹ s⁻¹. The shortest detonation lengthscale is associated with the consumption of ¹²C (0.2 cm), followed by the location of the maximum of the energy generation rate (0.3 cm), and trailed by the ¹⁶O consumption lengthscale (1.4 cm). Note that while the first reactions to take place behind the shock front are ¹²C+¹²C self-reactions, the burning material soon approaches quasi-NSE, consisting primarily of ²⁸Si and ³²S. Due to the much higher temperatures reached near the focal point, which is 2.7 × 10³ cm farther downstream and not shown in the figure, the composition of the detonation ashes is much closer to full NSE.

Zoom In Zoom Out Reset image size

Figure 7 shows a summary of these measures of the detonation lengthscale versus upstream density for calculations with successfully propagating detonations. The initial mass fractions are $X_{\rm ^{12}C}=X_{\rm ^{16}O}=0.5$. As in Figure 6, the shortest lengthscale is the distance to the location where $X_{\rm ^{12}C}$ has been halved, labeled $l_{\rm ^{12}C}$. The location of the maximum of the energy generation rate, labeled l, is at a distance that is a factor of 1–2 longer than $l_{\rm ^{12}C}$. The lengthscale for ¹⁶O to be halved, labeled $l_{\rm ^{16}O}$, is significantly longer than both $l_{\rm ^{12}C}$ and l. The values of these lengthscales are in good agreement with previous work (Khokhlov 1989; Gamezo et al. 1999).

Zoom In Zoom Out Reset image size

The success of a detonation depends on the initial strength of the imploding shock. Figures 8 and 9 compare a successfully propagating detonation and an unsuccessful detonation, respectively. Both figures show radial profiles of the pressure (black lines), normalized to 10²⁶ dyne cm⁻², and of the ¹²C mass fraction (red lines) at the same relative time in each simulation, measured from the moment when the imploding shock wave reaches the focal point. Both calculations begin with an initial density of 3.2 × 10⁷ g cm⁻³, but the outer boundary for the calculation in Figure 8 is at 2.75 × 10³ cm, while the simulation in Figure 9 has an outer boundary at 1.375 × 10³ cm. As a result, the lengthscales and timescales in the successful calculation are twice as large as in the unsuccessful run. The larger value of r(v_shock = v_CJ) in Figure 8 implies a larger initial shock strength, which is why it forms a successfully propagating detonation, as demonstrated by the superposition of the shock front and the compositional discontinuity. The simulation with the weaker initial shock fails to yield a detonation, as shown by the lack of coupling of the shock front and the compositional discontinuity, whose velocity in mass space has stalled.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These expanding detonations fail to propagate for two reasons. First, the velocity is constrained to be zero at the focal point, located behind the outwardly propagating detonation, while for a standard CJ detonation, the downstream ashes move at a finite velocity in the direction of the detonation. The zero-velocity boundary condition exerts a backward "pull" on the material behind the detonation, preventing it from reaching CJ conditions. This effect has been demonstrated for outwardly propagating explosions with a zero-velocity boundary condition in planar and spherical geometry by He & Clavin (1994) and Seitenzahl et al. (2009), among others. While a simple estimate of the lengthscale of the critical shocked volume necessary for a subsequent successful detonation might be ∼l, the actual critical lengthscale for successful planar C detonations with a zero-velocity boundary condition is a factor of 10³ to 3 × 10⁴ times larger.

The second major effect is due to the spherical curvature of the detonation front. As the post-shock, but pre-burned, material expands, the density and temperature are reduced and nuclear reactions proceed at a slower pace than for a planar detonation at the same velocity. Furthermore, when burning releases most of its energy at some finite distance ∼l behind the shock, it powers a shock front with a larger surface area than in the planar case. The effect of curvature increases the critical radius by an additional factor of three as compared to l, so that the critical radius found by Seitenzahl et al. (2009) was 3 × 10³ to 10⁵ times larger than the detonation lengthscale, l.

However, as we show in the next section, our critical radii are 300–10⁴ larger than the detonation lengthscales for material at the initial density. This somewhat smaller increase is due to the converging shock flow: the outwardly moving detonation propagates into previously shocked material, which consequently has a shorter detonation lengthscale than it would if it had the initial unshocked density.

By varying the size of the computational domain and thus the imploding shock strength, we can map the regions of parameter space that yield successful or unsuccessful detonation ignitions. Figure 10 shows this dividing line at different initial densities, characterized by the radius at which the imploding shock velocity would equal the CJ detonation velocity if nuclear reactions were neglected. The critical radius decreases as the initial density increases and ranges from 10²–10⁶ cm for typical WD densities.

Zoom In Zoom Out Reset image size

Also shown in Figure 10 are approximations of radii at which the converging shock wave reaches the CJ velocity from two-dimensional full star hydrodynamic simulations of double detonations. The diamond, circle, and triangle represent calculations with the PROMETHEUS code (Fryxell et al. 1989) of 0.45 M_☉ C/O core + 0.21 M_☉ He shell (Model L of Sim et al. 2012), 0.81 + 0.13 M_☉ (Model 1 of Fink et al. 2010), and 1.03 + 0.06 M_☉ (Model 3 of Fink et al. 2010) simulations, as labeled. The square represents a 1.02 + 0.03 M_☉ simulation (D. Townsley 2013, private communication) calculated with FLASH (Fryxell et al. 2000). The smallest resolution in these multi-dimensional simulations is 1–3 × 10⁶ cm. The He detonation in each of these simulations is initiated at a point.

Since the converging shock waves in the multi-dimensional simulations are somewhat aspherical, a direct mapping of their shock strengthening to our spherically symmetric calculations requires an approximation. We first estimate the radius of the imploding shock wave at a given time by calculating the enclosed volume within the shock, V_encl, and then inferring the spherically averaged radius, (3V_encl/4π)^1/3. The time evolution of this quantity is used to estimate the shock's inward velocity at two radii, from which the value of the shock strengthening scaling is deduced. These quantities are then used to estimate the radius at which the imploding shock reaches the CJ velocity, which is always >10 km.

The spherically averaged strength of the imploding shock in these multi-dimensional simulations is many orders of magnitude above our critical values for initial densities ≳ 10⁷ g cm⁻³, which is the central density of a 0.8 M_☉ WD. Thus, it appears that propagating He detonations can robustly ignite high-mass C/O cores via converging shock waves, even when the small 0.01–1 cm C detonation lengthscales are resolved. However, given the much higher critical r(v_shock = v_CJ) at lower densities, detonations in lower-mass C/O cores are not as certain. A complete analysis of double detonation ignition necessitates adequately resolved multi-dimensional simulations.

To verify that our results are converged, we performed a suite of simulations for each of our three values of ρ₀, using three ratios of ρ_pert/ρ₀ = 1.1, 1.2, and 1.5 and varying the minimum resolution. One representative convergence study is shown in Figure 11, which demonstrates the effect on the critical value for r(v_shock = v_CJ) when the minimum resolution is varied for an initial density of ρ₀ = 10⁷ g cm⁻³ and perturbed density ρ_pert = 1.5ρ₀. Black circles denote simulations with successful detonations; red crosses mark those with failed detonations. The necessity of resolving the burning lengthscale (∼0.3 cm for this example) is made clear in simulations with inadequate resolution, which can have successful detonations for imploding shocks that reach the CJ velocity at radii 20% smaller than the true critical value of 1000 cm.

Zoom In Zoom Out Reset image size

Figure 12 shows the upper and lower bounds on the critical values for r(v_shock = v_CJ) for various ratios of the perturbed density to the initial density, or alternatively, the ratio of r(v_shock = v_CJ) to the size of the computational domain, for our three initial densities as labeled. Smaller critical values of r(v_shock = v_CJ) are found for simulations with larger ρ_pert/ρ₀. This is due to the zero-gradient boundary condition at the outer edge of the grid. The density and pressure of the once-shocked imploding material at a fixed time should decrease with increasing radius. The zero-gradient condition thus implies densities and pressures of the inflowing material that are higher than they should be, and when the detonation propagates into this material, the burning lengthscales are correspondingly shorter. This in turn yields a smaller critical value for r(v_shock = v_CJ).

Zoom In Zoom Out Reset image size

It is thus important to ensure that the value of ρ_pert/ρ₀ is low enough to remove the effect of the outer boundary condition. The critical values of r(v_shock = v_CJ) do indeed converge by ρ_pert/ρ₀ = 1.1 for the higher initial densities of 10⁷ and 3.2 × 10⁷ g cm⁻³. In contrast, the derived critical value for the lowest initial density of 3.2 × 10⁶ g cm⁻³ is not fully converged even at our lowest value of ρ_pert/ρ₀ = 1.1, and the computational cost of probing lower values of ρ_pert/ρ₀ became prohibitive. However, the exhibited trend implies only a small change to this value for smaller perturbations. More importantly, the size of the domain for these lowest density simulations approaches the size of a WD, and our assumption of a constant initial density throughout the unperturbed region breaks down. The critical shock strength for this density should thus only be viewed as an approximation, and very likely a lower limit.

Kushnir et al. (2012) performed an earlier study of detonation initiation from imploding shock waves. While they apply a more generalized analysis to the problem and focused on reactions with Arrhenius-type temperature dependences, their order-of-magnitude estimate for C/O fuel with ρ₀ = 10⁷ g cm⁻³ is r(v_shock = v_CJ) = 10³ cm, which is in good agreement with our results.

The spontaneous initiation of spherical C detonations from regions with perturbed temperatures, but unperturbed velocities, has been previously explored by Arnett & Livne (1994), Niemeyer & Woosley (1997), Röpke et al. (2007), and Seitenzahl et al. (2009). These authors calculate the critical sizes of hot regions that ignite and yield propagating detonations for various parameterizations of the temperature profile. Direct comparison of our results to theirs is difficult, as the density ahead of the outwardly propagating detonation in our calculations changes with time, but for a density of 3 × 10⁷ g cm⁻³ and surrounding temperature of 10⁹ K, Seitenzahl et al. (2009) find critical radii ranging from 3 × 10³ to 10⁵ cm, depending on the temperature profile of the initiating volume. The detonation lengthscale at this density is ≃ 0.3 cm, so the ratio of their critical radius to detonation lengthscale is 10⁴ to 3 × 10⁵.

In our calculations, for an initial density of 3.2 × 10⁷ g cm⁻³, our critical radius is ≃ 100 cm, significantly smaller than the lower end of Seitenzahl et al. (2009)'s results. This is unsurprising because the outward detonation in our case initially propagates into previously shocked material, whose higher density yields smaller detonation lengthscales. For this initial density of 3.2 × 10⁷ g cm⁻³, the newly formed detonation propagates into material with density >10⁸ g cm⁻³ and temperature ≃ 10⁹ K. Extrapolating Seitenzahl et al.'s (2009) results to these high densities yields a range of critical radii that overlap our value of 100 cm.

We have also performed analogous calculations simulating the interiors of O/Ne WDs, with mass fractions of $X_{\rm ^{16}O}=0.7$ and $X_{\rm ^{20}Ne}=0.3$. However, the increased Coulomb barrier for O-burning results in a much longer detonation lengthscale, which is even larger by a factor of ∼10⁴ than the O-consumption lengthscale in a propagating C detonation shown in Figure 7, because O-burning is enhanced in that case by the presence of ⁴He nuclei liberated during C-burning. Thus, much higher imploding shock strengths and larger simulation volumes are required to achieve a successful detonation, which makes resolution of the detonation structure difficult. As a result, none of our O/Ne runs that spatially resolved the burning lengthscales yielded a successfully propagating detonation.

Seitenzahl et al. (2009) explored the effect of composition on the spontaneous initiation of detonations and found that decreasing the initial carbon mass fraction from $X_{\rm ^{12}C} = 0.5$ to $X_{\rm ^{12}C}=0.3$ for one of their simulations increased the critical radius for detonation ignition by a factor of 10–100. Extrapolating their results to $X_{\rm ^{12}C}=0$ suggests critical radii that are 5 × 10⁴ to 10⁸ times larger than for the case with $X_{\rm ^{12}C}=0.5$, and thus it is unsurprising that we have not resolved a successful O/Ne detonation in our calculations.

If the critical radii for O/Ne detonations are as much as 10⁶ times larger than for C/O detonations, Figure 10 suggests that double detonations do not occur if the WD core is C-deficient. This may explain why detonations of O/Ne WDs, which have masses ⩾1.2 M_☉ and would yield overluminous SNe Ia with relatively fast light curve evolution, have not been observed (Sim et al. 2010).