A publishing partnership

The following article is Free article

THE LAST MINUTES OF OXYGEN SHELL BURNING IN A MASSIVE STAR

, , , and

Published 2016 December 12 © 2016. The American Astronomical Society. All rights reserved.
, , Citation Bernhard Müller et al 2016 ApJ 833 124DOI 10.3847/1538-4357/833/1/124

0004-637X/833/1/124

ABSTRACT

We present the first  4π–three-dimensional (3D) simulation of the last minutes of oxygen shell burning in an 18 M supernova progenitor up to the onset of core collapse. A moving inner boundary is used to accurately model the contraction of the silicon and iron core according to a one-dimensional stellar evolution model with a self-consistent treatment of core deleptonization and nuclear quasi-equilibrium. The simulation covers the full solid angle to allow the emergence of large-scale convective modes. Due to core contraction and the concomitant acceleration of nuclear burning, the convective Mach number increases to ∼0.1 at collapse, and an  = 2 mode emerges shortly before the end of the simulation. Aside from a growth of the oxygen shell from 0.51 M to 0.56 M due to entrainment from the carbon shell, the convective flow is reasonably well described by mixing-length theory, and the dominant scales are compatible with estimates from linear stability analysis. We deduce that artificial changes in the physics, such as accelerated core contraction, can have precarious consequences for the state of convection at collapse. We argue that scaling laws for the convective velocities and eddy sizes furnish good estimates for the state of shell convection at collapse and develop a simple analytic theory for the impact of convective seed perturbations on shock revival in the ensuing supernova. We predict a reduction of the critical luminosity for explosion by 12%24% due to seed asphericities for our 3D progenitor model relative to the case without large seed perturbations.

Export citation and abstractBibTeXRIS

1. INTRODUCTION

It is well known that core and shell burning in massive stars typically drives convective overturn (Kippenhahn & Weigert 1990). Although convective heat transport and mixing are inherently multi-dimensional phenomena, the dynamical, convective, Kelvin–Helmholtz, and nuclear timescales are typically too disparate for modeling convection in three dimensions (3D) during most phases of stellar evolution. Spherically symmetric (one-dimensional; 1D) stellar evolution models therefore need to rely on mixing-length theory (MLT; Biermann 1932; Böhm-Vitense 1958) or some generalization thereof (Kuhfuss 1986; Wuchterl & Feuchtinger 1998; Demarque et al. 2008). Such an effective 1D treatment of convection in stellar evolution is bound to remain indispensable even with the advent of modern, implicit hydrodynamics codes (Viallet et al. 2011, 2016; Miczek et al. 2015) that permit multi-dimensional (multi-D) simulations over a wider range of flow regimes and timescales.

The final stages of a massive star before its explosion as a supernova (SN) are among the notable exceptions for an evolutionary phase where the secular evolution timescales are sufficiently short to remain within reach of multi-D simulations (see, e.g., Mocák et al. 2008, 2009; Stancliffe et al. 2011; Herwig et al. 2014 for other examples in the case of low-mass stars). There is also ample motivation for investigating these final stages in 3D. Aside from the implications of multi-D effects in convective shell burning for pulsar kicks (Burrows & Hayes 1996; Goldreich et al. 1997; Lai & Goldreich 2000; Fryer et al. 2004; Murphy et al. 2004) and their possible connection to pre-SN outbursts (Smith & Arnett 2014), they have recently garnered interest as a means for facilitating shock revival in the ensuing supernova (Couch & Ott 2013; Couch et al. 2015; Müller & Janka 2015), which has been the primary impetus for this paper. While the idea that progenitor asphericities arising from convective motions with Mach numbers ∼0.1 can aid shock revival by boosting turbulent motions in the post-shock regions appears plausible, it would be premature to claim that this new idea is a decisive component for the success of the neutrino-driven mechanism. Major questions about this so far undervalued ingredient remain unanswered, and in this paper we shall address some of them.

To evaluate the role of pre-SN seed perturbations in the explosion mechanism, we obviously need multi-D simulations of shell burning up to the onset of core collapse. As shown by the parametric study of Müller & Janka (2015), the typical Mach number and scale of the convective eddies at this stage determine whether the seed asphericities can effectively facilitate shock revival. None of the available multi-D models can reliably provide that information yet. While there is a large body of two-dimensional (2D) and 3D simulations of earlier phases of shell burning (Arnett 1994; Bazan & Arnett 1994, 1998; Asida & Arnett 2000; Kuhlen et al. 2003; Meakin & Arnett 2006, 2007b, 2007a; Arnett & Meakin 2011; Jones et al. 2016), a first, exploratory attempt at extending a model of silicon shell burning up to collapse has only been made recently by Couch et al. (2015), albeit based on a number of problematic approximations. Couch et al. (2015) not only assumed octant symmetry, which precludes the emergence of large-scale modes, but also artificially accelerated the contraction of the iron core due to deleptonization, which leads to a gross overestimation of the convective velocities in the silicon shell, as we shall demonstrate in this paper. Moreover, convective silicon burning often (though not invariably) terminates minutes before collapse in stellar evolution models (see, e.g., Figures 22 and 23 in Chieffi & Limongi 2013 and Figure 16 in Sukhbold & Woosley 2014), as it apparently also does in the 1D model of Couch et al. (2015) calculated with the MESA code (Paxton et al. 2011, 2013). Obviously, simulations covering the full solid angle (4π) with a more physical treatment of the core contraction are required as a next step.

Moreover, the efficiency of progenitor asphericities in triggering shock revival in supernova simulations varies considerably between different numerical studies. The models of Couch & Ott (2013, 2015) are compatible with a small or moderate reduction of the critical luminosity (Burrows & Goshy 1993) for runaway shock expansion. Couch et al. (2015) observed shock revival in their perturbed and non-perturbed model alike, i.e., the perturbations are not crucial for shock revival at all in their study (which appears somewhat at odds with their claims of a significant effect). On the other hand, a much stronger reduction of the critical luminosity of the order of tens of percent has been inferred by Müller & Janka (2015) for dipolar or quadrupolar perturbation patterns based on 2D models with multi-group neutrino transport. These claims may not be in conflict with each other, but could simply result from the different scale and geometry of the pre-collapse velocity/density perturbations, the different progenitor models, and the different treatment of neutrino heating and cooling in these works. A more quantitative theory about the impact of progenitor asphericities on shock revival that could provide a unified interpretation of these disparate findings is still lacking.

In this paper, we attempt to make progress on both fronts. We present the first full-4π–3D simulation of the last minutes of oxygen shell burning in an 18 M star. The model is followed up to collapse by appropriately contracting the outer boundary of the excised (non-convective) core as in the corresponding 1D stellar evolution model computed with the Kepler code (Weaver et al. 1978; Heger & Woosley 2010). By focusing on oxygen shell burning, we avoid the intricacies of deleptonization in the iron core and the silicon shell and the nuclear quasi-equilibrium during silicon burning, so that nucleon burning can be treated with an inexpensive α-network. Our simulation covers the last 293.5 s before collapse to keep the last three minutes (∼9 turnover timescales) free of the artificial transients.

In our analysis of the simulation, we single out the properties of the convective flow that are immediately relevant for understanding pre-collapse asphericities in supernova progenitors and their role in the explosion mechanism, while a more extensive analysis of the flow properties based on a Reynolds decomposition (as in Arnett et al. 2009; Murphy & Meakin 2011; Viallet et al. 2013; Mocák et al. 2014) is left to a future paper. The key question that we set out to answer in this paper is simply: can we characterize the multi-dimensional structure of supernova progenitors (and perhaps their role in the explosion mechanism) already based on 1D stellar evolution models? We argue that this question can be answered in the affirmative, and demonstrate that the typical velocity and scale of the convective eddies comport with the predictions of mixing-length theory (MLT) and linear stability analysis. In preparation for future core-collapse simulations using multi-D progenitors, we develop a tentative theory for the effects of pre-collapse seed perturbations on shock revival that allows one to single out promising models for such simulations. Aside from some remarks on convective boundary mixing, we largely skirt the much more challenging question, whether deviations from MLT predictions have a long-term effect on the evolution of supernova progenitors during earlier phases.

Our paper is structured as follows. In Section 2, we describe the numerical methods used for our 3D simulation of oxygen shell burning and briefly discuss the current version of the Kepler stellar evolution code and the 18 M supernova progenitor model that we consider. In Section 3, we present the results of our 3D simulation, compare them to the 1D stellar evolution model, and show that the key properties of the convective flow are nicely captured by analytic scaling laws. We point out that these scaling laws impose a number of requirements on 3D simulations of shell burning in Section 4. In Section 5, we formulate a simple estimate for the effect of the pre-collapse asphericities with a given typical convective velocity and eddy scale on shock revival. The broader implications of our findings and questions that need to be addressed by 3D stellar evolution models of supernova progenitors are summarized in Section 6. Two appendices address different formulations of the Ledoux criterion (Appendix A) and possible effects of resolution and stochasticity (Appendix B).

2. SETUP AND NUMERICAL METHODS

2.1. The Kepler Stellar Evolution Code

We simulate oxygen shell burning in a non-rotating 18 M solar metallicity star. This stellar model has been evolved to the onset of core collapse with an up-to-date version of the stellar evolution code Kepler (Weaver et al. 1978; Woosley et al. 2002; Heger & Woosley 2010). A 19-species nuclear network (Weaver et al. 1978) is used at low temperatures (up to oxygen burning); at higher temperatures, we switch to a quasi-equilibrium (QSE) approach that provides an efficient and accurate mean to treat silicon burning and the transition to a nuclear statistical equilibrium network after silicon depletion.

The mixing processes taken into account in this model include convective mixing according to MLT, thermohaline mixing according to Heger et al. (2005), and semiconvection according to Langer et al. (1983), but modified for a general equation of state as derived in Heger et al. (2005). All mixing is modeled as a diffusive process with appropriately determined diffusion coefficients. Since we will compare the predictions of MLT and the results of our 3D simulation in some detail, we elaborate further on the numerical implementation of MLT in fully convective (Ledoux-unstable) regions in Kepler, which has been outlined in a more compact form in previous papers (Woosley & Weaver 1988; Woosley et al. 2004). For the implementation of semiconvection and thermohaline convection (which are not immediately relevant for this paper), we refer the reader to Heger et al. (2000, 2005).

