Large-Scale Mixing in a Violent Oxygen-Neon Shell Merger Prior to a Core-Collapse Supernova
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.
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
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
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
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
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.
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
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 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 .
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
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:
, convectively inert Si shell,
, convective O shell,
, Zone I,
, convective Ne shell,
, Zone II,
, convective C shell,
, Zone III,
, 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 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
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 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 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 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 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).
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. \email@example.com