Large-Scale Mixing in a Violent Oxygen-Neon Shell Merger Prior to a Core-Collapse Supernova

Large-Scale Mixing in a Violent Oxygen-Neon Shell Merger Prior to a Core-Collapse Supernova

[    Bernhard Müller    [    Tobias Melson    [
00, 00, 201900, 00, 201900, 00, 2019
00, 00, 201900, 00, 201900, 00, 2019
00, 00, 201900, 00, 201900, 00, 2019

We present a seven-minute long -3D simulation of a shell merger event in a non-rotating supernova progenitor before the onset of gravitational collapse. The key motivation is to capture the large-scale mixing and asymmetries in the wake of the shell merger before collapse using a self-consistent approach. The geometry is crucial as it allows us to follow the growth and evolution of convective modes on the largest possible scales. We find significant differences between the kinematic, thermodynamic and chemical evolution of the 3D and the 1D model. The 3D model shows vigorous convection leading to more efficient mixing of nuclear species. In the 3D case the entire oxygen shell attains convective Mach numbers of , whereas in the 1D model, the convective velocities are much lower and there is negligible overshooting across convective boundaries. In the 3D case, the convective eddies entrain nuclear species from the neon (and carbon) layers into the deeper part of the oxygen burning shell, where they burn and power a violent convection phase with outflows. This is a prototypical model of a convective-reactive system. Due to the strong convection and the resulting efficient mixing, the interface between the neon layer and the silicon-enriched oxygen layer disappears during the evolution, and silicon is mixed far out into merged oxygen/neon shell. Neon entrained inwards by convective downdrafts burns, resulting in lower neon mass in the 3D model compared to the 1D model at time of collapse. In addition, the 3D model develops remarkable large-scale, large-amplitude asymmetries, which may have important implications for the impending gravitational collapse and the subsequent explosion.

stars:massive – convection – hydrodynamics – turbulence – supernovae: general
Corresponding author: Naveen Yadavny@MPA-Garching.MPG.DE

0000-0002-4107-9443]Naveen Yadav \move@AU\move@AF\@affiliationMax Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany


Astrophysics Research Centre, School of Mathematics and Physics, Queen’s University Belfast, Belfast, BT7 1NN, UK \move@AU\move@AF\@affiliationMonash Centre for Astrophysics, School of Physics and Astronomy, Monash University, Victoria-3800, Australia

0000-0002-0831-3330]Hans Thomas Janka \move@AU\move@AF\@affiliationMax Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany


Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany

0000-0002-3684-1325]Alexander Heger \move@AU\move@AF\@affiliationSchool of Physics & Astronomy, University of Minnesota, Minneapolis, MN-55455, USA


Center for Nuclear Astrophysics, Department of Physics and Astronomy, Shanghai Jiao-Tong University, Shanghai-200240, P. R. China \move@AU\move@AF\@affiliationJoint Institute for Nuclear Astrophysics, 1 Cyclotron Laboratory, National Superconducting Cyclotron Laboratory, Michigan State University, East Lansing, MI-48824-1321, USA

1 Introduction

Multi-dimensional effects in the late burning stages of massive stars have recently garnered considerable interests for a number of reasons. Whereas the mixing-length theory (MLT) of convection (biermann_32; boehm_58) as used in one-dimensional (1D) stellar evolution models provides a very good description of mixing in the interior of convective zones for the late, neutrino-cooled burning stages, it has long been speculated that additional phenomena such as turbulent entrainment (fernando_91; strang_01; meakin_07_a; spruit_15) and the excitation of internal waves could significantly affect shell growth, mixing, and angular momentum transport (cantiello_14; fuller_15) in a manner that is not captured by current 1D stellar evolution models.

Since the early 1990s, various groups have attempted to study the late convective burning stages using multi-dimensional simulations to investigate such effects. Since the seminal early work in two dimensions (2D; bazan_94; bazan_98; asida_00) and three dimensions (3D; kuhlen_03; meakin_06; meakin_07_a; meakin_07_b), additional convective boundary mixing has indeed been consistently observed in many modern simulations (mueller_16c; cristini_17; jones_17; andrassy_18) and appears to be well captured by the semi-empirical entrainment laws familiar from terrestrial settings. The long-term impact of such extra mixing on the evolution of massive stars remains more unclear, however. One major caveat concerns the duration of the simulations, which are currently limited to periods far shorter than the thermal adjustment time scale; and it has also been argued that the predicted extrapolation of entrainment rates might in many cases not qualitatively alter stellar structure over secular time scales (mueller_16b).

In some highly dynamical situations, however, multi-dimensional effects may result in qualitatively different behaviour compared to 1D stellar evolution models. Such situations often occur when material entrained across shell boundaries burns violently, which can lead to strong feedback on the dynamics of the flow; proton ingestion episodes as studied in 3D by several groups (herwig_11; stancliffe_11; herwig_14; woodward_15) are one noteworthy example. Such violent episodes of convective boundary mixing are also interesting because there is actually a potential fingerprint for the multi-dimensional dynamics in the stellar interior in the form of the nucleosynthesis enabled by turbulent entrainment (herwig_11; jones_16b; ritter_18).

Such dynamical mixing events may also occur during the late stages of burning in massive stars. 1D stellar evolution models already show that mergers between O, Ne, and C shells are quite commonplace (sukhbold_14, see, e.g., their model in their Figure 8); they are facilitated by relatively small buoyancy barriers between these shells. What is particularly intriguing is that these mergers very often occur only a few turnover time-scales before collapse, as pointed out by the systematic study of collins_18, who found such late mergers for about of all progenitors between and . This implies that the merged shells may be caught in a highly dynamical state at the time of the supernova. It is therefore conceivable that traces of such a merger remain imprinted in the explosion geometry. If this is the case, such mergers could help explain asymmetric features in supernova remnants such as the broad Si-Mg rich “jet”–like features in Cassiopeia A (delaney_10; isensee_10; grefenstette_14; grefenstette_17). 1D stellar evolution models with parameterized mixing also suggest that such shell mergers could result in very characteristic nucleosynthesis.