MLT assumes that the relative density contrast δρ/ρ between convective updrafts/downdrafts and the spherically averaged background state is related to the deviation of the spherically averaged stratification from convective neutrality and hence to the Brunt–Väisälä frequency ωBV. If the Ledoux criterion for convection is used (as in Kepler), one obtains

for δρ/ρ, where both entropy and composition gradients are implicitly taken into account (see Appendix A). Here, ρ, P, and cs denote the spherically averaged density, pressure, and adiabatic sound speed, and g denotes the local gravitational acceleration. ωBV is the Brunt–Väisälä frequency,7 and Λmix is the mixing length, which is chosen as one pressure scale height hP so that we have

under the assumption of hydrostatic equilibrium. The convective velocity in MLT can then be expressed in terms of ωBV, Λmix, g, and δρ/ρ, and a dimensionless parameter α1 as

Note that different normalizations and default values for α1 are used in the literature. Wherever a direct calibration against observations (as for the solar convection zone, Christensen-Dalsgaard et al. 1996) is not possible, physical arguments can only constrain α1 to within a factor of a few.

Together with the temperature contrast δT between the convective blobs and their surroundings, vconv determines the convective energy flux Fconv,

where cP is the specific heat at constant pressure and α2 is another dimensionless parameter. Note that the second and third lines in Equation (4) implicitly assume that the contribution of composition gradients to the unstable gradient can be neglected inside a convective zone, which is a good approximation for advanced burning stages.

For compositional mixing, Kepler uses a time-dependent diffusion model (Eggleton 1972; Weaver et al. 1978; Heger et al. 2000, 2005) for the evolution of the mass fractions Xi,

where is the diffusive partial mass flux for species i, and the diffusion coefficient is given by

where we have introduced another dimensionless parameter α3. If we introduce the composition contrast δXi = ΛmixXi/∂r between the bubbles and the background state, the symmetry to Equation (4) for the convective energy flux becomes manifest:

We note that only the products α1α2 and α1α3 enter the evolution equations, and we are therefore free to reshuffle an arbitrary factor between α1 and the other two coefficients. In Kepler, we choose α1α2 = 1/2 and α1α3 = 1/6, which is traditionally interpreted as the result of α1 = 1/2, α2 = 1 and α3 = 1/3, where the choice of α3 = 1/3 is motivated by the interpretation of convective mixing as a random-walk process in 3D with mean free path Λmix and an average total velocity (including the non-radial velocity components) vconv. Setting α3 = α2/3 arguably introduces an asymmetry in the equations, but we defer the discussion of its effect to Section 3.4. For extracting convective velocities from the Kepler model, we shall work with the alternative choice of α1 = 1, α2 = 1/2, α3 = 1/6, however, as this gives better agreement with the convective velocity field in our 3D simulation. This is equally justifiable; essentially, this choice amounts to a larger correlation length for velocity perturbations and less perfect correlations between fluctuations in velocity and entropy/composition.

For numerical reasons, ωBV is rescaled before computing the convective energy and partial mass fluxes according to Equations (4) and (7),

where f is an adjustable parameter that is set to f = 0.01 in our model. By rescaling ωBV, convective mixing and energy transport are suppressed until a reasonably large superadiabatic gradient has been established. This procedure avoids convergence problems due to zones switching too frequently between convective stability and instability. The repercussions and limitations of this numerical approach will be discussed in Section 3.1, where we compare the 1D stellar evolution model to our 3D hydrodynamic simulation.

2.2. 1D Supernova Progenitor Model

Entropy, density, and composition profiles of the 1D progenitor model at the onset of collapse are shown in Figure 1. The progenitor has an extended convective oxygen shell of about 0.5 M with a broader convective carbon burning shell directly on top of it. The inner and outer boundaries of the oxygen shell are located at 3000 and 8000 km at the beginning of our 3D simulation and contract considerably until collapse sets in. The entropy jump between the silicon and oxygen shell is relatively pronounced, so that no strong overshooting and/or entrainment at the inner convective boundary is expected because of the strong buoyancy barrier at the interface. The boundary between the oxygen and carbon shell is considerably "softer" with only a small jump of 0.5 kb/nucleon in entropy.

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Top panel: mass fractions Xi of relevant α-elements in the 1D progenitor model at the onset of collapse as a function of enclosed mass m. Bottom panel: profiles of entropy s and density ρ as a function of m. Dashed vertical lines indicate the boundaries of the region simulated in 3D.

Standard image High-resolution image

We note that the balance between energy generation by nuclear burning and neutrino cooling is broken during the final phase before collapse that we are considering here. This is due to the acceleration of shell burning induced by the contraction of the core on a timescale too short for thermal adjustment by neutrino cooling. Different from earlier phases, it is therefore sufficient to follow shell convection in multi-D merely for several overturn timescales to reach the correct quasi-steady state (instead of several Kelvin–Helmholtz timescales for earlier phases to ensure thermal adjustment).

2.3. 3D Simulation

At a time of 293.5 s before the onset of collapse, the stellar evolution model is mapped to the finite-volume hydrodynamics code Prometheus (Fryxell et al. 1989), which is an implementation of the piecewise parabolic method of Colella & Woodward (1984). An axis-free overset "Yin-Yang" grid (Kageyama & Sato 2004; Wongwathanarat et al. 2010), implemented as in Melson et al. (2015), using MPI domain decomposition and an algorithm for conservative advection of scalars (Wongwathanarat et al. 2010; Melson 2013), allows us to retain the advantages of spherical polar coordinates, which are best suited to the problem geometry, while avoiding excessive time-step constraints close to the grid axis. As in Kepler, nuclear burning is treated using a 19-species α-network. The simulations are performed in the so-called implicit large eddy simulations paradigm (Boris et al. 1992; Grinstein et al. 2007), in which diffusive processes (viscosity, mass diffusion, thermal diffusion) are not explicitly included in the equations. Instead, one relies on the truncation errors of the underlying numerical scheme to mimic the effects of irreversible processes taking place at unresolved scales (truncation errors act as an "implicit" sub-grid scale model).

Since there is no convective activity in the Fe core and the Si shell in the last stages before collapse in the Kepler model, we excise the innermost 1.68 M of the core and contract the inner boundary of the computational domain according to the trajectory of this mass shell in the Kepler run from an initial radius of 3000–1974 km at the onset of collapse. At both the inner and outer boundaries, we impose reflecting boundary conditions for the radial velocity, and use constant extrapolation for the non-radial velocity components. The density, pressure, and internal energy are extrapolated into the ghost zones assuming hydrostatic equilibrium and constant entropy. Excising the core not only reduces the computer time requirements considerably, but also allows us to circumvent the complications of deleptonization and Si burning in the QSE regime. The outer boundary is set to a mass coordinate of 4.07 M (corresponding to a radius of 50,000 km) so that the computational domain comprises the outer 0.08 M of the Si shell, the entire O and C shell, and a small part of the incompletely burnt He shell. On the other hand, using an inner boundary condition implies that we cannot address potential effects of shell convection on the core via wave excitation at the convective boundaries, such as the excitation of unstable g-mode (Goldreich et al. 1997; whose growth is likely too slow to be significant; see Murphy et al. 2004) or core spin-up due to angular momentum transport by internal gravity waves (Fuller et al. 2015).

We use a logarithmic radial grid with 400 zones, which implies a radial resolution of Δr/r = 0.7% at the beginning of the simulation. Equidistant spacing in logr is maintained throughout the simulation as the inner boundary is contracted. 56 × 148 angular zones are used on each patch of the Yin-Yang grid, which corresponds to an angular resolution of 2°. A limited resolution study based on two additional models with coarser meshes is presented in Appendix B.

3. SIMULATION RESULTS

In Figures 2 and 3, we show 2D slices depicting the evolution of the mass fraction XSi of silicon and the radial velocity vr to provide a rough impression of the multi-D flow dynamics in our 3D simulation. Convective plumes initially develop on small angular scales in the inner part of the oxygen shell (where the burning rate is high). After about 100 s we see fully developed convection with maximum plume velocities of ∼500 km s−1 that increase toward collapse, and large-scale modes dominate the flow. The latest snapshots at 286 s and 293.5 s suggest the emergence of a bipolar flow structure right before collapse. Large-scale structures are more clearly visible in the velocity field than in XSi. Indeed, the rising plumes enriched in silicon and the sinking plumes containing fresh fuel appear rather "wispy," an impression that is reinforced by the 3D volume rendering of XSi at the onset of collapse in Figure 4.

Figure 2. Refer to the following caption and surrounding text.

Figure 2. Slices showing the mass fraction XSi of silicon (left column) and the radial velocity vr (right column) at times of 20 s, 151 s, and 210 s after the beginning of the 3D simulation (top to bottom). vr is given in units of km s−1. The boundary between the patches of the Yin-Yang grid is located in the right half of the panels at 45° and 135° from the vertical direction. Note that convection initially develops on small angular scales after mapping from the 1D stellar evolution as a strongly superadiabatic gradient builds up in the narrow region of strongest nuclear burning at the base of the oxygen shell (top row). Once convection is fully developed, large-scale overturn emerges. The position of the updrafts and downdrafts shifts freely across the boundaries between the grid patches.

Standard image High-resolution image
Figure 3. Refer to the following caption and surrounding text.

Figure 3. Slices showing the mass fraction XSi of silicon (left column) and the radial velocity vr (right column) at times of 270 s, 286 s, and 293.5 s (onset of collapse) after the beginning of the 3D simulation (top to bottom). vr is given in units of km s−1. Note that wave breaking at the outer boundary of the oxygen shell and the global asymmetry of convective motions become more conspicuous at late times. At the onset of collapse, a bipolar flow pattern emerges (bottom row).

Standard image High-resolution image
Figure 4. Refer to the following caption and surrounding text.

Figure 4. Volume rendering of the mass fraction of silicon at the end of the 3D simulation at 293.5 s (onset of collapse) on one patch of the Yin-Yang grid, showing fuzzy silicon-rich updrafts of hot ashes (red) and silicon-poor downdrafts of fresh fuel. A global asymmetry in the updrafts is clearly visible. The inner boundary of the oxygen shell (cyan) is relatively "hard" due to the strong buoyancy jump between the silicon and oxygen shell and therefore remains almost spherical.

