text:
# Long-term variability in the X-ray emission of RX J0720.4-3125

## 1 Introduction

RX J0720.4-3125 (Haberl et al., 1997) belongs to the group of radio-quiet neutron stars discovered by ROSAT. As these objects are characterised by soft X-ray emission from the surface, it was hoped that high resolution spectroscopy of those sources with the Chandra’s Low Energy Grating Spectrometer (LETGS) or XMM-Newton’s Reflection Grating spectrometer (RGS) would reveal line features from a thin atmosphere that could provide information on the surface gravity of the neutron star and hence on the equation of state of the neutron star matter. However, deep observations with XMM-Newton and Chandra of RX J0720.4-3125 (Paerels et al., 2001; Kaplan et al., 2003) and RX J1856.5-3754 (Burwitz et al., 2001) have in that respect been somewhat of a disappointment, as the spectra are almost perfectly well described by a blackbody spectrum. The only unusual feature of the X-ray spectra is that the blackbody emission implies a radius too small for a neutron star and, moreover, underpredicts the observed optical and UV flux (Kaplan et al., 2003; Burwitz et al., 2001; Motch et al., 2003). This implies that another cooler component is present and that the X-ray emission only comes from part of the surface.

RX J0720.4-3125 has a low absorbing column density ($`N_\mathrm{H}10^{20}\mathrm{cm}^2`$), is at an approximate distance of $`300\mathrm{p}\mathrm{c}`$ (Kaplan et al., 2003), and shows sinusodial pulsations with a period of 8.39 s (Haberl et al., 1997; Cropper et al., 2001). The upper limit to the period derivative $`\dot{P}<n\times 10^{13}`$ ss-1, with $`n`$ ‘a few’ (Kaplan et al., 2002), implies a magnetic field strength $`B<3\sqrt{n}10^{13}`$ G and a characteristic age $`\tau _c0.5P/\dot{P}>(1.3/n)`$ Myr. As RX J0720.4-3125 was considered to be a perfect blackbody spectrum with a temperature of $`k\mathrm{T}_{\mathrm{BB}}86\mathrm{eV}`$ it was observed several times by XMM-Newton for calibration purposes. However, as we will demonstrate here, the spectrum of RX J0720.4-3125 has slowly hardened and cannot be described any more by a blackbody spectrum, while the pulse shape has become narrower and the phase dependence of the spectrum changed with time.

While we were finalising this paper, Haberl et al. (2003b) released a pre-print in which they show that the spectrum of RX J0720.4-3125 depends on the pulse phase. They also discuss changes in the spectrum as found with EPIC, but assign these to calibration inaccuracies.

## 2 Observations

Table 1 gives the log of XMM-Newton observations of RX J0720.4-3125. Since the source spectrum is rather soft, and most of the emission comes from below $`1.8`$ keV, we use RGS (den Herder et al., 2001) data for the spectral analysis; due to the nature of the gratings, the effective area of the RGS is rather insensitive to changes in CCD gain and charge transfer inefficiency. The RGS was operated in normal spectroscopic mode, yielding a time resolution of 4 s, making it unsuitable for a proper timing analysis. Therefore, for this purpose we use EPIC/PN (Strüder et al., 2001) since it has the largest effective area among the instruments on board XMM-Newton, and operates in a high-time resolution mode.

During the first four observations, PN was operated in full-frame mode, which produces an image of 378 $`\times `$ 384 pixels with a time resolution of 73.4 ms. During the last two observations the PN was switched to small-window mode, yielding images of 63 $`\times `$ 64 pixels with a time resolution of 6 ms. The thin filter was used in revolutions 0078, 0533, and 0534, the medium one in revolution 0175, and the thick one in revolution 0622. The observation of revolution 0711 started with the thin filter, and later switched to the medium one. All data were processed using the standard XMM-Newton software SAS version 5.4.1.

### 2.1 Spectral Analysis

Figure 1 shows the RGS spectra for the 6 revolutions. The indivual RGS’s agree with each other within the statistical uncertainties, so the data from both RGS’s have been combined. From the Figure it is apparent that the spectrum of RX J0720.4-3125 gradually changes, while the total flux appears to remain fairly constant. The change of the spectrum is most noticeable in orbit 0711, where there is an increase of the flux in the 15 Å to 25 Å range, and a drop in the flux above 30 Å.

To verify that the RGS instrument is stable and that its effective area does not change, we have also analyzed data from two sources that are stable on a time scale of years.

The first source is the supernova remnant 1E 0102-72.3 in the Small Magellanic Cloud, a calibration source that is regularly observed with the RGS. From the strength of the bright emission lines from this source in the range 13 to 23 Å, as measured in data from orbits 0065 to 0711, we find that the effective area of the instrument in this range changes by less than 3%. In the long wavelength range the SNR 0102-72.3 spectrum has a weak carbon line at 33.7 Å. The upper limit on possible changes in the effective area of the RGS instrument at this wavelength is about 10%, mainly determined by the limited statistics in this line. All data from 1E 0102-72.3 are consistent with no change in the RGS effective area. The second source is the calibration source Mrk 421. An upper limit to any change in the neutral Oxygen absorption edge around 22.8 Å, which is a tracer of possible instrument contamination (e.g. ice), is 5%.

Since all upper limits on possible changes in effective area are well below the changes seen in RX J0720.4-3125 we conclude that the X-ray spectrum of RX J0720.4-3125 itself is subject to change over the course of about 3 years.

For the spectral analysis we use the XSPEC spectral fitting package (Arnaud, 1996). As noted by Paerels et al. (2001) the data of orbit 0078 can be well fitted by an absorbed blackbody spectrum. With the improved calibration of the RGS we find a temperature of $`kT_{\mathrm{BB}}=86.7\pm 0.3`$ eV and an interstellar column density of $`N_\mathrm{H}=(1.41\pm 0.07)\times 10^{20}`$ cm-2, in reasonable agreement with the values obtained with the Chandra LETGS instrument, $`kT_{\mathrm{BB}}=81.4\pm 1.3`$ eV and $`N_\mathrm{H}=(1.32\pm 0.14)\times 10^{20}`$ cm-2, by Kaplan et al. (2003). The LETGS observation was made in February 2000, four months before the first XMM-Newton observation.

Fits with the same model to the remaining spectra yield a temperature increase with time, as well as changes as large as $`\mathrm{\Delta }N_\mathrm{H}5\times 10^{20}`$ cm-2 in the interstellar absorption. The model, however, does not provide satisfactory fits to all the spectra.

Next, we fit all the data simultaneously. Since large variations of the interstellar absorption over such a short period are unlikely, we assume that $`N_\mathrm{H}`$ remains constant; we therefore couple this parameter for all observations, while we let the blackbody temperature and normalization vary freely between observations (Table 2). Although these fits serve to emphasize the gradual hardening of the spectra, half of the spectra, especially that of revolution 0711, are poorly fit by this model.

Following Haberl et al. (2003a) and van Kerkwijk (2004), next we fit the data with a model that consists of a blackbody affected by additional absorption by a broad gaussian line, as might be expected in the case of cyclotron absorption. In these fits, we constrain the parameters of the blackbody model to be the same for all the observations, whereas the parameters of the line are allowed to vary independently. Although this model fits the data well, in most cases we find that the width of the line is larger than its central energy, or the central energy of the line lies outside the spectral range covered by the data. In those cases, the tail of the broad gaussian just serves to attenuate the emission of the source at long wavelengths. Furthermore, we find no clear trend in the central energy or the width of the line.

We find that the ratio of the spectra of the last and the first observation is close to a power law. We therefore fit the data to an empirical model that consists of a blackbody multiplied by a power law $`E^\mathrm{\Gamma }`$, all affected by interstellar absorption. While the index of the multiplicative power law is allowed to change between observations, for these fits we constrain the parameters of the blackbody and the interstellar absorption to be the same in all observations. While it is difficult to assign a physical interpretation to this model, it provides an acceptable description of the data in the RGS range ( 10 - 38 Å, see Table 2), it has fewer parameters than the gaussian absorption model and, in addition, the index of the power law increases steadily over the course of the observations.

### 2.2 Timing analysis

Starting from the raw data, we first produce a list of calibrated events. To reduce pile-up, in the next step we select only single events as well as events that are not affected by some of the imperfections (bad columns, hot pixels, etc.) of the CCDs. We extract events within a $`39`$ arcsec circle centred on the source. We barycenter these events using the SAS routine BARYCEN version 1.13.4, and we then separate the events according to their energy in 3 event lists; the bands that we use are 0.1 to 1.2 keV, 0.1 to 0.4 keV, and 0.4 to 0.8 keV, respectively.