The recent surge of interest in the 3D dynamics of the late burning stages has mostly focused on the potentially beneficial role of convective seed perturbations in the supernova mechanism (mueller_15a; mueller_17; couch_15) and on slow, steady-state entrainment (meakin_07_a; meakin_07_b; mueller_16c; jones_17; cristini_17; andrassy_18), but such dynamical shell mergers have not yet been investigated thoroughly. mueller_16b reported a late breakout of convection across shell boundaries in a thin O burning shell that was reignited shortly before collapse in an progenitor, but with insufficient time before collapse for a genuine merger of the O and Ne shell. More recently, mocak_18 observed the merging of shells in a 3D simulation of a star, but their model was restricted to a small wedge of and not evolved until core-collapse. Moreover, the merger occurred already within the initial numerical transient before a convective steady state developed.

Here, we present the first 3D simulation over the full solid angle of a fully developed large-scale merger of a silicon and neon shell up to the onset of core-collapse. Our simulation follows the last seven minutes of convective O and Ne shell burning in an progenitor (mueller_16a) from an early stage with clearly separated convective shells to the point where the buoyancy jump between the O and Ne shells is reduced sufficiently to trigger a violent merger that is still not fully completed at collapse. In this paper, we focus on analyzing the flow dynamics and mixing during the merger and the comparison of the 3D model to the corresponding 1D stellar evolution calculation. Implications of the shell merger for the ensuing supernova, though a major motivation for the current simulation, will be discussed in future work.

The layout of the paper is as follows. The initial model, which has been derived from a 1D simulation is described in Section 2. In Section 3, we provide important details of the numerical method and the microphysics used for the 3D simulation. In Section 4 we describe in detail the 3D model and compare it with the 1D model. In Section 5 we briefly discuss the pre-supernova model. In Section 6 we provide our conclusions and discuss the impact of our results on the pre-supernova structure and the core collapse.

2 1D Progenitor Model


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefProperties of the initial 1D model. The top panel shows the mass fractions of some relevant -elements, the specific nuclear energy generation rate, and the specific neutrino cooling rate with at the base of the O burning shell, whereas in the Fe/Si core. The middle panel shows the radial velocity and the convective velocity according to MLT. The Si/Fe core is slowly contracting whereas the O shell is slowly expanding. The O burning shell has larger convective velocity compared to the Ne and C burning shells, and the Fe/Si core is convectively inert. The bottom panel shows the density and entropy profiles. Entropy jumps correspond to composition interfaces; large relative jumps in entropy means high stiffness of such boundaries against convective entrainment. Vertical black dashed lines mark the shell which is simulated in 3D.

We consider a progenitor with solar metallicity and a zero-age main sequence mass of . The 1D model has been calculated using the Kepler stellar evolution code (heger_10) as outlined in mueller_16c. The 1D model is mapped to the 3D grid before the onset of core collapse. Figure 2 shows the relevant properties of the initial 1D model. The top panel shows the mass fraction profiles for key nuclei, the specific nuclear energy generation rate, and the specific neutrino cooling rate. The shells of principal interest in this paper are the O-rich shells between (). These include a convective O burning shell333Following common practice in stellar evolution literature, we label the convective shells by the fuel that burns at their base, unless we explicitly discuss the composition in more detail. with ashes of Si and S, a convective Ne burning shell between , which is separated by a non-convective layer from a convective C burning shell with ashes of O and Ne between . The O and Ne shells are separated by a thin interface at . The Si shell and the Fe core inside are practically inert, with neutrino cooling dominating over burning. The middle panel shows the profiles of the radial velocity, , and the convective velocity, , calculated using the same choice of dimensionless coefficients as (mueller_16b). The convective velocity near the base of the O shell reaches up to . The Ne and C shells are relatively quiet, having convective velocities close to and , respectively. The bottom panel displays the density and (specific) entropy profiles. The entropy profile shows the three aforementioned convective regions as flat sections. Significant entropy jumps of and (here is the Boltzmann constant) close to the bottom (just above the Fe/Si core) and close to the top (just below the He shell) of the three shells make the convective boundaries quite stiff, with the inner boundary of the oxygen shell even being visible as a discontinuity in the density profile. By contrast, the entropy jump between the oxygen and neon shells is small, and the neon and carbon shells are separated by a stable region with gradually increasing entropy, but not by a discontinuity.

3 Numerical Methods: 3D Model

\H@refstepcounter table \hyper@makecurrenttable\hb@xt@ Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSetup details for the 3D simulation.

Quantity Option/Value aRadial grid Geometric Cells in radial direction, 450 bCells in polar direction, 48 bCells in azimuthal direction, 140 cNo. of ghost cells (), 8 cNo. of buffer cells (), 2 Angular resolution dInner boundary at , dOuter boundary at , Inner boundary condition, radial refer to Section 3.1 Outer boundary condition, radial refer to Section 3.1 Gravitational potential Newtonian, Spherical Simulation time ePerturbation amplitude () Neutrino cooling Yes Equation of state Helmholtz EoS

, where . The number of angular cells are specified per patch of the Yin-Yang grid. For each grid patch (Yin and Yang). The inner and outer boundaries are moved according to the trajectory obtained from 1D model. , where is a uniformly distributed random number in the interval .