Standard image High-resolution image

Convection also develops in the overlying carbon shell. However, since the convective velocities in the carbon shell are lower, and since this shell extends out to a radius of 27,000 km, convection never reaches a quasi-steady state within the simulation time. We, therefore, do not address convection in the carbon shell in our analysis.

As in earlier studies of mixing at convective boundaries (Meakin & Arnett 2007b), the interface between the carbon and oxygen layer proves unstable to the Holmböe/Kelvin–Helmholtz instability8 with wave breaking leading to the entrainment of material from the carbon shell. The snapshots suggest that such entrainment events become more frequent and violent shortly before collapse.

The convective velocities and eddy scales thus fall roughly into the regime where the parametric study of Müller & Janka (2015) suggests a significant impact of pre-collapse asphericities on shock revival (a convective Mach number of order 0.05 or higher, corresponding to velocities of a few 100 km s−1, and dominant  = 1 or  = 2 modes).

3.1. Flow Dynamics for Quasi-Stationary Convection—Quantitative Analysis and Comparison with MLT

To analyze the flow dynamics more quantitatively, we consider the volume-integrated net nuclear energy generation rate (including neutrino losses) in the oxygen shell, , the volume-integrated turbulent kinetic energy, Er and , contained in the fluctuating components of radial and non-radial velocity components, and profiles of the root-mean-square (rms) averaged turbulent Mach number of the radial velocity fluctuations in Figures 5 and 6. Er, , and  are computed from the velocity field as follows.

where the domain of integration in Equations (9) and (10) extends from the inner boundary radius r to the outer boundary radius r+ of the oxygen shell. Angled brackets denote mass-weighted spherical Favre averages for quantity X,

We note that one does not expect any mean flow in the non-radial directions in the absence of rotation; therefore, only vθ and vφ appear in Equation (10). In Figure 5, we also show the results for and the kinetic energy in convective motions from the 1D Kepler run for comparison. MLT only predicts the radial velocities of rising and sinking convective plumes, so we only compute the 1D analog to Er,

where vconv is calculated according to Equation (3).

Figure 5. Refer to the following caption and surrounding text.

Figure 5. Top: volume-integrated net nuclear energy generation rate in the oxygen shell in the 3D simulation (black) and in Kepler (red). Bottom: kinetic energies (black) and Er (blue) contained in fluctuating non-radial and radial motions in the 3D simulation; see Equations (9) and (10) for definitions. The MLT estimate of the volume-integrated kinetic energy in radial convective motions in the oxygen shell for the Kepler model (red) is computed by using Equation (3) for the convective velocity assuming α1 = 1.

Standard image High-resolution image
Figure 6. Refer to the following caption and surrounding text.

Figure 6. Profiles of the turbulent Mach number of radial velocity fluctuations in the oxygen and carbon shell at different times during the 3D simulation. Note that there is a secular increase in the Mach number in the oxygen shell even after convection has reached a quasi-stationary state due to the contraction of the inner boundary. By contrast, the turbulent Mach number in the carbon shell merely increases because convection has not reached a quasi-stationary state in that shell.

Standard image High-resolution image

The volume-integrated nuclear energy generation rate increases by more than two orders of magnitude during the evolution toward collapse. Due to slight structural adjustments after the initial transient and slightly different mixing in the 3D model, is roughly 30%–50% higher in 3D than in the Kepler for most of the run (see the discussion in Section 3.4), but still parallels the Kepler run quite nicely and perhaps as closely as can be expected given the extreme dependence of the local energy generation on the temperature T during oxygen burning.

The convective kinetic energy oscillates considerably during the first 120 s, but exhibits a smooth secular increase reflecting the acceleration of nuclear burning. Equipartition between the radial and non-radial kinetic energy in convective motions as suggested by Arnett et al. (2009) does not hold exactly, instead we observe for most of the simulation, suggesting that there may not be a universal ratio between the non-radial and radial kinetic energy and that this ratio is instead somewhat dependent on the shell geometry (width-to-radius ratio, ratio of width and pressure scale height), which can vary across different burning shells, progenitors, and evolutionary phases. There may also be stochastic variations in the eddy geometry that the convective flow selects (see Appendix B). Anisotropic numerical dissipation might also account for different results in different numerical simulations. The turbulent Mach number in the oxygen shell (Figure 6) also increases steadily from about 0.04–0.05 after the initial transient to 0.1 at collapse.

Again, there is reasonable agreement between the MLT prediction for the convective kinetic energy and Er in the 3D simulation (Figure 5). and Er are in fact closer to each other than and Er in 3D. Somewhat larger deviations arise immediately prior to collapse when convection is no longer fast enough to adjust to the acceleration of nuclear burning as we shall discuss in Section 3.2.

Except for the last few seconds, the kinetic energy in convection scales nicely with the nuclear energy generation rate both in 1D and 3D. For a case where the convective luminosity Lconv and balance each other in the case of steady-state convection, MLT implies , where Mconv is the mass contained in the convective shell (Biermann 1932; Arnett et al. 2009; note that only the form of the equations is slightly different in these references). In Figure 7, we show the efficiency factors ηconv for the conversion of nuclear energy generation into turbulent kinetic energy9 Eturb,

for both the 3D model (using either the component Er or for Eturb) and the Kepler model (using ), with Λmix set to the pressure scale height at the inner boundary of the oxygen shell. Between 130 and 290 s, ηconv shows only small fluctuations around 0.35 and 0.5 for the kinetic energy in radial and non-radial convective motions in 3D. For the Kepler model, we find similar values around ηconv ≈ 0.37.

Figure 7. Refer to the following caption and surrounding text.

Figure 7. Top panel: efficiency ηconv for the conversion of nuclear energy generation into convective kinetic energy as defined in Equation (14) in the 3D run and the 1D Kepler model. In the 3D case, we compute ηconv both for the kinetic energy in radial motions (Equation (9), black curve) and transverse motions (Equation (10), blue); for the Kepler run (red), we use the energy contained in radial convective motions computed according to Equation (13). By default, we use the pressure scale height in Equation (14) as the mixing or damping length Λdamp. The efficiency is much lower if Λdamp is identified with the radial extent of the convective zone . Bottom panel: comparison of the Brunt–Väisälä frequency ωBV,max at the base of the oxygen shell (black), the reciprocal of the convective turnover time tconv (blue), and the logarithmic derivative of the volume-integrated nuclear energy generation rate (red). The freeze-out of convection (denoted by a dashed vertical line) occurs roughly when .

Standard image High-resolution image

The scaling law can also be understood as resulting from a balance of buoyant driving (or, equivalently, kinetic energy generation by a heat engine) and turbulent dissipation (see, e.g., Arnett et al. 2009 and in a different context Müller & Janka 2015). In this picture, the scaling law emerges if the mixing length is identified with the damping length Λdamp. This identification (Λdamp = Λmix = hP), however, has been criticized on the ground that Λdamp should correspond to the largest eddy scale, which can be considerably larger than hP if low- modes dominate the flow and the updrafts and downdrafts traverse the entire convection zone, which is precisely the situation that is realized in our 3D model. The disparity of the pressure scale height and the eddy scale can be quantified more rigorously by considering the radial correlation length Λcorr for fluctuations in the radial velocity, . Following Meakin & Arnett (2007b) and Viallet et al. (2013), we compute the vertical correlation length Λcorr as the full width at half maximum of the correlation function ,

The correlation function is computed at a radius of r = 4000 km in the inner half of the oxygen shell. Λcorr is shown in Figure 8 and compared to the pressure scale height Λmix = hP at the inner boundary of the oxygen shell and the extent Λconv of the convective region. Once convection is fully developed, we clearly have Λcorr > Λmix and Λcorr ≈ Λconv/2 (as expected for updrafts and downdrafts reaching over the entire zone).

Figure 8. Refer to the following caption and surrounding text.

Figure 8. Top: outer and inner boundary radius r+ and r of the oxygen shell as functions of time. Bottom: correlation length Λcorr for the radial velocity computed at r = 4000 km (black), pressure scale height hP at the base of the oxygen shell (red), and radial extent of the oxygen shell (blue) as a function of time.

Standard image High-resolution image

Arnett et al. (2009) argued that the damping length should be on the order of the width Λconv of the convective zone under such circumstance. However, if we compute the efficiency factor ηconv based on Λconv,

we obtain suspiciously low values ηconv ≲ 0.1. This suggests that the effective damping length is set by the pressure scale height (or a multiple thereof) after all. One could opine that the energetics of the flow might still be described adequately by Λdamp = Λcorr, and that the efficiency factor ηconv merely happens to be relatively low.

We argue, however, that there is a deeper reason for identifying Λdamp with a multiple of the pressure scale height in the final phases of shell convection when neutrino cooling can no longer balance nuclear energy generation. The crucial point is that the average distance after which buoyant convective blobs have to return their excess enthalpy h' to their surroundings cannot become arbitrarily large in a steady-state situation, and since enthalpy and velocity fluctuations are correlated (), this also limits the damping length.

During the final stages, nuclear energy generation, convective transport, and turbulent dissipation must balance each other in such a way as to avoid both a secular build-up of an ever-growing unstable entropy/composition gradient (although the spherically average stratification always remains slightly unstable) and a complete erasure of the superadiabatic gradient. Assuming that the Brunt–Väisälä frequency is primarily set by the gradient of the entropy s, this implies ∂2 s/∂rt ≈ 0, and hence roughly constant entropy generation,

throughout the convective region. In the late pre-collapse stages, we can relate to the local nuclear energy generation rate and the derivative of the "total" convective luminosity Lconv,

Here, Lconv denotes the net total energy flux resulting from fluctuations (denoted by primes) in the total energy density and velocity around their spherical Favre average,

where e is the specific internal energy, and P the pressure. Note that in formulating Equation (18), we implicitly assumed that is equal to the rate of change of the Favre average of the internal energy e (instead of the total energy density, which includes the contribution of the turbulent kinetic energy). This assumption is justified for steady-state convection in the late pre-collapse phase because of moderate Mach numbers and the minor role of neutrino cooling. These two factors imply that the energy that is generated by nuclear reactions and distributed throughout the unstable region by convection mostly goes into internal energy (whereas our argument cannot be applied to earlier phases where neutrino cooling and nuclear energy generation balance each other).

