Dependence of Convective Boundary Mixing on Boundary Properties and Turbulence Strength

TL;DR

This study uses 3D hydrodynamical simulations to explore how convective boundary mixing in stellar evolution depends on boundary properties and turbulence, supporting the entrainment law's applicability in modeling.

Contribution

It provides a controlled analysis of convective boundary mixing dependence on boundary properties and turbulence strength using 3D simulations and mean field analysis.

Findings

Entrainment inversely proportional to bulk Richardson number.

Vertical RMS velocity scales with driving luminosity.

Boundary positions estimated via composition profiles.

Abstract

Convective boundary mixing is one of the major uncertainties in stellar evolution. In order to study its dependence on boundary properties and turbulence strength in a controlled way, we computed a series of 3D hydrodynamical simulations of stellar convection during carbon burning with a varying boosting factor of the driving luminosity. Our 3D implicit large eddy simulations were computed with the PROMPI code. We performed a mean field analysis of the simulations within the Reynolds-averaged Navier-Stokes framework. Both the vertical RMS velocity within the convective region and the bulk Richardson number of the boundaries are found to scale with the driving luminosity as expected from theory. The positions of the convective boundaries were estimated through the composition profiles across them, and the strength of convective boundary mixing was determined by analysing the boundaries within the framework of the entrainment law. We find that the entrainment is approximately inversely proportional to the bulk Richardson number. Although the entrainment law does not encompass all the processes occurring at boundaries, our results support the use of the entrainment law to describe convective boundary mixing in 1D models, at least for the advanced phases. The next steps and challenges ahead are also discussed.

Figures

Tables

Table 1. Table 1: Summary of simulation properties. ϵ f ​ a ​ c subscript italic-ϵ 𝑓 𝑎 𝑐 \epsilon_{fac} : nuclear energy generation rate boosting factor, τ sim subscript 𝜏 sim \tau_{\rm sim} : simulated physical time, v rms subscript 𝑣 rms v_{\rm rms} : global RMS convective velocity, τ q subscript 𝜏 𝑞 \tau_{q} : quasi-steady state time (total simulation time minus initial transient time), τ c subscript 𝜏 𝑐 \tau_{c} : convective turnover time, Ri B,u subscript Ri B,u \textrm{Ri}_{\textrm{B,u}} , Ri B,l subscript Ri B,l \textrm{Ri}_{\textrm{B,l}} : bulk Richardson number for the upper and lower convective boundary regions, respectively, Ma: Mach number.

	eps1	eps33	eps100	eps333	eps1k	eps3k	eps10k	eps33k
${ϵ}_{fac}$	1	33	100	333	1000	3000	1$\times {10}^4$	3.3$\times {10}^4$
${\tau }_{sim}$	5000	5000	3600	1640	6000	3635	2000	1190
$v_{rms}$	6.80$\times {10}^5$	1.15$\times {10}^6$	1.59$\times {10}^6$	1.92$\times {10}^6$	4.56$\times {10}^6$	6.93$\times {10}^6$	1.03$\times {10}^7$	1.57$\times {10}^7$
${\tau }_q$	3908	3900	2489	636	5003	2743	1498	592
${\tau }_c$	3257	1789	1317	1025	465	316	219	153
${\tau }_q/{\tau }_c$	1.20	2.18	1.89	0.62	10.76	8.68	6.84	3.87
$R{i}_{B,u}$	1876	559	310	200	42	19	7	4
$R{i}_{B,l}$	2.93$\times {10}^4$	8089	4203	3101	435	188	101	44
Ma	0.0028	0.0048	0.0070	0.0074	0.0209	0.0321	0.0481	0.0727

Table 2. Table 2: 2 ​ σ r , l 2 subscript 𝜎 𝑟 𝑙 2\sigma_{r,l} and 2 ​ σ r , u 2 subscript 𝜎 𝑟 𝑢 2\sigma_{r,u} are twice the standard deviations from the horizontal means of the boundary surface positions (shaded envelopes in Fig. 12 ). These measurements can be associated with boundary surface height fluctuations, possibly due to plume penetration events. δ ​ r l 𝛿 subscript 𝑟 𝑙 \delta r_{l} and δ ​ r u 𝛿 subscript 𝑟 𝑢 \delta r_{u} are the approximate boundary widths determined from the composition profiles of the lower and upper convective boundaries respectively (see Fig. 16 ). These measurements of the boundary represent the size of the boundary region which is formed due to shear mixing at the boundaries. H p , l subscript 𝐻 𝑝 𝑙 H_{p,l} and H p , u subscript 𝐻 𝑝 𝑢 H_{p,u} are the average pressure scale heights across the lower and upper convective boundary regions, respectively.

	eps1k	eps3k	eps10k	eps33k
$2{\sigma }_{r,l}\left(\mathrm{cm}\right)$	3.88$\times {10}^6$	3.88$\times {10}^6$	5.96$\times {10}^6$	6.77$\times {10}^6$
$2{\sigma }_{r,u}\left(\mathrm{cm}\right)$	1.16$\times {10}^7$	2.26$\times {10}^7$	2.25$\times {10}^7$	1.04$\times {10}^8$
$\delta r_l\left(\mathrm{cm}\right)$	4.04$\times {10}^7$	5.14$\times {10}^7$	6.24$\times {10}^7$	5.88$\times {10}^7$
$\delta r_u\left(\mathrm{cm}\right)$	1.18$\times {10}^8$	1.29$\times {10}^8$	1.40$\times {10}^8$	2.02$\times {10}^8$
$\delta r_l/H_{p,l}$	0.141	0.181	0.218	0.201
$\delta r_u/H_{p,u}$	0.339	0.366	0.389	0.535

Table 3. Table 3: Truncated table of the re-zoned stellar model initial conditions used in the prompi simulations. From left to right the following variables are listed: zone number, radius, entropy, 12 C abundance, mean atomic mass, mean atomic number, density, temperature and pressure. Values are shown for zones within the bottom of the computational domain (top rows), lower convective boundary (middle rows) and the top of the computational domain (bottom rows). The full table of values is available in machine-readable form as supplementary online material.

i	r [cm]	s [erg K-1]	X${}_{{}_{}{}^{12}C}$	$\overline{A}$	$\overline{Z}$	$\rho$ [g cm-3]	T [K]	P [bar]
0001	4.161$\times {10}^8$	1.649$\times {10}^8$	6.099$\times {10}^{-3}$	18.17	9.072	1.737$\times {10}^6$	8.642$\times {10}^8$	9.594$\times {10}^{22}$
0002	4.165$\times {10}^8$	1.650$\times {10}^8$	6.099$\times {10}^{-3}$	18.17	9.072	1.733$\times {10}^6$	8.643$\times {10}^8$	9.570$\times {10}^{22}$
0003	4.169$\times {10}^8$	1.651$\times {10}^8$	6.099$\times {10}^{-3}$	18.17	9.072	1.729$\times {10}^6$	8.644$\times {10}^8$	9.546$\times {10}^{22}$
. . .
1249	8.987$\times {10}^8$	2.979$\times {10}^8$	5.056$\times {10}^{-2}$	17.63	8.807	1.870$\times {10}^5$	9.043$\times {10}^8$	9.805$\times {10}^{21}$
1250	8.990$\times {10}^8$	3.046$\times {10}^8$	6.852$\times {10}^{-2}$	17.42	8.699	1.826$\times {10}^5$	9.201$\times {10}^8$	9.868$\times {10}^{21}$
1251	8.994$\times {10}^8$	3.123$\times {10}^8$	8.877$\times {10}^{-2}$	17.18	8.578	1.773$\times {10}^5$	9.361$\times {10}^8$	9.895$\times {10}^{21}$
. . .
4998	2.348$\times {10}^9$	3.947$\times {10}^8$	3.571$\times {10}^{-1}$	14.38	7.183	8.171$\times {10}^3$	3.472$\times {10}^8$	1.717$\times {10}^{20}$
4999	2.349$\times {10}^9$	3.948$\times {10}^8$	3.571$\times {10}^{-1}$	14.38	7.183	8.163$\times {10}^3$	3.471$\times {10}^8$	1.715$\times {10}^{20}$
5000	2.349$\times {10}^9$	3.948$\times {10}^8$	3.571$\times {10}^{-1}$	14.38	7.183	8.155$\times {10}^3$	3.471$\times {10}^8$	1.713$\times {10}^{20}$

Table 4. Table 4: Convective turnover times, boundary migration time scales and the ratio of these two time scales for each boundary of the eps1k , eps3k , eps10k , eps33k models.

	eps1k	eps3k	eps10k	eps33k
