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

TL;DR

This study uses a detailed 3D simulation to explore large-scale mixing and asymmetries in a pre-supernova oxygen-neon shell merger, revealing significant differences from 1D models and implications for supernova dynamics.

Contribution

It provides the first self-consistent 3D simulation capturing large-scale mixing and asymmetries in a pre-supernova shell merger, highlighting the importance of 4π geometry.

Findings

3D model shows vigorous convection with Mach numbers ~0.1.

Efficient mixing causes silicon to spread into the oxygen/neon shell.

Large-scale asymmetries develop, affecting collapse and explosion dynamics.

Abstract

We present a seven-minute long $4\pi$-3D simulation of a shell merger event in a non-rotating $18.88\, M_\odot$ 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 $4\pi$ 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 $\mathord{\approx}\, 0.1$, 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.

Tables

Table 1. Table 1: Setup details for the 3D simulation.

Quantity	Option/Value
aa$r_{i+1}=\left(1+\alpha \right)r_i$, where $1+\alpha ={\left(r_{+}/r_{-}\right)}^{1/N_r}$.Radial grid	Geometric
Cells in radial direction, $N_r$	450
bbThe number of angular cells are specified per patch of the Yin-Yang grid.Cells in polar direction, $N_{\theta }$	48
bbThe number of angular cells are specified per patch of the Yin-Yang grid.Cells in azimuthal direction, $N_{\varphi }$	148
ccFor each grid patch (Yin and Yang).No. of ghost cells ($\theta ,\varphi$), $N_g$	8
ccFor each grid patch (Yin and Yang).No. of buffer cells ($\theta ,\varphi$), $N_b$	2
Angular resolution	$2^{\circ }$
ddThe inner and outer boundaries are moved according to the trajectory obtained from 1D model.Inner boundary at $t=0$, $r_{-}$	$3.30\times {10}^3\text{km}$
ddThe inner and outer boundaries are moved according to the trajectory obtained from 1D model.Outer boundary at $t=0$, $r_{+}$	$9.88\times {10}^4\text{km}$
Inner boundary condition, radial	refer to Section 3.1
Outer boundary condition, radial	refer to Section 3.1
Gravitational potential	Newtonian, Spherical
Simulation time	$420\text{s}$
ee$\delta \rho =f\rho U[-1,1]$, where $U[-1,1]$ is a uniformly distributed random number in the interval $[-1,1]$.Perturbation amplitude ($f$)	$2\times {10}^{-4}$
Neutrino cooling	Yes
Equation of state	Helmholtz EoS
ffCFL is the value of the Courant–Friedrichs–Lewy condition.CFL	0.6

Table 2. Table 2: Convectively stable zones according to the Ledoux criterion at the beginning of the 3D simulation.

Zone	Location (inner edge)		Width
	$m\left[{\text{M}}_{\odot }\right]$	$r\left[\text{km}\right]$	$\Delta m\left[{\text{M}}_{\odot }\right]$	$\Delta r\left[\text{km}\right]$
I	2.4	$7,500$	0.1	$1,000$
IIaaZone II in the 1D model shows a strong and narrow barrier () close to $\approx \mathrm{\hspace{0.17em}4.1}{\text{M}}_{\odot }$.	3.7	$21,500$	0.4	$10,000$
III	4.4	$49,000$	0.1	$5,000$

Equations

$${\omega }_{\text{BV}}^2=\left(\frac{\partial \ln \rho }{\partial r}-\frac{1}{{\Gamma }_1}\frac{\partial \ln p}{\partial r}\right)g,$$

$${\Gamma }_1={\left(\frac{\partial \ln p}{\partial \ln \rho }\right)}_{s,X_i},$$

$${\omega }_{\text{BV}}$$

$$C_{\mathrm{L}}$$

$${\overline{\omega }}_{\text{BV}}(r):=\frac{1}{4\pi }\underset{\Omega }{\int }{\omega }_{\text{BV}}(\mathbf{r})\mathrm{d}\Omega ,$$

$$\frac{{\rho }^{\prime }}{\overline{\rho }}:=\frac{\rho (\mathbf{r})-\overline{\rho }}{\overline{\rho }},$$

$$\mathbf{v}(\mathbf{r})={\tilde{v}}_r(r){\mathbf{e}}_{\mathbf{r}}+v_r^{\prime \prime }(\mathbf{r}){\mathbf{e}}_{\mathbf{r}}+v_{\theta }{\mathbf{e}}_{\theta }+v_{\varphi }{\mathbf{e}}_{\varphi },$$

$$v_r^{\prime \prime }(\mathbf{r}):=v_r(\mathbf{r})-{\tilde{v}}_r.$$

$$C(r,\delta r):=\frac{\overline{{v^{\prime }}_r(r,\theta ,\varphi ){v^{\prime }}_r(r+\delta r,\theta ,\varphi )}}{{\left\{\overline{{v^{\prime }}_r^2(r,\theta ,\varphi )}\overline{{v^{\prime }}_r^2(r+\delta r,\theta ,\varphi )}\right\}}^{1/2}},$$

$${\Lambda }_{\text{corr}}(r):=\underset{r_{-}}{\overset{r_{+}}{\int }}C(r,\delta r)\mathrm{d}{r}^{\prime },\delta r=r^{\prime }-r,$$

$${}^{20}\text{Ne}{+}^{20}\text{Ne}{\to }^{16}\text{O}{+}^{24}\text{Mg},$$

$$E_{\text{r}}(i):=\frac{1}{2\delta m_i}\underset{m_{i,-}}{\overset{m_{i,+}}{\int }}\rho v_r^{\prime \prime 2}\mathrm{d}m,$$

$$E_{\text{a}}(i):=\frac{1}{2\delta m_i}\underset{m_{i,-}}{\overset{m_{i,+}}{\int }}\rho v_{\text{a}}^2\mathrm{d}m,\text{a}\in \{\theta ,\varphi \}.$$

$$t_{\mathrm{conv}}:=\frac{H_P}{{\left(\widetilde{v_{\mathrm{r}}^{\prime \prime 2}}\right)}^{1/2}},t_{\mathrm{nuc}}:=\frac{e_{\mathrm{th}}}{{\dot{q}}_{\mathrm{net}}},t_{\exp }:={\left|\frac{\mathrm{d}\tilde{v}}{\mathrm{d}r}\right|}^{-1},$$

$$\ell :=\frac{\pi }{2}\times \frac{r_{+}+r_{-}}{r_{+}-r_{-}},$$

$${\tau }_{\mathrm{C}}:=\frac{2H_P}{{\left(\widetilde{v_{\mathrm{r}}^{\prime \prime 2}}\right)}^{1/2}},$$