Figure 9 shows the two terms contributing to in Equation (18) based on Favre averages over a few time steps around 210 s, and demonstrates that is indeed roughly constant throughout the convective region. Since strong nuclear burning is confined to a narrow layer at the bottom of the convective shell, we even have throughout a large part of the shell, and for a stratification with roughly ρ ∝ r−3 and T ∝ r−1, this leads to

For such an idealized case, one can directly compute that the energy transported by convective blobs from the lower boundary must be dissipated after an average distance of

where ξ = r+/r is the ratio of the outer and inner boundary radius. Evidently, Λdamp grows only moderately at large ξ, and is always smaller than .

Figure 9. Refer to the following caption and surrounding text.

Figure 9. Top: total convective luminosity (including the kinetic energy flux) at 210 ms in the 3D simulation as a function of enclosed mass m. Bottom: quantities determining the spherically averaged entropy production. The term T−1dLconv/dm stemming from the divergence of the total convective luminosity is shown in red, the entropy production due to the nuclear source term (neglecting terms in the chemical potential of the different nuclear species) is shown in blue, and the black curve denotes the sum of both terms. The curves are computed from averages over several time steps.

Standard image High-resolution image

It thus appears unlikely that large damping lengths can be realized in very extended convection zones in the final pre-collapse stage. This is decidedly different from earlier stages with strong neutrino cooling in the outer part of the convective zone, for which Arnett et al. (2009) found high values of . As outlined before, the different behavior is likely due to the specific physical conditions right before collapse; in the absence of strong cooling, the self-regulatory mechanism that we outlined above automatically ensures that Λdamp cannot be considerably larger than the pressure scale height. Thus, the implicit identification of Λdamp and hP (or a multiple thereof) in MLT is likely less critical for shell convection right before collapse than for earlier phases.

However, it still remains to be determined whether the damping length can reach considerably higher values in deep convection zones with ξ ≫ 1 during earlier stages when nuclear energy generation and neutrino cooling balance each other. Since neutrino cooling generally decreases with radius within a shell, it can still be argued that the convective luminosity must decay not too far away from the burning region. Thus, an analog to Equation (21) could still hold, and the damping length would only increase slowly with the width of the shell in the limit of large ξ. In that case, the difference between our simulation and the results of Arnett et al. (2009) would merely be due to a different depth of the convective zone, which is much deeper in our model (4–5 pressure scale heights as opposed to ∼2 pressure scale heights in Arnett et al. 2009) so that we approach a "saturation limit" for the damping length and can more conveniently distinguish the damping length from the width of the convective zone, since the different length scale are sufficiently dissimilar.

Radial profiles of the convective velocities also point to reasonable agreement between MLT and the 3D simulation. In the upper panel of Figure 10, we compare the convective velocities from Kepler to rms averages of the fluctuations of the radial velocity (δvr) and the transverse velocity component (δvt) at 210 s,

We also compare these to the MLT estimate vconv = ωBV Λmix computed from the Brunt–Väisälä frequency for the spherically averaged stratification of the 3D model. It is evident that the agreement especially between δvr and the convective velocity in Kepler is very good in the oxygen shell. In large parts of the shell, vconv = ωBVΛmix is also in very good agreement with δvr, which again demonstrates that the choice of the pressure scale height as the acceleration and damping length for convective blobs is a reasonable choice. However, no reasonable comparison can be made in the outer part of the oxygen shell, where ωBV is formally negative. This is due to the strong aspherical deformation of the shell boundary and the entrainment of light, buoyant material from the carbon shell; the fact that the outer part of the oxygen shell is formally stable if ωBV is computed from spherical averages of the density and pressure is thus merely a boundary effect and has no bearing on the validity of MLT in the interior of the shell.

Figure 10. Refer to the following caption and surrounding text.

Figure 10. Top panel: comparison of radial and transverse rms velocity fluctuations δvr (black) and δvt (blue) in the 3D model at 210 s to the convective velocity vconv computed in the Kepler model using MLT (red), and to ωBV hP (violet). Bottom panel: comparison of the Brunt–Väisälä frequency computed from spherical averages of the pressure, density, and sound speed in the 3D run (black) and in the 1D Kepler model. Note that there is a formally stable region around the boundary between the oxygen and carbon shell in the 3D model due to the aspherical deformation of the shell boundary and the entrainment of material from the carbon shell, which increases the spherically averaged entropy in the outer parts of the oxygen shell.

Standard image High-resolution image

The good agreement between the 3D simulation and the Kepler model may seem all the more astonishing considering the rescaling of the Brunt–Väisälä frequency according to Equation (8) for stability reasons. However, this procedure is justified by the fact that the convective luminosity automatically adjusts itself in such a way as to avoid a secular build-up of ωBV as discussed before. In a steady state, the convective luminosity in MLT in a shell will roughly balance the nuclear energy generation rate, , regardless of whether ωBV is rescaled or not. If a rescaling factor is introduced in Equation (4), the result is simply that a larger ωBV is maintained under steady-state conditions to balance the rescaling factor. Except for pathological situations, the convective energy flux and the convective velocities are thus essentially unaffected by this procedure. The superadiabaticity of the stratification is changed, however. For convection at low Mach numbers, it will be systematically overestimated. This trend is evident from the lower panel of Figure 10, which compares ωBV in Kepler and the 3D simulation. Since convection is not extremely subsonic in our case, the rescaling factor is only slightly smaller than unity, and the superadiabaticity in the 1D and 3D model remains quite similar.

3.2. Freeze-out of Convection

MLT in Kepler thus provides good estimates for the typical convective velocities in the final stages of oxygen shell burning as long as a steady-state balance between nuclear energy generation, convective energy transport, and turbulent dissipation is maintained. However, steady-state conditions are not maintained up to collapse. Figure 7 shows that the growth of the turbulent kinetic energy can no longer keep pace with the acceleration of nuclear burning in the last few seconds before collapse, where ηconv drops dramatically.

The time at which convection "freezes out" can be nicely determined by appealing to a timescale argument: freeze-out is expected once the nuclear energy generation rate (which sets the Brunt–Väisälä frequency and the convective velocity under steady-state conditions) changes significantly over a turnover timescale. More quantitatively, the efficiency factor ηconv drops abruptly once the freeze-out condition

is met, as shown in the bottom panel of Figure 7. Equivalently, the freeze-out condition can be expressed in terms of the convective turnover time tconv,

where is an appropriate global average of the convective velocity, e.g.,

Using these definitions, we find that freeze-out occurs roughly when

which may be even more intuitive than Equation (24).

Somewhat astonishingly, the Kepler run shows a similar drop of ηconv in the last seconds, although MLT implicitly assumes steady-state conditions when estimating the density contrast and the convective velocity. Kepler still overestimates the volume-integrated turbulent kinetic energy somewhat after freeze-out (Figure 5), but the discrepancy between the 1D and 3D models is not inordinate.

The key to the relatively moderate differences can be found in profiles of the turbulent convective Mach number vconv/cs in Kepler and in 3D at the onset of collapse in Figure 11. Evidently, MLT only overestimates the convective velocities in a narrow layer at the lower boundary of the oxygen shell, where the acceleration of nuclear burning greatly amplifies the superadiabaticity of the stratification (as quantified by ωBV). This immediately increases vconv, whereas the convective velocity field adjusts only on a longer timescale () in 3D. However, even in the Kepler run, the convective velocities in the middle and outer region of the oxygen shell remain unaffected by the increase of the ωBV close to the inner shell boundary. Different from the innermost region, where ωBV reacts instantaneously to the nuclear source term, ωBV (and hence the convective velocity in the outer region) responds to the accelerated burning on a diffusion timescale, which is again of order . For a slightly different reason (insufficient time for convective diffusion versus insufficient time for the growth of plumes), the Kepler run therefore exhibits a similar freeze-out of convection as the 3D model. We thus conclude that the volume-integrated turbulent kinetic energy and the average convective Mach number in 1D stellar evolution codes still provide a reasonable estimate for the state of convection even right at collapse. The spatial distribution of the turbulent kinetic energy, on the other hand, appears more problematic; it will be somewhat overestimated in the shell source at collapse due to the instantaneous reaction of ωBV to the increasing burning rate.

Figure 11. Refer to the following caption and surrounding text.

Figure 11. Profiles of the rms-averaged turbulent Mach number of radial velocity fluctuations at the onset of collapse in the 3D model (black) and the 1D Kepler model (red). Dashed lines denote the boundaries of the domain simulated in 3D in Prometheus. The turbulent Mach number in 1D and 3D agrees well in the bulk of the oxygen shell, but the acceleration of nuclear burning and the concomitant increase of the Brunt–Väisälä frequency artificially increases the convective velocities close to the base of the shell in Kepler, as MLT immediately translates the increase in ωBV into an increase in convective velocity. Note that high nuclear burning rates in individual zones inside the silicon core produce formally unstable zones in the Kepler run shortly before collapse, which do not affect the evolution of the model to any significant degree.

Standard image High-resolution image

The rescaling of ωBV in Kepler according to Equation (8) can also affect the time of freeze-out at a minor level. For a convective Mach number of ∼0.1, the rescaling procedure changes ωBV only by ∼30%, and given the very rapid increase of , this will not shift the time of freeze-out appreciably.

3.3. Scale of Convective Eddies

The role of progenitor asphericities in the explosion mechanism depends not only on the magnitude of the convective velocities in the burning shells, but also on the angular scale of the infalling eddies. MLT does not make any strong assumptions about the eddy scale; it assumes a radial correlation length for entropy and velocity perturbations, but such a correlation length can, in principle, be realized with very different flow geometries. Empirically, simulations of buoyancy-driven convection in well-mixed shells are usually characterized by eddies of similar radial and angular extent d that reach across the entire unstable zone (e.g., Arnett et al. 2009). The dominant modes are also typically close in scale to the most unstable modes in the linear regime (Chandrasekhar 1961; Foglizzo et al. 2006), which have . This correspondence between the linear and nonlinear regime has sometimes been justified by heuristic principles for the selection of the eddy scale based on maximum kinetic energy or maximum entropy production (Malkus & Veronis 1958; Martyushev & Seleznev 2006). Expressing the balance of kinetic energy generation due to the growth of an instability with a scale-dependent growth rate ω (d) and turbulent dissipation for the dominant mode in a shell with mass M yields

