Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany

Abstract
Type Ia supernovae are believed to be white dwarfs disrupted by a thermonuclear explosion. Here we investigate the scenario in which a rather low-mass, carbon-oxygen (C + O) white dwarf accumulates helium on its surface in a sufficient amount for igniting a detonation in the helium shell before the Chandrasekhar mass is reached. In principle, this can happen on white dwarfs accreting from a non-degenerate companion or by merging a C + O white dwarf with a low-mass helium one. In this scenario, the helium detonation is thought to trigger a secondary detonation in the C + O core. It is therefore called the "double-detonation sub-Chandrasekhar'' supernova model. By means of a set of numerical simulations, we investigate the robustness of this explosion mechanism for generic 1- $\ensuremath{M_\odot}$ models and analyze its observable predictions. Also a resolution dependence in numerical simulations is analyzed. Hydrodynamic simulations of the double-detonation sub-Chandrasekhar scenario are conducted in two and three spatial dimensions. The propagation of thermonuclear detonation fronts, both in helium and in the carbon-oxygen mixture, is computed by means of both a level-set function and a simplified description for nuclear reactions. The decision whether a secondary detonation is triggered in the white dwarf's core or not is made based on criteria given in the literature. In a parameter study involving different initial flame geometries for He-shell masses of 0.2 and 0.1 $\ensuremath{M_\odot}$ (and thus 0.8 and 0.9 $\ensuremath{M_\odot}$ of C + O), we find that a secondary detonation ignition is a very robust process. Converging shock waves originating from the detonation in the He shell generate the conditions for a detonation near the center of the white dwarf in most of the cases considered. Finally, we follow the complete evolution of three selected models with 0.2 $\ensuremath{M_\odot}$ of He through the C/O-detonation phase and obtain $\ensuremath{^{56}{\rm Ni}}$ -masses of about 0.40 to 0.45 $\ensuremath{M_\odot}$ . Although we have not done a complete scan of the possible parameter space, our results show that sub-Chandrasekhar models are not good candidates for normal or sub-luminous type Ia supernovae. The chemical composition of the ejecta features significant amounts of $\ensuremath{^{56}{\rm Ni}}$ in the outer layers at high expansion velocities, which is inconsistent with near-maximum spectra.

Key words: supernovae: general - nuclear reactions, nucleosynthesis, abundances - hydrodynamics - methods: numerical

1 Introduction

Therefore, a long standing question is whether other progenitor channels contribute to the observed SN Ia sample, and two alternatives have been suggested: the WD merger or double degenerate scenario (in contrast to single degenerate models, which consist of binaries with only one WD) and the sub-Chandrasekhar model. In the present work we explore the latter.

The basic idea of the sub-Chandrasekhar model is that if the accretion rate onto a WD is lower than about 1-4 $\times$ $10^{-8}~\ensuremath{M_\odot} ~{\rm yr}^{-1}$ the accreted He (or the H processed to He) is not fused steadily into C and O. Instead, after reaching a critical amount of He at relatively low densities, the He shell becomes unstable and detonates (cf. Woosley et al. 1986, and references therein). The detonation ignites most likely close to the bottom of the He shell, produces almost pure $\ensuremath{^{56}{\rm Ni}}$ , and can occur long before the WD reaches the Chandrasekhar limit. A second detonation may be triggered spontaneously when the He shell detonation shock wave hits the C + O core or with some delay after the shock has converged near the center of the WD (double detonation models). This way, in the sub-Chandrasekhar mass models the conditions for thermonuclear runaway are caused by the compression in the shock wave and not by the high degree of degeneracy as in the Chandrasekhar mass models.

Delayed double detonations have been studied extensively before. Woosley & Weaver (1994), and Livne (1990) carried out one-dimensional (1D) simulations and Livne & Glasner (1991,1990) considered two-dimensional (2D) setups. In 1D the He detonation ignites synchronously in a layer close to the bottom of the He shell and due to the spherical symmetry constraint the shock converges perfectly in the center. However, this is not a very realistic model. According to Livne (1990) an ignition most likely happens in a single point leading to an off-center convergence of oblique shock waves in the core.

Several other simulations predict successful directly ignited double detonations in two dimensions (Arnett 1997; Wiggins et al. 1998; Livne 1997; Wiggins & Falle 1997). The smoothed particle hydrodynamics simulations by Benz (1997) and García-Senz et al. (1999) were carried out in 3D. Livne & Glasner (1990) first reported an increased probability of the direct core ignition if the He detonation does not happen directly at the core-shell interface but at a certain distance above it. This way, the pressure jump can grow large enough before hitting the core.

