A review of direct numerical simulations of astrophysical detonations and their implications

Abstract

Multi-dimensional direct numerical simulations (DNS) of astrophysical detonations in degenerate matter have revealed that the nuclear burning is typically characterized by cellular structure caused by transverse instabilities in the detonation front. Type Ia supernova modelers often use one-dimensional DNS of detonations as inputs or constraints for their whole star simulations. While these one-dimensional studies are useful tools, the true nature of the detonation is multi-dimensional. The multi-dimensional structure of the burning influences the speed, stability, and the composition of the detonation and its burning products, and therefore, could have an impact on the spectra of Type Ia supernovae. Considerable effort has been expended modeling Type Ia supernovae at densities above 1\times {10}^7 g\cdotcm^-3 where the complexities of turbulent burning dominate the flame propagation. However, most full star models turn the nuclear burning schemes off when the density falls below 1\times {10}^7 g\cdotcm^-3 and distributed burning begins. The deflagration to detonation transition (DDT) is believed to occur at just these densities and consequently they are the densities important for studying the properties of the subsequent detonation. This work will review the status of DNS studies of detonations and their possible implications for Type Ia supernova models. It will cover the development of Detonation theory from the first simple Chapman–Jouguet (CJ) detonation models to the current models based on the time-dependent, compressible, reactive flow Euler equations of fluid dynamics.

Graphical abstract

Keywords

supernova / detonations / direct numerical simulations

1  Introduction

Type Ia supernovae are widely believed to be thermonuclear explosions of carbon-oxygen white dwarf (WD) stars in binary star systems, although the exact nature of these systems is still under debate. Current theories include single-degenerate progenitor systems, containing either a near-Chandrasekhar-mass or sub-Chandrasekhar-mass WD orbiting a common center of mass with a main sequence or red giant companion, or double-degenerate systems, containing two WDs of similar masses. Variations of these two categories give rise to the ever-increasing number of hypotheses that compose the current literature, each competing to fit some niche of Type Ia observations. Although discrepancies exist between competing theories, most of the modeling community will agree that detonations are a necessary component to describe their interpretation.

Detonations propagate through fuel by compressing it to burning conditions within a thin region known as a shock front. The energy that is released by the nuclear reactions compresses the adjacent fuel which burns releasing energy to compress the next segment of fuel, thus maintaining the detonation. Pure detonations were considered one of the first models of Type Ia supernovae but were shown to be inconsistent with observations. Detonations propagate at super-sonic velocities so a WD does not have sufficient time to expand as the shock proceeds through the star under high density, processing it completely into iron-group elements [1, 2]. However, a defining feature of Type Ia observations is the presence of intermediate-mass elements (elements between carbon and nickel) in the spectra [3], so the lack of this feature eliminated pure detonations as legitimate explanations. A current model that agrees well with observations is a deflagration to detonation transition (DDT). In this model the explosion begins as a sub-sonic deflagration which pre-expands the WD before transitioning into a detonation which then burns the remainder of the star at lower densities, giving rise to the intermediate-mass elements as witnessed in Type Ia supernovae spectra. Parameterized models of this hypothesis have been shown to reliably match light curves and elemental abundances necessary for a proper explanation [4], but they may not reflect reality as they are manually adjusted to fit the data. Thus, parameterized simulations do show DDT models to be plausible but a deeper understanding of the explosion mechanism is necessary to constrain current models. Transitions from deflagrations to detonations are known to occur in terrestrial experiments and simulations [5, 6], giving further reason to speculate that this scenario may be realized.

Simulations of Type Ia supernovae are notoriously plagued with the complication of a vast range of length scales that is impossible to resolve with current computing capabilities. The nuclear flame width is 10^-3 to10 cm, 8–12 orders of magnitude smaller than the stellar radius of approximately 10⁷ to 10⁸ cm [7]. In order to work around this problem the state-of-the-art technique is to employ sub-grid models to stand in for the small-scale flame physics within the overall full-star models [8]. These sub-grid models are developed from the results of direct numerical simulation (DNS) studies of detonations under various conditions encountered within a WD environment. DNS studies provide information such as propagation velocity, burning timescales, abundance data, and energy release which are unresolved at full-star lengths. Another consideration appropriate for DNS modeling is detonation initiation and survival. The collection of data from these DNS studies can then be used to inform sub-grid models and whole star Type Ia simulations, yielding realistic approximations of the large-scale explosions without the need to resolve all length scales. This article will focus on the results of such DNS studies which have been performed in order to discuss their possible implications in Type Ia SNe modeling.

2 Detonation theory

2.1 Simple detonation models

The first detonation theory was a simple one-dimensional model developed by Chapman [9] and Jouguet [10] at the turn of the twentieth century. The Chapman–Jouguet (CJ) theory models the shock front as a sharp discontinuity between the fuel and reaction products (ashes). The reactions occur instantaneously at the shock as the compressed front propagates through the fuel, leaving completely burned ashes behind. This model yields useful information such as the detonation velocity and thermodynamic values behind the front. Conservation equations for mass, momentum, and energy are used to determine the post-shock conditions from the pre-shock (fuel) conditions. These equations are, respectively,