${\tau }_{b,u}\left(s\right)$	7066	3377	1346	662
${\tau }_{b,l}\left(s\right)$	15843	7649	3391	975
${\tau }_c\left(s\right)$	465	316	219	153
${\tau }_{b,u}/{\tau }_c$	15.2	10.7	6.1	4.3
${\tau }_{b,l}/{\tau }_c$	34.1	24.2	15.5	6.4

Equations

$$f={\left(1+\delta t\omega f_d\right)}^{-1},$$

$$v_{\mathrm{x},\mathrm{rms}}\propto L^{1/3},$$

$$E=\frac{v_e}{v_{\mathrm{rms}}}=A{\text{Ri}}_{\mathrm{B}}^{-n},$$

$${\text{Ri}}_B\propto v_{\mathrm{rms}}^{-2}\propto L^{-2/3},$$

$${\dot{m}}_e=\frac{\partial m}{\partial r}{v}_e=4\pi r^2\rho v_{\mathrm{rms}}A{\text{Ri}}_B^{-n}.$$

$$\frac{\partial E_K}{\partial t}=-\frac{\partial s}{\partial t}=\frac{5}{4}\frac{\Delta v^3}{\Delta r}.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})$$

$$\rho \frac{\partial \mathbf{v}}{\partial t}+\rho \mathbf{v}\cdot \nabla \mathbf{v}$$

$$\rho \frac{\partial E_t}{\partial t}+\rho \mathbf{v}\cdot \nabla E_t+\nabla \cdot (p\mathbf{v})$$

$$\rho \frac{\partial X_i}{\partial t}+\rho \mathbf{v}\cdot \nabla X_i$$

$$\Delta t=C_{\max }\frac{\Delta x}{c_s},$$

$${ϵ}_{{}^{12}\mathrm{C}}\sim 4.8\times 10^{18}{Y}_{12}^2\rho {\lambda }_{12,12},$$

$${ϵ}_{\nu }=$$

$$+\frac{21.6{\rho }^2{\lambda }^2}{\rho +8.6\times 10^5}(1+2.215\times 10^{-6}{\xi }^2)e^{-4.5855\times 10^{-3}\xi }$$

$$+4.772\times 10^2{\lambda }^2\rho e^{-2.5817\times 10^{-5}{\rho }^{2/3}{\lambda }^{-1}},$$

$${ϵ}_{\ell }=\frac{v_{\mathrm{rms}}^3}{\ell }\text{and}{ϵ}_{\Delta x}=\frac{\Delta v^3}{\Delta x},\text{respectively},$$

$${\nu }_{\mathrm{eff}}=\Delta v\Delta x.$$

$${\nu }_{\mathrm{eff}}=v_{\mathrm{rms}}\ell {\left(\frac{\Delta x}{\ell }\right)}^{4/3}.$$

$${\text{Re}}_{\mathrm{eff}}={\left(\frac{\ell }{\Delta x}\right)}^{4/3}\sim N_x^{4/3},$$

$$⟨\omega ⟩=\frac{1}{\Delta A}{\int }_{\Delta A}\omega \text{d}A,$$

$$\omega =\overline{⟨\omega ⟩}+{\omega }^{\prime },$$

$$\overline{\omega }=\frac{1}{\Delta t}{\int }_{t_1}^{t_2}\omega (t)dt$$

$${\partial }_t\left(\rho E_k\right)+\nabla \cdot \left(\rho E_k\mathbf{v}\right)=-\mathbf{v}\cdot \nabla p+\rho \mathbf{v}\cdot \mathbf{g},$$

$$\overline{⟨\rho {\mathbf{D}}_{\mathbf{t}}{E}_k⟩}=$$

$${\tau }_c=\frac{2{\ell }_{cz}}{v_{\mathrm{rms}}},$$

$${\tau }_b=\frac{\delta r}{|v_e|},$$

$${\text{Ri}}_B=\frac{\Delta B\ell }{v_{\mathrm{rms}}^2},$$

$$\Delta B=\underset{r_c-\Delta r}{\overset{r_c+\Delta r}{\int }}{N}^{2}dr.$$

$$N^2=\frac{g}{\rho }\left(\frac{\partial \rho }{\partial s}\frac{ds}{dr}+\frac{\partial \rho }{\partial \bar{A}}\frac{d\bar{A}}{dr}+\frac{\partial \rho }{\partial \bar{Z}}\frac{d\bar{Z}}{dr}\right).$$

Full text

††pagerange: Dependence of Convective Boundary Mixing on Boundary Properties and Turbulence Strength–References††pubyear: 2017

1 Introduction

One-dimensional (1D) stellar evolution models are currently the only computational tool that can be used to simulate the structure and evolution of a star from the zero-age main sequence (ZAMS) up to their final stage. Such models are invaluable due to their ability to provide estimates for stellar masses (Bersten et al., 2014), nucleosynthesis yields (Pignatari et al., 2016), progenitor structures (Arnett and Meakin, 2011a) and evolutionary ages (Nieva and Przybilla, 2014). Due to their limited dimensionality, 1D stellar models must prescribe or approximate multi-dimensional phenomena. This often leads to a loss of predictability as such parameterisations usually contain free parameters which are tuned in order for the models to match observations (e.g. the mixing length theory; Böhm-Vitense 1958). In order to model these highly complex, non-linear processes such as convection, three-dimensional (3D) hydrodynamic models can be computed which solve the equations of fluid motion (e.g. the Euler equations). Such simulations can test the underlying physical assumptions in the stellar models under semi-realistic astrophysical, macroscopic, but not necessarily microscopic conditions. These simulations resolve the dynamical time-scales associated with the fluid flow, whereas stellar models resolve the secular (thermal) time-scales. Previous 3D hydrodynamic models of the carbon and oxygen burning shells in massive stars (Meakin and Arnett, 2007; Viallet et al., 2013; Cristini et al., 2017), have revealed enlightening results: the dynamic properties of the convective flow; the balance between turbulent driving and kinetic energy dissipation (viscosity); convective boundary mixing through turbulent entrainment111Turner (1973) defines entrainment as the transport of fluid across an interface between two bodies of fluid by a shear-induced turbulent flux.; propagation of g-mode waves throughout the stable regions; and some level of agreement with the entrainment law (Eq. 3). Building upon our previous work producing the first 3D hydrodynamic models of the carbon shell (Cristini et al., 2017), we extend this by studying the effects of varying the driving luminosity within the same initial set-up, provided by the same stellar evolution model of a 15 M⊙ star.

The structure of the paper is as follows. In Section 2, we briefly describe the set-up of the initial conditions for the hydrodynamical models. The results and their analyses are presented in Section 3. We conclude our findings in Section 4. The Appendices include more details about the computational tools used, various derivations and supplementary analyses using the RANS framework, with a particular emphasis on the turbulent kinetic energy.

2 Initial physical model

The computational tools used to calculate the stellar model providing the initial conditions and hydrodynamic simulations are described in Appendix A.

As in our previous study (Cristini et al., 2017, C17 hereinafter), we follow the “box-in-star” approach (e.g. Arnett et al., 2015), use a Cartesian coordinate system and a plane-parallel geometry. Details of the stellar model initial conditions can be found in Cristini et al. (2016) and in Appendix A. The computational domain represents a convective region of thickness, $t$, bounded either side by radiative regions of thickness, $t/2$. Initially, the computational domain and convective region span $5.7$ and $2.6$ pressure scale heights, respectively. By including parts of the surrounding radiative regions, in general (the eps33k model is an exception) this ensures that over several convective turnovers the convective boundaries will not interact with the vertical domain boundaries. The aspect ratio of the convective zone is 2:1 (width:height), and therefore the plane-parallel approximation is not ideal. As in C17 this choice was made to ease the difficult Courant time scale condition at the inner boundary of the grid, allowing longer run-times as well as better resolution near convective boundaries. Direct comparison with oxygen burning simulations (e.g. Meakin and Arnett, 2007), which use a spherical grid, suggest that no significant error results.

The computational domain uses reflective, stress-free boundary conditions in the vertical direction and periodic boundary conditions in the two horizontal directions. Although the material in the radiative regions is stable against convection it has oscillatory g-mode motions excited by the adjacent convection zone. Meakin and Arnett (2006) showed that such waves in their 2D model of the carbon and oxygen shell are well described by the linearised non-radial wave equation. If left to propagate freely and due to the reflective boundary conditions, standing waves will develop which will affect the hydrostatic balance. Therefore, in order to mimic the propagation of these waves out of the domain, a damping region is employed which extends radially between a radius of $0.6\times 10^{9}$ cm and the lower domain boundary at $0.42\times 10^{9}$ cm. The damping region covers the full horizontal extent of the computational domain between these radii. Within this region all velocity components are reduced by a common damping factor, $f$, resulting in damped velocities over the damping region, $\boldsymbol{v_{d}}=f\boldsymbol{v}$. The damping factor is defined as

