Idealised 4\pi simulations of O-shell convection

Idealised hydrodynamic simulations of turbulent oxygen-burning shell convection in geometry

S. Jones, R. Andrassy, S. Sandalski, A. Davis, P. Woodward and F. Herwig
Heidelberg Institute for Theoretical Studies, Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
Department of Physics & Astronomy, University of Victoria, P.O. Bos 3055 Victoria, B.C., V8W 3P6, Canada
Joint Institute for Nuclear Astrophysics, Center for the Evolution of the Elements, Michigan State University,
640 South Shaw Lane, East Lansing, MI 48824, USA
LCSE and Department of Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
Accepted 2016 October 27. Received 2016 October 25; in original form 2016 May 9

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 O-burning 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 spherically-averaged 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 spherically-averaged 3D abundance profiles.

stars: massive, evolution, interior — physical data and processes: turbulence, hydrodynamics, convection

1 Introduction

Modelling the long-term 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 (C-burning onwards), the structure of the star is generally said to be frozen-in. 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 core-collapse supernova of type II111however, note that note that a subset of type I supernovae (Ib and Ic) are also core-collapse supernovae, while the latter case (black hole formation) is predicted to produce a long-duration ( 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 core-collapse 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 hole-forming) or successful (neutron star-forming) 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 pre-supernova stellar structure, and different predictions for pre-SN 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.

Figure 1: Isotopes included in the nuclear reaction network in the 25  MESA model during oxygen shell burning (circles). Stable isotopes are outlined with black squares.

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 C-burning core and hence the location of the first subsequent C-burning shell. This ratio is set by the competition between the triple-alpha (3) and the C(O reactions at the end of the core He-burning 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 He-burning core, and the temperature is set by the mass of the H-depleted core, which is itself determined by the treatment of convective boundary mixing (CBM) during the main sequence (core H-burning). 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 He-burning 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 O-burning shell in a 25 star. What can be done, however, are multidimensional hydrodynamic simulations of the O-burning 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 He-burning 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 pre-supernova structure of core-collapse 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 in-falling core material with the neutrino-heated 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 shell-burning episode, the imprint of the convective velocity field is essentially frozen in. Fluctuations in the velocity amplitudes can impact the success of the simulated core-collapse supernova explosion of the progenitor star (Couch & Ott, 2013; Müller & Janka, 2015; Couch et al., 2015).

Figure 2: Kippenhahn (convective structure evolution) diagram of the core and shell oxygen-burning phases in the 25  MESA model with positive nuclear energy generation contours. In this figure, convective regions are shaded in grey, is the specific energy loss rate due to neutrino production by nuclear reactions and is the time until core collapse. The turquoise dot-dashed line marks the boundary of the C-free core, which is defined as the mass coordinate below which the mass fraction of C is lower than . The initial setup of our PPMstar simulations is based on the structure of the convective O shell in this MESA model at , i.e. about 10 days before core collapse (marked with a vertical dashed line). The solid portion of the vertical line at shows the region of the star that was simulated in PPMstar.

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 O-burning shell and the convective C-burning shell merged into one deep convective layer, engulfing in addition the Ne-burning shell. In their computations, this event was responsible for a late-time neutron burst that resulted in large overproduction of neutron-capture 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 O-burning 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.

Figure 3: Nuclear energy generation rate (), thermal neutrino energy loss rate () and luminosity profiles in the 25  MESA model during shell oxygen burning about 10 days before core collapse ( in Figs. 2 and 5). The oxygen shell reaches from about 1.15  () to about 1.85  (). The second peak at about 2.08  () is the base of the convective C shell. The extent of the PPMstar hydrodynamic simulations, in which the combined effects of nuclear burning and neutrino losses is approximated by the solid dark grey curve , can be seen in Fig. 4. The convectively unstable regions in the underlying MESA model are shaded (yellow: O-burning shell; purple: C-burning shell).

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 He-burning shell flash during a very late thermal pulse (convective He shell flash; VLTP) of a post-AGB star.

In this work, idealised simulations of convection in the first O-burning 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 He-shell 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 O-shell 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.

Figure 4: Contributions to the total pressure (top panel) and density stratification (bottom panel) during the shell O-burning phase of the representative MESA model about 10 days before core collapse ( in Figs. 2 and 5). The initial pressure and density stratifications of the PPMstar hydrodynamic simulations are shown with subscript PPM (solid dotted grey lines). The convectively unstable regions in the underlying MESA model are shaded (yellow: O-burning shell; purple: C-burning shell).

2 Simulations

Two types of simulations of the first convective O-burning 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 high-resolution grid and up to 55 minutes using our low-resolution 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 non-dimensionalised by putting  in the denominator, such a definition would be problematic for convective cores for which .

Figure 5: Time evolution of the defining quantities of the convective O-burning shell in a stellar model computed with the MESA stellar evolution code. The x-axis is the log of the time remaining until core collapse (cf. Fig. 2);  and  are the radius of the bottom and top of the convective O-burning shell, is the mean molecular weight in the convective region and is the peak luminosity driving the convection.

2.1 MESA simulations

Stellar evolution calculations of a non-rotating massive star with metallicity were performed using the MESA code222The 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 helium-burning. 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 oxygen-burning shell was approximated with an exponentially-decaying 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 e-folding 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 inner-most ) during the core and shell O-burning 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 O-burning 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 oxygen-burning shell is a deep convective C-burning 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 O-burning 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 O-burning 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 O-burning 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 O-shell, but in the convective C-burning 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 post-He 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 O-burning 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).