The 3D simulation uses the Prometheus code. The simulation employs a moving (but non-Lagrangian) grid in the radial direction and an axis-free Yin-Yang grid (kageyama_04; wongwathanarat_10a) as implemented in melson_15a with a uniform angular resolution () in the and directions. For the purpose of this simulation we changed the -order reconstruction originally used in the piecewise-parabolic method (PPM) of colella_84 to a -order extremum-preserving reconstruction (colella_08; sekora_09). We use the Helmholtz equation of state (EoS) described in timmes_99 and timmes_00, which is thermodynamically consistent to high accuracy. Nuclear burning is treated using a 19-species -network (weaver_78). An accurate modeling of silicon burning is avoided as its nuclear burning in the quasi-statistical equilibrium (QSE) regime involves tracking a very large number of nuclear species making it computationally expensive. Energy loss by neutrinos (itoh_96) is included as a local sink term.

3.1 Boundary Conditions

We excise the core inside as most of the nuclear energy generation is limited to a small region at the base of the O shell (described in Section 2) and the Si/Fe core is relatively inert. The region outside of is also excluded from the simulation as it remains dynamically disconnected from the interior during the short 3D simulation. In the 1D model the Si/Fe core cools via neutrino emission and contracts. For consistency, the radial boundaries of the computational domain are moved in step with the 1D model (movement of the outer boundary is negligible in practice). Thus, the 3D simulation covers 1) a small part of the Si shell as a “buffer” below the O shell, 2) the entire O, Ne, and C shells as the “region of interest” layers, and 3) a small part of the He shell. We emphasize that the steep entropy step at the Si/O interface is about away from the inner grid boundary (see Figure 2). This convectively stable interface ensures that the convective mass motions in the simulation volume remain radially separated from the grid boundary. For PPM reconstruction, we impose reflective boundary conditions for velocity at both the inner and outer radial boundary and extrapolate all thermodynamical quantities assuming adiabatic and hydrostatic stratification. We strictly enforce zero advective fluxes across the boundaries before updating the conserved variables.

3.2 1D-to-3D Mapping

The 3D simulation is initialized using density, temperature, and mass fraction profiles taken from the 1D model. In order to break the spherical symmetry we introduce small deviations from hydrostatic equilibrium by perturbing the density field such that , where is a uniformly distributed random number in the interval and the amplitude .

The Newtonian gravitational potential is computed in the monopole approximation; the appropriate contribution from the excised inner core is added. The details of the simulation setup are summarized in Table 3.

4 Simulation Results: 1D vs 3D

In the following sections we discuss the kinematic, thermodynamic, and chemical evolution of the 3D model and compare it to the 1D model. We start by describing the convective instability of the simulated shell as a function of time.

4.1 Convective stability


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefLeft panel: Space-time plot showing the spherically averaged Brunt-Väisälä frequency in the 3D run. Regions with negative are convectively stable and regions with positive are convectively unstable according to the Ledoux criterion. Initially the active O shell (), the active Ne shell (), and the active C shell () are clearly separated by stable layers with positive gradients in the initial entropy profile (blue curve for entropy vs. mass see axis on top of panel). The barrier between the O and Ne shells at disappears after due to the increasing entropy in the O shell. The result is a large merged convective layer with a slightly negative entropy gradient between at collapse (black curve for entropy vs. mass see axis on top of panel). Right panel: Space-time plot showing the Brunt-Väisälä frequency for the 1D model. The convectively stable interfaces in the 1D model are very narrow. The stability properties in most of the O shell remain unchanged during the course of evolution, except minor changes in convectively stable Zone I close to the end. The results are shown after the model has relaxed in the first . Note: The region below (the O burning shell) has in both 1D and the 3D case. It has been masked in the plots to get around the limited dynamic range of the color axis.

The Ledoux criterion can be used to test the stability of a mass element against convective overturn. The Brunt-Väisälä frequency is given by


where is the gravitational acceleration and


is the adiabatic index. We use the following equation from buras_06a


The conversion of Equation (1) to Equation (4) is exemplified in Appendix A of mueller_16b444heger_05 show an alternative approach (used in Kepler stellar evolution code) which can be used to extract the contribution of composition to without explicitly tracking composition.. For unstable modes (), is the growth rate, and for stable modes (), is the negative of the oscillation frequency. Therefore, according to the sign convention adopted in this paper corresponds to convective instability.

\H@refstepcounter table \hyper@makecurrenttable Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefConvectively stable zones according to the Ledoux criterion at the beginning of the 3D simulation.

Zone Location (inner edge) Width I 2.4 0.1 IIa 3.7 0.4 III 4.4 0.1

Zone II in the 1D model shows a strong and narrow barrier () close to .

Figure 4.1 (left panel) shows the spherically averaged Brunt-Väisälä frequency profiles as a function of time for the 3D model. is defined as


which is different from calculated using spherically averaged , etc. The plot shows three convectively stable layers (zones IIII at masses of , , and respectively) sandwiched between more voluminous convectively unstable regions. The location ( and ) and width ( and ) of these zones at are listed in Table 4.1 in order of increasing mass coordinate. The initial and final entropy profiles (blue and black lines in left panel of Figure 4.1) show the disappearance of entropy jumps between the convectively stable and unstable regions. Such entropy steps represent barriers that the convective eddies cannot cross. As burning at the bottom of the O shell increases the entropy, the stabilizing gradient in Zone I gradually becomes weaker, and disappears altogether after , allowing the O and Ne shell to merge into one large convection zone extending between . The disappearance of Zone I is mirrored by the disappearance of the entropy jump at . This increase in entropy towards collapse is characteristic of the late burning stages prior to collapse, when neutrino cooling is too slow to balance thermal energy deposition due to rapid burning in the contracting shells. Zone II also becomes narrower after , and internal gravity waves can transport more energy across it. Different from Zone I, this stable barrier is not completely eliminated prior to collapse, however. The right panel of Figure 4.1 shows for the 1D model. The convectively stable zones present in the 1D model are much thinner and remain unchanged over the course of the simulation. In the last there is a reduction in the strength of Zone I. At the same time, Zone II and Zone III are totally undisturbed in the 1D model.

4.2 Flow Dynamics