[TABLE]

where $\delta t$ is the time step of the simulation, $\omega=0.01$ is the damping frequency and is a free parameter chosen to correspond to a small fraction of the convective turnover. $f_{d}=0.5\;(\textrm{cos}\left(\pi r/r_{0}\right)+1)$, where $r$ is the radial position in the vertical direction and $r_{0}=0.6\times 10^{9}$ cm is the edge of the damping region in the vertical direction. Using this damping function, $f_{d}=0$ at $r=r_{0}$, where the damping region starts. This ensures a smooth transition between the non-damped and damped regions.

The energy generation rate is calculated using the same prescription as described by eq. 6 in C17, with the addition of a constant boosting factor, $\epsilon_{fac}$. This factor was set to 1000 in C17. In this study, $\epsilon_{fac}$ varies between $1$ (denoted model eps1) and $3.3\times 10^{4}$ (denoted model eps33k). Except for the energy generation rate, all of the models in this study are identical in set-up to the hrez model of C17. The radial profile of the nuclear energy generation rate is plotted for all models in Fig. 1, and are normalised to the maximum energy generation rate of the nominal luminosity model. The energy generation rate due to neutrino losses is unchanged for each model, as over such short (dynamical) time-scales it is not expected that changes in the energy generation rate due to neutrino losses will have any considerable effect on the structure of the shell.

In order to be computationally efficient while maximising the resolution, a mesh size of $512^{3}$ was chosen for all models. This resolution was shown to be sufficient to model the upper convective boundary for the hrez model of C17. Their vhrez model ($1024^{3}$) at $\epsilon_{fac}=10^{3}$ modelled the lower convective boundary more accurately, but it is currently too computationally expensive to compute many models at such a resolution. With a resolution of 5123, the effective Reynolds number of these simulations is Re${}_{\rm eff}\sim$ 4000, placing the convective shells within the turbulent regime. See Appendix B for a description of the effective Reynolds number.

Eight models are computed and are named according to the value used for $\epsilon_{fac}$, see Table 1. The model eps1k is an extension of the C17 hrez model up to 6000 s. All other models were computed from the same initial conditions. This allows one to study the effect that varying $\epsilon_{fac}$ has on the initial transient stage. In particular, the time required to reach a quasi-steady state, where an equilibrium is reached between the driving luminosity and numerical dissipation is of interest. Figure 2 shows the density (blue) and entropy (red) profiles for the 3D hydrodynamic model initial output and the corresponding 1D stellar model initial conditions (black). The highest energy model, eps33k, is difficult to evolve over long times as the energies are so extreme that the structure of the shell is disrupted, which results in the shell becoming dynamically unstable. As such, the physical simulation time for this model is relatively shorter than the other models. The global properties of each model are summarised in Table 1.

3 Results and discussions

3.1 General flow properties

The models first pass through an initial transient phase (whereby the turbulent velocity field is established) and the models eventually reach a quasi-steady turbulent state, where the initial conditions no longer influence the flow. The transition to the quasi-steady state for each model can be seen most easily in Fig. 3, where the specific total kinetic energy (KE; thick solid lines), $\mathrm{KE_{tot}}=(v_{x}^{2}+v_{y}^{2}+v_{z}^{2})/2$, for each model are plotted against the simulation time. The initial transient phase begins with a characteristic sharp rise in the KE up to a local maximum, which is then followed by a steady decrease. This marks the end of the initial transient, and the beginning of the quasi-steady state. From Fig. 3 it can be seen that the local maximum in KE during the transient phase increases with $\epsilon_{fac}$, and that the overall KE for each model increases as $\epsilon_{fac}$ is increased. This is intuitive as an increase in energy generation results in a larger flux of KE at the temperature peak near the bottom of the shell. The thin dotted lines in Fig. 3 show the specific total KE, $\mathrm{KE_{tot}}$ (thick solid lines) minus the specific total turbulent KE, $\mathrm{KE_{turb}}=(v_{x}^{\prime\,2}+v_{y}^{\prime\,2}+v_{z}^{\prime\,2})/2$ (see Appendix D for a description of the notation). These lines therefore represent the KE that is not associated with turbulent convection, i.e. KE due to waves propagating throughout the computational domain.

Specific KE spectra for each model are presented in Appendix C.

The logarithm of the specific KE evolution over the radius of the computational domain for the eps1k, eps3k, eps10k and eps33k models are shown as contour plots in Fig. 4. The colour-bars show values of log(KE) and are unique to each panel. Each panel is labelled for the respective model. Each model passes through an initial transient phase characterised by a strong pulse in KE near the start of each simulation. The models show semi-regular pulses in KE over their evolution, a characteristic of convective transport. Such strong turbulent motions also lead to turbulent entrainment. This is best seen at the upper boundary by the gradual migration of this boundary into the stable region above. In all models, gravity waves (excited by the convective flows hitting the boundary) can be seen in the upper stable region, identified by short, horizontal yellow streaks.

The entropy evolution over the radius of the computational domain for the eps1k, eps3k, eps10k and eps33k models are also shown as contour plots in Fig. 5. The colour-bars show values of the entropy in units of erg K*-1* and are unique to each panel. Each panel is labelled for the respective model. Turbulent entrainment is visible via the increasing extent of the high-entropy convective region.

The panels in Fig. 6 show contour plots for the final snapshot222The eps33k model snapshot is taken at 842 s rather than the final snapshot of the simulation, before the shell becomes disrupted. of the eps1k, eps3k, eps10k and eps33k models. Each contour plot shows the variation in entropy from the mean entropy in the convective region. This form of visualisation was chosen as the entropy in the convective region is nearly homogeneous. Entropy variations from the mean are seen in the convective region and at the boundaries, but are very small ($\sim 10^{-4}$, c.f. fig. 6 of Cristini et al. 2017). The magnitude and direction of the arrows represent the speed and direction of the $v_{x}$ and $v_{y}$ velocity fields within the x-y plane. With an increase in boosting factor it can be seen that the variation in entropy is increasing within the convective region, shown by darker red and blue hues. The comparative global distortion of the convective boundaries, in particular the upper boundary, clearly increases with an increasing boosting factor, as would be expected. Smaller scale distortions along the boundaries, such as Kelvin-Helmholtz instabilities, are also more abundant with increasing boosting factor.

Using dimensional analysis of a turbulent system, it can be shown that the energy dissipation rate is $\epsilon\sim v_{\rm rms}^{3}/\ell$, where $\ell$ is the integral length scale and $v_{\rm rms}$ represents the velocity of the largest energy-containing fluid elements (Kolmogorov, 1941). As the energy dissipation rate is set by the energy generation rate which is proportional to the luminosity, it can be shown that

[TABLE]

where $v_{\rm x,rms}$ is the vertical RMS velocity, assuming that the integral scale is constant and that $v_{\rm x,rms}\sim v_{\rm rms}$. Several other studies (Porter and Woodward, 2000; Arnett, Meakin and Young, 2009; Arnett et al., 2015; Jones et al., 2017) have also confirmed this proportionality.

For the eight models in this study, the radial (vertical) velocity does indeed have a positive correlation with the cube root of the nuclear energy generation rate. This is shown by Fig. 7, where the vertical RMS velocity of each model is plotted against the boosting factor, $\epsilon_{fac}$. A linear regression on these values reveals a line of best fit with a slope of $0.318\pm 0.034$. This is in agreement with Eq. 2.

The time-scale of the transient phase and the time required to establish the turbulent velocity field is approximately one convective turnover time (Meakin and Arnett, 2007, C17). The above scaling implies that, as the energy generation rate is increased while the initial structure is kept constant, the convective turnover time (given in Table 1) will decrease. This is confirmed in Fig. 3 by the relatively shorter initial transient times and shorter fluctuation333Such fluctuations are associated with a phase lag between the buoyant driving and dissipative damping (Arnett and Meakin, 2011b). periods during the quasi-steady phase for the more energetic models.

From above, assuming that the integral length scale is constant, it can be shown that the Mach number scales with the nuclear energy generation rate. So that, for constant sound speed the Mach number varies with the boosting factor as Ma $\propto\epsilon_{fac}^{1/3}$. The temporal evolution of the Mach number over the radius of the computational domain is plotted in Fig. 8 for models eps1k, eps3k, eps10k and eps33k. Again, each colour bar is unique to each model and shows the dimensionless Mach number. Similarly to Fig. 4, the episodic nature of turbulent convection and the generation of gravity mode waves in the upper stable region is well represented by the Mach number.

3.1.1 Model eps33k