In this work, the results of a (restricted) parameter study are presented investigating the possibility of triggering the second detonation in two and three dimensions. In all our models the WD has a total mass of 1 $\ensuremath{M_\odot}$ , but the (C + O)-core and the He-shell masses differ. We compute sequences of models with very different initial flame geometries in order to test the robustness of this explosion mechanism. Since this latter question is the main focus of our paper, most of the simulations are stopped once the conditions for a detonation are matched. However, for a few successful cases the energetics and nucleosynthesis of complete double detonations are computed. In the following chapter the model details will be described. Section 3 shows the simulation results, which are summarized and discussed in the last part of the paper.

2 Explosion model

2.1 General setup

Our hydrodynamic simulations were carried out on a uniform Eulerian grid with a finite volume scheme based on the PROMETHEUS implementation (Fryxell et al. 1989) of the "piecewise parabolic method'' (PPM) (Colella & Woodward 1984). Thus, discontinuities like shocks were automatically and accurately captured but, nonetheless, smeared out over several grid cells. A potential concern is the overheating effect observed in Godunov-type schemes when shocks hit walls or collide. Following the suggestion of Donat & Marquina (1996), an approximate Riemann solver (the so-called Marquina solver) was implemented but this did not lead to noticeable differences. Thus, all simulations presented in the following are based on the exact Riemann solver of Colella & Glaz (1985).

While the parameter study testing the possibility of igniting a detonation in the C + O core was carried out on a static grid, the complete double detonation models presented at the end of Sect. 3 employed a co-expanding grid as in Röpke & Hillebrandt (2005) in order to account for the expansion of the WD. The coordinates were (r, z)-cylindrical in 2D and Cartesian in 3D.

The equation of state we used includes contributions from an arbitrarily degenerate and arbitrarily relativistic electron gas, a photon gas, electron-positron pairs and nuclei with a Maxwell-Boltzmann distribution. The initial WD was assumed to be in hydrostatic equilibrium. Therefore its structure is uniquely determined once the density at the center and at the core-shell interface and the temperature profile are fixed. For simplicity a constant temperature of 5 $\times$ $10^5~{\rm K}$ was assumed throughout the star. Of course, this is an unrealistic choice, but because of the high degree of degeneracy the exact temperature has almost no effect on the state of the matter. Another simplifying assumption is the approximation of the gravitational potential in spherical symmetry. Only the monopole moment of the density distribution averaged over all angles was calculated at every time step.

2.2 Nuclear reactions

As we are mostly interested in the hydrodynamic evolution of the explosion, it is appropriate to use a strongly simplified scheme for nuclear reactions as in Reinecke et al. (2002a). On our numerical grid the structure of the flame cannot be resolved. Therefore in our simulations the detonation flame is approximated as a discontinuity propagating at an appropriate speed (see Sect. 2.3) with the fast reactions happening spontaneously at its passage and only the density dependent products, but no reaction rates being calculated. In order to approximate the energy release in the fast reactions and the correct molecular weight of fuel and ashes, it is sufficient to consider only one representative species for each group of elements involved which has a binding energy and mass typical for the group. Our scheme includes four species, namely $\ensuremath{^{4}{\rm He}}$ , a mixture of $\ensuremath{^{12}{\rm C}}$ and $\ensuremath{^{16}{\rm O}}$ with equal mass fractions, $\ensuremath{^{28}{\rm Si}}$ as a representative for the intermediate mass elements, and $\ensuremath{^{56}{\rm Ni}}$ as an iron group element. The following reactions are included in the model, depending on the density at the front:

In order to distinguish the $\ensuremath{^{4}{\rm He}}$ in the shell from that produced in NSE with $\ensuremath{^{56}{\rm Ni}}$ both are treated as different species in the code. With decreasing densities the $\ensuremath{^{4}{\rm He}}$ in NSE finally recombines to $\ensuremath{^{56}{\rm Ni}}$ . This is the only slow reaction included in our simulations.

The propagation of the detonation flame is modeled with the "level set technique'' (Osher & Sethian 1988) in the implementation of Golombek & Niemeyer (2005), which in turn is based on the deflagration flame model of Reinecke et al. (1999b). The zero level set of a signed distance function G represents the flame front. Information on whether matter has already been burnt is encoded in its sign. It is negative in the fuel and positive in the ashes. Most important, the location of the flame can be calculated at every time step simply by finding the roots of G. The propagation of the zero level set of G consists of two parts: the advection with the flow and self-propagation by burning, respectively. By means of an operator splitting, the advection of G can be treated in the same way as that of a passive scalar in the hydrodynamic scheme (Reinecke et al. 1999b; Mulder et al. 1992; Reinecke et al. 1999a). Burning is accounted for by (Reinecke 2001)

$\begin{displaymath}\left( \frac{\partial G}{\partial t} \right)_{\rm burn} = s \left\vert \vec{\nabla} G \right\vert. \end{displaymath}$

The changes in the mass fractions are calculated depending on $\rho$ , when the flame front passes a cell. Here and also in the case of the velocity s we have the problem that in general the state in such a cell is a mixture of a burnt and an unburnt matter. So we have e.g.:

$\begin{displaymath}\rho = \alpha \rho_1 + (1 - \alpha) \rho_2, \end{displaymath}$

The scheme described here is not in general appropriate to simulate detonations as, e.g., the burning speed s depends on the detailed geometry of the detonation front. The burning speed s is in general a function of $\frac{\delta}{R}$ , where $\delta$ is the width of the detonation front and R is its curvature radius. In our case, however, the front can be approximated as thin and planar in the density range that is relevant for the explosion dynamics.

2.3 Detonation physics

A detonation wave consists of a shock wave that compresses and heats the matter behind it so strongly that it is ignited. The shock wave, in turn, is maintained by the energy released due to burning. Depending on the properties of the reactions taking place, there exist different types of detonations. A very important case is the Chapman-Jouguet detonation, which propagates into the fuel at the lowest possible speed. According to Landau & Lifschitz (1991) a detonation that is created spontaneously by the combustion itself is a Chapman-Jouguet detonation in many cases. Therefore, in our simulations we assumed this type of detonation for the He burning. In the Chapman-Jouguet case there exists a simple law for the velocity of the burning front relative to the ashes, which is just equal to the sound speed there:

$\begin{displaymath}s_{\rm 2,He} = \sqrt{\gamma_2 \left( \frac{\partial p_2}{\partial \rho_2} \right)_T}\cdot \end{displaymath}$

According to Sharpe (1999), self supporting detonations in degenerate carbon/oxygen mixtures at densities between 2 $\times$ 10⁷and 1 $\times$ $10^9~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ are of pathological type. These arise in general from a non-monotonic heat release in the reactions. For them the front velocity cannot be determined independent of the reaction rates. Therefore, for the C/O detonations in this work, we used the values calculated by Sharpe (1999) in numerical simulations as a function of density:

$\begin{displaymath}s_{\rm 1,C/O} = f_{\rm Sharpe}(\rho)\footnotemark. \end{displaymath}$

2.4 Ignition conditions in the C + O core

Table 2: Critical temperatures for detonations at given density and mass according to Röpke et al. (2007b). Detonation is impossible for $\rho \le 1$ $\times$ $10^6~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ .

$\begin{figure} \par\includegraphics[width=11.5cm,clip]{8438fig1} \end{figure}$	Figure 1: Initial flame models. The He detonation flame fronts start at the borders of the black regions that already contain hot ashes at the beginning.
Open with DEXTER

As our code lacks a detailed nuclear network and reaction rates are not calculated, the onset of a C/O detonation can not be followed. Apart from that, in our full star simulations it would not be feasible to resolve the scales of the detonator (where the shock that supports the detonation is built up). Therefore, the decision on whether a detonation is triggered is made upon critical conditions given by Niemeyer & Woosley (1997) and Röpke et al. (2007b) (see Tables 1 and 2). They are based on the following arguments. If there is a sufficiently large region of supersonic spontaneous burning, a strong pressure wave must form and will synchronize with the burning wave if the compression and heating is strong enough to burn the material in the sound crossing time of the region under consideration (Khokhlov 1991; Blinnikov & Khokhlov 1986; Woosley 1990). Therefore a detonation wave is born if that region is big enough or, equivalently, provided that it consists of enough mass at a given density. Niemeyer & Woosley (1997) empirically determined those critical masses from hydrodynamic simulations for different densities and for a given linear radial temperature gradient across that mass. Röpke et al. (2007b) used the same code and setup, but determined the critical temperature at a given mass for different densities. They focussed on the low end of possible densities.

3 Simulations

With the numerical scheme described above, several simulations were performed. In Sect. 3.1, the results of a parameter study on the influence of the initial flame geometry for two different core and He-shell masses are presented, followed by a resolution study in Sect. 3.2 and, finally, the determination of explosion energies and nucleosynthesis yields in three successful double detonation simulations. We stress again that the main intention of the parameter study is not to scan the entire parameter space of possible WD properties but rather to investigate the conditions which may lead to a double detonation.

3.1 Initial flame geometry

The accretion of matter from the companion star onto the WD is a very complex process, and the formation and ignition of a He shell is still afflicted with large uncertainties. As many different initial flame geometries in the shell seem possible, their influence on the possibility and the properties of a double detonation is explored first.

The different initial flame models, which were treated in this study, are shown in Fig. 1 and the nomenclature of the simulations is as follows: the first letter in each model's name denotes the geometry, and the following number specifies it further. "s4.10'' means an ignition in a whole layer of the He shell at a radius of 4100 km. "bn'' stands for an ignition in n separate "bubbles'' (which are, of course, tori in 2D if located off the z-axis). Two special cases of the latter type, for n = 1, are "z5.00'' and "r5.00'', where the initial "bubble'' is located at the z- or the r-axis, respectively, at a distance of 5000 km from the origin. Finally the capital letters "A'' or "B'', which follow the above names, indicate the initial masses which were used:

Table 3: Results of the parameter study about the initial flame geometry. The last column states if the ignition conditions of a core detonation have been reached.

The most symmetric case with the simplest hydrodynamic flow pattern is the ignition in a whole layer of the He shell. In model s4.10A_2dq_256 starting with the initial flame geometry shown at the top left of Fig. 1, the inner and the outer spherical detonation front start to propagate inwards and outwards, respectively. While the combustion stops close to the border of the WD due to the low densities, it ceases at the core-shell interface because of the depletion of the He fuel. To investigate the chances of a core ignition, only the hydrodynamic evolution was considered afterwards and all the reactions taking place were stopped. The shock wave that was formerly driving the inner detonation front continues to propagate inwards until it converges at the center. There (at $t = 0.575~{\rm s}$ ) a strong shock collision takes place that causes a maximum temperature of 10.2 $\times$ 10⁹ K and densities up to 7.37 $\times$ $10^8~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ .

According to Table 1 these maximum values suffice for a C/O ignition, if there is at least a mass of about 2.0 $\times$ $10^{15}~{\rm g}$ with the same temperature and density or, more accurately, with a sufficiently small temperature decrease from the above mentioned value. A mass of 2.0 $\times$ $10^{15}~{\rm g}$ at this density corresponds to a spherical volume with a radius of the order of 1 km. This length was not resolved in this simulation, where one cell has a diameter of $23.2~{\rm km}$ . But as the high temperatures and densities in our simulations are not confined to one single cell, but distributed over several (which means that the region of interest is resolved), the ignition conditions seem to be clearly fulfilled. Consequently, a C/O detonation will occur at the center of the core. There is, however, the possibility that a detonation of the C/O material is triggered already when the shock wave enters at the edge of the core. About 0.1 s after the start of the simulation, temperatures of about 1 to 2 $\times$ $10^9~{\rm K}$ and densities of about $10^7~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ are reached. According to Table 2, this is within the range of conditions where a detonation might occur. But as the temperature gradient might possibly be too large across a radius of about 30 km, this question remains open here. As was already mentioned, other authors (Arnett 1997; Wiggins et al. 1998; Benz 1997; García-Senz et al. 1999; Wiggins & Falle 1997; Livne 1997) found successful direct ignitions and an increasing probability for the latter with increasing distance of the He ignition point from the C/O core (cf. Livne & Glasner 1990).

Our simulations confirm the earlier results. Helium shell detonations were ignited at smaller (s3.60A_2dq_256) and larger (s5.00A_2dq_256) radii than in the case described above. As was expected, the shock wave propagating through the core is weaker in the case of smaller ignition radii and stronger in the case of shocks starting further out. Accordingly the maximum temperatures and densities that are reached in the center are lower in the former and higher in the latter case (see Table 3).

Due to the high degree of symmetry, the models with shell ignition are the ones that can trigger a second detonation most easily. To test multi-dimensional configurations, the complete ignition layer was replaced by single spots (which effectively are tori in 2D) distributed over the He shell. Although this is probably not a realistic scenario, it is considered here as a numerical experiment. As far as their maximum densities and temperatures are concerned, the models with five initially ignited bubbles in the first quadrant differ only slightly from the complete shell ignition model (see Table 3). The values of $T_{\rm max}$ and $\rho_{\rm max}$ are only slightly lower than for the ignition of a shell, starting at the same radius, and they still match the conditions for a detonation. If the number of initially burning bubbles is reduced to three, peak temperature and density decrease further.

$\begin{figure} \par\includegraphics[width=12cm,clip]{8438fig2} \end{figure}$	Figure 2: Time evolution of the model b3A_2dq_256. Here and in the following figures the density is color coded and lines of constant pressure are drawn in blue. In cyan the solid line is the location of the He detonation flame (the zero level of G) and the dashed lines are at the border of the He shell.
Open with DEXTER

The evolution of simulation b3A_2dq_256 is shown in Fig. 2. Due to the coalescence of the parts of the detonation fronts moving inwards, an almost spherically symmetric shock front is entering the core. This resemblance of the dominant hydrodynamic processes is the reason for the similarity of the results given in Table 3.

In between the bubbles interactions of the detonation shocks take place. The collision process of two radial shocks can be described as a reflection at the plane of symmetry between the waves (see e.g. Courant & Friedrichs 1948). Generally, in this frame the incidence angle may differ from the angle of reflection. This can be seen best at $t = 0.15~{\rm s}$ in Fig. 2. The reflected waves are again reflected at the waves coming from the opposite side of the initial bubble. The repetition of this process should also have some influence on the geometry of the innermost shock front and thus on the conditions given in Table 3.