$$\widetilde{{\mathcal{M}}_r^2}:=\frac{\underset{\Omega }{\int }\rho (v_r-{\tilde{v}}_r{)}^2\mathrm{d}\Omega }{\underset{\Omega }{\int }\rho c_{\text{s}}^2\mathrm{d}\Omega },$$

$$\overline{f}:=⟨f⟩=\frac{1}{4\pi }\underset{\Omega }{\int }f(\mathbf{r})\mathrm{d}\Omega .$$

$$\tilde{g}:=⟨g⟩=\frac{\underset{\Omega }{\int }\rho (\mathbf{r})g(\mathbf{r})\mathrm{d}\Omega }{\underset{\Omega }{\int }\rho (\mathbf{r})\mathrm{d}\Omega }.$$

$$\overline{f^{\prime }}=\frac{1}{4\pi }\underset{\Omega }{\int }{f}^{\prime }(\mathbf{r})\mathrm{d}\Omega =0,\overline{g^{\prime \prime }}=\frac{1}{4\pi }\underset{\Omega }{\int }{g}^{\prime \prime }(\mathbf{r})\mathrm{d}\Omega \ne 0,\text{and},\overline{\rho g^{\prime \prime }}=\frac{1}{4\pi }\underset{\Omega }{\int }\rho (\mathbf{r})g^{\prime \prime }(\mathbf{r})\mathrm{d}\Omega =0.$$

$$c_{\ell }^m=\underset{\Omega }{\int }{Y_{\ell }^m}^{\ast }(\theta ,\varphi )f(\theta ,\varphi )\mathrm{d}\Omega ,$$

$$S_{\ell }:=c_{\ell }^2=\sum _{m=-\ell }^{\ell }{\left|\underset{\Omega }{\int }{Y}_{\ell }^m(\theta ,\varphi )f(\theta ,\varphi )\mathrm{d}\Omega \right|}^2.$$

Full text

stars:massive – convection – hydrodynamics – turbulence – supernovae: general

††journal: ApJ††software: NumPy (van der Walt et al., 2011), Scipy (Oliphant, 2007), IPython (Pérez & Granger, 2007), Matplotlib (Hunter, 2007), Visit (Childs et al., 2012), and SHTns  library (Schaeffer, 2013).

1 Introduction

Multi-dimensional effects in the late burning stages of massive stars have recently garnered considerable interests for a number of reasons. Whereas diffusive mixing models based on the mixing-length theory (MLT) of convection (Biermann, 1932; Böhm-Vitense, 1958) as used in one-dimensional (1D) stellar evolution models provide a reasonable estimate of mixing (Arnett et al., 2019, for a more quantitative discussion) in the interior of convective zones for the late, neutrino-cooled burning stages (Müller et al., 2016b, see Figure 10 & 14)111The upper panel of Figure 10 shows a comparison between the radial rms velocity fluctuations $\delta v_{\rm r}$ (black) in the 3D model to the convective velocity $v_{\rm conv}$ computed in the Kepler model using MLT (red). Inside the convection zone ($1.75\,\mathord{-}\,2.25\,\rm M_{\odot}$) the 1D and 3D velocities are in close agreement. Figure 14 shows a comparison between the mass fraction profiles in the 3D model and the 1D model. Although, the differences in the mass fractions in the interior of the oxygen shell are minute, the differences in gradients at shell boundaries are conspicuous., it has long been speculated that additional phenomena such as turbulent entrainment (Fernando, 1991; Strang & Fernando, 2001; Meakin & Arnett, 2007a; Spruit, 2015) and the excitation of internal waves could significantly affect shell growth, mixing, and angular momentum transport (Cantiello et al., 2014; Fuller et al., 2015) in a manner that is not captured by current 1D stellar evolution models. Convective-reactive systems (Dimotakis, 2005) are at the extreme end of mixing problems, where the feedback produced by nuclear burning strongly affects the flow and vice-versa, thereby making it extremely hard to describe chemical mixing in these systems using simplified 1D prescriptions.

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 & Arnett, 1994, 1998; Asida & Arnett, 2000) and three dimensions (3D; Kuhlen et al., 2003; Meakin & Arnett, 2006, 2007a, 2007b), additional convective boundary mixing has indeed been consistently observed in many modern simulations (Müller et al., 2016b; Cristini et al., 2017; Jones et al., 2017; Andrassy et al., 2018) 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 (Müller, 2016).

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 et al., 2011; Stancliffe et al., 2011; Herwig et al., 2014; Woodward et al., 2015) 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 et al., 2011; Jones et al., 2016; Ritter et al., 2018).

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 & Woosley, 2014, see, e.g., the $15\,\text{M}_{\odot}$ 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 et al. (2018), who found such late mergers in about $40$ percent of all progenitors between $16\,\text{M}_{\odot}$ and $26\,\text{M}_{\odot}$. 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 et al., 2010; Isensee et al., 2010; Grefenstette et al., 2014, 2017). 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 (Müller & Janka, 2015; Müller et al., 2017; Couch et al., 2015) and on slow, steady-state entrainment (Meakin & Arnett, 2007a, b; Müller et al., 2016b; Jones et al., 2017; Cristini et al., 2017; Andrassy et al., 2018), but such dynamical shell mergers have not yet been investigated thoroughly. Müller (2016) reported a late breakout of convection across shell boundaries in a thin O burning shell that was reignited shortly before collapse in an $12.5\,\text{M}_{\odot}$ progenitor, but with insufficient time before collapse for a genuine merger of the O and Ne shell. More recently, Mocák et al. (2018) observed the merging of shells in a 3D simulation of a $23\,\text{M}_{\odot}$ star, but their model was restricted to a small wedge of $27.5^{\circ}\,\mathord{\times}\,27.5^{\circ}$ 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 $4\pi$ 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 $18.88\,\text{M}_{\odot}$ progenitor (Müller et al., 2016a) 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 $18.88\,\text{M}_{\odot}$. The 1D model has been calculated using the Kepler stellar evolution code (Heger & Woosley, 2010) as outlined in Müller et al. (2016b, their Section 2.1)222 Kepler uses the Ledoux criterion for semiconvection as described in Weaver et al. (1978). We use $\alpha\,\mathord{=}\,1$ for MLT (larger [$\mathord{\approx}\,1.7$] and depth-dependent MLT values used in the solar convection zone may only be valid for the surface convection zone of the sun). Convective boundaries are ‘softened’ by one formal overshoot zone with $1/10$ mixing efficiency as semiconvection (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 or consensus about their validity in all evolution stages. . The 1D model is mapped to the 3D grid $420\,\text{s}$ before the onset of core collapse. Figure 1 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 $1.7\mathord{-}4.5\,\text{M}_{\odot}$ ($\mathord{\approx}\,2,000\mathord{-}42,000\,\text{km}$). 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 $2.5\mathord{-}3.7\,\text{M}_{\odot}$, which is separated by a non-convective layer from a convective C burning shell with ashes of O and Ne between $4.1\mathord{-}4.5\,\text{M}_{\odot}$. The O and Ne shells are separated by a thin interface at $\mathord{\approx}\,2.45\,\text{M}_{\odot}$. The Si shell and the Fe core inside $1.6\,\text{M}_{\odot}$ are practically inert, with neutrino cooling dominating over burning. The middle panel shows the profiles of the radial velocity, $v$, and the convective velocity, $v_{\text{conv}}$, calculated using the same choice of dimensionless coefficients as (Müller, 2016). The convective velocity near the base of the O shell reaches up to $\mathord{\sim}\,100\,\text{km}\,\text{s}^{-1}$. The Ne and C shells are relatively quiet, having convective velocities close to $20\,\text{km}\,\text{s}^{-1}$ and $10\,\text{km}\,\text{s}^{-1}$, 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 $\mathord{\approx}\,2.0\,k_{\text{b}}/\text{nucleon}$ and $\mathord{\approx}\,4.0\,k_{\text{b}}/\text{nucleon}$ (here $k_{\text{b}}$ 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.

To be clear, choosing “$18.88\,\text{M}_{\odot}$” (accurate to two decimal places) does not amount to “fine-tuning” of the model presented in Müller et al. (2016b). We would like to emphasize that the stellar model presented in this work and the stellar model used by Müller et al. (2016b) are $\mathord{\approx}\,1\rm\,\text{M}_{\odot}$ apart, and it is well known that the pre-supernova structure of a star is a very sensitive function of the zero-age main sequence mass (ZAMS) (Sukhbold et al., 2018, see upper left panel of their Figure 8). Collins et al. (2018) evaluated convective velocities and eddy scales in the oxygen and silicon burning shells of a large set (2353 solar-metallicity, non-rotating, calculated using the Kepler code) of progenitor models with ZAMS between $9.45$ and $35\rm\text{M}_{\odot}$. One of their key finding is that shell mergers are not rare; “about 40 percent of progenitors between 16 and 26 $\rm\text{M}_{\odot}$ exhibit simultaneous oxygen and neon burning in the same convection zone as a result of a shell merger shortly before collapse”. At the same time, there is no clear mapping between the outcomes in 1D and 3D, and therefore, it cannot be said a priori with certainty if a certain 1D model will exhibit a shell merger in 3D simulation. In a nutshell, the $18.88\,\text{M}_{\odot}$ model was chosen for this study as the dominant angular wave numbers ($\ell$) of convection in Si and O shell were small, and at the same time the convective Mach number was large.

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 & Sato, 2004; Wongwathanarat et al., 2010) as implemented in Melson et al. (2015) with a uniform angular resolution ($2^{\circ}$) in the $\theta$ and $\phi$ directions444Adequacy of $2^{\circ}$ angular resolution: A simulation using an equidistant spherical grid ($\Delta\theta\,\mathord{=}\,\Delta\phi\,\mathord{=}\,\alpha$) captures angular modes with $\ell$ ranging between zero and $180^{\circ}/\alpha$. Therefore, in our case with $\alpha\,\mathord{=}\,2^{\circ}$, we can extract modes up to $\ell\,\mathord{=}\,90$. The well developed inertial range (see left panel of Figure 15 in Section 5) is a proof that the simulation is in the regime of implicit large eddy simulations (Boris et al., 1992, ILES), and the energy transfer between different length scales is correctly modelled. Moreover, ignoring the smallest scale does not affect the evolution, as the mixing (and the related hot-spot burning as described in Section 4.3) is caused by the large-scale convective flow.. For the purpose of this simulation we changed the $4^{\mathrm{th}}$-order reconstruction originally used in the piecewise-parabolic method (PPM) of Colella & Woodward (1984) to a $6^{\mathrm{th}}$-order extremum-preserving reconstruction (Colella & Sekora, 2008; Sekora & Colella, 2009). We use the Helmholtz equation of state (EoS) described in Timmes & Arnett (1999) and Timmes & Swesty (2000), which is thermodynamically consistent to high accuracy. Nuclear burning is treated using a 19-species $\alpha$-network (Weaver et al., 1978). 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 et al., 1996) is included as a local sink term.