Figure 6: Radial profiles of the squared Brunt-Väisälä frequency at the beginning of the simulation and at the 22 minute in the d2 run. The MESA model is shown for comparison and to illustrate that the amplitude and slope of the profile at the upper boundary in the PPMstar model are very similar.

2.2 PPMstar simulations

Hydrodynamic simulations of O-shell convection were performed in 3D with geometry using the PPMstar code. The code as well as its performance for a similar application (He-shell flash convection in a low-mass 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 He-shell flash convection simulations. In summary, an ideal gas equation of state is adopted and a geometrically representative stratification is produced by combining three piece-wise 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 over-estimated 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 O-burning 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 O-burning shell to have when much more detailed microphysics have been considered.

Figure 7: Fractional volume profile for He-shell (Woodward, Herwig & Lin, 2015) at time step 808672 (AGB) and O shell (d2 from this work) runs at the same product of number of time steps and Courant number (i.e. the same computational effort). The much larger amount of material that has been entrained in the O shell simulation compared with the AGB simulation after expending the same amount of computational effort is the reason why O-burning shell is well suited for studying the convergence properties of PPMstar simulations of stellar convection.
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
Table 1: Summary of the hydrodynamic simulations that have been performed with the PPMstar code for this work, ordered by driving luminosity. d1 and d2 are both representative of the actual conditions in the first O-burning shell in a star (c.f. Fig. 5) and the other simulations constitute a heating series in which the properties of convection and entrainment are studied as a function of the luminosity driving the convection. The entrainment rate is given for the upper convective boundary. is the maximum tangential component of the velocity of the fluid in the upper half of the convection zone.
Figure 8: Fractional volumes of the entrained fluid in the 1536 PPMstar run d2 at 25.7 minutes of simulated time (dump number 154). The left image is a projection of the far hemisphere of the simulation; the near hemisphere has been cut away and the core has been made almost transparent, but is still visible as a faint purple sphere. The right image is a thin slice through the sphere. In both images, fluid that began inside the convection zone has been made transparent and therefore only entrained fluid is visible. That entrained fluid has a fractional volume that is indicated by the colour scale. See Section 3.1 of the text for a detailed description.

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 constant-volume 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 Brunt-Väisälä frequency333This 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.

Figure 9: A volume rendering of the vorticity of the fluid in the PPMstar run D2 at of simulated time (dump number 154). From the outside inward the domain boundary, a stable layer, the convectively unstable region, a rather thin stable region below the turbulent convection zone and the inert core can be distinguished.

The O-burning 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 O-burning 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 O-burning shell convection problem the same ratio is only 1.02. Another difference is that in these O-shell 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 He-shell 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.

Figure 10: Radial velocity of the entrained fluid in the 1536 PPMstar run d2. See Section 3.1 of the text for a detailed description.

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 He-shell 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.

Figure 11: Spherically-averaged total, radial and tangential fluid velocity components after 22 minutes of simulated time (dump number 132) in the d2 (thick lines) and d1 (thin lines) simulations (see Table 1).

The character of the flow above the O-burning shell in this problem is quite similar to that which we have seen in the shell-flash convection zone during the very late thermal pulse (VLTP) of post-AGB 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 O-burning shell convection zone compared to the AGB He shell convection zone after expending the same amount of computational effort makes the O-burning 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 semi-transparent, 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.

Figure 12: Temporal evolution of the maximum RMS radial velocity in the convection zone for all of the runs listed in Table 1. All of the runs experience an initial transient event after which the convection settles down into a steady state. The d1 and d2 simulations are in very good agreement following their respective initial transients.

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 Kelvin-Helmholtz 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.

Figure 13: Dependence of the maximum RMS radial velocity on the luminosity . The RMS velocity profiles were averaged over data dumps (spanning  s of simulation time) centred in the middle of the time interval, in which the mass entrainment rate was measured for the given run. The fitting relation, which assumes the same units as the axes of the plot, was constructed using the runs only. The run d2 and the 3D simulation of Meakin & Arnett (2007) are shown for comparison.

3.2 Velocity profiles

The spherically-averaged 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.

Figure 14: RMS tangential velocity profiles at the upper convective boundary in the d1 (left panel) and d8 (right panel) runs extracted at from 80 radial buckets (see Fig. LABEL:fig:bucket_map_r_ub and Section 3.3). The radius of the steepest decline in , which is used to define the position of the convective boundary, is marked for each profile separately by a short vertical line. The long, dashed and dotted vertical lines show the average value of and fluctuations in , respectively. The standard deviation was computed separately for positive and negative fluctuations. The lower part of each plot shows the distribution of values collected from 35 data dumps in the time interval from to . The width of the histogram bins corresponds to the cell size of the underlying PPMstar computational grid.
Figure 15: As Fig. 14, but is shown for the d1 run in the left panel and for the d8 run in the right panel.

The flow velocity is closely related to the luminosity of the heat source. The best-fit 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,


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 mixing-length 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.

Figure 16: Evolution of the radius of the upper boundary in the d1 (left panel) and d8 (right panel) runs as defined by the steepest gradient in the RMS tangential velocity (see Figs. 14 and 15). The full sphere is split into 80 radial buckets (see Fig. LABEL:fig:bucket_map_r_ub and Section 3.3). The standard deviation was computed separately for positive and negative fluctuations. The large-amplitude oscillations in the first minutes result from the initial transient (see Section 2.2).

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 space-filling 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 bucket-to-bucket 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 lowest-luminosity run d1 contains a distinct and persistent444Animated 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 large-scale 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 lowest-luminosity d1 and d2 runs it is weak enough not to influence mass entrainment rates as shown in Sect. 3.4. The flow in the higher-luminosity 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 renderings555See 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 ,


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 negligible666This also holds true for the high-luminosity 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 spherically-averaged 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 best-fit entrainment rates for these runs can be found in Fig. 22 in Appendix A. The entrainment in the two highest-luminosity 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 higher-luminosity 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 spherically-symmetric 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 spherically-averaged mixing in the 3D simulation?

Figure 17: Entrained mass as a function of time for the run d1 (left panel) and the run d2 (right panel). The two entrainment rates agree within 16 per cent. The dashed line in each of the panels shows the best-fit slope from the other panel for comparison. The values of have been adjusted for the dashed lines, because different amounts of mass were entrained during the initial transients in the two runs.

The diffusion equation was discretised to give



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 O-burning 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 time-scales. 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 piecewise-linear 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 by-eye 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).