4.2.1 Growth of Density and Velocity Fluctuations


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSlices ( plane) showing the density fluctuations at different times. The physical size of the region displays ranges from at early time to at late time. The time sequence demonstrates the transition from initially well-separated convective shells with strong density fluctuations at the boundaries from overshooting to a merging of the O and Ne shells (panels at and ), and eventually the C shell (panels at and ). Towards the end large-scale asymmetries dominate the convective flow (see text for details). The green circle marks the inner boundary of the computational domain. Minimum and maximum values in the plane are written in the top left corner of each panel. Each panel shows data only inside a spherical region bounded by a convectively stable layer (Zone I at , Zone II for , , and Zone III for and ) for clarity of presentation.

Consider the density fluctuations which are defined as


where is the average density evaluated over a spherical shell. Figure 4.2.1 shows the density fluctuations on an slice at different times. The panel at shows the development of plumes with density fluctuations at a level of few percent. These plumes are contained in the region below the first convectively stable layer (Zone I) interior to as described in Section 4.1; we refer to these as “primary plumes”. Primary plumes overshoot into the stably stratified layer above and are decelerated, creating “hot spots” in the density fluctuations as seen in both panels at and . They also excite internal gravity waves which transport energy across the stable zone and thus perturb the region directly above it, creating “secondary” plumes. In due course the mass entrainment caused by interfacial shear instabilities driven by convection scour material off the stable interface (strang_01, see Figure 1 and 3 there). As a result of the entrainment and mixing the stabilizing gradient in Zone I ceases to exist around (see Figure 4.1 left panel). This initiates the formation of a large convective region extending from the base of the O shell to the base of Zone II. The eddies overshooting into Zone II are decelerated, creating new hot spots close to , At this point, global asymmetries develop in the flow. The panel at in Figure 4.2.1 shows an asymmetric feature with a size of , which characterizes the phase of vigorous convective activity. Close to the end of simulation (panels at and ) the density fluctuations in the region within reach a level of and are again of smaller scale structure, while in the region between and they are as large as with a marked dipolar asymmetry.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefMinimum (dash-dotted line) and maximum (thick lines) values of density fluctuations as functions of mass coordinate. The arrows in the right panel mark the locations of convectively stable zones. The right panel shows the elimination of the convectively stable Zone I (seen by large density fluctuations due to overshooting) in the process of the shell merger. In the left panel, we see that the fluctuations in the Si shell travel inwards in mass coordinate (as a consequence of slow entrainment) with time and their amplitudes grow, but they never reach the inner grid boundary at . Blue arrows (right panel) mark the positions of Zone I and Zone II.

Figure 4.2.1 provides another perspective on the flow dynamics by showing the minimum and maximum of the fractional density fluctuations on spherical mass shells as functions of mass coordinate. One clearly sees (right panel) how the density fluctuations at shell interfaces initially decrease in magnitude and spread out towards larger as entrainment whittles down the stabilizing entropy gradients at the shell interfaces. After the shells have merged, the density fluctuations grow considerably, especially in the outer part of the convective region. Figure 4.2.1 (left) shows that the inner convective boundary moves towards lower by entrainment; here the density fluctuations in the boundary layer actually become stronger with time as the overall violence of convection increases. We also point out that the density fluctuations never touch the inner grid boundary, which confirms the proper choice of its location at .


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSlices in the - plane showing radial velocity fluctuations (in ) at different times. The panels till show the initially slow but steady build-up of the convective velocities. After the merger of the O and Ne shells, convection becomes considerably more violent, and large-scale flow patterns develop (see the text for details). The green circle marks the inner boundary of the computational domain. Minimum and maximum values in the plane are specified in the top left corner of each panel. Each panel shows data only inside a spherical region bounded by a convectively stable layer (Zone I at , Zone II for , , and Zone III for and ) for clarity of presentation.

We next consider the turbulent velocity fluctuations. Because the model is initially non-rotating, we do not include any non-radial components in the mean flow and decompose the velocity field as (see Appendix A)


where is the Favre-averaged radial velocity. The fluctuating component of the radial velocity is therefore given by


Figure 4.2.1 shows the radial velocity fluctuations corresponding to the snapshots shown in Figure 4.2.1. The panel at shows well-developed plumes extending from the O burning region at to the base of Zone I at . The panel at shows the secondary plumes extending from to . The primary and secondary plumes are physically separated from each other by the Zone I. Also, the secondary plumes have lower velocities (by a factor of ) compared to the primary plumes. These plumes are driven by the active Ne shell and develop later because of a lower Brunt-Väisälä frequency. The velocity difference is also apparent in panel at which is the stage when the stable Zone I has already disappeared.The panel at shows emergence of a large scale flow spanning the entire region from base of the oxygen burning layer to the base of Zone II. At this stage, convection becomes very vigorous and the magnitude of velocity fluctuations increases from to more than (panels at and ) within the next . The panel at shows the presence of a large-scale dipolar asymmetry in the flow. Convective downdrafts in the innermost region reach velocities in excess of aided by rapid contraction of the Fe/Si core.

4.2.2 Convection Length Scale


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefVelocity correlation functions (defined by Equation (9)) near the centers of layers ah (listed in Section 4.2.2) at various times. The mass and radius coordinates of these points are written in the corresponding panels. Vertical green lines mark the boundaries of these layers. The correlation function at is shown as filled grey curve in each panel. Placing the panels d and h together (aligning them along the green lines) shows two well developed convective layers with a contact somewhere inside Zone II. See the text for details.

During the course of evolution, the convective eddies gradually grow from small to large length scales. The longitudinal two-point correlation function gives a qualitative measure of the size of the convective region. It is defined as (Equation (15) of mueller_16c)


where we have used the Reynolds-averaged velocities (see Appendix A). The correlation length at a radius is defined by integrating the correlation function