This can be seen in model b2A_2dq_256, which is similar to model b3A_2dq_256, but with the bubble on the z-axis removed. It starts with only two initially ignited bubbles but reaches greater maximum densities and temperatures in the core. The reason may be a focusing effect on the z-axis in 2D cylindrical symmetry and the bigger density and pressure jump that can build up until the burning reaches the axis.

An analogous focusing effect in the plane of equatorial symmetry (the xy-plane) is not observed. On the contrary, removing the initial burning bubble from the r-axis (i.e. going from b2A_2dq_256 to b1A_2dq_256) leads to a strong decrease of the maximum values of density and temperature.

$\begin{figure} \par\includegraphics[width=8cm,clip]{8438fig3} \end{figure}$	Figure 3: Structure of a detonation propagating at the core-shell interface in the rest frame of the detonation front from Livne & Glasner (1990). (Here the core consists of carbon only, not of C/O like in our case.)
Open with DEXTER

Considering models with only a single initially burning bubble, the exact geometry of the appearing shock waves is of great importance. After the He detonation has reached the border of the core, it spreads around it in a wave propagating almost perpendicular to the core-shell interface. The detailed hydrodynamic structure has been analyzed in Livne & Glasner (1990), from which Fig. 3 is taken. An oblique shock wave that propagates into the core (tilted by an angle $\delta$ with respect to the interface) is associated with the detonation wave moving along the surface of the core. After the shock passage the core-shell interface is also tilted inwards by an angle $\beta$ . Thus in-flowing matter follows the oblique shock. The detonation wave is followed by a rarefaction wave. The angles $\beta$ and $\delta$ can be determined by applying the Rankine-Hugoniot conditions and depend on the strength of the detonation shock.

This structure can clearly be seen at $t = 0.40~{\rm s}$ in model r4.24A_2dq_256, which is the detonation of a torus in the z = 0 plane (see Fig. 4). Starting from the "equatorial'' plane the detonation wave propagates around the core until it reaches the "polar'' z-axis. There a collision of the waves coming from all azimuthal directions takes place. Equivalently, each wave is reflected at the wave coming from the opposite side of the star. This causes the shock wave visible in the last two snapshots of Fig. 4 to propagate to the lower right. In the following the oval inner shock front moves towards the center. At $t = 0.70~{\rm s}$ , even before the center is reached, the maximum density and temperature values are attained on the z-axis at about 850 km away from the center (fourth plot in Fig. 4). The density and temperature values of 1.54 $\times$ $10^8~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ and 5.62 $\times$ $10^9~{\rm K}$ are significantly lower than those in the spherically symmetric models, but according to Table 1 they still suffice for the ignition of an off-center core detonation.

In the following simulations, the symmetry of the first detonation was further reduced by departing from equatorial symmetry and simulating a whole rotationally symmetric star, i.e. now the negative z-range was considered too. Figure 5 shows the example simulation z4.24A_2d_256. Here the He detonation starts in a point on the z-axis and then spreads around the whole WD until it reaches the negative z-axis (see Fig. 5). In this sense, the setup is maximally asymmetric. The processes now taking place are analogous to the simulation r4.24A_2dq_256 described above. But here the maximum density and temperature are reached when the oval inner shock front converges almost perfectly in a "point'' on the negative z-axis at $z \approx -1500~{\rm km}$ . Thus at a time of $1.08~{\rm s}$ , $\rho_{\rm max} = 4.52$ $\times$ $10^8~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ and $T_{\rm max} = 8.53$ $\times$ $10^9~{\rm K}$ are reached. These maximum values are comparable with the spherically symmetric case discussed at the beginning of this section. The reason for this is the high degree of symmetry of the colliding shock fronts. Astoundingly those maximum values are reached comparatively far away from the center in a region of relatively low initial density.

The initial condition considered in the last paragraph seems more realistic than the more symmetric cases. Thus, this model was studied also in simulations with higher numerical resolution. Here, a resolution dependence of the peak temperatures and densities became apparent - an effect that is discussed in Sect. 3.2 below.

Finally, a series of models with a different initial mass distribution was investigated. Sequence "B'' corresponds to models with $0.9~\ensuremath{M_\odot}$ of C/O and $0.1~\ensuremath{M_\odot}$ of He and the parameters given in Table 3. They support the generality of the results discussed above. The dependence on the initial flame geometry is analogous to case "A'', but maximum densities and temperatures achieved are generally lower. This is consistent with the fact that the shock waves originating from the He shell have to cross a larger pressure gradient before they reach the center in model "B''. Still, most of the models reach the conditions for core ignition described in Sect. 2.4 (the "zr'' models did not ignite with mass distribution "A'' either).

3.2 A resolution study