Looking at Fig. 4, at first glance the eps33k model may appear to still be transitioning through the transient phase, but actually this transition is shown by the yellow region spanning the convective region up to around 200 s, followed by strong entrainment at the upper boundary. Towards the end of the model ($>700\,$s) a strong increase ($>10^{14}\,$erg g*-1*) in KE can be seen. At that time, the shell becomes dynamically unstable, and at later times ($>1000\,$s) the shell is completely disrupted. This can also be seen in the time series of the velocity magnitude snapshots for this model in Fig. 9, where the velocity magnitude is plotted at 340 s, 840 s, 1020 s and 1180 s. Comparing the velocity magnitudes at 340 s and 840 s (top two panels), it can be seen that the upper boundary has migrated, almost entirely encompassing the stable region above. By 1020 s (lower left panel) the upper boundary no longer exists and the entire previously upper stable region is now turbulent. These turbulent motions are no longer decelerated by approaching a stable region above, but still over-turn (as if these motions were approaching a boundary with an extremely high bulk Richardson number) due to the reflective boundary condition at the edge of the simulation domain. By 1180 s, the removal of the upper stable region, leads to the turbulent velocities increasing dramatically (lower right panel) as turbulence is still driven by a large nuclear energy generation rate at the bottom of the shell, and the up-flowing radial velocities produced from this turbulence are not reduced by a negative buoyancy force due to the presence of an upper boundary. Instead the strong radial deceleration and turning of fluid elements at the upper domain boundary results in a band of material ($\geqslant 2\times 10^{9}\,$cm) with visibly reduced velocities. Regardless, the removal of the boundary results in an un-damped, runaway acceleration within the convective region and the velocities become so high, that the fluid approaches the transonic regime ($|v|\sim 10^{8}\,$cm s*-1*) and the shell becomes dynamically unstable. This is shown in the bottom right panel of Fig. 8, where the Mach number can be seen to exceed $0.25$. Unsurprisingly, this demonstrates that for the given background stratification there is an upper limit to the amount of artificial boosting that can be applied to the model without complete shell disruption. Despite such dynamic and violent behaviour, the short time period between $450\,$s - $900\,$s of the eps33k model still represents turbulent entrainment at a greatly boosted rate, this time window is used in further boundary analysis for this model.

3.2 Temporal evolution

The importance of including several convective turnovers for temporal averaging of quantities is demonstrated in Fig. 10, where the RMS radial velocity profile (normalised to the velocity in the centre of the domain) is plotted using various temporal averaging windows from one to ten convective turnovers. An instantaneous profile (no time averaging) is also plotted as a black dashed line to illustrate the variance in velocity over a single snapshot in time. From Fig. 10 it can be seen that an averaging window of one convective turnover smoothes out the stochasticity of the instantaneous profile, but is qualitatively different from the remaining time windows.

The number of convective turnovers completed during the simulation time of each model is shown in Fig. 11, where it can be seen that the eps1, eps33, eps100 and eps333 models have all completed two or less convective turnovers. In such models, the turbulent velocity field may have developed but there is insufficient time during the quasi-steady state to provide statistically reliable results.

Therefore, the four models eps1, eps33, eps100 and eps333 have not completed a sufficient amount of time in the quasi-steady phase to warrant a full analysis. These models were thus excluded from the boundary analysis involving turbulent entrainment. It must be stressed that the poor temporal convergence of these models is not a result of the models themselves, but due to a lack of sufficient computational resources. We plan to evolve these models out to much later times ($>10$ convective turnovers) in the future.

The longer time-scales of the remaining models ($\gtrsim 4$ convective turnovers), leads to better converged properties over the longer quasi-steady states. An analysis of these models within the Reynolds-averaged Navier-Stokes framework and the effect that varying $\epsilon_{fac}$ has on the mean turbulent kinetic energy equation are presented in Appendix D.

The convective boundaries of these models are also in a quasi-equilibrated state, where the growth of the boundaries due to turbulent entrainment occurs within an equilibrium regime. This regime is known as the equilibrium entrainment regime (Fedorovich, Conzemius and Mironov, 2004; Garcia and Mellado, 2014) and is achieved when the time-scale for boundary growth exceeds the convective turnover time-scale. In such conditions the boundary evolution is quasi-equilibrated as the entrainment process samples the entire spectrum of turbulent motions, rather than strong, individual plumes. The definitions for the two time-scales mentioned here and the respective values for the various models are discussed in Appendix E.

3.3 Turbulent entrainment at convective boundaries

In order to determine the positions of the convective boundaries, we adopt the method of C17, who calculate the positions based on the difference in composition between the adjacent convective and radiative regions. This provides a two-dimensional boundary surface. The estimated radial position of the boundary is then obtained through a horizontal mean of the radial position over this surface. These boundary positions over the simulated times of the eps1k, eps3k, eps10k and eps33k models are plotted in Fig. 12. The positions of the upper (top panel) and lower (bottom panel) boundaries are shown by solid lines, and the shaded envelopes represent twice the standard deviation from the calculated means. The models show clear entrainment from the evolution of the boundary position over time.

The variance in the boundary positions increases as the driving luminosity is increased. This agrees with the notion that the boundaries become softer with an increased driving luminosity due to a higher TKE flux at the boundaries, and that a larger distortion in the boundary surface (relative to the mean) is characteristic of more plume penetration.

The evolution of the convective boundary positions for the eps1k (top left), eps3k (bottom left), eps10k (top right) and eps33k (bottom right) models are also shown separately in Fig. 13. These plots are very similar to Fig. 12, but allow an individual assessment of the entrainment and boundary migration for each model. The positions of the boundaries are shown by black solid lines. The best fit line from a linear regression of the boundary positions during the quasi-steady state for each model is shown by a solid coloured line. The corresponding best fit slope and relative error are shown in coloured text beside each fit. This slope is the entrainment velocity, and is remarkably close to linear in each case. The small errors in the best fit slope show this. The largest is 1.5% for the upper boundary of the eps33k model.

3.3.1 The entrainment law and boundary stiffness scaling

The time rate of change of the boundary position through turbulent entrainment (the entrainment velocity), $v_{e}$, is known to have a power law dependence on the bulk Richardson number (Eq. 24; see Appendix F for details). In this equation, the ratio of the entrainment velocity to the velocity representing the large-scale fluid elements, $E=v_{e}/v_{\rm rms}$, is considered (e. g. Garcia and Mellado, 2014). This relationship between the relative entrainment rate and the bulk Richardson number is referred to throughout the meteorological and atmospheric science communities as the entrainment law, and is typically written as:

[TABLE]

where $A$ and $n$ are constants. Simulations (e. g. Deardorff 1980) and laboratory (e. g. Chemel, Staquet and Chollet 2010) studies have found similar values for the coefficient, $A$, typically between 0.2 and 0.25. The exponent, $n$, is generally taken to be $1$. On the other hand, in a recent DNS study, Jonker et al. (2013) showed that $A\approx 0.35$ and $n=1/2$.

The eps1k, eps3k, eps10k and eps33k models have been interpreted within the framework of the entrainment law (Eq. 3). These models cover a significant period ($\gtrsim 4$ convective turnover times) in the quasi-steady state (see Fig. 11). The entrainment speed (normalised by the RMS turbulent velocity) is plotted as a function of the bulk Richardson number for these models in Fig. 14. The coloured circles and diamonds represent the values for the upper and lower boundaries, respectively, for the eps1k (yellow), eps3k (cyan), eps10k (orange) and eps33k (blue) models.

The solid line is the best fit line following a linear regression of all the data points in (logarithmic) $E-\textrm{Ri}_{B}$ space. The line of best fit is then of the form $y=mx+c$, where $m=-n$ and $c=\textrm{log}A$. The slope ($n$) and intercept ($A$) along with the respective errors for this best fit are noted to the left of the line as $n=0.74\pm 0.04$ and $A=0.05\pm 0.06$, respectively. The black points in Fig. 14 represent the entrainment rate, $E$ (as a function of the bulk Richardson number) obtained in the oxygen shell burning study by Meakin and Arnett (2007), and the dashed line is the best fit curve following a linear regression of their data points. The slope and intercept are also noted beside this fit as $n=1.05\pm 0.21$ and $A=1.06\pm 2.40$, respectively.

Comparing the values obtained for the constants in this luminosity study with those obtained in the resolution study by C17 (see their fig. 15), it can be seen that the values for $n$ are in reasonable agreement, and point towards a value in the range $1/2\leq n\leq 1$ for fusion driven (neutrino cooling dominated) turbulent convection in the carbon burning shell of massive stars. The value of $A$ is slightly smaller than the values quoted in other laboratory and numerical studies. The lower value of $A$ as compared with that of the oxygen shell burning simulation (black points) suggests that the efficiency of the work done by turbulent eddies on the stable stratification at the boundary of the carbon shell is less than that of the oxygen shell. This being said, the error in the calculation of $A$ for the oxygen shell is large, it stands to reason then that both fits could lead to the same value of $A$, within the error range of the simulations.