between and , which are usually taken to be and respectively. meakin_07_b define these points such that . In the present case, estimating the correlation length using the equation above is non-trivial because of the presence of multiple convectively stable zones. For this analysis we divide the simulation domain into multiple layers according to their initial convective stability and the burning processes:

  1. , convectively inert Si shell,

  2. , convective O shell,

  3. , Zone I,

  4. , convective Ne shell,

  5. , Zone II,

  6. , convective C shell,

  7. , Zone III,

  8. , He shell.

Figure 4.2.2 shows the correlation function evaluated at roughly the middle points of these layers. We have marked the boundaries of these layers by vertical lines. In panel a the correlation function drops rapidly from unity with , resulting in at all times (). Thus the Si layer underneath the O burning zone stays non-convective and dynamically disconnected from the overlying shells during the entire course of the evolution. The correlation function has a finite width () in panels b, c, and d such that its value approaches zero for . The correlation function approaches a symmetrical shape as we traverse the O and Ne shells. If the correlation function is evaluated at the centre of a convective region then it will have a symmetrical shape, which will change as we go to the bottom/top of the region. As one moves from the Ne shell into Zone II the function again becomes asymmetric, which implies that after there is a well defined convective region centered somewhere between . As we traverse Zone III and go into the He shell, the correlation length becomes large () and the function approaches a symmetrical form. This suggests that after there is a well defined convection layer centered somewhere between . The final state of the simulated shell has two well defined convective regions with a contact somewhere inside Zone II.

4.3 Thermodynamic Evolution

The continuous deposition of thermal energy by nuclear burning powers convection. In the following sections we compare nuclear energy production and entropy evolution in the 1D and 3D models.

4.3.1 Nuclear Energy Generation: 1D vs. 3D


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefLeft panel: Space-time plot showing net energy generation ( ,color-coded) rate for the 1D model. The solid curves (i) show initial mass fraction profiles and the dashed curves (f) show the final mass fraction profile for O (black), Ne (red) and Si (blue) with the corresponding scale on the top of the panel. Burning happens at the base of the O shell close to . Neutrino cooling is dominant in the central part of the Fe/Si core (gray). The burning rate is practically constant in time. The right panel shows the color-coded net energy generation rate for the 3D model. The upper panel shows the region inside . The energy generation is relatively constant in time except some changes close to the end resulting from Ne depletion. The lower panel shows a zoom of the region inside . The energy production here is dominated by O burning at base of the O shell. Material entrained from the Ne shell burns on its way down forming a secondary burning layer separated from the O burning zone. At there is a rapid increase in the spatial extent and energy production rate inside the secondary layer. This marks the episodic burning phase; a significant amount of Ne entering the O shell is ignited and consumed over a short duration. The colored lines represent the Ne mass fraction at different times(scale at the top of the panel).


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSequence of 3D volume renderings following the shell merger in progress. The Ne mass fraction is volume rendered and appears as nearly a transparent “cloud” (with the opacity function shown in the lower right corner of each panel). The inner opaque bubble with a convoluted morphology represents a color plot of the Si mass fraction (the color coding of latter is given in the upper left corner of each panel) on a radial velocity isosurface (at ). The choice of radial velocity for the isosurface is motivated by tracking the average value of Si mass fraction. The box is across, which is the extent of the main Ne shell (see Figure 2). The progression of states (from left to right and top to bottom) shows the violent O-Ne shell merger in progress.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefUpper panel: Specific kinetic energies in the radial velocity fluctuations as functions of time for the different layers. The evolution can be split into four phases: a) onset of convection, b) stationary phase, c) rapid rise, and d) saturation phase. The entire region between the inert Si shell and the Ne shell attains similar specific kinetic energy after the shell merger. Lower panel: Ratio of kinetic energy stored in the radial fluctuations () and the transverse flow () as a function of time. The horizontal line marks the value (1/2) when the kinetic energy is equally partitioned between all the components as in case of isotropic turbulence, showing that the flow is primarily dominated by large-scale radial motions due to convection.

Figure 4.3.1 (left panel) shows the net energy generation rate in the 1D model as a function of time. In the region inside (at higher density and temperature inside the Fe core) the neutrino cooling rate dominates over the nuclear energy generation rate leading to core contraction. In the Si layer energy deposition rate by nuclear burning of Si to Fe exceeds the energy loss rate due to neutrino cooling. Most of the nuclear energy generation happens at the base of the O shell between . In addition, there are contributions from Ne shell burning, C shell burning, and He shell burning. The energy production rate in the 1D model in all of these burning layers is relatively constant in time.

Figure 4.3.1 (right panel) shows the net energy generation rate ( everywhere) in the 3D model. The evolution between (upper half) does not show a lot of deviation from its initial behaviour, except that close to the Ne burning shell moves inwards in mass (in a Lagrangian sense). The lower half shows a zoom of the region between . Most of the energy production happens at the base of the O shell. Around another burning layer develops above the O burning zone but below the Ne shell. Ne entrained from the Ne shell (see colored lines representing the Ne mass fraction at different time) burns on its way to the base of the O shell. Both the extent of this layer and intensity of burning increases with time, and around there is a short but intense Ne burning episode. This phenomenon is marked as episodic burning, which is reminiscent of hot-spot burning seen by bazan_94 in their 2D study of oxygen burning in a progenitor (arnett_94).

Figure 4.3.1 shows the resulting shell merger marked by Si rich outflow rapidly moving into the Ne shell. The isosurfaces shown track the points corresponding to radial velocity. The Si-rich material moves rapidly outwards carving out an elongated cavity in the enclosing Ne shell. The velocity isosurface expands from to more than in a short span of . The episodic nuclear burning peaks at around (as shown in Figure 4.3.1 by contours equally spaced between ) and triggers/powers a phase of violent convective activity. A comparison between the Ne abundance at and shows that the total amount of Ne decreases by because of the episodic burning.