The numerical resolution of the parameter study presented above was relatively low with a spatial resolution of $23.2~{\rm km}$ only. Hence, for one of the most interesting models, z4.24A_2d_256, we performed a resolution study. Globally, we found a good agreement of the temporal evolution of the state variables. But the maximum densities and temperatures at the central shocks collision showed a significant dependence on the numerical resolution of the simulation (see Fig. 6). With decreasing cell size the maximum values increase almost exponentially. This trend continues down to the smallest tested cell width (7.4 km) and no indication for numerical convergence is found. Given that the shock collision is an almost perfect spherically symmetric problem, we reduced to one dimension in order to facilitate even higher resolutions. The state variables on the r-axis of simulation s4.10A_2dq_256 at $t = 0.25~{\rm s}$ linearly interpolated to higher resolutions were taken as initial condition. The results of the corresponding one-dimensional runs with increasing resolution are presented in Fig. 7. It shows an analogous exponential increase of the temperature and density maxima down to a cell width of 0.72 km. Below that value, however, the temperatures exceeded the limits of our equation of state. Here, more micro-physical processes would have to be taken into account and could restrict further growth of density and temperature.

But even without this restriction, from analytical considerations it is clear that the density growth ratio must be limited, although it can get very big. In the idealized case of infinitely strong shock reflection at the origin in spherical geometry (and constant initial conditions of $\rho_0(r) = 1$ , u₀(r) = -1, and p₀(r) = 0), its value is $\left(\frac{\gamma + 1}{\gamma - 1}\right)^3$ after the reflected shock (cf. Glaister 1988). Another problem that can be solved semi-analytically is that of a spherical shock wave that starts from infinity, converges to the origin, and is reflected there (cf. Ponchaut et al. 2006; Landau & Lifschitz 1991; Guderley 1942). While the shock is moving inwards it grows in strength because of its surface decrease. Close enough to the center for the strong shock approximation to hold, it can be shown that $p \sim r^{-\frac{2(1-n)}{n}}$ , $\rho = {\rm const.}$ , and $\vert u\vert \sim r^{-\frac{1-n}{n}}$ with n being a constant between 0 and 1 that has to be numerically determined and depends only on the ratio of specific heats $\gamma$ . According to Landau & Lifschitz (1991), at the time of perfect focusing the density has increased by a factor of 20.1 for $\gamma = \frac{7}{5}$ and the reflection is accompanied by another increase by about a factor of 145. Thus a maximum compression of about 2900 can be achieved and in the case of $\gamma = \frac{4}{3}$ , relevant in our case, those numbers are even bigger.

In our simulations the shock that propagates inwards into the WD due to the He detonation is followed by an extended wave of matter that is also flowing inwards. Therefore the case might lie somewhere in between these two idealized situations: The shock is amplified by the geometrical focusing effect and thus compresses a significant part of the core that is then burned at higher densities. On the other hand the matter succeeding the in-flowing shock leads to even higher densities after the shock is reflected. Taking all those properties into account, it becomes clear why the ignition conditions for a core detonation can be fulfilled so well in our simulations. Provided that the shock is symmetric enough also an off-center ignition at a point with much smaller initial density like in model z4.24A_2d_256 is possible in this scenario.

$\begin{figure} \par\includegraphics[width=8.35cm,clip]{8438fig6} \end{figure}$	Figure 6: Resolution study of the maxima in the central shock collision of model z4.24A_2d.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.35cm,clip]{8438fig7} \end{figure}$	Figure 7: One-dimensional resolution study of the maxima in the central shock collision of model s4.10A.
Open with DEXTER

As the minimum shock surface that we can reach depends on the grid resolution, the maximum temperature and density values that are achieved through the geometric amplification in the shock implosion naturally also do so. In our simulations a successful second detonation ignition thus depends on the resolution and it seems likely that all the models in Table 3 would explode provided that they were well enough resolved.

Therefore, the results given in Table 3 have to be interpreted as a conservative estimate of the possibility of a detonation ignition. If at a certain grid resolution the ignition conditions are fulfilled, a second core detonation must happen. However, if the ignition conditions are not fulfilled a core detonation can still not be excluded and might be achieved with increasing resolution.

3.3 Three double-detonation supernova simulations

As a first step to compare the simulations of sub-Chandrasekhar mass models with actual SN Ia observations, measurable quantities like the nucleosynthetic abundances of an explosion have to be determined. To this end, two complete double detonations of model z4.24A until $t = 2.0~{\rm s}$ have been simulated. One of them was set up in 3D with 256³ cells for the full WD. For comparison with previous results, the second, 2D rotationally symmetric simulation was also done with this low resolution (128 $\times$ 256 cells).

Table 4: Results of three complete double detonation simulations. All nucleosynthetic abundances are given in solar masses. $M_{\rm loss}$ is the mass loss over the grid boundaries.