For each observation we find the best period in the full band using an epoch folding technique; in all cases we find a period of 8.391 s, consistent with the value previously found for this source by Kaplan et al. (2003). We then produce folded light curves in the three bands, and we also compute a folded hardness-ratio light curve from the ratio of the 0.4–0.8 keV and the 0.1–0.4 keV light curves. In Figure 2 we show the 0.1-1.2 keV and the hardness-ratio light curves. For each observation we define the phase such that the maximum of the full-band light curve occurs at phase zero; the phase of the hardness ratio light curves is the same as for the full-band light curves.

The pulse profile in the 0.1–1.2 keV band, as well as the hardness-ratio pulse profile, change from one observation to the other. The first panel in Figure 2 shows a sinusoidal fit to the pulse profile during the first observation; the same sine function is overplotted to the full-band pulse profiles obtained from the other observations. It is apparent that the pulse profile becomes narrower with time.

At the same time, the hardness-ratio pulse profile also changes. In the first observation there is a clear modulation, and the hardness-ratio profile leads the full-band light curve by $`0.064\pm 0.017`$ in phase. In the following observations the amplitude of the hardness-ratio modulation decreases and the phase difference between the full-band and the hardness-ratio light curves is consistent with zero. Eventually, in the last observation the modulation increases again, but now the hardness-ratio light curve lags the full band-light curve by $`0.126\pm 0.010`$ in phase.

## 3 Conclusions and discussion

The XMM-Newton data of RX J0720.4-3125 show that the spectrum of the source changes on a time scale of years, the first time ever that the X-ray spectrum of an isolated neutron star, other then soft gamma-ray repeaters or anomalous X-ray pulsars, is seen to change. Whereas the changes are most pronounced in the last observation, we think that the actual change is gradual, as witnessed by a gradual increase in the temperatures derived from the blackbody fits; or by a gradual increase in the index of the powerlaw in the fits with a blackbody multiplied with a power law (Table 2 and Figure 1). The spectral changes are accompanied by an energy-dependent change in the pulse shape; in particular the pulse phase where the spectrum is hardest has moved with respect to the phase of maximum flux (Figure 2).

The phase (i.e. angle) dependent spectrum of single neutron stars is currently not explained. The broad absorption features have been interpreted as a proton-cyclotron absorption feature (Haberl et al., 2003a). In pulsars with a strong field (probably stronger than the limit for RX J0720.4-3125) the absorption feature may be weakened by the strong-field quantum electrodynamics effect of vacuum resonance mode conversion (Lai & Ho, 2003). The neutron star spectra have also been interpreted as due to cyclotron-resonance scattering of the spectrum from the surface of the neutron star, by electron-positron pairs in the magnetosphere (Ruderman, 2003). For both interpretations the spectrum is likely to be angle and energy dependent, in accordance with the variation of the X-ray spectrum (as measured by hardness ratio) with pulse phase.

To explain the gradual, long-term variation we consider two general possibilities: either the intrinsic spectrum of the neutron star changes, or our view of the neutron star changes. The intrinsic spectrum of the neutron star could change, due to energy release deep in the neutron star, due to a glitch for example, causing the surface to gradually become hotter. We consider this unlikely as the explanation for the changing spectrum of RX J0720.4-3125, because it does not explain the change in pulse shape. Also this model would predict the total flux to increase with the temperature, in contrast to what is observed. Another possibility, valid for Ruderman’s model, would be that the electron-positron plasma surrounding the neutron star changes. So far, there is no specific prediction in the Ruderman model for changes in the magnetospheric plasma on a year-long time scale.

Therefore, we suggest that the variation in the spectrum of RX J0720.4-3125 is caused by a change in the angle under which we see the emitting region and/or the covering electron-positron plasma, caused by precession of the neutron star. Precession arises when the form of the neutron star deviates from a perfect sphere and its rotation is not around a principal axis (as reviewed by e.g. Link 2003). The changes in pulse form and phase of the radio pulsar B1828$``$11 are successfully described with free precession, and have a time scale of years (Stairs et al., 2000), i.e. comparable to the time scale on which the variations in the spectrum and pulses of RX J0720.4-3125 occur. Whether a change in viewing angle can produce the observed hardening of the spectrum, with only a modest change in overall flux, will depend on the details of the emission model, and requires further investigation.

To further investigate the cause of the long-term variability it will be useful to investigate the phase-dependence of the spectra in more detail; to obtain a reliable period derivative of the pulses, and look for the phase changes expected from precession; and finally to see whether the changes continue.

###### Acknowledgements.

The authors whish to thank Frank Haberl for a very useful discussion on the effect of EPIC-pn observation modes and Jan-Willem den Herder for his continuous stimulation.


format: text

text:
# Direct Numerical Simulations of Type Ia Supernovae Flames II: The Rayleigh-Taylor Instability

## 1 INTRODUCTION

Accelerating a thermonuclear flame to a large fraction of the speed of sound (possibly supersonic) is one of the main difficulties in modeling Type Ia supernovae (see the recent review by Hillebrandt & Niemeyer 2000 and references therein). Numerical results have shown that a prompt detonation does not yield the right abundance pattern to account for the observations (Arnett et al., 1971); consequently, the burning must begin as a subsonic flame. One dimensional simulations show that the flame needs to reach approximately one third of the speed of sound to properly account for the explosion energetics and nucleosynthetic yields (Woosley et al., 1984; Nomoto et al., 1984), although, a transition from a deflagration to a detonation at some late stage of the explosion (Khokhlov, 1991; Niemeyer & Woosley, 1997) is not ruled out. Laminar flame speeds are too slow by several orders of magnitude (Timmes & Woosley, 1992), so instabilities and the interaction with flame-generated turbulence are turned to in hopes of providing a mechanism to considerably accelerate the flame.

In this paper, we focus on the effects of the Rayleigh-Taylor (RT) instability (Taylor, 1950; Chandrasekhar, 1981) on the flame. After ignition, the subsonic flame moves outward from the center of the star. The pressure is essentially continuous across the flame. Furthermore, since the star has time to expand in response to the burning, this pressure also remains constant in time. As the flame propagates, it leaves behind hot ash that is less dense than the cool fuel. Gravity points toward the center of the star, so this flame front is RT unstable. Bubbles of hot ash will rise and try to exchange places with spikes of cool fuel, increasing the surface area of the flame and leading to an enhancement of the flame speed. It is possible that this acceleration, operating on the stellar scale, can bring the flame speed up to a significant fraction of the speed of sound. Furthermore, the shear layer between the fuel and ash which develops as the instability evolves is Kelvin-Helmholtz unstable, creating turbulence in the region of the flame. As the flame surface is wrinkled by the RT instability, it will interact with this turbulence. We look in detail at the interaction between the RT instability and the flame on small scales through 2-dimensional spatially resolved simulations in conditions appropriate to the late stages of a Type Ia supernova explosion.

The difference in scale between the white dwarf and the flame width is enormous (up to 12 orders of magnitude), making direct numerical simulations of the whole explosion process impossible. Nevertheless, simulations of RT unstable flames have been attempted for some time, both on the full stellar scale and as smaller, constant density micro-physical studies. Various techniques and approximations are used to follow an unresolved flame on these scales. Some of the earliest such calculations (Müller & Arnett, 1982, 1986) used donor-cell advection to follow a combustion front through the star, burning most of it, and releasing enough energy to produce an explosion. In their model, the flame was propagated essentially through numerical diffusion. The RT instability was resolved only on the very largest scales of the star, missing the turbulence generated on the small scales. Livne (1993) used a two-dimensional implicit Lagrangian-remap hydrodynamics code with a mixing-length subgrid model to follow a RT unstable flame. The diffusion of reactants in the Lagrangian-remap scheme is much lower than that used in Müller & Arnett (1982, 1986), so a more physical flame speed could be used. The flame speed was chosen to be either the conductive (laminar) or turbulent speed—which ever was larger. Again, the whole star was modeled, with coarse resolution, showing that the RT instability sets up and greatly convolutes the flame front, but the resulting acceleration left the flame too slow to match observations. Livne (1993) suggested that more resolution is needed to get the small scale physics right.