for the change of the kinetic energy Ekin in a given mode. The dominant mode(s) in the nonlinear regime will be the one(s) for which

is maximal, which actually suggests a bias toward slightly larger scales than in the linear regime.

A superficial inspection of Figures 2 and 3 already reveals that our 3D models conform to the typical picture with . More quantitatively, the dominance of large-scale modes is shown by a decomposition of the radial velocity in the inner half of the oxygen shell (at a radius of 4000 km) into spherical harmonics (for more sophisticated decompositions of the flow field, see Chatzopoulos et al. 2014; Fernández et al. 2014). In Figure 12, we plot the total power for each multipole order ,

which shows a clear peak at low that slowly moves from  = 4 down to  = 2 over the course of the simulation. The tail at high above the typical eddy scale roughly exhibits an −5/3 slope as expected for a Kolmogorov-like turbulent cascade (Kolmogorov 1941) because of the rough proportionality between and the wave number (Peebles 1993).

Figure 12. Refer to the following caption and surrounding text.

Figure 12. Power in different multipoles for the decomposition of the radial velocity at r = 4000 km into spherical harmonics in the 3D model at different times computed according to Equation (30). The dominant angular wave number shifts from  = 3–5 to  = 2 over the course of the simulation. The dashed line indicates a slope of −5/3, which is roughly expected for a Kolmogorov spectrum above the injection scale in wave number (i.e., at smaller spatial scales).

Standard image High-resolution image

The dominant eddy scale is consistent with the crude estimate that the dominant is given by the number of convective eddies of diameter that can be fitted into one hemisphere of the convective shell (Foglizzo et al. 2006),

This estimate for the dominant multipole order is plotted in Figure 13. It agrees well with spectra of the radial velocity, although it may not clearly predict the emergence of the dominant quadrupole at the end (which is compatible with our argument that the dominant angular scale for fully developed convection is slightly larger than in the linear regime). The slowly changing geometry of the shell evidently accounts nicely for the secular trend toward modes of lower . Figure 8 reveals that both the contraction of the inner boundary of the shell by about one-third in radius and a secular expansion of the (somewhat ill-defined) outer shell boundary contribute to this trend. The fast change of right before collapse is clearly due to the contraction, however, as the outer boundary radius decreases again shortly before collapse.

Figure 13. Refer to the following caption and surrounding text.

Figure 13. Estimate of the typical angular wave number of convection in the linear regime from the inner and outer boundary radius of the oxygen shell according to Equation (31) and Foglizzo et al. (2006). The rapid drop at the end of the simulation is evident and suggests that the emergence of a global  = 2 mode is due to the rapid contraction of the iron and silicon core shortly before collapse.

Standard image High-resolution image

The expansion of the outer boundary is not seen in the Kepler model and is the result of entrainment of material from the carbon shell (see Section 3.4 below). If the amount of entrainment is physical, this is another reason to suspect that estimates of the dominant angular scale based on stellar evolution models using Equation (31) will slightly overestimate the dominant . Considering uncertainties and progenitor variations in the shell structure, Equation (31) nonetheless furnishes a reasonable zeroth-order estimate of the typical eddy scale.

3.4. Comparison of Convective Mixing in 1D and 3D

Although the properties of the velocity field are more directly relevant for the potential effect of progenitor asphericities on supernova shock revival, some remarks about convective mixing in our 3D model are still in order.

In Figure 14, we compare spherically averaged profiles of the entropy s, and the mass fractions of oxygen, silicon, and sulphur (XO, XSi, and XS) from the 3D model to the Kepler run at a time of 210 ms. Although the treatment of convective mixing as a diffusive process in 1D has sometimes been criticized (Arnett et al. 2009), the differences in the interior of the oxygen shell remain minute; the most conspicuous among them are the somewhat steeper gradients in the mass fractions in Kepler. These could potentially contribute (on a very modest level) to the lower total nuclear energy generation rate in Kepler, since the nuclear energy generation rate is roughly proportional to the square of the mass fraction XO of oxygen in the burning region. Even if we account for spatial fluctuations in the composition by computing , the compositional differences do not appear to be sufficiently large to explain the different burning rates; temperature changes due to hydrostatic adjustment thus seem to be the major cause of the somewhat higher total nuclear energy generation rate in Prometheus.

Figure 14. Refer to the following caption and surrounding text.

Figure 14. Spherically averaged profiles of the entropy (violet curves, top panel) and the mass fractions of oxygen (black), silicon (red), and sulfur (blue) in the 3D run (solid curves) at 210 s compared to profiles from the 1D Kepler model (dashed) at the same time. Note that the slope in the mass fractions is somewhat steeper in the Kepler model, which we ascribe to the use of an extra factor of 1/3 in the diffusion equation for compositional mixing.

Standard image High-resolution image

It is unclear whether the composition gradients are really an artifact of MLT; we find it equally plausible that they simply stem from the choices of different coefficients α2 and α3 for energy transport and compositional mixing in Equations (4) and (7). The introduction of an additional factor of 1/3 in Equation (7) is typically justified by interpreting turbulent mixing as a random-walk process of convective blobs with a mean free path Λmix and a total velocity vconv with random orientation, which translates into a radial correlation length and an rms-averaged radial velocity of . However, the mixing length and MLT velocity are implicitly identified with the radial correlation length and in Equation (4) already, so that the choice α3 = α2 rather than α3 = α2/3 is arguably more appropriate. With such a (more parsimonious) choice of parameters, the composition gradients would be flattened considerably.

Figure 14 also shows evidence of boundary mixing (entrainment; Fernando 1991; Strang & Fernando 2001; Meakin & Arnett 2007b) that is not captured in the Kepler run. The fact that the entropy and composition gradients are smeared out at the boundaries (especially at the outer boundary) is mostly due to the aspherical deformation of the shell interface by Kelvin–Helmholtz/Holmböe waves; the shell boundary remains relatively well defined in the multi-D snapshots in Figures 2 and 3. However, the oxygen shell is clearly expanding in m at the outer boundary. To capture the increase of the total mass Mconv in the convective oxygen shell, we integrate the mass in all zones with entropies between 3.6kb/nucleon and 5.2kb/nucleon (Figure 15). Mconv increases by about 0.05 M over the course of the simulation with some evidence for higher toward the end, corresponding to an entrainment rate of 1.4 × 10−4 M, which is also roughly the maximum value of the turbulent mass flux that is reached in the formally stable region around the outer boundary (Figure 16).

Figure 15. Refer to the following caption and surrounding text.

Figure 15. Mass Mconv contained in the convective oxygen shell in the 3D simulation as a function of time. The mass increases by about 0.05 M due to entrainment, which accelerates slightly toward the end of the simulation as a result of higher convective velocities and Mach numbers. Note that the small changes in the first ∼30 s are simply due to the advection of the entropy discontinuities over the grid in the wake of hydrostatic adjustment, as a result of which cells can jump around the threshold entropies of 3.6kb/nucleon and that we use to define the shell boundaries. "Physical" entrainment begins once the first convective plumes reach the boundary between the carbon and oxygen shell around ∼30 s (denoted by a vertical line).

Standard image High-resolution image
Figure 16. Refer to the following caption and surrounding text.

Figure 16. Turbulent mass flux in the 3D model at a time of 210 s as a function of enclosed mass m. Positive values around the outer boundary of the oxygen shell at m ≈ 2.3 M indicate entrainment of material from the carbon shell. The peak value of 1.4 × 10−4 M s−1 roughly corresponds to the average entrainment rate over the course of the simulation.

Standard image High-resolution image

Higher resolution is ultimately required to decide whether this entrainment rate is physical or partially due to numerical diffusion, which could lead to an overestimation of the amount of entrained mass in wave breaking events (see Appendix B). Our simulations are, however, consistent with semi-empirical entrainment laws found in the literature. Laboratory experiments and simulations (Fernando 1991; Strang & Fernando 2001; Meakin & Arnett 2007b) suggest

for the entrainment rate in the relevant regime of the bulk Richardson number RiB and a dimensionless proportionality constant A. RiB is defined in terms of the density contrast δρ/ρ at the interface, the gravitational acceleration g, the typical convective velocity vconv, and the eddy scale Λ as

If we identify Λ with the pressure scale height, this amounts to

In our case, we have δρ/ρ = 0.1, and with vconv = 2. 5 × 107 cm s−1 (corresponding to the non-radial velocities near the boundary, which are relevant for the dynamics of interfacial Holmböe/Kelvin–Helmholtz waves), we obtain RiB = 17, indicating a very soft boundary. Together with an average convective velocity of ∼200 km s−1 and an average entrainment rate of 1.4 × 10−4 M, this points to a low A ∼ 0.1 in the entrainment law (32), although the ambiguities inherent in the definition of RiB can easily shift this by an order of magnitude, which may account for the higher value A ≈ 1 obtained by Meakin & Arnett (2007b). It is obvious that the calibration of the entrainment law is fraught with ambiguities: if we calibrate Equation (32) by using a global average for vconv,

and the initial values for δρ/ρ, and the density ρ at the outer boundary radius r in (32) and (34), the time-dependent entrainment rate is well fitted by A = 0.37 (Figure 17). If anything, relatively low values of A merely demonstrate that entrainment in our 3D model is no more affected by numerical diffusion than in comparable simulations. Considering the low value of the bulk Richardson number and the small entropy jump of ∼0.5kb/nucleon, which should be conducive to entrainment effects, the dynamical impact of boundary mixing in our simulation is remarkably small, but its long-term effect warrants further investigation.

Figure 17. Refer to the following caption and surrounding text.

Figure 17. Comparison of the measured entrainment rate in the 3D simulation (black) and a fit based on Equation (32) (red) computed using A = 0.37 and a global average for the convective velocity (see the text for details). Overall, the time-dependent entrainment rate nicely follows Equation (32). Note that no data is shown later than 290 s, as the detection of the outer boundary of the oxygen shell becomes problematic due to increasingly violent boundary mixing shortly before collapse. As in Figure 15, the dashed vertical line indicates the time when convective plumes first encounter the outer boundary and physical entrainment begins.

