Large-Scale Mixing in a Violent Oxygen-Neon Shell Merger Prior to a Core-Collapse Supernova
Naveen Yadav, Bernhard Müller, Hans Thomas Janka, Tobias Melson, Alexander Heger
DD RAFT VERSION A PRIL
14, 2020Typeset using L A TEX twocolumn style in AASTeX62
Large-Scale Mixing in a Violent Oxygen-Neon Shell Merger Prior to a Core-Collapse Supernova N AVEEN Y ADAV ,
1, 2 B ERNHARD
M ¨
ULLER ,
3, 4 H ANS T HOMAS J ANKA , T OBIAS M ELSON , AND A LEXANDER H EGER
3, 5, 61
Max Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany Exzellenzcluster ORIGINS, Boltzmannstr. 2, D-85748 Garching, Germany Monash Centre for Astrophysics, School of Physics and Astronomy, Monash University, Victoria-3800, Australia Astrophysics Research Centre, School of Mathematics and Physics, Queen’s University Belfast, Belfast, BT7 1NN, UK Center for Nuclear Astrophysics, Department of Physics and Astronomy, Shanghai Jiao-Tong University, Shanghai-200240, P. R. China Joint Institute for Nuclear Astrophysics, 1 Cyclotron Laboratory, National Superconducting Cyclotron Laboratory, Michigan State University, East Lansing,MI-48824-1321, USA (Received 00, 00, 2019; Revised 00, 00, 2019; Accepted 00, 00, 2019)
Submitted to ApJABSTRACTWe present a seven-minute long π -3D simulation of a shell merger event in a non-rotating . M (cid:12) su-pernova progenitor before the onset of gravitational collapse. The key motivation is to capture the large-scalemixing 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 largestpossible scales. We find significant differences between the kinematic, thermodynamic and chemical evolutionof the 3D and the 1D model. The 3D model shows vigorous convection leading to more efficient mixing ofnuclear species. In the 3D case the entire oxygen shell attains convective Mach numbers of ≈ . , whereas inthe 1D model, the convective velocities are much lower and there is negligible overshooting across convectiveboundaries. In the 3D case, the convective eddies entrain nuclear species from the neon (and carbon) layersinto the deeper part of the oxygen burning shell, where they burn and power a violent convection phase withoutflows. This is a prototypical model of a convective-reactive system. Due to the strong convection and theresulting efficient mixing, the interface between the neon layer and the silicon-enriched oxygen layer disappearsduring the evolution, and silicon is mixed far out into merged oxygen/neon shell. Neon entrained inwards byconvective downdrafts burns, resulting in lower neon mass in the 3D model compared to the 1D model at timeof collapse. In addition, the 3D model develops remarkable large-scale, large-amplitude asymmetries, whichmay have important implications for the impending gravitational collapse and the subsequent explosion. Keywords: stars:massive – convection – hydrodynamics – turbulence – supernovae: general INTRODUCTIONMulti-dimensional effects in the late burning stages ofmassive stars have recently garnered considerable interestsfor a number of reasons. Whereas diffusive mixing mod-els based on the mixing-length theory (MLT) of convec-tion (Biermann 1932; B¨ohm-Vitense 1958) as used in one-dimensional (1D) stellar evolution models provide a reason-able estimate of mixing (Arnett et al. 2019, for a more quan-titative discussion) in the interior of convective zones for thelate, neutrino-cooled burning stages (M¨uller et al. 2016b, seeFigure 10 & 14) , it has long been speculated that additional Corresponding author: Naveen [email protected] The upper panel of Figure 10 shows a comparison between the radialrms velocity fluctuations δv r (black) in the 3D model to the convective ve-locity v conv computed in the Kepler model using MLT (red). Inside theconvection zone ( . − .
25 M (cid:12) ) the 1D and 3D velocities are in close phenomena such as turbulent entrainment (Fernando 1991;Strang & Fernando 2001; Meakin & Arnett 2007a; Spruit2015) and the excitation of internal waves could significantlyaffect shell growth, mixing, and angular momentum trans-port (Cantiello et al. 2014; Fuller et al. 2015) in a mannerthat is not captured by current 1D stellar evolution models.Convective-reactive systems (Dimotakis 2005) are at the ex-treme end of mixing problems, where the feedback producedby nuclear burning strongly affects the flow and vice-versa,thereby making it extremely hard to describe chemical mix-ing in these systems using simplified 1D prescriptions.Since the early 1990s, various groups have attemptedto study the late convective burning stages using multi- agreement. Figure 14 shows a comparison between the mass fraction pro-files in the 3D model and the 1D model. Although, the differences in themass fractions in the interior of the oxygen shell are minute, the differencesin gradients at shell boundaries are conspicuous. a r X i v : . [ a s t r o - ph . H E ] D ec Y ADAV ET AL .dimensional simulations to investigate such effects. Sincethe seminal early work in two dimensions (2D; Bazan & Ar-nett 1994, 1998; Asida & Arnett 2000) and three dimensions(3D; Kuhlen et al. 2003; Meakin & Arnett 2006, 2007a,b),additional convective boundary mixing has indeed been con-sistently observed in many modern simulations (M¨uller et al.2016b; Cristini et al. 2017; Jones et al. 2017; Andrassy et al.2018) and appears to be well captured by the semi-empiricalentrainment laws familiar from terrestrial settings. The long-term impact of such extra mixing on the evolution of massivestars remains more unclear, however. One major caveat con-cerns the duration of the simulations, which are currentlylimited to periods far shorter than the thermal adjustmenttime scale; and it has also been argued that the predictedextrapolation of entrainment rates might in many cases notqualitatively alter stellar structure over secular time scales(M¨uller 2016).In some highly dynamical situations, however, multi-dimensional effects may result in qualitatively different be-haviour compared to 1D stellar evolution models. Suchsituations often occur when material entrained across shellboundaries burns violently, which can lead to strong feed-back on the dynamics of the flow; proton ingestion episodesas studied in 3D by several groups (Herwig et al. 2011; Stan-cliffe et al. 2011; Herwig et al. 2014; Woodward et al. 2015)are one noteworthy example. Such violent episodes of con-vective boundary mixing are also interesting because thereis actually a potential fingerprint for the multi-dimensionaldynamics in the stellar interior in the form of the nucleosyn-thesis enabled by turbulent entrainment (Herwig et al. 2011;Jones et al. 2016; Ritter et al. 2018).Such dynamical mixing events may also occur during thelate stages of burning in massive stars. 1D stellar evolutionmodels already show that mergers between O, Ne, and Cshells are quite commonplace (Sukhbold & Woosley 2014,see, e.g., the M (cid:12) model in their Figure 8); they are fa-cilitated by relatively small buoyancy barriers between theseshells. What is particularly intriguing is that these mergersvery often occur only a few turnover time-scales before col-lapse, as pointed out by the systematic study of Collins et al.(2018), who found such late mergers in about percent ofall progenitors between M (cid:12) and M (cid:12) . This implies thatthe merged shells may be caught in a highly dynamical stateat the time of the supernova. It is therefore conceivable thattraces of such a merger remain imprinted in the explosiongeometry. If this is the case, such mergers could help ex-plain asymmetric features in supernova remnants such as thebroad Si-Mg rich “jet”–like features in Cassiopeia A (De-Laney et al. 2010; Isensee et al. 2010; Grefenstette et al.2014, 2017). 1D stellar evolution models with parameter-ized mixing also suggest that such shell mergers could resultin very characteristic nucleosynthesis.The recent surge of interest in the 3D dynamics of the lateburning stages has mostly focused on the potentially bene-ficial role of convective seed perturbations in the supernovamechanism (M¨uller & Janka 2015; M¨uller et al. 2017; Couchet al. 2015) and on slow, steady-state entrainment (Meakin & Arnett 2007a,b; M¨uller et al. 2016b; Jones et al. 2017;Cristini et al. 2017; Andrassy et al. 2018), but such dynami-cal shell mergers have not yet been investigated thoroughly.M¨uller (2016) reported a late breakout of convection acrossshell boundaries in a thin O burning shell that was reignitedshortly before collapse in an . M (cid:12) progenitor, but withinsufficient time before collapse for a genuine merger of theO and Ne shell. More recently, Moc´ak et al. (2018) observedthe merging of shells in a 3D simulation of a M (cid:12) star, buttheir model was restricted to a small wedge of . ◦ × . ◦ and not evolved until core-collapse. Moreover, the mergeroccurred already within the initial numerical transient beforea 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 sili-con and neon shell up to the onset of core-collapse. Our sim-ulation follows the last seven minutes of convective O andNe shell burning in an . M (cid:12) progenitor (M¨uller et al.2016a) from an early stage with clearly separated convec-tive shells to the point where the buoyancy jump between theO and Ne shells is reduced sufficiently to trigger a violentmerger that is still not fully completed at collapse. In thispaper, we focus on analyzing the flow dynamics and mixingduring the merger and the comparison of the 3D model to thecorresponding 1D stellar evolution calculation. Implicationsof the shell merger for the ensuing supernova, though a ma-jor motivation for the current simulation, will be discussed infuture work.The layout of the paper is as follows. The initial model,which has been derived from a 1D simulation is describedin Section 2. In Section 3, we provide important details ofthe numerical method and the microphysics used for the 3Dsimulation. In Section 4 we describe in detail the 3D modeland compare it with the 1D model. In Section 5 we brieflydiscuss the pre-supernova model. In Section 6 we provideour conclusions and discuss the impact of our results on thepre-supernova structure and the core collapse.
1D PROGENITOR MODELWe consider a progenitor with solar metallicity and a zero-age main sequence mass of . M (cid:12) . The 1D modelhas been calculated using the K EPLER stellar evolutioncode (Heger & Woosley 2010) as outlined in M¨uller et al.(2016b, their Section 2.1) . The 1D model is mapped tothe 3D grid s before the onset of core collapse. Fig-ure 1 shows the relevant properties of the initial 1D model.The top panel shows the mass fraction profiles for key nu-clei, the specific nuclear energy generation rate, and thespecific neutrino cooling rate. The shells of principal inter- K EPLER uses the Ledoux criterion for semiconvection as described inWeaver et al. (1978). We use α = 1 for MLT (larger [ ≈ . ] and depth-dependent MLT values used in the solar convection zone may only be validfor the surface convection zone of the sun). Convective boundaries are‘softened’ by one formal overshoot zone with / mixing efficiency assemiconvection (Weaver et al. 1978). We emphasize that we do not use the f -parameter model for overshoot in our 1D models, as the physics behind“overshooting prescriptions” is very uncertain, and there is no evidence orconsensus about their validity in all evolution stages. HELL MERGER IN
3D 3 . . . . . m a ss f r a c t i o n HeC ONe MgSi SAr
Si ONe C . . . v [ k m s ] . . . . . . m [ M ] ⇢ [ g c m ] . . . . . . r [10 km] . . . . ˙ q [ e r gg s ] ˙ q nuc ˙ q ⌫ vv conv s [ k b / nu c l e o n ] ⇢s v c o n v [ k m s ]
Properties of the initial 1D model. The top panel shows the mass fractions of some relevant α -elements, the specific nuclear energygeneration rate, and the specific neutrino cooling rate with | ˙ q nuc | (cid:29) | ˙ q ν | at the base of the O burning shell, whereas | ˙ q ν | (cid:29) | ˙ q nuc | in the Fe/Sicore. The middle panel shows the radial velocity and the convective velocity according to MLT. The Si/Fe core is slowly contracting whereasthe O shell is slowly expanding. The O burning shell has larger convective velocity compared to the Ne and C burning shells, and the Fe/Sicore is convectively inert. The bottom panel shows the density and entropy profiles. Entropy jumps correspond to composition interfaces; largerelative jumps in entropy means high stiffness of such boundaries against convective entrainment. Vertical black dashed lines mark the shellwhich is simulated in 3D. est in this paper are the O-rich shells between . − . M (cid:12) ( ≈ , − , km). These include a convective O burn-ing shell with ashes of Si and S, a convective Ne burningshell between . − . M (cid:12) , which is separated by a non-convective layer from a convective C burning shell with ashesof O and Ne between . − . M (cid:12) . The O and Ne shells areseparated by a thin interface at ≈ . M (cid:12) . The Si shell andthe Fe core inside . M (cid:12) are practically inert, with neutrinocooling dominating over burning. The middle panel showsthe profiles of the radial velocity, v , and the convective veloc-ity, v conv , calculated using the same choice of dimensionlesscoefficients as (M¨uller 2016). The convective velocity nearthe base of the O shell reaches up to ∼ km s − . The Neand C shells are relatively quiet, having convective velocitiesclose to km s − and km s − , respectively. The bottompanel displays the density and (specific) entropy profiles.The entropy profile shows the three aforementioned con-vective regions as flat sections. Significant entropy jumps of ≈ . k b / nucleon and ≈ . k b / nucleon (here k b is the Boltz-mann constant) close to the bottom (just above the Fe/Sicore) and close to the top (just below the He shell) of thethree shells make the convective boundaries quite stiff, withthe inner boundary of the oxygen shell even being visible as Following common practice in stellar evolution literature, we label theconvective shells by the fuel that burns at their base, unless we explicitlydiscuss the composition in more detail. a discontinuity in the density profile. By contrast, the entropyjump between the oxygen and neon shells is small, and theneon and carbon shells are separated by a stable region withgradually increasing entropy, but not by a discontinuity.To be clear, choosing “ . M (cid:12) ” (accurate to two deci-mal places) does not amount to “fine-tuning” of the modelpresented in M¨uller et al. (2016b). We would like to em-phasize that the stellar model presented in this work and thestellar model used by M¨uller et al. (2016b) are ≈ M (cid:12) apart,and it is well known that the pre-supernova structure of a staris a very sensitive function of the zero-age main sequencemass (ZAMS) (Sukhbold et al. 2018, see upper left panel oftheir Figure 8). Collins et al. (2018) evaluated convectivevelocities and eddy scales in the oxygen and silicon burningshells of a large set (2353 solar-metallicity, non-rotating, cal-culated using the K EPLER code) of progenitor models withZAMS between . and M (cid:12) . One of their key findingis that shell mergers are not rare; “about 40 percent of pro-genitors between 16 and 26 M (cid:12) exhibit simultaneous oxygenand neon burning in the same convection zone as a result of ashell merger shortly before collapse”. At the same time, thereis no clear mapping between the outcomes in 1D and 3D, andtherefore, it cannot be said a priori with certainty if a certain1D model will exhibit a shell merger in 3D simulation. Ina nutshell, the . M (cid:12) model was chosen for this study asthe dominant angular wave numbers ( (cid:96) ) of convection in Si Y ADAV ET AL . Table 1.
Setup details for the 3D simulation.Quantity Option/Value a Radial grid GeometricCells in radial direction, N r b Cells in polar direction, N θ b Cells in azimuthal direction, N φ c No. of ghost cells ( θ, φ ), N g c No. of buffer cells ( θ, φ ), N b ◦ d Inner boundary at t = 0 , r − . × km d Outer boundary at t = 0 , r + . × kmInner boundary condition, radial refer to Section 3.1Outer boundary condition, radial refer to Section 3.1Gravitational potential Newtonian, SphericalSimulation time s e Perturbation amplitude ( f ) × − Neutrino cooling YesEquation of state Helmholtz EoS f CFL 0.6 ar i +1 = (1 + α ) r i , where α = ( r + /r − ) /N r . b The number of angular cells are specified per patch of theYin-Yang grid. c For each grid patch (Yin and Yang). d The inner and outer boundaries are moved according to thetrajectory obtained from 1D model. e δρ = fρ U [ − , , where U [ − , is a uniformly dis-tributed random number in the interval [ − , . f CFL is the value of the Courant–Friedrichs–Lewy condition. and O shell were small, and at the same time the convectiveMach number was large. NUMERICAL METHODS: 3D MODELThe 3D simulation uses the P
ROMETHEUS code. The sim-ulation employs a moving (but non-Lagrangian) grid in theradial direction and an axis-free Yin-Yang grid (Kageyama &Sato 2004; Wongwathanarat et al. 2010) as implemented inMelson et al. (2015) with a uniform angular resolution ( ◦ ) inthe θ and φ directions . For the purpose of this simulation wechanged the th -order reconstruction originally used in the Adequacy of ◦ angular resolution: A simulation using an equidistantspherical grid ( ∆ θ = ∆ φ = α ) captures angular modes with (cid:96) ranging be-tween zero and ◦ /α . Therefore, in our case with α = 2 ◦ , we can extractmodes up to (cid:96) = 90 . The well developed inertial range (see left panel of Fig-ure 15 in Section 5) is a proof that the simulation is in the regime of implicitlarge eddy simulations (Boris et al. 1992, ILES), and the energy transfer be-tween different length scales is correctly modelled. Moreover, ignoring thesmallest scale does not affect the evolution, as the mixing (and the relatedhot-spot burning as described in Section 4.3) is caused by the large-scaleconvective flow. piecewise-parabolic method (PPM) of Colella & Woodward(1984) to a th -order extremum-preserving reconstruction(Colella & Sekora 2008; Sekora & Colella 2009). We use theHelmholtz equation of state (EoS) described in Timmes &Arnett (1999) and Timmes & Swesty (2000), which is ther-modynamically consistent to high accuracy. Nuclear burn-ing is treated using a 19-species α -network (Weaver et al.1978). An accurate modeling of silicon burning is avoided asits nuclear burning in the quasi-statistical equilibrium (QSE)regime involves tracking a very large number of nuclearspecies making it computationally expensive. Energy loss byneutrinos (Itoh et al. 1996) is included as a local sink term.3.1. Boundary Conditions
We excise the core inside . M (cid:12) as most of the nuclearenergy generation is limited to a small region at the base ofthe O shell (described in Section 2) and the Si/Fe core is rel-atively inert. The region outside of . M (cid:12) is also excludedfrom the simulation as it remains dynamically disconnectedfrom the interior during the short 3D simulation. In the 1Dmodel the Si/Fe core cools via neutrino emission and con-tracts. For consistency, the radial boundaries of the compu-tational domain are moved in step with the 1D model (move-ment 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 shellsas the “region of interest” layers, and 3) a small part of the Heshell. We emphasize that the steep entropy step at the Si/Ointerface is about . M (cid:12) away from the inner grid bound-ary (see Figure 1). This convectively stable interface ensuresthat the convective mass motions in the simulation volumeremain radially separated from the grid boundary. For PPMreconstruction, we impose reflective boundary conditions forvelocity at both the inner and outer radial boundary and ex-trapolate all thermodynamical quantities assuming adiabaticand hydrostatic stratification. We strictly enforce zero ad-vective fluxes across the boundaries before updating the con-served variables. 3.2. The 3D simulation is initialized using density, temperature,and mass fraction profiles taken from the 1D model. Hy-drodynamic and thermodynamic variables are mapped usingsecond order polynomial interpolation. Mass fractions aremapped using linear interpolation (to preserve the disconti-nuities at shell boundaries). The question of initial transients(mapping errors) is discussed in Section 4.5. In order to breakthe spherical symmetry we introduce small deviations fromhydrostatic equilibrium by perturbing the density field suchthat δρ = U [ − , f ρ , where U [ − , is a uniformly dis-tributed random number in the interval [ − , and the am-plitude f = 2 . × − .The Newtonian gravitational potential is computed in themonopole approximation; the appropriate contribution fromthe excised inner core is added. The details of the simula-tion setup are summarized in Table 1. We have not donea convergence study for the model presented in the paper. HELL MERGER IN
3D 5Please refer to the discussion on the effect of angular and ra-dial resolution on a simulation of oxygen burning in case ofthe M (cid:12) model by M¨uller et al. (2016b, see Appendix B).The simulations are converged in terms of the coupling be-tween nuclear energy generation rate and growth of turbulentmotions (M¨uller et al. 2016b, Figure 18); and dominant mul-tipoles (M¨uller et al. 2016b, (cid:96) -modes in Figure 19) emergingin the flow. SIMULATION RESULTS: 1D VS 3DIn the following sections we discuss the kinematic, thermo-dynamic, and chemical evolution of the 3D model and com-pare it to the 1D model. We start by describing the convectiveinstability of the simulated shell as a function of time.4.1.
Convective stability
The Ledoux criterion can be used to test the stability ofa mass element against convective overturn. The Brunt-V¨ais¨al¨a frequency is given by ω BV = (cid:18) ∂ ln ρ∂r − ∂ ln p∂r (cid:19) g, (1)where g = − dΦ / d r < is the gravitational acceleration and Γ = (cid:18) ∂ ln p∂ ln ρ (cid:19) s,X i , (2)is the adiabatic index. For unstable modes ω BV > is thegrowth rate; for stable modes ω BV is imaginary and | ω BV | isthe oscillation frequency. Following Buras et al. (2006), weredefine the Brunt-V¨ais¨al¨a frequency such that ω BV := sign( C L ) (cid:114) − gρ | C L | , (3) C L := (cid:18) ∂ρ∂s (cid:19) X i ,p d s d r + (cid:88) i (cid:18) ∂ρ∂X i (cid:19) s,p d X i d r , (4)which is more convenient for visualization purposes. Theconversion of Equation (1) to Equation (4) is exemplified inAppendix A of M¨uller et al. (2016b) . Note that our defi-nition changes only affects the sign convention, but not theabsolute value of ω BV . For unstable modes ( C L > ), ω BV isthe growth rate, and for stable modes ( C L < ), ω BV is thenegative of the oscillation frequency. Therefore, accordingto the sign convention adopted in this paper ω BV > corre-sponds to convective instability.Figure 2 (left panel) shows the spherically averaged Brunt-V¨ais¨al¨a frequency ω BV profiles as a function of time for the3D model. ω BV is defined as ω BV ( r ) := 14 π (cid:90) Ω ω BV ( r ) dΩ , (5) Heger et al. (2005) show an alternative approach (used in K
EPLER stel-lar evolution code) which can be used to extract the contribution of compo-sition to ω BV without explicitly tracking composition. Table 2.
Convectively stable zones according to theLedoux criterion at the beginning of the 3D simulation.Zone Location (inner edge) Width m [ M (cid:12) ] r [ km ] ∆ m [ M (cid:12) ] ∆ r [ km ] I 2.4 , , II a , , III 4.4 , , a Zone II in the 1D model shows a strong and narrow bar-rier ( < . M (cid:12) ) close to ≈ . M (cid:12) . which is different from ω BV calculated using spherically av-eraged ρ , p etc. The plot shows three convectively stable lay-ers (zones I − III at masses of . − . M (cid:12) , . − . M (cid:12) , and . − . M (cid:12) respectively) sandwiched between more volu-minous convectively unstable regions. The location ( m and r ) and width ( ∆ m and ∆ r ) of these zones at t = 0 are listedin Table 2 in order of increasing mass coordinate. The initialand final entropy profiles (green and red lines in left panelof Figure 2) show the disappearance of entropy jumps be-tween the convectively stable and unstable regions. Suchentropy steps represent barriers that the convective eddiescannot cross. As burning at the bottom of the O shell in-creases the entropy (see right panel of Figure 10), the stabi-lizing gradient in Zone I gradually becomes weaker, and dis-appears altogether after ≈ s, allowing the O and Ne shellto merge into one large convection zone extending between . − . M (cid:12) . The disappearance of Zone I is mirrored by thedisappearance of the entropy jump at . M (cid:12) . This increasein entropy towards collapse is characteristic of the late burn-ing stages prior to collapse, when neutrino cooling is too slowto balance thermal energy deposition due to rapid burning inthe contracting shells. Zone II also becomes narrower after ≈ s, and internal gravity waves can transport more en-ergy across it. Different from Zone I, this stable barrier is notcompletely eliminated prior to collapse, however. The rightpanel of Figure 2 shows ω BV for the 1D model. The convec-tively stable zones present in the 1D model are much thinnerand remain unchanged over the course of the simulation. Inthe last ≈ s there is a reduction in the strength of Zone I.The contraction of the Si core after the end of the last Si shellburning phase, and deleptonisation of the Fe core leading tocontraction, drives more burning (at the base of) the O shell.This raises the entropy and weakens the ‘buoyancy gap’ be-tween the two shells (in Zone I). At the same time, Zone IIand Zone III are totally undisturbed in the 1D model.4.2. Flow Dynamics
Growth of Density and Velocity Fluctuations
Consider the density fluctuations ρ (cid:48) which are defined as ρ (cid:48) ρ := ρ ( r ) − ρρ , (6) Y ADAV ET AL . t [10 s] . . . . . . m [ M ] -0.4-0.3-0.2-0.10.00.10.20.30.4 ! BV [s ]4 . . . . . . . s [k b / nucleon] t [10 s] . . . . . . m [ M ] -0.4-0.3-0.2-0.10.00.10.20.30.4 ! BV [s ]4 . . . . . . . s [k b / nucleon] Figure 2.
Left panel: Space-time plot showing the spherically averaged Brunt-V¨ais¨al¨a frequency in the 3D run. We use the definition ofBrunt-V¨ais¨al¨a frequency from Buras et al. (2006). According to this definition, regions with negative ω BV are convectively stable and regionswith positive ω BV are convectively unstable according to the Ledoux criterion. Initially the active O shell ( . − . M (cid:12) ), the active Ne shell( . − . M (cid:12) ), and the active C shell ( . − . M (cid:12) ) are clearly separated by stable layers with positive gradients in the initial entropy profile(green curve for entropy vs. mass see axis on top of panel). The barrier between the O and Ne shells at . M (cid:12) disappears after s due to theincreasing entropy in the O shell. The result is a large merged convective layer with a slightly negative entropy gradient between ≈ . − . M (cid:12) at collapse (red curve for entropy vs. mass see axis on top of panel). Right panel: Space-time plot showing the Brunt-V¨ais¨al¨a frequency for the1D model. The convectively stable interfaces in the 1D model are very narrow. The stability properties in most of the O shell remain unchangedduring the course of evolution, except the changes in convectively stability at the base of Zone I close to the end (sign reversal). where ρ is the average density evaluated over a sphericalshell. Figure 3 shows the density fluctuations on an x − y slice at different times. The panel at s shows the devel-opment of plumes with density fluctuations at a level of fewpercent. These plumes are contained in the region below thefirst convectively stable layer (Zone I) interior to ≈ kmas described in Section 4.1; we refer to these as “primaryplumes”. Primary plumes overshoot into the stably stratifiedlayer above and are decelerated, creating “hot spots” in thedensity fluctuations as seen in both panels at s and s.They also excite internal gravity waves which transport en-ergy across the stable zone and thus perturb the region di-rectly above it, creating “secondary” plumes. In due coursethe mass entrainment caused by interfacial shear instabilitiesdriven by convection scour material off the stable interface(Strang & Fernando 2001, see Figure 1 and 3 there). As aresult of the entrainment and mixing the stabilizing gradi-ent in Zone I ceases to exist around s (see Figure 2 leftpanel). This initiates the formation of a large convective re-gion extending from the base of the O shell to the base ofZone II. The eddies overshooting into Zone II are deceler-ated, creating new hot spots close to ≈ , km, At this point, global asymmetries develop in the flow. The panel at s in Figure 3 shows an asymmetric feature with a sizeof ≈ , km, which characterizes the phase of vigorousconvective activity. Close to the end of simulation (pan-els at s and s) the density fluctuations in the regionwithin , km reach a level of ≈ − and are again ofsmaller scale structure, while in the region between ≈ , and , km they are as large as ≈ − with a markeddipolar asymmetry.The cautious reader may have also noticed the sharp dis-continuities in all the panels, and the associated “waves” inthe panels from s onward. The reason is that the fluc-tuation field renders the convective boundaries sharp; theconvective flow rapidly decelerates due to buoyancy brak-ing over a short length scale ( ∼ . H P , where H P is thelocal pressure scale height). The flow deposits energy andmomentum and produces density fluctuations which markthe convective boundaries. The convective boundaries arestiff; or in other words have a large bulk Richardson Number(Turner 1979, see Chapter 10), and act like a “stretched mem-brane”. These boundaries reflect/transmit the internal wavesgenerated by turbulent convection, and penetrative convec- HELL MERGER IN
3D 7 . . . . . . . . . . y [ k m ]
35 s .
03 +0 .
05 -0.03-0.01+0.00+0.01+0.03 .
03 +0 .
02 -0.02-0.01+0.00+0.01+0.02 .
03 +0 .
03 -0.02-0.01+0.00+0.01+0.02 x [10 km] y [ k m ]
323 s .
05 +0 .
05 -0.03-0.01+0.00+0.01+0.03 x [10 km] .
11 +0 .
16 -0.05-0.03+0.00+0.03+0.05 x [10 km] .
19 +0 .
38 -0.16-0.08+0.00+0.08+0.16 ⇢ /⇢
Slices ( x − y plane) showing the density fluctuations ρ (cid:48) /ρ at different times. The physical size of the region displays ranges from , km at early time to , km at late time. The time sequence demonstrates the transition from initially well-separated convectiveshells with strong density fluctuations at the boundaries from overshooting to a merging of the O and Ne shells (panels at s and s), andeventually the C shell (panels at s and s). 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 above the topleft corner of each panel. Each panel shows data only inside a spherical region bounded by a convectively stable layer (Zone I at s, Zone IIfor , , s and Zone III for and s) for clarity of presentation. Please see the animation provided as a supplementary material(Movie A). tion. On comparing the panels at s and s, we see thatthe waves are absent outside , km, and are present inthe convectively active region only. This indicates that thewaves are not spurious and are certainly linked to the con-vective activity. Figure 4 provides another perspective on theflow dynamics by showing the minimum and maximum ofthe fractional density fluctuations on spherical mass shells asfunctions of mass coordinate. One clearly sees (right panel)how the density fluctuations at shell interfaces initially de-crease in magnitude and spread out towards larger m as en-trainment whittles down the stabilizing entropy gradients atthe shell interfaces. After the shells have merged, the den-sity fluctuations grow considerably, especially in the outerpart of the convective region. Figure 4 (left) shows that theinner convective boundary moves towards lower m by en-trainment; here the density fluctuations in the boundary layeractually become stronger with time as the overall violence ofconvection increases. We also point out that the density fluc-tuations never touch the inner grid boundary, which confirmsthe proper choice of its location at . M (cid:12) . We next con-sider the turbulent velocity fluctuations. Because the model isinitially non-rotating, we do not include any non-radial com-ponents in the mean flow and decompose the velocity field as (see Appendix A) v ( r ) = ˜ v r ( r ) e r + v (cid:48)(cid:48) r ( r ) e r + v θ e θ + v φ e φ , (7)where ˜ v r is the Favre-averaged radial velocity. The fluctuat-ing component v (cid:48)(cid:48) r of the radial velocity is therefore given by v (cid:48)(cid:48) r ( r ) := v r ( r ) − ˜ v r . (8)Figure 5 shows the radial velocity fluctuations correspond-ing to the snapshots shown in Figure 3. The panel at sshows well-developed plumes extending from the O burn-ing region at ≈ , km to the base of Zone I at , km.The panel at s shows the secondary plumes extendingfrom ≈ , km to ≈ , km. The primary and sec-ondary plumes are physically separated from each other bythe Zone I. Also, the secondary plumes have lower veloci-ties (by a factor of − ) compared to the primary plumes.These plumes are driven by the active Ne shell and developlater because of a lower Brunt-V¨ais¨al¨a frequency. The ve-locity difference is also apparent in panel at s which isthe stage when the stable Zone I has already disappeared.Thepanel at s shows emergence of a large-scale flow span-ning the entire region from base of the oxygen burning layerto the base of Zone II. At this stage, convection becomes veryvigorous and the magnitude of velocity fluctuations increases Y ADAV ET AL . .
75 1 .
80 1 . m [ M (cid:12) ] − . − . . . . ρ / ρ . . . . . . m [ M (cid:12) ] − . − . − . − . . . . . .
59 s118 s179 s238 s 295 s359 s420 s
Figure 4.
Minimum (dash-dotted line) and maximum (thick lines) values of density fluctuations ρ (cid:48) /ρ as functions of mass coordinate. Thearrows 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 Sishell travel inwards in mass coordinate (as a consequence of slow entrainment) with time and their amplitudes grow, but they never reach theinner grid boundary at . M (cid:12) . Blue arrows (right panel) mark the positions of Zone I and Zone II. from ≈ km s − to more than , km s − (panels at s and s) within the next s. The panel at sshows the presence of a large-scale dipolar asymmetry in theflow. Convective downdrafts in the innermost region reachvelocities in excess of , km s − aided by rapid contrac-tion of the Fe/Si core.4.2.2. Convection Length Scale
The initial seed perturbations are applied as random cell-by-cell variations and hence on a very short length scale. Inthe absence of a buoyancy barrier convective flow naturallygrows to the largest possible unimpeded length scale. Thelongitudinal two-point correlation function gives a qualitativemeasure of the size of the convective region. It is defined as(Equation 15 of M¨uller et al. 2016b) C ( r, δr ) := v (cid:48) r ( r, θ, φ ) v (cid:48) r ( r + δr, θ, φ ) (cid:110) v (cid:48) r ( r, θ, φ ) v (cid:48) r ( r + δr, θ, φ ) (cid:111) / , (9)where we have used the Reynolds-averaged velocities (seeAppendix A). The correlation length Λ corr at a radius r isdefined by integrating the correlation function Λ corr ( r ) := r + (cid:90) r − C ( r, δr ) d r (cid:48) , δr = r (cid:48) − r, (10)between r − and r + , which are usually taken to be −∞ and + ∞ respectively. Meakin & Arnett (2007b, Appendix B) de- fine these points such that C ( r, r ± − r ) ≥ . . In the presentcase, estimating the correlation length using the equationabove is non-trivial because of the presence of multiple con-vectively stable zones. For this analysis we divide the sim-ulation domain into multiple layers according to their initialconvective stability and the burning processes: a. . − . M (cid:12) , convectively inert Si shell, b. . − . M (cid:12) , convective O shell, c. . − . M (cid:12) , Zone I, d. . − . M (cid:12) , convective Ne shell, e. . − . M (cid:12) , Zone II, f. . − . M (cid:12) , convective C shell, g. . − . M (cid:12) , Zone III, h. . − . M (cid:12) , He shell.Figure 6 shows the correlation function evaluated at roughlythe middle points of these layers. We have marked theboundaries of these layers by vertical lines. In panel a thecorrelation function C ( r, δr ) drops rapidly from unity with δr , resulting in Λ corr < , km at all times ( Λ corr /r (cid:28) ).Thus the Si layer underneath the O burning zone stays non-convective and dynamically disconnected from the overlyingshells during the entire course of the evolution. The corre-lation function has a finite width ( Λ corr (cid:46) r ) in panels b, c,and d such that its value approaches zero for δr ≈ , km. HELL MERGER IN
3D 9 . . . . . . . . . . y [ k m ]
35 s+509
557 -497-248+0+248+497
529 -327-163+0+163+327
571 -443-222+0+222+443 x [10 km] y [ k m ]
323 s+741
798 -564-282+0+282+564 x [10 km] x [10 km] v r [km s ]
Slices in the x - y plane showing radial velocity fluctuations v (cid:48)(cid:48) r (in km s − ) at different times. The panels s till s show theinitially slow but steady build-up of the convective velocities. After the merger of the O and Ne shells, convection becomes considerablymore violent, and large-scale flow patterns develop (see the text for details). The green circle marks the inner boundary of the computationaldomain. Minimum and maximum values in the plane are specified in the top left corner of each panel. Each panel shows data only inside aspherical region bounded by a convectively stable layer (Zone I at s, Zone II for , , s and Zone III for and s) for clarityof presentation. Please see the animation provided as a supplementary material (Movie A). The correlation function approaches a symmetrical shape aswe traverse the O and Ne shells. If the correlation functionis evaluated at the centre of a convective region then it willhave a symmetrical shape, which will change as we go to thebottom/top of the region. As one moves from the Ne shellinto Zone II the function again becomes asymmetric, whichimplies that after ≈ s there is a well defined convectiveregion centered somewhere between . − . M (cid:12) . As wetraverse Zone III and go into the He shell, the correlationlength becomes large ( Λ corr ∼ r ) and the function approachesa symmetrical form. This suggests that after ≈ s there isa well defined convection layer centered somewhere between . − . M (cid:12) . The final state of the simulated shell has twowell defined convective regions with a contact somewhereinside Zone II.4.3. Thermodynamic Evolution
The continuous deposition of thermal energy by nuclearburning powers convection. In the following sections wecompare nuclear energy production and entropy evolution inthe 1D and 3D models.4.3.1.
Nuclear Energy Generation: 1D vs. 3D
Figure 7 (left panel) shows the net energy generation ratein the 1D model as a function of time. In the region inside . M (cid:12) (at higher density and temperature inside the Fe core) the neutrino cooling rate dominates over the nuclear energygeneration rate leading to core contraction. In the Si layerenergy deposition rate by nuclear burning of Si to Fe exceedsthe energy loss rate due to neutrino cooling. Most of thenuclear energy generation happens at the base of the O shellbetween . − . M (cid:12) . In addition, there are contributionsfrom Ne shell burning, C shell burning, and He shell burning(not shown in the figure). Note that the C shell burning andthe He shell burning are not relevant for the dynamics. Also,the total energy generation in the He shell is negligible. Noenergy from it enters or leaves the simulation domain, as itsmain purpose is to provide ‘boundary pressure’.Figure 7 (right panel) shows the net energy generation rate( | ˙ q nuc | (cid:29)| ˙ q ν | everywhere) in the 3D model. The evolutionbetween . − . M (cid:12) (upper half) does not show a lot of de-viation from its initial behaviour, except that close to ≈ sthe Ne burning shell moves inwards in mass (in a Lagrangiansense). The lower half shows a zoom of the region between . − . M (cid:12) . Most of the energy production happens at thebase of the O shell. Around s another burning layer de-velops above the O burning zone but below the Ne shell. Neentrained from the Ne shell (see colored lines representingthe Ne mass fraction at different time) burns on its way to thebase of the O shell. Both the extent of this layer and intensityof burning increases with time, and around s there is ashort but intense Ne burning episode. This phenomenon is0 Y ADAV ET AL .
101 convectively inert Si shell1 . M , . ⇥ km59 s118 s179 s 238 s295 s359 s
101 convective O shell2 . M , . ⇥ km
101 zone I2 . M , . ⇥ km
101 convective Ne shell3 . M , . ⇥ km C ( r , r ) zone II3 . M , . ⇥ km
101 convective C shell4 . M , . ⇥ km
101 zone III4 . M , . ⇥ km r [10 km]
101 He shell4 . M , . ⇥ km abcdefgh Figure 6.
Velocity correlation functions (defined by Equation (9)) near the centers of layers a − h (listed in Section 4.2.2) at various times. Themass 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 s is shown as filled grey curve in each panel. Placing the panels d and h together (aligning them along the greenlines) shows two well developed convective layers with a contact somewhere inside Zone II. See the text for details. marked as episodic burning , which is reminiscent of hot-spotburning seen by Bazan & Arnett (1994, see Fig. 3) in their2D study of O burning in a M (cid:12) progenitor (Arnett 1994). Note that the term “hot-spot burning” is being used here in a figurativesense. It is meant to describe a place of vigorous nuclear burning activity;large-scale mixing of Ne deep in to the O shell and its subsequent burning atthe higher ambient density and temperature. We think that our usage of theterm ”hot-spot burning” confirms with the phenomenon described by Bazan& Arnett (1994, see their Figure 3). Quoting Bazan & Arnett (1998), “
Sig-nificant mixing beyond convective boundaries determined by mixing-lengththeory brings fuel ( C ) into the convective region, causing hot spots ofnuclear burning .” Figure 8 shows the resulting shell merger marked by Sirich outflow rapidly moving into the Ne shell. The isosur-faces shown track the points corresponding to km s − radial velocity. The Si-rich material moves rapidly outwardscarving out an elongated cavity in the enclosing Ne shell.The velocity isosurface expands from ≈ , km to morethan , km in a short span of ≈ s. The episodic nu-clear burning peaks at around s (as shown in Figure 7 by log ˙ q contours) and triggers/powers a phase of violent con-vective activity. A comparison between the Ne abundance at s and s shows that the total amount of Ne decreases HELL MERGER IN
3D 11 . . . . . l og ˙ q n e t , ˙ q n e t > . . . . . . . . t [10 s] . . . . . m [ M ] episodic burning . . . . l og ˙ q n e t [ e r gg s ] , ˙ q n e t > .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 . mass fraction 295 s359 s420 s t [10 s] m [ M ] l og ˙ q n e t [ e r gg s ] , ˙ q n e t > ˙ q net [erg g s ] , ˙q net < . . . . . . mass fraction ifif i f Figure 7.
Left panel: Space-time plot showing net energy generation ( ˙ q nuc − ˙ q ν ,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 thecorresponding scale on the top of the panel. Burning happens at the base of the O shell close to . M (cid:12) . Neutrino cooling is dominant in thecentral part of the Fe/Si core (top colorbar). The burning rate is practically constant in time. The right panel shows the color-coded net energygeneration rate for the 3D model. The upper panel shows the region inside . − . M (cid:12) . The energy generation is relatively constant in timeexcept some changes close to the end resulting from Ne depletion. The lower panel shows a zoom of the region inside . − . M (cid:12) . The energyproduction here is dominated by O burning at base of the O shell (see the contour at log ˙ q net = 15 ) till ≈ s. Material entrained from theNe shell burns on its way down forming a secondary burning layer separated from the O burning zone. At ≈ s there is a rapid increase inthe spatial extent and energy production rate inside the secondary layer. This marks the episodic burning phase; a significant amount of Neentering the O shell is ignited and consumed over a short duration. The colored lines represent the Ne mass fraction at different times(scale atthe top of the panel). by ∆ M Ne ≈ . M (cid:12) ( ≈
70 percent) because of the episodicburning.A rough estimate of the energy produced can be obtainedby considered the principle energy producing reactions in Neburning which effectively convert two Ne nuclei to O and Mg nuclei (Iliadis 2015, Equation 5.108) Ne + Ne → O + Mg , (11)with the average energy production (per nucleons) of Q Ne ≈ . MeV ( ≡ . × erg g − ) near T ≈ . GK (Il-iadis 2015, Equation 5.109). This implies a total re-lease of Q Ne ∆ M Ne ≈ . × erg of thermal energy inside ≈ . M (cid:12) of stellar plasma within ≈ s, which furthertranslates into ≈ × erg g − in each component ( E r , E θ and E φ ) of kinetic energy.The upper panel of Figure 9 shows the growth of the spe-cific kinetic energy (per unit mass) in the fluctuating compo- nent in all of the layers, individually defined as E r ( i ) := 12 δm i m i, + (cid:90) m i, − ρv (cid:48)(cid:48) r d m, (12)where m i, − , m i, + , δm i are the inner boundary, outer bound-ary and mass, respectively, of the i th layer (defined in Sec-tion 4.2.2). Kinetic energy starts building up with the onset ofconvection and reaches a plateau phase by ≈ − s exceptfor the Ne shell. The stationary phase lasts till ≈ s, whenepisodic Ne burning begins. In the next ≈ − s there isa rapid rise in kinetic energy by a factor of ∼ . In the endthe growth saturates , and the entire merged shell has uni-form specific kinetic energy. The lower panel of Figure 9shows the ratio ψ of kinetic energy in the fluctuating radialcomponent and the transverse components. The transverse2 Y ADAV ET AL . X Ne X Si X Si X Si X Si X Si X Si X Si X Si X Si X Ne X Ne X Ne X Ne X Ne X Ne X Ne X Ne X Si κκκκκκκκκ
230 s
250 s
270 s
290 s
310 s
330 s
350 s
370 s
390 s
000 km
000 km
000 km
000 km
000 km
000 km
000 km
000 km
000 km
Sequence of 3D volume renderings following the shell merger in progress. The Ne mass fraction is volume rendered and appearsas nearly a transparent “cloud” (with the opacity function κ shown in the lower right corner of each panel). The inner opaque bubble with aconvoluted 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 km s − ). The choice of radial velocity for the isosurface is motivated by tracking the average value of Simass fraction. The box is × , km across, which is the extent of the main Ne shell (see Figure 1). The progression of states (from left toright and top to bottom) shows the violent O-Ne shell merger in progress. components of the specific kinetic energy are defined as E a ( i ) := 12 δm i m i, + (cid:90) m i, − ρv a d m, a ∈ { θ, φ } . (13)The ratio fluctuates with time and hardly stays close to theequipartition value ( E r = E θ = E φ ) for the convectivelyactive shells. For the stable Zone I the value deviates awayfrom / when the shell is finally eroded away. Interestingly,transverse motions dominate over radial motion in Zone IIafter ≈ s. 4.3.2. Entropy: 1D vs 3D
Regions inside a star tend to become convective when ra-diation and neutrinos are unable to carry away the thermalenergy deposited by nuclear burning. In a stably stratifiedregion the entropy gradient is positive (entropy increases ra-dially outwards). Figure 10 shows the entropy profiles forthe 1D model (left panel) and the 3D model (right panel).The differences between the entropy profiles of the 1D and3D models are conspicuous. The negative entropy gradientat the base of the O layer facilitates convection.
HELL MERGER IN
3D 13 E r [ e r gg ] convectively inert Si shellconvective O shellzone Iconvective Ne shellzone IIconvective C shellzone IIIOnset Stationary phase Rapid rise Saturate0 . . . . . . . . . t [10 s] = E r / ( E ✓ + E ) isotropic turbulence onset phase
Upper panel: Specific kinetic energies in the radial velocity fluctuations as functions of time for the different layers. The evolutioncan be split into four phases: a) onset of convection, b) stationary phase, c) rapid rise, and d) saturation phase. The entire region between theinert Si shell and the Ne shell attains similar specific kinetic energy after the shell merger. Lower panel: Ratio of kinetic energy stored in theradial fluctuations ( E r ) and the transverse flow ( E θ + E φ ) as a function of time. The horizontal line marks the value (1/2) when the kineticenergy is equally partitioned between all the components as in case of isotropic turbulence, showing that the flow is primarily dominated bylarge-scale radial motions due to convection. In case of the 1D model, the entropy in the region in-side . − . M (cid:12) increases gradually until the first s( ∆ s ≈ . k b / nucleon), keeping the gradient unchanged. Inthe last s, the inner Fe/Si core experiences significantcontraction ( ≈ leads to more efficient mixing. In case of the3D model the entropy evolution is remarkably different. Theoverall entropy increase is mostly due to the heat depositedby oxygen burning and hot-spot burning. In the first s,the average entropy for the matter inside . − . M (cid:12) de-creases by ∆ s . − . ≈ . k b / nucleon. This high-entropymaterial is carried inwards by convective downdrafts. Dur-ing the same time, the average entropy for the matter inside Note: In K
EPLER , convective boundaries are ’softened’ by one formalovershoot zone with / mixing efficiency as semiconvection. . − . M (cid:12) increases by ∆ s . − . ≈ . k b / nucleon.The ratio of ∆ s . − . and ∆ s . − . is > , which clearlyshows that the entropy gain of the inner region cannot beexplained as resulting from entrainment of high-entropy ma-terial and must be a result of nuclear burning. The entropyprofile inside . − . M (cid:12) is levelled in the first s, but theadditional Ne burning preserves the small negative entropygradient between . − . M (cid:12) . At the same time the steepentropy barrier close to . M (cid:12) is softened (the entropy jumpis reduced). With the complete disappearance of the stabiliz-ing gradient provided by Zone I at ≈ s, thermal energycan be convected into a much larger volume. We define threetimescales relevant for the growth of entropy: the convec-tive timescale ( t conv ), the nuclear timescale ( t nuc ), and theexpansion timescale ( t exp ) which are defined as: t conv := H P (cid:16) (cid:102) v (cid:48)(cid:48) (cid:17) / , t nuc := e th ˙ q net , t exp := (cid:12)(cid:12)(cid:12)(cid:12) d˜ v d r (cid:12)(cid:12)(cid:12)(cid:12) − , (14)4 Y ADAV ET AL . m [ M ] s [ k b / nu c l e o n ]
57 s118 s179 s237 s 297 s358 s417 s 2.0 2.5 3.0 3.5 4.0 m [ M ] s [ k b / nu c l e o n ]
59 s118 s179 s238 s 295 s359 s420 s
Figure 10.
Entropy evolution in the 1D model (left panel) and the 3D model (right panel). The dotted blue curve represents the initial ( t = 0 s)entropy profile in each of the panels. In the 1D model, most of the change is confined to the region between . − . M (cid:12) . Energy depositedby O burning causes a gradual increase in entropy until ≈ s. At the end, the sharp jump in entropy at ≈ . M (cid:12) disappears. In the 3Dmodel efficient heat transport quickly levels the entropy gradient inside . − . M (cid:12) . In the last ≈ s rapid energy deposition by episodicNe 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 inthe right panel is to guide the eye to the inversion of entropy gradient, which happens between s and s. where H P = P d r/ d P is the pressure scale height, and e th is the specific internal energy. Figure 11 shows a com-parison between the nuclear and convective timescales, ξ := t nuc /t conv , in the region inside . M (cid:12) during thelast s before the onset of collapse. The value of ξ ismuch larger than unity everywhere (except close to the innerboundary), which means convection can carry away all thenuclear energy released. However, this can not explain theentropy peak seen at the base of the O shell (see the rightpanel of Figure 10). Now, let us take the contraction of theinner core into account. During the last s, the inner bound-ary of the O shell recedes inwards rapidly, while at the sametime the outer layers move inward at lower velocities. There-fore, a given mass shell located between r and r + δr expandsin size. A convective plume which starts at the base of sucha shell needs to cover a longer distance to leave the shell.This would result in trapping of the convective flow leadingto a reduction in the convective efficiency. Figure 11 showsa comparison of the expansion and convective timescales, ξ := t exp /t conv . During the last seconds before collapsethis ratio approaches unity, which implies that the convectiveflow is trapped. The inability of convection to transfer thereleased nuclear energy leads to a rapid buildup of entropy( ∆ s ≈ . k b / nucleon during the last s) at the base of the Oshell, leading to an overall negative entropy gradient extend-ing up to ≈ . M (cid:12) , which is the base of Zone II. In the endthe 3D model has higher entropy compared to the 1D model.4.4. Mixing and Shell Merger
Figure 12 shows the spatial distribution of Ne (upper row)and Si (lower row) confirming the picture of large-scale mix-ing in the 3D model. The spatial distributions of both Neand Si at ≈ s (panels at s) show significant deviationsfrom spherical symmetry. The dipolar feature (panel at s) in the Ne distribution is a cavity created by a Si-rich bubble(panel at s). The downdrafts of Ne-rich material propa-gating into the Si-enriched inner volume (panel at s) pro-vide the nuclear fuel for hot-spot burning. As Ne mixes in-wards and burns, the Si-rich matter expands outwards result-ing in a complete merger of the Ne shell and the Si-enrichedoxygen shell as seen in panels at s. The Ne mass fractiondecreases considerably ( ≈ . M (cid:12) of Ne is burnt in total) af-ter the merger and the leftover Ne is thoroughly mixed withSi. Figure 13 shows an overview of the merger and expan-sion of the Si and Ne shells. We also show the evolution ofthe typical angular mode number characterizing convectionin the merging layers evaluated using the following relation(Chandrasekhar 1981) (cid:96) := π × r + + r − r + − r − , (15)where r + , r − are the inner and outer radii of the convectiveshell. The value of (cid:96) decreases rapidly as the ongoing shellmerger leads to the formation of a radially extended convec-tive region reaching from the base of the O shell to the outeredge of the convective Ne shell.Figure 14 compares the mass fractions of Si, Ne, and Oin the 1D and 3D models. Although both 1D and 3D mod-els show smoothing of the sharp features, the abundances inthe 1D model are practically unaffected outside of ≈ . M (cid:12) .Abundance changes in the 1D model take place during thelast ≈ s, while in the 3D model they come about slowlyover the course of the simulation. In the 1D model thesechanges results from convective mixing (MLT), overshoot-ing, thermohaline mixing and semiconvection (treated as dif-fusive processes, for details refer to M¨uller et al. 2016b, Sec-tion 2.1). Specifically, during the last s, the core experi-ences significant contraction ( ≈ HELL MERGER IN
3D 15 . . . . . . . . . m [ M (cid:12) ] − − ξ ξ := t nuc /t conv ξ := t exp /t conv “freezeout”“trapping” h ξ i
390 s360 s h ξ i
390 s360 s
420 s420 s
Figure 11.
A comparison between nuclear ( t nuc ), convective( t conv ), and expansion ( t exp ) timescales during the last s lead-ing to collapse. Thick solid curves show the average value of ξ := t nuc /t conv and ξ := t exp /t conv over t ∈ [360 , s, and thethick dashed curves show the value of ξ and ξ just before the on-set of collapse. The ratio of nuclear timescale and the convectiveturnover timescale ξ (cid:29) everywhere, with the exception of the in-ner boundary. Close to the onset of collapse ξ ∼ in the regioninside . M (cid:12) , pointing to a rapid decline in the efficiency of con-vective heat transfer. This effective trapping of the convective flowproduces a peak in entropy (during the last s entropy changes by ≈ . k b / nucleon in the innermost region) at the base of the O burn-ing shell. sity and temperature increment leads to an enhancement inthe neon burning rate, which in turn leads to an increase inthe convective velocity resulting in faster mixing (accordingto the MLT model). In the 3D model they result from massentrainment at interfaces separating convectively stable andunstable regions, the large-scale convective flow, and the ris-ing Si-rich bubbles powered by rapid Ne burning.Panels a and d compare the O abundance profiles. In the3D model, the O mass fraction inside . M (cid:12) increases dur-ing the last s. The increase results from a combinationof large-scale mixing (convective downdrafts carrying O-richmaterial inwards, ≈ . M (cid:12) ) and rapid Ne burning (produc-ing O via Ne ( γ, α ) O), and negates the change caused byburning of O to Si. In the end, the total mass of O in thesimulated region does not change.The 1D model fails to capture the evolution of the O abun-dance altogether. Panels b and e compare the Ne abundanceprofiles. The outer part of the neon-carbon shell locatedbetween . − . M (cid:12) remains unaffected in both cases. In the 3D model, the gradual rise in the Ne mass fraction in-side . M (cid:12) continues till ≈ s only to be reversed by aphase of rapid Ne burning (described in Section 4.3.1). Incontrast, Ne is virtually unaffected in the 1D model. Pan-els c and f compare the Si abundance profiles. In the 3Dmodel, the evolution of silicon mass fraction is a result oftwo processes: 1) silicon production by oxygen burning via O ( O , α ) Si, and 2) outward mixing caused by convec-tion. During the course of evolution ≈ . M (cid:12) of Si isproduced and ≈ . M (cid:12) of Si is gradually transported outfrom the region between . − . M (cid:12) into the region be-tween . − . M (cid:12) . As a result the outer boundary of theSi-enriched layers from ≈ , km to ≈ , km (closeto the C/He interface).4.5. Initial Transients and Relaxation
Mapping data from one code to another invariably leads totransients. At this point, we must delineate transients causedby mapping errors ( spurious ), and the transients caused byrapid changes in the thermodynamical properties of a system( genuine ).First, we provide two arguments why our results are notaffected by mapping errors:1.
Convective Turnover Time : In order to smooth (damp)mapping related errors, one aims at simulating a suf-ficiently large number of convective turnovers of theflow. This approach to eliminate the effects of ini-tial transients (spurious) in a simulation is useful whenconvection is supposed to lead to a steady-state. Theconvective turnover timescale is conventionally de-fined as τ C := 2 H P (cid:16) (cid:102) v (cid:48)(cid:48) (cid:17) / , (16)where H P is the local pressure scale height, v (cid:48)(cid:48) r is thefluctuating component of the radial velocity. Figure 16shows the value of the convective turnover timescaleaveraged over a period of s at various times. Thetimescale evolves rapidly over the course of the sim-ulation. In the innermost convectively active region,which is the Si-rich part of the O shell (see Fig-ure 1), the convective turnover timescale is between ≈ − s during the entire simulation, which is suf-ficient for − full turnovers during the simulation.In the Ne shell, the turnover timescale becomes rel-atively stable at around s and its value is (cid:46) sthereafter. Therefore, even after the onset of the shellmerger, which begins at ≈ s (see section 4), the re-gion has enough time to go through − convectiveturnovers. Furthermore, the convectively stable layerssandwiched between the convectively unstable layers(see Section 4.1) play an effective role in absorbingand dissipating the initial transients because they dampthe propagation of acoustic waves.2. Growth of Kinetic Energy In the Merging Shells : Dur-ing the onset phase, between − s (Figure 9), the6 Y ADAV ET AL . y [ k m ]
268 sNe +0.05+0.10+0.15+0.20
202 323 sNe +0.05+0.10+0.15+0.20 x [10 km] y [ k m ]
268 sSi +0.05+0.10+0.15+0.20+0.25+0.30+0.35 x [10 km]
202 323 sSi +0.05+0.10+0.15+0.20+0.25+0.30+0.35 x [10 km] (a)
Slices ( x − y plane) showing the spatial distribution of Ne (upper row) and Si (lower) at various times compared to the volume-rendered 3D views in Figure 8. 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 insidea spherical region bounded by a convectively stable layer (Zone II at , s and Zone III for s) for clarity of presentation. Please seethe animation provided as a supplementary material (Movie B). kinetic energy in radial velocity fluctuations grows andreaches a steady value in the layers of interest, whichare the Si-rich part of the O shell (labelled as “con-vective O shell”) and the convective Ne shell, imply-ing that most of the initial transients (mapping errors)have already dissipated away and do not affect the fur-ther evolution.In the case of a genuine transient, the system may undergoan irreversible change in a short time, and may never ap-proach a steady-state, which is also the case with the shell-merger described above. In this case, the numerical solutionis stable (the CFL value is 0.6; the code preserves the hydro-static state; PPM is well understood), but the physical systemitself is on the verge of a thermodynamical instability. Thefact that the system under consideration here manifests sucha behaviour is evident from the nature of the transition be-tween the ‘stationary’ phase and the ‘rapid rise’ phase (seeFigure 9). Before this transition, the kinetic energy in theconvective O shell and Zone-I stays at a plateau between ≈ − s, and in the Ne-shell energy is slowly depositedby convective activity underneath and gently reaches a steadystate value at ≈ s. For the next ≈ s the change in ki-netic energy in the Ne shell is negligible. Its value starts ris-ing again only at the time of the transition in response to theelevated energy production in the O burning shell, associatedwith the shell merger beginning at ≈ s. Therefore, the ensuing shell merger represents a “genuinely transient phe-nomenon” consistent with the thermodynamical evolution. Itis obvious that the conditions emerging from the shell mergerdo not represent a steady state situation, and we cannot claimthat the same conditions or the same dynamical situationwould be present at the onset of collapse if we had startedour simulation significantly earlier. However, we considerthe shell merger starting ≈ s before collapse as an exem-plary case of such a possibility, which could happen in otherstars even if not in the one considered here. In Section 6, wefurther discuss the effects of our choice of the starting timeon the onset and dynamics of the shell merger and its impacton the collapse conditions. PRE-SUPERNOVA MODEL PROPERTIESLarge-scale modes with (cid:96) ∼ − are most effective forshock revival (M¨uller & Janka 2015; Abdikamalov et al.2016; M¨uller et al. 2016a). In order to gain a quantitative un-derstanding of the modes, we have done spherical harmonicdecomposition of the density and velocity perturbations. Fig-ure 15 shows the multipole power distribution at collapse forthe density fluctuations (left panel) and the velocity fluctua-tions (right panel). The red markers indicate the multipolewith the largest power ( (cid:96) max ) at a given mass coordinate. Thevalue of (cid:96) max changes with radius, but there are continuoussections where it remains stationary. In the case of densityfluctuations inside the merged O-Ne shell ( . − . M (cid:12) ) a HELL MERGER IN
3D 17 t [10 s] r [ k m ] m = 2 . r Ne , − r Ne , + r Si − Ne interface r Si , − r Si , + ‘ Figure 13.
Time evolution of the position of the inner boundaries(solid lines, symbol ‘ − ’) and outer boundaries (dashed lines, sym-bol ‘ + ’) of the Ne and Si layers. The boundaries have been de-termined using angular averaged mass fraction profiles. The innerboundary of the Ne shell is defined at X Ne = 0 . and the innerboundary of the Si layer is defined at X Si = 0 . . The outer bound-aries for both Ne and Si are defined at X Ne = 0 . or X Si = 0 . ,respectively. The evolving radius corresponding to the initial masscoordinate of the Ne/Si interface ( m ≈ . ) is shown as the dot-ted green curve. In the linear regime, the outer boundary of Sishell moves outwards slowly, whereas in the non-linear phase (shellmerger) the motion is much faster. The late expansion of the outerNe shell boundary is caused by an expansion of the whole O-Nelayer due to the violent energy release from enhanced burning. Thesolid black line shows an estimate of the typical angular mode num-ber of convection cells according to Equation 15 (evaluated for theentire Si-enriched part of the O layer; see last panel of Figure 12). mixture of dipole ( (cid:96) = 1 ) and quadrupole mode ( (cid:96) = 2 ) dom-inates, and in the region outside of Zone II the dipole modedominates (similar to the global non-spherical oscillation de-scribed by Herwig et al. (2014)). This confirms our visual im-pression obtained from Figure 3. In the case of radial velocityfluctuations a mixture of dipole mode and (cid:96) = 3 mode dom-inates in the merged O-Ne shell, whereas the dipole modedominates in the C and He shells. The normalized probabil-ity distribution of (cid:96) max (shown as histograms) confirms that (cid:96) ∼ − dominate both the density fluctuations and the ve-locity fluctuations. Figure 17, left panels, show another viewof the spectra at the middle points of the layers defined inSection 4.2.2. In the lower half-panel the spectra for the re-gion inside ≈ . M (cid:12) (layers c and d) are shallow comparedto the region outside and including Zone II. In the upper half-panel the spectra for (cid:96) (cid:46) are shallower in the merged O-Ne shell including Zone II. The spectra above (cid:96) (cid:38) agree quitewell with the ‘Kolmogorov slope’.Figure 17 (right panel) shows the profiles of the turbulentMach number corresponding to the radial velocity fluctua-tions at multiple times. The quantity is defined followingM¨uller et al. (2016b) as (cid:103) M r := (cid:82) Ω ρ ( v r − ˜ v r ) dΩ (cid:82) Ω ρc s dΩ , (17)where c s is the adiabatic sound speed. The Mach number inthe O and Ne shells increases gradually until the shell merger.During the first s, the value of M r is almost negligiblein the Ne shell outside . M (cid:12) , which is due to the smallerBrunt-V¨ais¨al¨a in the Ne shell. The region outside . M (cid:12) stays quiescent until s as the stable Zone II located be-tween . − . M (cid:12) does not allow convective plumes to pen-etrate further. After the merger, between − s, M r increases rapidly by a factor of reaching a peak value of ≈ . at the onset of collapse. Convection is so violentthat the plumes penetrating deep into Zone II excite stronginterfacial waves that reach as far as the He shell, visibleby the non-negligible Mach numbers outside . M (cid:12) . Fig-ure 18 shows pseudocolor maps of the radial velocity ondifferent Si mass-fraction isosurfaces prior to collapse. Aswe go from small to large values of X Si , we effectivelyprobe the Si distribution in successively deeper regions ofthe star. The left panel shows large-scale Si-rich bubbles ris-ing at moderate speeds of ≈ km s − somewhere around ≈ , − , km. In the panel for X Si = 0 . (mid-dle), we see a highly asymmetrical and clumpy distributionof Si with some of the clumps moving outwards at veloc-ities in excess of , km s − . This may have an impacton the observed asymmetries in the element distributions insome young supernova remnants such as Cassiopeia A. Inthe panel for X Si = 0 . (right) we probe an even higher Simass fraction which is present deep inside the O shell (scaleof ≈ , km), where we again see fast outward moving Si-rich plumes. In short, the Si distribution in bulk is clumpyand has asymmetric velocities (although it looks as if it hadisotropized in panel f of Figure 12). CONCLUSIONSIn this work, we present the first 3D simulation (full π solid angle) of a fully developed large-scale O-Ne shellmerger prior to core-collapse in an . M (cid:12) progenitor. Thework builds upon the previous study of M¨uller et al. (2016b)in which they simulated the last minutes of O-shell burningin an M (cid:12) supernova progenitor up to the onset of corecollapse. In the present case we find significant differencesbetween the kinematic, thermodynamic, and chemical evolu-tion of the 3D and the 1D models. Ne is mixed deep intothe O-shell, leading to a fuel ingestion episode which re-leases significant amounts of thermal energy by nuclear burn-ing on a short timescale ( ∼ s), further boosting convec-tion: Convective downdrafts transport Ne close to the base8 Y ADAV ET AL . (a) (b) (c)(d) (e) (f) m [ M ] m a ss f r a c t i o n , X O m [ M ] m a ss f r a c t i o n , X N e m [ M ] m a ss f r a c t i o n , X S i
57 s118 s179 s237 s 297 s358 s417 s2.0 3.0 4.0 5.0 m [ M ] m a ss f r a c t i o n , X O
59 s118 s179 s238 s 295 s359 s420 s 2.0 3.0 4.0 5.0 m [ M ] m a ss f r a c t i o n , X N e m [ M ] m a ss f r a c t i o n , X S i Figure 14.
Mass 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 . M (cid:12) , whereas such a feature is absent in the 1D model(panel a). Panel e shows the sudden decrease in the Ne abundance ( ≈ . M (cid:12) ) in the last s 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 ≈ , km in the 3D model over the course of simulation. In contrast, the 1D model (panel c) fails to capture the large-scale mixing of Sientirely. Note: The initial abundance profiles are shown in Figure 1 (top panel). . . . . . . m [ M ] log S ` -5.0-4.0-3.0-2.0-1.00 5 10 15 20 25 30 ` . . . f r e q u e n c y . . . . . . m [ M ] log S ` ` . . . f r e q u e n c y multipole moment , `
Power in various multipoles ( S (cid:96) ) of the density fluctuations ( ρ (cid:48) /ρ , left panel) and the radial velocity fluctuations ( v (cid:48)(cid:48) r , 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 thebottom of each panel shows the probability density of multipoles with maximum power in the entire simulated shell. Note: Tick-marks on (cid:96) axis are placed on left edge of a bin, such that the first bin on (cid:96) axis corresponds to (cid:96) = 0 . Mathematically, the spherical average of a fluctuation(Reynolds averaged) should be zero. In this case, we get a small but non-zero value because of numerical reasons. HELL MERGER IN
3D 19 . . . . . . m [ M (cid:12) ] τ C [ s ] onset of shell mergerO shell Ne shell C s h e ll Z o n e − I Zone − II Zone − III h −
060 s ih −
120 s ih −
180 s ih −
240 s i h −
300 s ih −
360 s ih −
420 s i Figure 16.
Time averages of convective turnover times (see Equation 16) at different instants. The boundaries of the O shell, the Ne shell, andthe C shell are marked by black, red, and gray vertical lines, respectively. The curve at s shows that the convective turnover timescale in thelayer between the O shell and Zone-I (see Table 2 for the definition) is ≈ − s, which gives a conservative estimate of ≈ − convectiveturnovers in the remaining s before the collapse. In the Ne-shell, the convective turnover time is (cid:46) s, which provides enough time for atleast − convective turnovers in the remaining time before the collapse. The horizontal dashed line marks s, which is the total duration ofthe simulation. The shell merger (described in more detail in Section 4) begins at ≈ s (shown here by the solid cyan curve). . . . . . . . m [ M ] . . . . . . . . ⇣ g M r ⌘ / S i / N e i n t e r f a c e O / H e i n t e r f a c e
59 s118 s179 s238 s 295 s359 s420 s10 S ` / S , v r ` / multipole moment , ` S ` / S , ⇢ / ⇢ cde fgh p M r
Left panels: Multipole power spectra ( S ( (cid:96) ) /S (0) ) for density fluctuations (lower half-panel) and radial velocity fluctuations (upperhalf-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 thecomputed 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 squareMach 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 isat boundary of layer g and layer h. Up to about s, distinct convective O and Ne shells can be recognized, both with strongly subsonicconvection. After the shell merger, the convective Mach numbers increase considerably. ADAV ET AL . v r [100 km s ]
000 km
000 km
000 km
Pseudocolor maps of radial velocity v r on Si mass fraction isosurfaces for X Si = 0 . (left), X Si = 0 . (middle) and X Si = 0 . (right),which effectively probe the Si distribution at successively greater depths. The snapshots are taken at s when the Fe-core begins to undergogravitational collapse. “Scale” refers to the size of the y -axis (vertically upwards). The left panel shows the large-scale Si-rich bubbles. Themiddle panel shows the asymmetry and clumpy Si distribution prior to collapse, and the right panel shows Si-rich plumes and downdrafts in thedeep interior moving outward or inward, respectively, at ≈ , km s − . Please see the animation provided as a supplementary material forthe case with X Si = 0 . (Movie C). of the O shell, where it ignites and powers convection in re-turn, leading to a positive feedback loop for rapid Ne en-trainment. The entire simulated shell attains a convectiveMach number of ≈ . . The maximum convective velocitiesare ∼ km s − for updrafts and ∼ km s − for down-drafts. In contrast the 1D model is relatively quiescent withconvective shell burning constrained to layers. The specifickinetic 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 mergedshell at ( ≈ s ) , the distribution of elements becomes moreisotropic in the following time leading up to collapse. How-ever, the velocity field still shows a large global asymmetryat collapse. In short, the 3D pre-supernova model exhibitsstrong density perturbations and large-scale velocity asym-metries (with dominant (cid:96) ∼ − modes) in the flow as well asthe distribution of nuclear species.Another question worth investigating is that, “how the choice of starting time affects the onset and dynamics of theshell merger?” Although, we can definitely say that the shellmerger reported in this study is a robust phenomenon, the ex-act time at which it happens may perhaps change if we startat a different time. This question can only be answered bydoing a longer simulation, which is at the moment unfeasibledue to limited computing resources. Nevertheless, the shellmerger must occur quite late during the evolution of the Oshell because this is when the critical conditions (a reductionof the entropy jump between the O shell and Ne shell, andfast O shell convection that allows fast entrainment) are es-tablished. This point has been discussed in detail in Collinset al. (2018): For the entropy of the O shell to rise and be-come similar to that of the Ne shell, O burning and neutrinocooling must have already dropped out of balance due to thecontraction of the O shell. The contraction of the O shell in the wake of core contraction allows convection to reach suf-ficiently high Mach numbers and power strong entrainment.Simulations by Couch et al. (2015) and M¨uller et al.(2016b) have already shown that the progenitor structureis genuinely three-dimensional at collapse, but only consid-ered quasi-steady state convection. Recently, Yoshida et al.(2019) have presented a suite of 2D and 3D models of O-shell burning for massive stars between − M (cid:12) . They findconvection with large-scale eddies and the turbulent Machnumber ∼ . in their models despite the simulations beingshort (lasting ≈ s before collapse), which neverthelessshows that large-scale perturbations may be a generic featurein pre-supernova stars. Our work demonstrates that asym-metries can become exceptionally large when the convectiveflow becomes non-steady because of a shell merger. In sucha case, the nuclear timescale becomes comparable to theconvective timescale, and hence nuclear burning becomesstrongly coupled with the flow dynamics and the resultantmixing is considerably faster and can be highly asymmet-ric. 1D stellar evolution calculations are severely limited incapturing these phenomena self-consistently as they occurprimarily due to a combination of turbulent entrainment, fuelingestion, and the excitation/propagation of internal waves,which are inherently 3D effects. Despite the inherent limi-tations of 1D calculations, such shell mergers may still takeplace in 1D (as seen by Rauscher et al. (2002, model S20) andin recent studies by Sukhbold & Woosley (2014) and Daviset al. (2018, case C3), and merely be delayed as in modelS25 of Rauscher et al. (2002), where the shell merger hap-pens only s before collapse. However, shell mergers in 1Dmodels cannot capture the highly dynamical and asymmetricflow during and after the merger. Therefore, the dynamics ofshell mergers during the late burning stages can only be cap-tured using self-consistent multi-dimensional calculations. HELL MERGER IN
3D 21Such calculations are still in their infancy and the work re-ported here attempts to make progress in our understandingof the internal structure of supernova progenitors.Also, recent studies show that progenitor asymmetries mayhelp in shock revival (Couch & Ott 2013; M¨uller & Janka2015; M¨uller et al. 2017) by aiding the growth of convectionand/or instabilities like the standing accretion shock instabil-ity (SASI) (Takahashi et al. 2016; M¨uller et al. 2017). Ourwork furnishes further evidence that some supernova pro-genitors have convective seed perturbation that are dynam-ically important. For a shell merger like the one described inthis work, we see even more violent convection and strongerlarge-scale asymmetries than in the . M (cid:12) model of M¨ulleret al. (2017).Shell mergers will also affect the compactness parameter(O’Connor & Ott 2011) and thus the explodability with pos-sible consequences for the location and extent of “islands ofexplodability” (e.g. Sukhbold & Woosley (2014)). Due tothe high convective Mach numbers and pronounced globalasymmetries, we expect a strong impact on the explosion dy-namics for our . M (cid:12) model.This work constitutes only the first step towards under-standing the importance of late-stage shell mergers. For ex-ample, mergers may also impact the pre-supernova nucle-osynthesis yields, as already suggested by 1D stellar evo-lution models. The S20 model of Rauscher et al. (2002)showed a strong overproduction of several elements betweenSi and V and an underproduction of Cr, Mn, and the light Feisotopes, which they ascribed to the stellar structure result-ing from the C/O shell merger. Ritter et al. (2018) showedthat even moderate C-ingestion events during O burning cansignificantly overproduce odd-Z elements such as K, Sc, Cl,etc. Some of their models with full C-O shell merger (albeitartificially increased O burning luminosity) have overproduc-tion factors greater than 10. Therefore, we expect to see sub-stantial effects on the nucleosynthetic yields due to the atyp-ical way in which Ne burns during the shell merger event described in our work. During the explosion phase, the sub-stantial asymmetries in the element distribution and the ve-locity field may have repercussions for observable asymme-tries in supernovae and nucleosynthesis (Hughes et al. 2000;DeLaney et al. 2010; Lopez & Fesen 2018). Such questionswill be addressed in the future.ACKNOWLEDGEMENTSThis project was supported by the European ResearchCouncil through grant ERC-AdG No. 341157-COCO2CASAand by the Deutsche Forschungsgemeinschaft (DFG, Ger-man Research Foundation) under Germany’s ExcellenceStrategy through Excellence Cluster ORIGINS (EXC-2094)–390783311. This work was also supported by the Aus-tralian Research Council through ARC Future FellowshipsFT160100035 (BM), Future Fellowship FT120100363 (AH).This research was undertaken with the assistance of resourcesobtained via NCMAS and ASTAC from the National Com-putational Infrastructure (NCI), which is supported by theAustralian Government and was supported by resources pro-vided by the P AWSEY S UPERCOMPUTING C ENTRE withfunding from the Australian Government and the Govern-ment of Western Australia. On the Garching side, computingresources from the G
AUSS C ENTRE FOR S UPERCOMPUT - ING (at the L
EIBNIZ S UPERCOMPUTING C ENTRE (LRZ)under SuperMUC project ID: pr53yi) are acknowledged.The analysis of the simulation was done using the MPGSupercomputer H
YDRA , C
OBRA , and D
RACO provided byT HE M AX P LANCK C OMPUTING AND D ATA F ACILITY atGarching.
Software:
NumPy (van der Walt et al. 2011), Scipy(Oliphant 2007), IPython (P´erez & Granger 2007), Mat-plotlib (Hunter 2007), Visit (Childs et al. 2012), and SHTnslibrary (Schaeffer 2013). This research has made use ofNASA’s Astrophysics Data System.APPENDIX A. ANGLE AVERAGING METHODAn intensive physical quantity (for e.g., density) can be written as a sum of its average and the fluctuating component (notationfor mean and fluctuating components follows from Oertel (2010, chapter 5)) as f ( r ) = f + f (cid:48) ( r ) , where the angle-averagedvalue (Reynolds) f is defined as f := (cid:104) f (cid:105) = 14 π (cid:90) Ω f ( r ) dΩ . (A1)An extensive physical quantity (for e.g., specific energy e ) can be written as g ( r ) = ˜ g + g (cid:48)(cid:48) ( r ) , where the angle-averaged value(Favre) is defined as ˜ g := (cid:104) g (cid:105) = (cid:82) Ω ρ ( r ) g ( r ) dΩ (cid:82) Ω ρ ( r ) dΩ . (A2)The fluctuating components f (cid:48) and g (cid:48)(cid:48) satisfy the following relations f (cid:48) = 14 π (cid:90) Ω f (cid:48) ( r ) dΩ = 0 , g (cid:48)(cid:48) = 14 π (cid:90) Ω g (cid:48)(cid:48) ( r ) dΩ (cid:54) = 0 , and , ρg (cid:48)(cid:48) = 14 π (cid:90) Ω ρ ( r ) g (cid:48)(cid:48) ( r ) dΩ = 0 . (A3)2 Y ADAV ET AL . B. SPHERICAL HARMONIC TRANSFORMSThe forward spherical harmonic transform of a function f ( θ, φ ) involves computing the coefficients c m(cid:96) according to c m(cid:96) = (cid:90) Ω Y m(cid:96) ∗ ( θ, φ ) f ( θ, φ ) dΩ , (B4)where Y m(cid:96) ( θ, φ ) is the orthonormalized spherical harmonic function of degree (cid:96) and order m . The total power c (cid:96) for a mode withmultipole order l is defined as S (cid:96) := c (cid:96) = (cid:96) (cid:88) m = − (cid:96) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω Y m(cid:96) ( θ, φ ) f ( θ, φ ) dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (B5)We have used the routines from the SHTns library (Schaeffer 2013) to calculate the spherical harmonic decomposition.REFERENCES Abdikamalov, E., Zhaksylykov, A., Radice, D., & Berdibek, S.2016, MNRAS, 461, 3864, doi: 10.1093/mnras/stw1604Andrassy, R., Herwig, F., Woodward, P., & Ritter, C. 2018, arXive-prints, arXiv:1808.04014. https://arxiv.org/abs/1808.04014Arnett, D. 1994, ApJ, 427, 932, doi: 10.1086/174199Arnett, W. D., Meakin, C., Hirschi, R., et al. 2019, ApJ, 882, 18,doi: 10.3847/1538-4357/ab21d9Asida, S. M., & Arnett, D. 2000, ApJ, 545, 435,doi: 10.1086/317774Bazan, G., & Arnett, D. 1994, ApJL, 433, L41,doi: 10.1086/187543—. 1998, ApJ, 496, 316, doi: 10.1086/305346Biermann, L. 1932, Zeitschrift fur Astrophysik, 5, 117B¨ohm-Vitense, E. 1958, Zeitschrift fur Astrophysik, 46, 108Boris, J. P., Grinstein, F. F., Oran, E. S., & Kolbe, R. L. 1992, FluidDynamics Research, 10, 199,doi: 10.1016/0169-5983(92)90023-PBuras, R., Rampp, M., Janka, H.-T., & Kifonidis, K. 2006, A&A,447, 1049, doi: 10.1051/0004-6361:20053783Cantiello, M., Mankovich, C., Bildsten, L., Christensen-Dalsgaard,J., & Paxton, B. 2014, ApJ, 788, 93,doi: 10.1088/0004-637X/788/1/93Chandrasekhar, S. 1981, Hydrodynamic and HydromagneticStability, Dover Books on Physics Series (Dover Publications)Childs, H., Brugger, E., Whitlock, B., et al. 2012, in HighPerformance Visualization–Enabling Extreme-Scale ScientificInsight (United States. Department of Energy. Office ofScience), 357–372Colella, P., & Sekora, M. D. 2008, Journal of ComputationalPhysics, 227, 7069, doi: 10.1016/j.jcp.2008.03.034Colella, P., & Woodward, P. R. 1984, Journal of ComputationalPhysics, 54, 174, doi: 10.1016/0021-9991(84)90143-8Collins, C., M¨uller, B., & Heger, A. 2018, MNRAS, 473, 1695,doi: 10.1093/mnras/stx2470 Couch, S. M., Chatzopoulos, E., Arnett, W. D., & Timmes, F. X.2015, ApJL, 808, L21, doi: 10.1088/2041-8205/808/1/L21Couch, S. M., & Ott, C. D. 2013, ApJL, 778, L7,doi: 10.1088/2041-8205/778/1/L7Cristini, A., Meakin, C., Hirschi, R., et al. 2017, MNRAS, 471,279, doi: 10.1093/mnras/stx1535Davis, A., Jones, S., & Herwig, F. 2018, MNRAS,doi: 10.1093/mnras/sty3415DeLaney, T., Rudnick, L., Stage, M. D., et al. 2010, ApJ, 725,2038, doi: 10.1088/0004-637X/725/2/2038Dimotakis, P. E. 2005, Annual Review of Fluid Mechanics, 37,329, doi: 10.1146/annurev.fluid.36.050802.122015Fernando, H. J. S. 1991, Annual Review of Fluid Mechanics, 23,455, doi: 10.1146/annurev.fl.23.010191.002323Fuller, J., Cantiello, M., Lecoanet, D., & Quataert, E. 2015, ApJ,810, 101, doi: 10.1088/0004-637X/810/2/101Grefenstette, B. W., Harrison, F. A., Boggs, S. E., et al. 2014,Nature, 506, 339, doi: 10.1038/nature12997Grefenstette, B. W., Fryer, C. L., Harrison, F. A., et al. 2017, ApJ,834, 19, doi: 10.3847/1538-4357/834/1/19Heger, A., & Woosley, S. E. 2010, ApJ, 724, 341,doi: 10.1088/0004-637X/724/1/341Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350,doi: 10.1086/429868Herwig, F., Pignatari, M., Woodward, P. R., et al. 2011, ApJ, 727,89, doi: 10.1088/0004-637X/727/2/89Herwig, F., Woodward, P. R., Lin, P.-H., Knox, M., & Fryer, C.2014, ApJL, 792, L3, doi: 10.1088/2041-8205/792/1/L3Hughes, J. P., Rakowski, C. E., Burrows, D. N., & Slane, P. O.2000, ApJ, 528, L109, doi: 10.1086/312438Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90,doi: 10.1109/MCSE.2007.55Iliadis, C. 2015, Nuclear Physics of Stars (Wiley-VCH Verlag,Wenheim, Germany, 2015.), doi: 10.1002/9783527692668
HELL MERGER IN
3D 233D 23