Table Of ContentPulsating reverse detonation models of Type Ia supernovae. II:
Explosion
Eduardo Bravo1,2, Domingo Garc´ıa-Senz1,2, Rub´en M. Cabez´on1, 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`ecnica ∼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 det-
onation 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 pul-
sation 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 snap-
shot. 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´on 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 tempera-
tures 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