\begin{array}{}\text{(1)}& {\rho }_{1}D={\rho }_2(D-u_2)\end{array}

\begin{array}{}\text{(2)}& P_2-P_1={\rho }_1{u}_{2}D\end{array}

\begin{array}{}\text{(3)}& E_2-E_1=P_1{V}_1+{\displaystyle \frac{1}{2}}{D}^2-P_2{V}_2-{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\rho }_1^2}{{\rho }_2^2}}{D}^2+q\end{array}

where {\rho }_i is the density, u_i is the fluid velocity, P_i is the fluid pressure, E_i is the internal energy, V_i is the specific volume V=\frac{1}{\rho } and i=1 represents the pre-shock value and i=2 represents the post-shock. The values D and q are the detonation velocity and the energy released during the burning, respectively. In these equations it is assumed that the fuel is at rest, u_1=0 [11, 12].

Combining Eqs. (1) and (2) to eliminate the fluid velocity u₂ yields the Rayleigh line

\begin{array}{}\text{(4)}& R={\rho }_1^2{D}^2-{\displaystyle \frac{P_2-P_1}{V_1-V_2}}=0\end{array}

which relates the pre-shock and post-shock fluid states. Combining (1), (2), and (3) to eliminate both the fluid velocity u₂ and the detonation speed D gives the Hugoniot curve

\begin{array}{}\text{(5)}& H=E_2-E_1-q-{\displaystyle \frac{1}{2}}(P_1+P_2)(V_1-V_2)=0\end{array}

which represents all possible post-shock states for the given initial conditions and value of q.

Plotting the Rayleigh line and Hugoniot curve in PV-space, the post-shock conditions may be analyzed. The intersection of these two curves gives the point(s) where the conservation equations (1)–(3) are satisfied, defining the post-shock state of the system. There are three general classes of solutions differing by the slope of the Rayleigh line (see Fig. 1), which is defined by the detonation velocity, D. In the first case, D is small enough that the slope of the Rayleigh line does not allow intersection of the two curves, so there is no stable detonation possible under this system. In the second case, there is a single intersection point at which the Rayleigh line is tangent to the Hugoniot curve. This defines the stable detonation condition for the system referred to as the CJ-velocity, D_{\mathrm{CJ}}. In the third case, the detonation velocity is greater than the stable CJ-velocity, and the curves intersect in two places. The upper and lower points are denoted as strong and weak detonations, respectively. Only strong detonations are believed to be physical [12], and they are unstable to fluctuations behind the shock [11]. Which class of solutions is found depends entirely on the choice of detonation speed used for the Rayleigh line.

2.2 Zeldovitch, von Neumonn, and Döring

In the middle of the 20\mathrm{t}\mathrm{h} century Zeldovitch [13], von Neumann [14], and Döring [15] developed the next advance in detonation theory. In the Zeldovitch, von Neumann, Döring (ZND) theory the shock front is also modeled as a discontinuity, but the reactions do not take place instantaneously. Instead the reactions are initiated by the shock and then proceed according to the proper reaction rates, resulting in a finite reaction zone behind the shock. The amount of burning which has occurred can be tracked by replacing the quantity q in Eq. (5) by the quantity \lambda q, where λ is a dimensionless parameter denoting the fraction of energy released, 0⩽\lambda ⩽1. This allows the width of the reaction zone to be determined and makes it possible to track thermodynamic values within the reaction zone.

2.3 Fluid dynamics

The simple detonation models discussed in Sections 2.1 and 2.2 are limited to one-dimensional symmetry and handle the nuclear energy release in one-step calculations where all the energy from the burned fuel is deposited into the detonation region. While their solutions serve as a gauge for classifying detonation behavior, they lack the complexity needed for realistic detonation models. Real detonations are multi-dimensional and time-dependent with feedback occurring between the nuclear energy generation and the hydrodynamics. More realistic models of astrophysical detonations employ the time-dependent, compressible, reactive flow Euler equations of fluid dynamics:

\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}\rho }{\mathrm{d}t}}+\nabla \cdot \rho v=0\end{array}

\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}\rho v}{\mathrm{d}t}}+\nabla \cdot (\rho vv)+\nabla P=\rho g\end{array}

\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{d}\rho E}{\mathrm{d}t}}+\nabla \cdot [(\rho E+P)v]=\rho v\cdot g\end{array}

\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}\rho X_i}{\mathrm{d}t}}+\nabla \cdot \rho X_{i}v=0\end{array}

where ρ is fluid density, v is the fluid velocity, P is the pressure, g is the gravitational acceleration. X_i is the mass fraction of the ith nuclear species such that the sum of all X_i is one. E is the sum of the internal energy, e, and kinetic energy per unit mass such that

\begin{array}{}\text{(10)}& E=e+{\displaystyle \frac{1}{2}}{|v|}^2\end{array}