It is not clear at this stage if a single set of parameters ($A$ and $n$) in the entrainment law would apply to all burning stages. In order to reduce the errors in the calculation of both the bulk Richardson number and the entrainment rate, models which span larger fractions of the quasi-steady state are required. The variance of $n$ between the carbon and oxygen shells suggests that there may be additional mixing processes occurring besides entrainment in the carbon shell. Uncertainties in the ‘correct’ value of $n$ for turbulent entrainment are also present throughout terrestrial simulations, where $n=1/2$ was found by Jonker et al. (2013) while $n=1$ was found by Fernando (1991), for example.

Furthermore, an approximate scaling relation can be obtained between the bulk Richardson number and the driving luminosity. This can then allow the determination of the stiffness of convective boundaries in 1D stellar models. This relation can be obtained by starting with the formula for the bulk Richardson number (Eq. 24) and substituting for $v_{\rm rms}$ using Eq. 2 leading to

[TABLE]

if the integral scale and buoyancy jump ($\Delta B$) are assumed to be constant (reasonable assumptions given an initial hydrostatic stratification and short dynamical time-scales).

This scaling between the driving luminosity and bulk Richardson number can easily be tested within the current luminosity study. As previously mentioned, only the four models eps1k - eps33k are used to test this relation, as the remaining models do not contain a sufficient number of convective turnovers over which entrainment can occur. The bulk Richardson numbers for the upper (bottom line) and lower (top line) boundaries of these four models are plotted in Fig. 15 as a function of the boosting factor, $\epsilon_{fac}$. A linear regression of the data points for each boundary is performed over logarithmic space in order to determine a best fit power law and scaling exponent, $\beta$, assuming Ri${}_{B}\propto L^{\beta}$. This best fit power law is shown by the grey line for each boundary and the corresponding value of the scaling exponent is noted beside each one, $\beta=-0.686$ for the upper boundary and $\beta=-0.648$ for the lower boundary. For both boundaries, the value obtained is close to the expected value of $\beta=-0.667$. Unfortunately, due to such a sparse data set it is not possible to obtain errors in the values of the slope. Nevertheless, such an agreement between theory and simulation is encouraging and with more data points and longer time sampling for each model the confidence in these estimates can be improved.

From the best fit power laws calculated in Fig. 15, the bulk Richardson numbers can be extrapolated back to the original carbon burning luminosities of the 1D input model. These are Ri${}_{B,u}=4750$ and Ri${}_{B,l}=3.77\times 10^{4}$ for the upper and lower boundaries, respectively. These values are in agreement (to within a factor of around 3) with the values calculated in the eps1 model (see Table 1), although arguably, the eps1 model values are not well constrained either due to the short number of convective turnovers simulated (e.g. see Fig. 11). The above values also agree (also to within a factor of around 3) with the values calculated from the stellar model initial conditions (see table A1 of C17), which also have inherent uncertainties, due to the difficulty of calculating the bulk Richardson numbers over coarse stellar model grids.

The entrainment velocities for the eps1 model can also be estimated using the above values for the bulk Richardson numbers. Inserting the above values for the bulk Richardson numbers into the entrainment law (Eq. 3) with the values for $A$ and $n$ determined in Fig. 14 and the RMS velocity from Table 1 yields entrainment velocities of around 14 cm s*-1* and 65 cm s*-1* for the lower and upper boundaries, respectively. The corresponding mass entrainment rates can be estimated using, $\dot{m}_{e}=4\pi r^{2}\rho\,v_{e}$, where the density is that of the boundary region averaged over the quasi-steady state. The mass entrainment rates for the lower and upper boundaries of the eps1 model are $1.35\times 10^{-8}$ M${}_{\odot}\,$s*-1* and $3.77\times 10^{-8}$ M${}_{\odot}\,$s*-1*, respectively. Such mass entrainment rates imply that the time required for a mass comparable to that of the carbon shell ($\sim 0.75\,$M⊙) to be entrained is roughly half a year. This time-scale seems rather short though considering that the evolutionary time-scale for the carbon shell used as input is roughly 30 years, suggesting that the inferred mass entrainment rate is too large compared to what is expected from stellar evolution models.

It is very important at this point to remember that the initial 1D model used in this study corresponds to the start of the carbon shell burning phase, during which the convective region also grows significantly in the 1D calculation. We do not expect the mass entrainment rate derived in this study to apply for the entire carbon burning stage nor to other burning stages. The Bulk Richardson number is expected to evolve during a burning stage and between burning stages. It is expected to increase significantly as the entropy and mean molecular weight gradients at the convective boundaries build up throughout the burning stage. At the very end of a burning stage, the driving luminosity will decrease and entrainment is expected to cease altogether. This is why a theoretical framework such as the entrainment law is crucial in order to correctly apply results of 3D hydrodynamic simulations back into 1D stellar evolution models and to avoid unrealistic growths of convective zones.

Despite the difficulties discussed above, the reliable estimation of the bulk Richardson number is still best pursued through long time-scale (ideally $\gtrsim 10$ convective turnovers), 3D hydrodynamic simulations of turbulent convection. The results of which can populate the entrainment law parameter space ($E-\textrm{Ri}_{B}$) and constrain the free parameters $A$ and $n$. The agreement with the scaling relation, Ri${}_{B}\propto L^{-2/3}$, for several models over a range of different driving luminosities, demonstrates the applicability of this scaling relation to stellar flows. The discrepancy between the extrapolated values of the bulk Richardson number for natural carbon burning luminosities and the stellar model initial conditions, could be used to fit free parameters within the formulation of the bulk Richardson number for stellar models, namely the integration length, $\Delta r$, and the integral scale, $\ell$ (see Appendix F).

With better constrained values of $A$ and $n$ as well as more accurate estimates of the bulk Richardson numbers within stellar models, new mixing prescriptions can be incorporated into stellar models, such as the one suggested by Meakin and Arnett (2007)

[TABLE]

3.4 Effects of driving luminosity on boundary widths

The composition profiles at both boundaries for the final time-step of each of the four models (eps1k - eps33k) are shown in Fig. 16. Each composition profile has been shifted in radius such that the boundary position coincides with that of the eps1k model (for ease of comparison). Boundary region widths are calculated through the adoption of the method used by C17 (see their §4.5.3). This method uses the jump in composition, $\bar{A}$, between the convective and stable regions. The boundary region is estimated to include all but 1% of the composition contained within the convective region and surrounding stable region.

In Fig. 16, the boundary width is denoted for each model by the distance between two filled squares at the edges of the boundary region, and is shown by the respectively coloured shaded region. These boundary widths and their fractions of the pressure scale height are given in Table 2.

For the upper boundary (right panel of Fig. 16) the composition profiles are similar except for the eps33k model, which has a shallower slope over the boundary region. In this model the upper boundary is much broader, and the increased driving leads to a dramatic smoothing in the abundance slope; this is not seen in the other models. As noted in §3.1.1, the strong increase in driving luminosity, is likely leading to a structural readjustment of the shell at this time-step and will eventually result in the complete disruption of the shell.

At the lower boundary (left panel) the abundance profiles are very similar in slope between all of the different cases, with an increase in driving luminosity leading to a clear increase in boundary width. As entrainment becomes more effective with stronger driving luminosity, the composition in the convective region near the boundary also increases slightly. The increased amount of heavier material is due to mixing of material from the stable region below into the turbulent region.

The smoothing of the horizontally averaged abundance profile due to the boundary deformation can be measured by the standard deviation of the boundary position (Fig. 12). The deviation can be associated with the vertical fluctuations in the boundary surface, possibly due to plume penetration events. The extent of the boundary region (its width) can be associated with the strength of shear mixing (or Kelvin-Helmholtz instability) at the boundary caused by the U-turning of turbulent fluid elements. The boundary widths (Table 2) have a weaker dependence on the boosting factor, $\epsilon_{fac}$, than originally expected, but it can be said that the boundary width generally increases with boosting factor, leading to a better spatial sampling of the boundaries. The eps33k model with a boosting factor of $3.3\times 10^{4}$ is the only model where both convective boundaries are spatially resolved at a resolution of $512^{3}$ (see the bottom left panel of Fig. 18 and Fig. 19). This is encouraging as it implies that there exists a resolution where the ILES method adopted here can provide a spatially resolved model of the carbon shell and its boundaries at nominal luminosity.

4 Conclusions

