Table Of Content

Pulsating reverse detonation models of Type Ia supernovae. II: Explosion Eduardo Bravo1,2, Domingo Garc´ıa-Senz1,2, Rubén M. Cabezón1, Inmaculada Dom´ınguez3 9 0 0 2 n ABSTRACT a J ObservationalevidencespointtoacommonexplosionmechanismofTypeIasupernovaebased 0 on a delayed detonation of a white dwarf. However, all attempts to find a convincing ignition 2 mechanism based on a delayed detonation in a destabilized, expanding, white dwarf have been ] elusive so far. One of the possibilities that has been invoked is that an inefficient deflagration R leads to pulsation of a Chandrasekhar-mass white dwarf, followed by formation of an accretion S shock that confines a carbon-oxygen rich core, while transforming the kinetic energy of the h. collapsing halo into thermal energy of the core, until an inward moving detonation is formed. p This chain of events has been termed Pulsating Reverse Detonation (PRD). In this work we - present three dimensional numerical simulations of PRD models from the time of detonation o initiation up to homologous expansion. Different models characterized by the amount of mass r t burnedduringthe deflagrationphase,M ,giveexplosionsspanninga rangeofkinetic energies, s defl a K ∼ (1.0−1.2)×1051 erg, and 56Ni masses, M(56Ni) ∼ 0.6−0.8 M⊙, which are compatible [ withwhatisexpectedfortypicalTypeIasupernovae. Spectraandlightcurvesofangle-averaged 1 sphericallysymmetricversionsofthePRDmodelsarediscussed. TypeIasupernovaspectrapose v the most stringent requirements on PRD models. 3 1 Subjectheadings: Supernovae: general–hydrodynamics–nuclearreactions,nucleosynthesis,abundances 0 3 1. Introduction nation as required, for example, to determine the 1. equation of state of the dark energy component Type Ia supernovae (SNIa) are the most en- 0 of our Universe, it is necessary to understand the 9 ergetic transient phenomena in the Universe dis- physics of SNIa explosions. From the theoreti- 0 playing most of its energy in the optical band cal point of view, the accepted model of SNIa : of the electromagnetic spectrum, where they ri- v consists of a carbon-oxygen white dwarf (WD) i val in brightness with their host galaxies during near the Chandrasekhar mass that accretes mat- X several weeks. The importance of SNIa in as- ter from a companion in a close binary system. r trophysics and cosmology is highlighted by their a This model accounts for the SNIa sample homo- use as standard (or, better, calibrable) candles to geneity, the lack of prominent hydrogen lines in measurecosmicdistancesandrelatedcosmological their spectra, and its detection in elliptical galax- parameters(Perlmutter et al.1998;Schmidt et al. ies. Moreover, massive WDs are extremely un- 1998; Riess et al. 2001). However, in order to stable bodies in which a modest release of en- achieve a high precision in the distance determi- ergy can produce a huge change of size (i.e. the ratio of binding to gravitational energy is only 1Dept. F´ısica i Enginyeria Nuclear, Univ. Politècnica ∼15%,whileinanormalstaritis∼50%). There de Catalunya, Diagonal 647, 08028 Barcelona, Spain; are two main ingredients of the standard model [email protected] [email protected] [email protected] that are still poorly known: the precise configura- 2Institut d’Estudis Espacials de Catalunya, Barcelona, tion of the stellar binary system and its evolution Spain priortothermalrunawayoftheWD(Langer et al. 3Depto. F´ısica Teo´rica y del Cosmos, Univ. Granada; 2000; Nomoto et al. 2003; Piersanti et al. 2003; [email protected] 1 Han & Podsiadlowski 2004; Badenes et al. 2007; SNIa are remarkably homogeneous. In particular, Hachisu et al.2008),andtheexplosionmechanism the maximum velocity reached by intermediate (Hillebrandt & Niemeyer 2000). Fortunately, as mass elements, ∼ 11000 kms−1, indicates a sim- long as a carbon-oxygen WD reaches the Chan- ilar extent of thermonuclear burning in all SNIa. drasekhar mass, its previous evolution is not crit- Admitting thatpartofthe diversityshownby the ical for the explosion because the WD structure SNIa sample is directly related to the explosion is determined by the state of a degenerate gas of mechanism of Chandrasekhar-mass white dwarfs fermions that can be described with a single pa- (Hatano et al. 2000), multidimensional hydrody- rameter: the mass of the star. This fact leaves namic simulations of the phenomenon should be the explosionmechanismas the mostrelevantun- able to explain the wide range of nickel masses known concerning SNIa. that are seemingly synthesized in different events. In spite of continued theoretical efforts dedi- Lackingmultidimensionalcalculationsofthe light cated to understand the mechanism behind SNIa, curve and spectra of these models, it is certainly realistic simulations are still unable to provide riskyto comparedirectly the outcome ofthree di- a satisfactory description of the details of these mensional explosion simulations with SNIa obser- thermonuclear explosions. Nowadays, there is vations, in order to elucidate the more promising consensus that the initial phases of the explo- explosion mechanism and discriminate between sion involve a subsonic thermonuclear flame (de- models. However,even acknowledgingthat a per- flagration), whose propagationcompetes with the fect match between the observed luminosity and expansion of the WD. After a while, the cor- spectralevolutionofSNIaandthosededucedfrom rugation of the flame front induced by hydro- angle-averaged versions of the three dimensional dynamic instabilities culminates in an accelera- explosionmodelsisnottobeexpected,therearea tion of the effective combustion rate. However, few elementary requirements that a model of typ- the huge difference of lengthscales between the ical SNIa should pass: white dwarf (∼ 108 cm) and the flame width • There should be no more than a few tenths (.1 cm) prevents that large-scale numerical sim- of unburnt carbon moving at low veloci- ulations can resolve the deflagration front, mak- ties (Marion et al. 2003; Baron et al. 2003; ing it necessary to implement a model of the sub- Kozma et al. 2005). sonic flame. The precise value of the effective combustion rate is currently under debate, as it • The upper limit of the ejected mass of 56Ni depends on details of the flame model. Recent should reach ∼ 1 M⊙, after neutronization three-dimensional (3D) models calculated by dif- due to electron captures is properly taken ferent groups have shown that pure deflagrations into account (Mazzali et al. 2007). always give final kinetic energies that fall short of 1051 ergs, while leaving too much unburnt car- • Intermediate-mass elements should be gen- bonandoxygenclosetothecenter(Reinecke et al. erated in amounts exceeding ∼ 0.2 M⊙ 2002; Gamezo et al. 2003; Garc´ıa-Senz & Bravo (Garc´ıa-Senz et al. 2007). 2005; Ro¨pke et al. 2007). Because both signa- • The ejecta should be chemically stratified, tures are at odds with observational constraints, as demanded by observations of supernovae Gamezo et al. (2003, 2004, 2005) concluded that and remnants (e.g. Badenes et al. 2006; the only way to reconcile three-dimensional simu- Gerardy et al. 2007; Mazzali et al. 2008). lations with observations is to assume that a detonation ignites after a few tenths of a solar mass • Large clumps of radioactive 56Ni and other have been incinerated subsonically (delayed deto- Fe-group elements should not be present at nation). the photosphere at the time of maximum Through a systematic analysis of a large sam- brightness (Thomas et al. 2002). ple of well observed SNIa, Mazzali et al. (2007) showedthatinspiteofthesignificantdispersionof • In order to reproduce the observed width- supernova luminosities, implying different ejected luminosity relationship the kinetic energy masses of 56Ni, other characteristic features of and the total mass burned have to be in 2 the ranges K ∼ (1.0−1.4) × 1051 ergs Hewouldimprovethematchofcalculatedandob- and Mburnt ∼ 1.0 − 1.2 M⊙, respectively served light curves. There have been some claims (Woosley et al. 2007)1. ofdetectionofHeIlinesat2.04µmand1.052µm in the spectra of the spectroscopically normal SN Baron et al. (2008) calculated detailed non- 1994D and other SNIa (Cumming et al. 1994b,a; LTE synthetic spectra of angle-averaged spher- Nomoto et al.2003;Pignata et al.2008),although ically symmetric versions of Pulsating Reverse its identification is very uncertain (Wheeler et al. Detonation (PRD) models. As compared with 1998; Mazzali & Lucy 1998; Hatano et al. 1999). early, maximum, and post-maximum spectra of Inafewtestcalculationsshowninthepresentpa- wellobservedSNIa,the PRDmodelsproducedan per we have assumed the presence of such a small excessively prominent C II feature, and the spec- cap of He. tra appeared too red due to the presence of iron- In the present work we give details of the det- peak elements in the outer layers. Baron et al. onation phase of PRD models computed in three (2008)pointedoutthatfittingofrealSNIaspectra dimensions, together with an analysis of the nu- mightimproveifsomefractionoftheiron-peakel- cleosynthesis and the bolometric light curve ob- ements moving athigh velocities weresubstituted tained from angle-averagedspherically symmetric by intermediate-mass elements. versionsofsuchmodels. Wealsodiscusstheconse- In a companion paper (Bravo & Garc´ıa-Senz, quences that departures from spherical symmetry 2009; hereafter paper I) we have analyzed the of the ejecta can have over the spectra. In Sec- conditions for detonation ignition after a pulsat- tion 2 we explain the methodology used to com- ing phase of a white dwarf. A confined detona- pute the three-dimensionalsupernovamodels pre- tion ignition (CDI) can be induced by an accre- sented in this paper. In Section 3 the evolution tion shock, formed around a carbon-oxygen rich of the white dwarf during the detonation phase is core after pulsation following subsonic burning of described, while in Section 4 the properties of the a mass Mdefl < 0.3 M⊙. The mechanism of CDI homologous expanding ejecta are analyzed. Fi- is robust, and detonation ignition conditions are nally, our conclusions are given in Section 5. reached non-marginally within the shocked flow for a wide range of initial configurations. How- 2. Methods and models ever,if the nuclearenergyreleasedbefore the pulsation is close to the white dwarf binding energy, The simulations were performed with a three- Mdefl ∼ 0.29 M⊙, the accretion shock is not effi- dimensional Lagrangian Smoothed Particle Hy- cient enough to initiate a detonation in a mixture drodynamics (SPH) code suited to the modelling of carbon and oxygen. In that paper we have also of Type Ia supernovae (Garcia-Senz et al. 1998; shown that if C-O matter is contaminated by a Garc´ıa-Senz et al. 1999). Contrary to Eulerian small amount of He (X(He) ∼ 0.05− 0.10) the methods in which usually the spatial resolution conditions for a successful CDI are relaxed sub- is a free numerical parameter and the mass reso- stantially, making even easier the detonation of lution is determined by the local density, in SPH the core. A small cap made of He might rest atop it is the mass resolution of the calculation which of the white dwarfat the moment of thermal run- is fixed2. In SPH, the number of interpolating away as a result of previous accretion from the neighbors is what determines the spatial resolu- companion star in the binary system. This he- tion of the calculation. In our models, the num- liumcouldeitherhavebeenaccreteddirectlyfrom ber of neighbors was variable in the range 40-90, a He star or result from the nuclear processing whilethetotalnumberofequal-massparticleswas of accreted hydrogen. Ho¨flich & Khokhlov (1996) N = 250,000. The maximum spatial resolution suggestedthatthepresenceofa∼0.01M⊙ capof reached during the simulations of the detonation 1Woosleyetal.(2007)alsorequiredthattheinnerlayersof 2Although it is possible to dynamically increase the mass ejecta (. 0.8 M⊙) should be fairly well mixed, but such resolution during the SPH calculation using a technique extentofmixingcanbeinconflictwithwhatisrequiredto knownas particle splitting,that istheequivalent to AMR explainSNIaspectrabothintheopticalandintheinfrared methods inEuleriancodes,wehavenotuseditinthecal- (Motoharaetal. 2006; Gerardyetal. 2007; Mazzalietal. culationsreportedinthepresentpaper. 2008) 3 phaseofthePRDmodelswas9−12km,depending proach was used to make the calculation feasible. onthe model(the radiusofthe staratthatepoch As soon as the characteristic nuclear time-step of was∼40,000−95,000km). Furtherdetailsofthe a given particle became smaller than a prescribed numericaltechniques appliedin the SPHcode are fraction of the dynamical time-step its chemical provided in the appendix. evolution was followed isochorically,and both the The physical ingredients included in the hy- nuclear network and the equation of state were drocode allow to describe accurately all the phe- decoupled from the hydrodynamic evolution until nomena relevant to thermonuclear supernovae. thedynamicaltime-stepwasrecovered. Whenthe Gravity was calculated by means of the hierar- temperatureofamasselementreached5.5×109K, chical tree method (Hernquist 1987), retaining NSEwasassumed. Onceattained,NSEwasmain- up to the quadrupole terms in the multipolar ex- tained as long as the temperature stayed above pansion. In order to improve the accuracy of 2×109 K, providing a nuclear energy generation the calculation, the gravitational interaction be- rateaccurateenoughforthe hydrodynamicalsim- tween close particles was computed directly. The ulations. Weak interactions during NSE deter- equation of state includes contributions of par- mine the electron mole number, Ye, whose evo- tially degenerate and relativistic electrons with lution was calculated at each timestep by solving pair corrections, ions considered as an ideal gas the equation: dYe/dt = ΣiλiYi, where λi stands plus Coulomb corrections and radiation. The nu- forallkindofweakinteractions: λi >0forβ−dis- clear binding energy and electron capture rates integrations and e+ captures, and λi < 0 for β+ (Fuller et al. 1982; Mart´ınez-Pinedo et al. 2000) disintegrations and e− captures, while the molar onmatterinnuclearstatisticalequilibrium(NSE) fractions, Yi, were set by the NSE equations. were interpolated from a table that uses density, 2.1. Testing the code temperature and electron mole number as input. Weak interactions in general, and electron cap- We have tested the performance of our SPH tures in particular, are necessary for an accurate hydrocode simulating two three-dimensional de- description of thermonuclear supernovae because layed detonation scenarios available in the scien- they affect the dynamical evolution by decreas- tific literature, which had been computed using ing the total number of electrons that support quitedifferentnumericalmethods(bothforresolv- the white dwarf structure. The mean electron ingthe hydrodynamicalequationsandtodescribe mole number of the ashes, ∆Ya, determined by e the nuclear flames). We have chosen delayed det- electron captures during the deflagration phase onation models that bear many similarities with is ∆Ya ∼ 0.484 mol g−1. This neutronization e the PRD scenario, i.e. those in which a delayed translates into a decrease of electron pressure of detonation is formed at a low density and propa- asmuchas∆pa/p∼(8/3)∆Ya ∼4%,ofthe same e e gatesinwardsthrougha carbon-oxygenrichwhite order as the increase in pressure due to incinera- dwarfcore. OurreferencemodelsaretheGravita- tion of carbon and oxygen at a density of 109 g tionally Confined Detonation (GCD) model Y100 cm−3, ∆p /p∼10%. burn fromPlewa(2007),andtheturbulentdelayeddet- The nuclear network consumed a huge part of onation model DD 005 from (Ro¨pke & Niemeyer thecomputationalresourceschieflyduetothevery 2007). The first one begins with the incinera- small time-steps required to integrate accurately tion of a small off-centered bubble which floats to thenuclearkineticequations,andbecauseatsome low-density regions consuming a few hundredths point of the calculation there were tens of thou- ofasolarmasstofinallybreak-upthewhitedwarf sands of particles undergoing simultaneously nu- surface. The combustion ashes slide through the clear combustion. A small network of 9 isotopes: surface to finally concur at the opposite pole ig- n, 1H, 4He, 12C, 16O, 20Ne, 24Mg, 28Si, and 56Ni niting a detonation which processes most of the (Woosley 1986) allows to calculate the approxi- white dwarf to NSE, mainly in the form of 56Ni. matenuclearenergyinputandaroughnucleosyn- The turbulent delayed detonation model DD 005 thesis with a moderate load of CPU time. The is based on a physical criteria of detonation igni- performanceofthisnetworkwasextensivelytested tion,namelythetransitionfromalaminarflameto in Timmes et al. (2000). An operator-split ap- the distributed regime of nuclear combustion. We 4 have computed models GCD1 and TURB7 using of the neutronization on the explosion dynamics similarinitialmodelsandphysicalassumptionsas should be also larger in TURB7. The final mean forY1003andDD 005,respectively. ModelGCD1 electron mole number in GCD1 and TURB7 was started from a single incinerated bubble of radius Y =0.4964andY =0.4932,respectively,i.e. the e e 53 km located 100 km off-center. Model TURB7 neutronization in the turbulent model was about startedfromthesameinitialmodelasmodelDF29 double than in the other one. in paper I, but it incorporated a prescription for turbulent accelerationof the burning rate accord- 2.2. Pulsating reverse detonation simula- ingtothelocalvelocityfluctuations. Althoughthe tions initial configuration of TURB7 is not the same as Our working hypotheses are that during the inDD 005,bothmodelsburntsimilarmassessub- deflagration phase it is burnt an insufficient sonically,that were distributed in a smallnumber amount of mass to unbind the white dwarf, of ash clumps, and used the same criteria for the and that the outcome of the explosion is pre- deflagration-to-detonationtransition. As we want dominantly determined by M . The simula- defl to test the ability of the code to deal with det- tions of the deflagration phase started from ini- onation propagation, the similarity of the global tial models consisting on spherically symmetric propertiesofthe whitedwarfattheendofthede- cold (T = 108 K) isothermal white dwarfs built flagration era is what makes meaningful the com- through a fourth order Runge-Kutta integration. parison of models DD 005 and TURB7. The one-dimensional profile was mapped to a The outcome of the additional models is sum- three-dimensional distribution of N particles and marized in Table 1, where we give the amount of afterwards relaxed with the procedure explained mass burned during the initial deflagration phase inGarcia-Senz et al.(1998). Onceastableenough (when available), the final kinetic energy, and the configuration was obtained, the temperature of a ejected mass of 56Ni, NSE elements, carbon plus small subset of particles was artificially increased oxygen, and intermediate-mass elements (details up to the self-consistent value of NSE tempera- of these simulationsand the numericalsetup used ture atsuchdensity (assumingisochoricadiabatic will be presented in a forthcoming paper, Garc´ıa- burning)andtheensuingdynamicalevolutionwas Senz et. al. 2009). The mass of 56Ni in the ref- followed by means of the SPH hydrocode. erence models is just an upper limit because elec- The main properties of our PRD models are tron captures were not included in the hydrody- givenin Table 2. For clarity,we havechangedthe namic calculations. As can be seen in the Table, name of the models with respect to previous pub- the matching between our models and the refer- lications. In the present work, PRD models are ence ones is excellent, especially when one com- designated with the acronym PRDnn, where ’nn’ pares GCD1 with Y100, which are nearly twin gives the hundredths of a solar mass burned dur- models. Model TURB7 reproduces quite well the ing the initial deflagrationphase. For instance, in mainpropertiesofmodelDD 005,but somesmall modelPRD18themassburnedsubsonicallybefore differences remain in the synthesized masses of the pulsation was Mdefl = 0.18 M⊙ (this model intermediate-mass elements (IME) and NSE ele- was called PRD6 in Bravo & Garc´ıa-Senz 2006 ments. The reason for these discrepancies may and Baron et al. 2008, while the present model reside in different approaches to the description PRD14 was called PRD5.5 in Baron et al. 2008). of the flame as well as on the effect of the elec- The evolution of the central density and the tron captures on the dynamics of the explosion. chemical composition of model PRD14 is shown Even though neither DD 005 nor Y100 included in Fig. 1. At the beginning combustion propa- electron captures, the amount of mass burned at gated subsonically as a deflagration wave, with high densities during the deflagration phase was an effective combustion rate low enough to en- much larger in the turbulent delayed detonation sure that the flame quenched before the white model(0.38M⊙ vs ∼0.06M⊙), hence the impact dwarf became unbind. During the first seconds of the ensuing pulsation the incinerated bubbles 3Note that Y100 is a two dimensional model, while GCD1 floated to the surface, leading to composition in- wascomputedin3D version, i.e. the internal volume, which we call 5 hereafterthecore,wasplentyofcoldfuel,i.e. car- two times: just after the nuclear timescale be- bonandoxygen,whiletheashesoftheinitialcom- comesshorterthanthehydrodynamicaltimescale, bustion,mostlyhotironandnickel,werescattered and 0.1 s later (top images). The maximum tem- around (see Fig. 2). Shortly after, the expansion peratures achieved in the region are high enough of the core came to anend while partof the enve- (T & 5 × 109 K) for explosive silicon burning, lope re-collapsed. Eventually, an accretion shock whereas carbon and oxygen are destroyed in a was born because of the impact of the in-falling shortertimescaleoncethetemperaturerisesabove material onto the carbon-oxygen core. Upon re- ∼ 3×109 K as highlighted by the fuel contours. collapse nuclear reactions light againa few tenths The thermal gradient in the radial direction is ofa secondbefore the time ofmaximumcompres- large,ontheorderof100−200Kcm−1,althoughit sion (∼ 4.4 s). In this paper, our interest is fo- isinfluencedbytheresolutionofthenumericalcal- cused on the events that ensue after a detona- culation. Shock fronts propagate ahead and close tion is ignited due to the energy transferred to to the thermal waves, as revealed by the sculpt- the core through the accretion shock (the defla- ingoftheisobarsfollowingthecontourofthehigh gration phase has been thoroughly described in temperatureregionsandbythelargepressuregra- Garc´ıa-Senz & Bravo 2005). dient in front of the thermal front (bottom left image): the detonation is born. 3. Description of the detonation phase of The performance of the detonation capture al- PRD models gorithm is illustrated in the bottom right image on Fig. 3, where it is shown the region where the In paper I we have shown that after a pulsa- adaptiveellipticalinterpolatorgivesasharperdet- tion the white dwarf reaches conditions favorable onation front. As can be seen, a detonation is for a CDI if Mdefl is well below ∼ 0.3 M⊙ and detected all around the thermal wave except the nuclear reactions are turned off. However, when outermostpartcharacterizedbyashallowthermal nuclear reactions are included in the calculation gradient. Figure4displaysacomplementaryview the release of nuclear energy contributes not only showing the location (in fact, the projection onto totemperatureincreasebutalsotoraisingthegas the yz plane) of the actual SPH particles within pressure, thus disturbing the dynamical evolution the detonation front, i.e. those for which the det- during the collapsing phase of the pulsation. In onation capture condition has been verified. The the present simulations we have allowed the det- color map in the same Fig. shows the smoothing onation to ignite spontaneously at the time and length, whose value at the detonation front is of place determined by the SPH code as a result order 30 km. of the hydrodynamical evolution coupled to the nuclear network. This approach must be looked 3.2. Burning of the white dwarf core with caution as we cannot resolve the detonation front thickness, hence our simulations might not In the PRD scenario, the last phase of the ex- describe accurately the process of detonation for- plosionstartswhentheconvergingreversedetona- mation. However, as we will see later in Sect. 4, tionwaveislaunched. Onceinitiated,the detona- the explosion properties are quite independent of tionis self-sustainedby the nuclearenergyrelease the details of the detonation initiation, hence our from burning of the carbon-oxygen fuel that per- approach is justified. vades the hydrostatic core. The evolution of the radii of spherical shells of 3.1. Initial steps of detonation propaga- givenLagrangianmassisshowninFig.5,together tion with the evolution of the detonation front. Be- cause the model does not have spherical symme- It is interesting to analyze the initial steps of try, the location of the detonation front is only detonation propagation, keeping in mind the cau- approximate. In that figure, the inner and outer tionary remark given above about the capabili- boundaries of the detonated region have been de- ties of our code to describe accurately the birth fined as the radius of the shells for which at least of a detonation front. Figure 3 shows the evolu- half ofthe carbonand oxygenremainingafter the tion of a zoomed slice across model PRD14 for deflagrationphasehadbeenprocessedby the det- 6 onation wave (see Fig. 6 and the next paragraph unburntmass,∼0.08M⊙,andthecentraldensity, for a complementary view of the evolution of the ρ ∼4×108gcm−3,giveavolumeof∼4×108km3 c detonationwave). It canbe seenthatthe inwards or a typical lengthscale of the unburned spot of propagation of the thermonuclear burning is very ∼ 450 km. Our resolution at the center at that fast, as corresponds to a detonation front. This time was of order 11 km, implying that the det- detonation front processes most of the inner core onation failure at the center of the white dwarf of the white dwarf, but it stalls before arriving was not due to insufficient numerical resolution. to the very central regions and fails to burn the The physical reason for detonation failure at the innermost few hundredths of a solar mass. center is the expansion of the external, already Fig. 6 illustrates the evolution of the white detonated, core matter. As the detonation wave dwarf core during the propagation of the reverse is accessible to the shockedmatter throughsound detonation for model PRD14, with a view com- waves, the rarefaction wave originated by expan- plementary to the one in Fig. 5. At the time of sion of the core caught and weakened the shock the firstsnapshot,t=3.8s,the accretionshockis andfinally quenchedthe burning closeto the cen- about to form. The C+O mass fraction contours ter. It took a few tenths ofa secondsince detona- maintain a high degree of spherical symmetry in tion initiationto releaseenoughnuclear energyto the core, increasing their value steeply towards overcome the impact pressure of accreting matter the center (lighter contours belong to largerC+O and allow the detonated matter to expand appre- mass fractions). In the first snapshot, a hot spot ciably. As a result of such delay most of the core is visible in the yz plane at r ∼ 2000 km, which matterwasdetonated,leavingonlyasmallmassof is a result of a slight asymmetry in the collapsing unburnedfuelnearthe center. Theamountofun- flow (the temperature and chemical composition burned carbonandoxygenis sensitive to the den- maps a few tenths of a second before do not show sitystructureatthemomentofdetonationforma- any sign indicating the imminent formation of a tion, beingthe largestforthe lessenergeticmodel hot spot). The combination of high temperatures PRD26 (see Table 2). and abundant fuel culminated in the formationof 3.3. Sensitivity and resolution a detonation already visible in the second snapshot. At the same time, the accretion shock had Wehavecheckedthesensitivityoftheresultsof heated the fuel in a hot spherical shell located at our simulations to the numerical method, in par- r ∼ 3000 km up to T & 1.5×109 K. The lateral ticulartotheimplementationofartificialviscosity, spreading of burning was not due to detonation and to the spatial resolution. As it is well known, propagation, but it was caused by the accretion artificialviscosityintroducesnumericaldissipation shock, that continuously fed mechanical energy in the evolution equations, which can generate into the core sparking detonations all around the enough entropy at shocks to induce an artificial hot shell. During the first 0.3 s after detonation detonation. We have computed two additional ignition the accretion shock remained stationary, models starting from Mdefl = 0.18 M⊙, one with close to the core. Afterwards, the overpressure the artificial viscosity parameters enlarged by a generated by the nuclear energy released pushed factor two, the other with the parameters divided out the accretion shock, that detached from the by the same factor. In both cases, the code in- dense core (last snapshot). In the outer layers of cluded an algorithm to inhibit combustion inside thecorethedensitywaslowenoughtoallowforin- the artificially widened shock fronts. The detona- complete silicon burning and leave a composition tion developed in both cases in the same way as rich in intermediate mass elements, mainly Si, S, in our standard calculation, irrespectively of the Ca and Ar, while in the inner regions the burning valuestakenbytheartificialviscosityparameters. proceededallthe wayupto56Ni(the centralden- The finalkinetic energiesand56Nimassessynthe- sityatthe momentofformationofthe detonation sizeddifferredbylessthan8%withrespecttoour was in the range (0.5−4)×108 g cm−3). standard model PRD18. In our simulations the detonation propagated Wealsocheckedthe sensitivityofthe resultsto inwardsburningallthefuelexceptforasmallvol- thenumericalresolutionbyrunningmodelPRD18 ume close to the center of the white dwarf. The with the same number of particles but reducing 7 the number ofneighbors,n . The smallernumber wouldprobablygiverisetoexplosionswithslightly 0 of neighbors translated into a smaller smoothing smallerkineticenergies,substantiallysmaller56Ni length, and thus in an increased capability of the masses and much larger intermediate-mass ele- codetorepresentsmallscalefeatures. Suchproce- ments yields. Unfortunately, due to the failure dure allowedus to reacha maximumresolutionof to maintain a steady detonation in model PRD29 4 km during detonation propagation. Once more, we have not been able to verify this last point. thedetonationdevelopedinthesamewayasinour As seen in the last snapshot of Fig. 6, the det- standard calculation. As a drawback, the small onation smooths the thermal and chemical struc- number of neighbors originated numerical insta- ture of the white dwarf, wiping out the imprints bilitiesintheoutermostlow-densityregionsofthe left by the deflagration phase in the central vol- ejecta. These instabilities prevented the model to ume of the ejecta. The detonation also endows reach the homologous expansion phase, making it themechanicalstructureoftheejectawithalarge impossibletocomparethefinaloutcomewiththat degree of spherical symmetry. Fig. 7 shows the of our standard PRD18 model. homologous density profile of model PRD14 at a Finally, we checked if the outcome of the ex- timeof34ssincetheinitialthermalrunaway. The plosion was sensitive to the presence of a small solidlinerepresentstheangle-averageddensityob- mass of He atop of the white dwarf. This helium tained in 1000 shells evenly distributed in mass, mixed with the underlying layers of carbon and whiletheerrorbarsgivethestandarddeviationof oxygen during the pulsation and was present at the particles density. The small dispersion is an the detonated layer at ignition time. By adding indication of the high degree of spherical symme- either a total He mass of 0.005 M⊙ (i.e. a 0.36%) try. The high-velocity low-density external mate- or 0.010 M⊙ (i.e. a 0.72%) to model PRD18 the rial is best approximated by an exponential den- kinetic energy and the mass of 56Ni synthesized sity law, whereas the central region can be ap- in the explosion changed negligibly (. 4%) with proximatelydescribedasanearlyuniformdensity respect to the standard calculation of the same core whose mass is ∼ 0.82 M⊙. The exponential model. Thus, we conclude that the eventual pres- tail begins at a velocity ∼ 10000 km s−1, where ence of He does not introduce additional uncer- the intermediate-masselements(Si, S,Ca,...) are tainties in the outcome of the PRD models. abundant: XIME & 0.2. The innermost 0.08 M⊙, composedmainlybyunburnedC-O,isdenserthan 4. Nucleosynthesis, light curves and spec- the uniform core by up to an order of magnitude. tra: asymmetries, inhomogeneities A comparisonof the density profiles of the last computed models (normalized by t3) is presented In Table 2 there are given the properties of in Fig. 8. The larger kinetic energy of model the PRD models we have computed, covering PRD14 with respect to PRD18 produces a more the range of subsonically burned masses M = defl extended tail in the former. The beginning of the 0.14 − 0.26 M⊙. The kinetic energy at infinity exponential tail matches quite approximately the spans the range (1.0−1.2) × 1051 erg, in good location where the detonation was ignited in each agreement with requirements deduced from the model(seevaluesofv inTable2). Whereasthe SNIa light curve shape (see Sect. 1). In principle, deto transitionfromthe centraluniformdensity region it should be possible to obtain PRD models with to the exponential tail is quite smooth in model stillsmallervaluesofMdefl,downto∼0.10M⊙ or PRD14, it becomes increasingly abrupt as M even less, which according to the trend displayed deto and v decrease. In model PRD26 there is not in Table 2 would translate into larger kinetic en- deto a uniform density core, which is substituted by a ergies and 56Ni masses. However, as the accre- density inversion below M . In models PRD18 tion shock in model PRD14 already forms quite deto and PRD26 the density of the innermost region close to the total mass of the white dwarf, it is rich in fuel is ∼4 times larger than the density of to be expected that the outcome of the explosion the nearlyhomogeneouscore,afactormuchlower for smaller M will not be much different from defl than that found in model PRD14. thatofmodelPRD14. Ontheotherside,valuesof M slightly larger than those in model PRD26 defl 8 4.1. Nucleosynthesis belongingtothelastgroupwascomputedbyinte- gration of the nuclear kinetic equations along the Figs. 9 and 10 depict three-dimensional views thermodynamic trajectory recorded for that par- ofthe chemicalstructureofthe ejecta. InFig.9a ticle during the SPH simulations. The method renderingofthefinaldistributionof56Niinmodel followed for the integration of the nuclear equa- PRD14 is shown. The central volume is homoge- tions was described in Cabezón et al. (2004). Fi- neously filled with 56Ni, except for a small ball nally, the particles belonging to the NSE group made of carbon and oxygen that is not visible in were again divided in two categories according to this image, surrounded by intermediate-mass ele- the density at which they left NSE. Those par- ments and unburned carbon and oxygen. Small ticles that cooled below 2 × 109 K at densities and medium size clumps of 56Ni and other Fe- in excess of 108 g cm−3were assigned the chem- group elements pervade the outer layers of the icalcompositioninNSEatthesetemperatureand ejecta. Theseclumpshaveadoubleorigin: someof density, and at the electron mole number given themwereproducedduringtheinitialdeflagration by the SPH model. Conversely, particles that left phase and migrated to the external layers of the NSE at densities below 108 g cm−3were assumed whitedwarfduringtheexpandingphaseofthepul- to experience an alpha-rich freeze-out, and their sation,whileothersdetachedfromtheincinerated final chemical composition was obtained by inte- core after the end of the detonation phase. These grating the nuclear network following the cooling clumps are the end product of the impact of the path of the particle, starting from T = T and detonated core into the surrounding low-velocity NSE the chemical composition of NSE matter at that expanding layers. Thus the final distribution of temperature. 56Ni comes from two different regimes: a smooth inner region, rich in 56Ni and other Fe-group ele- The mass of 56Ni synthesized in each PRD model can be found in Table 2, while Table 3 ments, and a clumpy outer region where radioac- gives the final abundances after radioactive de- tiveandstableFe-groupandintermediate-massel- cays,wherewehavesummeduptheabundancesof ementssharethevolumewithcarbonandoxygen. theisotopesofeachelementandweshowonlythe A complementary picture of the chemical struc- most abundant species in the ejecta. The amount ture is given in Fig. 10, in which the mean molar weight in an octant of model PRD14 is displayed. of 56Ni produced is in the range 0.6 − 0.8 M⊙, whichis enoughtoexplainthe lightcurvesoftyp- It shows the same gross features seen in Fig. 9, icalSNIa. Aspointedbefore,modelswithM in but here more details are visible, as for instance defl the centralball made of carbonand oxygen. Nev- excessof∼0.27M⊙ mightproduceloweramounts of56Nithatmightaccountforsub-luminousSNIa, ertheless, carbon and oxygen are predominantly butwehavenotbeenabletoverifyitinthepresent located at the outermost layers of the ejecta. work. The total mass of IME synthesized in the In order to compute even coarse light curves three computed models spans a similarly narrow andspectra ofSNIa models andto performa pre- range, M(IME)∼0.20−0.26 M⊙. liminary comparison with observations, it is nec- Figures 11-13 show the final (after radioactive essary to obtain the isotopic yields up to A∼60. decays)chemicalprofilesofthethreePRDmodels AsourSPHcodeimplementsjusta9-isotopesnu- as a function of the ejecta velocity. The location clearnetwork,wehavepost-processedthethermo- of the accretion shock, v , and the CDI, v , dynamic history of the particles to derive the nu- acc deto left no imprint on the chemical profile. The dis- cleosynthetic output of eachexplosionmodel. For persionofthemassfractionsagreeswiththequal- eachcomputedmodel,theSPHparticleswereclas- itative view obtained with Figs. 9 and 10: the in- sified into three groups: 1) those whose temperatures neverexceeded109 K,thatwere assumedto ner volumeis chemicallystratifiedwhile the outer layers display large scatters. The transition be- keeptheirinitialcompositionunchanged,2)those tween both regimes takes place at an ejecta ve- that reached temperatures in excess of T = 5.5 × 109 K, that were assumed to attainNSNESE locity of v ≃ 10000 km s−1 in model PRD14 and v ≃ 8000 km s−1 in model PRD18, while composition and, 3) the remaining particles, that in model PRD26 the scatter is already high at a were candidates to experience incomplete explo- few1000kms−1. Theabundanceofintermediate- sive burning. The nucleosynthesisof each particle 9 masselementsincreasestowardsthesurfaceofthe ric light curves from spherically symmetric angle- ejecta in all models. In model PRD14 the IME averagedversionsofthree-dimensionalSNIamod- amount to a significant mass fraction at velocities els. Bolometric light curves should be less sensi- in excess of v ≃10000 km s−1, and are dominant tive to the averaging procedure than broadband abovev ≃20000kms−1. The largemassfraction light curves, as the first ones are not affected ofFe-groupelementsathighvelocitiesisaconcern by shifting of photon wavelengths in and out of because its expected spectral signatures are not monochromatic bands, as due for instance to iron seeninSNIa. The sameistrueformodels PRD18 line blocking. and PRD26. In model PRD18 the abundance of Bolometric light curves have been computed carbonplus oxygenincreasesmonotonically up to from angle-averaged versions of the PRD models thesurfacewhereitamountstoX(C+O)∼0.30. by means of the one-dimensional code described TheC+OprofileinmodelPRD26showsapeakat in Bravo et al. (1993, 1996). In general, the bolo- thelocationoftheCDI,vdeto =8000kms−1,and metric light curves obtained with this code dur- declines at larger velocities. Although this model ing the pre-maximum and maximum phases are has the largest mass of IME, Fe-group elements in fairly good agreement with those computed by are still dominant through the outermost layers. directly solving the radiative transfer equations InFig.14itisshownthedetailedchemicalcom- (seeHo¨flich et al.1993). AlthoughourPRDmod- position of model PRD14 after radioactivedecays els reach the peak of luminosity slightly earlier as a function of the velocity. We stress that the than expected (∼ 12 days), the maximum lumi- poorchemicalstratificationoftheejectaismainly nosity, L , and 56Ni mass, M 56Ni , match al- max duetothemixingbetweenthedetonatedcorema- most perfectly the empirical relationship derived (cid:0) (cid:1) terialandthesurroundingmatterthatisproduced by Stritzinger & Leibundgut (2005, their Eq. 6). during the initial phases of expansion following In Fig. 16 we show a plot of M(56Ni) vs the bolo- the end of the detonation era. Such a mixing is metric light curve decline parameter ∆m , both 15 expected to happen in any model in which a det- fromthe presentthree-dimensionalmodels aswell onation propagates inwards from a point close to as from observational estimates for a set of well the surface of the star, as is the case in the two observed SNIa (Stritzinger et al. 2006). Models test models presented in Sect. 2.1. This is illus- PRD14 and PRD26 reproduce nicely the obser- trated in Fig. 15, that shows the chemical profiles vational trend, while model PRD18 declines too of model TURB7 at two times: just after the end fast in comparison to SNIa having a similar mass of detonation and when the ejecta is in the ho- of 56Ni. Models GCD1 and TURB7 synthesize a mologous expansion phase. It is plain to see that largeamountof56Niand,duetotheirlargekinetic the chemicalspecies areredistributedthroughthe energyandthepresenceofradioactivenucleiupto outer ∼ 0.6 M⊙ between both times, with the re- the surface of the ejecta, decline much faster than sult that Fe-groupelements replaceC-O andIME observations seem to indicate (the main parame- in the outer layers. Thus, the final chemical pro- ter determining the width of the bolometric light file of TURB7 displays also a large mass of Fe-Ni curve is the photon diffusion time, t , that is in d moving at high velocities. turn affected by the kinetic energy: t ∝ K−1/4, d Woosley et al. 2007). 4.2. Light curves 4.3. Spectra Multidimensional thermonuclear supernova models assuming different explosion mechanisms Baron et al. (2008) obtained detailed spectra havebeenshowntoreproducereasonablywellthe ofangle-averagedone-dimensionalversionsofsev- gross features of broadband light curves of nor- eral PRD models, from slightly before up to a mal Type Ia supernovae (Blinnikov et al. 2006; few days after maximum light. They concluded Kasen & Plewa 2007; Ro¨pke et al. 2007). The that the synthetic spectra of the PRD models model light curves in these works were obtained did not match the spectra of a few typical well- from one and two-dimensional averaged versions observedSNIa. The maindisagreementsfoundby of two and three-dimensional supernova models. Baron et al. (2008) were due to the presence of An alternative approach is to compute bolomet- Fe-groupelementsintheouterlayersoftheejecta, 10