For astrophysical detonations, these equations are closed by an equation of state capable of describing the pressure changes in degenerate matter, for example see Timmes and Swesty [16], and a scheme for tracking the evolution of the composition and the nuclear energy released. For a review and derivation of the Euler equations see Laudau and Lifhshitz [17], Arnett [18], and Fryxell et al. [19].

There are many historical numerical methods for solving these equations. For a detonation it is important to use a scheme that can adequately resolve the steep gradients in temperature, pressure, density, and composition associated with the shock front. One of the key challenges is finding a solution that will remain stable without introducing a large amount of numerical diffusion across steep gradients. Two methods that meet this challenge are the flux-corrected transport methods (FTC), see for example Oran and Boris [20] and Boris et al. [21], and the piecewise parabolic method (PPM) of Colella and Woodward [22]. Both methods, FTC a finite differencing scheme and PPM a finite volume Godunov scheme, correct for the effects of numerical diffusion by imposing monotonicity constraints. To achieve this, the FTC methods moderate hydrodynamic fluxes using a positivity preserving algorithm and the PPM solves a Riemann problem at each grid cell interface. The low computational cost and good shock capturing ability of these methods, especially the PPM, make them well suited for use in astrophysical simulations, which are often challenged by the presence of strong shocks. The FTC method has a history of use in terrestrial detonation modeling as well. All the multi-dimensional DNS studies of Type Ia detonations that we discussed in Section 3.2 use variants of these two methods to solve the Euler equations.

2.4 Considerations for nuclear kinetics

For astrophysical detonations, the coupling between the hydrodynamics and the nuclear burning can be quite complex. For example, reactions behind the shock front can become endothermic, effectively causing portions of the reaction zone to fall out of sonic contact with the detonation front. This leaves some of the energy from the reactions unavailable to sustain the detonation, for example see Sharpe [23] and Fickett & Davis [12]. In this situation the detonation is said to be pathological. The equations of fluid dynamics must be coupled with a realistic nuclear reaction network to create a reasonable model for the speed and structure of the pathological detonation front. A careful handling of the nuclear kinetics is also required for models attempting to accurately capture the multi-dimensional nature of the detonation because hydrodynamic effects may lead to complex incomplete burning.

Nuclear reaction networks involve a set of cross-coupled differential equations that represent the different reactions for all the species, which must be solved simultaneously to evolve the composition. The reaction rates are determined by the temperature and density. There are hundreds of species and thousands of nuclear reactions involved in realistic nuclear kinetics, so modeling their nucleosynthesis is computationally expensive when coupled with hydrodynamics. For the nuclear energy release, most astrophysical hydrodynamics simulations use a parameterization based on the output of separate detailed nuclear kinetics models. This method neglects the feedback between the nuclear reactions and the hydrodynamics. The best multi-dimensional DNS studies of detonations to date use an in situ alpha network composed of a limited set of 13 to 17 nuclei that describe the reactive flow in terms of alpha particles (⁴He) and nuclei that can be composed of integer multiples of alpha particles. While alpha networks neglect many of the important reactions channels, they are capable of following the basic endothermic and exothermic reaction processes that are key to modeling a detonation, see Hix et al. [24] and Timmes et al. [25]. Currently there are a few nuclear reaction network schemes, see Hix et al. [26] and Guidry [27], under development that run with an efficiency approaching that of the alpha network but also include the needed reaction channels to give them an accuracy approaching that of the more complete networks.

Due to the matter of numerical resolution an important complication when modeling astrophysical detonations is whether or not to allow nuclear burning to occur within the shock. Conventional wisdom states that this should not be allowed based on the findings of Ref. [11]. They performed a study comparing several different Eulerian hydrodynamics methods; the Piecewise Parabolic Method (PPM) [22], Godunov’s Method [28], SADIE [29], Lax–Wendroff [30, 31], Donor Cell [32]. Tests were run in one-dimension to quantify how well each method was able to model a sharp hydrodynamic shock traveling down a tube both with and without shock burning. These were then compared to results from the Lagrangian version of the PPM method [22]. A Lagrangian code was chosen for comparison due to its ability to maintain a perfectly sharp discontinuity. In all tests it was found that the accuracy with which a method was able to propagate the shock improved with resolution, and at most resolutions the PPM method performed better than the others. With regard to the question of shock burning, the paper made several important findings. In an ideal hydrodynamics model, the shock is assumed to be an infinitely thin discontinuity that is completely governed by a set of conservation equations describing the jump conditions across the discontinuity. In such a model, the fuel would spend no time within the shock itself, and does not begin to react until it has passed through the discontinuity. In all of the Eulerian methods tested, the shock is spread across some number of zones which can cause problems if the length of time the fuel spends within the shock becomes comparable to the reaction time. This can then result in a significant amount of the burning taking place inside the shock, possibly causing numerical errors in the hydrodynamics. The errors can quickly grow to yield unphysical results such as a bulge in the shock profile (if there is enough energy production the front of the shock) that may propagate as a separate wave ahead of the real shock. Fryxell et al. [11] found that by increasing the resolution such that the width of the shock was small compared to the burning length, they could eliminate these errors in all of the models. It should be noted that the PPM method was tested using both a single reaction network consisting solely of ¹²C and ⁵⁶Ni and the 13 isotope alpha network [33], while the other methods were only tested using the single reaction network.