The effects of turbulence on the burning were investigated by Niemeyer & Hillebrandt (1995b); Khokhlov (1995). Two-dimensional, small scale RT simulations using a thickened-flame model, coupled to an explicit Lagrange-remap PPM (Colella & Woodward, 1984) implementation to model a thin flame front were performed by Khokhlov (1993) and later extended to 3-d (Khokhlov, 1994). In this model, a reaction-diffusion equation for the fuel concentration is solved along with the Euler equations, with the conductivity and reaction rate modified to yield the desired laminar flame speed and a flame that was a few computational zones wide. The thermal structure of the true flame is stretched greatly, and, as a result, curvature effects of the flame are not properly included. Energy is deposited behind the flame to couple this reaction-diffusion model to the hydrodynamics. Similar techniques have also been used for terrestrial flames (see for example O’Rourke & Bracco 1979). Simulations of the RT instability at a density of $`10^8\mathrm{g}\mathrm{cm}^3`$ lead the author to conclude that while the flame accelerates, it is too small of an acceleration (when scaled up to the size of the star) to account for the observations, and that a deflagration-detonation transition may be more promising. Owing to the difference in behavior between two-dimensional and three-dimensional turbulence, Khokhlov (1994) concludes that the turbulence was more effective in increasing the flame speed in 3-d. Khokhlov (1995) demonstrated scaling relations for the 3-d RT unstable flame and concluded that the turbulent flame speed does not depend on the small scale burning or the laminar flame speed.

The independence of the turbulent flame speed on the laminar speed was also postulated by Niemeyer & Hillebrandt (1995b), who performed a full-stellar model based on Kolmogorov scaling of turbulence from the large stellar scales down to the flame scales. Here, a turbulent scaling model using the Gibson length (the scale at which an eddy is completely burned during a single turnover) as the small scale cutoff yields a turbulent flame velocity prediction for the grid scale. This was coupled with a subgrid model for the turbulent kinetic energy, and implemented within the PPM algorithm. Their simulations of the flame propagating through the whole star showed the RT instability dominating the flow, but the flame did not accelerate enough to produce an explosion.

More recently, large scale simulations of pure deflagration explosion have been performed using two different methods, level-sets (Reinecke et al., 1999a, b, 2002) and a thickened-flame model (Gamezo et al., 2003) with a Sharp-Wheeler (Sharp, 1984; Glimm & Li, 1988) RT subgrid model. Both groups using PPM as their underlying explicit hydrodynamics method. Different numbers and locations of the ignition points are used, but the 3-d simulations by both groups release enough energy to explode the star, with the burning proceeding solely as a deflagration.

These studies have demonstrated the tremendous effect the RT instability has as the dominant acceleration mechanism of the flame—to the point where, in some simulations, a deflagration alone can unbind the star. However, in all the cases discussed above, some model for the flame or turbulence on scales smaller than the grid resolution was required. In this paper, we present fully resolved simulations of RT unstable flames at the densities appropriate to the late stages of the explosion—no flame model is used. Fully resolving the flame means that the effects of curvature and strain on the flame are implicitly accounted for. Curvature effects on flames in conditions close to those studied here were found to be important in Dursi et al. (2003), where laminar flames were propagated in diverging and converging spherical geometries. Here, the effects can be even larger as the flame will experience sharp kinks as it is distorted by the RT instability. To achieve converged steady-state laminar flames, approximately five grid point resolution in the thermal width is required. This limits the size of the domain that can be modeled significantly (although the adaptive mesh refinement employed here helps greatly), but frees us from the need to deal with subgrid models. Resolved simulations at the smallest scales complement the full star flame-model calculations, and it is hoped that progress from the small scales up will allow us in the future to develop and calibrate more accurate subgrid models.

Fully resolved multidimensional Type Ia-like flame simulations are rare (Niemeyer & Hillebrandt, 1995a; Zingale et al., 2001; Bell et al., 2003a), mainly because the flame moves so slowly, requiring excessive numbers of timesteps for compressible codes. Resolved calculations of the reactive RT instability have been performed before (Vladimirova & Rosner, 2003), but using model flames (a KPP-type reaction, Kolmogoroff et al. 1937) and a Boussinesq approximation, under conditions that are not directly applicable to Type Ia SNe. When extrapolated to the astrophysical regime, they predict a flame speed independent of the laminar speed, similar to that found by Khokhlov (1995). However, compressibility was neglected in these calculations. Here we use a low Mach number formulation of the equations of motion, which retains the compressibility, but frees us from the restrictive acoustic timestep constraint, and we use realistic input physics (reaction rate, EOS, etc.) appropriate to the conditions in a Type Ia SNe explosion.

Flame instabilities stretch and wrinkle the flame, increasing the bulk burning rate. The SNe flame is subject to the Landau-Darrieus (Landau, 1944; Darrieus, 1938) and RT instabilities, as well as the interaction with turbulence. In a previous paper, (Bell et al. 2003a, henceforth referred to as paper I), we looked at Landau-Darrieus unstable flames, and validated the low Mach number method for astrophysical flames, confirming the results of earlier studies that the acceleration of the laminar flame is quite small, while finding no evidence for the breakdown of the cusps in the non-linear regime. In the present paper, we look in detail at the interaction between the burning and the RT instability. A RT driven flame behaves very differently from its laminar counterpart. Aside from looking at the increase in velocity over the laminar speed the wrinkling provides, we seek to determine whether there is a simple scaling relation that matches these results, and whether such a relation holds over the whole range of densities, where the character of the burning is expected to change dramatically.

In §2, we discuss the effects of the burning on the RT instability when we enter the distributed burning regime. In §3, the low Mach number numerical method and input physics are described, followed by the results of our simulations in §4. We conclude in §5 by discussing the implications of the RT burning process on the explosion mechanism for Type Ia SNe.

## 2  THE EFFECTS OF BURNING ON THE RAYLEIGH-TAYLOR INSTABILITY

The RT instability has been well studied analytically (Taylor, 1950; Layzer, 1955; Chandrasekhar, 1981) and numerically, both as an incompressible (see recent results by Young et al. 2001; Cook & Dimotakis 2001) and compressible (see for example Glimm & Li 1988; Calder et al. 2002) fluid. In the absence of burning, the growth of the multimode RT instability in the non-linear regime is characterized by two parameters, one describing the terminal velocity, and the other describing bubble merger (Glimm & Li, 1988). Bubbles are assumed to merge when the difference in their heights exceeds some critical multiple of the smaller bubbles size. This leads to an expression for the average position of the interface,

$$h=\alpha Agt^2,$$
(1)
(Sharp, 1984) where $`\alpha `$ is a function of the two parameters, and $`A=(\rho _{\mathrm{fuel}}\rho _{\mathrm{ash}})/(\rho _{\mathrm{fuel}}+\rho _{\mathrm{ash}})`$ is the Atwood number. Much theoretical, computational, and experimental work aims to measure $`\alpha `$, and some evidence suggests that it may be constant for a wide range of initial conditions (George et al., 2002). The velocity of the merged bubbles can be written as

$$v=c\sqrt{Agh},$$
(2)
where $`c1/2`$ (Davies & Taylor, 1950). These expressions were all formulated for the non-reactive RT instability and studies of rising bubbles. Numerical experiments of reactive RT instabilities in SNe conditions (Khokhlov, 1995) using a thickened-flame model found agreement with the above velocity expression, with $`c=0.5`$. Expressions of this form have been used as a subgrid model for the flame speed in full-star, pure deflagration calculations of Type Ia explosions in 1-d (Niemeyer & Woosley, 1997) and 3-d (Gamezo et al., 2003), resulting in a successful explosion of the star. These results suggest that the RT instability is the dominant acceleration mechanism for the flame.

Woosley (1990) and Timmes & Woosley (1992) argued that, unlike a non-reactive RT instability, where all wavelengths can grow (in the absence of surface tension), burning introduces a small scale cutoff to the instability. The growth rate for the pure RT instability is determined by the dispersion relation

$$\omega ^2=gkA$$
(3)
(Chandrasekhar, 1981), where $`\omega `$ is the angular frequency and $`k`$ is the horizontal wave number for the dominant mode. We note that we have ignored contributions due to surface tension. By equating the period for the growth of the instability to the laminar flame propagation timescale, Timmes & Woosley (1992) showed that there is a critical wavelength below which any perturbations will be washed out by burning. This wavelength is called the fire-polishing wavelength, $`\lambda _{\mathrm{fp}}`$:

$$\lambda _{\mathrm{fp}}=4\pi \frac{v_{\mathrm{laminar}}^2}{g_{\mathrm{eff}}},$$
(4)
where $`g_{\mathrm{eff}}=gA`$ is the effective gravitational acceleration. We note that this is a simple estimate, nevertheless, it provides a useful measure of the effectiveness of the RT instability on disrupting the flame. Table 1 lists the flame properties, as computed by the code described in §3, including $`\lambda _{\mathrm{fp}}`$ for the densities we are considering.