3.1 Boundary Conditions

We excise the core inside $1.7\,\text{M}_{\odot}$ 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 $5.0\,\text{M}_{\odot}$ 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 $0.1\,\text{M}_{\odot}$ away from the inner grid boundary (see Figure 1). 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. Hydrodynamic and thermodynamic variables are mapped using second order polynomial interpolation. Mass fractions are mapped using linear interpolation (to preserve the discontinuities at shell boundaries). The question of initial transients (mapping errors) is discussed in Section 4.5. In order to break the spherical symmetry we introduce small deviations from hydrostatic equilibrium by perturbing the density field such that $\delta\rho=U[-1,1]f\rho$, where $U[-1,1]$ is a uniformly distributed random number in the interval $[-1,1]$ and the amplitude $f=2.0\times 10^{-4}$.

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 1. We have not done a convergence study for the model presented in the paper. Please refer to the discussion on the effect of angular and radial resolution on a simulation of oxygen burning in case of the $18\,\rm\text{M}_{\odot}$ model by Müller et al. (2016b, see Appendix B). The simulations are converged in terms of the coupling between nuclear energy generation rate and growth of turbulent motions (Müller et al., 2016b, Figure 18); and dominant multipoles (Müller et al., 2016b, $\ell$-modes in Figure 19) emerging in the flow.

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

[TABLE]

where $g\,\mathord{=}\,-\mathrm{d}\Phi/\mathrm{d}r\,\mathord{<}\,0$ is the gravitational acceleration and

[TABLE]

is the adiabatic index. For unstable modes $\omega_{\text{BV}}\,\mathord{>}\,0$ is the growth rate; for stable modes $\omega_{\text{BV}}$ is imaginary and $|\omega_{\text{BV}}|$ is the oscillation frequency. Following Buras et al. (2006), we redefine the Brunt-Väisälä frequency such that

[TABLE]

which is more convenient for visualization purposes. The conversion of Equation (1) to Equation (4) is exemplified in Appendix A of Müller et al. (2016b)555Heger et al. (2005) show an alternative approach (used in Kepler stellar evolution code) which can be used to extract the contribution of composition to $\omega_{\text{BV}}$ without explicitly tracking composition.. Note that our definition changes only affects the sign convention, but not the absolute value of $\omega_{\text{BV}}$. For unstable modes ($C_{\mathrm{L}}\,\mathord{>}\,0$), $\omega_{\text{BV}}$ is the growth rate, and for stable modes ($C_{\mathrm{L}}\,\mathord{<}\,0$), $\omega_{\text{BV}}$ is the negative of the oscillation frequency. Therefore, according to the sign convention adopted in this paper $\omega_{\text{BV}}\,\mathord{>}\,0$ corresponds to convective instability.