Figure 18: Dependence of the entrainment rate on the luminosity (left panel) and on the shear velocity (right panel). 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). The RMS velocity profiles were averaged over data dumps (spanning 100 s of simulation time) centred in the middle of the time interval, in which the mass entrainment rate was measured for the given run. The fitting relations, which assume the same units as the axes of the plots, were constructed using the runs only. The run d2 and the 3D simulation of Meakin & Arnett (2007) are shown for comparison.

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 times777The 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 spatially-varying 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 e-folding 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.

Figure 19: Spherical averages of the fractional volumes of the upper (initially stable) fluid as a function of radial coordinate at different times over the duration of the d2 simulation.

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 fall-off 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


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 re-scaled 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


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 e-folding length of the exponentially-decaying 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 convection-induced temperature fluctuations could have on the energy generation rate in a real star. The temperature sensitivity of reactions,


with the reaction rate given by Eq. 18.75 of Kippenhahn, Weigert & Weiss (2012) is


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 luminosity-velocity 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:

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

  2. 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 large-scale 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.888Although 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.

Figure 20: Time-averaged radial diffusion coefficient profile calculated from the spherically-averaged abundance profiles by the method described in Section 3.5 (brown solid line; black solid line is a fit to the noisy region). The convective velocities computing using MLT agree with the spherically-averaged 3D velocities to within about a factor of 2 inside the convection zone but are too large in the vicinity of the convective boundary, resulting in an overestimation of the diffusion coefficient there. Limiting the mixing length to the distance from the convective boundary reproduces the fall-off of the diffusion coefficient inside the convection zone approaching the boundary that is seen in the spherically averaged 3D simulation results.

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 luminosity-entrainment-rate 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 velocity-entrainment-rate 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 set-up (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 large-scale convective cells meet and draw the partially-mixed 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:


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:


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


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 helium-carbon fluids in the core-helium-burning problem investigated by Spruit (2015). The temperature at the top of the O-shell 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 helium-burning core.

5 Summary and conclusions

In this work, 3D hydrodynamic simulations representing the first O-burning 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 shell-burning 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.

Figure 21: Results of the 3D–1D diffusion analysis at the upper convective boundary of the d1 ( grid) and d2 ( grid) simulations (see Table 1). The vertical dotted lines represent the upper boundary of the convection zone. B is where the entropy gradient becomes positive in our PPMstar setup (equivalent to the Schwarzschild criterion); C is where the radial gradient of the tangential component of the fluid velocity is steepest after 46.7 (16.0) minutes of simulated time for simulation d1 (d2). We also give the MESA model upon which these simulations were based; it has been aligned so that the convective boundary according to the Schwarzschild criterion is located at B. is the spherically-averaged mass fraction of the overlying fluid and is plotted at a simulated time indicated by the subscript in tens of seconds. is the diffusion coefficient computed in the framework of mixing length theory with . (solid black line) is the derived diffusion coefficient that gives the same net mixing as the 3D hydrodynamic simulation when its output is spherically averaged. is the recommended diffusion coefficient to use in a 1D code given by , where is the radial coordinate of the Schwarzschild boundary at B, as described in Section 3.6 of the text, with an exponentially decaying convective boundary mixing from radius with , as formulated by Freytag, Ludwig & Steffen (1996, see Section 3.6).

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 radially-averaged 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 exponentially-decaying 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 spherically-averaged abundance profile evolution. For the case of O shell burning the parameter in Eq. 5 at the upper convective boundary that best reproduces the spherically-averaged 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 shell-burning 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 sought-after in the stellar physics community.


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 PHY-1430152 (JINA Center for the Evolution of the Elements). The 3-D 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.


  • 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., Christensen-Dalsgaard 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 e-prints
  • 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., Vangioni-Flam 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 e-prints
  • 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.

Figure 22: Entrained mass as a function of time in the runs d5, d6, d8, d9, and d10. The entrained mass is fit with a straight line, whose gradient gives the entrainment rate at the upper convective boundary of the simulation. The properties of the simulations are given 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 spherically-averaged 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 spherically-averaged 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 .

Figure 23: As Fig. 21. Comparison of the diffusion coefficient profile that reproduces the mixing of spherically-averaged 3D hydrodynamic simulation data from the PPMstar d1 simulation (, black solid line) with the exponentially decaying convective boundary mixing model of Freytag, Ludwig & Steffen (1996) with different CBM parameters. While using results in slight underestimation of the diffusion coefficient, it gives the best match to the gradient of above 8 Mm.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description