In this paper, we studied the dependence of convective boundaries on their stiffness and the turbulence strength. To study this dependence in a controlled way, we performed a series of simulations of the carbon burning shell, in which the nuclear energy generation rate was artificially boosted. The initial conditions were taken from the 15 M⊙ stellar model described in C17. Each model was computed within a Cartesian cube consisting of 5123 computational zones. The nuclear energy generation rate was calculated using the same prescription described in C17, with the addition of a boosting factor, $\epsilon_{fac}$, which we varied between 1 and $3\times 10^{4}$. The results showed that over the short dynamical time-scales, which these models were run over (due to limited computational resources), the shell structure was not adversely affected by the boosting (except in the most extreme case) and served to accelerate the evolution of the shell through increased rates of entrainment and mixing. As the structure was unchanged between each model, the vigour of turbulence was also increased for the same stable stratification and so the representation of entrainment in the most energetic cases ($\epsilon_{fac}\geq 10^{4}$) was far from been physically realistic during the adopted evolutionary state.

Such a luminosity study allowed us to populate the parameter space ($E-\textrm{Ri}_{B}$) of the entrainment law (Eq. 3) and constrain the values of $A$ and $n$ for high Péclet number, fusion driven convection. We found values of $A=0.05\pm 0.06$ and $n=0.74\pm 0.04$. The central result of this study is the strikingly clear dependence of the entrainment rate on the stiffness of the convective boundary. This result was found for both the top and bottom boundaries. The bottom boundary being stiffer than the top one, entrainment at the bottom boundary is slower. While this was already observed in our previous study C17, the present luminosity study places this finding on a firm footing.

Comparing to the values found in the resolution study by C17, a similar value for $n$ was found, and further points to a value of $n$ for fusion driven (neutrino cooling dominated) turbulent convection in the range $1/2\leq n\leq 1$. The variance of $n$ between the carbon and oxygen burning shells (Meakin and Arnett, 2007) suggests that there may be additional mixing processes occurring besides entrainment in the carbon shell. Comparison to the simulations of Jonker et al. (2013) suggests that these entrainment values may depend upon the Péclet number. Further, the structure of the convective region changes as evolution proceeds. While the entrainment law seems valid for times of the order of the turnover time, it may not be universal for evolutionary time scales and conditions.

Further simulations of stellar convection are needed to help populate the entrainment law parameter space. Longer time-scale simulations ($>10\,\tau_{c}$) are also needed to provide better statistics during the quasi-steady state and reduce the errors in the calculated values of $A$ and $n$. As computing power increases, and the porting of existing codes to GPU technologies is becoming more common, more high-resolution, longer time-scale simulations will be possible.

This luminosity study also confirmed the scaling relation between the vertical RMS velocity and the driving luminosity, $v_{\rm x,rms}\propto L^{1/3}$, as found in previous studies. This study also confirmed the expected scaling between the bulk Richardson number (stiffness) of the boundary and the driving luminosity of the shell, Ri${}_{B}\propto L^{-2/3}$. Such a relation will prove useful in developing new CBM prescriptions for entrainment in stellar evolution models. For example, it helps to extrapolate values of the bulk Richardson number back to the nominal carbon burning luminosities. Comparing these extrapolated values to those from the initial 1D stellar evolution model calculations in C17 reveals that they are similar to within a factor of around 3.

This study represents a few steps along the very long road of integrating the results of 3D hydrodynamics simulations into a theoretical framework that can be used in 1D stellar evolution models; framework coined “321D” by John Lattanzio and others. Our simulations give a plausible representation for the physics of turbulent stellar convection, but only for a few turnover times. Longer term behaviour, which is crucial for stellar evolution properties such as the size of mixing regions, remains a challenge. The evolution of the convective region may modify the entrainment law, which would depend upon the properties of the convection zone (e.g., luminosity, entropy jumps, composition jumps, etc.). The problem may be a coupled one, so that there might not be a universal entrainment law as suggested by the different values derived for $A$ and $n$ between carbon and oxygen burning. Despite this potential non-universality, 3D simulations at various points throughout the evolution of stars may still be able to constrain the values of $A$ and $n$ to be used with the entrainment law in 1D models.

Including entrainment (or the penetrative overshoot commonly used in GENEC and other 1D models) moves the convective boundaries but does not change their shape. As shown in Fig. 16, abundance profiles at convective boundaries in 3D simulations are smooth, sigmoid-like functions rather than the step-like functions found in GENEC and other 1D codes. These smoother shapes are due to convection-induced shear mixing (probably Kelvin-Helmholtz instability) occurring at the boundary as the fast upward convective flow U-turns. Note that rotation-induced shear mixing (Edelmann et al., 2017) plays a similar role as the convection-induced shear and both lead to smoother profiles at boundaries. Using a diffusion coefficient, which is exponentially decaying from the convective boundary in 1D codes following multi-D simulations of outer convective regions (Herwig, 2000; Freytag, Ludwig and Steffen, 1996) commonly used in MESA (Paxton et al., 2011), enables the reproduction of asteroseismic constraints and leads to abundance profiles at boundaries, which are comparable to those of 3D simulations (Arnett and Moravveji, 2017). This exponentially-decaying mixing leads to larger convective cores like penetrative overshoot or entrainment but at a rate/extent which may not be correct. We may thus have to include the growth rate of convective regions (e.g. via entrainment) and the process that determines the shape of their boundaries (e.g. shear mixing) as separate processes in 1D codes in order to reproduce all the properties of convective boundaries found in 3D hydrodynamic simulations. Further steps along the 321D road include replacing the mixing-length theory (Böhm-Vitense, 1958) with a theory able to predict convective properties such as velocities in the entrained region (see discussion in Arnett et al., 2015).

Acknowledgements