Figure 2 (left panel) shows the spherically averaged Brunt-Väisälä frequency $\overline{\omega}_{\text{BV}}$ profiles as a function of time for the 3D model. $\overline{\omega}_{\text{BV}}$ is defined as

[TABLE]

which is different from $\omega_{\text{BV}}$ calculated using spherically averaged $\rho$, $p$ etc. The plot shows three convectively stable layers (zones I$-$III at masses of $2.4\mathord{-}2.5\,\text{M}_{\odot}$, $3.7\mathord{-}4.1\,\text{M}_{\odot}$, and $4.4\mathord{-}4.5\,\text{M}_{\odot}$ respectively) sandwiched between more voluminous convectively unstable regions. The location ($m$ and $r$) and width ($\Delta m$ and $\Delta r$) of these zones at $t=0$ are listed in Table 2 in order of increasing mass coordinate. The initial and final entropy profiles (green and red lines in left panel of Figure 2) 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 (see right panel of Figure 10), the stabilizing gradient in Zone I gradually becomes weaker, and disappears altogether after $\mathord{\approx}\,250\,\text{s}$, allowing the O and Ne shell to merge into one large convection zone extending between $1.7\mathord{-}3.6\,\text{M}_{\odot}$. The disappearance of Zone I is mirrored by the disappearance of the entropy jump at $2.4\,\text{M}_{\odot}$. 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 $\mathord{\approx}\,350\,\text{s}$, 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 2 shows $\omega_{\text{BV}}$ 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 $\mathord{\approx}\,60\,\text{s}$ there is a reduction in the strength of Zone I. The contraction of the Si core after the end of the last Si shell burning phase, and deleptonisation of the Fe core leading to contraction, drives more burning (at the base of) the O shell. This raises the entropy and weakens the ‘buoyancy gap’ between the two shells (in 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 $\rho^{\prime}$ which are defined as

[TABLE]

where $\overline{\rho}$ is the average density evaluated over a spherical shell. Figure 3 shows the density fluctuations on an $x\mathord{-}y$ slice at different times. The panel at $35\,\text{s}$ 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 $\mathord{\approx}\,8000\,\text{km}$ 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 $35\,\text{s}$ and $175\,\text{s}$. 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 & Fernando, 2001, see Figure 1 and 3 there). As a result of the entrainment and mixing the stabilizing gradient in Zone I ceases to exist around $250\,\text{s}$ (see Figure 2 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 $\mathord{\approx}\,20,000\,\text{km}$, At this point, global asymmetries develop in the flow. The panel at $323\,\text{s}$ in Figure 3 shows an asymmetric feature with a size of $\mathord{\approx}\,20,000\,\text{km}$, which characterizes the phase of vigorous convective activity. Close to the end of simulation (panels at $375\,\text{s}$ and $420\,\text{s}$) the density fluctuations in the region within $30,000\,\text{km}$ reach a level of $\mathord{\approx}\,5\mathord{-}10\%$ and are again of smaller scale structure, while in the region between $\mathord{\approx}\,30,000$ and $40,000\,\text{km}$ they are as large as $\mathord{\approx}\,30\mathord{-}40\%$ with a marked dipolar asymmetry.

The cautious reader may have also noticed the sharp discontinuities in all the panels, and the associated “waves” in the panels from $323\,\text{s}$ onward. The reason is that the fluctuation field renders the convective boundaries sharp; the convective flow rapidly decelerates due to buoyancy braking over a short length scale ($\mathord{\sim}\,0.1H_{P}$, where $H_{P}$ is the local pressure scale height). The flow deposits energy and momentum and produces density fluctuations which mark the convective boundaries. The convective boundaries are stiff; or in other words have a large bulk Richardson Number (Turner, 1979, see Chapter 10), and act like a “stretched membrane”. These boundaries reflect/transmit the internal waves generated by turbulent convection, and penetrative convection. On comparing the panels at $323\,\text{s}$ and $375\,\text{s}$, we see that the waves are absent outside $38,000\,\text{km}$, and are present in the convectively active region only. This indicates that the waves are not spurious and are certainly linked to the convective activity.

Figure 4 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 $m$ 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 (left) shows that the inner convective boundary moves towards lower $m$ 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 $1.75\,\text{M}_{\odot}$.

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)

[TABLE]

where $\tilde{v}_{r}$ is the Favre-averaged radial velocity. The fluctuating component $v_{r}^{\prime\prime}$ of the radial velocity is therefore given by

[TABLE]

Figure 5 shows the radial velocity fluctuations corresponding to the snapshots shown in Figure 3. The panel at $35\,\text{s}$ shows well-developed plumes extending from the O burning region at $\mathord{\approx}\,3,000\,\text{km}$ to the base of Zone I at $8,000\,\text{km}$. The panel at $175\,\text{s}$ shows the secondary plumes extending from $\mathord{\approx}\,8,000\,\text{km}$ to $\mathord{\approx}\,20,000\,\text{km}$. 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 $2\mathord{-}3$) 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 $268\,\text{s}$ which is the stage when the stable Zone I has already disappeared.The panel at $323\,\text{s}$ 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 $\mathord{\approx}\,800\,\text{km}\,\text{s}^{-1}$ to more than $1,600\,\text{km}\,\text{s}^{-1}$ (panels at $375\,\text{s}$ and $420\,\text{s}$) within the next $120\,\text{s}$. The panel at $420\,\text{s}$ shows the presence of a large-scale dipolar asymmetry in the flow. Convective downdrafts in the innermost region reach velocities in excess of $1,600\,\text{km}\,\text{s}^{-1}$ aided by rapid contraction 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. In the absence of a buoyancy barrier convective flow naturally grows to the largest possible unimpeded length scale. The longitudinal two-point correlation function gives a qualitative measure of the size of the convective region. It is defined as (Equation 15 of Müller et al., 2016b)

[TABLE]

where we have used the Reynolds-averaged velocities (see Appendix A). The correlation length $\Lambda_{\text{corr}}$ at a radius $r$ is defined by integrating the correlation function

[TABLE]

between $r_{-}$ and $r_{+}$, which are usually taken to be $-\infty$ and $+\infty$ respectively. Meakin & Arnett (2007b, Appendix B) define these points such that $C(r,r_{\pm}-r)\,\mathord{\geq}\,0.5$. In the present case, estimating the correlation length using the equation above is non-trivial because of the presence of multiple convectively stable zones. For this analysis we divide the simulation domain into multiple layers according to their initial convective stability and the burning processes:

1. a.