2.5 Realistic detonation models

Experiments with detonations in terrestrial gases exhibit a multi-dimensional nature [12, 34, 35]. These experiments show that an initially planar detonation front will become perturbed (or wrinkled) by transverse and longitudinal instabilities, creating a complex cellular pattern in its wake. These real detonations have a complicated reaction zone as opposed to the simple one-dimensional models discussed previously, leading to nonuniform energy release and ash composition. The one-dimensional CJ and ZND models are however still useful in these more realistic simulations. The post-shock hydrodynamic variables in multi-dimensional detonations can deviate appreciably from these one-dimensional models along the detonation front due to the formation of the cellular structure, but the detonation velocity only differs by approximately 1%–2%. Furthermore, the ZND conditions for a corresponding one-dimensional detonation can be used as a good approximation for the initial conditions in the more realistic case. However, Dominguez and Khokhlov [36] warns that the correct assumption of final ash composition must be used in order to obtain the closest approximation to the initial conditions.

A schematic of the cellular structure is illustrated in Fig. 2, showing the detonation front (red curves in the figure) at three times [12, 37]. Instead of a planar detonation front, the real front is composed of strong burning regions called mach stems and weaker burning regions called incident shocks separated by triple-points. As the detonation front propagates into the fuel, the stronger mach stems create over-reacted regions of ash while the weaker incident shocks create under-reacted regions, leaving a nonuniform distribution of ashes behind the front. The triple-points (shown as blue dots in Fig. 2) are “hot spots” of increased temperature and density which initiate small spherical detonations in the front. The high-temperature, high-pressure environments left behind become the strong mach stem regions of the detonation front as these small detonations expand. The transverse shock waves (shown in green in Fig. 2) that emanate from these miniature explosions propagate perpendicular to the detonation velocity. As they propagate, they burn through the lateral regions of under-reacted ash created behind the weak incident shocks. The mach stems weaken as these transverse waves carry energy away from them, and it is these regions that become the new incident shocks. When transverse waves from neighboring hot spots travel toward each other into a common incident shock, their eventual collision causes the temperature and density in this region to increase, giving rise to the next generation of triple-points. So a particular region of the detonation front is constantly evolving, alternating between mach stem and incident shock, as the triple-points move back and forth across the front. As the detonation front proceeds into the fuel, a time-dependent pattern evolves which is created by the paths of the triple points (represented as black lines in Fig. 2). These paths are high-pressure regions which trap pockets (or cells) of under-reacted or over-reacted ash as the front propagates, creating the cellular structure. This structure was discovered while performing detonation experiments in terrestrial gases [12]. In these “smoke foil” experiments, soot covered plates were placed in square tubes where detonations were allowed to propagate. Cellular patterns were observed on these plates, where the diamond-shaped cells were traced out by the removal of ash along the paths of the high-pressure triple-points.

3 Direct numerical simulations

3.1 Structure and instabilities

DNS studies are perhaps the best tool available for probing the detailed structure and stability of astrophysical detonations. As evident by the multi-dimensional structure of real detonations described in Section 2.5, multi-dimensional effects play an important role in the evolution and fate of the detonation. Much progress has been made toward developing realistic multi-dimensional numerical models that produce the gross properties observed in the chemical flame experiments, for example see Oran et al. [5], Shepherd [38], and Stewart & Kasimov [39]. These works are a good foundation for the development of realistic DNS studies of astrophysical detonations. For astrophysical environments only a handful of multi-dimensional studies have been published in the literature. We will give a brief overview of a few of the key one-dimensional DNS having motivated our focus on the multi-dimensional studies.

The question of the size of the burning region that follows the leading shock of the detonation is important for Type Ia supernova because when the density is low key nuclear burning length scales can become comparable or greater than the characteristic scale of an exploding WD (1\times {10}^8 cm). Khokhlov [40] was able to obtain complex descriptions of the structure of the one-dimensional detonation front by conducting DNS studies of detonations propagating in carbon–oxygen mixtures and in helium. Nuclear burning in carbon–oxygen fuel proceeds in stages: C-burning \to O-burning \to Si-burning. When a detonation occurs a shock wave propagates into the carbon-oxygen fuel, resulting in the sequence of burning stages. Khokhlov [40] explained that this creates three distinct burning regions within the reaction zone, with each type of burning corresponding to a different burning width: x_{\mathrm{C}}\ll x_{\mathrm{O}}\ll x_{\mathrm{Si}}. These widths increase with decreasing density, and can lead to incomplete burning when the width becomes comparable to the radius of the WD. Khokhlov [40] has shown that this incomplete burning occurs at \rho \sim 1\times {10}^7 g\cdotcm^-3 and \rho \sim 1\times {10}^6 g\cdotcm^-3 for Si and O, respectively, while C-burning is always thought to continue to completion. Khokhlov [40] used a full nuclear reaction network and considered only supported (overdriven) detonations. Sharpe [23] extended this work by considering self-sustaining pathological detonations, which are more like those expected to form in Type Ia supernovae. Sharpe [23] used an alpha nuclear reaction network, but demonstrated that the resulting structure of tests with supported detonations matched that of Khokhlov [40], which used a more complete nuclear reaction network. Sharpe [23] concluded that the thickness of self-sustaining detonations was a few times larger than that predicted by the supported detonations of Khokhlov [40]. As will be shown in Section 3.2, multi-dimensional effects are also known to increase the size of the reaction zone.