A rough estimate of the energy produced can be obtained by considered the principle energy producing reactions in Ne burning which effectively convert two nuclei to and nuclei (iliadis_2015, Equation (5.108))


with the average energy production (per nucleons) of () near (iliadis_2015, Equation (5.109)). This implies a total release of of thermal energy inside of stellar plasma within , which further translates into in each component (, and ) of kinetic energy.

The upper panel of Figure 4.3.1 shows the growth of the specific kinetic energy (per unit mass) in the fluctuating component in all of the layers, individually defined as


where are the inner boundary, outer boundary and mass, respectively, of the th layer (defined in Section 4.2.2). Kinetic energy starts building up with the onset of convection and reaches a plateau phase by except for the Ne shell. The stationary phase lasts till , when episodic Ne burning (the shell merger) begins. In the next there is a rapid rise in kinetic energy by a factor of . In the end the growth saturates, and the entire merged shell has uniform specific kinetic energy. The lower panel of Figure 4.3.1 shows the ratio of kinetic energy in the fluctuating radial component and the transverse components. The transverse components of the specific kinetic energy are defined as


The ratio fluctuates with time and hardly stays close to the equipartition value () for the convectively active shells. For the stable Zone I the value deviates away from when the shell is finally eroded away. Interestingly, transverse motions dominate over radial motion in Zone II after .

4.3.2 Entropy: 1D vs 3D


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefEntropy evolution in the 1D model (left panel) and the 3D model (right panel). The dotted blue curve represents the initial () entropy profile in each of the panels. In the 1D model, most of the change is confined to the region between . Energy deposited by O burning causes a gradual increase in entropy until . At the end the sharp jump in entropy at is softened slightly. In the 3D model efficient heat transport quickly levels the entropy gradient inside . In the last rapid energy deposition by episodic Ne burning leads to a negative entropy gradient extending from base of the O shell to the base of Zone II. Note: The dashed horizontal line in the right panel is to guide the eye to the inversion of entropy gradient, which happens between and .

Layers inside a star tend to become convective when radiation and neutrinos are unable to carry away the thermal energy deposited by nuclear burning. In a stably stratified region the entropy gradient is positive (entropy increases radially outwards). Figure 4.3.2 shows the entropy profiles for the 1D model (left panel) and the 3D model (right panel). The differences between the entropy profiles of the 1D and 3D models are conspicuous. The negative entropy gradient at the base of the O layer facilitates convection.

In case of the 1D model, the entropy in the region inside increases gradually until the first (), keeping the gradient unchanged. During the last minute the jump in the entropy profile close to is reduced, making this region more prone to overshooting. Accelerated oxygen burning cased by Fe/Si core contraction leads to a rapid increase in entropy at the base of the O layer, which further steepens the entropy gradient.

In case of the 3D model the entropy evolution is remarkably different. The entropy profile inside is levelled in the first , but the additional Ne burning preserves the small negative entropy gradient between and . At the same time the steep entropy barrier close to is softened. With the complete disappearance of the stabilizing gradient provided by Zone I at , thermal energy is efficiently convected into a larger volume. However, convection is unable to carry away the excess thermal energy deposited by episodic Ne burning. This leads to a rapid rise in entropy at the base of the O shell and an overall negative entropy gradient extending up to , which is the base of Zone II. In the end the 3D model has higher entropy compared to the 1D model.

4.4 Mixing and Shell Merger


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefSlices ( plane) showing the spatial distribution of Ne (upper row) and Si (lower) at various times compare to the volume-rendered 3D views in Figure 4.3.1. The inner cavity in the Ne distribution (panels a and b) is shaped by the expanding Si shell (panels d and e). Panel b shows prominent Ne downdrafts extending deep into the underlying region and panel e shows the large-scale Si-rich dipolar plumes. Panels c and f show the well merged/mixed Ne and Si abundances prior to the onset of gravitational collapse. Each panel shows data only inside a spherical region bounded by a convectively stable layer (Zone II at , and Zone III for ) for clarity of presentation.

Figure 4.4 shows the spatial distribution of Ne (upper row) and Si (lower row) confirming the picture of large-scale mixing in the 3D model. The spatial distributions of both Ne and Si at (panels at ) show significant deviations from spherical symmetry. The dipolar feature (panel at ) in the Ne distribution is a cavity created by a Si-rich bubble (panel at ). The downdrafts of Ne-rich material propagating into the Si-enriched inner volume (panel at ) provide the nuclear fuel for hot-spot burning. As Ne mixes inwards and burns, the Si-rich matter expands outwards resulting in a complete merger of the Ne shell and the Si-enriched oxygen shell as seen in panels at . The Ne mass fraction decreases considerably ( of Ne is burnt in total) after the merger and the leftover Ne is thoroughly mixed with Si.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefMass fraction profiles for O, Ne, and Si (from left to right) at equally spaced epochs for the 1D (upper row) and 3D (lower row) models. The O abundance in the 3D model (panel d) shows an increase inside , whereas such a feature is absent in the 1D model (panel a). Panel e shows the sudden decrease in the Ne abundance () in the last caused by rapid burning. On the other hand, the Ne abundance in the 1D model (panel b) is virtually unaffected. Panel f shows the outer boundary of the Si-enriched shell reaching out to in the 3D model over the course of simulation. In contrast, the 1D model (panel c) fails to capture the large-scale mixing of Si entirely. Note: The initial abundance profiles are shown in Figure 2 (top panel).


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefPower in various multipoles () of the density fluctuations (, left panel) and the radial velocity fluctuations (, right panel) at the end of simulation. Red marks indicate the mode number with the maximum amplitude at a given mass coordinate. The histogram at the bottom of each panel shows the probability density of multipoles with maximum power in the entire simulated shell. See text for more details.