$1.7\mathord{-}1.8\,\text{M}_{\odot}$, convectively inert Si shell, 2. b.

$1.8\mathord{-}2.4\,\text{M}_{\odot}$, convective O shell, 3. c.

$2.4\mathord{-}2.5\,\text{M}_{\odot}$, Zone I, 4. d.

$2.5\mathord{-}3.6\,\text{M}_{\odot}$, convective Ne shell, 5. e.

$3.6\mathord{-}4.1\,\text{M}_{\odot}$, Zone II, 6. f.

$4.1\mathord{-}4.4\,\text{M}_{\odot}$, convective C shell, 7. g.

$4.4\mathord{-}4.5\,\text{M}_{\odot}$, Zone III, 8. h.

$4.5\mathord{-}5.0\,\text{M}_{\odot}$, He shell.

Figure 6 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 $C(r,\delta r)$ drops rapidly from unity with $\delta r$, resulting in $\Lambda_{\text{corr}}\,\mathord{<}\,1,000\,\text{km}$ at all times ($\Lambda_{\text{corr}}/r\,\mathord{\ll}\,1$). 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 ($\Lambda_{\text{corr}}\,\mathord{\lesssim}\,r$) in panels b, c, and d such that its value approaches zero for $\delta r\,\mathord{\approx}\,10,000\,\text{km}$. 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 $\mathord{\approx}\,180\,\text{s}$ there is a well defined convective region centered somewhere between $3.0\mathord{-}3.8\,\text{M}_{\odot}$. As we traverse Zone III and go into the He shell, the correlation length becomes large ($\Lambda_{\text{corr}}\mathord{\sim}r$) and the function approaches a symmetrical form. This suggests that after $\mathord{\approx}\,300\,\text{s}$ there is a well defined convection layer centered somewhere between $4.5\mathord{-}4.8\,\text{M}_{\odot}$. 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 7 (left panel) shows the net energy generation rate in the 1D model as a function of time. In the region inside $0.6\,\text{M}_{\odot}$ (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 $1.8\mathord{-}2.1\,\text{M}_{\odot}$. In addition, there are contributions from Ne shell burning, C shell burning, and He shell burning (not shown in the figure). Note that the C shell burning and the He shell burning are not relevant for the dynamics. Also, the total energy generation in the He shell is negligible. No energy from it enters or leaves the simulation domain, as its main purpose is to provide ‘boundary pressure’.

Figure 7 (right panel) shows the net energy generation rate ($|\dot{q}_{\text{nuc}}|\,\mathord{\gg}|\,\dot{q}_{\nu}|$ everywhere) in the 3D model. The evolution between $2.3\mathord{-}4.5\,\text{M}_{\odot}$ (upper half) does not show a lot of deviation from its initial behaviour, except that close to $\mathord{\approx}\,300\,\text{s}$ the Ne burning shell moves inwards in mass (in a Lagrangian sense). The lower half shows a zoom of the region between $1.8\mathord{-}2.3\,\text{M}_{\odot}$. Most of the energy production happens at the base of the O shell. Around $100\,\text{s}$ 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 $300\,\text{s}$ there is a short but intense Ne burning episode. This phenomenon is marked as episodic burning, which is reminiscent of hot-spot burning 666Note that the term “hot-spot burning” is being used here in a figurative sense. 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 at the higher ambient density and temperature. We think that our usage of the term ”hot-spot burning” confirms with the phenomenon described by Bazan & Arnett (1994, see their Figure 3). Quoting Bazan & Arnett (1998), “Significant mixing beyond convective boundaries determined by mixing-length theory brings fuel (${}^{12}\rm C$) into the convective region, causing hot spots of nuclear burning.” seen by Bazan & Arnett (1994, see Fig. 3) in their 2D study of O burning in a $20\,\text{M}_{\odot}$ progenitor (Arnett, 1994).

Figure 8 shows the resulting shell merger marked by Si rich outflow rapidly moving into the Ne shell. The isosurfaces shown track the points corresponding to $250\,\text{km}\,\text{s}^{-1}$ radial velocity. The Si-rich material moves rapidly outwards carving out an elongated cavity in the enclosing Ne shell. The velocity isosurface expands from $\mathord{\approx}\,10,000\,\text{km}$ to more than $20,000\,\text{km}$ in a short span of $\mathord{\approx}\,160\,\text{s}$. The episodic nuclear burning peaks at around $350\,\text{s}$ (as shown in Figure 7 by $\log\dot{q}$ contours) and triggers/powers a phase of violent convective activity. A comparison between the Ne abundance at $360\,\text{s}$ and $420\,\text{s}$ shows that the total amount of Ne decreases by $\Delta M_{\text{Ne}}\mathord{\approx}\,0.2\,\text{M}_{\odot}$ ($\mathord{\approx}$ 70 percent) 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 ${}^{20}\text{Ne}$ nuclei to ${}^{16}\text{O}$ and ${}^{24}\text{Mg}$ nuclei (Iliadis, 2015, Equation 5.108)

[TABLE]

with the average energy production (per $40$ nucleons) of $\overline{Q}_{\text{Ne}}\mathord{\approx}\,6.2\,\text{MeV}$ ($\mathord{\equiv}\,1.5\mathord{\times}10^{17}\,\text{erg}\,\text{g}^{-1}$) near $T\,\mathord{\approx}\,1.5\,\text{GK}$ (Iliadis, 2015, Equation 5.109). This implies a total release of $\overline{Q}_{\text{Ne}}\Delta M_{\text{Ne}}\,\mathord{\approx}\,6.0\mathord{\times}10^{49}\,\text{erg}$ of thermal energy inside $\mathord{\approx}\,2.5\,\text{M}_{\odot}$ of stellar plasma within $\mathord{\approx}\,100\,\text{s}$, which further translates into $\mathord{\approx}\,4\mathord{\times}10^{15}\,\text{erg}\,\text{g}^{-1}$ in each component ($E_{\text{r}}$, $E_{\theta}$ and $E_{\phi}$) of kinetic energy.

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

[TABLE]

where $m_{i,-},m_{i,+},\delta m_{i}$ are the inner boundary, outer boundary and mass, respectively, of the $i$th layer (defined in Section 4.2.2). Kinetic energy starts building up with the onset of convection and reaches a plateau phase by $\mathord{\approx}\,70\mathord{-}80\,\text{s}$ except for the Ne shell. The stationary phase lasts till $\mathord{\approx}\,320\,\text{s}$, when episodic Ne burning begins. In the next $\mathord{\approx}\,20\mathord{-}30\,\text{s}$ there is a rapid rise in kinetic energy by a factor of $\mathord{\sim}10$. In the end the growth saturates, and the entire merged shell has uniform specific kinetic energy. The lower panel of Figure 9 shows the ratio $\psi$ of kinetic energy in the fluctuating radial component and the transverse components. The transverse components of the specific kinetic energy are defined as

[TABLE]

The ratio fluctuates with time and hardly stays close to the equipartition value ($E_{r}=E_{\theta}=E_{\phi}$) for the convectively active shells. For the stable Zone I the value deviates away from $1/2$ when the shell is finally eroded away. Interestingly, transverse motions dominate over radial motion in Zone II after $\mathord{\approx}\,250\text{s}$.

4.3.2 Entropy: 1D vs 3D

Regions 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 10 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 $1.8\mathord{-}2.4\,\text{M}_{\odot}$ increases gradually until the first $360\,\text{s}$ ($\Delta s\,\mathord{\approx}\,0.1\,k_{\text{b}}/\text{nucleon}$), keeping the gradient unchanged. In the last $60\rm\,\text{s}$, the inner Fe/Si core experiences significant contraction ($\approx$50%). The associated density and temperature increment leads to a rapid increase in the burning rate of oxygen as well as the neon (1D mixing). This raises the overall entropy in the oxygen shell and weakens the ‘buoyancy gap’ between the O and Ne shells. The softening of convective boundary 777Note: In Kepler , convective boundaries are ’softened’ by one formal overshoot zone with $1/10$ mixing efficiency as semiconvection. leads to more efficient mixing.

In case of the 3D model the entropy evolution is remarkably different. The overall entropy increase is mostly due to the heat deposited by oxygen burning and hot-spot burning. In the first $360\,\rm\text{s}$, the average entropy for the matter inside $2.5\,\mathord{-}\,3.7\,\rm\text{M}_{\odot}$ decreases by $\Delta s_{2.5-3.7}\,\mathord{\approx}\,0.03\,k_{\text{b}}/\text{nucleon}$. This high-entropy material is carried inwards by convective downdrafts. During the same time, the average entropy for the matter inside $1.8\,\mathord{-}\,2.5\,\rm\text{M}_{\odot}$ increases by $\Delta s_{1.8-2.5}\,\mathord{\approx}\,0.54\,k_{\text{b}}/\text{nucleon}$. The ratio of $\Delta s_{1.8-2.5}$ and $\Delta s_{2.5-3.7}$ is $\mathord{>}\,10$, which clearly shows that the entropy gain of the inner region cannot be explained as resulting from entrainment of high-entropy material and must be a result of nuclear burning. The entropy profile inside $1.8\mathord{-}2.4\,\text{M}_{\odot}$ is levelled in the first $60\,\text{s}$, but the additional Ne burning preserves the small negative entropy gradient between $2.4\mathord{-}3.6\,\text{M}_{\odot}$. At the same time the steep entropy barrier close to $2.4\,\text{M}_{\odot}$ is softened (the entropy jump is reduced). With the complete disappearance of the stabilizing gradient provided by Zone I at $\mathord{\approx}\,250\,\text{s}$, thermal energy can be convected into a much larger volume. We define three timescales relevant for the growth of entropy: the convective timescale ($t_{\rm conv}$), the nuclear timescale ($t_{\rm nuc}$), and the expansion timescale ($t_{\exp}$) which are defined as:

[TABLE]

where $H_{P}\,\mathord{=}\,P\mathrm{d}r/\mathrm{d}P$ is the pressure scale height, and $e_{\rm th}$ is the specific internal energy. Figure 11 shows a comparison between the nuclear and convective timescales, $\xi_{1}\,\mathord{\coloneqq}\,t_{\rm nuc}/t_{\rm conv}$, in the region inside $2.5\,\rm\text{M}_{\odot}$ during the last $60\,\text{s}$ before the onset of collapse. The value of $\xi_{1}$ is much larger than unity everywhere (except close to the inner boundary), which means convection can carry away all the nuclear energy released. However, this can not explain the entropy peak seen at the base of the O shell (see the right panel of Figure 10). Now, let us take the contraction of the inner core into account. During the last $60\,\text{s}$, the inner boundary of the O shell recedes inwards rapidly, while at the same time the outer layers move inward at lower velocities. Therefore, a given mass shell located between $r$ and $r\mathord{+}\delta r$ expands in size. A convective plume which starts at the base of such a shell needs to cover a longer distance to leave the shell. This would result in trapping of the convective flow leading to a reduction in the convective efficiency. Figure 11 shows a comparison of the expansion and convective timescales, $\xi_{2}\,\mathord{\coloneqq}\,t_{\rm exp}/t_{\rm conv}$. During the last seconds before collapse this ratio approaches unity, which implies that the convective flow is trapped. The inability of convection to transfer the released nuclear energy leads to a rapid buildup of entropy ($\Delta s\,\mathord{\approx}\,0.2\,k_{\text{b}}/\text{nucleon}$ during the last $5\,\rm\text{s}$) at the base of the O shell, leading to an overall negative entropy gradient extending up to $\mathord{\approx}\,3.6\,\text{M}_{\odot}$, 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 12 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 $\mathord{\approx}\,250\,\text{s}$ (panels at $268\,\text{s}$) show significant deviations from spherical symmetry. The dipolar feature (panel at $323\,\text{s}$) in the Ne distribution is a cavity created by a Si-rich bubble (panel at $323\,\text{s}$). The downdrafts of Ne-rich material propagating into the Si-enriched inner volume (panel at $323\,\text{s}$) 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 $420\,\text{s}$. The Ne mass fraction decreases considerably ($\mathord{\approx}\,0.2\,\text{M}_{\odot}$ of Ne is burnt in total) after the merger and the leftover Ne is thoroughly mixed with Si. Figure 13 shows an overview of the merger and expansion of the Si and Ne shells. We also show the evolution of the typical angular mode number characterizing convection in the merging layers evaluated using the following relation (Chandrasekhar, 1981)

[TABLE]

where $r_{+}$, $r_{-}$ are the inner and outer radii of the convective shell. The value of $\ell$ decreases rapidly as the ongoing shell merger leads to the formation of a radially extended convective region reaching from the base of the O shell to the outer edge of the convective Ne shell.

Figure 14 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 $\mathord{\approx}\,2.7\,\text{M}_{\odot}$. Abundance changes in the 1D model take place during the last $\mathord{\approx}\,60\,\text{s}$, 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 Müller et al. 2016b, Section 2.1). Specifically, during the last $60\,\text{s}$, the core experiences significant contraction ($\approx$50%). The associated density and temperature increment leads to an enhancement in the neon burning rate, which in turn leads to an increase in the convective velocity resulting in faster mixing (according to the MLT model). 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 $2.5\,\text{M}_{\odot}$ increases during the last $60\,\text{s}$. The increase results from a combination of large-scale mixing (convective downdrafts carrying O-rich material inwards, $\mathord{\approx}\,0.05\,\text{M}_{\odot}$) and rapid Ne burning (producing O via ${}^{20}\text{Ne}(\gamma,\alpha)^{16}\text{O}$), and negates the change caused by burning of O to Si. In the end, the total mass of O in the simulated region does not change.

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 $4.0\mathord{-}4.5\,\text{M}_{\odot}$ remains unaffected in both cases. In the 3D model, the gradual rise in the Ne mass fraction inside $2.4\,\text{M}_{\odot}$ continues till $\mathord{\approx}\,360\,\text{s}$ 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, the evolution of silicon mass fraction is a result of two processes: 1) silicon production by oxygen burning via ${}^{16}\text{O}(^{16}\text{O},\alpha)^{28}\text{Si}$, and 2) outward mixing caused by convection. During the course of evolution $\mathord{\approx}\,0.15\,\text{M}_{\odot}$ of Si is produced and $\mathord{\approx}\,0.30\,\text{M}_{\odot}$ of Si is gradually transported out from the region between $1.8\mathord{-}2.5\,\text{M}_{\odot}$ into the region between $2.5\mathord{-}4.0\,\text{M}_{\odot}$. As a result the outer boundary of the Si-enriched layers from $\mathord{\approx}\,10,000\,\text{km}$ to $\mathord{\approx}\,35,000\,\text{km}$ (close to the C/He interface).