Standard image High-resolution image
Figure 18. Refer to the following caption and surrounding text.

Figure 18. Efficiency factors ηconv for the conversion of the total nuclear energy generation rate into turbulent energy in radial motions (Er, solid lines) and non-radial motions (, dashed lines) for the baseline run with an angular resolution of 3° degrees and Nr = 400 radial zones (black curves) and for two low-resolution runs with an angular resolution of 2° and a radial resolution of 400 zones (red) and 286 zones (blue).

Standard image High-resolution image

4. REQUIREMENTS FOR 3D PRE-SUPERNOVA SIMULATIONS

If 3D simulations of shell burning in massive stars are to be used as input for core-collapse simulations, it is essential that the typical convective velocities and eddy scales are captured accurately. The analysis of our model in the preceding section provides guidelines about the approximations that can (or cannot) be justified in such simulations.

The emergence of large-scale motions ( = 2 modes) during the final phase of our model implies that pre-SN models generally need to cover the full solid angle (which has been done previously for oxygen shell burning only by Kuhlen et al. 2003, albeit for an earlier phase). However, for sufficiently narrow convective shells, simulations restricted to a wedge or octant may still cover the flow geometry accurately notwithstanding that such symmetry assumptions remain questionable in the ensuing SN phase. Thus, for the pre-SN phase, the assumption of octant symmetry in Couch et al. (2015) may be adequate for their model of silicon shell burning, which has toward the end of the simulation. The eddies should then remain of a moderate scale with a preferred of .

An accurate treatment of nuclear burning is even more critical because of the scaling of convective velocities with . Since the nuclear generation rates in the silicon and oxygen shell are sensitive to the contraction of the deleptonizing iron core, this not only applies to the burning shell in question itself, but also to the treatment applied for the iron core. If the contraction of the core is artificially accelerated as in Couch et al. (2015), this considerably reduces the nuclear timescale in the outer shells as well. For example ∼0.2 M of intermediate mass-elements in the silicon shell are burned to iron group elements within 160 s in the 3D model of Couch et al. (2015), i.e., silicon burning, on average, proceeds 6.25 times faster than in the corresponding stellar evolution model, where this takes 1000 s. This suggests an artificial increase of the convective velocities by 84% in their 3D model.

Approximations that affect the nuclear burning timescale are also problematic because they change the ratio , which plays a crucial role in the freeze-out of convective motions shortly before the onset of collapse (see Section 3.2). If the nuclear burning is artificially accelerated and continues until collapse, then the freeze-out will occur somewhat earlier, which may compensate the overestimation of convective velocities discussed before. However, the simulation of Couch et al. (2015) suggests that the opposite may also occur: in their 3D model, silicon burning slows down toward the end of their simulation as the shell almost runs out of fuel. In the corresponding 1D stellar evolution model, convection in the original silicon shell has already died down completely, as can be seen from their Figure 2, which shows non-zero convective velocities only in regions with Ye = 0.5. While it is conceivable that convection subsides more gradually in 3D as the available fuel is nearly consumed—probably over a few turnover timescales—increasing the ratio τconv/τnuc by more than a factor of ≳3 evidently introduces the risk of artificially prolonging convective activity in almost fully burned shells.

Other worries about the feasibility of multi-D simulations of supernova progenitors include the problem of thermal adjustment after mapping from a 1D stellar evolution model as well as artificial boundary mixing. We have largely circumvented the problem of thermal adjustment in this study by focusing on the final stages. The somewhat higher nuclear burning rate in the 3D model (by up to ∼50% compared to Kepler), which may be due to physical multi-D effects or transients after the mapping, such as an adjustment to a new hydrostatic equilibrium, suggests that even for a setup where the problem of hydrostatic and thermal adjustment is rather benign, we still face uncertainties of the order of 15%—because of —in the final convective velocity field at collapse. The slight expansion of the outer boundary of the oxygen shell, which may be the result of an adjustment effect or driven by (physical) entrainment, also deserves attention, because it plays some role in fostering the emergence of an  = 2 mode right before collapse. It appears less worrisome, however, since there are natural variations in shell geometry anyway, and since the emergence of the  = 2 mode may still be primarily driven by the contraction of inner shell boundary. There is no evidence for artificial boundary mixing at this stage, although further high-resolution tests remain desirable.

5. EFFECT OF CONVECTIVE SEED PERTURBATIONS ON SUPERNOVA SHOCK REVIVAL

With typical convective Mach numbers of ∼0.1 and a dominant  = 2 mode at collapse, the progenitor asphericities fall in the regime where they may be able to affect shock revival in the ensuing core-collapse supernova, as has been established by the parameter study of Müller & Janka (2015). Considering that several recent works have shown that the conditions for shock revival in multi-D can be captured with good accuracy by surprisingly simple scaling laws (Müller & Janka 2015; Janka et al. 2016; Müller et al. 2016; Summa et al. 2016) that generalize the concept of the critical luminosity (Burrows & Goshy 1993) to multi-D, it is reasonable to ask whether the effect of progenitor asphericities can also be predicted more quantitatively by simple analytic arguments. Given the good agreement between our 3D model and MLT, such a theory could help to better identify progenitors for which convective seed asphericities play a major role in the explosion before investing considerable computer time into multi-D simulations.

The key ingredient to accomplish this consists in a first quantitative theory for the interaction of asymmetries in the supersonic infall region with the shock, which Müller & Janka (2015) only described qualitatively as "forced shock deformation." The starting point is the translation of initial radial velocity perturbations into density perturbations at the shock due to differential infall (Müller & Janka 2015),

which can also be understood more rigorously using linear perturbation theory (Goldreich et al. 1997; Takahashi & Yamada 2014). Note that we now designate the typical convective Mach number during convective shell burning simply as Ma to avoid cluttered notation. The perturbations in the transverse velocity components are amplified as r−1 (Goldreich et al. 1997) and are roughly given by

where and rini are the initial sound speed and radius of the shell before collapse and rsh is the shock radius. Radial velocity perturbations only grow with r−1/2 (Goldreich et al. 1997) and can therefore be neglected.

5.1. Generation of Turbulent Kinetic Energy by Infalling Perturbations