Also important to Type Ia supernovae is the question of detonation stability. Instabilities in the burning region behind the detonation front can generate large variations in the temperature which will be reflected in the composition of the nuclear products. Terrestrial experiments have shown that instabilities increase the size of the reaction zone and also make initiation of the detonation more difficult [41]. Khokhlov [42], explored the one-dimensional stability of detonations to induced longitudinal perturbations by observing the growth or decay of small amplitude perturbations in one-dimensional simulations of steady state detonations under a range of density conditions. These studies show that thermonuclear detonations traveling in degenerate carbon-oxygen mixtures were unstable to longitudinal perturbations for densities above 2\times {10}^7 g\cdotcm^-3. This was due to positive feedback between the hydrodynamic fluctuations and the ¹²C+¹²C reaction. The instability developed in a thin layer behind the detonation front. For the lower fuel density of 2\times {10}^7 g\cdotcm^-3, the layer of burning which generated the instability was on the order of 1 cm thick, and for the higher fuel density of 3\times {10}^9 g\cdotcm^-3 it was only {10}^{-2} cm thick. The study concluded that this small-scale instability would not impact the nucleosynthesis, however, it could change the conditions for detonation initiation and extinction.

3.2 Multi-dimensional DNS of detonations

The first numerical study to show the multi-dimensional structure of a detonation in degenerate carbon and oxygen was Boisseau et al. [43]. They developed a one- and two-dimensional Cartesian hydrodynamics code that solved the continuity equations of fluid dynamics using the flux-corrected transport algorithm of Boris et al. [21]. The code also used an equation of state that was suitable for degenerate matter and employed an alpha-network to follow the nuclear burning. They tested their setup relative to Khokhlov [40] and Khokhlov [42] with one-dimensional simulations at densities of 3\times {10}^7 g\cdotcm^-3 and 1\times {10}^7 g\cdotcm^-3. Indeed they found close agreement with those studies for the size of the carbon burning induction zone and the instability of the detonation to longitudinal perturbations. They used conditions derived from the one-dimensional studies as initial conditions for two-dimensional studies which featured a fuel density of 3\times {10}^7 g\cdotcm^-3. The size of their computational grid was designed to resolve only the carbon burning scale and was too small in extent to capture the oxygen and silicon burning structures. Transverse perturbations were added to the leading edge of the detonation to provoke the formation of transverse waves. Like real detonations discussed in Section 2.5, the detonation cells left pockets of incompletely burned material at higher densities and lower temperatures than the surrounding material. The structure extended well behind the front. The cell size predicted by their two-dimensional study was approximately three times longer than that carbon burning zone predicted by their own one-dimensional study and that of Khokhlov [40]. This was an indication that the carbon induction zone had been lengthened by the pattern of burning in the cell structure. Like Khokhlov [42], Boisseau et al. [43] reasoned that small increases in the size of the carbon induction zone would have little impact on the results of the detonations. However, if the same scale of increase was present for the silicon burning zone, the burning zone might be larger than the white dwarf or larger than the scale of the white dwarf’s density gradient. This could impact the completeness of the energy released and products produced.

Gamezo et al. [44] directly tested multi-dimensional effects on the predicted lengths of the carbon, silicon, and oxygen burning zones. They chose to test each burning scale separately over a range of densities from 1\times {10}^6 to 3\times {10}^7 g\cdotcm^-3. They used one-dimensional calculations of pathological detonations to find the reaction zone widths and time scales for carbon, oxygen, and silicon burning and chose the resolutions of analogous two-dimensional models accordingly. This ensured that the main cellular features of each burning scale could be resolved, however, it prevented the two-dimensional trials from completely capturing the interactions between the different burning scales. In each case the detonation was at first overdriven from the starting shock but then relaxed to a velocity near the CJ velocity and a few trials were run with overdriven detonations. All trials were able to produce cells for the carbon, silicon, and oxygen length scales over a range of densities.