4.5 Initial Transients and Relaxation

Mapping data from one code to another invariably leads to transients. At this point, we must delineate transients caused by mapping errors (spurious), and the transients caused by rapid changes in the thermodynamical properties of a system (genuine).

First, we provide two arguments why our results are not affected by mapping errors:

Convective Turnover Time: In order to smooth (damp) mapping related errors, one aims at simulating a sufficiently large number of convective turnovers of the flow. This approach to eliminate the effects of initial transients (spurious) in a simulation is useful when convection is supposed to lead to a steady-state. The convective turnover timescale is conventionally defined as

[TABLE]

where $H_{P}$ is the local pressure scale height, $v_{r}^{\prime\prime}$ is the fluctuating component of the radial velocity. Figure 16 shows the value of the convective turnover timescale averaged over a period of $10\,\text{s}$ at various times. The timescale evolves rapidly over the course of the simulation. In the innermost convectively active region, which is the Si-rich part of the O shell (see Figure 1), the convective turnover timescale is between $\mathord{\approx}\,10\mathord{-}20\,\text{s}$ during the entire simulation, which is sufficient for $20\mathord{-}40$ full turnovers during the simulation. In the Ne shell, the turnover timescale becomes relatively stable at around $230\,\text{s}$ and its value is $\mathord{\lesssim}\,80\,\text{s}$ thereafter. Therefore, even after the onset of the shell merger, which begins at $\mathord{\approx}\,250\,\text{s}$ (see section 4), the region has enough time to go through $2\mathord{-}3$ convective turnovers. Furthermore, the convectively stable layers sandwiched between the convectively unstable layers (see Section 4.1) play an effective role in absorbing and dissipating the initial transients because they damp the propagation of acoustic waves. 2. 2.