In order to test the influence of symmetry in the initial flame configuration in the He shell on the explosion abundances, besides these maximally asymmetric cases a 2D double detonation of the spherically symmetric model s4.10A was also simulated. The results of all three simulations are given in Table 4. For model z4.24A_2d_128 also the time evolution of the mass fractions of the species is plotted in Fig. 8 and will be discussed in the following.

At first, the $\ensuremath{^{4}{\rm He}}$ in the shell is burnt completely to $\ensuremath{^{56}{\rm Ni}}$ . But the mass of the latter is somewhat reduced due to some mass loss through the boundaries of the computational domain. This is clearly visible also in the decrease of the total mass. The reason for this is that the co-expanding grid was set up to track a mass shell at $0.85~\ensuremath{M_\odot}$ . Since the matter above this shell expands with a larger velocity, some of it leaves the numerical grid. The tracking of a bigger mass shell would reduce this mass loss, but it would also reduce the resolution of the colliding shocks in the core.

After about $1.2~{\rm s}$ the collision of the shock waves resulting from the He detonation gives rise to conditions for a detonation in the C + O core. In the corresponding grid zones, this detonation was artificially triggered by setting up a second level set. At the high densities near the center the C/O mixture is burnt completely to $\ensuremath{^{56}{\rm Ni}}$ in NSE with a small amount of $\ensuremath{^{4}{\rm He}}$ . But the $\ensuremath{^{4}{\rm He}}$ recombines to $\ensuremath{^{56}{\rm Ni}}$ after a short time due to the density decrease caused by expansion. Farther away from the center at smaller densities the burning product is $\ensuremath{^{28}{\rm Si}}$ and some low-density matter at the border of the core remains unburnt. In this way $0.40~\ensuremath{M_\odot}$ of $\ensuremath{^{56}{\rm Ni}}$ and $0.51~\ensuremath{M_\odot}$ of $\ensuremath{^{28}{\rm Si}}$ are produced after about $1.5~{\rm s}$ and $0.01~\ensuremath{M_\odot}$ of C/O remain unburnt.

The results of the 3D case z4.24A_3d_128 are very similar to the 2D one. However, $0.05~\ensuremath{M_\odot}$ more $\ensuremath{^{56}{\rm Ni}}$ were produced. The main reason for this difference is the loss of different amounts of mass over the boundaries of the computational domain, which amounted to $0.08~\ensuremath{M_\odot}$ in 2D and only $0.05~\ensuremath{M_\odot}$ in 3D. The explanation for the smaller mass loss is a simple geometric effect: In 3D Cartesian coordinates there is more volume outside the star that lies on the numerical grid than in 2D cylinder coordinates. For example, the largest distance from the center is the length of the space diagonal ( $\!\sqrt{3}~a$ , when a is the side length), whereas in 2D it is the length of the diagonal in the rz-plane ( $\!\sqrt{2}~a$ ).

Due to the spherical symmetry in the model s4.10A_2dq_128 He burning and shock convergence are faster. Therefore all the reactions finish already after about one second, but according to Table 4, almost the same nucleosynthetic abundances are produced. This indicates that the density distribution in the core in case of shocks converging asymmetric to the density distribution differs only slightly from the perfect spherical symmetric case. As the now following core detonation propagates supersonic, the final abundances are independent of the position of the ignition point.

In all three cases considered the $\ensuremath{^{56}{\rm Ni}}$ -mass and with it the explosion energy are rather large. In fact, they reach the range of (weaker) normal SNe Ia. However, these numbers have to be taken with care since the density limit of 3 $\times$ $10^7~\ensuremath{{\rm g}~{\rm cm}^{-3}}$ adopted here for separating burning of C/O to $\ensuremath{^{56}{\rm Ni}}$ and $\ensuremath{^{28}{\rm Si}}$ , respectively, is quite uncertain, but will have an impact on the final abundances.

$\begin{figure} \par\includegraphics[width=8.4cm,clip]{8438fig9} \end{figure}$	Figure 9: Mass fractions of the species (averaged over all angles) of models z4.24A_2d_128 (solid line) and s4.10A_2dq_128 (dashed line) in radial velocity space at $t = 10~{\rm s}$ .
Open with DEXTER

We could not calculate synthetic spectra or light curves of our simulations yet. First order spectra and light curves are given, e.g., in Livne & Arnett (1995) for models very similar to our z4.24A case. For comparison with their results the distributions of the species of models z4.24A_2d_128 and s4.10A_2dq_128 in radial velocity space (averaged over all angles) are presented in Fig. 9. Both are very similar except for a minor difference in the unburnt C/O, which is more abundant in the spherically symmetric case. The velocity ranges are comparable to model 6 of Livne & Arnett (1995), which had a similar setup ( $M_{\rm C/O} = 0.8~\ensuremath{M_\odot}$ , $M_{\rm He} = 0.17~\ensuremath{M_\odot}$ ), but resulted in a larger $\ensuremath{^{56}{\rm Ni}}$ mass ( $0.648~\ensuremath{M_\odot}$ ). Also, the maximum velocity of $\ensuremath{^{56}{\rm Ni}}$ and the minimum velocity of $\ensuremath{^{28}{\rm Si}}$ from core burning are shifted to significantly larger values. This could be caused by a difference in the transition density to NSE that might also be responsible for the difference in $\ensuremath{^{56}{\rm Ni}}$ mass.