Stratification effects are negligible on the scales that we are considering, since we are modeling a tiny fraction of the star. The relative change in pressure over our domain is

$$\frac{\mathrm{\Delta }p}{p}=\frac{\rho g\mathrm{\Delta }r}{p}\frac{v_{\mathrm{SW}}^2}{c_s^2}=M_{\mathrm{SW}}^21,$$
(5)
where we have expressed this ratio in terms of the Mach number of the Sharp-Wheeler velocity for a RT flame encompassing our entire domain. Since we are resolving the structure of the flame, the domain size restrictions limit the speed-ups we see to about 10 times the laminar speed, which gives a maximum $`\mathrm{\Delta }p/p`$ of $`O(10^8)`$.

As the density of the white dwarf decreases, the flame speed decreases sharply while the flame width, $`l_f`$, increases. As Table 1 shows, at the densities we are considering, we pass from the regime where $`\lambda _{\mathrm{fp}}l_f`$ (low densities) to $`\lambda _{\mathrm{fp}}l_f`$ (high densities). When $`\lambda _{\mathrm{fp}}l_f`$, the RT instability will dominate the flame, to the point where it may no longer be continuous, and one can no longer draw a simple curve from one side of the domain to the other marking the flame position. This increase in complexity of the flame surface signifies that we are leaving the flamelet regime.

Niemeyer & Woosley (1997) made estimates of the density where the burning transitions from the flamelet regime to the distributed burning regime, in which a distinguished local flame surface is not apparent. Using the flame properties tabulated in Timmes & Woosley (1992) and equating the Gibson scale to the flame width, they found that for a 0.5 12C/0.5 16O flame, the burning will enter the distributed burning regime at a density of $`1`$$`5\times 10^7\mathrm{g}\mathrm{cm}^3`$. Niemeyer & Kerstein (1997) refined this estimate by suggesting that since the turbulence present in a Type Ia explosion is generated by the RT instability, it should follow Bolgiano-Obukhov rather than Kolmogorov statistics. They found a slightly lower transition density of $`O(10^7)\mathrm{g}\mathrm{cm}^3`$. At these densities, the turbulence is intense enough to disrupt the flame. We note that this transition density is right in the range where the fire polishing length becomes equal to the flame width, as discussed above.

This density range is also interesting in the context of deflagration to detonation transitions (DDT). Arguments based on nucleosynthesis and one-dimensional modeling suggest that if the burning in a Type Ia supernova were to transition from a deflagration to a detonation, this transition would need to occur at a density around $`10^7\mathrm{g}\mathrm{cm}^3`$ (see for instance Hoeflich & Khokhlov 1996). The correlation between this DDT density and the density marking the transition to the distributed burning regime was pointed out by Niemeyer & Woosley (1997) and Khokhlov et al. (1997). One-dimensional turbulence studies of the flame-turbulence interaction at these densities were performed by Lisewski et al. (2000), where it was observed that the turbulence disrupts the flame sufficiently in some cases to create local explosions as the hot ash comes in contact with pockets of cool fuel. If these regions of fuel are large enough, these local explosions may be able to initiate a DDT.

Direct numerical simulations of thermonuclear flames at these densities do not exist. Niemeyer et al. (1999) performed simplified flame model calculations in 3-d using a source term relevant to flames in Type Ia supernovae. They concluded that for densities larger than $`O(10^7)\mathrm{g}\mathrm{cm}^3`$, scaling the flame speed in direct proportion to the increase in surface area is a valid subgrid model. For densities below $`10^7\mathrm{g}\mathrm{cm}^3`$ however, they state that the effects of turbulence become important.

As Table 1 shows, the laminar flame speed in our density range is extremely subsonic—our highest laminar flame Mach number is $`2\times 10^5`$. Simulating a flame at these densities with a fully compressible code would be prohibitively expensive, requiring millions of timesteps. This in turn leads to a large accumulation of error. Instead, a low Mach number hydrodynamics code (Bell et al., 2003b) is employed for the present study. We discuss the numerical method in §3.

Together, the flame width and the fire-polishing length set strict limits on the densities and size of domains that can be addressed through fully resolved simulations. Steady state laminar flame simulations show that we need about 5–10 zones inside the flames thermal width,

$$l_f\frac{T_{\mathrm{ash}}T_{\mathrm{fuel}}}{\mathrm{max}\{T\}}.$$
(6)
We note that this width measurement is a factor of 2–3 narrower than the alternate measure defined as the width of the region where the temperature is 10% above the fuel temperature to 90% of the ash temperature that was used in Timmes & Woosley (1992). In order for the RT instability to develop, we need the computational domain to be at least one fire-polishing length wide, but in practice, we want 10 or more fire-polishing lengths in the domain for bubble merger to take place and the RT instability to become well developed. These restrictions limit us to consider densities $`1.5\times 10^7\mathrm{g}\mathrm{cm}^3`$ in the present study. The lower limit of densities we can consider is set only by what is attained in the supernova explosion itself—we stop at $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ here.

When the RT instability dominates the burning, it is the fluid motion mixing the fuel with the hot ash that controls the burning of the new fuel, rather than thermal conduction. Thus, the flame has a very large thermal width, governed by the mixing, at the lowest densities. To initiate a detonation, the temperature in a sufficiently large region needs to rise in unison, such that the over-pressure generated by the burning is strong enough to drive a shock just ahead of the burning layer. The critical size of this region is called the matchhead size. If this RT thermal width becomes as large as the critical detonation matchhead size, as tabulated in Niemeyer & Woosley (1997), it is possible that a DDT can occur. Furthermore, the size of the matchhead increases dramatically as the density decreases, and at low densities, it may be larger than the entire star. We note that recent simulations (Gamezo et al., 2003; Reinecke et al., 2002) indicate that a pure deflagration can be successful in unbinding the star. We will explore the possible of a DDT, using the scaling behavior we find for flames in the distributed burning regime.

## 3 NUMERICAL METHODS

The simulations were performed using an adaptive, low Mach number hydrodynamics code, as described in Bell et al. (2003b). The state variables in the inviscid Navier-Stokes equations are expanded in powers of Mach number, following Majda & Sethian (1985). The result is that the pressure is decomposed into a dynamic and thermodynamic component, the ratio of which is of order Mach number squared. Only the dynamic component appears in the momentum equation. These equations are solved using an approximate projection formalism (Almgren et al., 1998; Day & Bell, 2000), breaking the time evolution into an advection and a projection step. The advection is computed using an unsplit, second-order Godunov method that updates the species and enthalpy to the new time level and finds provisional velocities. A divergence constraint on the velocities is provided by the thermodynamics, and enforced in the projection step, where the provisional velocities are projected onto the space of vectors satisfying this constraint. Burning is handled through operator splitting. The low Mach number formulation allows us to follow the evolution of these slow moving flames—in this density range, the laminar speeds are $`\mathrm{M}10^5`$—without the constraint of the sound speed on the timestep. This affords us timesteps that are $`O(1/\mathrm{M})`$ larger than a fully compressible code would take. We are still free however to have large density jumps as we cross the flame. This method was validated for astrophysical applications by comparing solutions of 1-d laminar flames to fully compressible results (Bell et al., 2003b), and in our study of the Landau-Darrieus instability (paper I), where the growth rate computed from our calculations matched the theoretical predictions across a range of wave numbers. The code is unchanged from the description in Bell et al. (2003b). We review the input physics used in the present simulations below.

This method contrasts sharply with that used in the only previous study of the RT instability in Type Ia SNe, as presented in Khokhlov (1995). There, the fully compressible PPM algorithm was used, along with a model for a thickened flame to represent it on the small scales. The laminar flame speed matched the correct physical value, but the flame thickness was much larger than the physical value. In the present case, resolving the thermal structure frees us from the need for a flame model.

The equation of state consists of a Helmholtz free-energy tabular component for the degenerate/relativistic electrons/positrons, an ideal gas component for the ions, and a blackbody component for radiation, as described in Timmes & Swesty (2000). The conductivities contain contributions from electron-electron and electron-ion processes, and are described in Timmes (2000). All of the flames are half carbon/half oxygen, but only the carbon is burned. A single reaction, 12C(12C,$`\gamma `$)24Mg, is followed, using the unscreened rate from Caughlan & Fowler (1988). Since we do not expect the burning to proceed up to the iron-group elements at these low densities, this reaction alone is sufficient to model the nuclear reactions. The reaction rate for oxygen burning is several orders of magnitude slower than the carbon rate at the ash temperatures we reach, so neglecting oxygen burning is valid on our scales.