Growth of Kinetic Energy In the Merging Shells: During the onset phase, between $0\,\mathord{-}\,80\,\text{s}$ (Figure 9), the kinetic energy in radial velocity fluctuations grows and reaches a steady value in the layers of interest, which are the Si-rich part of the O shell (labelled as “convective O shell”) and the convective Ne shell, implying that most of the initial transients (mapping errors) have already dissipated away and do not affect the further evolution.

In the case of a genuine transient, the system may undergo an irreversible change in a short time, and may never approach a steady-state, which is also the case with the shell-merger described above. In this case, the numerical solution is stable (the CFL value is 0.6; the code preserves the hydrostatic state; PPM is well understood), but the physical system itself is on the verge of a thermodynamical instability. The fact that the system under consideration here manifests such a behaviour is evident from the nature of the transition between the ‘stationary’ phase and the ‘rapid rise’ phase (see Figure 9). Before this transition, the kinetic energy in the convective O shell and Zone-I stays at a plateau between $\mathord{\approx}100\,\mathord{-}\,300\,\text{s}$, and in the Ne-shell energy is slowly deposited by convective activity underneath and gently reaches a steady state value at $\mathord{\approx}\,200\,\text{s}$. For the next $\mathord{\approx}\,125\,\text{s}$ the change in kinetic energy in the Ne shell is negligible. Its value starts rising again only at the time of the transition in response to the elevated energy production in the O burning shell, associated with the shell merger beginning at $\mathord{\approx}\,250\,\text{s}$. Therefore, the ensuing shell merger represents a “genuinely transient phenomenon” consistent with the thermodynamical evolution. It is obvious that the conditions emerging from the shell merger do not represent a steady state situation, and we cannot claim that the same conditions or the same dynamical situation would be present at the onset of collapse if we had started our simulation significantly earlier. However, we consider the shell merger starting $\mathord{\approx}\,200\,\text{s}$ before collapse as an exemplary case of such a possibility, which could happen in other stars even if not in the one considered here. In Section 6, we further discuss the effects of our choice of the starting time on the onset and dynamics of the shell merger and its impact on the collapse conditions.

5 Pre-Supernova Model Properties

Large-scale modes with $\ell\mathord{\sim}1\mathord{-}2$ are most effective for shock revival (Müller & Janka, 2015; Abdikamalov et al., 2016; Müller et al., 2016a). In order to gain a quantitative understanding of the modes, we have done spherical harmonic decomposition of the density and velocity perturbations. Figure 15 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 ($\ell_{\text{max}}$) at a given mass coordinate. The value of $\ell_{\text{max}}$ changes with radius, but there are continuous sections where it remains stationary. In the case of density fluctuations inside the merged O-Ne shell ($1.8\mathord{-}3.1\,\text{M}_{\odot}$) a mixture of dipole ($\ell\,\mathord{=}\,1$) and quadrupole mode ($\ell\,\mathord{=}\,2$) dominates, and in the region outside of Zone II the dipole mode dominates (similar to the global non-spherical oscillation described by Herwig et al. (2014)). This confirms our visual impression obtained from Figure 3. In the case of radial velocity fluctuations a mixture of dipole mode and $\ell=3$ mode dominates in the merged O-Ne shell, whereas the dipole mode dominates in the C and He shells. The normalized probability distribution of $\ell_{\text{max}}$ (shown as histograms) confirms that $\ell\,\mathord{\sim}\,1\mathord{-}3$ dominate both the density fluctuations and the velocity fluctuations. Figure 17, 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 $\mathord{\approx}\,3.6\,\text{M}_{\odot}$ (layers c and d) are shallow compared to the region outside and including Zone II. In the upper half-panel the spectra for $\ell\,\mathord{\lesssim}\,15$ are shallower in the merged O-Ne shell including Zone II. The spectra above $\ell\,\mathord{\gtrsim}\,15$ agree quite well with the ‘Kolmogorov slope’.