Fluctuations in the nuclear species concentrations due to the burning cells were clearly seen in the reaction zones at all burning scales. Gamezo et al. [44] posits that these fluctuations are a result of two mechanisms. The first occurs when the large induction zone behind the weak part of the leading shock is cut off by the collision of two triple points, and gives an under-reacted pocket of material. The second mechanism is the formation of over-reacted regions behind the strong shock in the vicinity of the triple-point collisions. The over-reacted regions contain less of the lower mass number reactants and more of the higher mass number products than the surrounding material. The net effect of these concentration fluctuations is to increase the width of the cell averaged reaction lengths for carbon, oxygen, and silicon compared to that given by one-dimensional models. The authors compared the reaction lengths of their one- and two-dimensional studies at all scales for a density of 6.5\times {10}^6 g\cdotcm^-3 and found that the two-dimensional carbon and oxygen lengths were up to 1.6 times longer than their one-dimensional counterparts when averaged over many cells. The silicon length was no more than 1.3 times its one-dimensional counterpart. Gamezo et al. [44] also shows that the shock velocity within the cells varies greatly, ranging on average from 0.75 to 1.7 times the CJ velocity for the carbon and oxygen scales and between 0.85 and 1.3 times the CJ velocity for the silicon cells. Coupled with the changes in reaction zone width and species concentration, these velocity fluctuations may influence the spectra of Type Ia supernova.

Prior to this study modelers wondered if the cellular nature of the detonation could revive the pure detonation model by allowing enough pockets of less reacted fuel to persist within the silicon cells to account for the quantity of intermediate-mass elements seen in Type Ia spectra. Gamezo et al. [44] shows that for the higher densities characteristic of pure detonation models, the silicon cells probably would not be large enough to account for the required quantities of intermediate mass elements. While it does not impact this conclusion, Gamezo et al. [44] cautions that the detonation cell sizes resulting from their calculations were still resolution-dependent, and 2 to 3 times larger than the “true” cell sizes expected for higher numerical resolution.

A study performed by Timmes et al. [37] differs with Gamezo et al. [44] on the issue of the resolution dependence of the flow features in the multi-dimensional burning. Timmes et al. [37] utilized the FLASH hydrodynamics code which solves the problem of adequate resolution with Adaptive Mesh Refinement (AMR) [19]. AMR, which employs different adaptive levels of resolution on the same grid, potentially allows all the detonation burning scales to be considered at once. However, Timmes et al. studied only the carbon scale [37]. Spatial resolutions of 1\times {10}^{-1}, 5\times {10}^{-2}, 2.5\times {10}^{-2} and 1.25\times {10}^{-2} cm were examined for detonations propagating near the CJ velocity in pure ¹²C fuel at a density 1\times {10}^{-7} g\cdotcm^-3. Timmes found that the cell size did not change with increasing resolution. Timmes et al. [37] also found that the cells on average have a longer length than width, whereas Gamezo et al. found a more even width to length ratio [44]. Since Timmes et al. was using a pure carbon fuel [37], the difference may be rooted in the initial energy source for the burning, but further studies are needed to verify this. Timmes et al. [37] helped to define the minimum resolution required by a multi-dimensional DNS where cellular detonation is a key feature of the model.

Although the true structure of a detonation is three-dimensional there are few DNS studies of detonations performed in three-dimensions due to the computational cost being at least an order of magnitude greater than that of a two-dimensional simulations. Calder et al. [45] ran the first three-dimensional DNS of astrophysical detonations which were published in a methods paper that won the 2000 Gordon Bell Prize. The test problems were the direct three-dimensional analog of the 0.1 cm resolution two-dimensional DNS shown in Ref. [37]. The paper featured the method of using the FLASH code’s AMR to make the problem computationally efficient and able to handle all the scales required for exploring the three-dimensional carbon burning cell structure. Calder et al. [45] states that AMR allowed them a savings of a factor of forty in the number of grid points needed to run the three-dimensional simulation compared to what would have been required by a fixed grid. Three-dimensional DNS studies of terrestrial detonations, such as Williams et al. [46], show that the fine structure within cells and the vorticity at flame cell borders are influenced by the extra dimension but the spacing of the transverse waves which create the burning cell are very similar to that produced by two-dimensional models. Gamezo et al. [44] interprets this as meaning that the cell sizes and the width of the reaction zones produced by two- and three-dimensional studies may be very similar in size.

Parete-Koon et al. [47] carried out a small number of two- and three-dimensional DNS studies of detonations at a density of 5\times {10}^7 g\cdotcm^-3 using the FLASH code to test the ability of detonations to survive passage through thin funnels of fuel separated by neighboring regions of ash. The study was motivated by the idea that the Type Ia detonation may follow a deflagration phase that would leave the detonation’s fuel riddled with regions of ash. They also explored the results of Maier and Niemeyer [48], which in a similar set of FLASH based one- and two-dimensional studies found that the detonation was quenched when passing through fuel funnels of less that 16 cm. However, limited computational resources restricted the grid resolution of the Maier and Niemeyer [48] two-dimensional studies to 0.5 cm, which was not sufficient to resolve the carbon cellular structure detonation. Parete-Koon et al. [47] concurred that the detonation could survive a fuel funnel of no less than 16 cm in width when the simulation had a maximum grid resolution of 0.5 cm. However, with a maximum grid resolution greater 0.5 cm, the detonation was quenched in its passage through the funnel for widths as great as 20 cm. They also tested the problem at 0.125 cm grid resolution in a three-dimensional tube with semi-cylinders of ash used to constrict the fuel. They found that the results were similar to the analogous two-dimensional tests in terms of detonation survival; however, there were significant differences in the smaller scale structure of the burning front and its ash. While neither of these studies directly explored the impact of the cellular structure, Parete-Koon et al. [47] shows in Fig. 3 that it is present in their studies with grid resolutions of 0.125 cm and greater.

