Distributed Burning in Type Ia Supernovae: A Statistical Approach A.M. Lisewski, W. Hillebrandt Max-Planck-Institut fu¨r Astrophysik, Karl-Schwarzschild-Str. 1, 85740 Garching, Germany 0 0 S.E. Woosley 0 2 UCO/Lick Observatory,University of California Santa Cruz, Santa Cruz, CA 95064,USA n a J J.C. Niemeyer 4 2 University of Chicago, Department of Astronomy and Astrophysics, 5640 S. Ellis Avenue, Chicago, IL 2 60637,USA v 8 and 0 5 9 A.R. Kerstein 0 9 Combustion Research Facility, Sandia National Laboratories,Livermore,CA 94551-0969,USA 9 / h p - Received ; accepted o r t s a : v i X r a – 2 – ABSTRACT We present a statistical model which shows the influence of turbulence on a thermonuclear flame propagating in C+O white dwarf matter. Based on a Monte Carlo description of turbulence,itprovidesamethodforinvestigatingthephysicsintheso-calleddistributedburning regime. Using this method we perform numerical simulations of turbulent flames and show that in this particular regime the flamelet model for the turbulent flame velocity loses its validity. In fact, at high turbulent intensities burning in the distributed regime can lead to a deceleration of the turbulent flame and thus induces a competing process to turbulent effects that cause a higher flame speed. It is also shown that in dense C+O matter turbulent heat transport is described adequately by the Peclet number, Pe, rather than by the Reynolds number, which means that flame propagation is decoupled from small-scale turbulence. Finally, at the onset of our results we argue that the available turbulent energy in an exploding C+O white dwarf is probably too low in order to make a deflagration to detonation transition possible. Subject headings: methods: statistical – nuclear reactions – stars: supernovae: general – turbulence – 3 – 1. Introduction The thermonuclear explosion of a C+O Chandrasekhar-mass white dwarf, which is believed to be the underlying process of a type Ia supernova (SN Ia), has been subject to numerous investigations. However, despite the fact that there are different plausible models explaining the history of the explosion (Nomoto et al. 1984; Woosley & Weaver 1986; Mu¨ller & Arnett 1986; Livne 1993; Khokhlov 1991a; Arnett and Livne 1994; Niemeyer & Hillebrandt 1995; Ho¨flich, Khokhlov & Wheeler 1995; Ho¨flich 1995; Wheeler et al. 1995; Ho¨flich & Khokhlov 1996), many important details still remain unclear. Thermonuclear reactions provide the source of energy which possibly unbinds the white dwarf. Thus, a SN Ia is characterizedby the physics of thermonuclearflames that propagatethroughthe star. The physicalconditions ofthis flame vary drastically during the different temporal and spatial stages of this process. A better understanding of these conditions, their interaction with the flame and finally their consequences regarding the explosion itself are still important issues from the theoretical point of view. Here, we focus on the interaction between thermonuclear burning and turbulence taking place during the burning process. Turbulence is caused by different kinds of instabilities (like shear instabilities or the Rayleigh-Taylor instability) that occur on certain length and time scales (for an overview, see Niemeyer & Woosley 1997). Rough estimates give a turbulent Reynolds number Re 1014 at an integral scale of ≈ L 106cm (Hillebrandt & Niemeyer 1997). Consequently, the Kolmogorov scale l , i.e. the scale where k ≈ microscopic dissipation becomes important, is about 10−4cm. This enormous dynamical range makes a directrepresentation– atleastbymeans ofnumericalmethods –practicallyimpossible. Onthe otherhand, turbulence is a characteristic feature during the explosion process and it must be consideredin any realistic model of a SN Ia. The simultaneous coupling to energy generation due to nuclear reactions leads directly to the physics of turbulent combustion, where a few crucial problems, even for terrestrial conditions, still remain unsolved. These problems include the prescription of the effective turbulent flame speed or the existence of different modes in turbulent combustion along with their physical properties. Global features of turbulent combustion can be systematically classified by a small set of dimensionless parameters. This classification, which already found a wide utilization among the chemical combustion community, can be used in the field of SNe Ia in order to emphasize universality in the physics of turbulent combustion. In this work we study the properties of turbulent flames in a certain state, the so-called distributed flame regime. The physical characterizationof the latter is given by the situation where turbulent motions are fast enough to disturb the flame on microscopic scales. Therefore a distributed flame does locally not – 4 – look like a laminar flame anymore which is a basic difference to turbulent flames in the flamelet regime. Turbulent burning fronts in an exploding white dwarf are flamelets at high and intermediate densities, but below 5 107g cm−3 onecannotexpectthatthe flameletpictureisstillvalid(Niemeyer&Woosley1997). ∼ × Until now, in the context of SNe Ia, burning in the distributed regime has only been considered qualitatively. Therefore we wish to present first quantitative results regarding this issue. Furthermore, there are good reasons to believe that this combustion mode plays an important role in the deflagration to detonation transition (DDT), which in turn is a promising model, for empirical and theoretical reasons, for the explosive stage of a SN Ia (e.g. Khokhlov1991b;Niemeyer & Woosley 1997;Niemeyer & Kerstein1997; Khokhlov et al. 1997). In order to attempt a representation turbulent dynamics in the distributed flame regime, we use a new model, formulated in one spatial dimension, which nevertheless provides essential features of three dimensional homogeneous turbulence. It consists of a statistical description of turbulent mixing and a deterministic evolution of the underlying microphysics. This method allows a systematical investigation of turbulence phenomena. When it is coupled to a nuclear reaction network, it gives first insights of how the flame structure is affected by turbulence on scales, which have not been resolved in direct numerical simulations. Inparticular,we investigatethe flame propertiesin caseswhere the Gibsonscaleis comparable to the thickness of an undisturbed conductive flame. The Gibson scale l is defined as the length scale G on which the turbulent velocity fluctuations equal the laminar flame velocity. In the case of flames in degenerate white dwarf matter l becomes comparable to the thickness of the flame only for densities G around 2 107g cm−3 and below (Khokhlov et al. 1997; Niemeyer & Woosley 1997). × The outline of the paper is as follows: first we introduce and describe the statistical method that we use to model turbulence. Then we couple the latter to physically relevant microscopicaldiffusion processes, such as temperature and viscous diffusion, and to external energy sources coming from nuclear reactions. Theresultingmethodiseventuallyusedtoinvestigatesomepropertiesofturbulentflamesinthe distributed regime, like their effective flame speed. Finally we discuss the relevance of our results to the DDT problem. 2. One-Dimensional Turbulence Since fundamental aspects of turbulence can be recovered from the knowledge of the statistical moments and correlations of the velocity flow, the statistical approach to turbulence is particulary – 5 – appealing. Therefore, the history of statistical methods in turbulence theory is rather long. In this context we present a novel model of turbulence (Kerstein 1999). It is a stochastic method, realized as a Monte Carlo simulation, which allows to compute statistical properties of the flow velocity and of passive scalars in stationary and decaying homogeneous turbulence. One-dimensional turbulence (ODT) represents many aspects of three-dimensional turbulence, but it is formulated in only one spatial dimension. This model provides the temporal evolution of a characteristic transverse velocity profile u(y,t) of the turbulent medium, where y is the spatial location on a finite domain [0,Y] and t is the elapsed time. This is done in a two-fold way: u(y,t) is subject to a molecular diffusion process and to a random sequence of profile rearrangementsrepresenting turbulent eddies. Reflecting the typical behavior of turbulence kinematics, the event rate of these profile rearrangements (so-called eddy mappings) is proportional to a locally averaged shear of the velocity profile u. Given the transverse velocity profile at a certain time, u(y,t), the mapping which models the action of an turbulent eddy of size l and at the position y reads 0 u(3y 2y ,t) y y y +f l 0 0 0 1 − ≤ ≤ uˆ(y,t)=u( 3y+4y +2l,t) y +f l y y +f l (1)  − 0 0 1 ≤ ≤ 0 2 u(3y 2y 2l,t) y +f l y y +l 0 0 2 0 − − ≤ ≤ where it is f =1 f = 1/ anduˆ(y,t)=u(y,t) for all y / [y ,y +l]. This three-valued map represents 1 2 3 0 0 − ∈ the typical features of a turbulent vortex, namely rotation and compression. Each random map defines an eddy time scale, τ(y ,l,t), via 0 l τ(y ,l,t)= . (2) 0 2A u(y ,t) u (y +l/2,t) l 0 l 0 | − | Here, u is the boxcar averaged profile of u over a length scale l. A is the only model independent and l dimensionless parameter which has to be fixed empirically. Vortical kinetic energy is fed by the kinetic energy of the local shear. Thus equation (2) can be interpreted as an energy balance. It then is 1 1 ρ (lτ−1)2 = ρ (2A u (y ,t) u (y +l/2,t))2, (3) 0 0 l 0 l 0 2 2 | − | where ρ is the fluid density. Equation (2) resp. (3) is used to introduce the statistical hypothesis of ODT. 0 It assumes that the occurrence of eddies with size l and location y (with respective tolerances dl, dy) is 0 governedby a Poisson random process with mean event rate A dydl =:λ(y ,l,t)dydl, (4) l2τ(y ,l,t) 0 0 – 6 – where the processes for different values of y and l are statistically independent. In this context A can be 0 viewed as a factor which scales the event rate λ. Motivated by the results in Kerstein (1999) we choose A= 0.23. In addition, the microscopic evolution is modelled by a diffusion process of the kind u =νu , t yy with the kinematic viscosity ν. The temporal evolution of turbulence, stationary or decaying, depends on the choice of boundary conditions for u. A decaying turbulent intensity is obtained by periodic boundary conditions, whereas the stationary case is given by the choice of jump periodic boundary conditions: u(y+Y,t)=u(y)+u , u :=u(Y,t=0) and u(y,t=0) being strictly monotonic on [0,Y]. 0 0 The numerical implementation of ODT reproduces many typical features of three-dimensional homogeneous turbulence (Kerstein 1999). For instance, power spectra of the kinetic energy show the self similar k−5/3 power law within the inertialrangedownto the scale where the transitionto dissipationtakes place, c.f. Fig. 1. An obvious advantage of this ansatz is the high spatialresolutionof turbulence compared withmultidimensionalnumericalmodels. Incombinationwiththe relativelymoderatecomputationaleffort, ODT appears as a useful tool for performing parameter studies in turbulence theory. On the other hand, ODT does not consider any pressure fluctuations (dynamical or external) in the temporal evolution of the velocity profile u. The reason for this artefact is the inherent conservation of kinetic energy in only one spatial velocity component due to equation (3). Note that the pressure gradient term in the Navier-Stokes equations redistributes energy among the different spatial components. No such redistributionis considered here. However, for isobaric flows ODT appears to be an appropriate model of turbulence. 3. Turbulent Flames in Dense C+O Matter Although the aim of this work is to investigate some properties of turbulent burning in dense C+O matter, we begin this section by recalling some properties of laminar, conductive flames. Since the physics of undisturbed flames inside a white dwarf is relatively well understood, the first reason for it is to verify established results, such as those given in the work of Timmes & Woosley (1992). The other reason is that thermal conduction and viscous diffusion will be employed as the underlying microscopical processes for ODT in order to model turbulent burning fronts in white dwarfs. The state of unburned matter, i.e. density and nuclear composition, uniquely defines the propagation velocitys ofthe conductiveflame. Timmes&Woosley(1992)calculatedtheflamevelocitiesaswellastheir l thickness for various fuel compositions and densities. At lower densities, for 108gcm−3 > ρ > 107gcm−3, ∼ ∼ the speed of the laminar flame decreases rapidly. This behavior is accompanied by a strong increase of – 7 – the flame thickness δ, which is essentially the size of the nuclear reactive zone (Timmes & Woosley 1992; l Khokhlov et al. 1997). As already mentioned in the introduction, we are especially interested in densities of the order (1 3) 107gcm−3. Thus, we set up a conductive flame propagating into unburned matter ∼ − × consisting of half 12C and half 16O at densities of ρ= 1.3 107gcm−3 and of ρ= 2.3 107gcm−3. This × × is done by solving the equations for the conservationofmole fractionsand enthalpy in planargeometry,viz. dY i = Y Y λ (i)+Y Y λ (i), (5) i k jk i k kj dt − j,k X ∂T 1 ∂ ∂T ∂ 1 1 = σ P + S˙, (6) ∂t ρc ∂x ∂x − ∂tρ c p (cid:18) (cid:19) p dY S˙ =N iB . (7) A i dt i X Herein S˙ denotes the local specific energy generation rate, σ the thermal conductivity and N is Avogadro A number and B is the nuclear binding energy of the nucleus considered. The nuclear reaction network i consists of seven species, viz. 4He, 12C, 16O, 20Ne, 24Mg, 28Si, and 56Ni. Equation (5) already makes use of the factthat the Lewis number,thatis the ratioofheatdiffusion tomass diffusion, inwhite dwarfmatter is around 107 and consequently microscopic transport of the element species can be neglected. Furthermore, in white dwarf matter the laminar flame speed is much smaller than the speed of sound, s < 0.001a . l s Therefore pressure does not change significantly across the flame front. Finally, the typicaltime scale forfree collapseorfree expansionofthe star canbe estimatedby (Fowler & Hoyle 1964) 446 τ s. (8) h ≈ ρ¯1/2 For a white dwarf of radius R = 108 cm and of mass M = 1.4M , equation (8) gives τ 0.02 s. wd wd ⊙ h ≈ Timmes & Woosley (1992) stated that gravitational influence can be safely dropped since the diffusion timescale of the flame, τ δ2/D, where δ is the flame thickness and D is the temperature diffusion d ≈ l l coefficient, is several orders of magnitude smaller than τ . However, for very low densities (ρ<107gcm−3) h ∼ these two timescales can become comparable and a possible expansion of the star could directly affect the laminar flame. To avoid additional complexity we assume in this work that the hydrodynamical timescales is much longer than the relevant timescales of the burning front. Thus we consider the small-scale burning front to be not affected by stellar expansion. Thenumericalsolutionofequations(5)and(6)isgivenexplicitlyintimecombinedwithanappropriate equation of state representingions, black body radiation,electrons and positrons. For ρ=2.3 107gcm−3, × – 8 – Figure 2 shows the propagation of the laminar flame which is indicated by the moving jump in the temperature profile. In the steady state evolution, the flame moves into the unburned matter at a velocity s = 2.1 104cm s−1. The flame thickness δ is characterized by the width of the temperature jump l l × caused by carbon burning. It is δ =1.1 cm. After 12C destruction the burned material has a temperature l 3.2 109 K and a density of 1.45 107g cm−3. For these conditions the subsequent oxygen burning has a × × destructiontimescale ofτ 0.1s. That means that as the carbonflame propagateswith a steady-state des,O ≈ velocity s, oxygen burning ignites at a distance of 2.1 103 cm behind the flame front. This distance is l × about ten times larger than our whole computational domain. Consequently, we are not able to track the whole spatialregionwhere nuclearreactionstake place thataretriggeredbycarbonburning. The extension of the reactive zone is even larger if one takes burning of heavier elements, such as silicon, into account. Thus, the resolution of the flame with all its reactions up to nuclear statistical equilibrium would require an immense spatial resolution. However,inourcasethe fuelconsists halfofcarbonandthe destructionofthe latter contributesnearly all of the total energy release by the laminar flame. In addition, the speed of the laminar flame is governed by carbon destruction. Therefore it is fair to concentrate only on the latter. We also solved the above equations for a density of 1.3 107 g cm−3 with the same nuclear composition as before. The resulting × flame velocity for this case is s =0.85 104cm s−1, while the flame’s thickness is 3.8 cm. Our results agree l × well with the ones given in the work of Timmes & Woosley (1992). As already mentioned in the introduction, the Gibson scale l can be used to measure the influence G of turbulence on a flame. Since it is directly related to the turbulent fluctuation velocity u′, the Gibson scale depends strongly on the kind of turbulence that is considered. Essential features of turbulence, like structure functions, the temporal evolution of the turbulent intensity and the energy cascade, depend on the absence or presence of certain geometrical (homogeneity, isotropy) and physical (external forces, energy sources, state of the turbulent matter) conditions. On large scales in the interior of a white dwarf, turbulence is mainly caused by the Rayleigh-Taylor instability of the flame which takes place on a length scaleofl 106cm. This instabilityproducesaturbulentenergycascadedownto dissipativelengthscales RT ≈ l (10−4 10−3)cm (Khokhlov 1995; Niemeyer & Woosley 1997). On all scales l <l l turbulence k k RT ≈ − ≪ is believed to be decoupled from gravitationaleffects and thus can be described by the Kolmogorovtheory. The latter sets the scaling law for u′ to be u′(l) l1/3. Inbetween these scale bounds nuclear burning is ∼ affected by turbulent motion and different kinds of burning regimes reveal (For an overview, see Niemeyer & Woosley 1997). Here, we put our emphasis on the special situation where homogeneous, isotropic – 9 – turbulence interacts with the microscopic structure of the laminar flame, i.e. where l becomes comparable G to δ. The relation l δ marks the transition into the distributed burning regime, where the smallest l G l ≈ turbulent vortices can enter the interior of the laminar flame and can carry away reactive material before it is completely burned (Peters 1986; Niemeyer & Woosley 1997). Taking into account the small Prandtl numbers in the dense matter of a white dwarf, we can describe this situation by a constraint on the turbulent Karlovitznumber. This number is the ratio of the diffusion timescale and the eddy-turnover-time at the smallest length scale of turbulent heat transport. The latter turns out to be the Kolmogorov scale, l , only in case of Pr 1. For Prandtl numbers much smaller than one turbulent heat transport is subject k ≈ to significant diffusion effects already at a scale of l Pr−3/4. Thus we have k u′(l )δ Ka k l Pr1/2 >1. (9) ≡ s l l k ∼ Using Kolmogorovtheory an equivalent representation of the latter condition reads u′(L)δ Ka = l Pr1/2Re1/2 s L l u′(L) 2 = Re−1/2Pr−1/2 >1. (10) s (cid:18) l (cid:19) ∼ There is a simple relation between the Gibson length and the Karlovitz number as follows l δ Ka−2. (11) G l ≈ Therefore, l < δ is equivalent to Ka > 1. The actual value Ka for a given turbulent flame plays an G l outstanding role in turbulent combustion physics. We now combine ODT with the structure of a laminar flame to show how certain properties of a nuclear flame in dense degenerate matter of white dwarf depend on the Karlovitz number. The obvious step is to incorporate equation (6) together with a viscous transport evolution for the transversal velocity profile u(x,t). Then turbulent advection is modelled by the random eddy mappings generatedby ODT, c.f. equation (1). This system reads as ∂uˆ Pr ∂ ∂(uˆ) = σ , (12) ∂t ρˆc ∂x ∂x p (cid:18) (cid:19) ∂Tˆ 1 ∂ ∂Tˆ ∂1/ρˆ 1 = σ P + S˙, ∂t ρˆcp∂x ∂x!− ∂t cp where uˆ,ρˆ,Tˆ are the rearranged profiles of the velocity, the density and the temperature according to the eddy-mapping, c.f. equation (1). – 10 – We choose the amplitude of velocity fluctuations of u to be 107 cm s−1 at 106 cm, and to ∼ ∼ obey Kolmogorov scaling. These values are motivated through the expected speed of buoyant unstable hot bubbles of size L 106 cm, which become Rayleigh-Taylor unstable and eventually give the main ≈ contribution to turbulent energy on large scales. Having fixed this velocity at a certain lengthscale, we can use Kolmogorov scaling to estimate u′(l) on smaller scales l. Thus the amplitude of turbulent velocity fluctuations in our model is comparable to the expected small-scale velocity fluctuations within a SN Ia. However, due to computational limitations we cannot consider u to be the actual fluid velocity of the white dwarf matter. Since its Prandtl number is around Pr > 10−5, the numerical resolution of both, the ∼ viscous and the heat conductive diffusion scale, would cost too much numerical effort even in one spatial dimension. To get out of this dilemma, we use an ’artificial’ Prandtl of 2 10−2 and of 10−1. This choice × underestimates the inertial range of turbulent motions within an exploding Chandrasekhar-mass white dwarf. But we believe that the physics of distributed burning in this kind of matter is decoupled from scales much smaller than l Pr−3/4, because turbulent eddies of sizes well below this limit will be smeared k instantly out by temperature diffusion. Therefore, we do not essentially change the flame properties in the distributed burning regime by increasing Pr. Of course, this is only valid for sub-unity Prandtl numbers, otherwise one possibly leaves the distributed regime. The random velocity field u is generated by stationary ODT, thus u′ becomes constant in time after an initial transient phase. For our studies we perform simulations where u′ (at a length scale l =1 cm) is equalorlargerthanagivenvalues. Thenweestimate someflame propertiessuchastheeffective turbulent l flame speed, s . So far, analytical results do not predict the function s (u′;Ka) correctly. For instance, T T experimentally observed bending effects or even complete flame quenching are still mostly unexplained problems (Ronney 1995). An overview of the various numerical models is shown in Table 1. Each realization is uniquely parametrized by its Reynolds, Karlovitz and Prandtl number. To get the turbulent Reynolds numbers, Re =u′(L)L/ν, we have estimated the integral scale L for each model from the wavenumber k for which the L relevant energy power spectrum E(k) has a maximum. Together with the r.m.s. turbulent velocity and the kinematic viscosity, ν = Pr(ρc )−1σ, we eventually obtain the Reynolds number Re and the Karlovitz p number Ka. An example of a numerical realization is illustrated in Figure 3, where a laminar and a turbulent flame are shown in a space-time diagram. Table 1 also shows the measured turbulent flame speeds s as a function of Ka, and Pr. These T