Figure 4.4 compares the mass fractions of Si, Ne, and O in the 1D and 3D models. Although both 1D and 3D models show smoothing of the sharp features, the abundances in the 1D model are practically unaffected outside of . Abundance changes in the 1D model take place during the last , while in the 3D model they come about slowly over the course of the simulation. In the 1D model these changes results from convective mixing (MLT), overshooting, thermohaline mixing and semiconvection (treated as diffusive processes, for details refer to mueller_16c). In the 3D model they result from mass entrainment at interfaces separating convectively stable and unstable regions, the large-scale convective flow, and the rising Si-rich bubbles powered by rapid Ne burning.

Panels a and d compare the O abundance profiles. In the 3D model, the O mass fraction inside increases during the last . The increase results from a combination of large-scale mixing (convective downdrafts carrying O-rich material inwards) and rapid Ne burning (producing O via ). The 1D model fails to capture the evolution of the O abundance altogether. Panels b and e compare the Ne abundance profiles. The outer part of the neon-carbon shell located between remains unaffected in both cases. In the 3D model, the gradual rise in the Ne mass fraction inside continues till only to be reversed by a phase of rapid Ne burning (described in Section 4.3.1). In contrast, Ne is virtually unaffected in the 1D model. Panels c and f compare the Si abundance profiles. In the 3D model of Si is gradually transported out from the region between into the region between . As a result the outer boundary of the Si-enriched layers from to (close to the C/He interface).

5 Pre-Supernova Model Properties

Large-scale modes with are most effective for shock revival (mueller_15a; abdikamalov_16; mueller_16a). In order to gain a quantitative understanding of the modes, we have done spherical harmonic decomposition of the density and velocity perturbations. Figure 4.4 shows the multipole power distribution at collapse for the density fluctuations (left panel) and the velocity fluctuations (right panel). The red markers indicate the multipole with the largest power () at a given mass coordinate. The value of changes with radius, but there are continuous sections where it remains stationary. In the case of density fluctuations inside the merged O-Ne shell () a mixture of dipole () and quadrupole mode () dominates, and in the region outside of Zone II the dipole mode dominates (similar to the global non-spherical oscillation described by herwig_14). This confirms our visual impression obtained from Figure 4.2.1. In the case of radial velocity fluctuations a mixture of dipole mode and mode dominates in the merged O-Ne shell, whereas the dipole mode dominates in the C and He shells. The normalized probability distribution of (shown as histograms) confirms that dominate both the density fluctuations and the velocity fluctuations. Figure 5, left panels, show another view of the spectra at the middle points of the layers defined in Section 4.2.2. In the lower half-panel the spectra for the region inside (layers c and d) are shallow compared to the region outside and including Zone II. In the upper half-panel the spectra for are shallower in the merged O-Ne shell including Zone II. The spectra above agree quite well with the ‘Kolmogorov slope’.


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefLeft panels: Multipole power spectra () for density fluctuations (lower half-panel) and radial velocity fluctuations (upper half-panel) at the center of the various layers (c, d, e, f, g, and h) defined in Section 4.2.2 at the end of simulation. Black circles mark the computed harmonics (integer values only) and the smooth curves are cubic splines. The black lines in the upper half-panel indicate a slope of expected for the inertial range in a fully developed turbulent flow according to Kolmogorov’s theory. Right panel: Root mean square Mach number of the radial velocity fluctuations. The Si/Ne interface is at the boundary of layer c and layer d, and the O/He interface is at boundary of layer g and layer h. Up to about , distinct convective O and Ne shells can be recognized, both with strongly subsonic convection. After the shell merger, the convective Mach numbers increase considerably.

Figure 5 (right panel) shows the profiles of the turbulent Mach number corresponding to the radial velocity fluctuations at multiple times. The quantity is defined following mueller_16c as


where is the adiabatic sound speed. The Mach number in the O and Ne shells increases gradually until the shell merger. During the first , the value of is almost negligible in the Ne shell outside , which is due to the smaller Brunt-Väisälä in the Ne shell. The region outside stays quiescent until as the stable Zone II located between does not allow convective plumes to penetrate further. After the merger, between , increases rapidly by a factor of reaching a peak value of at the onset of collapse. Convection is so violent that the plumes penetrating deep into Zone II excite strong interfacial waves that reach as far as the He shell, visible by the non-negligible Mach numbers outside .


figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefPseudocolor maps of radial velocity on Si mass fraction isosurfaces for (left), (middle) and (right), which effectively probe the Si distribution at successively greater depths. The snapshots are taken at when the Fe-core begins to undergo gravitational collapse. “Scale” refers to the size of the -axis (vertically upwards). The left panel shows the large-scale Si-rich bubbles. The middle panel shows the asymmetry and clumpy Si distribution prior to collapse, and the right panel shows Si-rich plumes and downdrafts in the deep interior moving outward or inward, respectively, at .

Figure 5 shows pseudocolor maps of the radial velocity on different Si mass-fraction isosurfaces prior to collapse. As we go from small to large values of , we effectively probe the Si distribution in successively deeper regions of the star. The left panel shows large-scale Si-rich bubbles rising at moderate speeds of somewhere around . In the panel for (middle), we see a highly asymmetrical and clumpy distribution of Si with some of the clumps moving outwards at velocities in excess of . This may have an impact on the observed asymmetries in the element distributions in some young supernova remnants such as Cassiopeia A. In the panel for (right) we probe an even higher Si mass fraction which is present deep inside the O shell (scale of ), where we again see fast outward moving Si-rich plumes. In short, the Si distribution in bulk is clumpy and has asymmetric velocities (although it looks as if it had isotropized in panel f of Figure 4.4).

6 Conclusions