While the fine detail of the carbon cells produced by a realistic detonation is on a scale that is several orders of magnitude below that of the finest resolution of the current full star Type Ia explosion models, the silicon cells and in lower density regions even the oxygen cells may be resolvable. Meakin et al. [49] developed a full Type Ia model using the FLASH code that featured a deflagration followed by a parameterized single point detonation. The nuclear reactions in the study were parameterized. According to the cell sizes observed by Gamezo et al. [44], the resolution of their grid was capable of capturing the silicon cells of the detonation. However, they state that it was difficult to discern if the cellular structure had formed because the layer behind the detonation wave that would have contained it quickly mixed with the turbulent layer of deflagration ash immediately above it. Models with finer resolution or lower density conditions would be needed to determine the possible impact of the cells.

4 Conclusions

The Type Ia modeling community has only begun to exploit the information gained from the current DNS studies of detonations. These studies illustrate that the burning in the detonation is quite complex and requires a careful treatment especially on scales that are comparable to the density gradient of the exploding white dwarf. Considerable effort has been expended modeling supernovae at higher densities where the complexities of turbulent burning dominate the model, for example see Khokhlov [50] and Zingale et al. [51]. However, most full star models turn the nuclear burning schemes off when the density falls below 1\times {10}^7 g\cdotcm^-3 and distributed burning begins. This is especially problematic given that the DDT is believed to occur at just these densities. A recent notable exception are the two-dimensional simulations of Jackson et al. [52] and Townsley et al. [53] where the front tracking algorithm in the FLASH code has been extended to facilitate tracking both deflagrations and detonations. Though this method was calibrated using DNS, its inclusion in the front tracking algorithm is somewhat heuristic. DNS studies of detonations can provide a good test bed for detonation stability at the key densities where DDT is believed to occur and they serve as a motivator for careful treatment of the nuclear burning in full star models.

There is still much work to do in developing DNS studies of detonations. The question of how a complete nuclear reaction network will impact the predicted structure of the multi-dimensional detonation remains open. A complete consistent systematic study of the effect of dimensionality on the structure of the detonation has not been published in the astrophysical literature. The question of how the cell sizes are related to the size of the burning zone has only begun to be addressed. Multi-dimensional studies with all burning scales represented on the same simulation grid have not been done. The knowledge gained from existing DNS studies provides a good foundation for exploring these questions. The publicly available simulation codes, such as FLASH, coupled with the currently available computational resources, leave modelers in a good position for making progress.

References

[1] W. Arnett, Astrophysics and Space Sciences, 1969, 5: 180

[2] S. Woosley and T. Weaver, Ann. Rev. Astron. Astrophys., 1986, 24: 205

[3] A. Filippenko, Ann. Rev. Astron. Astrophys., 1997, 35: 309

[4] D. Kasen, F. Röpke, and S. Woosley, Nature, 2009, 460: 869

[5] E. S. Oran, V. N. Gamezo, and D. A. Kessler, Tech. Rep. RL/MR/6400–11-9332, Naval Research Laboratory, 2011

[6] A. Y. Poludnenko, T. A. Gardiner, and E. S. Oran, Phys. Rev. Lett., 2011, 107(5): 054501

[7] F. Timmes and S. Woosley, Astrophys. J., 1992, 396: 649

[8] A. Calder, D. Townsley, I. Seitenzahl, F. Peng, O. Messer, N. Vladimirova, E. Brown, J. Truran, and D. Lamb, Astrophys. J., 2007, 656: 313

[9] D. Chapman, Philosophical Magazine, 1899, 47: 90

[10] J. Jouguet, J. Math. Pure Appl., 1905, 1: 347

[11] B. Fryxell, E. A. Muller, and D. Arnett, Max Plank Institute for Astrophysics, Pre-print, 1989: 449

[12] W. Fickett and W. Davis, Detonation: Theory and Experiment, Courier Dover, 1979

[13] Y. Zeldovitch, Zh. Eksp. Teor. Fiz., 1940, 10: 524

[14] J. von Neumann, OSRD Reports, 1942: 549: 1

[15] W. Döring, Annalen der Physik, 1943: 435

[16] F. X. Timmes and F. D. Swesty, Astrophys. J. Suppl., 2000, 126: 501