The authors would first like to thank the referee for their continued questioning, concerns and suggestions, resulting in a greatly improved manuscript. The authors acknowledge support from EU-FP7-ERC-2012-St Grant 306901. RH acknowledges support from the World Premier International Research Centre Initiative (WPI Initiative), MEXT, Japan. This article is based upon work from the “ChETEC” COST Action (CA16117), supported by COST (European Cooperation in Science and Technology). AC acknowledges partial support from NASA Grant NNX17AG24G. This research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575. CM and WDA acknowledge support from NSF grant 1107445 at the University of Arizona. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin (http://www.tacc.utexas.edu) for providing HPC resources that have contributed to the research results reported within this paper. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. The authors acknowledge PRACE for awarding access to the resource MareNostrum 4 based in Spain at Barcelona Supercomputing Center. The support of David Vicente and Janko Strassburg from Barcelona Supercomputing Center, Spain, to the technical work is gratefully acknowledged. CG, RH, and CM acknowledge ISSI, Bern, for its support in organising a collaboration meeting which sparked many fruitful discussions. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

Appendix A Computational tools

Initial conditions for the 3D hydrodynamic simulations are provided by a stellar evolution model of a 15 M⊙ star evolved up to the point of carbon shell burning. This model was calculated using the Geneva stellar evolution code (genec). The 3D hydrodynamic simulations presented here were calculated using the Prometheus MPI (prompi; where MPI is an acronym for message passing interface) code. Both codes are described below.

The 1D stellar model structure and composition were mapped onto a finer grid for use in the prompi code. These finely zoned stellar model data can be found in full in machine-readable format as supplementary online material. Several zones from these initial conditions are also shown in Table 3.

A.1 The Geneva stellar evolution code

genec (Eggenberger et al., 2008; Ekström et al., 2012) solves the stellar evolution equations (see e.g. Eqs. 10.1 - 10.5 of Kippenhahn, Weigert and Weiss 2013) using a finite difference, time implicit method within a Lagrangian framework.

The chemical composition is homogeneously mixed in convective regions up to oxygen burning. The structure equations are de-coupled from the abundance equations; changes in abundances due to nuclear burning and diffusive mixing are calculated separately.

Reaction rates are calculated for a nuclear reaction network of 23 isotopes. Energy losses due to the production and loss of neutrinos is included. A perfect gas equation of state including radiation and partial degeneracy is used. Opacities are interpolated from tables provided by the OPAL group (Alexander and Ferguson, 1994; Rogers, Swenson and Iglesias, 1996). Mass loss estimates are a function of metallicity, and calculated according to the prescriptions by Vink, de Koter and Lamers (2001) and de Jager, Nieuwenhuijzen and van der Hucht (1988). Energy transport due to convection is assumed to be adiabatic for deep internal convection, and for convective envelopes, energy transport is modelled using the mixing length theory with $\alpha_{ml}=1.6$ (Schaller et al., 1992). The extent of convectively unstable regions, is determined by the Schwarzschild criterion, along with penetrative convection of up to 10% of the local pressure scale height (for core hydrogen and helium burning phases), where the pressure scale height is $H_{p}=-dr/d\,\textrm{ln}\,p=p/g\rho$.

A.2 The Prometheus MPI code

prompi (Meakin and Arnett, 2007) is a finite-volume, time explicit, Eulerian code which uses the piecewise-parabolic method of Colella and Woodward (1984).

A.2.1 Piecewise-Parabolic Method

The piecewise-parabolic method (PPM) is a higher order extension of Godunov’s method (Godunov, 1959; Godunov, Zabrodin and Prokopov, 1962), which utilises a forward Euler method for time integration444Note that, PPM is a special method in the sense that it has a higher effective resolution than other differencing schemes, e.g. DNS (see e.g. Sytine et al., 2000, for a discussion). In multiple dimensions PPM is second order accurate in space and time (Colella, 1990). PPM provides a Kolmogorov description of the turbulent cascade down to the dissipation level at the sub-grid level. It was developed with shock capturing in mind, and relates the change in turbulent kinetic energy and traversal time across a shock to the specific entropy production rate (Kolmogorov, 1962; Grinstein, Margolin and Rider, 2007), or

[TABLE]

This is the definition for the energy dissipation rate along a turbulent cascade (in a statistically steady state, see also Appendix B). Without the use of an explicitly defined viscosity PPM matches motions on the turbulent cascade to those at the grid scale. In the ILES paradigm then, the physical scale at which dissipation occurs is replaced by the grid scale. The conservation laws are enforced to machine accuracy, and so the turbulent cascade is well represented (see Fig. 17 and fig. 8 of Cristini et al. 2017), providing that the inertial range is sufficiently sampled on the grid.

A.2.2 prompi design

prompi is derived from the legacy astrophysics code prometheus (Fryxell, Müller and Arnett, 1989). prompi utilises domain decomposition in order to be parallelised over multiple processors. Communication between computational nodes is handled by MPI.

The Euler equations for fluid motion (inviscid approximation) are solved within the ILES paradigm (e.g. Grinstein, Margolin and Rider, 2007),

[TABLE]

where $p$ is the pressure, $\mathbf{g}$ the gravitational acceleration, $E_{t}$ the total energy, $X_{i}$ the mass fraction of nuclear species $i$ and $R_{i}$ the rate of change of nuclear species $i$.

A.2.3 Model set-up

Kolmogorov (1941) showed that the rate of TKE dissipation at all scales does not depend on the details of the dissipative process. This implies that it may be unnecessary to resolve the dissipation sub-range of the cascade, meaning that the sub-grid dissipation becomes the effective Kolmogorov length, through its implied viscosity. In a turbulent flow the momentum diffusion (viscosity) is negligible in comparison to the advection of KE at all scales except near the Kolmogorov scale and below it. In ILES, the minimum amount of dissipation required in order to maintain monotonicity and energy conservation is adopted.

The implicit sub-grid scale model used in ILES is the leading order term in the truncation error, which is a result of discretising the Euler equations through the use of PPM, and is therefore of second order accuracy555Unlike other LES methods, ILES does not require the use of an explicit sub-grid model.. The qualitative effects of dissipation at the grid scale are implemented in ILES. Energy conservation is built into the ILES scheme; KE that cascades down to the sub-grid scales from the resolved scales is damped as the velocity fluctuations are dissipated. The internal energy of the fluid is suitably increased in such a way that it mimics viscosity at the Kolmogorov scale, which would dissipate the structures at this scale into heat.

The specific form of this numerical diffusion is the same as the simplest case described in the appendix of Colella and Woodward (1984). That is the interpolation functions are flattened based on the steepness of the pressure jump across a zone and an extra explicit diffusive flux is also added. In multiple dimensions, this method is valid providing that the same flattening is applied to the derivatives in each direction (Colella, 1990).

A Cartesian geometry is used, and the boundary conditions are reflective in the vertical direction and periodic in the horizontal directions. Velocities are also damped near the domain boundary in the vertical direction using a sinusoidal function. Explicit time-stepping is used to ensure that all structural and acoustic changes within cells are temporally resolved.

The time-step, $\Delta t$, is restricted such that it is in accordance with the CFL condition (Courant, Friedrichs and Lewy, 1928), specifically

[TABLE]

where $\Delta x$ is the cell width, $c_{s}$ is the local sound speed and $C_{\rm max}$ is the Courant factor set to be 0.8.

The Helmholtz EOS (Timmes and Arnett, 1999; Timmes and Swesty, 2000) is used to describe the thermodynamic state of the plasma. Convective heat transport is initiated in the model through small, random and equal perturbations in temperature and density. Energy generation from nuclear fusion is parameterised as a function of composition, density and temperature, using a slightly modified version of the parameterisation given by Audouze, Chiosi and Woosley (1986); Maeder (2009):

[TABLE]

where $Y_{12}=X_{{}^{12}\rm C}/12$, $X_{{}^{12}\rm C}$ is the mass fraction of 12C, $\lambda_{12,12}=5.2\times 10^{-11}\,T_{9}\hskip 0.14226pt^{30}$ and $T_{9}=T/10^{9}$. The loss of energy due to escaping neutrinos is parameterised using the analytical formula provided by Beaudet, Petrosian and Salpeter (1967),

[TABLE]

where $\lambda=k_{B}T/mc^{2}$ and $\xi=\rho^{1/3}\lambda^{-1}$.

Appendix B Effective Reynolds number derivation

A useful dimensionless number for determining the extent of turbulence in a simulation is the effective Reynolds number. This is the discrete analogue of the Reynolds number, $\textrm{Re}=v\,\ell/\nu$, where $v$ and $\ell$ are the velocity and length scales and $\nu$ is the viscosity. The effective Reynolds number is defined using the following arguments. The rate of energy dissipation at a length scale $\lambda$ is $\epsilon\sim v^{3}/\lambda$ (Kolmogorov, 1941). One can then approximate the rates of energy dissipation at the extreme scales of the simulation, i.e. at the integral scale and the grid scale. These energy dissipation rates are

[TABLE]

where $\Delta v$ is the flow velocity across a grid cell; this velocity can also be used to define an effective numerical viscosity at the grid scale

[TABLE]

For a turbulent system within a statistically steady state, Kolmogorov (1962) showed that the rate of energy dissipation is equal at all scales; applying this equality to Eq. 13 yields (with the use of Eq. 14)

[TABLE]

Assuming the integral scale is the size of the convective region, $\ell_{cz}$, then the effective Reynolds number can be expressed as (to within a factor of 2)

[TABLE]

where $N_{x}$ is the number of grid points in the vertical direction. In these simulations this is a slight over-estimate as in the vertical direction only half of the grid points represent the convective region. Therefore, for these models which have a vertical resolution of $N_{x}=512$ zones (256 in the convective zone), the effective Reynolds number is, Re${}_{\rm eff}\sim 1625-4000$, which is comfortably within the turbulent regime.

Appendix C Specific kinetic energy spectra

The properties of the inertial sub-range of scales within the turbulent cascade are explored by comparing velocity spectra for each model at different luminosities. This is achieved by performing a horizontal 2D fast Fourier transform666Using the Python package numpy.fft.fft2. (FFT) of the vertical velocity at a constant height, within the centre of the convection zone. The results of this transform are presented in Fig. 17, where the transform squared, $\hat{V}(k)^{2}$, is plotted as a function of the wave-number. These spectra are time-averaged over several convective turnover times, where the convective turnover time is given by, $\tau_{c}=2\ell_{cz}/v_{\rm rms}$, with $\ell_{cz}$ being the height of the convective region (see Table 1 for values of $\tau_{c}$). The 1D profile is obtained by binning the 2D transform within the $k_{y}-k_{z}$ plane, where $k_{y}$ and $k_{z}$ are the wave-numbers in the y and z directions, respectively ($k_{y},k_{z}=0,2\pi,4\pi,...,2\pi(N/2),$ where $N$ is the number of grid points in one dimension, i.e. the resolution). A scaling of $(1/N_{y}N_{z})$ is applied, where $N_{y}$ and $N_{z}$ are the two horizontal resolutions. The spectra have also been normalised by $k^{5/3}$ in order to depict the regions of the spectra that follow Kolmogorov’s $k^{-5/3}$ power law (Kolmogorov, 1941) as roughly horizontal. This region is roughly located between the two vertical dashed lines at approximate wave numbers of $k=20$ and $k=250$. This wave-number range represents the scales of the models that are no longer affected by the initial or boundary conditions, driving force or dissipation. We also attribute the leftmost region as the integral range and the rightmost region as the dissipation range. The magnitude of the specific kinetic energy (velocity squared) increases with driving luminosity, as would be expected.

Velocity power spectra at various resolutions up to $1024^{3}$ are calculated by Sytine et al. (2000) for both PPM simulations solving the Euler equations and up to $512^{3}$ solving the Navier-Stokes (NS) equations with explicit viscosity terms. They showed that by comparing the kinetic energy time dependence that the PPM solutions required around 4 times less resolution than the NS solutions to model the same system.

Appendix D Reynolds-averaged Navier Stokes framework

It is common, when studying turbulent flows, to average the governing equations both spatially, to obtain a mean turbulent state over two dimensions, and temporally, to smooth out the stochastic nature of turbulence and provide a statistical average. When taking a Reynolds average of the Euler equations, a new ‘mean evolution’ of the fluid flow can be represented by averaging the horizontal components into a radial one.

Reynolds decomposition by construction separates the mean flow component from the fluctuating component. One can then construct physically relevant terms, which represent competing processes for a given conservation law (Ch.5 of Chassaing, 2002). The mean field is calculated by averaging over the horizontal plane, i.e. for a quantity $\omega$

[TABLE]

where $\textrm{d}A=\textrm{d}y\,\textrm{d}z$ and $\Delta A=\Delta y\,\Delta z$ is the area of the computational domain perpendicular to the radial axis. The original quantity can then be expressed as the sum of the mean and fluctuating components

[TABLE]

where the Reynolds average of the fluctuation is by definition zero, i.e. $\overline{\left<\omega^{\prime}\right>}=0$. The overbar notation denotes a temporal average over the quasi-steady state. This provides a better, statistically more meaningful representation of the flow. It also smooths out many of the instantaneous fluctuations provided the number of convective turnovers simulated is large enough ($\gtrsim 4$). This temporal averaging is defined as:

[TABLE]

for an averaging window $\Delta t=t_{2}-t_{1}$.

D.1 Turbulent Kinetic Energy Equation

The Eulerian equation of turbulent kinetic energy can be written as (eq. A12 of Meakin and Arnett 2007):

[TABLE]

where $v$ is the velocity and $E_{k}=\frac{1}{2}(v\cdot v)$ is the specific kinetic energy.

Applying horizontal and temporal averaging to the above equation yields the mean turbulent kinetic energy equation, which can be written as,

[TABLE]

where $\mathbf{D_{t}}=\partial_{t}+\boldsymbol{v}\cdot\boldsymbol{\nabla}$ is the material derivative;

$\mathbf{F_{p}}=p^{\prime}\boldsymbol{v^{\prime}}$ is the turbulent pressure flux;

$\mathbf{F_{k}}=\rho E_{k}\boldsymbol{v^{\prime}}$ is the TKE flux;

$\mathbf{W_{p}}=p^{\prime}\boldsymbol{\nabla}\cdot\boldsymbol{v^{\prime}}$ is the pressure dilatation;

$\mathbf{W_{b}}=\rho^{\prime}\mathbf{g}\cdot\boldsymbol{v^{\prime}}$ is the work due to buoyancy; and

$\epsilon_{k}$ is the implied numerical dissipation of KE.

D.2 Effects of driving luminosity on the turbulent kinetic energy budget

The turbulent kinetic energy (TKE) budget of the models in this driving luminosity study are interpreted within the Reynolds-averaged Navier-Stokes (RANS) framework. All of the terms in Eq. D.1 were calculated for the eps1k - eps33k models and are shown in Fig. 18. The radial profiles of each term are presented in the left panels of this figure, with the inferred viscous dissipation shown by a dashed line. This numerical dissipation is calculated from the residual TKE ($\epsilon_{k}$ in Eq. D.1). For a well resolved turbulent system, it should be a smooth profile. Each profile is time averaged over multiple convective turnover times and normalised by the surface area of the domain. The right panels show bar charts of the radial integration of each term shown in the left panel. It can be seen that all profiles increase in magnitude as the driving luminosity is increased, simply because more nuclear energy is being put into the system.

The main driving term for the TKE is the buoyancy work, $\mathbf{W_{b}}$, which for the eps1k, eps3k and eps10k models is balanced by the numerical dissipation at the grid scale, $\epsilon_{k}$. Due to such a balance, the shell is in a statistically steady state, as shown by the very small values of the material time derivative of the TKE, $\rho\mathbf{D_{t}}E_{k}$. For the eps33k model, the driving term is so large that while most of the TKE is dissipated, some of the energy affects the shell dynamically. Towards the end of the simulation, it is no longer in a statistically steady state and the shell itself is eventually completely disrupted.

The peak in the residual TKE (dashed curve in Fig. 18) at the bottom of the convective shell decreases in amplitude with increased driving luminosity, and the profile also broadens with increased driving luminosity. This is mainly due to the broadening of the boundary itself as the boundary is more easily overwhelmed by turbulent motions with an increased driving luminosity. For the most energetic model, eps33k, the lower boundary is broad enough that the spatial resolution is sufficient to model this boundary physically and without any adverse numerical effects due to poor resolution.

To investigate this in more detail, the residual TKE is compared between these four models for the upper and lower boundary in Fig. 19. In this figure the numerical dissipation is normalised by a value at a common position within the convective region close to the boundary. This normalisation value is numerically converged at these spatial resolutions (5123) and so poor spatial resolution at the boundary is revealed by peaks in the dissipation curves. For the upper boundary (right panel in Fig. 19), the absence of any peaks suggests that the adopted resolution is adequate to resolve this boundary of these four models. With the exception of the eps33k model, all of the models steadily decrease in relative dissipation towards the boundary, followed by a steeper slope beyond the boundary, and eventually plateauing to a very small value within the stable upper region. For the eps33k model, however, the TKE is so great that the reduction in dissipation across the boundary region occurs at an almost constant slope. A possible explanation for such behaviour is that the TKE of this model is so great that it destroys the boundary layer altogether (see Fig. 9).

The lower boundary is much narrower than the upper boundary (see e.g. Table 2). This is also apparent in Fig. 19 from the appearance of spurious peaks in the relative dissipation curves in all models except eps33k. The increase in driving luminosity broadens the lower boundary, resulting in an increase in the effective resolution for the higher energy models. This broadening reduces the amplitude of the dissipation peaks. In the eps10k model the dissipation peak is small and in the eps33k it is almost non-existent. As expected, a lower spatial resolution is needed to resolve a broader boundary. We see that an increase in the luminosity at a fixed resolution explored in this study has the same positive effect on resolving the boundary as an increase in resolution at a fixed luminosity (as done in C17).

Appendix E Equilibrium entrainment regime

Turbulent mixing is considered to occur within the equilibrium entrainment regime if the time scale for the boundary migration due to entrainment, $\tau_{b}$, is comparable to or larger than the convective turnover time scale, $\tau_{c}$ (Fedorovich, Conzemius and Mironov, 2004; Garcia and Mellado, 2014). These time scales are defined as

[TABLE]

[TABLE]

where $\delta r$ is the boundary width (see Table 2) and $v_{e}$ is the entrainment velocity (see Fig. 13). The time scales for each boundary and the ratio $\tau_{b}$/$\tau_{c}$ are given for the eps1k, eps3k, eps10k and eps33k models in Table 4. This ratio is moderately larger than one for all of these boundaries suggesting that turbulent entrainment occurs within the equilibrium entrainment regime.

Appendix F Bulk Richardson number

The bulk Richardson number is a useful diagnostic for the susceptibility of a boundary region to the entrainment of material through turbulent motions (more broadly categorised as convective boundary mixing). We refer to this susceptibility as the stiffness of the boundary, i.e. a stiffer boundary is less susceptible to convective boundary mixing. The bulk Richardson number is defined as the ratio of the specific stabilisation potential (analogous to the work done against convective motions by the boundary) to the specific TKE within the convective region. It is written as

[TABLE]

where $\Delta B$ is the buoyancy jump across the boundary. Based on the results of Meakin and Arnett (2007), we take the integral length scale, $\ell$, to be half of the local pressure scale height at the boundary. The buoyancy jump is estimated by integrating the square of the buoyancy (Brunt-Väisälä) frequency, $N$, over a suitable distance ($\Delta r$) either side of the boundary centre, $r_{c}$,

[TABLE]

The integration distance $\Delta r$ is not well defined theoretically but it should be large enough to capture the dynamics of the boundary region and the distance over which fluid elements are decelerated. In these simulations we take the integration distance, $\Delta r$, to be a quarter of the local pressure scale height at the boundary.

The buoyancy frequency is the frequency with which a perturbed fluid element will oscillate at if it is surrounded by a stably stratified medium. This frequency is imaginary for a convectively unstable fluid element and is defined as:

[TABLE]

Bulk Richardson numbers less than 10 are associated with relatively soft convective boundaries, whereas bulk Richardson numbers greater than 100 are associated with relatively stiff convective boundaries.