In this work, we present the first 3D simulation (full solid angle) of a fully developed large-scale O-Ne shell merger prior to core-collapse in an progenitor. The work builds upon the previous study of mueller_16c in which they simulated the last minutes of O-shell burning in an supernova progenitor up to the onset of core collapse. In the present case we find significant differences between the kinematic, thermodynamic, and chemical evolution of the 3D and the 1D models. Ne is mixed deep into the O-shell, leading to a fuel ingestion episode which releases significant amounts of thermal energy by nuclear burning on a short timescale (), further boosting convection: Convective downdrafts transport Ne close to the base of the O shell, where it ignites and powers convection in return, leading to a positive feedback loop for rapid Ne entrainment. The entire simulated shell attains a convective Mach number of . The maximum convective velocities are for updrafts and for downdrafts. In contrast the 1D model is relatively quiescent with convective shell burning constrained to layers. The specific kinetic energy in the merged shell increases by a factor of during the course of evolution. Although the Si-rich material forms a strongly dipolar structure in the merged shell at , the distribution of elements becomes more isotropic in the following time leading up to collapse. However, the velocity field still shows a large global asymmetry at collapse. In short, the 3D pre-supernova model exhibits strong density perturbations and large scale velocity asymmetries (with dominant modes) in the flow as well as the distribution of nuclear species.

Recent simulations (couch_15; mueller_16c) have already shown that the progenitor structure is genuinely three-dimensional at collapse, but only considered quasi-steady state convection. Our work demonstrates that asymmetries can become exceptionally large when the convective flow becomes non-steady because of a shell merger. In such a case, the nuclear timescale becomes comparable to the convective timescale, and hence nuclear burning becomes strongly coupled with the flow dynamics and the resultant mixing is considerably faster and can be highly asymmetric. 1D stellar evolution calculations are severely limited in capturing these phenomena self-consistently as they occur primarily due to a combination of turbulent entrainment, fuel ingestion, and the excitation/propagation of internal waves, which are inherently 3D effects. Despite the inherent limitations of 1D calculations, such shell mergers may still take place in 1D (as seen by rauscher_02 and in recent studies by sukhbold_14 and davis_2018, and merely be delayed as in model S25 of rauscher_02, where the shell merger happens only before collapse. However, shell mergers in 1D models cannot capture the highly dynamical and asymmetric flow during and after the merger. Therefore, the dynamics of shell mergers during the late burning stages can only be captured using self-consistent multi-dimensional calculations. Such calculations are still in their infancy and the work reported here attempts to make progress in our understanding of the internal structure of supernova progenitors.

Also, recent studies show that progenitor asymmetries may help in shock revival (couch_13; mueller_15a; mueller_17) by aiding the growth of convection and/or instabilities like the standing accretion shock instability (SASI) (takahashi_16; mueller_17). Our work furnishes further evidence that some supernova progenitors have convective seed perturbation that are dynamically important. For a shell merger like the one described in this work, we see even more violent convection and stronger large-scale asymmetries than in the model of mueller_17.

Shell mergers will also affect the compactness parameter (oconnor_11) and thus the explodability with possible consequences for the location and extent of “islands of explodability” (e.g. sukhbold_14). Due to the high convective Mach numbers and pronounced global asymmetries, we expect a strong impact on the explosion dynamics for our model.

This work constitutes only the first step towards understanding the importance of late-stage shell mergers. For example, mergers may also impact the pre-supernova nucleosynthesis yields, as already suggested by 1D stellar evolution models. The S20 model of rauscher_02 showed a strong overproduction of several elements between Si and V and an underproduction of Cr, Mn, and the light Fe isotopes, which they ascribed to the stellar structure resulting from the C/O shell merger. ritter_18 showed that even moderate C-ingestion events during O burning can significantly overproduce odd-Z elements such as K, Sc, Cl, etc. Some of their models with full C-O shell merger (albeit artificially increased O burning luminosity) have overproduction factors greater than 10. Therefore, we expect to see substantial effects on the nucleosynthetic yields due to the atypical way in which Ne burns during the shell merger event described in our work. During the explosion phase, the substantial asymmetries in the element distribution and the velocity field may have repercussions for observable asymmetries in supernovae and nucleosynthesis (hughes_2000; delaney_10; lopez_2018). Such questions will be addressed in the future.


This project was supported by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA and by the Deutsche Forschungsgemeinschaft through the Excellence Cluster “Universe” (EXC-153). This work was also supported by the Australian Research Council through ARC Future Fellowships FT160100035 (BM), Future Fellowship FT120100363 (AH). This research was undertaken with the assistance of resources obtained via NCMAS and ASTAC from the National Computational Infrastructure (NCI), which is supported by the Australian Government and was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. On the Garching side, computing resources from the Gauss Centre for Supercomputing (at the Leibniz Supercomputing Centre (LRZ) under SuperMUC project ID: pr53yi) are acknowledged. The analysis of the simulation was done using the MPG Supercomputer Hydra, Cobra, and Draco provided by The Max Planck Computing and Data Facility at Garching. Software: NumPy & Scipy (oliphat_2007), IPython (PER-GRA:2007), Matplotlib (Hunter:2007), Visit (Childs_11), and SHTns  library (schaeffer_2013). This research has made use of NASA’s Astrophysics Data System. \onecolumngrid APPENDIX

A Angle Averaging Method

An intensive physical quantity (for e.g., density) can be written as a sum of its average and the fluctuating component (notation for mean and fluctuating components follows from Oertel_2010) as , where the angle-averaged value (Reynolds) is defined as


An extensive physical quantity (for e.g., specific energy ) can be written as , where the angle-averaged value (Favre) is defined as


The fluctuating components and satisfy the following relations


B Spherical Harmonic Transforms

The forward spherical harmonic transform of a function involves computing the coefficients according to


where is the orthonormalized spherical harmonic function of degree and order . The total power for a mode with multipole order is


We have used the routines from the SHTns library (schaeffer_2013) to calculate the spherical harmonic decomposition. \bibliography@latexpaper.bbl

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description