[17] L. D. Landau and E. M. Lifhshitz, Fluid Mechanics, Pergamon Press, 1959

[18] D. Arnett, Supernovae and Nucleosynthesis: An Investigation of the History of Matter, from the Big Bang to the Present, Princeton: Princeton University Press, 1996

[19] B. Fryxell, K. Olsen, P. Ricker, F. Timmes, M. Zingale, D. Lamb, P. MacNeice, R. Rosner, J. Truran, and H. Tufo, Astrophys. J. Suppl., 2000, 131: 273

[20] E. S. Oran and J. P. Boris, Numerical Simulation of Reactive Flow, Cambridge: Cambridge University Press, 1987

[21] J. P. Boris, A. M. Landsberg, E. S. Oran, and J. H. Gardner, LCPFCT-A Flux-Corrected Transport Algorithm for Solving Generalized Continuity Equations, Tech. Rep., Naval Research Lab, 1993

[22] P. Colella and P. Woodward, J. Comput. Phys., 1984, 54: 174

[23] G. Sharpe, Mon. Not. R. Astron. Soc., 1999, 310: 1039

[24] W. Hix, A. M. Khokhlov, J. C. Wheeler, and F. Thielemann, Astrophys. J., 1998, 503: 332

[25] F. X. Timmes, R. D. Hoffman, and S. E. Woosley, Astrophys. J. Suppl., 2000a, 129: 377

[26] W. R. Hix, S. T. Parete-Koon, C. Freiburghaus, and F.-K. Thielemann, Astrophys. J., 2007, 667: 476

[27] M. W. Guidry, J. Comput. Phys., 1012, 232: 5266

[28] S. Godunov, Matematicheskii Sbornik, 1959, 47: 165

[29] C. Arnold, Ph.D. thesis, University of Michigan, 1985

[30] P. Lax and B. Wendroff, Communications on Pure and Applied Mathematics, 1960, 13: 217

[31] P. Lax and B. Wendroff, Communications on Pure and Applied Mathematics, 1964, 17: 381

[32] R. A. Gentry, R. E. Martin, and B. J. Daley, J. Comput. Phys., 1966, 1: 87

[33] F. Thielemann, 1985 (private communication)

[34] Y. N. Denisav and Y. K. Troshin, Dokl. Akad. Nauk SSSR (Phys.-Chem. Sec.), 1959, 125: 110

[35] B. V. Voitsekhovsky, V. V. Mitrofanov, and M. E. Topchian, Izd. Akad. Nauk SSSR, 1963

[36] I. Dominguez and A. Khokhlov, Astrophys. J., 2011, 730: 87

[37] F. Timmes, R. Hoffman, and S. Woosley, Astrophys. J., 2000b, 129: 377

[38] J. E. Shepherd, in: Proceedings of the Combustion Institute, 2009

[39] D. Stewart and A. Kasimov, Journal of Propulsion and Power, 2006, 22: 1230

[40] A. Khokhlov, Mon. Not. R. Astron. Soc., 1989, 239: 785

[41] J. H. Lee, Annual Review of Fluid Mechanics, 1984, 16: 311

[42] A. M. Khokhlov, Astrophys. J., 1993, 419: 200

[43] J. Boisseau, J. Wheeler, E. Oran, and A. Khokhlov, Astrophys. J., 1996, 471: L99

[44] V. N. Gamezo, J. Wheeler, A. Khokhlov, and E. Oran, Astrophys. J., 1999, 512: 827

[45] A. C. Calder, B. C. Curtis, L. J. Dursi, B. Fryxell, G. Henry, P. MacNeice, K. Olson, P. Ricker, R. Rosner, F. X. Timmes, et al., in: Super Computing Gordon Bell Prize Paper, 2000

[46] D. N. Williams, L. Bauwens, and E. S. Oran, Symposium on Combustion and Flame, and Explosion Phenomena, 1996, 26: 2991

[47] S. Parete-Koon, C. Smith, M. Guidry, R. Hix, and O. Messer, J. Phys., 2012 (in press)

[48] A. Maier and J. Niemeyer, Astron. Astrophys., 2006, 451: 207

[49] C. A. Meakin, I. Sietenzahl, D. Townsley, G. C. J. IV, J. Truran, and D. Lamb, Astrophys. J., 2009, 693: 1188

[50] A. M. Khokhlov, Astrophys. J., 1995, 449: 695

[51] M. Zingale, S. Woosley, C. A. Rendleman, M. Day, and J. B. Bell, Astrophys. J., 2005, 632: 1021

[52] A. P. Jackson, A. C. Calder, D. M. Townsley, D. A. Chamulak, E. F. Brown, and F. X. Timmes, Astrophys. J., 2010, 720: 99

[53] D. M. Townsley, A. P. Jackson, A. C. Calder, D. A. Chamulak, E. F. Brown, and F. X. Timmes, Astrophys. J., 2009, 701: 1582, 1604

Rights & permissions

Higher Education Press and Springer-Verlag Berlin Heidelberg 2013