The interaction of the pre-shock perturbations with the shock can then be interpreted as an injection of additional turbulent kinetic energy into the post-shock region. While this problem has not yet been addressed in the context of spherical accretion onto a neutron star, the interaction of planar shocks with incident velocity and density perturbations has received some attention in fluid dynamics (Ribner 1987; Andreopoulos et al. 2000; Wouchuk et al. 2009; Huete Ruiz de Lira 2010; Huete et al. 2012). The perturbative techniques that allow a relatively rigorous treatment in the planar case cannot be replicated here, and we confine ourselves to simple rule-of-thumb estimate for the generation of turbulent energy by the infalling perturbations: if we neglect the deformation of the shock initially, we can assume transverse velocity perturbations δvt and density fluctuations δρ/ρ compared to the spherically averaged flow are conserved across the shock as a first-order approximation. The anisotropy of the ram pressure will also induce pressure fluctuations δP/P ∼ δρ/ρ downstream of the shock. In a more self-consistent solution, these pressure fluctuations would induce lateral flows and modify the shape of the shock, and larger vorticity perturbations would arise if the shock is asymmetric to begin with (which is important in the context of the SASI (Foglizzo et al. 2007; Guilet & Foglizzo 2012). As a crude first-order estimate, such a rough estimate is sufficient for our purpose; it is not incompatible with recent results about shocks traveling in inhomogeneous media (Huete Ruiz de Lira 2010; Huete et al. 2012).

From the density and pressure perturbations δρ/ρ ∼ δP/P ∼ Ma and transverse velocity perturbations downstream of the shock, we can estimate fluxes of transverse kinetic energy (Ft), acoustic energy (Fac), and an injection rate of kinetic energy due to the work done by buoyancy during the advection of the accreted material through down to the gain radius rg. Ft is roughly given by

where we approximated the initial sound speed as , which is a good approximation for the shells outside the iron core. Note that we use a typical ratio rsh/rg = 3/2 during the pre-explosion phase to express Ft in terms of the gravitational potential of the gain radius; the reason for this will become apparent when we compare the injection rate of turbulent kinetic energy at the shock to the contribution from neutrino heating.

Following Landau & Lifshitz (1959), the acoustic energy flux can be estimated by assuming that the velocity fluctuations in acoustic waves are roughly δv ∼ cs δP/P (where cs is the sound speed behind the shock and δP ≈ Ma P). The post-shock pressure P can be determined from the jump conditions,

where ρpre and ρ are the pre- and post-shock density, β ≈ 7 is the compression ratio in the shock, and vpre is the pre-shock velocity, which we approximate as . The acoustic energy flux is thus

Here, is the spherical average of the post-shock velocity, and has been used following Müller & Janka (2015).

Finally, the gravitational potential energy corresponding to density fluctuations δρ will be converted into kinetic energy by buoyancy forces at a rate of10

Especially for moderate Mach numbers, Fpot is clearly the dominating term, as the flux of acoustic and transverse kinetic energy scale with Ma2.

In the absence of infalling perturbations, Müller & Janka (2015) established a semi-empirical scaling law that relates transverse kinetic energy stored in the post-shock region to the volume-integrated neutrino heating rate , the mass in gain region, Mg, and the shock and gain radius

At least for convection-dominated models, this scaling law can be understood as the result of a balance between kinetic energy generation by buoyancy and turbulent dissipation with a dissipation length Λ = rsh − rg (see also Murphy et al. 2013). Assuming a local dissipation rate of , this leads to

from which Equation (42) immediately follows.

In the presence of infalling perturbations, it is natural to add another source term to Equation (43),

To keep the calculation tractable, we only include the dominant contribution Fpot arising from infalling perturbations and discard Fac and Ft.

However, this obviously poses the question about the appropriate choice for Λ, which can now no longer assumed to be simply given by rsh − rg. To get some guidance, we can consider the limit in which neutrino heating is negligible; here, the appropriate choice for Λ is clearly given by the scale of the infalling perturbations, i.e., Λ ≈ πrsh/ in terms of their typical angular wave number . Hence we find

in this limit. The general case can be accommodated by simply interpolating between the two limits,

We emphasize that a different dissipation length enters in both terms: in the limit of neutrino-driven convection with small seed perturbations, the dissipation length is given by the width of the gain layer, whereas the dissipation length πrsh/ can be considerably larger for "forced" convection/shock deformation due to infalling perturbations with small .

For deriving the modification of the critical luminosity, it will be convenient to express in terms of its value in the limit of small seed perturbations (Equation (42)) and a correction term ψ,

where ψ is defined as

Different from the case of negligible seed perturbations, it is hard to validate Equation (46) in simulations. In the 2D study of Müller & Janka (2015), the amplitudes of the infalling perturbations change significantly over relatively short timescales, and the phase during which they have a significant impact on the turbulent kinetic energy in the post-shock region, but have not yet triggered shock revival and was therefore too short to detect any deviations from Equation (42), especially since the turbulent kinetic energy fluctuates considerably around its saturation value in 2D.

5.2. Effect on the Heating Conditions and the Critical Luminosity

Conceptually, the steps from Equation (46) to a modified critical luminosity are no different from the original idea of Müller & Janka (2015), i.e., one can assume that the average shock radius can be obtained by rescaling the shock radius for the stationary 1D accretion problem with a correction factor that depends on the average rms Mach number in the gain region (Equation (42) in Müller & Janka 2015),

which then leads to a similar correction factor for the critical values for the neutrino luminosity and mean energy Lν and Eν (Equation (41) in Müller & Janka 2015).

In the presence of strong seed perturbations, we can express at the onset of an explosive runaway in terms of its value at shock revival in the case of small seed perturbations and a correction factor (1 + ψ)2/3 as in Equation (46),

Equation (50) obviously hinges on the proper calibration (and validation) of Equation (46), which needs to be provided by future core-collapse supernova simulations. Nonetheless, it already allows some crude estimates.

The ratio of the critical luminosity with strong seed perturbation to the critical luminosity value in multi-D for small seed perturbations is found to be

where we linearized in ψ. In order not to rely on an increasingly long chain of uncertain estimates, it is advisable to use the known multi-D effects without strong seed perturbations as a yardstick; they bring about a reduction of the critical luminosity by about 25% compared to 1D (Murphy & Burrows 2008; Hanke et al. 2012; Couch 2013; Dolence et al. 2013; Müller & Janka 2015). This reduction is obtained by setting at the onset of runaway shock expansion, which is also the value derived by Müller & Janka (2015) based on analytic arguments.

Using this value, we estimate a reduction of the critical luminosity by ∼0.15 ψ relative to the the critical luminosity in multi-D without perturbations, which remains only a very rough indicator for the importance of perturbations in shock revival, barring any further calibration and a precise definition of how and where Ma is to be measured.

It is illustrative to express ψ in terms of the heating efficiency ηheat, which is defined as the ratio of the volume-integrated neutrino heating rate and the sum of the electron neutrino and antineutrino luminosities and ,

and the accretion efficiency ηacc,

We then obtain

Using Equation (54), we can verify that the estimated reduction of the critical luminosity due to seed perturbations by ∼0.15ψ is in the ballpark: if we take Ma to be half the maximum value of the Mach number in the infalling shells in the models of Müller & Janka (2015) and work with reasonable average values of ηacc = 2 and ηheat = 0.05, we obtain a reduction of 11% for their model p2La0.25 (Ma = 0.045), 24% for p2La1 (Ma = 0.1), and 36% for p2La2 (Ma = 0.15), which agrees surprisingly well with their inferred reduction of the critical luminosity (Figure 12 in their paper). It also explains why their models with  = 4 require twice the convective Mach numbers in the oxygen shell to explode at the same time as their corresponding  = 2 models. For the models of Couch & Ott (2013, 2015) with ηheat ≈ 0.1 and  = 4, our estimate would suggest a reduction in critical luminosity by 6%. This prediction cannot be compared quantitatively to the results of Couch & Ott (2013, 2015), since an analysis of the effect on the critical luminosity in the vein of Müller & Janka (2015) and Summa et al. (2016) would require additional data (e.g., trajectories of the gain radius). Qualitatively, such a moderate reduction of the critical luminosity seems consistent with their results: the effect of infalling perturbations corresponds to a change in the critical heating factor11 (by which they multiply the critical luminosity to compute the neutrino heating terms) by only 2%–3%, and their inferred reduction of the "critical heating efficiency" by ∼10% due to infalling perturbations is much smaller than the reduction of the critical heating efficiency by a factor of ∼2 in 3D compared to 1D. For the simulations of Couch et al. (2015), for which we estimate the convective Mach number in the silicon shell as roughly 0.02, the expected reduction in the critical luminosity (again for ηheat ≈ 0.1 and a dominant  = 4 mode) is roughly 1%, which is consistent with the development of an explosion in both the perturbed and the unperturbed model.

For our 18 M progenitor model with a typical convective Mach number Ma ≈ 0.1 in the middle of the oxygen shell, we expect a much more sizable reduction of the critical luminosity by 12%–24% if we assume ηacc = 2,  = 2 and ηheat = 0.05–0.1, although this crude estimate still needs to be borne out by a follow-up core-collapse supernova simulation. Because the importance of the infalling perturbations relative to the contribution of neutrino heating to non-radial instabilities is determined by ηacc and ηheat, reasonably accurate multi-group transport is obviously required; the inaccuracy of leakage-based models like Couch & Ott (2013) and Couch et al. (2015) that has been pointed out by Janka et al. (2016) evidently does not permit anything more than a proof of principle.

6. SUMMARY AND CONCLUSIONS

In this paper, we presented the first 3D simulation of the last minutes of oxygen shell burning outside a contracting iron and silicon core in a massive star (ZAMS mass 18 M) up to the onset of collapse. Our simulation was conducted using a 19-species α-network as in the stellar evolution code Kepler (Weaver et al. 1978) and an axis-free, overset Yin-Yang grid (Kageyama & Sato 2004; Wongwathanarat et al. 2010) to cover the full solid angle and allow for the emergence of large-scale flow patterns. To circumvent the problem of core deleptonization and nuclear quasi-equilibrium in the silicon shell without degrading the accuracy of the simulation by serious modifications of the core evolution, a large part of the silicon core was excised and replaced by a contracting inner boundary with a trajectory determined from the corresponding Kepler run. The model was evolved over almost 5 minutes, leaving ample time for transients to die down and roughly 3 minutes or 9 turnover time-scales of steady-state convection for a sufficiently trustworthy analysis of the final phase before collapse.

For the simulated progenitor, an 18 M star of solar metallicity with an extended oxygen shell, our 3D simulation shows the acceleration of convection from typical Mach numbers of ∼0.05 to ∼0.1 at collapse due to the increasing burning rate at the base of the shell. The contraction of the core also leads to the emergence of larger scales in the flow, which is initially dominated by  = 3 and  = 4 modes before a pronounced quadrupolar ( = 2) mode develops shortly before collapse. As a result of a small buoyancy jump between the oxygen and carbon shell, the oxygen shell grows from 0.51 M to 0.56 M due to the entrainment of material from the overlying carbon shell over the course of the simulation, which appears compatible with empirical scaling laws for the entrainment rate at convective boundaries (Fernando 1991; Strang & Fernando 2001; Meakin & Arnett 2007b).

The comparison with the corresponding Kepler model shows that—aside from entrainment at the boundaries—convection is well described by mixing-length theory (MLT) in the final stage before collapse in the model studied here. MLT at least captures the bulk properties of the convective flow that matter for the subsequent collapse phase quite accurately: if properly "gauged," the convective velocities predicted by MLT in Kepler agree well with the 3D simulation, and the time-dependent implementation of MLT even does a reasonable job right before collapse when the nuclear energy generation rate changes significantly within a turnover timescale, which results in a "freeze-out" of convection. The good agreement with MLT is also reflected by the fact that the kinetic energy in convective motions obeys a scaling law of the expected form (Biermann 1932; Arnett et al. 2009). The kinetic energy Ekin,r in radial convective motions can be described to good accuracy in terms of the average nuclear energy generation rate per unit mass , the pressure scale height hP at the base of the shell, and the mass Mconv in the shell as

and the convective velocities are not too far from

where ωBV is the Brunt–Väisälä frequency and hP is the local value of the pressure scale height. Our results are consistent with the assumption that convective blobs are accelerated only over roughly one pressure scale height, and there appears to be no need to replace the pressure scale height in Equation (55) with the extent of the convective zone as the arguments of Arnett et al. (2009) suggest. We surmise that this may be a specific feature of the final phases of shell convection before collapse that requires the dissipation of turbulent energy within a limited distance: since neutrino cooling no longer balances nuclear energy generation, the convective flow will adjust such as to maintain a constant rate of entropy generation throughout the shell to avoid a secular build-up or decline of the unstable gradient. During earlier phases with appreciable neutrino cooling in the outer regions of convective shells, Equation (55) may no longer be adequate.

Similarly, the dominant scale of the convective eddies agrees well with estimates based on linear perturbation theory (Chandrasekhar 1961; Foglizzo et al. 2006). In terms of the radii r and r+ of the inner and outer shell boundary, the dominant angular wave number is roughly

Our findings already allow some conclusions about one of the primary questions that has driven the quest for 3D supernova progenitors, i.e., whether progenitor asphericities can play a beneficial role for shock revival. We suggest that Equations (55) and (57) can be used to formulate an estimate for the importance of convective seed perturbations for shock revival in the ensuing supernova (Couch & Ott 2013; Couch et al. 2015; Müller & Janka 2015). To this end, these two equations (or alternatively, the convective velocities obtained via MLT in a stellar evolution code) need to be evaluated at the time of freeze-out of convection to obtain the typical convective Mach number Ma and angular wave number at collapse. The time of freeze-out can be determined by equating the typical timescale for changes in the volume-integrated burning rate with the turnover timescale tconv, which results in the condition

or

Relying on an estimate for the extra turbulent energy generated in the post-shock region in the supernova core by the infall of seed perturbations, and using the reduction of the energy-weighted critical neutrino luminosity crit for explosion by ∼25% in multi-D (Murphy & Burrows 2008; Hanke et al. 2012; Müller & Janka 2015) as a yardstick, one finds that strong seed perturbations should reduce crit further by

relative to the control value in multi-D simulations without strong seed perturbations. Here ηacc and ηheat are the accretion and heating efficiency in the supernova core, Ma is the typical convective Mach number in the infalling shell at the onset of collapse, and includes the proper weighting of the neutrino luminosity Lν with the square of the neutrino mean energy (Janka 2012; Müller & Janka 2015). This estimate appears to be roughly in line with recent multi-D studies of shock revival with the help of strong seed perturbations and nicely accounts for the range of effect sizes from Couch et al. (2015; no qualitative change in shock revival) to Müller & Janka (2015; reduction of by tens of percent for  = 2 modes with sufficiently strong perturbations). For our 3D progenitor model, we expect a reduction of the critical luminosity by 12%–24%.

Considering these numbers, the prospects for a significant and supportive role of progenitor asphericities in the supernova explosion mechanism seem auspicious. Yet caution is still in order. Because the relative importance of seed perturbations is determined by the ratio Ma/(ηaccηheat), a reliable judgment needs to be based both on a self-consistent treatment of convective burning in multi-D before collapse (which determines Ma) and accurate multi-group neutrino transport after bounce (which determines ηacc and ηheat). Again, first-principle models of supernovae face a curious coincidence: as one typically finds Ma ∼ ηheat, progenitor asphericities are just large enough to play a significant role in the explosion mechanism, but not large enough to provide a clear-cut solution for the problem of shock revival. The danger of a simplified neutrino treatment has already been emphasized repeatedly in the literature (see Janka et al. 2016 for a recent summary), but pitfalls also abound in simulations of convective burning: for example, the recipes employed by Couch et al. (2015) can be shown to considerably affect the convective Mach number at collapse by appealing to scaling laws from MLT and timescale considerations.

Our method of excising the core seems to be a viable avenue toward obtaining 3D initial conditions in the oxygen shell (which is the innermost active convective shell in many progenitor models) without introducing inordinate artifacts due to initial transients or artificial changes to the nuclear burning. Nonetheless, the model presented here is only another step toward a better understanding of the multi-D structure of supernova progenitors. In particular, the effects of resolution and stochasticity on the convective flow need to be studied in greater depth, though a first restricted resolution study (Appendix B) suggests that the predicted convective velocities are already accurate to within 10% or less.

Future simulations will also need to address progenitor variations in the shell geometry, shell configuration, and the burning rate; in fact, the 18 M was deliberately chosen as an optimistic case with strong oxygen burning at the base of a very extended convective shell, and may not be representative of the generic situation (if there is any). Moreover, massive stars with active convective silicon shells at collapse also need to be explored even if they form just a subclass of all supernova progenitors. Treating this phase adequately to avoid the artifacts introduced by an approach like that of Couch et al. (2015) is bound to prove a harder challenge due to the complications of nuclear quasi-equilibrium. Finally, the long-term effects of entrainment and other phenomena that cannot be captured by MLT need to be examined: if such effects play a major role in the evolution of supernova progenitors, capturing them with the help of exploratory 3D models and improved recipes for 1D stellar evolution in the spirit of the 321D approach (Arnett et al. 2015) will be much more challenging than 3D simulations of the immediate pre-collapse stage, where the problems of extreme timescale ratios (e.g., of the thermal adjustment and turnover timescale), numerical diffusion, and energy/entropy conservation errors are relatively benign. It is by no means certain that supernova progenitor models will look fundamentally different once this is accomplished; but there is little doubt that groundbreaking discoveries will be made along the way.

We thank T. Foglizzo, E. Müller, and S. Woosley for useful discussions and T. Melson for support and discussions concerning the Yin-Yang grid. We acknowledge support by the Australian Research Council through a Discovery Early Career Researcher Award DE150101145 (B.M.) and an ARC Future Fellowship FT120100363 (A.H.) by the Deutsche Forschungsgemeinschaft through the Excellence Cluster universe EXC 153 (T.J.) and by the European Research Council through grant ERC-AdG No. 341157-COCO2CASA (M.V., T.J.). This research was undertaken with the assistance of resources 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. This material is based upon work supported by the National Science Foundation under grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).

