MNRAS
496,
2039–2084 (2020)
doi:10.1093/mnras/staa1691
Advance Access publication 2020 June 15
Three-dimensional models of core-collapse supernovae from low-mass
progenitors with implications for Crab
G. Stockinger,
1
,
2
H.-T. Janka
,
1
‹
D. Kresse,
1
,
2
‹
T. Melson,
1
T. Ertl,
1
M. Gabler,
3
A. Gessner,
4
,
5
A. Wongwathanarat,
1
A. Tolstov,
6
S.-C. Leung,
7
K. Nomoto
8
and A. Heger
9
,
10
,
11
,
12
,
13
Affiliations are listed at the end of the paper
Accepted 2020 June 10. Received 2020 June 4; in original form 2020 May 6
ABSTRACT
We present 3D full-sphere supernova simulations of non-rotating low-mass (
∼
9M
)pro-
genitors, covering the entire evolution from core collapse through bounce and shock revival,
through shock breakout from the stellar surface, until fallback is completed several days later.
We obtain low-energy explosions (
∼
0.5–1.0
×
10
50
erg) of iron-core progenitors at the low-
mass end of the core-collapse supernova (LMCCSN) domain and compare to a super-AGB
(sAGB) progenitor with an oxygen–neon–magnesium core that collapses and explodes as
electron-capture supernova (ECSN). The onset of the explosion in the LMCCSN models is
modelled self-consistently using the
VERTEX
-
PROMETHEUS
code, whereas the ECSN explosion
is modelled using parametric neutrino transport in the
PROMETHEUS
-HOTB code, choosing
different explosion energies in the range of previous self-consistent models. The sAGB and
LMCCSN progenitors that share structural similarities have almost spherical explosions with
little metal mixing into the hydrogen envelope. A LMCCSN with less second dredge-up results
in a highly asymmetric explosion. It shows efficient mixing and dramatic shock deceleration in
the extended hydrogen envelope. Both properties allow fast nickel plumes to catch up with the
shock, leading to extreme shock deformation and aspherical shock breakout. Fallback masses
of
5
×
10
−
3
M
have no significant effects on the neutron star (NS) masses and kicks. The
anisotropic fallback carries considerable angular momentum, however, and determines the
spin of the newly born NS. The LMCCSN model with less second dredge-up results in a
hydrodynamic and neutrino-induced NS kick of
>
40 km s
−
1
and a NS spin period of
∼
30 ms,
both not largely different from those of the Crab pulsar at birth.
Key words:
hydrodynamics – neutrinos – stars: massive – stars: neutron – supernovae: gen-
eral – supernovae: individual: Crab.
1 INTRODUCTION
According to current understanding, stars with initial masses of
8–9 M
end their lives in a core-collapse supernova (CCSN). The
explosion is powered by gravitational energy, which is released
when the core of the star collapses to a compact remnant (a
neutron star or black hole), and a fraction of which is transferred to
the ejecta by neutrino-energy deposition (Colgate & White
1966
;
Bethe & Wilson
1985
). In the past six decades, numerous studies
have focused on the collapse phase and subsequent evolution using
1D simulations. Over the past three decades, multidimensional
simulations with successively improved treatment of the micro-
E-mail:
thj@mpa-garching.mpg.de
(HTJ);
danielkr@mpa-garching.mpg.de
(DK)
physics have driven our understanding of the explosion mechanism
(see e.g. Janka
2012
; Janka et al.
2012
; Burrows
2013
; Janka,
Melson & Summa
2016
;M
̈
uller
2016
). A close analysis of these
models led to the discovery of new hydrodynamic instabilities
such as the standing accretion shock instability (SASI) (Blondin,
Mezzacappa & DeMarino
2003
; Ohnishi, Kotake & Yamada
2006
;
Blondin & Mezzacappa
2007
; Foglizzo et al.
2007
;Fern
́
andez
2010
). The increase of computational capabilities in the recent years
along with new developments for neutrino transport methods (see
e.g. Buras et al.
2006
; Takiwaki, Kotake & Suwa
2012
; O’Connor &
Couch
2018a
; Glas et al.
2019a
; Skinner et al.
2019
and references
therein) have enabled full 3D simulations of the early explosion
phase (see e.g. Couch & O’Connor
2014
; Takiwaki, Kotake & Suwa
2014
; Lentz et al.
2015
; Melson, Janka & Marek
2015a
;Melson
et al.
2015b
; Roberts et al.
2016
; Summa et al.
2016
;M
̈
uller et al.
C
The Author(s) 2020.
Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Cre
ative
Commons Attribution License (
http://creativecommons.org/licenses/by/4.0/
), which permits unrestricted reuse, distribution, and reproduction in any medium,
provided the original work is properly cited.
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
2040
G. Stockinger et al.
2017
,
2018
; O’Connor & Couch
2018b
; Ott et al.
2018
;Vartanyan
et al.
2018
; Burrows, Radice & Vartanyan
2019
; Vartanyan et al.
2019
; Melson, Kresse & Janka
2020
).
Motivated by the historical SN1987A and its progenitor de-
tection, a variety of studies in two and three dimensions also
investigated the propagation of the shock wave from its initiation
to its breakout from the stellar surface in 15–20 M
blue supergiant
(BSG) models, which are suitable as progenitors of SN1987A.
In first studies, the blast wave was launched through artificial
energy deposition near the center (see e.g. Nomoto, Shigeyama &
Hashimoto
1987
; Nomoto et al.
1988
;Arnett,Fryxell&M
̈
uller
1989b
; Fryxell, M
̈
uller & Arnett
1991
;M
̈
uller, Fryxell & Arnett
1991b
), later the explosion was initiated with the neutrino-driven
mechanism (see Kifonidis et al.
2003
; Hammer, Janka & M
̈
uller
2010
; Wongwathanarat, Janka & M
̈
uller
2013
; Wongwathanarat,
M
̈
uller & Janka
2015
). M
̈
uller et al. (
2018
) followed the long-time
evolution of the explosion of ultra-stripped progenitors by 3D simu-
lations, motivated by the importance of such stars in understanding
the progenitor systems of the recent detections of NS–NS mergers
by LIGO/Virgo (GW170817, Abbott, LIGO Scientific Collabora-
tion & VIRGO Collaboration
2017
, and GW190425, Abbott, LIGO
Scientific Collaboration & VIRGO Collaboration
2020
).
These theoretical works showed that supernova (SN) explosions
are by far not spherical events as previously thought. Three-
dimensional instabilities facilitate the explosion (Herant et al.
1994
; Burrows, Hayes & Fryxell
1995
; Janka & M
̈
uller
1996
)and
are a necessary ingredient to explain the clumpiness and mixing
found in photospheric emission (Utrobin et al.
2015
,
2017
)and
spectral analyses of the nebular phase of core-collapse events
(Jerkstrand
2017
). Wongwathanarat et al. (
2015
) showed that the
final ejecta distribution carries imprints of the asphericities pro-
duced during the onset of the explosion (
t
∼
1 s), which are further
modified during later phases. Depending on the detailed progenitor
structure, hydrodynamic instabilities arising at the composition
interfaces, such as the Rayleigh–Taylor (RT) instability and the
Richtmeyer–Meshkov (RM) instability, shape the final spatial and
velocity distributions of nucleosynthetic products. The resulting
ejecta morphology ranges from quasi-spherical ejecta (M
̈
uller
et al.
2018
) to strongly pronounced RT-fingers including cases
that resemble the geometry found in Cas A (Grefenstette et al.
2017
; Wongwathanarat et al.
2017
). Due to the highly non-linear
and stochastic behaviour of non-radial instabilities and turbulence
during the onset of the explosion phase and the subsequent evolution
of RT instabilities, which depend on the progenitor structure, a
clear connection between the asymmetries, and thus the degree
of mixing, and the progenitor properties has still to be worked
out.
In this paper, we consider CCSN progenitors with initial masses
near the low-mass end of about 9–10 M
, where around 20 per cent
of all CCSNe are thought to occur (assuming a Salpeter initial mass
function Salpeter (
1955
) and an upper mass limit of
∼
20 M
for
CCSNe), to study the differences in their development of mixing
instabilities during the explosion. To this end, we compare ECSNe
from a super-AGB progenitor with an ONeMg core and CCSNe
from red supergiant (RSG) progenitors with iron cores, all in a
zero-age main sequence (ZAMS) mass range around 9 M
.
The evolution of stars with masses
12 M
is very sensitive
to the initial stellar mass, various pulsational instabilities, and
mass-loss phenomena (Woosley & Heger
2015
). Iron-core CCSN
progenitors with initial masses around 9–10 M
ignite oxygen
burning off-centre in contrast to their more massive (
M>
15 M
)
counterparts. After oxygen burning, silicon ignites in a degenerate
flash which might, in some cases, lead to additional mass-loss in
the last decade of evolution or is speculated to even eject parts of
the hydrogen envelope (Woosley & Heger
2015
).
The 8
.
8M
progenitor of Nomoto (
1984
) is even more peculiar. It
experiences several thermal pulses and off-centre ignition of fusion
material. In the end, it has a degenerate ONeMg-core surrounded
by a dilute and extended hydrogen envelope.
When the core approaches its Chandrasekhar mass, electron cap-
tures on
24
Mg and
20
Ne via the reaction chains
24
Mg(
e
−
,
ν
e
)
24
Na(
e
−
,
ν
e
)
24
Ne and
20
Ne(
e
−
,
ν
e
)
20
F(
e
−
,
ν
e
)
20
O (Miyaji et al.
1980
) destabi-
lize the core due to a reduction of the effective adiabatic index of the
electron-degeneracy dominated gas pressure. Continuous electron
capture on
20
Ne, which further reduces the pressure support, works
against the now beginning oxygen-burning as temperatures increase
during the collapse. Simulations in 1D (see Kitaura, Janka &
Hillebrandt
2006
; Fischer et al.
2010
;H
̈
udepohl et al.
2010
) and 2D
(see Janka et al.
2008
;Radiceetal.
2017
) suggest that the collapse
proceeds despite the oxygen burning. Jones et al. (
2016
) simulated
the deflagration of oxygen in ONeMg cores with different core den-
sities. At log
10
(
ρ
c
/
gcm
−
3
)
=
9
.
95 and lower densities their cores do
not collapse but get partly unbound due to the inefficient semicon-
vective mixing during the electron-capture phase and the resulting
strong thermonuclear runaway. Only when the central densities are
higher than this threshold value of
ρ
c
, the core is found to collapse to
a proto-neutron star (PNS). Recently, Kirsebom et al. (
2019
)inves-
tigated the influence of a newly measured strong transition between
the ground states of
20
Ne and
20
F on the electron-capture rate and
thus on the evolution of ONeMg cores. Adding the new transition
increases the likelihood that the star is (partially) disrupted by a
thermonuclear explosion (termed tECSN) rather than collapsing
to form a PNS. However, Zha et al. (
2019
), using state-of-the-art
electron-capture rates including the latest rate for the second forbid-
den transition of
20
Ne(
e
−
,
ν
e
)
20
F from Suzuki et al. (
2019
), found
that the oxygen deflagration starting from log
10
(
ρ
c
/
gcm
−
3
)
>
10
.
01
(
<
10
.
01) leads to collapse (thermonuclear explosion). Their es-
timate of the central density when the oxygen deflagration
is initiated in an evolving ONeMg core exceeds this critical
value. Therefore they conclude that ONeMg cores are likely to
collapse.
For this reason, in our study we assume that the ONeMg core
collapses to a PNS, leading to a ‘collapse ECSN’ (cECSN).
This assumption receives additional motivation by the fact that
recent studies considering the galactic chemical evolution of the
Milky Way stress the importance of cECSNe to reproduce the
solar abundances of several important and problematic isotopes
including, e.g.
48
Ca,
50
Ti, and several of the isotopes from Zn to Zr
(Jones et al.
2019
).
Despite the narrow range of central densities for which cECSNe
are expected to occur (Leung, Nomoto & Suzuki
2020
), and
despite the open questions associated with a variety of competing
processes that decide about collapse or thermonuclear explosion
and that depend strongly on many uncertain aspects of the employed
physics, connections to cECSNe have been made for observations of
SN1994N, 1997D, 1999br, 1999eu, 2011dc, and 2005cs (Stevenson
2014
). However, comparisons of the nebular spectra of some of
these cases with 1D neutrino-driven SN models are ambiguous or
disfavour the link to cECSNe (Jerkstrand et al.
2018
). Also SN1054
(the Crab) has been speculated to be a cECSN (Hillebrandt
1982
;
Nomoto et al.
1982
; Smith
2013
; Tominaga, Blinnikov & Nomoto
2013
), although such an interpretation is in conflict with results by
Gessner & Janka (
2018
) for the maximum kick velocity of PNSs
produced by cECSNe.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
3D models of supernovae from low-mass progenitors
2041
With the help of full-sphere three-dimensional simulations we
aim at investigating the following questions:
(i) What are the differences in the early stages (first seconds) of
the explosion in low-mass Fe-core and ONeMg-core progenitors?
(ii) What is the influence of the different progenitor structures on
the long-time evolution of the explosion? In particular, what is the
influence on the formation of reverse shocks and the efficiency of
outward mixing of neutrino-heated material?
(iii) Are CCSNe of low-mass progenitors able to produce highly
asymmetric ejecta and strong radial mixing of metals similar to
findings for more massive RSG and BSG stars in previous studies?
(iv) How do the properties of the compact remnants change on
long time-scales due to the fallback of matter? Are there significant
changes to the remnants’ mass, kick, and angular momentum?
The structure and contents of our paper are the following:
In Section 2, the basic properties of the considered models of
non-rotating, low-mass (super-AGB and RSG) progenitors are
introduced. Section 3 provides a brief description of the numerical
methods and input physics used in our simulations. Section 4
contains our results for the first second(s) of the explosion, focusing
on shock dynamics, explosion energies, neutrino emission, PNS
properties and the chemical composition of the ejecta. In Section 5,
for the first time in 3D explosion modelling, the SN evolution
of low-mass super-AGB and RSG progenitors is described until
and beyond shock breakout concerning the development of mixing
instabilities and ejecta asymmetries, the spatial distribution of
chemical elements, and the effects of fallback on the properties of
the newly formed NSs. In Section 6, we briefly compare our results
with previous studies, and in Section 7, we conclude with a summary
and discussion. Several appendices contain basic information on
more technical aspects concerning the simulation inputs and the
analysis methods.
2 PROGENITORS
In this paper, we study an ECSN of a non-rotating 8.8 M
super-
AGB star (e8.8), which is constructed from the envelope model
of Jones et al. (
2013
) and a collapsing core model (Leung et al.
2020
; Tolstov, Leung, and Nomoto, private communication), and
two CCSNe resulting from non-rotating low-mass RSGs with iron
cores (z9.6, s9.0), evolved to the onset of collapse by Heger (private
communication) and by Woosley & Heger (
2015
), respectively.
The considered ECSN progenitor is explored here for the first time,
whereas the explosions of the iron core progenitors were simulated
in 3D with the
VERTEX
-
PROMETHEUS
code and some results were
publishedinMelsonetal.(
2015a
,
2020
).
2.1 ONeMg-core progenitor
Model e8.8is a solar-metallicity progenitor with a ZAMS mass of
M
ZAMS
=
8
.
8M
. Its degenerate ONeMg-core undergoes electron
capture which ignites the O-Ne deflagration at the center. The
central density (
ρ
c
) and the electron fraction (
Y
e
), and thus the
core mass
M
(ONeMg) at the ignition, depend on the convective
stability criterion, extent of the convective mixing, and the electron
capture rate (Zha et al.
2019
). Leung et al. (
2020
) adopted the core
with
ρ
c
∼
10
9
.
975
gcm
−
3
,
Y
e
=
0
.
496, and
M
(ONeMg)
=
1
.
39 M
,
and calculated the propagation of the convective deflagration.
The deflagration incinerates ONeMg-core material into nuclear
statistical equilibrium (NSE). In the NSE region, electron capture
on iron-group nuclei and free protons dominates nuclear energy
release, thus inducing collapse. When the central density reaches
ρ
c
∼
10
10
.
64
gcm
−
3
, the NSE region extends to
∼
0.45 M
.TheSN
simulations started from this progenitor condition. The structure of
the core is well approximated by a spherical model.
When transferring the progenitor model of e8.8 to the SN
simulations, the mass of the ONeMg-core was reduced to 1.34 M
in
the course of providing an easy-to-handle fit to the complex density
structure. This has no relevant influence on the dynamical evolution
of the explosion, as can be concluded from the close similarity
of the explosion behaviour of model e8.8 with that of a previous
version of the progenitor (Nomoto
1984
,
1987
; Nomoto, private
communication), which is commonly termed e8.8 in the literature
and e8.8
n
in this publication. Model e8.8
n
had a degenerate core
of 1.375 M
, but a smaller hydrogen envelope (see Fig.
1
). It is a
reference case for simulations of ECSNe and was used in various
studies, focusing on explosion properties (Kitaura et al.
2006
;
Fischer et al.
2010
), nucleosynthesis (Janka et al.
2008
;Wanajo,
Janka & M
̈
uller
2011
;Wanajoetal.
2018
), effects of microphysics
(H
̈
udepohl et al.
2010
; von Groote
2014
;Radiceetal.
2017
),
and PNS kicks (Gessner & Janka
2018
). In order to compensate
for possible uncertainties in our ECSN simulations with the new
progenitor model of e8.8, we vary its explosion energy in a set
of 2D simulations, denoted by e8.8
3
, e8.8
6
, e8.8
10
, and e8.8
15
for
E
exp
=
(3
,
6
,
10
,
15)
×
10
49
erg.
Progenitor model e8.8 has the very sharp density gradient at
1
.
34 M
near the edge of its compact ONeMg-core (
ξ
2
.
5
=
5
.
7
×
10
−
6
,
ξ
1
.
5
=
8
.
0
×
10
−
6
)
1
that is characteristic of such progenitors
prior to the onset of core collapse. This steep gradient is a prominent
feature of the density (
ρ
) profile in Fig.
1
. In the same figure
we also show the
ρ
r
3
-profile as well as the electron fraction
Y
e
and the temperature
T
. All of these quantities are displayed as
functions of enclosed mass and radius. The inner
∼
0.45 M
of the
degenerate core contain iron-group nuclei and
α
-particles in NSE.
The dash–dotted, dashed, and dotted lines indicate the positions
of the NSE/ONeMg, CO/He, and He/H composition interfaces,
respectively. We define the locations of the composition interfaces,
similar to Wongwathanarat et al. (
2015
), as those positions at the
bottom of the respective layers of the star where the mass fractions
X
i
drop below half of their maximum values in the layer. The
radial positions of the composition interfaces are summarized in
Ta b l e
1
. In the top panel of Fig.
2
we present the composition of
model e8.8, where we combine all elements with mass numbers
greater than 28 into ‘iron-group’ (IG) material or nuclei in NSE.
The ONeMg-core is surrounded by very thin carbon and helium
shells (
M
C
≈
8
.
1
×
10
−
3
M
,
M
He
∼
2
.
1
×
10
−
6
M
) and a hydrogen
(H
+
He) envelope (
M
H
≈
4
.
49 M
). The total masses of the differ-
ent nuclei present in the entire pre-SN model are listed in Table
5
.
For collapse and post-bounce evolution of an ONeMg-core
progenitor, we employ in all of this paper the new progenitor model
e8.8 in 1D, 2D, and 3D simulations with the
PROMETHEUS
-HOTB
code as detailed in Section 3.4. The profile of the old progenitor
e8.8
n
is shown in Fig.
1
merely for illustration and reference.
2.2 Fe-core progenitors
As
a
second
progenitor
we
employ
a
zero-metallicity
M
ZAMS
=
9
.
6M
star, termed z9.6. It was first used by Janka
1
Because the core structure of higher mass progenitors can be characterized
by the compactness parameter
ξ
M
≡
M/
M
R
(
M
bary
=
M
)
/
10
3
km
(O’Connor & Ott
2011
), we provide values for this quantity here. However, due to the steep
density gradient outside of the core in the considered progenitors, the value
of this parameter in the context of our work is limited.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
2042
G. Stockinger et al.
Figure 1.
Profiles of the temperature (
T
), density (
ρ
), electron fraction (
Y
e
), and
ρ
r
3
for the pre-collapse progenitor models as functions of radial coordinate
(left-hand panels) and mass coordinate (right-hand panels). Indicated by dash–dotted, dashed, and dotted lines are the outer boundaries of the dege
nerate (iron
or NSE), CO, and He cores, respectively. Note the huge differences in the density and
ρ
r
3
-profiles between the progenitors with iron and ONeMg-cores, in
particular just outside of the CO core. We show the difference in the core structures of our ONeMg-core models in a zoom of the
ρ
versus
M
(
r
) profiles in the
rightmost panel.
et al. (
2012
) and was also considered in other studies such as
M
̈
uller, Janka & Marek (
2013
), Radice et al. (
2017
), and M
̈
uller
et al. (
2019
). This iron core progenitor is structurally similar to
the ECSN model. It also shows a sharp decline of the density at
the edge of its iron core, enabling low-energy explosions in 1D
(Melson et al.
2015a
;Radiceetal.
2017
). Evolved by Heger (private
communication) as an extension to Heger & Woosley (
2010
), the
pre-SN model develops an iron core of about 1.30 M
. The iron
core is surrounded by a 0
.
061 M
Si-layer, a 0
.
016 M
CO-layer,
and has a hydrogen-free helium layer of about 0.004 M
below a
0.33 M
convective H/He-layer with a hydrogen mass fraction of
∼
6 per cent, which is surrounded by a massive hydrogen envelope
of nearly 8 M
(
ξ
2
.
5
=
7
.
66
×
10
−
5
,
ξ
1
.
5
=
2
.
38
×
10
−
4
). The H/He-
layer is in the process of extended, but incomplete, second dredge-
up of the He-layer that is typical for AGB/sAGB stars and those
at the low-mass end of the CCSN domain. The star is basically
an iron-core AGB star, maybe should be called a hyper-AGB
(hAGB) star.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
3D models of supernovae from low-mass progenitors
2043
Table 1.
Shell structure of the pre-collapse progenitor models.
Model
Interface
R
shell
M
shell
t
3D
sh
,
max
E
bind
(cm)
(M
)(s)
(10
49
erg)
e8.8
NSE/ONeMg
2
.
57
×
10
7
0.45
–
ONeMg/C
8
.
48
×
10
7
1.33
0.19
E
bind
(
m
>
M
map
)
=−
5.99
C/He
1
.
09
×
10
8
1.34
0.20
He/H
1
.
21
×
10
8
1.34
0.21
Surface
8
.
43
×
10
13
5.83
4
.
1
×
10
5
z9.6
IG/Si
1
.
10
×
10
8
1.30
0.09
Si/CO
1
.
45
×
10
8
1.36
0.11
E
bind
(
m
>
M
map
)
=−
5.82
CO/He
6
.
48
×
10
8
1.37
0.40
He/H (H
<
1per
cent)
6
.
24
×
10
9
1.38
2.64
He/H (H
<
10 per
cent)
1
.
40
×
10
12
1.70
1
.
8
×
10
3
Surface
1
.
50
×
10
13
9.60
1
.
1
×
10
5
s9.0
IG/Si
1
.
24
×
10
8
1.30
0.08
Si/CO
1
.
55
×
10
8
1.33
0.30
E
bind
(
m
>
M
map
)
=−
2.63
CO/He
1
.
34
×
10
9
1.40
1.30
He/H
1
.
22
×
10
11
1.57
124.0
Surface
2
.
86
×
10
13
8.75
1
.
8
×
10
5
Notes
. The radii of the composition interfaces,
R
shell
, are defined as those positions at the bottom of the stellar layers
(e.g. CO) where the mass fractions (e.g. C
+
O) drop below half of their maximum values in the respective layer.
Progenitor model z9.6 is in the process of deep second dredge-up with hydrogen reaching basically to the bottom of the
former helium layer; hence two values are provided for the He/H interface. In all figures, we refer to the ‘high’ value,
since it is often used also in other context, e.g. for the ‘
α
-parameter’ in common-envelope (CE) studies. We also show
the mass
M
shell
contained within the corresponding radius and the post-bounce time when the outermost radius of the
forward shock of our 3D models reaches the interface.
E
bind
is the binding (i.e. internal
+
kinetic
+
potential) energy
in the progenitor star outside of
M
map
, which is the location of the final mass cut. For values see Table
3
.
As the envelope is not rich of metals, mass-loss is expected to
play only a negligible role during the star’s evolution. This leaves
the total mass of the star almost unchanged (see Table
1
). Due to its
structure it was one of the first iron-core progenitors that exploded in
fully self-consistent 3D simulations by Melson et al. (
2015a
) with
VERTEX
-
PROMETHEUS
and was also investigated by Radice et al.
(
2017
)andM
̈
uller et al. (
2019
). The result of Melson et al. (
2015a
)
provides the initial state for our investigation.
Moreover, we investigate a solar-metallicity
M
ZAMS
=
9
.
0M
star, termed s9.0, of Sukhbold et al. (
2016
). Its 1.30 M
iron core
is surrounded by a silicon shell of 0.03 M
, a carbon–oxygen
layer of 0.068 M
, and a helium shell of 0.169 M
(see Table
1
).
The hydrogen envelope extends from 1.57 M
up to 8.75 M
(
ξ
2
.
5
=
3
.
83
×
10
−
5
,
ξ
1
.
5
=
5
.
24
×
10
−
3
). In comparison to model
z9.6, model s9.0 is just slightly less evolved on its track to second
dredge-up of the He core, however, the convection is better driven by
the iron-peak opacity from the outside-in. The s9.0 progenitor was
chosen to be representative for low-mass CCSNe by Jerkstrand et al.
(
2018
), who studied the late-time nebular spectra of the supernova
(SN), by Glas et al. (
2019a
) focusing on the neutrino emission
during the explosion, and by Burrows et al. (
2019
)inalargeset
of 3D simulations. The three-dimensional exploding model for our
investigation is provided by Melson et al. (
2020
) and has also been
modelled with
VERTEX
-
PROMETHEUS
.
Although the progenitors considered in this study have very
similar ZAMS masses, their pre-collapse core structures differ
significantly (see Fig.
1
). We stress that the
ρ
r
3
-profiles of the
progenitors are decisive for the long-time evolution of the explosion
(Kifonidis et al.
2003
; Wongwathanarat et al.
2015
). The behaviour
of the
ρ
r
3
-profile yields important information on the propagation
of the shock through the stellar structure, because, according to
Sedov et al. (
1961
), positive gradients of
ρ
r
3
cause shock deceler-
ation, whereas negative gradients cause the opposite. Additionally,
variations of the shock velocity produce crossing pressure and
density gradients behind the shock front near the composition-
shell interfaces (Chevalier & Klein
1978
). Such conditions are
essential for RT instabilities as detailed in Section 5.1, assigning
them a crucial role for explaining high-velocity metal-rich ejecta
and radial mixing of heavy elements (e.g. Arnett et al.
1989a
;
Nomoto, Filippenko & Shigeyama
1990
; Wongwathanarat et al.
2015
).
The ECSN progenitor exhibits an extremely sharp drop of the
ρ
r
3
-profile just outside of the ONeMg-core. It is this drop in
density that enables fast explosions due to an early, rapid decline
of the mass accretion rate and of the corresponding ram pressure
at the SN shock (Kitaura et al.
2006
). Outside of the core, the
ρ
r
3
-profile grows monotonically as no other composition interfaces
are encountered. The z9.6 progenitor falls into the class of ECSN-
like progenitors also in this respect: Similar to the electron-capture
progenitor, the z9.6 model shows a monotonic growth of the
ρ
r
3
-
profile exterior to the CO core, where only a small step in the density
profile can be seen at the He/H interface. Model s9.0on the other
hand exhibits strong variations in its
ρ
r
3
-profile. Each interface
of different composition layers is accompanied by a negative
ρ
r
3
-
gradient close to the interface. Of particular interest are the CO/He
and He/H interfaces. These interfaces have an impact on the long-
time evolution of the explosion as will be discussed in Section 5.1.
For this paper we performed spherically symmetric (1D), ax-
isymmetric (2D), and fully three-dimensional (3D) simulations for
all progenitors beyond the moment when the shock reaches the
surface of the star. The different set-ups and approaches for these
simulations will be described in the following section.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
2044
G. Stockinger et al.
Figure 2.
Pre-collapse composition of models e8.8 (top), z9.6 (middle), and s9.0 (bottom) as a function of radius (left) and enclosed mass (right). Note
the broken horizontal axis of the right-hand panels. We combine all chemical elements with mass numbers greater than 28 into the ‘iron-group’ (IG). Th
e
dash–dotted, dashed, and dotted lines indicate the outer boundaries of degenerate (iron or NSE), CO, and He cores, respectively.
3 NUMERICAL METHODS AND PHYSICAL
SET-UP
In order to cover collapse, shock revival, and shock propagation
through the envelope and circumstellar material we employ a step-
wise approach. Core collapse, bounce, and post-bounce evolution
until the explosion is well on its way, i.e. the neutrino-dominated
phases of the first second(s) around and after core bounce, of the
iron-core progenitors are simulated with the
VERTEX
-
PROMETHEUS
code. The corresponding long-time simulations covering the shock
propagation through the star, after the onset of the explosion until
and beyond shock breakout, are conducted with the
PROMETHEUS
-
HOTB code. The explosion of the ECSN progenitor is simulated
entirely with
PROMETHEUS
-HOTB.
2
2
The reason for these different treatments is mainly technical: The
VERTEX
-
PROMETHEUS
version used for ONeMg core collapse by Kitaura et al. (
2006
),
Janka et al. (
2008
), and H
̈
udepohl et al. (
2010
) has not yet been updated
with the 3D developments and parallelization optimization applied to the
In the following sections, we describe the numerical and physical
features of the codes and the set-ups applied during the different
stages of the evolution.
3.1
VERTEX
-
PROMETHEUS
code
VERTEX
-
PROMETHEUS
is a hydrodynamics code based on an imple-
mentation of the piecewise parabolic method (PPM) of Colella &
Woodward (
1984
), coupled with a three-flavour, energy-dependent,
ray-by-ray-plus (RbR
+
) neutrino transport scheme that iteratively
solves the neutrino energy and momentum equations with a closure
determined from a tangent-ray Boltzmann solver (Rampp & Janka
version used for Fe-core collapse. However, for suitable choices of parameter
values, the explosion dynamics of ECSNe computed with
PROMETHEUS
-
HOTB is very similar to the fully self-consistent 2D
VERTEX
-
PROMETHEUS
and
COCONUT
-
VERTEX
explosion models discussed by Janka et al. (
2008
,
2012
), respectively.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
3D models of supernovae from low-mass progenitors
2045
2002
). It employs the full set of neutrino reactions and microphysics
presented in Buras et al. (
2006
) and the high-density equation of
state (EoS) of Lattimer & Swesty (
1991
) with a nuclear incompress-
ibility of K
=
220 MeV. At low densities (
ρ
≤
ρ
HD
=
10
11
gcm
−
3
),
VERTEX
-
PROMETHEUS
uses the EoS of Janka & M
̈
uller (
1995
),
which includes the contributions of photons, arbitrarily degen-
erate and arbitrarily relativistic
e
+
/
e
−
, and non-degenerate, non-
relativistic nucleons and nuclei. The relative abundances of 23
nuclear species (including some neutron-rich nuclei) are determined
by an NSE solver in regions with temperatures above
T
NSE
.Below
T
NSE
a flashing scheme is used to approximately treat nuclear
burning(seeRampp&Janka
2002
). For unshocked, collapsing
stellar matter we choose
T
NSE
=
0
.
5 MeV, and for neutrino-heated
postshock matter we take
T
NSE
=
0
.
5 MeV for the simulation of
model z9.6 and
T
NSE
=
0
.
34 MeV for model s9.0.
3
The simulations presented here are performed with a 1D gravita-
tional potential including general relativistic corrections, Case A
of Marek et al. (
2006
). The neutrino transport solver contains
corrections for general relativistic redshift and time dilation effects.
V
ERTEX
-P
ROMETHEUS
makes use of the axis-free Yin-Yang grid
(Kageyama & Sato
2004
) based on the implementation of Melson
(
2016
).
3.2
PROMETHEUS
-HOTB code
PROMETHEUS
-HOTB is based on the same hydrodynamics mod-
ule as
VERTEX
-
PROMETHEUS
and uses the implementation of the
Yin-Yang grid presented in Wongwathanarat, Hammer & M
̈
uller
(
2010a
). It employs the EoS of Lattimer & Swesty (
1991
) for high
densities above a threshold value
ρ
HD
(usually 10
11
gcm
−
3
)and
the ‘Helmholtz’ EoS of Timmes & Arnett (
1999
) for densities
below
ρ
HD
, which takes into account arbitrarily degenerate and
relativistic electrons and positrons, photons, and a set of non-
degenerate, non-relativistic nuclei. The set of nuclei consists of
neutrons
n
, protons
p
,13
α
-nuclei (
4
He,
12
C,
16
O,
20
Ne,
24
Mg,
28
Si,
32
S,
36
Ar,
40
Ca,
44
Ti,
48
Cr,
52
Fe,
56
Ni), and an additional tracer
nucleus Tr, which tracks the production of neutron-rich nuclei and
replaces
56
Ni in environments with low electron fraction,
Y
e
<
0
.
49.
These nuclear species are described as non-relativistic Boltzmann
gases. The advection of the species is treated with the consistent
multifluid advection (CMA) scheme of Plewa & M
̈
uller (
1999
).
NSE is assumed above
T
NSE
=
9
×
10
9
K and accounted for by an
NSE table including the nuclei listed above (Kifonidis, private
communication). Nuclear burning is considered at temperatures
below
T
NSE
with a 13-species
α
-network, which is consistently
coupled to the hydrodynamic modelling. At the boundary between
network and NSE we assume that all free neutrons and protons
recombine to yield
4
He. We thus add the mass fractions of
p
and
n
on
to the mass fraction of
4
He, accounting for the corresponding energy
release.
4
The P
ROMETHEUS
-HOTBcode uses a 3D gravitational
potential with the general-relativistic (GR) monopole correction
of Marek et al. (
2006
) as discussed by Arcones, Janka & Scheck
(
2007
), while higher multipoles are obtained from a solution of
Poisson’s equation as described in M
̈
uller & Steinmetz (
1995
).
3
V
ERTEX
-P
ROMETHEUS
does not apply a nuclear network for
T<T
NSE
.The
choice of a lower value of
T
NSE
permits us to follow the ejection of mass
through neutrino heating for a longer time period in model s9.0, where
the expansion velocity of the expelled matter is smaller than in z9.6, thus
facilitating nucleon recombination to
α
particles and heavy nuclei.
4
In a newer version of the code, we allow only paired free neutrons and
protons to recombine to
4
He, thus also satisfying charge conservation.
Table 2.
Summary of the PNS core parameter values used for the e8.8
model and resulting explosion energies from 1D simulations. Definitions of
these parameters in the context of our modelling approach are given in Ap-
pendix A. The explosion energy is essentially independent of dimensionality
(1D, 2D, 3D).
Model
E
exp
paR
c,f
t
0
R
ib,f
(foe)
(index)
(factor)
(km)
(s)
(km)
e8.8
3
0.03
−
31.0
×
10
−
2
27
0.1
40
e8.8
6
0.06
−
31.2
×
10
−
2
22
0.1
40
e8.8
10
0.10
−
34.0
×
10
−
1
20
0.1
40
e8.8
15
0.15
−
35.8
×
10
−
1
18
0.1
40
Different from
VERTEX
-
PROMETHEUS
,
PROMETHEUS
-HOTB uses
a three-flavour grey neutrino transport scheme as presented in
Scheck et al. (
2006
),
5
which is applicable at low and moderate
optical depths. Therefore, the high-density core of the PNS, with
amassof
M
c
=
1
.
1M
and densities well above those of the
neutrinospheric layer, is replaced by a closed (Lagrangian) inner
grid boundary at radius
R
ib
. The excised 1.1 M
PNS core is taken
into account in the gravitational potential as a central point mass.
The contraction of the PNS is mimicked by an inward movement
of the boundary radius
R
ib
, whose motion is followed by all grid
points in the computational domain. We use the contraction of
R
ib
as prescribed by Ertl et al. (
2016
). The time-dependent neutrino
luminosities at
R
ib
are imposed as boundary conditions as provided
by an analytic one-zone cooling model following Ugliano et al.
(
2012
), Sukhbold et al. (
2016
), Ertl et al. (
2016
), and Ertl et al.
(
2020
). This time-dependent treatment of the central-core region
employs five parameters (
p
,
a
,
R
c, f
,
t
0
,
R
ib,f
; see Appendix A for
definitions), which can be calibrated to yield explosions that fulfil
the constraints set by observed SNe or by fully self-consistent 3D
simulations of CCSNe. The reader is referred to Appendix A for a
more detailed description of the parametric approach and to Table
2
for the parameters of the core model employed in our work.
3.3 Collapse and post-bounce set-up in
VERTEX
-
PROMETHEUS
The collapse of the iron core progenitors is computed in 1D using the
full set of neutrino interactions until 10 ms after bounce. Thereafter,
the simulations are mapped on to the three-dimensional Yin-Yang
grid and random cell-to-cell density perturbations are imposed
with an amplitude of 0.1 per cent. The simulations employ a non-
equidistant radial grid with initially 400 zones extending to 10
9
cm,
which is refined in steps to more than 600 zones. This guarantees
a resolution
r
/
r
of better than 1 per cent at the gain radius. The
innermost 1.6 km are calculated in spherical symmetry to avoid
time-stepping constraints at the grid centre. The angular resolution
of the z9.6 model is 2
◦
. The post-bounce evolution of model s9.0is
computed with a newly implemented static mesh refinement (SMR)
scheme presented in Melson et al. (
2020
), which increases the
angular resolution to 1
◦
outside of the gain radius and to 0.5
◦
exterior
to a radius of 160 km.
The simulations with full neutrino transport are too expensive
to continue them to late post-bounce times. At
t
pb
0
.
5 s the
neutrino transport is therefore switched off and replaced by a
simplified scheme for neutrino heating and cooling, which ensures
an essentially seamless continuation with a minimum of transient
5
Some improvements to the neutrino transport module are described in
Appendix B.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
2046
G. Stockinger et al.
Table 3.
Initial positions of inner (
R
ib
) and outer (
R
ob
) grid boundaries,
baryonic mass contained within the inner boundary (
M
map
) and mapping
times
t
map
in seconds after bounce for our long-time simulations.
Model
Dim.
R
ib
R
ob
M
map
t
map
(km)
(km)
(M
)(s)
e8.8
3
2D
1000
8.4
×
10
8
1.334
2.515
e8.8
6
2D
1000
8.4
×
10
8
1.327
2.515
e8.8
10
2D
1000
8.4
×
10
8
1.319
2.515
e8.8
15
2D
1000
8.4
×
10
8
1.309
2.515
e8.8
3D
500
8.4
×
10
8
1.326
0.470
z9.6
3D
600
1.5
×
10
8
1.353
1.440
s9.0
3D
1000
2.9
×
10
8
1.351
3.140
artefacts. Details of this scheme are given in Appendix E. During
this phase of simplified neutrino treatment both model z9.6 and s9.0
are simulated with uniform angular resolution of 2
◦
.
3.4 Collapse and post-bounce set-up in
PROMETHEUS
-HOTB
During the spherically symmetric simulation of the collapse up
to core bounce,
PROMETHEUS
-HOTB uses the parametrized delep-
tonization scheme described in Liebend
̈
orfer (
2005
). The necessary
Y
e
(
ρ
)-trajectory was provided by H
̈
udepohl (private communi-
cation) from his core-collapse simulations of the ONeMg-core
progenitor e8.8
n
with
VERTEX
-
PROMETHEUS
.
For the simulation of the ECSNe progenitor we take
ρ
HD
=
10
11
gcm
−
3
and assume NSE in regions where the tem-
perature exceeds
T
NSE
=
9
×
10
9
K and apply the
α
-network for
temperatures lower than this value. After bounce
PROMETHEUS
-
HOTB employs the grey neutrino transport scheme and modelling
approach as presented in Scheck et al. (
2006
). Thus, the neutrino-
opaque central core of the PNS is excised from the computational
domain and replaced by the analytic core model of Ugliano et al.
(
2012
). In Table
2
we list the parameter values of the PNS core model
used for a set of simulations of model e8.8. We perform 1D and
2D simulations for all four sets of parameter values and choose the
e8
.
8
10
calibration as our reference case for a 3D simulation. The 1D
and 2D simulations possess a non-equidistant radial grid with 2000
zones up to a radius of
R
ob
=
2
×
10
10
cm. The 3D run has only 1400
radial zones for computational efficiency. The multidimensional
simulations are conducted with an angular resolution of 2
◦
, and the
3D simulation makes use of the Yin-Yang grid. We restrict ourselves
to a 1D gravitational potential with GR corrections (Arcones et al.
2007
) because the explosions of model e8.8 are nearly spherical.
3.5 Set-up for the long-time simulations
The simulations of the long-time evolution of all models are
computed with
PROMETHEUS
-HOTB. For this we map the final state
of a post-bounce simulation at time
t
map
on to a new computational
grid within
PROMETHEUS
-HOTB, similar to the procedure described
in Wongwathanarat et al. (
2015
). We also add the low-density
extensions to the Helmholtz EoS described therein. In Table
3
we list the times of mapping and the inner and outer radii of the
new computational domain,
R
ib
and
R
ob
, respectively. The mass
contained within
R
ib
is treated as a point mass
6
and is called
M
map
.
6
We ensure that matter at radii smaller than
R
ib
has velocities smaller than
the local escape velocity and will thus eventually contribute to the final NS
mass.
The time
t
map
is chosen such that the explosion energy has
effectively converged to its asymptotic value and a neutrino-driven
wind region has developed around the PNS, where the outflow is
essentially spherical and reaches supersonic velocity.
All long-time simulations are computed with an angular reso-
lution of 2
◦
. We use a non-equidistant (geometrically increasing)
radial grid from the inner to the outer boundary. In order to guarantee
sufficient resolution where needed, the radial grid is allowed to
move with the ejecta starting from
t
pb
∼
10 s. Gravity is accounted
for by a 1D GR-corrected potential and nuclear reactions are still
considered. When mapping the iron core models, z9.6 and e9.0, into
PROMETHEUS
-HOTB, we recombine free
n
and
p
from the freeze-
out of NSE into
4
Heunder the condition of charge conservation
and account for the energy release. Moreover, we combine all
neutron-rich nuclei formed in neutrino-heated ejecta into tracer (Tr)
material.
When mapping from the simulations of the onset of the explosion
to the follow-up simulations, the central region interior to
R
ib
(Table
3
) is removed from the computational domain. Similar to
Wongwathanarat et al. (
2015
) we prescribe an inflow boundary
condition at
R
ib
, which corresponds to the neutrino-driven baryonic
mass-loss (‘neutrino-wind’; e.g. Qian & Woosley
1996
) generated
by ongoing neutrino-energy deposition in the surface layers of
the cooling PNS. In contrast to Wongwathanarat et al. (
2015
),
we employ neutrino-wind results adopted from 1D simulations of
the explosions, seamlessly connected to the fully multidimensional
explosion simulations by choosing
R
ib
to be in the supersonic wind
region (ensuring that perturbations cannot propagate back to the
inner boundary creating artefacts) and by applying the wind data
at times when the outflow properties match closely between the
1D and the (angle-averaged) multidimensional models. This is
possible because the PNSs in 1D and multidimensional models are
extremely similar and the neutrino-emission and thus the neutrino-
driven winds also have very similar properties. For the long-time
run of model z9.6, we therefore employ neutrino-wind conditions
of a 1D simulation of this model with the
VERTEX
-
PROMETHEUS
code. This model treats PNS convection with an approach based on
mixing-length theory and exhibits neutrino-emission properties that
are hardly distinguishable from the multidimensional calculation
(see Mirizzi et al.
2016
). The time dependence of the radial velocity
v
r
, density
ρ
, and total (i.e. kinetic
+
internal) energy density
e
,
which are needed for setting the boundary condition, are shown in
Fig.
3
(solid lines) as extracted from the 1D explosion simulation of
the z9.6model at a radius of 600 km (in the supersonic wind domain).
For the 3D simulation of model e8.8we use the neutrino-wind results
from the corresponding 1D run with the same explosion energy (see
dashed lines in Fig.
3
), whereas we do not impose a wind boundary
condition in the long-time simulation of model s9.0, since the
neutrino-driven wind in this model is already very weak at the time
of mapping. Using time dependences of the boundary conditions
normalized by the initial value at the mapping time
t
map
guarantees
a smooth, seamless transition from the earlier evolution to the long-
time evolution of the explosion. Transient artefacts are thus kept
minimal.
Additionally, for the long-time simulations we include the decay
of radioactive nickel, which becomes a relevant source of energy
during late phases of the explosion. Radioactive
56
Ni (half-life
t
1
/
2
=
6
.
077 d) decays to
56
Covia electron capture (EC) decay. The
resulting
56
Conucleus is unstable (
t
1
/
2
=
77
.
23 d) and decays to
56
Feby means of electron capture and via positron decay (
β
+
). We
thus add
56
Co and
56
Fe to our set of nuclei. The respective decay
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
3D models of supernovae from low-mass progenitors
2047
Figure 3.
Time-dependent behaviour of neutrino-wind density
ρ
, radial
velocity
v
r
, and total energy density
e
, normalized to their initial values at
t
map
(see Table
3
), which defines the start of the long-time simulations
of model e8.8and z9.6. The data for model z9.6 are extracted from a
1D simulation of the explosion of the 9
.
6M
progenitor with
VERTEX
-
PROMETHEUS
(see Mirizzi et al.
2016
), evaluated at a radius of 600 km
(solid lines). For model e8.8, we use data of the respective 1D simulation
with
PROMETHEUS
-HOTB (dashed lines).
reactions are given by
EC :
e
−
+
56
28
Ni
→
56
27
Co
+
ν
e
+
γ,
EC :
e
−
+
56
27
Co
→
56
26
Fe
+
ν
e
+
γ
(81 per cent)
,
β
+
:
56
27
Co
→
56
26
Fe
+
e
+
+
ν
e
+
γ
(19 per cent)
.
The above reactions provide an energy source for the surrounding
plasma in the form of gamma radiation (
E
γ
) and kinetic energy
(
E
kin
,
e
+
) of the positrons that are produced in the
β
+
decays. We
include the annihilation energy of the positrons with electrons
(
E
ann
)in
E
γ
. The produced neutrinos escape freely. The average
energy available (including the kinetic energy of the positron in the
cobalt decay) per decay is
E
γ,
Ni
=
1
.
72 MeV,
E
γ,
Co
=
3
.
735 MeV
(Nadyozhin
1994
). A fraction of the
γ
’s may escape depending on
the (radial) optical depth
τ
(
r
) of the gas up to the stellar surface at
radius
R
∗
. This optical depth is defined as
τ
(
r
)
=−
r
R
∗
κ
γ
Y
e
(
r
)
ρ
(
r
)d
r
,
(1)
where
Y
e
the electron fraction and
κ
γ
the optical opacity. In the prac-
tical application the integral boundary
r
in equation (1) is the radial
location of a considered grid cell, and we assume that the trapped
fraction of the locally produced
γ
radiation, (1
−
exp[
−
τ
(
r
)]),
deposits its energy locally in the same cell of the computational
grid. Assuming Compton-scattering is the dominant opacity source,
we adopt a constant value of
κ
γ
=
6
.
0
×
10
−
2
cm
2
g
−
1
(Swartz,
Sutherland & Harkness
1995
). Therefore, the energy per mass
E
i
/
M
deposited by each species
i
, with mass fraction
X
i
and
nuclear mass
m
i
, into the surrounding plasma during a time-step
t
is given by
E
i
M
=
X
i
m
i
E
γ
1
−
e
−
τ
(
r
)
+
E
kin
,
e
+
,
(2)
where
X
i
=
1
−
e
−
t/t
0
,i
is the change of
X
i
during
t
with
t
0
,i
=
t
1
/
2
,i
(ln(2))
−
1
being the lifetime of species
i
,and
E
kin
,
e
+
is the
kinetic energy of the positron in the cobalt decay. The energy is
assumed to be deposited locally, and thermodynamic quantities are
self-consistently updated.
4 EVOLUTION DURING THE FIRST SECONDS
4.1 Shock propagation and explosion energetics
In the following we provide a brief overview of the most important
features of the early post-bounce evolution in the 3D simulations
of models e8.8, z9.6, and s9.0. The reader is referred to Melson
et al. (
2015a
,
2020
) for a detailed analysis of the post-bounce phase
of the iron-core progenitors in
VERTEX
-
PROMETHEUS
simulations.
Generic properties of the explosion of ECSNe are given in Kitaura
et al. (
2006
), Janka et al. (
2008
), H
̈
udepohl et al. (
2010
), Fischer
et al. (
2010
), Radice et al. (
2017
), and Gessner & Janka (
2018
).
The dynamics of our simulations for the 8.8 M
progenitor closely
resemble these previous findings. In Fig.
4
we show the evolution
of the angle-averaged radius of the SN shock and the diagnostic
explosion energy of our three-dimensional simulations during the
first three seconds after bounce. The angle-averaged shock radius
is calculated as
R
sh
=
1
4
π
R
sh
(
θ,φ
)d
,
(3)
where d
=
sin
θ
d
θ
d
φ
. The diagnostic explosion energy at all times
is given by the integral of the total (i.e. kinetic plus internal plus
gravitational) energy density in the postshock region, defined as
e
b
=
e
int
+
e
kin
+
e
grav
, over volume elements where it has a positive
value
E
exp
=
V
postshock
(
e
b
>
0)
dV
e
b
.
(4)
The sharp drop in density outside of the ONeMg-core in model
e8.8leads to an early and strong decrease of the accretion rate
and hence ram pressure at the shock. Consequently, the SN
shock expands rapidly and reaches the core/envelope boundary
(
R
He
/
H
=
1210 km), at 0
.
21 s after core bounce. This is in stark
contrast to the typical CCSN, where the high ram pressure stalls the
shock expansion at a small radius for several 100 ms. The explosion
energy starts rising steeply, fuelled by the onset of a neutrino-driven
wind, as soon as the shock leaves the ONeMg-core.
The acceleration of the shock at the core/envelope boundary is
followed by a drastic switch to deceleration after the shock passes
the lower boundary of the H-envelope (see upper right panel of
Fig.
4
). This is caused by a sudden change of the density gradient
at the core/envelope transition. Consequently, the neutrino-heated
ejecta pile up in a dense, compressed and decelerated shell behind
the SN shock.
The ECSN-like structure of model z9.6 is reflected in the
evolution of the SN shock front and in the growth of the diagnostic
explosion energy. The shock radii remain almost perfectly spherical
in both the e8.8 and z9.6 models. The acceleration of the blast wave
outside of the iron core, however, ends earlier in model z9.6 than
in model e8.8 because of the more gradual changes in the density
profile. In both cases the explosion energies also start to rise early
and saturate just after
t
pb
∼
1s.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
2048
G. Stockinger et al.
Figure 4.
Shock radius (upper left panel), shock velocity (upper right panel), and diagnostic explosion energy (lower left panel) versus post-bounce time for
all
of our 3D models. The wiggles in the shock velocity of model z9.6 at
∼
0.45 s are a consequence of small-amplitude neutron-star vibrations when the
VERTEX
neutrino transport is switched off and the heating/cooling scheme of Appendix E is switched on. Owing to an improved treatment this numerical transie
nt is
much reduced in model s9.0. We also show the total PNS kick velocities (thick lines) and the hydrodynamically induced PNS kick velocities (thin lines)
in
the lower right panel. For better visibility, the inset of this panel displays the PNS kick velocities caused by asymmetric neutrino emission due to th
e LESA
phenomenon (see text for details). The total velocities are the vector sum of the hydrodynamic PNS kick and the neutrino-induced kick. Despite nearly
equal
neutrino-induced kicks in z9.6 and s9.0 and lower hydrodynamic kick in z9.6, this latter model has a higher total kick velocity for some time, because i
ts
hydrodynamic and neutrino-induced kick directions are essentially parallel, in contrast to the situation in model s9.0. For the iron-core progenit
ors we can track
only hydrodynamic contributions to the kick after the transport module of
VERTEX
-
PROMETHEUS
is switched off at
t
neut
=
0
.
45 s and
t
neut
=
0
.
49 s for models
z9.6 and s9.0, respectively. For the ONeMg case our simplified neutrino treatment with the excised core of the PNS does not allow us to monitor the LESA
induced kick. The kick of model e8.8is scaled by a factor of 5 for better visibility.
In contrast to the z9.6 and e8.8 models, model s9.0lacks the very
steep density gradient at the edge of the Fe-core, as can be seen in
Fig.
1
. The shock can expand initially up to 180 km at
t
pb
∼
130 ms,
but then it enters a phase of recession. The arrival of the Si/O
interface at the shock and the decreasing mass-accretion rate within
the oxygen shell eventually lead to shock expansion at
∼
0
.
32 s after
bounce. Shock expansion is aided by strong convection behind the
shock front (see also Melson et al.
2020
). Similar to the results
presented in Glas et al. (
2019a
), who used the same progenitor, the
model remains convection-dominated and does not exhibit any sign
of the oscillatory growth of the SASI (Blondin et al.
2003
; Foglizzo
et al.
2007
). However, although SASI does not develop in the
simulation, the forward shock experiences large-scale deformation
with a dipole amplitude of
∼
10 per cent (compared to the angle-
averaged shock radius). This deformation is driven by big plumes
that form in the post-shock layer. Contrary to the ECSN-like models,
strong anisotropic, non-radial mass flows persist around the PNS
during several seconds after bounce. Continuous mass accretion
on to the PNS through narrow funnels delays the emergence of
the spherical neutrino-driven wind, which is why we needed to
continue the simulations including PNS and neutrino treatment for
more than three seconds after bounce. Explosion energy and shock
velocity in model s9.0remain considerably lower than in the other
two progenitors (see Fig.
4
).
4.2 Neutrino emission properties
In the following we present the neutrino emission properties of
our two 3D simulations performed with the
VERTEX
-
PROMETHEUS
code. Fig.
5
displays the neutrino luminosities (defined as 4
π
-
integrated energy fluxes) and mean neutrino energies (defined as
ratio of angle-averaged energy density to number density) for
ν
e
, ̄
ν
e
and heavy-lepton neutrinos
ν
x
, as well as the three lowest order
multipoles (monopole, dipole, and quadrupole) of the electron-
neutrino lepton-number flux for models z9.6 and s9.0 as functions
of post-bounce time (z9.6 left column, s9.0 right column). All
quantities are evaluated at 400 km (transformed to an observer
frame at rest at infinity). The formulas for the spherical harmonics
decomposition are provided in Appendix C.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
3D models of supernovae from low-mass progenitors
2049
Figure 5.
Neutrino luminosities (top), mean energies (middle), all averaged over all viewing directions, and lowest order multipole moments (
A
0
,
A
1
,
A
2
for monopole, dipole, and quadrupole) of the electron-neutrino lepton-number flux (bottom) as functions of time for models z9.6 (left) and s9.0 (righ
t). All
quantities are evaluated at 400 km, transformed to an observer frame at rest at infinity. With
ν
x
we denote one species of the heavy-lepton neutrinos. Note that the
dipole moment,
A
1
, plotted here is one third of the dipole amplitude of the lepton-number flux defined by Tamborra et al. (
2014a
) (see Appendix C for details).
In the case of z9.6the luminosities of all three species become
very similar after only
∼
180 ms, signalling the end of PNS accretion
caused by the quick onset of the explosion and the rapid shock
expansion. In contrast, PNS accretion continues at a significant rate
until roughly 350 ms in model s9.0. Only afterwards the luminosities
in this model converge to nearly the same level, mirroring the
characteristic trend when the cooling emission of the newly formed
NS begins. Overall, the time evolution and the values of the neutrino
luminosities and mean energies of both models are very similar to
each other, consistent with the nearly equal masses of the PNSs in
both cases (see Table
4
). Model s9.0 exhibits additional accretion
emission between
∼
150 and
∼
350 ms, which enhances the
ν
e
and
̄
ν
e
luminosities slightly and drives a continuous increase of the
mean energies of
ν
e
and ̄
ν
e
. However, the trend does not persist for
long enough (and the PNS does not gain enough mass) to reach a
crossing of the electron antineutrino energy with the heavy-lepton
neutrino energy as reported by Marek, Janka & M
̈
uller (
2009
)for
more massive progenitors with more massive PNSs.
The bottom panels of Fig.
5
demonstrate that both models develop
a lepton-number emission dipole that dominates the higher order
multipoles and can reach and even exceed the monopole for longer
periods of time (i.e. hundreds of milliseconds, see also Tamborra
et al.
2014b
). Note that the dipole amplitude displayed in Fig.
5
is, by normalization, one-third of the dipole amplitude considered
by Tamborra et al. (
2014b
) (see also Appendix C). A long-lasting
lepton-number emission dipole, whose direction is nearly stable
or migrates only very slowly, is a characteristic feature of the
LESA (lepton-emission self-sustained asymmetry) phenomenon
that was first witnessed in 3D
VERTEX
-
PROMETHEUS
simulations
with neutrino transport by Tamborra et al. (
2014b
) and Janka et al.
(
2016
). This striking phenomenon has meanwhile been confirmed
with fully multidimensional instead of ray-by-ray neutrino transport
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020
2050
G. Stockinger et al.
Table 4.
Overview of PNS properties in our multidimensional models at
t
map
.
Model
Dim.
t
map
E
exp
v
tot
NS
v
hyd
NS
v
ν
NS
J
NS
/10
45
θ
vJ
α
ej
α
ν
M
b
R
NS
M
map
M
g
P
NS
(s)
(10
50
erg) (km s
−
1
)(kms
−
1
)(kms
−
1
)(cm
2
gs
−
1
)(
◦
)
(per
cent)
(per
cent)
(M
)(km)(M
)(M
)(s)
e8
.
8
3
2D
2.515
0.3
1.55
1.55
–
1.78
–
1.154
–
1.323
49.85
1.334
1.216
4.18
e8
.
8
6
2D
2.515
0.6
0.94
0.94
–
2.75
–
0.041
–
1.316
50.22
1.327
1.210
2.69
e8
.
8
10
2D
2.515
1.0
0.13
0.13
–
4.19
–
0.004
–
1.308
50.57
1.319
1.203
1.75
e8
.
8
15
2D
2.515
1.5
0.59
0.59
–
1.67
–
0.011
–
1.299
50.86
1.309
1.195
4.37
e8.8
3D
0.470
1.0
0.44
0.44
–
0.70
90.0
0.004
–
1.307
50.57
1.326
1.210
10.58
z9.6
3D
1.440
0.86
34.90
10.16
24.89
2.55
45.4
4.623
1.354
1.340
20.98
1.353
1.231
2.96
s9.0
3D
3.140
0.48
40.87
28.45
26.46
8.05
31.3
10.02
1.178
1.350
19.58
1.351
1.230
0.94
Notes
. All values are given at the end of our explosion simulations with neutrino treatment (
t
map
).
E
exp
is the diagnostic explosion energy, which is essentially
identical to the final explosion energy because of the small envelope binding energy (see Table
5
).
v
tot
NS
is the total NS kick velocity resulting from the
hydrodynamic (
v
hyd
NS
) plus the neutrino-induced (
v
ν
NS
) contributions. We measure the contribution of neutrinos to the total kick until the neutrino transport is
switched of at
t
neut
=
0
.
45 s and
t
neut
=
0
.
49 s after core bounce for models z9.6 and s9.0, respectively. Hydrodynamic contributions and total kick velocities are
monitored until
t
map
.
J
NS
is the total angular momentum transported to the PNS through a radius of 100 km until
t
map
.
θ
vJ
is the angle between the direction of
the total kick velocity
v
v
v
and the direction of
J
J
J
.
α
ej
and
α
ν
are the final hydrodynamic and neutrino anisotropy-parameters, respectively (see Appendix D). For
the iron-core progenitors we average
α
ν
over the time when the emission dipole is largest until the end of the simulation.
M
b
is the baryonic PNS mass, which
is defined as the enclosed mass within the radius
R
NS
, at which the density drops below 10
11
gcm
−
3
. Note that for the simulations of the e8.8 progenitor,
R
NS
is
determined by the chosen parameters of the inner grid boundary (see Table
2
).
M
map
(see also Table
3
) is the central mass contained within the excised region,
from which we calculate the gravitational mass
M
g
for a PNS radius of 12 km (see Appendix D).
P
NS
is the PNS spin period at the end of our post-bounce
simulations, assuming a final NS radius
R
NS
of 12 km, angular momentum conservation, and a gravitational mass of
M
g
.
by O’Connor & Couch (
2018b
), Glas et al. (
2019b
), and Vartanyan
et al. (
2019
).
Neither in model z9.6 nor in s9.0 does SASI play a role and,
in addition, model z9.6 explodes after a very brief period of post-
bounce accretion. In both models the PNS emission is therefore
not masked by asymmetric accretion due to SASI shock sloshing
or spiral motions. Because accretion in particular in model z9.6
does not contribute to the neutrino emission at any significant level
after the onset of the explosion, its lepton-number flux asymmetry
is merely determined by asymmetric convection inside of the PNS,
where outward transport of lepton number is strongly suppressed
in the anti-LESA direction. Therefore the lepton-number flux can
even be negative in the anti-LESA direction (see Fig.
6
, two upper
left panels). This means that one hemisphere of the PNS radiates a
greater number of ̄
ν
e
than
ν
e
, whereas there is the usual excess of
ν
e
number loss from the other hemisphere. The models, in particular
z9.6with no long-lasting accretion, therefore confirm that LESA is
a phenomenon primarily generated by hemispherically asymmetric
convection in the interior of the PNS, well below the neutrinosphere;
the reader is referred to the more detailed discussion by Glas et al.
(
2019b
), where also a 3D simulation with successful explosion of
the 9.0
M
progenitor is evaluated.
7
It is also interesting to note that
the lepton-emission multipole (
=
0) that grows fastest and thus de-
velops high amplitudes of asymmetry first is the one of order
=
4,
which is also in line with the results reported by Glas et al. (
2019b
).
Interestingly, there are phases in both of our models when the
dipole and quadrupole amplitudes become similar (Fig.
5
). This
is also suggested by the Aitoff projections for time
t
pb
=
0
.
45 s
near the end of our simulation of z9.6 (Fig.
6
). At this late time
a quadrupole pattern is superimposed on a hemispheric dipole
asymmetry, whereas at
t
pb
=
0
.
30 s a much cleaner dipole is present.
This reflects the evolution of the multipole amplitudes visible for
7
In contrast to the models of Glas et al. (
2019b
), in which a 1D core of
10 km was used, the 3D simulations discussed here were computed with a
very small central 1D core of only 1.6 km radius. Both sets of simulations
show very similar LESA features, which means that the 10 km core had no
relevant influence on the previous results.
z9.6 in the bottom left panel of Fig.
5
. The late-time rapid growth
of the dipole in the lepton-number emission of model s9.0 at
t
pb
>
0
.
43 s (bottom right panel of Fig.
5
) is an apparent feature
caused by a transient phase of a reduced lepton-emission dipole
between
∼
0.35 and
∼
0.45 s. This reduction is a consequence of
an accretion asymmetry that channels matter towards the PNS
predominantly in one hemisphere, whereas a mighty outflow in the
form of a huge, rising bubble develops in the opposite hemisphere
(see the bottom panels of Figs
7
and
8
). The main downflow
direction is misaligned with the LESA dipole direction by roughly
90
◦
and fuels enhanced lepton-number emission by the one-sided
accretion. This enhanced lepton-number emission combined with
the displaced LESA dipole strengthens the quadrupole component
of the electron lepton-number emission and at the same time
weakens its dipole. The remaining dipole level in the interval of
[0.35 s, 0.45 s] signifies that the LESA dipole is stronger than the
emission of lepton number by accretion. The dipole amplitude
recovers to the full LESA strength after
t
pb
∼
0
.
45 s because the
mass-accretion rate on to the PNS declines continuously.
The Aitoff projections show that not only the electron-neutrino
lepton-number flux exhibits large-scale (low-order multipolar)
asymmetries, but also the
ν
e
plus ̄
ν
e
energy flux as well as the total
neutrino energy flux (the summed energy fluxes of
ν
e
plus ̄
ν
e
plus
four times that of
ν
x
). However, while the directional variations of
the lepton-number flux can be several times bigger than the angular
average of this quantity, the
ν
e
plus ̄
ν
e
energy flux varies only
within roughly 6–8 per cent and the total neutrino energy flux even
less within about 4–6 per cent (Fig.
6
). These directional variations,
in particular the hemispheric asymmetries, have consequences for
NS kicks and the neutron-to-proton ratio in the ejecta. This will be
discussed in the following Sections 4.3 and 4.4.
4.3 Neutron star kicks by asymmetric mass ejection and
neutrino emission
In order to compare the final state of the post-bounce simulations
of our model set, we present in Fig.
7
planar slices showing the
entropy at the times
t
map
when we start our long-time simulations.
MNRAS
496,
2039–2084 (2020)
Downloaded from https://academic.oup.com/mnras/article/496/2/2039/5857660 by California Institute of Technology user on 03 September 2020