The calculations are initiated by mapping a one-dimensional steady state laminar flame (computed with the same code) onto our grid, shifting the zero point according to a random phase, 10 frequency sinusoidal perturbation. Beginning with a steady-state flame in pressure equilibrium virtually eliminates any transients from disturbing the flow ahead of the flame. In all of the results presented, the flame is moving downward in our domain. The transverse boundaries are periodic, the upper boundary is outflow, and the lower boundary is inflow, with the inflow velocity set to the laminar flame speed. This keeps an unperturbed flame stationary on our grid. The solutions presented here were computed with an advective CFL number of $`0.225`$ and typically required 8000 timesteps. For all the runs we present, adaptive mesh refinement was used, with a base grid plus one finer level that has twice the resolution. This allows us to focus the resolution on the interesting parts of the flow, allowing larger domains to be studied. Refinement triggered on the temperature gradient and the gradient of carbon mass fraction, ensuring that the flame’s reaction zone was always at the finest resolution.

## 4 RESULTS

Table 2 lists the parameters used for all the simulations we present. For the lowest density case, $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$, we actually present results for three different domain sizes and two resolutions and for the $`10^7\mathrm{g}\mathrm{cm}^3`$ case, we present an additional simulation with the burning disabled. Except where noted otherwise, all discussion of the $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ results refer to the widest domain run. In the discussion below, we first concentrate on the change in character of the burning as the density is changed. Next we look at integral quantities that help characterize the flame, followed by a measurement of the RT growth rate and a comparison of the RT instability with and without burning. Finally, we consider resolution and domain size effects on the results we present.

### 4.1 Varying the Density and the Transition to Distributed Burning

Figures 1 to 3 show the results for the three densities we consider in this paper, $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$, $`10^7\mathrm{g}\mathrm{cm}^3`$, and $`1.5\times 10^7\mathrm{g}\mathrm{cm}^3`$. The fuel is at the bottom of the domain (and colored red), and the ash is at the top of the domain. Gravity points up, toward increasing $`y`$, in all simulations. For all runs, we set $`g=10^9\mathrm{cm}\mathrm{s}^2`$. This value is appropriate for the outer regions of the star, and assumes that some pre-expansion of the white dwarf has taken place. The domain width was chosen based on the flame width and fire-polishing length (and considerations of how many zones we can reasonably compute) to contain several unstable wavelengths. We see from the figures that as the density increases, the small scale structure is greatly diminished, as the burning is more effective in containing it. The calculations were run until the instability reached the top or bottom of the computational domain, or the size of the mixed region grew to be much larger than the width of the domain, leading to saturation of the instability which will be discussed later. In the descriptions below, ‘bubble’ refers to the hot ash buoyantly rising into the cool fuel, and ‘spike’ refers to the cool fuel falling into the hot ash.