APPENDIX A: LEDOUX CRITERION IN TERMS OF ENTROPY AND COMPOSITION GRADIENTS

The Ledoux and Schwarzschild criteria for convective instability are often expressed in terms of entropy and composition gradients in the stellar evolution literature. It is straightforward to show that Equation (1) can equally be expressed in terms of these gradients:

APPENDIX B: EFFECTS OF RESOLUTION AND STOCHASTICITY

Limited resolution is always a cause of concern for turbulent flow as in convective shell burning. Current models of shell burning during the final stages like Couch et al. (2015) and ours are limited to ∼1.5 million zones to cover the region of interest, and a significant increase in resolution is currently not affordable. To gauge potential resolution effects, we have, however, computed two additional runs with lower resolution, and compare key quantities of those models to the baseline run. We both consider the case of a reduced angular resolution of 3° with the same radial grid as well as a simulation with a resolution of 3° and a coarser radial grid of 286 zones (with Δr/r = 1%).

The two complementary low-resolution simulations also help to address another important effect to a limited extent, i.e., stochastic variations in the convective flow geometry. With just three different models computed at different resolutions, it is obviously impossible to distinguish cleanly between the effects of resolution and stochasticity. The three models merely define a vague "band of uncertainty," and only permit limited conclusions about the underlying cause for the variations between the models.

While all three models exhibit very similar nuclear energy generation rates, the dynamics of the convective flow is slightly different. Figure 18 compares the efficiency factor ηconv (Equation (14)) for the conversion of nuclear energy into the components Er and of the turbulent kinetic energy (which is essentially the same as comparing Er and directly because of the extremely similar burning rate). ηconv does not differ substantially for the two low-resolution runs. Variations between the different simulations are typically on the level of 10% or less, i.e., of the same order as the stochastic fluctuations within each run, and there is thus no clear evidence for a strong resolution dependence. This is in line with findings from a different context (convection in core-collapse supernova explosions), where Handy et al. (2014) and Radice et al. (2016) found that the global energetics of the flow are well captured, even with a modest resolution of 2°–3° in angle and ∼100 radial zones across the convective region.

The non-radial kinetic energies in the low-resolution runs differ more strongly from the baseline run, especially after 190 s. It is not clear, however, whether this is a resolution effect or if it is due to stochastic variations. The higher non-radial kinetic energy in the baseline run appears to be connected to a slightly different eddy geometry, which could suggest stochasticity rather than differences in resolution as the culprit for the differences between the runs. Figure 19 shows the coefficients for the decomposition of the radial velocity into spherical harmonics (computed according to Equation (30)) at a time of 210 s and demonstrates a slightly stronger preponderance of low- modes in the baseline run compared to the low-resolution models. It is thus conceivable that the lateral flows from the updrafts to the downdrafts at the top and bottom of the convective shell must be slightly faster in the baseline run merely to ensure mass conservation in a flow that is (almost) anelastic: in the baseline run, the lateral flow must transport mass from the updrafts to the downdrafts (and vice versa at the bottom of the shell) at the same rate as in the low-resolution runs, but over a larger distance; hence, the ratio of non-radial to radial turbulent velocities is larger.

Figure 19. Refer to the following caption and surrounding text.

Figure 19. Power in different multipoles for the decomposition of the radial velocity at r = 4000 km into spherical harmonics at a time of 210 s for the baseline run with an angular resolution of 3° degrees and Nr = 400 radial zones (black curve) and for two low-resolution runs with an angular resolution of 2° and a radial resolution of 400 zones (red) and 286 zones (blue). Although the dominant angular wave numbers are similar, the baseline run shows more power in  = 2 and  = 3 modes.

Standard image High-resolution image

Whatever the underlying cause, the lesson of our minimal resolution study is that we should anticipate uncertainties in the convective velocities, which scale with , of ≲10%, and that the dominant angular wave number of the convective eddies may also vary slightly.

Even the growth of the oxygen shell by entrainment is not too dissimilar for the three different runs, as shown by Figure 20. Again, it is unclear to what extent the variations in entrainment are stochastic or due to resolution effects. The resolution requirements for the problem at hand may be mitigated by the "softness" of the convective boundary: interfacial waves develop on relatively large scales, so that the individual breaking billows are typically covered by ∼10 zones or more. This is still well below the resolution of almost 50 zones in the transition of the best-resolved global simulations of convective boundary mixing in the context of hydrogen ingestion (Herwig et al. 2014; Woodward et al. 2015), and hence, more caution is in order concerning the convergence of the entrainment rate than for the convective velocities and eddy scales.

Figure 20. Refer to the following caption and surrounding text.

Figure 20. Growth of the mass Mconv of the oxygen shell as a function of time for the baseline run with an angular resolution of 3° degrees and Nr = 400 radial zones (black curve) and for two low-resolution runs with an angular resolution of 2° and a radial resolution of 400 zones (red) and 286 zones (blue). Note that Mconv starts at a slightly different value in the run with Nr = 286 because the initial model is mapped to a different grid than for Nr = 400.

Standard image High-resolution image

Footnotes

  • Note the sign convention used in this paper: corresponds to convective instability.

  • We do not attempt to classify the precise type at instability at play, since this is immaterial for our purpose.

  • Note that ηconv does not correspond to the "convective efficiency" as often used in stellar evolution, i.e., it is not the ratio of the convective luminosity to the radiative luminosity.

  • 10 

    If we assume that the density perturbations in the post-shock region adjust on a dynamical timescale as pressure equilibrium between over- and underdensities is established, then this estimate might be lower, but pressure adjustment itself would involve the generation of lateral flows and hence generate turbulent kinetic energy, so that our estimate is probably not too far off.

  • 11 

    The change in the critical heating factor is related to but not necessarily identical to the change in the generalized critical luminosity as introduced by Müller & Janka (2015) and Summa et al. (2016), which also depends, e.g., on the relative change of the gain radius and the specific binding energy in the gain region.

10.3847/1538-4357/833/1/124
undefined