Figure 17 (right panel) shows the profiles of the turbulent Mach number corresponding to the radial velocity fluctuations at multiple times. The quantity is defined following Müller et al. (2016b) as

[TABLE]

where $c_{\text{s}}$ is the adiabatic sound speed. The Mach number in the O and Ne shells increases gradually until the shell merger. During the first $120\,\text{s}$, the value of $\mathcal{M}_{r}$ is almost negligible in the Ne shell outside $2.5\,\text{M}_{\odot}$, which is due to the smaller Brunt-Väisälä in the Ne shell. The region outside $4.1\,\text{M}_{\odot}$ stays quiescent until $360\,\text{s}$ as the stable Zone II located between $3.6\mathord{-}4.1\,\text{M}_{\odot}$ does not allow convective plumes to penetrate further. After the merger, between $300\mathord{-}360\,\text{s}$, $\mathcal{M}_{r}$ increases rapidly by a factor of $\mathord{3}$ reaching a peak value of $\mathord{\approx}\,0.14$ 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 $4.5\,\text{M}_{\odot}$.

Figure 18 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 $X_{\text{Si}}$, 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 $\mathord{\approx}\,100\,\text{km}\,\text{s}^{-1}$ somewhere around $\mathord{\approx}\,30,000\mathord{-}40,000\,\text{km}$. In the panel for $X_{\text{Si}}\,\mathord{=}\,0.2$ (middle), we see a highly asymmetrical and clumpy distribution of Si with some of the clumps moving outwards at velocities in excess of $1,000\,\text{km}\,\text{s}^{-1}$. 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 $X_{\text{Si}}\,\mathord{=}\,0.3$ (right) we probe an even higher Si mass fraction which is present deep inside the O shell (scale of $\mathord{\approx}\,6,000\,\text{km}$), 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 12).

6 Conclusions

In this work, we present the first 3D simulation (full $4\pi$ solid angle) of a fully developed large-scale O-Ne shell merger prior to core-collapse in an $18.88\,\text{M}_{\odot}$ progenitor. The work builds upon the previous study of Müller et al. (2016b) in which they simulated the last minutes of O-shell burning in an $18\,\text{M}_{\odot}$ 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 ($\mathord{\sim}100\,\text{s}$), 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 $\mathord{\approx}\,0.1$. The maximum convective velocities are $\mathord{\sim}\,1300\,\text{km}\,\text{s}^{-1}$ for updrafts and $\mathord{\sim}\,1700\,\text{km}\,\text{s}^{-1}$ 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 $\mathord{\approx}\,50\mathord{-}100$ during the course of evolution. Although the Si-rich material forms a strongly dipolar structure in the merged shell at $(\mathord{\approx}\,330\,\text{s})$, 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 $\ell\mathord{\sim}1\mathord{-}2$ modes) in the flow as well as the distribution of nuclear species.

Another question worth investigating is that, “how the choice of starting time affects the onset and dynamics of the shell merger?” Although, we can definitely say that the shell merger reported in this study is a robust phenomenon, the exact time at which it happens may perhaps change if we start at a different time. This question can only be answered by doing a longer simulation, which is at the moment unfeasible due to limited computing resources. Nevertheless, the shell merger must occur quite late during the evolution of the O shell because this is when the critical conditions (a reduction of the entropy jump between the O shell and Ne shell, and fast O shell convection that allows fast entrainment) are established. This point has been discussed in detail in Collins et al. (2018): For the entropy of the O shell to rise and become similar to that of the Ne shell, O burning and neutrino cooling must have already dropped out of balance due to the contraction of the O shell. The contraction of the O shell in the wake of core contraction allows convection to reach sufficiently high Mach numbers and power strong entrainment.

Simulations by Couch et al. (2015) and Müller et al. (2016b) have already shown that the progenitor structure is genuinely three-dimensional at collapse, but only considered 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 $9\mathord{-}40\,\text{M}_{\odot}$. They find convection with large-scale eddies and the turbulent Mach number $\mathord{\sim}\,0.1$ in their models despite the simulations being short (lasting $\mathord{\approx}\,100\,\text{s}$ before collapse), which nevertheless shows that large-scale perturbations may be a generic feature in pre-supernova stars. 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 et al. (2002, model S20) and in recent studies by Sukhbold & Woosley (2014) and Davis et al. (2018, case C3), and merely be delayed as in model S25 of Rauscher et al. (2002), where the shell merger happens only $5\,\text{s}$ 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 & Ott, 2013; Müller & Janka, 2015; Müller et al., 2017) by aiding the growth of convection and/or instabilities like the standing accretion shock instability (SASI) (Takahashi et al., 2016; Müller et al., 2017). 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 $18.0\,\text{M}_{\odot}$ model of Müller et al. (2017).

Shell mergers will also affect the compactness parameter (O’Connor & Ott, 2011) and thus the explodability with possible consequences for the location and extent of “islands of explodability” (e.g. Sukhbold & Woosley (2014)). Due to the high convective Mach numbers and pronounced global asymmetries, we expect a strong impact on the explosion dynamics for our $18.88\,\text{M}_{\odot}$ 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 et al. (2002) 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 et al. (2018) 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 et al., 2000; DeLaney et al., 2010; Lopez & Fesen, 2018). Such questions will be addressed in the future.

Acknowledgements

This project was supported by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy through Excellence Cluster ORIGINS (EXC-2094)–390783311. 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.

This research has made use of NASA’s Astrophysics Data System.

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, chapter 5)) as $f(\mathbf{r})=\overline{f}+f^{\prime}(\mathbf{r})$, where the angle-averaged value (Reynolds) $\overline{f}$ is defined as

[TABLE]

An extensive physical quantity (for e.g., specific energy $e$) can be written as $g(\mathbf{r})=\tilde{g}+g^{\prime\prime}(\mathbf{r})$, where the angle-averaged value (Favre) is defined as

[TABLE]

The fluctuating components $f^{\prime}$ and $g^{\prime\prime}$ satisfy the following relations

[TABLE]

Appendix B Spherical Harmonic Transforms

The forward spherical harmonic transform of a function $f(\theta,\phi)$ involves computing the coefficients $c_{\ell}^{m}$ according to

[TABLE]

where $Y_{\ell}^{m}(\theta,\phi)$ is the orthonormalized spherical harmonic function of degree $\ell$ and order $m$. The total power $c_{\ell}^{2}$ for a mode with multipole order $l$ is defined as

[TABLE]

We have used the routines from the SHTns library (Schaeffer, 2013) to calculate the spherical harmonic decomposition.