Idealised hydrodynamic simulations of turbulent oxygenburning shell convection in geometry
Abstract
This work investigates the properties of convection in stars with particular emphasis on entrainment across the upper convective boundary (CB). Idealised simulations of turbulent convection in the Oburning shell of a massive star are performed in geometry on and grids, driven by a representative heating rate. A heating series is also performed on the grid. The simulation exhibits an entrainment rate at the upper CB of . The simulation with the same heating rate agrees within 17 per cent. The entrainment rate at the upper convective boundary is found to scale linearly with the driving luminosity and with the cube of the shear velocity at the upper boundary, while the radial RMS fluid velocity scales with the cube root of the driving luminosity, as expected. The mixing is analysed in a 1D diffusion framework, resulting in a simple model for CB mixing. The analysis confirms previous findings that limiting the MLT mixing length to the distance to the CB in 1D simulations better represents the sphericallyaveraged radial velocity profiles from the 3D simulations and provides an improved determination of the reference diffusion coefficient for the exponential diffusion CB mixing model in 1D. From the 3D simulation data we adopt as the convective boundary the location of the maximum gradient in the horizontal velocity component which has spatial fluctuations of . The exponentially decaying diffusion CB mixing model with reproduces the sphericallyaveraged 3D abundance profiles.
keywords:
stars: massive, evolution, interior — physical data and processes: turbulence, hydrodynamics, convection1 Introduction
Modelling the longterm evolution of stars requires the assumption of spherical symmetry and hydrostatic equilibrium. Such models are crucial for predicting the characteristics of stars with a given initial mass, metallicity, composition and rotation rate, which are important for many areas of astronomy (e.g. distance, age and mass determinations, population synthesis, galactic chemical evolution, etc.). The predictive power of these models is somewhat restricted by their dependence on approximate and often parametrised treatment of underlying physical processes that become necessary when enforcing symmetries and equilibria upon the models. Two such processes are convective mixing and the mixing at convective boundaries, for which either an instantaneous mixing or diffusion approximation is typically used (e.g. Freytag, Ludwig & Steffen, 1996; Herwig et al., 1997; Eldridge & Tout, 2004; Young et al., 2005; Eggenberger et al., 2008; Brott et al., 2011; Limongi & Chieffi, 2012; Ekström et al., 2012; Rauscher et al., 2002; Woosley & Heger, 2007; Paxton et al., 2013, 2015). A key quantity of interest is the rate at which material is entrained into convection zones, that is the rate at which material is transported across the convective boundaries from neighbouring stable layers by hydrodynamic instabilities in the vicinity of the convective boundary.
The varied ways in which these approximations are formulated and then implemented into stellar evolution codes result in different codes making different predictions for, e.g., nuclear burning lifetimes, core masses, evolution in the Hertzsprung–Russell diagram, weak process production (e.g. Martins & Palacios, 2013; Jones et al., 2015). In the advanced burning stages of massive stars (Cburning onwards), the structure of the star is generally said to be frozenin. However, as Sukhbold & Woosley (2014) have shown, seemingly insignificant changes in the masses of convective regions within the stellar core during the advanced burning stages can have a marked effect on the interior core structure. These changes to the core structure can shift the compact remnant that is predicted to be produced when the star dies, from a neutron star to a black hole (Ertl et al., 2016). The former case would typically result in a corecollapse supernova of type II^{1}^{1}1however, note that note that a subset of type I supernovae (Ib and Ic) are also corecollapse supernovae, while the latter case (black hole formation) is predicted to produce a longduration ( year) low luminosity ( erg s) transient often referred to as a failed supernova (Lovegrove & Woosley, 2013). To date, the only progenitor stars of type II supernovae that have been detected fall in the luminosity range , corresponding to an initial progenitor mass of (Smartt, 2015). However, red supergiant stars have indeed been observed with luminosities in the range (Levesque et al., 2006, 2009) and from a statistical point of view based on the current theory of stellar evolution, stars this luminous/massive should have contributed to the number of direct detections mentioned previously. This suggests that red supergiants with () always result in weak or failed supernovae that form black holes. This is a picture that is not, however, corroborated by recent spherically symmetric supernova simulation efforts (O’Connor & Ott, 2011; Ugliano et al., 2012; Ertl et al., 2016; Sukhbold et al., 2016) although the explosion mechanism of corecollapse supernovae is still not completely understood, and is certainly strongly influenced by asymmetries in the structure of the progenitor star and in the explosion itself (see Müller, 2016, for a recent review). With the outcome of such simulated explosions being so intrinsically linked to the structure of the progenitor star, the accuracy of stellar models is arguably as important as that of the supernova simulations when distinguishing between a failed (black holeforming) or successful (neutron starforming) supernova explosion for a star of given initial mass, metallicity and rotation rate. The mass of the compact remnant produced during core collapse is also strongly related to the presupernova stellar structure, and different predictions for preSN structures as a function of initial stellar mass and metallicity can lead to rather different conclusions. A timely example is the metallicity upper limit placed on the progenitor system of the binary black hole merger whose detection was the first of its kind using gravitational wave telescopes (Abbott et al., 2016b). The binary population synthesis code StarTrack placed a limit of (Belczynski et al., 2016) while the limit from another such code BPASS was (Eldridge & Stanway, 2016), in agreement with the prediction of Abbott et al. (2016a) which was obtained independently. Binary population synthesis models, however, due to their very nature, inherit the significant uncertainties of stellar models, and so to better constrain stellar models is also to improve the predictions of population synthesis models, through which the stellar models may be validated in the future with further gravitational wave measurements from compact binary mergers.
The structure of the core during the advanced burning stages of massive stars is determined by many different aspects, which are themselves products of other physical processes that have in some cases been described in a stellar model by their own parametrisation. For example, the C/O ratio in the core at the onset of C burning is a major factor setting the mass of the convective Cburning core and hence the location of the first subsequent Cburning shell. This ratio is set by the competition between the triplealpha (3) and the C(O reactions at the end of the core Heburning phase, which are sensitive to both the number of particles and the temperature. The particles are introduced into the core at different rates depending upon the treatment of mixing at the boundary of the convective Heburning core, and the temperature is set by the mass of the Hdepleted core, which is itself determined by the treatment of convective boundary mixing (CBM) during the main sequence (core Hburning). So, the core structure of a massive star during O shell burning, for example, is not only determined by the mixing assumptions in the O shell but in every convective episode prior.
Asteroseismological measurements of mixed modes in core Heburning stars provide a constraint for the extent and behaviour of convective (boundary/overshoot) mixing in stellar models (e.g. Constantino et al., 2015; Bossini et al., 2015). However, it is rather optimistic to expect that the same diagnostics will be measured for an Oburning shell in a 25 star. What can be done, however, are multidimensional hydrodynamic simulations of the Oburning shell that describe the convective properties very well on short timescales (Arnett, 1994; Bazán & Arnett, 1998; Asida & Arnett, 2000; Young et al., 2005; Meakin & Arnett, 2007). These kinds of simulations are much more difficult during core Heburning because of the much lower Mach number of the flow and the much longer evolutionary time scales.
Mixing processes at the convective boundaries of advanced (C, O and Si) burning shells that can be studied in multidimensional hydrodynamic simulations (see, e.g., Arnett, 1994; Bazán & Arnett, 1998; Asida & Arnett, 2000; Young et al., 2005; Meakin & Arnett, 2007) are known to have an impact on the presupernova structure of corecollapse supernova (CCSN) progenitor models (Young et al., 2005; Frey, Fryer & Young, 2013; Sukhbold & Woosley, 2014), which in turn affects the dynamics of the supernova explosion via an alteration of the competition between the ram pressure of the infalling core material with the neutrinoheated material behind the stalled shock. Because the time scale of the collapse is much shorter than the time scale of convection in the final O shellburning episode, the imprint of the convective velocity field is essentially frozen in. Fluctuations in the velocity amplitudes can impact the success of the simulated corecollapse supernova explosion of the progenitor star (Couch & Ott, 2013; Müller & Janka, 2015; Couch et al., 2015).
Nucleosynthesis in massive stars is also sensitive to the convective structure of the core during the advanced burning stages. Rauscher et al. (2002), reporting on their computations of the complete nucleosynthesis in massive stars, highlighted an event in a 20 model in which the convective Oburning shell and the convective Cburning shell merged into one deep convective layer, engulfing in addition the Neburning shell. In their computations, this event was responsible for a latetime neutron burst that resulted in large overproduction of neutroncapture nuclei. The impact or likelihood of such a shell interaction has not yet, to our knowledge, been studied or reported on in any great detail. However, Bazán & Arnett (1994, 1998) and Asida & Arnett (2000) have shown that the entrainment of from the overlying stable layer into the convective Oburning shell produced hot spots of nuclear burning which do feed back and affect the flow. Meakin & Arnett (2006) and Arnett & Meakin (2011), performing 2D simulations of multiple stratified convective burning shells, reported on the interaction of the shells by wave propagation in the intershell cavity.
The temporally and spatially stochastic fluctuations in the energy generation rate due to the burning of entrained fuel as it is advected toward the flame at the base of the convection zone in the simulations of Arnett & Meakin (2011) results in an enhanced entrainment rate of the fuel. The enhanced entrainment rate boosts the nuclear energy generation rate and towards the end of the calculation, large sloshing motions appear. A similar phenomenon was reported by Herwig et al. (2014) in a PPMstar simulation of H ingestion into the Heburning shell flash during a very late thermal pulse (convective He shell flash; VLTP) of a postAGB star.
In this work, idealised simulations of convection in the first Oburning shell of a massive star are performed with the PPMstar code in geometry. The original motivation of this work was to investigate the numerical convergence properties of mass entrainment at the top convection boundary. A similar investigation for Heshell flash convection was performed by Woodward, Herwig & Lin (2015), who found a very low entrainment rate at a seemingly very stiff convective boundary. In that work it was shown that the numerical simulations approach convergence for this property. One goal of this paper is therefore to establish convergence properties for the entrainment rate for this alternative setup. Following encouraging initial results, the analysis was extended to develop a mixing model in the diffusion framework that describes convective mixing in the region of the upper boundary of Oshell convection.
The concept and setup of the idealised simulations are described in Section 2 along with the stellar evolution calculations on which they are based. The results of the simulations are described in Section 3, including a 1D diffusion analysis connecting the 3D hydrodynamic simulations back to 1D stellar models. Section 4 contains a comparison to other works. In Section 5 the results are summarized and a brief outlook is given.
2 Simulations
Two types of simulations of the first convective Oburning shell were performed in the present work: 1D stellar evolution calculations using the MESA code (Paxton et al., 2011, 2013, 2015) and explicit 3D hydrodynamic simulations using the PPMstar code (Woodward, Herwig & Lin, 2015). As is described in the following sections, the MESA models boast an impressively detailed picture of the microphysics while the PPMstar simulations are performed with idealised physics assumptions. These assumptions are but one of the reasons that it is possible to use very high numerical resolution in geometry with the PPMstar code and still follow the simulation for 27 minutes with the highresolution grid and up to 55 minutes using our lowresolution grid. Additionally, we would like to be able to compare the O shell simulations directly with the AGB thermal pulse entrainment simulations by Woodward, Herwig & Lin (2015), and this requires being able to separate the effects of macrophysics and microphysics. The latter will be investigated in future work. This compromise on the details of the microphysics means that the initial setup of the PPMstar simulations will not exactly match the spherically symmetric MESA models that they are based on. It is of course reasonable to expect that the 3D hydrodynamic and 1D hydrostatic models will differ in many aspects, and in the present work it is assumed that the important aspects to optimize for are the aspect ratio and driving luminosity of the convection zone, and the mean molecular weights of the fluids interior and exterior to the convection zone. Furthermore, when comparing 1D and 3D simulations, it can actually become quite challenging to determine metrics by which the two can actually be compared (see Section 3.3). We define the aspect ratio of a convection zone to be , where and are the radii of the upper and lower boundary of the convection zone, respectively. While perhaps it may seem more intuitive to have the depth nondimensionalised by putting in the denominator, such a definition would be problematic for convective cores for which .
2.1 MESA simulations
Stellar evolution calculations of a nonrotating massive star with metallicity were performed using the MESA code^{2}^{2}2The initial chemical composition of the stellar models was , and with the metal distribution after Grevesse & Noels (1993). This stellar model is a continuation of the m25 model from Jones et al. (2015) beyond core heliumburning. The MESA model assumes hydrostatic equilibrium and uses the mixing length theory of convection (MLT) with mixing length , where is the pressure scale height. Convective stability is determined on the basis of the Schwarzschild criterion and convective mixing of nuclear species is modelled as a diffusive process. Inside the convection zone, the diffusion coefficient is assumed to be , where is the convective velocity predicted by MLT. Mixing at the convective boundaries of the oxygenburning shell was approximated with an exponentiallydecaying diffusion coefficient as proposed by Freytag, Ludwig & Steffen (1996, see also , ). The convective boundary mixing parameter used in this approximation was , which is used to give the efolding length of the diffusion coefficient, . The nuclear species included in the nuclear reaction network in the MESA simulations performed as part of this work are given in Fig. 1.
A Kippenhahn (convective structure evolution as a function of time) diagram of the MESA model with positive nuclear energy generation contours is shown in Fig. 2. The figure shows only the central region of the star (the innermost ) during the core and shell Oburning phases. A deep overlying convective C shell is present throughout, although only the bottom portion can be seen in the figure. The initial setup of our PPMstar simulations is based on the structure of the convective O shell in this MESA model at , i.e. roughly 10 days before core collapse. The Oburning shell in the stellar model sits atop a convectively stable silicon core of about in radius (about ). The shell is convective and is driven by the input of heat into the gas, which is predominantly from the fusion of nuclei with a small contribution made by other nuclear reactions and gravitational contraction. The majority of this energy generation is confined to a thin shell located at a radius of (about ) where both the temperature and the concentration of fuel are high enough for fusion to occur. The convection zone extends out from the shell at to about (). Above the convective oxygenburning shell is a deep convective Cburning shell. The two convective shells are separated by a stable layer of roughly (about ) of material, mostly composed of O, Ne, Mg and a small amount of Si (i.e. products of incomplete Ne burning). The aspect ratio of the convective Oburning shell is thus about 0.5. In a survey performed as part of this work, massive star models were computed with the MESA code with initial masses ranging from 12 to . The aspect ratio of this first Oburning shell is quite similar in all of the models across the range of initial masses that were simulated.
The dominant source of pressure in the convective Oburning shell comes from the gas, which contributes about 75 per cent of the total pressure, while degenerate electrons contribute at the per cent level. Radiation pressure accounts for the remaining fraction of the total pressure (Fig. 4). The contributions in the overlying stable layer are similar to the Oshell, but in the convective Cburning shell (with its base at a mass coordinate of about ), the contribution to the pressure is split equally between the gas and radiation.
Several mechanisms lead to the rapid production of neutrinos in the deep interior of massive stars during their advanced burning stages (see, e.g., Itoh et al., 1996). The densities are, however, low enough that the interaction of the neutrinos with the stellar plasma is negligible and they stream freely from the stellar core. The production of neutrinos accelerates the evolution of massive stars during their postHe core burning phases: nuclear binding energy must be released at a rate that not only provides the necessary luminosity to support the star but also at a rate that compensates for the neutrino energy losses. Radial profiles of the nuclear energy generation rate , thermal neutrino loss rate (i.e. neutrinos produced by the plasma via pair production, which is by far the dominant source of neutrinos in the Oburning shell) and luminosity in the vicinity of the oxygen shell are shown in Fig. 3. Despite the rapid neutrino energy sink, the net power generated due to the nuclear burning is still some ( erg s).
2.2 PPMstar simulations
Hydrodynamic simulations of Oshell convection were performed in 3D with geometry using the PPMstar code. The code as well as its performance for a similar application (Heshell flash convection in a lowmass star) are described in detail in Woodward, Herwig & Lin (2015). The setup procedure, as well as the physics assumptions in this work are the same as in the Heshell flash convection simulations. In summary, an ideal gas equation of state is adopted and a geometrically representative stratification is produced by combining three piecewise polytropes. Particular attention is paid to reproducing the pressure and density stratification. These are the important quantities for describing the flow; the temperature, if it were considered, would be overestimated since the contribution of the radiation to the total pressure (see Fig. 4) is assumed to be provided solely by the gas. In Fig. 5 defining quantities of the convective Oburning shell in the 25 MESA model – the mean molecular weight in the convection zone , the radial coordinate of the convective boundaries and , and the peak luminosity inside the convection zone – are given as a function of time until the onset of core collapse (cf. Fig. 2). This figure is intended to demonstrate the similarity of the initial setup of the PPMstar simulations presented in this work to the structure one would expect the Oburning shell to have when much more detailed microphysics have been considered.
id  grid  

(min)  ()  ()  (km/s)  (km/s)  (km/s)  
d1  55.2  11.8  1.15  0.016  32.0  38.0  1.21  
d2  27.3  11.8  1.33  0.018  32.6  36.8  1.19  
d8  36.2  29.5  3.60  0.035  44.4  50.6  1.62  
d5  37.2  59.1  8.07  0.060  53.5  64.1  2.08  
d6  41.3  118  16.8  0.116  69.4  87.7  2.82  
d9  43.5  295  38.3  0.257  92.6  116  3.77  
d10  43.7  591  79.4  0.524  123  151  4.85 
grid resolution  simulated time  driving luminosity 
entrainment rate  upper boundary velocity  maximum radial velocity 
maximum tangential velocity  maximum Mach number 
The convective shell is represented by a polytropic stratification with . The layers above and below have stable stratifications with and , respectively. The convection is driven by a constantvolume heating rate of , which is equivalent to the rate of energy generation due to oxygen burning in the MESA model to within a factor of two (see Fig. 5 and Section 2.1) accounting for the neutrino energy losses from both nuclear reactions and pair production. The heated shell has a thickness of cm (Fig. 3) and the transition region between the top of the convective shell and the overlying stably stratified layer is cm (about ), centred at about cm. The thickness of the transition layer was chosen to give the best agreement between the profiles of and in the piecewise polytropic PPMstar setup and the MESA model on which it was based. It is important to also note that the thickness of the transition layer should be such that it can be adequately resolved on both the and resolutions grids of PPMstar. At the base of the convective shell, the density and gravitational acceleration are g cm and , respectively. At the top convection boundary a pressure scale height is .
At the start of the simulation, the stable layer consists of fluid with a mean molecular weight , and the rest of the simulation domain is filled with fluid with a mean molecular weight . The resulting profile of the squared BruntVäisälä frequency^{3}^{3}3This expression is equivalent to the more familiar expression for an ideal gas, in which . as compared with that given by MESA is shown in Fig. 6. The MESA model in Fig. 6 matches the position of the convective boundary in PPMstar, although the mean molecular weight of the convective fluid in the MESA model is slightly too low at . The profile in a somewhat more evolved MESA model (not shown in Fig. 6 for clarity) that matches of the PPMstar model has a similar structure and amplitude, but the top of the convection zone in that model is located at a radius larger than in the PPMstar model.
The Oburning shell convection setup is very similar in important respects to the VLTP convection in Sakurai’s object, which has also been simulated using PPMstar (Herwig et al., 2011, 2014). The aspect ratios (see the beginning of Section 2 for a definition) of the Oburning convective shell and the VLTP shell flash convection zone are 0.49 and 0.69, respectively. Mach numbers in the convection flows are also similar. An appreciable difference is the ratio of the mean molecular weight of the convective fluid to that of the overlying stable fluid. The VLTP in Sakurai’s object has such a mean molecular weight ratio of 2.26, while in the Oburning shell convection problem the same ratio is only 1.02. Another difference is that in these Oshell simulations we do not take into account the burning of any entrained material from above. In that regard the simulation approach and goals are the same as those of the entrainment simulations for Heshell flash convection of Woodward, Herwig & Lin (2015). In this paper we wish to study the entrainment and general flow properties without the additional effect of possible feedback from nuclear burning of the entrained material. A summary of the 3D simulations performed for this work is given in Table 1.
3 Results
3.1 General properties of the flow
The 3D simulation begins with a constant entropy in the convection zone that is driven into an unstable convection flow by introducing a continuous injection of heat with a spherically symmetric profile near the base of the convection zone (see Fig. 3). Rather than stir the flow initially with a spectrum of disturbances, we instead allow the very slight differences in the numerical representation of the flow on our fine Cartesian grid to provide the initial perturbations that permit some of the gas to rise and some to sink in the convection zone. This approach is not only convenient, but it also provides a natural way to test that the grid that we use does not ultimately affect the disturbance spectrum that develops, except of course to truncate it at its short wavelength end. This same approach was used in Woodward, Herwig & Lin (2015), where visualisations near the beginning of the simulation (e.g. their Fig. 6) illustrate the initial grid imprint effects as well as the process by which they ultimately are overwhelmed by the flow’s natural development of its turbulent convection spectrum. One might think that waiting for the initial perturbations to be overwhelmed in this way could be shortened by beginning with a full spectrum of disturbances, but we would argue that one must wait this same length of time in the simulation in any case to be sure that the spectrum is the one the star prefers rather than just the one that has been put in initially by fiat.
In the simulation here (d2), the imprint of the grid on the properties of the flow persists at a detectable level through to about 1300 s (21.7 min), by which time the flow seems to have been fully turbulent already for a considerable period. In this setup grid imprints on the disturbance spectrum persist for a longer time compared to the Heshell convection simulations. The reason seems to be that the aspect ratio of the convection zone (see Section 2), about one half, prefers the largest convection cells that correspond to the Cartesian grid’s natural mode. This can be seen at early times by the prominence of 4 large convection cells in thin sections taken through the center of the star and aligned with planes of the grid (see movies at URLs given in Section 5). These develop at the outset and gradually fade in importance. Ultimately a still larger mode with just 3 rather than 4 cells visible in such sections develops. This mode cannot possibly be preferred by the Cartesian grid, and its development and ultimate dominance indicates that the flow has taken on a character by 1300 s that is determined by the physics and the properties of the stellar convection setup rather than by our numerical treatment. We choose to use a Cartesian grid because this makes the design and execution of the simulation code fit the properties of the machine especially well. The ultimate overwhelming of grid imprints indicates that our grid is fine enough that we are justified in designing the numerics to fit the machine rather than to fit the geometry of the star. This analysis of the evolution and eventual disappearance of grid imprints complements the analysis of simulation properties under grid refinement (Section 3.4) and provide some level of code verification for this problem.
The character of the flow above the Oburning shell in this problem is quite similar to that which we have seen in the shellflash convection zone during the very late thermal pulse (VLTP) of postAGB star Sakurai’s object (Herwig et al., 2014) and that of a thermal pulse in an AGB star (Woodward, Herwig & Lin, 2015). The most significant difference from a fluid dynamics perspective is the relatively large entrainment rate of gas from above the convection zone that results here for reasons that are explained in Section 3.4. This is illustrated in Fig. 7. Fractional volume profiles of the entrained fluid are shown for the O shell and AGB thermal pulse simulations at equivalent computational effort, which is defined here as the product of the number of time steps and the Courant number. The much greater amount of material that has been entrained into the Oburning shell convection zone compared to the AGB He shell convection zone after expending the same amount of computational effort makes the Oburning shell better suited for studying the properties of convergence in PPMstar simulations of stellar convection. Otherwise, in all these three cases (VLTP, AGB thermal pulse and now O shell) we have a dominance of very large convection cells, demanding a global treatment, and low, but not too low, flow Mach numbers, as has been pointed out earlier. In the left panel of Fig. 8, we see a hemisphere of the star, with the hemisphere in the foreground cut away. The plane facing us is the plane in the grid of the simulation. The gas of the degenerate core has been made transparent, and only mixtures of convection zone gas with entrained gas from above the convection zone are visible. These mixtures at the bottom of the convection zone have been made to appear dark blue and semitransparent, so that the location of this spherical surface is apparent. High concentrations of entrained gas are red, and as the concentrations become smaller colors go from red to yellow, white, turquoise, and dark blue, with differing degrees of opacity. The pervasive low concentrations of entrained gas in the convection zone have been made transparent, so that the convection cell pattern and downward moving entrainment lanes near the top of the convection zone can be easily seen. In this hemisphere, there are clearly three large areas where updrafts in large convection cells scrape against the bottom of the stably stratified layer to reveal the red and yellow colors of what is a relatively sharp transition between convection zone and stably stratified gases. There are also three large regions of descending gas that include higher concentrations of entrained gas from above the convection zone. These appear as turquoise and blue. Where this layer of entraining gas is cut by our slicing plane at , its internal strong gradient of entrained concentration can be made out. As discussed in Woodward, Herwig & Lin (2015), there is resistance to entrainment where the gas of the convection zone flows roughly tangentially to the top surface of the convection zone. Those regions show up as red and yellow in this image. In the regions that show up as turquoise and blue, convection cells meet and the flow is forced to descend. It is here that the entrainment happens. It is resisted by the relative buoyancy of the entrained gas, but its concentration is not large, and this resistance is overcome even by the low Mach numbers that apply in this region.
The image on the right in Fig. 8 is the concentration of entrained gas in the same plane at the same time in the simulation, but the rendering has been altered to show the entrained concentration not only near the top of the convection zone but all the way down to its bottom. This is a thin slice, just 1 per cent of the diameter of the simulation volume in thickness. It is clear from this image that globs of relatively enriched gas descend right to the bottom of the convection zone in the very large, dominant convection cells that develop due to the convection zone’s large depth as a fraction of the radius of its top surface. In the real star this entrained material would of course be burned by nuclear reactions, which we ignore in this set of simulations. The active entrainment near this top surface can also be seen in the KelvinHelmholtz trains of vortices that develop where the flow near the top of the convection zone turns downward, as was discussed in detail in Woodward, Herwig & Lin (2015). These trains of vortices show up mainly as red and yellow features in this image. The images in Figs. 9 and 10 are renderings of the vorticity and the radial fluid velocity, respectively, in this same thin slice through the star. The vorticity image emphasizes the smallest vortices and shear layers, making clear that a full, turbulent spectrum of vortices has developed by this time. The radial velocity image shows the three major upwellings in shades of red and yellow, while descending gas appears in shades of grey tending to blue in the most rapidly descending plumes.
3.2 Velocity profiles
The sphericallyaveraged RMS velocity profiles in the d2 run are shown in Fig. 11. The typical RMS total velocity is km s with about equal contributions from the radial and tangential components at Mm where the radial component reaches a maximum. The broad local maxima of the tangential velocity correspond to the fact that the speed becomes mostly tangential near the convection zone boundary as the flow turns around. While both the tangential and the radial velocity profiles have prominent features the total velocity profile is rather flat. Velocity profiles in the rest of the runs listed in Table 1 have the same structure. However, the detailed shape of profiles depends on the shape of the heat source, see the discussion in Sect. 4. Fig. 12 shows that all of the runs reach a steady state after a short initial transient.
The flow velocity is closely related to the luminosity of the heat source. The bestfit scaling relation shown in Fig. 13 is very close to , in agreement with Porter & Woodward (2000). This scaling can be motivated as follows.
In a stationary state, the luminosity can be decomposed to the fluxes and of enthalpy and kinetic energy, respectively,
(1) 
because we neglect the small radiative contribution. Both terms on the right hand side of Eq. 1 can be shown to scale in proportion to . The kinetic energy flux scales with the kinetic energy density and velocity, hence . The enthalpy flux is , where is a typical temperature fluctuation between the upflows and downflows. Assuming that the relative temperature fluctuations are of the same order as the relative dynamic pressure fluctuations , we have , because by the ideal gas law. The enthalpy flux is then .
Another way to see why is to consider the kinetic energy generation and dissipation rates in a convection zone. Per unit of volume, the turbulent dissipation rate is , where is an integral length scale of the turbulence. Therefore also the total dissipation rate in the convection zone, which in a stationary state is of the same order as (see Viallet et al., 2013), scales in proportion to . A similar argument was made by Biermann (1932).
3.3 Where is the convective boundary?
Standard stellar evolution models, such as our MESA models, are spherically symmetric. The local Schwarzschild or Ledoux criterion is used to determine convectively unstable regions. The mixinglength theory (MLT) is then applied to these regions to estimate the temperature gradient, convective velocity and convective energy flux. According to MLT, the convective velocity vanishes at convective boundaries.
In this section we will describe how convection in the 3D simulations presented here differs from the 1D picture. The entropy gradient averaged spherically and over a few overturning time scales is weakly positive in the upper half of the convection zone. If interpreted as a 1D representation of the stratification it would imply this upper part to be formally stable. However, convection there is just as vigorous as that in the lower part of the convection zone (see Figs. 9, 10, and 11). This suggests that instead of the average stratification the velocity field itself should be used to delineate the boundaries of a 3D convection zone. Figure 11 shows that the tangential velocity component reaches a local maximum both near the bottom and near the top of the convection zone and then drops suddenly at the convective boundary, where the steep entropy gradient forces the flow to turn over. The flow velocity does not vanish completely because of the presence of internal gravity waves and sound waves in the stable stratification. We adopt the convective boundary as the location where the decline in is the steepest.
This criterion can be applied to the spherical average as well as to different radial directions. This shows that in contrast to the 1D MLT picture, the boundary of a 3D convection zone is not located on a perfect sphere. The variation of the location of the convective boundary in different directions translates into finite thickness of the boundary when spherical averages are formed. To illustrate this, the data has been split into 80 spacefilling radial tetrahedra (which we call buckets), each covering approximately the same solid angle (see Fig. LABEL:fig:bucket_map_r_ub). Figures 14 and 15 show that there are indeed significant buckettobucket fluctuations when the radius of the upper boundary is defined by the steepest decline in .
From the spatial point of view, the left panel of Fig. LABEL:fig:bucket_map_r_ub shows that the lowestluminosity run d1 contains a distinct and persistent^{4}^{4}4Animated versions of the plots from Fig. LABEL:fig:bucket_map_r_ub show this effect clearly, see the movies at the URLs given in Section 5 pattern in the form of two largescale circular concentrations of mostly negative values and a lack of negative values on the rest of the sphere. This effect results from the grid imprint discussed in Sect. 3.1, but even in the lowestluminosity d1 and d2 runs it is weak enough not to influence mass entrainment rates as shown in Sect. 3.4. The flow in the higherluminosity runs is stronger and it quickly destroys the grid imprint as illustrated in the right panel of Fig. LABEL:fig:bucket_map_r_ub on the example of the d8 run, which has the second lowest luminosity considered ( the luminosity of d1, see Table 1).
The temporal evolution of the upper boundary is shown in Fig. 16. One can see that the boundary trembles stochastically and has a thickness, as defined by the spatial fluctuations shown in Fig. 16, of Mm () and Mm in the d1 and d8 runs, respectively. Bucket data are not available for the rest of our runs, but volume renderings^{5}^{5}5See movies at the URLs given in Section 5 suggest that the fluctuations at higher luminosities are larger still as one would expect. As the material from above the convection zone is entrained, the mean radius of the boundary statistically increases. This effect is difficult to see in the left panel of Fig. 16, because the average velocity of the boundary is small in the d1 run, but the mean radius of the boundary is clearly seen to increase with time in the right panel of Fig. 16, which shows the d8 run with the luminosity of d1. Measured values of are given in Table 1. The temporal fluctuations in the mean boundary radius around the linear trend are about Mm in the d1 and d2 runs after the initial transient has passed. When extrapolated to , the boundary radius in the d1 (d2) run is Mm ( Mm). These radii coincide within Mm (; less than of the spatial fluctuations in , see above) with the radius Mm of the formal Schwarzschild boundary in the PPMstar setup.
3.4 Mass entrainment at the upper convective boundary
The entrained mass is defined in the present work to be the total amount of the initially stable fluid (see Sect. 2.2) present in the convection zone at time ,
(2) 
where is the density of fluid alone. The lower integration limit is in practice taken to be (i.e. the bottom of the simulation box), because the amount of fluid getting below the convection zone is utterly negligible^{6}^{6}6This also holds true for the highluminosity runs, including d10. (see Fig. 19). The upper integration limit would ideally be taken equal to the radius of the upper boundary as defined in Sect. 3.3 using the gradient of the sphericallyaveraged RMS tangential velocity . However, experience has shown that the fractional volume of fluid is so high at as compared to the rest of the convection zone that the resulting becomes very sensitive to small fluctuations in . Therefore we use , where the scale height of is evaluated at .
The profiles of for runs d1 and d2 are shown in Fig. 17. The initial stage of flow development, which is influenced by the initial transient and grid imprints, is not considered. The entrainment rate is calculated by performing a least squares fit to the time series in , which was constructed by applying the formula in Eq. 2 at each available output dump (every s of simulated time). The entrainment rate obtained at a resolution of (run d2) is , which is only 16 per cent higher than the entrainment rate of obtained at (run d1).
The remaining runs listed in Table 1 were analysed in the same way. The profiles of and the associated bestfit entrainment rates for these runs can be found in Fig. 22 in Appendix A. The entrainment in the two highestluminosity runs (d9, d10) was found to be so fast that also a certain period of time at the end of the simulation had to be removed from the entrainment analysis, because the flow started to become influenced by the presence of the upper boundary of the simulation box well before the end of the simulation.
The only parameter changed in the series of simulations d1, d5, d6, d8, d9, and d10 is the luminosity of the heat source, which spans almost two orders of magnitude. The left panel of Fig. 18 shows that the entrainment rate is almost directly proportional to the luminosity . This scaling predicts an entrainment rate of for the luminosity of the actual stellar model considered ().
Likewise, there is a tight relation between the entrainment rate and the shear velocity at the upper boundary as shown in the right panel of Fig. 18. The shear velocity is approximated by the value of the RMS tangential velocity at the point where reaches a local maximum just below the upper boundary of the convection zone (see Fig. 11).
In Sect. 3.3, a potential concern was raised about the presence of a persistent grid imprint in the d1 and d2 runs. Figure 18 shows that these two runs do not deviate in any significant way from the scaling relations established by the higherluminosity runs. The grid imprint is apparently weak enough not to influence the entrainment rates even at the low luminosity of the d1 and d2 runs.
3.5 Modelling convection as a diffusive process in spherically symmetric models
One of the important goals of simulating stellar convection in 3D hydrodynamic codes such as PPMstar is to inform the treatments of convection and convective boundary mixing in 1D stellar models. In such sphericallysymmetric 1D models, convection is most commonly approximated as a diffusive process. Hence, the approach taken here to reconnect the mixing in the 3D simulation with 1D simulations is to try to answer the question: what diffusion coefficient profile would have been needed in 1D in order to reproduce the sphericallyaveraged mixing in the 3D simulation?
The diffusion equation was discretised to give
(3) 
where
is the mass fraction of a particular component of the fluid, is the diffusion coefficient, is the radius, and are spatial and temporal indices, respectively, and is the time between temporal indices and .
After having performed the 3D simulation in PPMstar, the mass fraction is known from spherical averages of the results, and of course the time between spherically averaged profiles is also a known quantity. Applying the boundary condition allows equation 3 to be solved for and, hence, provide an answer to the question posed above.
3.5.1 General shape of the diffusion coefficient
The method described above has been used to calculate the 1D diffusion coefficient profile across the whole depth of the Oburning convection zone that represents the spherically and temporally averaged mixing of the 3D hydrodynamic PPMstar d2 simulation (see Table 1) between s and s. The radial abundance profiles of the upper fluid (as are shown in Fig. 19) used for the initial and final conditions of the diffusion problem are 400 s (on the order of a convective turnover time scale) time averages centred at and . This averaging smooths out some of the noise that is inherent when stochastic fluctuations in the direction of inclination or azimuth appear on dynamical timescales. The inverted discretised diffusion equation is then solved for over the time step of . The resulting diffusion coefficient is shown as the brown line in Fig. 20. When the gradient of the abundance approaches 0 (essentially flat; see Fig. 19), there is not a unique solution to the inverse diffusion problem, only a minimum diffusion coefficient that would completely mix two neighboring zones when left for a time , above which any value would be a valid solution. This is the origin of the noise seen in Fig. 20: a loss of sensitivity of the method as approaches 0. In the noisy region, a piecewiselinear representation of a downsampled diffusion coefficient was constructed, to which a cubic function was fit, resulting in the solid black line in Fig. 20. This represents a kind of byeye fit to the lower envelope of the diffusion coefficient profile in the very noisy region between radial coordinates of about 5 and 6 Mm.
There are a few other interesting quantities plotted in Fig. 20 in addition to the diffusion coefficient derived as described above. Firstly, the diffusion coefficient computed using MLT in a representative MESA simulation is shown, which compares within an order of magnitude, underestimating in the middle of the convection zone (where the method used to derive loses sensitivity) and overestimating near the convective boundary. The light blue dashed line marked with circular glyphs is a diffusion coefficient given by the same formula as the one used in MLT () but is the average radial velocity of the convective fluid from the 3D PPMstar simulation. It is within a factor of of the value predicted using MLT. When one considers that the luminosity of the PPMstar simulation is generally within a factor of the MESA model (see Fig. 5) together with the scaling of velocity with luminosity (Fig. 13), the PPMstar and MLT velocities agree within a factor of inside the convection zone. Interestingly, Herwig et al. (2006) also found in their 2D simulations that the MLT velocities were too low by a factor of about 2–3 (their Fig. 24). Such a difference between 3D hydro and 1D MLT velocities can be due to neglect of kinetic energy flux in MLT or the different geometric and averaging approximations made in various flavours of MLT (see, e.g., Tassoul, Fontaine & Winget, 1990; Salaris & Cassisi, 2008).
In any case, using overestimates in the upper and lower regions of the convective layer, too. The blue starred line looks somewhat different. This line assumes instead that the mixing length is the minimum of and the distance to the upper convective boundary. This reproduces better the shape of the derived D as it falls off in the upper half of the convection zone. The reduction of such a mixing length as the flow approaches the upper convective boundary was also suggested by Eggleton (1972) and is analogous to the decrease of the horizontal scale of convection cells over a similar domain, as has been shown by Porter & Woodward (2000, their Fig. 5).
3.5.2 Convective boundary mixing
The results of the determination of a radial diffusion coefficient profile are shown in Fig. 21 for the upper convective boundary of the d1 and d2 simulations. The vertical dotted lines mark the upper boundary of the convection zone between the convectively unstable (on the left of the boundary) and stable (on the right) fluids. The two lines represent different definitions of the convective boundary. B is the radius at which the entropy gradient becomes positive in the initial stratification of the 3D simulations and C is where the radial gradient of the tangential component of the fluid velocity is steepest in the 3D simulation (see Section 3.3) after 40 (10) minutes of simulated time for the d1 (d2) run. This is where the convective boundary is in 3D. The underlying MESA model has been aligned so that its Schwarzschild boundary lies exactly on top of B, so that a better comparison can be made between the 3D and 1D simulation results.
The abundance of the upper (initially stable) fluid is shown as a function of radius at two different times^{7}^{7}7The actual subscript of in Fig. 21 denotes the dump number of the profile quantities. In this simulation, the dump interval was 10 seconds of simulated time. by (brown line corresponding to time ) and (grey line corresponding to time ). The abundance of the upper fluid changes over the domain of the figure as the fluid from the overlying stable layer is entrained into the convection zone. The solid light blue curve shows the diffusion coefficient (with ) used in the underlying MESA model using mixing length theory to estimate the convective velocities. The diffusion coefficient resulting from our calculation is shown as a solid black curve. The result was verified by solving the diffusion equation over the d2 problem domain using as the initial abundance profile, as the time step and our spatiallyvarying diffusion coefficient. The abundance profile at was obtained exactly. The method used here can only distinguish mixing in which some of the upper and lower fluids have been exchanged. Because of this, we are not able to constrain the mixing above roughly cm (for d1; cm for d2) using the present methods with the results from the simulations presented in this work. Combining the exponentially decaying mixing model of Freytag, Ludwig & Steffen (1996) and the mixing length modification described in Section 3.5.1, the shape of the diffusion coefficient can be very well reproduced in a 1D model. This is plotted as and will be explained in greater detail in Section 3.6. The difference between the two panels in Fig. 21 is the numerical resolution of the PPMstar simulations: in the top panel and in the bottom panel. The exponential decay of across the convective boundary in both panels is best fit with an efolding length of (i.e., in Eq. 5). It is a very encouraging result that the shape and scale height of the derived diffusion coefficient changes so little upon refinement of the computational grid. In Appendix B, the same plots are shown with different values of , illustrating the suitability of to this particular convective boundary.
3.6 A mixing prescription for 1D codes
In Fig. 21 a recommended diffusion coefficient profile is drawn that, based on our PPMstar simulations, is a better representation of the space and time averaged mixing during O shell burning than using MLT with the Schwarzschild or Ledoux criteria. It is a combination of two simple considerations that one can easily adopt in a 1D stellar evolution code. Firstly, instead of using simply for the mixing length, where is a constant of order unity, using the minimum of and the distance to the upper convective boundary as the mixing length, as suggested by Eggleton (1972), much better reproduces the falloff of the diffusion coefficient inside the convection zone approaching the convective boundary (see Fig. 20 and Section 3.5.1). The recommended diffusion coefficient is then written as
(4) 
where is the radius of the Schwarzschild convective boundary. contains a constant factor of 3, which has cancelled the factor of one third that usualy appears in Eq. 4. This factor accounts for the remaining discrepancy between the RMS radial velocities from the 3D hydrodynamic simulation and those predicted using MLT (see Section 3.5.1) after the velocities have been rescaled by the relation shown in Fig. 13 for better agreement between the driving luminosities in the PPMstar and MESA simulations. Secondly, applying the exponentially decaying diffusive convective boundary mixing model of Freytag, Ludwig & Steffen (1996) in combination with the modified mixing length seems to reproduce the shape of the derived diffusion coefficient rather well. In Fig. 21 this has been plotted, and is formulated for the upper convective boundary as
(5) 
where is the radial coordinate (i.e. inside the convection zone, with the scale height of pressure at the Schwarzschild boundary. Lastly, is the parameter that sets the efolding length of the exponentiallydecaying diffusion coefficient—which is in units of the presure scale height —outward from the radius . In Fig. 21 a value of was used. Note that this mixing model is already implemented in the MESA stellar evolution code, and the recommendation made in this section would be implemented by setting the parameters and with an additional minor modification to multiply the diffusion coefficient by the factor 3.
3.7 Temperature fluctuations and their feedback on the energy generation rate
Our use of a static heat source instead of a nuclear network raises the question of what effect the convectioninduced temperature fluctuations could have on the energy generation rate in a real star. The temperature sensitivity of reactions,
(6) 
with the reaction rate given by Eq. 18.75 of Kippenhahn, Weigert & Weiss (2012) is
(7) 
where ; the temperature sensitivity of the electron screening factor has been neglected. Eq. 7 gives for , which is the typical temperature at the bottom of the O shell convection zone. The relative RMS temperature fluctuations reach in the heating layer in the d1 run. The expected relative RMS fluctuations in the energy generation rate are then . With such a small value of , it seems safe to neglect the fluctuations in the energy generation rate even after having allowed for some spread around the RMS values.
4 Comparison to other works
As shown in Fig. 13, the luminosityvelocity scaling established by our set of simulations agrees within dex with the maximum radial velocity reached in the 3D simulation of O shell convection by Meakin & Arnett (2007, MA07) (see their Fig. 6). This agreement is surprisingly good given the significant differences in physics and geometry between the simulations. However, a more detailed comparison of our Fig. 11 with Fig. 6 of MA07 reveals that in the work of MA07:

the decrease in the flow speed towards the top of the convection zone is much more pronounced and

the local maxima in the tangential velocity at the convective boundaries are narrower.
The velocity field in the simulation of MA07 may be influenced by their constraining the flow to a narrow ( deg) wedge whereas Figs. 8 and 10 show that the flow develops largescale structures if it is allowed to. Another contribution to the decrease of the flow speed with radius could be provided by neutrino cooling, which is neglected in this work, but included in the work of MA07. The velocity field at the bottom of the convection zone is likely influenced by the structure of the heating source, which drives the flow. MA07 use a nuclear reaction network comprising 25 nuclei to compute the heating rate. Figure 3 shows that the artificial heat source used in our work is more extended and centred at a slightly larger radius compared to the energy source that would be provided by the actual nuclear reactions.^{8}^{8}8Although the detailed nuclear network used in the MESA simulations shown in Fig. 3 likely differs from that used by MA07, one would expect qualitatively the same behaviour. Preliminary results from a PPMstar run with a simple nuclear network indeed show a flatter profile with narrower maxima in at the convective boundaries.
The flow in the 3D simulations of a similar O shell by Kuhlen, Woosley & Glatzmaier (2003) reaches an RMS velocity of km s. Based on the velocity scaling shown in Fig. 13 and the RMS velocity of km s observed in the d2 run (see Sect. 3.2), one would expect and RMS velocity of km s at the luminosity of used by Kuhlen, Woosley & Glatzmaier (2003). These authors simulate the whole sphere, although it is unclear how the heating is modelled. They obtain a rather flat velocity profile. Unfortunately, the information they give about the velocity components does not allow a detailed comparison.
The left panel of Fig. 18 illustrates that the 3D simulation of MA07 follows the luminosityentrainmentrate scaling established by our simulations surprisingly well. However, the right panel of the same figure shows that the entrainment rate of s is measured by MA07 at a tangential velocity times lower compared to our velocityentrainmentrate scaling. The stability of the upper convective boundary, as measured by the profile of the squared buoyancy frequency , in the simulation of MA07 is very similar to our setup (compare their Fig. 2 with our Fig. 6). However, as discussed in Sect. 3.1, entrainment in our simulations is mostly observed at places where two neighbouring largescale convective cells meet and draw the partiallymixed material from the convective boundary to the interior of the convection zone. There are only a few such places in our simulations, but there could be many more (as extrapolated to the full sphere) in the simulation of MA07, in which the artificial lateral boundaries of their simulation box constrain the largest spatial scales in the velocity field (discussed above). This effect could compensate the velocity disadvantage mentioned above.
Spruit (2015) argues that the entrainment rate is limited by the buoyancy arising from the relatively lower mean molecular weight of the entrained material. His simple model predicts a linear dependence of the entrainment rate on the convective luminosity, which is close to the relation measured in this work (see the left panel of Fig. 18). In the case of the entrainment of a fluid with a mean molecular weight into a convective fluid with a mean molecular weight , Eq. 10 of Spruit (2015) becomes:
(8) 
where is the mass fraction of the lighter fluid, and the mean molecular weight of the mixture. The equivalent of Spruit’s Eq. 12 is then:
(9) 
where the subscript refers to “ingestion” (which is called “entrainment” in this work), is a dimensionless efficiency parameter, and erg g K the gas constant. In the O shell simulations presented in this work, and , so Eq. 9 becomes
(10) 
which contains a numerical factor times larger than that in Eq. 12 of Spruit (2015), because the two fluids in the O shell problem are much more similar in terms of than the helium and heliumcarbon fluids in the coreheliumburning problem investigated by Spruit (2015). The temperature at the top of the Oshell convection zone is K. At the luminosity of erg s (runs d1 and d2), Eq. 10 gives s and the entrainment rate measured in d2 corresponds to the efficiency parameter . For comparison, Spruit (2015) estimates for the ingestion of helium into a heliumburning core.
5 Summary and conclusions
In this work, 3D hydrodynamic simulations representing the first Oburning shell following core O extinction in a star were performed using the PPMstar code. The initial setup of the stratification was constructed using three piecewise polytropes calibrated to match the O shellburning structure of a MESA model with initially solar metallicity (). The luminosity from O burning was approximated using a constant volume heating rate that also accounted for energy losses due to thermal neutrino production and emission. The simulations were performed on and computational grids, including a heating series at in which several simulations were performed, each with a heating rate scaled from the original value by a constant factor.
At the representative heating rate of , the run d2 entrained material into the convection zone at the top boundary at a rate of s, while the run d1 did so at a rate of s. This is an increase of only 16 per cent with twice the numerical grid resolution in each spatial dimension. The upper convective boundary is best defined at the radius with the steepest radial gradient in the tangential velocity and the fluctuations in the radius of the convective boundary over the full domain are about cm. Larger fluctuations are present for higher driving luminosities. The entrainment rate at the upper convective boundary scales almost linearly with the driving luminosity of the convection and approximately with the cube of the peak RMS shear velocity at the upper convective boundary. The heating series also reproduces the expected scaling of velocities with the cube root of the driving luminosity.
In order to inform approximate mixing models for 1D stellar models, the 3D simulation data was folded down to give the average values of quantities on a sphere. Averaging the resulting radial profiles over approximately one convective turnover timescale then gives spherically symmetric profiles that may be comparable to those obtained from 1D stellar evolution codes that take time steps much longer than a dynamical time. Diffusion coefficient profiles that accurately represent the mixing of the radiallyaveraged 3D simulation data are comparable with those computed using MLT. However, MLT overestimates the diffusion coefficients near the convective boundary and likely underestimates them inside the convection zone. This can be corrected by limiting the mixing length to the distance to the convective boundary, as has been suggested by Eggleton (1972), and multiplying the diffusion coefficnent by a factor of 3. Combining this simple modification with the exponentiallydecaying diffusion model of Freytag, Ludwig & Steffen (1996) at and across the convective boundary gives a mixing model for O shell burning in 1D codes that reproduces the sphericallyaveraged abundance profile evolution. For the case of O shell burning the parameter in Eq. 5 at the upper convective boundary that best reproduces the sphericallyaveraged mixing in the 3D hydrodynamic simulation is .
High resolution simulations of stellar convection with an idealised approach to microphysics appear to be a very useful technique with which to study stellar convection in the O shellburning episodes of massive stars. As others have found, a key challenge is the interpretation of 3D simulation data in the context of spherically symmetric 1D stellar models. This must of course involve both the geometric averaging of quantities that fluctuate over spheres and the temporal averaging of quantities over long enough timescales to be meaningful in the context of the assumptions that are made for hydrostatic 1D models. The analysis presented in Section 3.5 is a somewhat pragmatic step in a useful direction, but it is not intended to be the predictive model that is so highly soughtafter in the stellar physics community.
Movies of the simulations performed in this work are available at http://www.lcse.umn.edu and http://csa.phys.uvic.ca/research/stellarhydrodynamics/movies/oshellconvectionmovies
Acknowledgments
SJ is a fellow of the Alexander von Humboldt Foundation and acknowledges support from the Klaus Tschira Stiftung (KTS). RA, a CITA national fellow, acknowledges support from the Canadian Institute for Theoretical Astrophysics. This work has been supported by NSF grant PHY1430152 (JINA Center for the Evolution of the Elements). The 3D hydrodynamical simulations reported here were carried out in part on the Compute Canada WestGrid Orcinus cluster at UBC, Canada and in part on the Blue Waters machine at NCSA under the support of NSF PRAC awards 0832618 and 1440025 and CDS&E award 1413548 to the University of Minnesota. FH acknowledges support from a NSERC Discovery grant. The development of the PPM code used in this work has also been supported by contracts with the Los Alamos National Laboratory and the Sandia National Laboratory through the University of Minnesota. FH would like to thank Roman Baranowski from UBC/WestGrid for his kind support and SJ and FH thank Michael Bennett for his preliminary work on the diffusion analysis presented in this paper. We would like to thank Bernhard Müller for several helpful comments.
References
 Abbott et al. (2016a) Abbott B. P. et al., 2016a, ApJL, 818, L22
 Abbott et al. (2016b) —, 2016b, Physical Review Letters, 116, 061102
 Arnett (1994) Arnett D., 1994, ApJ, 427, 932
 Arnett & Meakin (2011) Arnett W. D., Meakin C., 2011, ApJ, 733, 78
 Asida & Arnett (2000) Asida S. M., Arnett D., 2000, ApJ, 545, 435
 Bazán & Arnett (1994) Bazán G., Arnett D., 1994, ApJL, 433, L41
 Bazán & Arnett (1998) —, 1998, ApJ, 496, 316
 Belczynski et al. (2016) Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016, Nat., 534, 512
 Biermann (1932) Biermann L., 1932, Zeitschrift fuer Astrophysik, 5, 117
 Bossini et al. (2015) Bossini D. et al., 2015, MNRAS, 453, 2290
 Brott et al. (2011) Brott I. et al., 2011, A&A, 530, A115
 Constantino et al. (2015) Constantino T., Campbell S. W., ChristensenDalsgaard J., Lattanzio J. C., Stello D., 2015, MNRAS, 452, 123
 Couch et al. (2015) Couch S. M., Chatzopoulos E., Arnett W. D., Timmes F. X., 2015, ApJL, 808, L21
 Couch & Ott (2013) Couch S. M., Ott C. D., 2013, ApJL, 778, L7
 Eggenberger et al. (2008) Eggenberger P., Meynet G., Maeder A., Hirschi R., Charbonnel C., Talon S., Ekström S., 2008, Astrophysics and Space Science, 316, 43
 Eggleton (1972) Eggleton P. P., 1972, Monthly Notices of the Royal Astronomical Society, 156, 361
 Ekström et al. (2012) Ekström S. et al., 2012, A&A, 537, A146
 Eldridge & Stanway (2016) Eldridge J. J., Stanway E. R., 2016, ArXiv eprints
 Eldridge & Tout (2004) Eldridge J. J., Tout C. A., 2004, MNRAS, 353, 87
 Ertl et al. (2016) Ertl T., Janka H.T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124
 Frey, Fryer & Young (2013) Frey L. H., Fryer C. L., Young P. A., 2013, ApJL, 773, L7
 Freytag, Ludwig & Steffen (1996) Freytag B., Ludwig H. G., Steffen M., 1996, Astronomy and Astrophysics, 313, 497
 Grevesse & Noels (1993) Grevesse N., Noels A., 1993, in Origin and evolution of the elements, Pratzo N., VangioniFlam E., Casse M., eds., p. 15
 Herwig et al. (1997) Herwig F., Blöcker T., Schönberner D., El Eid M. F., 1997, A&A, 324, L81
 Herwig et al. (2006) Herwig F., Freytag B., Hueckstaedt R. M., Timmes F. X., 2006, ApJ, 642, 1057
 Herwig et al. (2011) Herwig F., Pignatari M., Woodward P. R., Porter D. H., Rockefeller G., Fryer C. L., Bennett M., Hirschi R., 2011, The Astrophysical Journal, 727, 89
 Herwig et al. (2014) Herwig F., Woodward P. R., Lin P.H., Knox M., Fryer C., 2014, ApJL, 792, L3
 Itoh et al. (1996) Itoh N., Hayashi H., Nishikawa A., Kohyama Y., 1996, ApJS, 102, 411
 Jones et al. (2015) Jones S., Hirschi R., Pignatari M., Heger A., Georgy C., Nishimura N., Fryer C., Herwig F., 2015, MNRAS, 447, 3115
 Kippenhahn, Weigert & Weiss (2012) Kippenhahn R., Weigert A., Weiss A., 2012, Stellar Structure and Evolution
 Kuhlen, Woosley & Glatzmaier (2003) Kuhlen M., Woosley W. E., Glatzmaier G. A., 2003, in Astronomical Society of the Pacific Conference Series, Vol. 293, 3D Stellar Evolution, Turcotte S., Keller S. C., Cavallo R. M., eds., p. 147
 Levesque et al. (2006) Levesque E. M., Massey P., Olsen K. A. G., Plez B., Meynet G., Maeder A., 2006, ApJ, 645, 1102
 Levesque et al. (2009) Levesque E. M., Massey P., Plez B., Olsen K. A. G., 2009, AJ, 137, 4744
 Limongi & Chieffi (2012) Limongi M., Chieffi A., 2012, ApJS, 199, 38
 Lovegrove & Woosley (2013) Lovegrove E., Woosley S. E., 2013, ApJ, 769, 109
 Martins & Palacios (2013) Martins F., Palacios A., 2013, A&A, 560, A16
 Meakin & Arnett (2006) Meakin C. A., Arnett D., 2006, ApJL, 637, L53
 Meakin & Arnett (2007) —, 2007, ApJ, 665, 690
 Meakin & Arnett (2007) Meakin C. A., Arnett D., 2007, The Astrophysical Journal, 667, 448
 Müller (2016) Müller B., 2016, ArXiv eprints
 Müller & Janka (2015) Müller B., Janka H.T., 2015, MNRAS, 448, 2141
 O’Connor & Ott (2011) O’Connor E., Ott C. D., 2011, ApJ, 730, 70
 Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
 Paxton et al. (2013) Paxton B. et al., 2013, The Astrophysical Journal Supplement Series, 208, 4
 Paxton et al. (2015) Paxton B. et al., 2015, ApJS, 220, 15
 Porter & Woodward (2000) Porter D. H., Woodward P. R., 2000, The Astrophysical Journal Supplement Series, 127, 159
 Rauscher et al. (2002) Rauscher T., Heger A., Hoffman R. D., Woosley S. E., 2002, ApJ, 576, 323
 Salaris & Cassisi (2008) Salaris M., Cassisi S., 2008, A&A, 487, 1075
 Smartt (2015) Smartt S. J., 2015, PASA, 32, e016
 Spruit (2015) Spruit H. C., 2015, A&A, 582, L2
 Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H.T., 2016, ApJ, 821, 38
 Sukhbold & Woosley (2014) Sukhbold T., Woosley S. E., 2014, ApJ, 783, 10
 Tassoul, Fontaine & Winget (1990) Tassoul M., Fontaine G., Winget D. E., 1990, ApJS, 72, 335
 Ugliano et al. (2012) Ugliano M., Janka H.T., Marek A., Arcones A., 2012, ApJ, 757, 69
 Viallet et al. (2013) Viallet M., Meakin C., Arnett D., Mocák M., 2013, The Astrophysical Journal, 769, 1
 Woodward, Herwig & Lin (2015) Woodward P. R., Herwig F., Lin P.H., 2015, ApJ, 798, 49
 Woosley & Heger (2007) Woosley S. E., Heger A., 2007, Physics Reports, 442, 269
 Young et al. (2005) Young P. A., Meakin C., Arnett D., Fryer C. L., 2005, ApJL, 629, L101
Appendix A Mass entrainment in the runs d5, d6, d8, d9, and d10
Table 1 gives entrainment rates for several simulations that have been performed as part of this work. In Fig. 17 the entrained mass is plotted as a function of time for the d1 and d2 simulations with a linear fit whose gradient is the entrainment rate. Fig. 22 shows the same kind of plot for the other runs that are listed in Table 1.
Appendix B Determination of convective boundary mixing scale height
In Section 3.5.2 it is described how a diffusion coefficient can be derived that represents the 1D sphericallyaveraged mixing properties of the 3D PPMstar hydrodynamic simulations. It is shown in Fig. 21 that combining a restricted mixing length with an exponential decay of the diffusion coefficient (, solid blue starred line) with CBM parameter reproduces the sphericallyaveraged mixing properties very well indeed. The value is not the result of a fitting procedure, however to the eye it is the best fit to one decimal place. This is illustrated in Fig. 23, which shows the same plot as Fig. 21 for the d1 simulation with different values of .