The results for $`\rho =6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$, 768 cm wide simulation (Figure 1) shows the flow quickly becoming mixed and the initial perturbations rapidly forgotten. A mixed region 500–1000 cm wide forms, separating the fuel and ash. Some smaller features begin to burn away as the flame evolves and the ash gets entrained by the fluid motions. At late times, it becomes impossible to draw a well defined flame surface separating the fuel and ash, as the RT mixing clearly dominates. In §4.4, we look at the effect of different domain widths on this growth.

The $`10^7\mathrm{g}\mathrm{cm}^3`$ run (Figure 2) shows a balance between the RT growth and the burning. Spikes of fuel push well into the hot ash, forming well defined mushroom caps. Before these caps can travel too far into the ash, they burn away. This is seen for several spikes in Figure 2—the mushrooms gradually burn away, almost uniformly, turning yellow/green, and finally leaving a ephemeral outline in the ash as the burning completes. As the instability evolves, we see longer wavelength modes beginning to dominate, but they also begin to burn away. We are prevented from watching the final dominant mode that grows burn away because of it reaching the top of our domain.

The highest density case, $`1.5\times 10^7\mathrm{g}\mathrm{cm}^3`$ (Figure 3) behaves qualitatively different than the two lower densities. A defining characteristic of this density is that the flame front is always sharply defined—a trait that we expect to carry over to higher densities as well. At this density, the burning proceeds rapidly enough that it can suppress the RT instability significantly. This is as expected, since the fire-polishing length is larger than the flame width at this density (see Table 1). As the spikes of fuel extend into the hot ash, we see them rapidly burning away from the outside inward. The RT mushroom caps hardly have time to form before their ends are burned away, giving them a hammerhead-like appearance.

Figure 4 shows the carbon destruction rate, $`|\dot{\omega _C}|`$, at the midpoint of each calculation for the three densities. At the low density end, the burning is concentrated in small regions whereas at the high end, the burning is smoother. We would expect this trend to continue for densities outside the range we consider here. At higher densities, the burning will dominate even more, and we would expect the flame surface to be well defined on the small scales. At densities lower than $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$, the mixed region will grow even larger, as the burning becomes even less effective at suppressing the RT instability. We will look at how the width of the reactive region grows in proportion to that of the mixed region at these low densities in §4.3.

### 4.2  Wrinkling and Acceleration of the Flame

Having described the overall evolution of the flames, we now present some quantitative measures of flame behavior. We define the effective flame speed, $`V_{\mathrm{eff}}`$, in terms of the carbon consumption in the domain; namely,

$$V_{\mathrm{eff}}(t)=\frac{_\mathrm{\Omega }\rho \dot{\omega }_C𝑑𝐱}{W(\rho X_C)^{\mathrm{in}}}$$
(7)
where $`\mathrm{\Omega }`$ is the spatial domain of the burning region, $`W`$ is the width of inflow face, $`(\rho X_C)^{\mathrm{in}}`$ is the inflow carbon mass fraction, and $`\rho \dot{\omega }_C`$ is the rate of consumption of carbon due to nuclear burning. Accurate direct evaluation of this integral is problematic because of the operator split treatment of reactions in the numerical method. A robust and accurate estimate of $`V_{\mathrm{eff}}`$ can be obtain by integrating the conservation equation for carbon mass fraction $`(\rho X_C)/t+(\rho uX_C)=\rho \dot{\omega }_C`$ over $`\mathrm{\Omega }`$ and a time interval $`[T_1,T_2]`$ to obtain

$$V_{\mathrm{eff}}=\frac{_\mathrm{\Omega }(\rho X_C)|_{T_1}^{T_2}d𝐱}{(T_2T_1)W(\rho X_C)^{\mathrm{in}}}u^{\mathrm{in}}$$
(8)
$`u^{\mathrm{in}}`$ is the inflow velocity. Note that this formula is only valid when essentially all carbon is being consumed in the domain.

Another quantity of interest is the area of the flame, which for our two-dimensional studies is a length. Computing a flame length for the lower density cases is problematic since there is often not a distinct flame surface. In this paper we define a flame length as the number of zone edges where the carbon mass fraction passes through 0.25. We normalize the length to the initial flame length so we can measure the growth of the burning surface. This approach, while crude, is robust and defines a reasonable method for computing the length of a well mixed flame where other approaches are infeasible. In particular, this approach systematically overestimates the flame length; however, this can be compensated for by normalizing by the initial flame length. We have validated that for well-defined flames, this method agrees well with the length of a contour as computed with the commercial IDL package, after normalization.

Figures 5 to 10 show the speeds and lengths of the flames at the three densities as a function of time (the additional curves on the $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ plots will be discussed in §4.4). The maximum speedups we observe are 3 times the laminar speed for the $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ run, 5 times for the $`10^7\mathrm{g}\mathrm{cm}^3`$ run, and 2.5 times for the $`1.5\times 10^7\mathrm{g}\mathrm{cm}^3`$ run. These speedups are much larger than any seen in our companion study on the Landau-Darrieus instability (paper I), as expected, based on the outcome of the large scale simulations of the explosion. The flame length increases by a factor of 50 throughout the simulation for the two lowest densities, but only increases to by a factor of 8 for the $`1.5\times 10^7\mathrm{g}\mathrm{cm}^3`$ case, due to the strong moderation of the RT instability by the burning.

Comparing the velocity curves to the length curves, we see that, in general, the velocity is not strictly proportional to the flame surface area. This is not unexpected, since the stretch and curvature that the flame experiences as well as the complicated motions of the RT instability pushing the fuel and ash together will modify the local burning rate. Evidence of this variability is seen in Figure 4 where considerable variation in the carbon destruction rate is observed, even for the highest density case. To further quantify this effect, in Figures 11 to 13 we plot the velocity divided by flame surface area normalized to $`v_{\mathrm{laminar}}/W`$, where $`W`$ is the width of the box. With the caveat that there are ambiguities in defining flame surface area for the “distributed” flames, we observe that, after an initial transient, the flame speed versus area relaxes to an essentially statistically stationary value. For the highest density case, the effective flame velocity is approximately 30% of the laminar flame speed $`\times `$ the area enhancement. For the lower density cases, the value drops to approximately 10%. Thus, at these densities, we find

$$v(t)=cv_{\mathrm{laminar}}\frac{A(t)}{A_0},$$
(9)
where $`c`$ is a proportionality constant not equal to 1. This proportionality constant increases with density, reflecting the fact that the fire-polishing wavelength, which sets the smallest scale on which we can wrinkle the flame, grows as well. Thus at higher densities, the localized curvature is smaller, and as a result, the localized velocity is less affected. This has strong implications for subgrid models, particularly in the flamelet regime, where it is normally assumed that the velocity scales directly as the increase in area of the flame surface. It appears that this is not the case, and is further complicated by the density dependence of the proportionality constant.

Woosley (1990) introduced a fractal model to describe the growth of the flame surface in a supernova when the burning is in the flamelet regime. The self-similar range through which a fractal description applies is bounded by the fire-polishing length on the small scales and the Sharp-Wheeler description for the RT instability growth on the large scales. In two-dimensions, a fractal model for the growth of the flame surface would be

$$L=L_0\left(\frac{\lambda _{\mathrm{max}}}{\lambda _{\mathrm{min}}}\right)^{D1}=L_0\left(\frac{\alpha g_{\mathrm{eff}}^2(tt_0)^2}{4\pi v_{\mathrm{laminar}}^2}\right)^{D1},$$
(10)
where we expect $`1<D<2`$ in the self-similar scaling regime. The offset in time, $`t_0`$, is chosen such that $`L=L_0`$ at $`t=0`$. The value of $`\alpha `$ is taken to be $`0.05`$. As we will see in §4.3, the Sharp-Wheeler model may not be the proper description of the reactive RT instability. With that caveat, Figures 8 and 10, we plot the predictions of the fractal scaling model using $`D=1.5`$ and $`1.7`$, which give area growing as $`t`$ and $`t^{1.4}`$ respectively. From these figures, it seems that, if this fractal model is correct, the growth in the flame surface is closely represented by a fractal dimension of 1.7. This agrees with the value computed from two-dimensional RT studies by Hasegawa et al. (1996). Some evidence suggests that the fractal dimension of RT induced turbulence may increase with time (Dimotakis et al., 1998). We caution however that we only span about a decade between $`\lambda _{\mathrm{min}}`$ and $`\lambda _{\mathrm{max}}`$ in these simulations. If such a scaling holds, it could serve as a subgrid model for large scale simulations, however, the evolution of the proportionality constant in Equation (9) needs to be understood. Additional calculations that encompass a larger range of spatial scales are needed to further validate this model.

Not all of the energy release by the flame goes into driving the expansion of the star. The RT unstable flame and associated Kelvin-Helmholtz instabilities also generate turbulence. Turbulence diagnostics can be difficult to define—here we use an integral quantity, the Favre average turbulent kinetic energy which plays a role in the $`kϵ`$ subgrid model for turbulence. The Favre average turbulent kinetic energy (for a discussion of which see Peters 2000), can be written

$$\stackrel{~}{k}(y,t)=\frac{1}{2}(\rho 𝐯\rho 𝐯)^2/\rho ^2$$
(11)
where $``$ denotes the horizontal spatial average of the quantity enclosed in the brackets. Figure 14 shows $`\stackrel{~}{k}`$ as a function of height at several instances in time for the $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ flame. We see that, within the mixed region, the turbulent kinetic energy is roughly constant, and this plateau rises quadratically with time, as shown in Figure 15. The integral of $`\stackrel{~}{k}(y,t)`$ over the vertical extent of the domain is the total turbulent kinetic energy $`k(t)`$. The generation of turbulent kinetic energy as a function of time is shown for the three densities in Figures 16 to 18. The data suggest that the increase in kinetic energy as a function of time is well approximated by a power law of the form

$$k(t)=at^b,$$
(12)
are shown in Figures 16 to 18, and summarized in Table 3. The power law fits, summarized in Table 3 match the data quite well, giving an exponent, $`b=2.3`$$`2.8`$, depending on the density, with the power increasing with increasing density. In the lowest density run, there is considerable turbulent energy, and, as Figure 14 shows, $`10^{11}\mathrm{erg}\mathrm{g}^1`$ is reached after $`6.4\times 10^3`$ s. This is on a scale of $`O(10^3)`$ cm. Niemeyer & Hillebrandt (1995b) presumed pre-existing turbulent kinetic energy of this level on scales of $`10^6`$ cm in their main calculation. Assuming a one-third power scaling of the energy cascade, on their length scales the RT generate turbulence already overwhelms this pre-existing turbulence, and it is continuing to grow, apparently quadratically with time. These results suggest that the RT instability alone can provide the turbulent motions presumed to exist in the star.

### 4.3  The Growth Rate of RT Unstable Flames

In the absence of any reactions, the extent of the mixed region should increase as $`t^2`$, as predicted by the Sharp-Wheeler model, Equation (1). This relation is only expected to hold during the phase where bubbles are merging, which excludes the initial linear growth range. Also, because we are resolving the flame structure, the interface between the fuel and ash is not infinitesimally thin, which complicates the definition of the Atwood number. Furthermore, at late times, once the bubble merging in our domain has stopped and the size of the mixed region rivals the width of the domain, this relation will also break down. This relation was derived for the purely hydrodynamical RT instability. Burning can have large consequences here. The Sharp-Wheeler model describes the growth of this mixed region, which will contain fuel and ash, but that fuel is burning, and as it turns into ash, we would expect to find the mixed region smaller than that of the purely hydrodynamical case. For the present simulations, since we have a limited number of modes in our box (typically 10), we are only going to be able to do the fits on a small subset of the data. A future study will focus on scaling in the flamelet regime, using much larger domains.

Figures 19 to 21 show the extrema of the mixed region, computed by laterally averaging the carbon mass fractions and finding the positions where it first exceeds 0.05 and 0.45. This definition is consistent with that used in Khokhlov (1995). The curves are measured with respect to the initial position of the interface. The top curve measures the position the spikes of fuel pushing into the region of hot ash. The lower curve is the position of the bubbles of hot ash floating into the fuel. The sharp kinks in the spike curves represent the instances where the plume of fuel, having pushed far into the ash region, burned away. At the highest density (Figure 21), the bubble curve is very smooth.

We can make fits of the growth of the mixed region to time by fitting to a power law,

$$w(t)=ct^n.$$
(13)
We do not attempt to fit to the Sharp-Wheeler scaling and find a value of $`\alpha `$ because of the limited range of wavelengths we follow, and the uncertainty as to whether that relation holds in the reactive case. We exclude the initial linear phase of the growth of the RT instability in our fits. Figures 22 to 24 shows the results. In all cases, we find $`n`$ less than 2.0, indicating that we are not in the Sharp-Wheeler regime. Furthermore, $`n`$ increases with increasing density. We find $`n=1.16`$ for $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$, $`n=1.2`$ for $`10^7\mathrm{g}\mathrm{cm}^3`$, and $`n=1.55`$ for $`1.5\times 10^7\mathrm{g}\mathrm{cm}^3`$. Larger scale studies, with a greater range of wavelengths are needed to determine whether this is a general result for reactive Rayleigh-Taylor, or if this is because of the small size domains considered here. One way we can get some insight into this is to rerun one of these simulations without burning. (Note that in this case, we shift to a non-moving frame of reference whereas the reacting cases are performed in a reference frame moving at the laminar flame speed.) The growth rate for the $`10^7\mathrm{g}\mathrm{cm}^3`$ run with burning disabled is over-plotted in Figure 20. We see that it is much smoother than the corresponding reactive case, since the spikes of fuel which push into the hot ash never burn away. The only wiggles in the mixing region growth curves result from bubble mergers. The fit to a power law for this non-reactive simulation is presented in Figure 25. The carbon mass fraction for several points in time are shown in Figure 27, which compares directly to Figure 2. Here, $`n=1.36`$, higher than the reactive case, but still not equal to 2.0. The difference between this value and the reactive case shows that the burning does have some influence. In the non-reactive case we really do expect $`n=2.0`$, and therefore the difference we see is very likely because we also are probably not following enough modes to see the Sharp-Wheeler scaling.

For the $`1.5\times 10^7\mathrm{g}\mathrm{cm}^3`$ run, we were closest to growing as $`t^2`$. The burning here is quite vigorous, compared to the RT growth, so as argued above, we would expect to see its influence on the extent of the mixed region. However, we would expect it to affect the position of the spikes of fuel only, since the bubbles of ash do not react. Therefore, we can fit the growth of the bubbles and spikes from the initial interface separately. This is shown in Figure 26. Again, we exclude the initial linear growth phase of the instability. We find the bubble position growing as $`n=1.85`$ and the spike region growing as $`n=1.28`$. Based on the trends we see here and above, we would expect that once we have enough unstable modes in the box, the width of the mixed region during the bubble merger phase will scale as slightly smaller $`t^2`$, with the bubble position itself growing as $`t^2`$. Therefore, more studies, in larger domains with more unstable modes are needed.

It is interesting to look at how the size of the reactive region scales with that of the mixed region, especially at the lowest density. As discussed in §2, if the reactive region grows to the size of the detonation matchhead, then it may be possible for a deflagration-detonation transition to occur. Measuring the scaling of the reactive region to the mixed region is not easy, because the flame is so wrinkled. We will use a volume weighted definition here. The collection of zones, $`A`$, where the carbon mass fraction, $`X_C`$, is between $`\eta _{\mathrm{min}}`$ and $`\eta _{\mathrm{max}}`$ is

$$A\left\{X_{C(i,j)}\right|\eta _{\mathrm{min}}X_{C(i,j)}\eta _{\mathrm{max}}\}.$$
(14)
Similarly, we can define the collection of zones, $`B`$, where the density-weighted carbon destruction rate, is within $`\gamma `$ of the peak,

$$B\left\{\left(\rho \dot{X}_C\right)_{(i,j)}\right|\left(\rho \dot{X}_C\right)_{(i,j)}\gamma \mathrm{max}\left(\left(\rho \dot{X}_C\right)_{(i,j)}\right)\}.$$
(15)
Then, the ratio of the size of the reactive region to the mixed region is

$$\mathrm{\Gamma }=\frac{V_B}{V_A},$$
(16)
where $`V_A`$ is the volume of the zones in set $`A`$. Once the mixed region develops, we expect this quantity to be always less than unity. Because of the freedom to choose the parameters $`\eta _{\mathrm{min}}`$, $`\eta _{\mathrm{max}}`$, and $`\gamma `$, the numerical value of $`\mathrm{\Gamma }`$ itself does not have much meaning, but the trend with time does. We pick $`\eta _{\mathrm{min}}=0.05`$, $`\eta _{\mathrm{max}}=0.45`$, and two values of $`\gamma `$: 0.1 and 0.8. We note that the trend is relatively insensitive to the choice of these parameters, as illustrated by Figure 28, which shows $`\mathrm{\Gamma }`$ as a function of time for the $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ RT unstable flame for the two choices of $`\gamma `$. After an initial transient, it appears that these curves level off, indicating that the reactive region continues to growth in size proportional to the mixed region as the flame evolves. We note that this shows that the reactive region is a small fraction of the mixed region, but this absolute scaling is dependent on the choice of the three parameters. This is illustrated in Figure 29, where the mixed and reactive regions are shown, midway through the calculation. At late times, the peak nuclear energy generation comes from regions where the carbon mass fraction is $`0.15`$, as shown in Figure 30. This is also apparent by looking at the laterally averaged flame profile and comparing to the laminar state—see Figure 31. If this reactive region can growth to a critical matchhead size, before the entire star is consumed, it may be possible for a deflagration to detonation transition to proceed. However, since at late times the peak burning is occurring at mass fractions of 0.15 rather than 0.5, this will greatly increase the required matchhead size, further complicating the possibility of a deflagration-detonation transition.

### 4.4  Resolution and Domain Size Dependence on the Results

In this subsection, we examine the role of domain size and resolution on the computational results. Figures 32 and 33 show the evolution of the $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ RT simulation in narrower domains (96 cm and 384 cm wide respectively). Other than the domain size, the parameters are identical to those used in the results described shown in Figure 1. In the narrowest case, once the size of the mixed region reaches a fixed size, comparable to the width of the domain, the mixing and reacting processes saturate and the flow enters a quasi-steady-state, with vortices that are remnants of the early mixing driving transverse shear flow in the mixed region. Once this saturation is reached, the burning rate begins to slowly decay. This behavior is not seen in the corresponding wider domain calculation, although, we would expect a similar pattern if we use a taller domain and ran out to longer times.

Figure 5 shows three additional curves corresponding to these narrower domains, the narrowest at two different resolutions (see Table 2). In the 96 cm wide domain, the peak velocity is much smaller, reaching only about 3 times the laminar speed. This is to be expected from Equation (2)—the wider domain allows longer wavelength modes to go unstable. The decrease in velocity at late times in the narrow domain begins when the size of the mixed layer becomes comparable to the width of the domain, preventing any further modes from growing. This is reflected in the later panels of Figure 32. In the 384 cm wide domain, the velocity continues to grow, beyond 6 times the laminar speed—we would expect our widest domain run to continue to accelerate as well, but it began to interact with the top of our computational domain and was not run out as long.

The velocity of the low resolution run tracks that of the high resolution case very closely, suggesting that, for the large scale diagnostics, we have converged. This lowest resolution run has approximately 5 points in the thermal width (Equation 6), a value found acceptable in the convergence study performed in Bell et al. (2003b) for Landau-Darrieus unstable astrophysical flames. We note again that this definition of thermal width is smaller than some other commonly used definitions (see, for example Dursi et al. 2003). Finally, in all calculations, the initial perturbations are well resolved, with typically 50–100 zones per wavelength (with the only exception being the 96 cm wide, $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ run, where there are only 20 zones per wavelength). This zoning is at or exceeds the resolution determined to give acceptable convergence of single-mode RT growth rates in the 3-d compressible study presented in Calder et al. (2002).

We also can look at how the scaling of the turbulent kinetic energy depends on the domain size. For the 384 cm wide, $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$ run, this is shown in Figure 34. Two fits are shown. When all of the data is included, the fit is quite poor, since it is biased by by data at the end of the run when the mixing has saturated and the turbulent kinetic energy stops increasing. Excluding this data, the fit to the first $`0.005`$ s is quite good, and with the kinetic energy $`t^{2.325}`$ quite close to the 2.267 power obtained from the 768 cm wide run. Thus, it appears that the scaling of the turbulent kinetic energy prior to saturation is insensitive to the domain size.

## 5 DETONATION IN THE DISTRIBUTED REGIME?

Based upon the observed systematics of flame propagation in the distributed regime, we can begin to speculate on the potential for a transition to detonation. As discussed by Niemeyer & Woosley (1997) and Niemeyer (1999), there are two possible modes for a delayed transition to detonation: a) a “local” transition to detonation because a fluid element of critical mass burns with a supersonic phase velocity and b) a “macroscopic” transition because some appreciable fraction of the white dwarf volume develops such complex topology that burning briefly consumes more fuel than could a spherical front encompassing that region and moving at supersonic speed. The former, also known as the Zeldovich mechanism, is the basis for the “delayed detonation model” by Khokhlov (1991); the latter is the basis for a similar model by Woosley & Weaver (1994). Niemeyer (1999) has argued that the volume detonation model requires special preconditioning, but because we study only the small scales on which the flame is resolved, we cannot comment in a meaningful way on this issue. Also, because we have not included turbulence cascading down from scales larger than our grid, we cannot conclusively argue about the local transition to detonation either (a large amount of turbulence could, for example, induce a transition to distributed burning at a higher density). Still some scaling relations are observed that, if they can be generalized by larger scale studies, argue against a transition to detonation.

As is well recognized (e.g. Niemeyer, 1999), transition to detonation can never occur so long as a well defined flame exists. The width of a steady flame is always thinner than the critical mass required for detonation. However, a mixture of hot fuel and cold ash is potentially explosive as calculations in the SN Ia context (Lisewski et al., 2000) have suggested. Figure 4 shows the existence, in the distributed regime, of localized hot spots of unsteady burning. The high temperature sensitivity of the 12C + 12C reaction itself makes it difficult for one of these hot spots to become supercritical in mass. For the conditions we consider, carbon fusion dominates the energetics, and the energy production is

$$S_{\mathrm{nuc}}X_C^2T^n,$$
(17)
where $`X_C`$ is the mass fraction of carbon, and $`n`$ given by the Coulomb barrier between two carbons nuclei,

$$n=\frac{\tau 2}{3},$$
(18)
$`\tau `$ $`=`$ $`4.248\left({\displaystyle \frac{Z_1^2Z_2^2\overline{A}}{T_9}}\right)^{1/3}`$ (19)
$`=`$ $`84.13T_9^{1/3},`$ (20)
with $`T_9`$ the temperature in 109 K. For our calculations with $`\rho `$ = $`6.67\times 10^6\mathrm{g}\mathrm{cm}^3`$, the temperature of the hot ash is $`T_9=2.4`$. Anticipating that the temperature of a combustible fuel-ash mixture will not be far from that, we find $`n20`$. The carbon mass fraction in the fuel is 0.5 and in the ash it is zero, and the temperature of the unburned fuel is negligible.

In any mixed mass, $`M`$, composed of fractions $`f`$ of ash and $`1f`$ of fuel, the temperature will be

$$T_{\mathrm{mix}}=f\frac{C_\mathrm{P}(T_{\mathrm{ash}})}{C_\mathrm{P}(T_{\mathrm{mix}})}T_{\mathrm{ash}}$$
(21)
and the carbon mass fraction

$$X_{C,\mathrm{mix}}(1f)X_{C,\mathrm{fuel}}.$$
(22)
The dependence of the heat capacity on temperature is important and mitigates somewhat the extreme sensitivity of the reaction rate. In the vicinity of $`T_9`$ = 2.5, the heat capacity is predominantly due to the electron gas (with $`C_{\mathrm{P},\mathrm{e}}T`$), but with a non-negligible contribution from radiation ($`C_\mathrm{P}T^3`$) (Woosley et al., 2003). To good approximation, near $`T_9`$ = 2.5, $`C_\mathrm{P}T^{7/4}`$.

The maximum in energy generation will then occur for

$$\frac{d}{df}[(1f)^2f^{80/11}]=0,$$
(23)
or for $`f40/51=0.784`$—that is, in the burning mixture, $`X_C=0.11`$ and $`T_9=2.20`$. A factor of two variation in the energy generation occurs for $`0.60<f<0.91`$, or $`X_C`$ from $`0.045`$ to $`0.20`$. This agrees well with what is observed in the numerical simulation (Fig. 30).

Coupled with the fact that burning to magnesium or silicon releases less energy than burning to the iron group, this low carbon mass fraction implies that the overpressure from burning will be small in the regions where ash and fuel are mixed. Burning from $`T_9=2.2`$ to $`2.4`$, the typical ash temperature, will only raise the pressure by 8%. Since this overpressure determines the temperature reached in a detonation and the burning rate depends on this to a high power, the critical mass for initiating a detonation will be large. The size will certainly be greater than the distance a sound wave can travel during the time it takes the mixture to burn. At $`T_9`$ = 2.2 and $`X_C=0.11`$, the nuclear time scale is

$$\tau _\mathrm{T}\frac{C_\mathrm{P}}{\dot{S}_{\mathrm{nuc}}},$$
(24)
or about 0.07 s. The speed of sound is approximately 5000 km s-1. Hence a critical mass for these conditions must be at least 300 km in size.

On the other hand, Figure 28 shows that, in steady state, the size of the reactive region is approximately a small constant times the size of the mixed region. Even allowing a range of burning time scales of a factor of 10 across the burning region, the explosive mixture only has dimensions about 10% that of the region that has been mixed by the RT instability. Taking $`0.05g_{\mathrm{eff}}t^2`$ as an upper bound to that (see §4.3), and realizing that the temperature will decline still further in a few tenths of a second because of expansion, the largest mixed, potentially explosive region is less than 5 km. Thus, though a transition to detonation is in principal possible, there simply may not be room, nor time enough to develop one in an exploding white dwarf.

## 6 CONCLUSIONS

We presented direct numerical simulations of reactive RT instabilities in conditions appropriate to the late stages of a Type Ia supernova explosion. Fully resolving the flame frees us from the need to specify any model parameters, such as the flame speed, Markstein length, etc. In the density range we consider, the flame transitions from having a distinct, well defined interface at the high density end to being a chaotic, well mixed burning region at the lower density end, with a transition density of $`10^7\mathrm{g}\mathrm{cm}^3`$. This transition is expected based on the arguments presented in §2. Furthermore, this would suggest that at the even higher densities characteristic of the early stage of the explosion, the flame surface continues to be well defined, consistent with the conclusions of Niemeyer et al. (1999) and the 2-d flame-model simulations at $`10^8\mathrm{g}\mathrm{cm}^3`$ presented in Khokhlov (1994). The effects of the burning were further made concrete by presenting the $`10^7\mathrm{g}\mathrm{cm}^3`$ case both with and without reactions.

At the lowest density, the mixed region becomes very large, although, still smaller than the critical carbon detonation matchhead size (Niemeyer & Woosley, 1997). It appears that the reactive region grows in direct proportion to the size of the mixed region, making it possible in theory that such a transition can occur, however, at these low densities, much of the star is already consumed, and, based on scaling predictions from the present simulations, the time required for the reactive region to grow to the matchhead size seems to be longer than the total explosion time. In effect, the flame runs out of star before a critical mass can be built up. Future simulations will explore this in more depth.

In all cases presented here, the flames accelerated considerably, reaching speeds between 2 and 6 times the laminar flame speed. This is significantly larger than accelerations seen in small-scale/resolved Landau-Darrieus instability studies (paper I). The growth of the flame surface appeared to be well described by a fractal model, with a fractal dimension of $`1.7`$. In the cases considered, the effective flame speed was asymptotically proportional to the increase in flame surface area but the constant of proportionality was about 0.1 for the low density cases and 0.3 for the higher density case. Thus, the flame acceleration was considerably less than would be predicted by the increase in flame surface area alone. This proportionality constant may approach unity as the density is increased, because the smallest scale on which the RT instability can bend the flame (the fire-polishing length) increases dramatically with density, and therefore the magnitude of the curvature that the flame experiences locally decreases. We would expect that on larger domains the flame would continue to accelerate, further supporting the already widely held view that the RT instability provides most of the acceleration to the flame in a pure deflagration Type Ia SNe explosion.

The growth of the mixed region for the reactive RT instability appears to be slower than the non-reactive case. The present results do not seem to support the Sharp-Wheeler model, but larger scale studies are needed, with more bubbles merging over a longer period of time to better understand the growth of the mixed region. Understanding the growth of the mixed region is critical to providing an accurate subgrid model.

Future studies will focus on the three-dimensional counterparts to the flames we discussed here. Turbulence behaves very differently in two and three dimensions, so we expect to see some differences in the evolution of the instability. Numerical results indicate that the growth rate of the pure RT instability is faster in 3-d than in 2-d (Kane et al., 2000; Calder et al., 2002). Furthermore, in 3-d the surface to volume ratio of the spikes of fuel is larger than in 2-d, exposing more fuel to the hot ash, so we may expect them to burn away more quickly. We still expect, based on the arguments presented in §2 and the results of the simulations, that the flame will undergoing the same transition in behavior as the density decreases. Larger two-dimensional studies in the flamelet regime are also needed, moving toward higher densities. Capturing a greater range of scales on the grid will allow for validation of the fractal scaling model and a better understanding of the moderation of the Sharp-Wheeler prediction for the growth of the mixed region. Finally, these fully resolved simulations can serve as the basis for testing various models to represent the flame on subgrid scales, which become necessary to address the larger scale Type Ia physics.

The authors thank F. X. Timmes for making his equation of state and conductivity routines available online. Support for this work was provided by the DOE grant No. DE-FC02-01ER41176 to the Supernova Science Center/UCSC and the Applied Mathematics Program of the DOE Office of Mathematics, Information, and Computational Sciences under the U.S. Department of Energy under contract No. DE-AC03-76SF00098. SEW acknowledges NASA Theory Award NAG5-12036. Some calculations were performed on the IBM SP (seaborg) at the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the DOE under Contract No. DE-AC03-76SF00098, the IBM Power4 (cheetah) at ORNL, sponsored by the Mathematical, Information, and Computational Sciences Division; Office of Advanced Scientific Computing Research; U.S. DOE, under Contract No. DE-AC05-00OR22725 with UT-Battelle, LLC, and the UCSC UpsAnd cluster supported by an NSF MRI grant AST-0079757.