In order to test possible asphericities of the ejecta in model z4.24A_2d_128 the abundances are plotted separately for the z-, r-, and -z-directions in Fig. 10. Despite the off-center explosion, the differences in the velocity distribution are minor. Only in -z-direction a shift in the maximum velocities of the $\ensuremath{^{56}{\rm Ni}}$ from the core detonation and of the $\ensuremath{^{28}{\rm Si}}$ is visible.

$\begin{figure} \par\includegraphics[width=8.35cm,clip]{8438figa} \end{figure}$	Figure 10: Mass fractions of the species of model z4.24A_2d_128 in velocity space along the z-, r-, and negative z-direction at $t = 10~{\rm s}$ .
Open with DEXTER

4 Summary and conclusions

We have studied the sub-Chandrasekhar model of SNe Ia by means of a series of two-dimensional and a few three-dimensional simulations with different initial conditions. The numerical scheme used is based on the PPM algorithm, and the propagation of the detonation front is modeled applying the level set technique. This novel implementation allowed for an accurate treatment of the hydrodynamic features of the sub-Chandrasekhar model (such as shock waves and thin detonation fronts) in multiple dimensions.

With our generic 1- $\ensuremath{M_\odot}$ model we performed a parameter study involving different He masses and ignition geometries of the initial He detonation. We find that a detonation in the He shell can clearly trigger a second detonation in the core and at least for the mass configurations studied, the double detonation seems to be a very robust process, which works without any "fine-tuning'' of our model parameters. In almost all of the simulations performed, the ignition conditions for a core detonation were reached at or near the center of the WD as a result of the convergence of the shock from the He shell detonation and in the few other cases the failure is likely a result of the finite numerical resolution. The high maximum densities and temperatures that are needed for detonation ignition and the fact that the shock is not diluted too much when propagating inwards against the pressure gradient are made possible through a geometrical shock amplification effect that appears in spherical symmetry: As the surface of the shock decreases its strength has to increase. The high compression that can be achieved theoretically in this way, especially if full convergence is reached and the shock is reflected (cf. Sect. 3.2), could even allow successful double detonations with considerably smaller He masses. The shock amplification that is reached, however, turns out to be resolution dependent on the numerical grid as it is coupled to the smallest resolvable shock surface. Taking this into account it is very likely that all our models would explode, if they were simulated with sufficiently high resolution.

The question of whether an incineration could also happen directly at the core-shell interface, was not addressed here. This could of course prevent the spherical shock convergence mechanism from playing a role, at least in parts of the parameter space. Therefore, edge-lit detonations are postponed to a separate study.

The complete double detonation simulations that were performed for the models z4.24A and s4.10A resulted in $\ensuremath{^{56}{\rm Ni}}$ -masses of about 0.40 to $0.45~\ensuremath{M_\odot}$ . Are the studied events thus good candidates for normal SNe Ia? Most likely not. Starting with a He-shell detonation, all of our models show a significant amount of rapidly expanding $\ensuremath{^{56}{\rm Ni}}$ in the outer layers. In the observed spectra of normal and sub-luminous events, however, this has never been observed (cf. Arnett 1997; Branch et al. 1982; Woosley et al. 1986; Branch 1984; and also the discussion in Livne & Arnett 1995). The only exception is the super-luminous SN 1991T. Thus our solar-mass models most likely will not be able to explain normal or sub-luminous SNe Ia and, given the robustness of the explosion mechanism of the model, the only conceivable explanation for the lack of observations of corresponding SN Ia events is that the progenitors are not realized in nature - or are very rare.

For a better agreement with observations a reduction of the He-shell mass would be the most obvious choice. However, then a significantly larger core mass would most likely be required, because otherwise the He would not detonate. In this case the core density would be much higher resulting in a very large Ni mass. This would again not be a candidate for normal SNe Ia, but it might be a promising candidate for objects like SN 1991T.

Further work could cover model improvements like a more realistic treatment of nuclear reactions including the calculation of real reaction rates depending on the actual thermodynamic state of the burnt matter. This would make an investigation of the onset and the explosion dynamics of the core detonation possible. Also, more extended parameter studies, especially towards the lower end of possible He masses, would be interesting. For a comparison with observations the calculation of synthetic light curves and spectra would also be desirable.

Double-detonation supernovae of sub-Chandrasekhar mass white dwarfs