3D simulations of a neon burning convective shell in a massive star

ABSTRACT

The treatment of convection remains a major weakness in the modelling of stellar evolution with one-dimensional (1D) codes. The ever-increasing computing power makes now possible to simulate in three-dimensional (3D) part of a star for a fraction of its life, allowing us to study the full complexity of convective zones with hydrodynamics codes. Here, we performed state-of-the-art hydrodynamics simulations of turbulence in a neon-burning convective zone, during the late stage of the life of a massive star. We produced a set of simulations varying the resolution of the computing domain (from 128³ to 1024³ cells) and the efficiency of the nuclear reactions (by boosting the energy generation rate from nominal to a factor of 1000). We analysed our results by the mean of Fourier transform of the velocity field, and mean-field decomposition of the various transport equations. Our results are in line with previous studies, showing that the behaviour of the bulk of the convective zone is already well captured at a relatively low resolution (256³), while the details of the convective boundaries require higher resolutions. The different boosting factors used show how various quantities (velocity, buoyancy, abundances, and abundance variances) depend on the energy generation rate. We found that for low boosting factors, convective zones are well mixed, validating the approach usually used in 1D stellar evolution codes. However, when nuclear burning and turbulent transport occur on the same time-scale, a more sophisticated treatment would be needed. This is typically the case when shell mergers occur.

convection, hydrodynamics, nuclear reactions, nucleosynthesis, abundances, turbulence, stars: evolution, stars: interiors

1 INTRODUCTION

To study the physics and evolution of stars, the most efficient and complete software tools available are one-dimensional (1D) stellar evolution models (Heger, Langer & Woosley 2000; Paxton et al. 2011; Ekström et al. 2012). Not only can these models be employed to study the physics of stars, but they also represent a key aspect for interpreting stellar observations. Without an up-to-date and reliable grid of 1D stellar models, it would not be possible to obtain important information e.g. from asteroseismic measurements (e.g. Aerts et al. 2003), or for isochrone fitting and age determination (e.g. Jørgensen & Lindegren 2005; Bossini et al. 2019). Despite the great progress made in recent years to improve stellar evolutionary models, several uncertainties still affect the outcome of the models, undermining the accuracy of the predictions. These uncertainties arise from the multiphysical and multidimensional processes that occur in stars, often included only through simplified prescriptions. Among the most concerning uncertainties in stellar modelling are determining the extent of convective regions and the mixing that occurs near the boundaries (convective boundary mixing, CBM).

For decades, 1D models have implemented the so-called mixing length theory (MLT; Böhm-Vitense 1958), a purely one-dimensional and local treatment of convection that, while simple and easy to implement, fails to appreciate the multidimensionality of the problem. This reflects in the well-known problem of having to add ad hoc mixing beyond the convective boundary, often referred to as ‘overshoot’. Over the years, many prescriptions have been suggested and improved to account for an extension of the convective region (e.g. Canuto & Mazzitelli 1991; Zahn 1991; Freytag, Ludwig & Steffen 1996; Herwig 2000; Gabriel et al. 2014), but their validation and calibration has been difficult due to the impossibility of directly measuring the size of convective regions in observed stars. Indirect information can be deduced from asteroseismology (Pedersen et al. 2021), the Hertzsprung–Russell diagram (Castro et al. 2014), and eclipsing binaries (Claret & Torres 2016), but these methods are still heavily based on 1D models for interpreting the observations.

In the last two decades, multidimensional hydrodynamic simulations of stellar interiors have been used to constrain and parametrize CBM (e.g. Freytag et al. 1996; Meakin & Arnett 2007). This approach is sometimes called ‘321D-guided’, referring to the fact that 1D models are improved with the results obtained from multidimensional models. In fact, multidimensional simulations are started from initial conditions assumed from the 1D models they are trying to validate, therefore they are subjected to the uncertainties of the 1D models. In addition, three-dimensional (3D) models are not able to reproduce the long time-scales typical of the stellar evolution, due to the excessive computing resources required, therefore their results need to be generalized into more extensive prescriptions.

Running multidimensional hydrodynamic simulations of stars is very challenging due to the great amount of computing resources required. Consequently, an important limiting factor in the simulations is the Mach number or the convective velocity of the fluid. The late burning phases of massive stars are normally characterized by large convective velocities, which make it easier and cheaper to perform hydrodynamic simulations of these environments using time explicit methods. They normally include oxygen- and silicon-burning phases (Meakin & Arnett 2007; Arnett, Meakin & Young 2009; Viallet et al. 2013; Couch et al. 2015; Müller et al. 2016; Jones et al. 2017; Yoshida et al. 2019), and occasionally also neon burning (Rizzuti et al. 2022, 2023). Recently, these convective phases have also been explored with rotation (Yoshida et al. 2021; Fields 2022; McNeill & Müller 2022), magnetic fields (Varma & Müller 2021; Leidi et al. 2023), and both rotation and magnetic fields (Varma & Müller 2023). Earlier phases, such as the main-sequence burning, helium- and carbon-burning, are characterized by slow convective velocities, which require large boosting in luminosity to make the hydrodynamic computation affordable (Gilet et al. 2013; Woodward, Herwig & Lin 2015; Cristini et al. 2017; Horst et al. 2021; Herwig et al. 2023; Andrassy et al. 2024) or in two dimensions (Baraffe et al. 2023). However, it is still not clear whether such boosting also introduces additional effects on the physics of the star.

In this work, we present a large grid of hydrodynamic simulations reproducing a neon-burning shell in a 15 M_⊙ star; this grid explores a wide range of resolution and luminosity. Some of these simulations have been used in Rizzuti et al. (2022) to study entrainment. Here, we perform a detailed analysis of the dynamics and nucleosynthesis and their interplay. The presence of a nominal-luminosity run allows us to validate the results normally obtained through luminosity extrapolation; having simulations with different boosting factors allows us to study how different quantities scale with the luminosity. We also perform a detailed mean-field decomposition of some key equations, to study the statistical properties of the fluid. Finally, an explicit network of isotopes has been included for realistically reproducing nuclear burning, allowing us to study the evolution of the abundances and their variance. The results we present here contribute to the general understanding of stellar evolution through our detailed analysis and critical comparison to the initial conditions.

We organize the paper as follows: In Section 2, we present the initial conditions and a general overview of the hydrodynamic simulations. In Section 3, we present the results concerning the dynamics and nucleosynthesis of the convective flow and we present the detailed mean-field analysis performed with the Reynolds-averaging framework. In Section 4, we analyse the time evolution and transport of the chemical composition. Finally, we discuss our results and draw conclusions in Section 5.

2 INITIAL CONDITIONS AND OVERVIEW OF 3D HYDRODYNAMIC SIMULATIONS

In this paper, we simulate a neon-burning shell of a 15 M_⊙ star. The initial conditions for our 3D models are mapped from a 1D genec stellar evolution model run at solar metallicity (Z = 0.014) and with the physical ingredients described in Ekström et al. (2012). The Schwarzschild criterion is used, and penetrative overshoot is included for core hydrogen and core helium-burning phases only, with an overshooting distance |$\ell _{\text{over}}=0.1\, H_\mathrm{ P}$|⁠. From carbon-burning onward, an α-chain network is used (see Hirschi, Meynet & Maeder 2004, for details). This input stellar model has been described in more detail in Rizzuti et al. (2022). The evolution of the structure of the model is presented in Fig. 1 (left) with a close-up of the neon-burning shell used as initial conditions for the 3D simulations presented in Fig. 1 (right). Compared to other convective episodes, neon-burning convective zones are relatively small and short-lived with typical spatial extents of the order of 10^8–9 cm and temporal extents of weeks to months. These time-scales are still too computationally expensive to be simulated completely in 3D (though see Rizzuti et al. 2023), so in this study we opt instead to simulate a statistically significant number of convective turnovers.

Figure 1.

Left: Structure evolution (also known as Kippenhahn) diagram of the 1D genec 15 M_⊙ stellar model used as input for this study. The top solid black line shows the time evolution of the total mass. The other solid black lines are radial contours (values of log₁₀(r) in cm are indicated on the curves). The Mach number of the flow inside convective zones (shaded areas) is colour coded. The extent of the computational domain covered by the 3D simulations corresponds to the red vertical bar (indicated by the red arrow). Right: Zoom-in on the neon shell convective region simulated in this study. The age_start-hydro corresponds to the time at which the 3D simulations start. Black solid lines are isomass contours (⁠|$M_\mathrm{ r}=0.9, \ 1.0,\ 1.1, \ 1.2\, \mathrm{ M}_\odot$|⁠). These isomass contours show that the shell studied does not undergo any significant contraction or expansion during the Ne-burning phase. The 3D simulations presented in this paper cover a physical time of an hour or less, thus not even covering the width of the red vertical line in this figure. These diagrams are taken from Rizzuti et al. (2022).

Using the initial conditions described above, we produced a series of 3D hydrodynamic simulations using the prompi code (Meakin & Arnett 2007; Cristini et al. 2017). We used a plane-parallel geometry and our computing domain is a cube of side equal to |$0.65\times 10^8\, \text{cm}$| encompassing the convective shell, which occupies half of the domain, as well as stable (radiative) layers both below and above (see e.g. Fig. 2). The gravity is fixed according to the initial stellar model, and nuclear burning is followed with a minimal nuclear network of four species (¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si), α-particles being considered to be at equilibrium (Arnett 1974). As already done in our previous works (Cristini et al. 2017; Rizzuti et al. 2022, 2023), for a subset of our simulations, we multiplied the reaction rates and the neutrino generation rate by a so-called ‘boosting’ factor taking the following values: 1 (nominal case), 10, 100, and 1000. This energy generation boosting provides two main advantages. First, the increased energy generation speeds up the nuclear burning and convection (convective velocities scale with the cubic root of the energy generation rate, see below). Higher convective velocities mean that sufficient statistics (i.e. a reasonable number of convective turnovers) can be obtained with a smaller computing budget, as time-steps are limited by the sound speed in prompi (time explicit integration). Second, and most important, having several simulations with different energy generation rates enable us to study the dependence of the results on the strength of turbulence. It is worth noting that this study includes a nominal case, i.e. exactly the same nuclear energy generation and neutrino losses as in the 1D input model. This enables a direct comparison to the 1D input model without the need for extrapolation.

Figure 2.

Vertical cross-sections of the velocity magnitude (values in colour scale) for the Ex1 model, taken at 1000, 1700, 2400, and 3001 s (last time-step) through the simulation. The progression shows an increase in the velocity magnitude of fluid elements. It is possible to see in the upper panels the formation of the convective cell, and in the lower ones its growth into the upper stable region as expected from turbulent entrainment. Gravity waves are also visible in the stable region.

We summarize the models presented in this study and their most important properties in Table 1. One can estimate an effective Reynolds number for our simulations as Re ∼(n/2)^4/3, where n is the resolution, considering that the radial extent of the convective zone is half of the domain and thus covers n/2 cells (see Arnett et al. 2019). This is a reasonable assumption, which is confirmed by the cross-sections of Figs 2 and 6 during the quasi-steady state. In this way, we obtain effective Reynolds numbers of 256, 645, 1625, and 4096 for the resolutions n = 128, 256, 512, and 1024, respectively. For large Reynolds numbers, we can assume the regime of turbulent cascade (Biermann 1932; Kolmogorov 1941), where kinetic energy cascade is driven from the large scales. If we use Re ≥ 1000 as the condition for the turbulent regime, we see that the highest resolutions 512³ (and 1024³) reach the numerically turbulent regime. We thus focus our analysis on simulations at the 512³ resolution throughout the paper and shorten their code name to include only their boosting factor for convenience (Ex1, Ex10, Ex100, and Ex1000). This being said, we analyse simulations with different resolutions to test the dependence of the results on the grid size in Section 3. Finally, we confirm from Table 1 that increasing the boosting factor by a factor of 10 increases the convective velocity roughly by a factor of 10^1/3 ≃ 2.15 (v_rms ∼ ϵ^1/3, where ϵ is the energy generation rate; Arnett et al. 2018). This in turn decreases the turnover, τ_c, as expected. The general thermodynamic properties of the simulated region of the star remain mostly unchanged by changing the boosting factor, thus the sound speed also remains the same. The Mach number scales therefore with v_rms, and is small in any case: about 2 × 10⁻³ for the Ex1 case to about 20 × 10⁻³ for the Ex1000 case. We also confirm that simulations at different resolutions with the same boosting factor have approximately the same v_rms, τ_c, and Ma.

Table 1.

Summary of the models presented in this study, listed under names indicating their boosting factor and resolution. The properties are: resolution N_xyz, boosting factor of the nuclear energy generation rate ϵ, physical time simulated τ_sim (s), global rms convective velocity v_rms (cm s⁻¹), convective turnover time τ_c (s), time spent in quasi-steady state τ_q (s), number of convective turnovers simulated in quasi-steady state phase n_c, and computational cost in CPU core hours.

	Nxyz	ϵ	τsim	vrms	Ma	τc	τq	nc	Cost
			(s)	(⁠|$10^6 \,\rm{cm\,s}^{-1}$|⁠)	(10−3)	(s)	(s)		(106 h)
Ex10_128	1283	10	1500	1.58	4.53	45	1250	27	0.04
Ex10_256	2563	10	1832	1.70	4.44	44	1600	36	0.36
Ex1a	5123	1	3037	0.70	1.93	100	700	7	11.4
Ex10	5123	10	1004	1.49	3.96	46	700	15	4.66
Ex100	5123	100	291	3.72	9.95	23	150	6	1.15
Ex1000	5123	1000	73b	8.00	23.2	13	10	1	0.28
Ex10_1024	10243	10	310c	1.49	3.94	47	310	6	48.2

Notes. ^aModels with resolution 512³ are indicated by their boosting factor only, since they are the most studied in this work.

^bTime of the entire simulation, although the upper domain is reached at about 30 s, as it can be seen in Fig. 3.

^cModel Ex10_1024 was restarted from the Ex10 simulation at 500 s.

3 GENERAL PROPERTIES OF THE CONVECTIVE FLOW

In this section, we present our analysis of the convective flow with a particular emphasis on the velocity field and the turbulent kinetic energy (TKE).

3.1 Flow velocity and specific kinetic energy

Fig. 2 shows vertical cross-sections of model Ex1 taken at key phases of the simulation, with the velocity magnitude represented in colour scale. These cross-sections allow us to see the time evolution of the velocity field in this simulation, which is representative of all our models. At the beginning of the simulation (t = 1000 and 1700 s), the nuclear burning at the bottom of the neon shell drives turbulent motions (initially plumes and later on eddies), which gradually fill the region that was convective in the 1D input stellar model. After an initial transient (the duration of which strongly depends on the boosting factor), convection is fully developed and the simulation enters a quasi-steady state (at around |$2400\, \text{s}$| for the Ex1 model). As time proceeds, the convective region grows and entrains material from the stable regions, as can be seen at the top boundary by comparing the cross-sections at t = 2400 and 3001 s.

Another way to follow the time evolution of the velocity in our simulations is via the specific total kinetic energy |$(v_x^2+v_y^2+v_z^2)/2$| (see Fig. 3). The initial and sudden rise in kinetic energy characterizes the initial transient. Its length depends on the boosting factor, being considerably longer for the non-boosted model (Ex1). After the transient, the simulations enter a quasi-steady state, where the average kinetic energy remains constant or slowly increases. This phase is also characterized by periodic pulses in kinetic energy of approximately the same time-scale as the convective turnover time. These pulses are related to the time delay between the formation and rise of large-scale eddies in the convective zone and their dissipation (Meakin & Arnett 2007; Arnett & Meakin 2011; Viallet et al. 2013; Arnett et al. 2015).

Figure 3.

Time evolution of the specific kinetic energy integrated over the computational domain, for the four models Ex1, Ex10, Ex100, and Ex1000. The trends cover the entire simulated time range. After an initial transient, whose duration depends on the boosting, the simulations enter a quasi-steady state (starting time indicated by the vertical dashed lines), which lasts until the upper boundary of the domain is reached (Ex100 and Ex1000) or the simulation is otherwise terminated (Ex1 and Ex10). As expected, models with larger boosting factors reach higher kinetic energies.

The roughly constant or slow rise of convective velocities is confirmed by comparing the radial profiles of the root-mean-square velocity between the start (red solid line in Fig. 4 for Ex1 model) and the end of the quasi-steady state phase (blue solid line). In addition to a mild increase in magnitude, the outward shifting of the upper boundary of the convective zone expected from entrainment is visible. This figure also confirms findings from previous 3D simulations (see e.g. Cristini et al. 2017; Jones et al. 2017; Mocák et al. 2018). The radial velocity component (dotted line) peaks in the centre of the convective zone, while the horizontal component peaks at the boundaries. This reflects the convective motions and the U-turning of fluid elements at the boundaries. The non-negligible velocities outside the convective zone are produced by gravity waves. Furthermore, using our nominal case simulation (Ex1), we can compare convective velocities directly between 3D and 1D models. We see that convective velocities in our 3D simulations are about twice as large as the velocity predicted by MLT in the 1D genec model (around 3.6 × 10⁵ cm s⁻¹, black solid line) consistent with the results of Jones et al. (2017).

Figure 4.

Radial profiles of different velocity components. In black, the MLT velocity v_MLT of the 1D stellar model was used as initial conditions for all simulations. In red, the root-mean-square velocity v_start at the beginning of the quasi-steady state in Ex1 simulation, averaged over one convective turnover. In blue, different components of the convective velocity at the end of Ex1 simulation, averaged over one convective turnover: root-mean-square v_end (solid), radial v_r (dotted), and horizontal v_hor (dashed) velocity components.

The boosting of the nuclear energy generation rate has a strong impact on the evolution of our simulations. This is clearly visible in Fig. 5, where we compare the velocity fields between our four 512³ simulations with different boosting factors (see Table 1). Since the time-scale of the evolution of the models is affected by the boosting, as can be seen in Fig. 3, we choose here to compare the different simulations when their upper convective boundary has reached a radial location of approximately 3.89 × 10⁸ cm, which is the maximum extension of the convective zone in the non-boosted model. This corresponds to a time of 3001, 400, 60, and 16 s for the Ex1, Ex10, Ex100, and Ex1000 models, respectively (during the quasi-steady state phase of these simulations). Looking at the highest velocity in the snapshot shown in Fig. 5 (coloured in red), we can see that larger boosting factors produced a higher kinetic energy of the fluid, which is confirmed in Fig. 3.

Figure 5.

Same as Fig. 2, but for the four models Ex1, Ex10, Ex100, and Ex1000 taken at 3001, 400, 60, and 16 s, respectively. Time-steps were chosen so that the upper convective boundary is located at 3.89 × 10⁸ cm in all simulations, for comparison. It is clear that one major effect of boosting the nuclear energy generation rate is an increase in the velocity magnitude of fluid elements, therefore in kinetic energy, as confirmed by Fig. 3.

In order to test whether there is a dependence of the flow velocity on resolution, we present in Fig. 6 the velocity fields of four simulations with the same boosting factor and initial conditions, but different resolutions. It is clear from the figure that the small-scale features of convection depend on the mesh size we choose for our simulations. In particular, according to the Implicit Large Eddy Simulation (ILES) paradigm, the grid scale sets the limits for the numerical dissipation of kinetic energy, which mimics the effects of viscosity. For this reason, if we increase the resolution of our models, the dissipation scale decreases, allowing the simulations to produce eddies on a smaller scale, which are closer to the real case scenario, as visible in Fig. 6 and discussed further in Section 3.2.

Figure 6.

Same as Fig. 2, but for the four different resolutions 128³, 256³, 512³, and 1024³ with boosting factor 10 (see Table 1), taken at 800 s from the beginning of each simulation. Since a higher resolution is linked to a smaller dissipation range, it is expected that the simulations predict eddies on a smaller scale when the resolution is increased.

At large scale, however, the structures look similar at different resolutions so we do not expect the bulk properties of the convective region to depend on resolution. This is confirmed in Fig. 7 showing the time evolution of the specific total kinetic energy for the four simulations with different resolution and the same boosting factor (Ex10_128, Ex10_256, Ex10, and Ex10_1024, see Table 1). As expected, the simulations have a very similar evolution, in agreement with the fact that global properties (like v_rms and τ_c estimated in Table 1) are comparable for all simulations with the same boosting factor and thus do not depend on the resolution.

Figure 7.

Same as Fig. 3, but for the four models with different resolutions 128³, 256³, 512³, and 1024³ with the same boosting factor 10. The simulations share a very similar evolution. The vertical dashed line is the beginning of the quasi-steady state.

High resolution is nevertheless needed to better resolve convective boundaries as discussed in Section 3.3.

3.2 Turbulent kinetic energy spectra

A different approach to study the kinetic energy of the simulations is to compute the turbulent kinetic energy (TKE) spectra. In order to do so, we performed 2D fast Fourier transforms of different velocity components on horizontal planes at constant height, within the convective zone. We also normalized the spectra by dividing them by k^−5/3, which is the power-law scaling for the inertial range (Kolmogorov 1941). In this way, the regions where the spectra have a horizontal slope correspond to the inertial range (see e.g. Cristini et al. 2017, and references therein).

We show in Fig. 8 a comparison between the spectra of the radial velocity squared, for different resolutions (left panel) and different boosting factors (right panel). The spectra are averaged through the entire quasi-steady state phase of each simulation (see Fig. 3). As expected, increasing the resolution has the effect of extending the inertial range plateau toward higher k, because dissipation occurs at smaller scales r in ILES. The vertical dashed lines indicate the grid size for each resolution. We can see that the dissipation range is 5–10 times larger than the grid size. On the other hand, if we increase the boosting factor while keeping the resolution fixed, we witness a rise in the velocity magnitude without changes in the length of the plateau, as visible in the right panel of Fig. 8.

Figure 8.

Two-dimensional Fourier transforms of the radial velocity multiplied by k^5/3, for the four different resolutions 128³, 256³, 512³, and 1024³ with a boosting factor 10 (left panel), and for the four different boosting factors of the 512³ resolution models (right panel), taken during the quasi-steady state phase, at the centre of the convective zone. The horizontal axes show both the Fourier space |$k=\sqrt{k_y^2+k_z^2}$| (bottom axis) and the real space r (top axis). The vertical dashed lines show in the left plot the grid size for each resolution, and in the right plot they show 32, 16, and 8 times the grid size for the 512³ models.

The spectra presented in Fig. 8 are taken in the bulk of the convective region and during the quasi-steady-state phase. It is interesting to find out if the spectra vary with location and time. In Fig. 9, we study the time evolution of the velocity spectra for the Ex10 model from the beginning of the simulation to the quasi-steady state, with time-averaging windows of 50 s, which is approximately the convective turnover time for this model (see Table 1). The spectra are taken near the bottom of the convective zone (at a radius of 3.58 × 10⁸ cm), and we present results for the radial velocity (left panel), the horizontal velocity (central panel), and the total root-mean-square velocity (right panel). As can be seen from the plots, at the beginning of the simulation (during the initial transient), the spectra have a peak at high k, i.e. at small scales. This happens because convection is at a very early stage, and eddies on the smallest scales are dominant. As time passes, the velocity magnitude increases and the peaks are shifted toward smaller k as the turbulent flow fills the entire region that was convective in the 1D input stellar model. After the initial transient (lasting about 300 s), the spectra do not vary significantly and they assume the more familiar shape of Fig. 8, which is characteristic of homogeneous and isotropic turbulence.

Figure 9.

Same as Fig. 8, but for the Fourier transforms of the radial velocity (left panel), horizontal velocity (central panel), and total root-mean-square (rms) velocity (right panel), taken at different times through the Ex10 simulation, with an averaging window of 50 s. It can be clearly seen the progression toward the quasi-steady state.

In a similar way, we compare in Fig. 10 the spectra of the velocity and its components at different radial locations inside the convective zone, from 3.58 to 3.88 × 10⁸ cm, during the quasi-steady state phase (500–550 s, Ex10). The spectra at different heights are very similar, which means that convection remains turbulent and generally isotropic throughout the convective zone. The main exception is the radial velocity at the lowest k values (largest scales) near convective boundaries. Indeed, both the uppermost and lowermost spectra have a lower radial velocity magnitude around k = 10 compared to the bulk of the convective region. This is easily explained if we consider that the convective boundaries are limiting the velocity in the radial direction, while the horizontal components of velocity are not affected by this restriction. Instead, they are stronger near the convective boundaries and present a peak around k = 11, which is less important for central regions.

Figure 10.

Same as Fig. 9, but for the Fourier transforms taken during the quasi-steady-state phase (500 s, Ex10) at different radial locations inside the convective zone.

3.3 Mean-field (RA-ILES) analysis of the turbulent kinetic energy

In order to gain a better understanding of the key processes taking place in our simulations, we performed a mean-field analysis, called Reynolds-Averaged analysis of Implicit Large Eddy Simulations (RA-ILES; see Meakin & Arnett 2007; Viallet et al. 2013; Cristini et al. 2017; Arnett et al. 2018). This allows us to disentangle the contributions and interplay between nuclear reactions and turbulence and to determine the dominant processes at play. It also provides insightful quantitative information. The basics of this analysis rely on a time- and space averaging of all the quantities. The time averaging is done on a time window T, which is long enough to be statistically meaningful, but short enough so that the main properties of the fluid do not change significantly. Usually, we perform this averaging over a time corresponding to two convective turnover time-scales. The spatial averaging is performed over a given volume V (in our simulation, it is usually the ‘volume’, or rather surface of a horizontal slice of cells). One can therefore define the ‘Reynolds average’ of a quantity q as

$$\begin{eqnarray} \bar{q}(x,t) = \frac{1}{T\Delta S}\int ^{t+\frac{T}{2}}_{t-\frac{T}{2}}\int _{\Delta S} q(t^{\prime },x,y,z) \,\text{d}S \, \text{d}t^{\prime }, \end{eqnarray}$$

(1)

dS being an infinitesimal surface centred on the point (x, y, z). The quantity q can then be decomposed as

$$\begin{eqnarray} q = \bar{q} + q^{\prime }. \end{eqnarray}$$

(2)

Another useful type of averaging is a density-weighted average, also called ‘Favre average’,

$$\begin{eqnarray} \widetilde {q} = \frac{\overline{\rho q}}{\overline{\rho }}. \end{eqnarray}$$

(3)

Again, we can decompose any quantity q as

$$\begin{eqnarray} q = \widetilde {q}+ q^{\prime \prime }. \end{eqnarray}$$

(4)

Note that most of the time, q′ ≠ q″.

With these notations, we can perform averages of the Euler equations governing the fluid we are simulating. Without providing details (which can be found in¹ Viallet et al. 2013; Mocák et al. 2014, 2018; Arnett et al. 2015), this leads to the following results for the kinetic energy equation:

$$\begin{eqnarray} \overline{\rho }\, \widetilde {{D}}_{t}\widetilde {\epsilon }_\text{k} = -\nabla _x\overline{\rho }\widetilde {u_x^{\prime \prime }\epsilon _\text{k}}-\nabla _x f_\text{P} + W_\mathrm{ P} + W_\text{B}. \end{eqnarray}$$

(5)

|$\overline{\rho }\, \widetilde {{D}}_{t}$| is an operator similar to the Lagrangian derivative in the RA-ILES framework: |$\overline{\rho }\widetilde {{D}}_{t} \widetilde {q} = \partial _{t}\left(\overline{\rho }\widetilde {q}\right) + \nabla _x \left(\overline{\rho }\widetilde {u}_x\widetilde {q}\right)$|⁠, with ∇_x the x component of the divergence operator,² and |$f_\text{P} = \overline{P^\prime u^{\prime }_x}$| is the acoustic flux (i.e. the flux of pressure variations). |$W_\text{P} = \overline{P^{\prime }\partial _x u_x^{\prime }}$| represents the turbulent pressure dilatation, and |$W_\text{B} = -\overline{\rho }\overline{u_x^{\prime \prime }}\widetilde {g}$| is the buoyancy work.

In this section, we go into the details of the RA-ILES decomposition of the TKE by discussing the importance and behaviour of each term in equation (5) for different resolutions (128³ to 1024³ for the Ex10 boosting factor, see Fig. 11) and boosting factors (nominal luminosity to 10³ times nominal luminosity at a resolution of 512³, see Fig. 12).

Figure 11.

Left column: Mean kinetic energy equation terms as a function of radius for various resolutions. From top to bottom: 128³, 256³, 512³, and 1024³. The time averaging is made over two turnover time-scales centred on |$t=650\,\text{s}$| for each resolution. The meaning of each curve is described in the caption in the first row. The grey area shows the convective region, and the red area the region where we applied a damping. Right column: Radial integration of each term of the mean kinetic equation for the same resolutions as in the left column.

Figure 12.

Same as Fig. 11, but for various boosting factors instead of various resolutions. From top to bottom: nominal nuclear energy generation rate, energy generation rate boosted by a factor of 10, 100, and 1000. The time averaging is made over two turnover time-scales centred on a time where the convective zone has about the same radial extent for each case (see text). The resolution of the simulations shown here is 512³.

3.3.1 Time evolution

|$\overline{\rho }\widetilde {D_t}\widetilde {\epsilon _{\rm k}}$| represents the Lagrangian time derivatives of the kinetic energy. A negligible time derivative within the chosen time average implies that the convection is in a statistically steady state. This can seen to be true in all resolutions (Fig. 11) and most boosting factors. The only exception is the Ex1000 model shown in Fig. 12. In this case, the high boosting factor leads to rapid growth of the convective region, causing it to interact with the domain boundary before a steady state can be achieved. As such, we will not analyse the Reynolds-averaged quantities of this model further.

3.3.2 Transport terms

The transport of kinetic energy throughout the convective region is defined by the TKE flux, |$\nabla _x(\overline{\rho } \widetilde {u^{\prime \prime }_{x}\epsilon ^{\prime \prime }_{\rm k})}$|⁠, and the acoustic flux, ∇_xf_P, where |$f_\mathrm{ P} = \overline{P^\prime u^{\prime }_x}$|⁠. We see that the acoustic flux often opposes the TKE flux within the convection zone, except at the convective boundaries, where the acoustic flux spikes where it launches gravity waves. The general behaviour observed in the Ex10 case in Fig. 11 can also be seen when the boosting factor varies (Fig. 12). The general pattern is preserved, but its amplitude is increased for increasing boosting factors, since the higher velocity of the fluid makes the transport more efficient (note the change in the scale on the y-axis by a factor of 1000 from the top to the bottom panel).

3.3.3 Source terms

Turbulence in the convective region is driven by the work due to turbulent pressure fluctuations, W_P, and the buoyancy work due to density fluctuations, W_B. As in previous work by Viallet et al. (2013) and Cristini et al. (2017), W_P appears to be negligible in shell convection in deep interiors. |${W}_\mathrm{ B} \gt 0$| is generally seen in the convective zone as expected, since it is the main driving term. Regions near the convective boundary have |${W}_\mathrm{ B} \lt 0$|⁠, meaning that the flow decelerates in these regions. These regions are usually referred to as ‘overshooting’ regions in 1D stellar modelling but we prefer the more general term of CBM regions. We notice a minor systematic decrease in the integrated buoyancy work as resolution is increased. We suspect this effect is caused by the lower resolution models having a larger CBM extent. This leads to a higher rate of entrainment of new material to burn into the convective Ne shell, which in turn, slightly increases the rate of energy generation and hence the buoyancy work done. The difference, however, is quite minor even over the relatively long time-scales that we simulate. Looking at the different boosting factors (Fig. 12), we see that W_B scales linearly with them: the work done by buoyancy force is about 1000 times higher in the Ex1000 simulation compared to the nominal case Ex1. This is due to the fact that nuclear processes drive convection by heating up the plasma at the bottom of the shell.

3.3.4 Dissipation

Due to the finite size of our grid, our code is not able to perfectly reproduce the behaviour of the fluid at spatial sizes that are smaller than the grid size and our simulations do not include explicit viscosity (ILES). The numerical dissipation (taking place at the grid scale) is nevertheless tracked by the |$\mathcal {N}_{\epsilon _{\rm k}}$| term, which is the difference between the left-hand side and right-hand side terms of equation (5). Similar to the results found in Cristini et al. (2017) for a convective carbon-burning shell, we find that the work due to buoyancy, W_B, is closely balanced by numerical dissipation |$\mathcal {N}_{\epsilon _{\rm k}}$| as expected from Kolmogorov’s turbulence theory as numerical dissipation in ILES replaces the dissipation due to physical viscosity. One key feature that depends on resolution is the numerical dissipation peak at the lower boundary. This confirms that high resolution is needed to fully resolve the lower boundary of convective region. Another feature that affects low-resolution simulation is the leap-frog instability that creates zigzags in the W_B and thus also the |$\mathcal {N}_{\epsilon _{\rm k}}$| profiles. This analysis shows that a resolution of 256 or higher is best to resolve the key processes taking place in the bulk of convective regions and not too steep boundaries like the top boundaries in this study. A resolution of 1024 or higher is needed to perfectly resolve sharp boundaries like the lower ones but measurements of entrainment do not necessarily require the lower boundary to be perfectly resolved (Rizzuti et al. 2022, 2023).

4 INTERPLAY BETWEEN TURBULENCE AND NUCLEAR PROCESSES

Our simulations follow explicitly the evolution of the four key nuclides for neon burning: ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si, which enables us to study the interplay between turbulence and nuclear processes in great details. These nuclides are linked by a small tailored nuclear reaction network including the following reactions: ²⁰Ne(γ, α)¹⁶O, ¹⁶O(α, γ)²⁰Ne, ²⁰Ne(α, γ)²⁴Mg, and ²⁴Mg(α, γ)²⁸Si (α particles are considered to be at nuclear equilibrium, which is a reasonable assumption for the neon shell studied here).

We can see the effects of the above-mentioned neon-burning reactions in the vertical cross-section snapshots presented in Fig. 13, where we show in colour scale the mass fraction of these four species after 1500 s from the start of the Ex1 model (note that the simulation has not reached the quasi-steady state by that time and that the turbulent region has not reached the upper boundary from the initial 1D model). In particular, the reference colour (white) indicates the average mass fraction inside the neon shell at the beginning of the simulation. The colour represents the same range of deviations from the initial composition in the convective region in all panels. The results fully reflect what is expected for neon burning: neon is consumed inside the convective region, while oxygen, magnesium, and silicon are produced via the reactions listed above. These abundance plots also show that the species are not always completely homogeneously mixed and that they can be used as tracers of the turbulent motion of the fluid. Particularly interesting is the ²⁰Ne case: the bottom part of the convective region shows an underabundance with respect to the initial one. This is a sign of the nuclear burning taking place at the bottom boundary. On the contrary, the top part of the convective region shows an overabundance of neon. This is due to entrainment, which brings neon from the stable region above the convective zone (where abundance is higher) into the convective zone, where it is slowly mixed. Note that the strong overabundance in Ne between 3.5 and |$3.6 \times 10^8\, \text{cm}$| comes from the initial 1D structure (see Fig. 14, dashed line), which has not been erased by turbulent mixing at the time when the snapshot is taken here, but it is erased soon afterward.

Figure 13.

Vertical cross-sections of the mass fraction (values in colour scale) of the different isotopic nuclei ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si, taken at 1500 s in the Ex1 model. The reference value of the colour bars (white) is the average mass fraction of each nuclide in the convective zone at the beginning of the simulation. For this reason, the overabundance of ¹⁶O, ²⁴Mg, and ²⁸Si (towards the red) and the underabundance of ²⁰Ne (towards the blue) inside the convective zone reflect the nuclear reactions that are burning neon to produce oxygen and magnesium, and burning magnesium to produce silicon.

Figure 14.

Profiles of the abundances of the four chemical elements followed in our set of simulations at the beginning (dashed lines) and at the end (solid lines) of the simulated time, and for Ex1 (top-left), Ex10 (top-right), Ex100 (bottom-left), and Ex1000 run (bottom-right).

4.1 Evolution of the chemical composition

An overview of the time evolution of the chemical composition in our simulations is presented in Fig. 14, which shows the radial profiles of the four chemical elements followed in our simulations (i.e. the horizontal average of the abundance) at the start of the simulation (dashed lines, the same in each case), and at the end (solid lines), for the four boosting factors. The convective region corresponds to the plateaus in the middle of the computational domain. A few general observations can already be made from this figure. The effects of entrainment are clearly visible, particularly at the top boundary: the boundary of the convective zone has moved outwards, and the steep initial profiles have been smoothed by CBM. Interestingly, in the current set of simulations, entrainment is more efficient at bringing ²⁰Ne inside the convective zone than nuclear burning at destroying it, leading to a slight increase in ²⁰Ne abundance (except for the Ex1000 case). Qualitatively, the Ex1 and Ex10 cases are very similar, except that, as expected, entrainment is more efficient with the higher boosting factor. This is exacerbated for the Ex100 and Ex1000 cases where the convective boundary has moved up to the top of the computing domain over the course of the simulation.

4.2 RA-ILES analysis of the composition transport equation

As was done for the TKE in the previous section, our RA-ILES can also be applied to the chemical composition transport equation, resulting in the following equation:

$$\begin{eqnarray} \overline{\rho }\widetilde {{D}}_{t}\widetilde {X}_i = -\nabla _x f_i + \overline{\rho }\widetilde {\dot{X}}^{\text{nuc}}_i. \end{eqnarray}$$

(6)

X_i is the mass fraction of the chemical element i, |$f_i = \overline{\rho }\widetilde {X_i^{\prime \prime }u_x^{\prime \prime }}$| is the turbulent flux of the element i, and |$\dot{X}_i$| is the rate of creation/annihilation of the element i due to nuclear reactions. The right-hand side of this equation includes the two sources of the modification of the mass fraction at a given location inside the computing domain: turbulent motion transport and nuclear processes.

Before discussing the composition transport equation, it is useful to understand the turbulent composition flux, f_i, plotted in Fig. 15, for ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si. These fluxes represent the rate at which each species is transported. Negative values indicate that the fluxes are oriented inward. We see that large quantities of neon are being entrained from the upper boundary, and transported towards the inner parts of the convective burning zone, where the neon is burnt. The other species, on the other hand, are being transported outward (to the right of the convective zone). The behaviour at the inner boundary is more complex, due to the interplay between entrainment through the boundary, and the burning occurring slightly above. ¹⁶O and ²⁴Mg are more abundant inside the convective zone than outside (see Fig. 14). There is thus a slightly negative flux near the bottom boundary for these elements, showing that they are partly transported outside the convective region. On the other hand, both elements are produced through the burning of ²⁰Ne. Convective motions transport the freshly produced elements through the convective zone, as indicated by the positive flux inside the bulk of the convective region. The case of ²⁸Si is different still: there is more of it below the convective zone than inside. It is thus transported through the bottom boundary by entrainment into the convective region, hence the positive flux at the boundary. Furthermore, this element is also synthesized by neon burning (by a double α-capture), and redistributed inside the whole convective region from bottom to top, hence the positive flux for ²⁸Si everywhere.

Figure 15.

Turbulent composition fluxes for the nuclides explicitly followed in our work. These fluxes are taken from the 512³ resolution simulation with a boosting factor of 10. The time averaging is done as in Fig. 11, and the coloured areas have the same meaning. The flux for each element is defined as |$f_i = \overline{\rho }\widetilde {X_i^{\prime \prime }u_x^{\prime \prime }}$| (cf. Mocák et al. 2018).

The RA-ILES composition transport equation profiles for neon are shown in Fig. 16. Each term of equation (6) is shown by a different curve in the plot and is described below.

Figure 16.

Left panel: RA-ILES decomposition of the mean ²⁰Ne abundance equation. The meaning of each curve is specified in the legend. The coloured regions have the same meaning as in Fig. 11, and the time average is performed over the same time window. The results come from the 512³ resolution simulation, with a boosting factor of 10. Right panel: Radial integration of each term shown on the left-hand side of the figure.

4.2.1 Time evolution

|$\bar{\rho }\widetilde {D_t}\widetilde {X}_{i}$| and |$\partial _t(\bar{\rho }\widetilde {X}_{i})$| represent the Lagrangian and Eulerian time derivatives of the composition. While the convective region is in a steady state, we see that significant mixing occurs at the convective boundaries, particularly at the upper boundary, which changes the size of the convection zone over time.

4.2.2 Transport terms

The mean composition flux term is given by |$\nabla _x(\bar{\rho }\widetilde {u_x}\widetilde {X}_{i})$|⁠, which we see is non-negligible at the upper convective boundary, implying that there is a change in the composition profile due to the overall growth of the convective zone. This is due to a mean flux of Ne into the convective zone from the upper boundary, as seen in Fig. 15. The divergence of the turbulent flux, |$\nabla _x(\overline{\rho }\widetilde {X_i^{\prime \prime }u_x^{\prime \prime }})$|⁠, shows how material mixed into the convective zone is transported by turbulent velocity variations. This term is positive where the abundance is decreased by the element flux, and negative where there is an accumulation of a chemical element due to transport. We see in Fig. 16 that the turbulent fluxes accumulate the neon abundance near the top boundary. The time rate of change of neon balances this term. At the bottom boundary, the destruction of neon by the nuclear burning is compensated by neon being transported here.

4.2.3 Nuclear burning

The rate at which compositions change due to nuclear burning, |$\bar{\rho }\widetilde {\dot{X}}_{i}$|⁠, shows how the rate of neon burning increases with depth, and hence density and temperature. As it also depends on the abundance of the fuel, this term peaks at the bottom of the convective zone.

4.2.4 Numerical residual

|$\mathcal {N}_{X_i}$| term highlights the loss of information of our code at the grid level. It is small throughout most of our simulated domain, but we see that it is important near the edge of the convective zone, where strong turbulence develops. Increasing the resolution of the simulation makes this term smaller (Arnett et al. 2018, and Fig. 11).

Integrated over the whole computational domain (right panel), we see that the global change of the ²⁰Ne abundance is, as expected, mostly due to the nuclear reactions, which are responsible for the decrease in the abundance of this chemical element.

4.3 RA-ILES analysis of the composition variance

The variance, σ, of the composition traces the variation of the abundance of an element with respect to its average value at a given radius inside the computational domain. In our RA-ILES framework, it is defined as |$\sigma _X = \widetilde {X^{\prime \prime }X^{\prime \prime }}$|⁠. In Fig. 17, we show the relative standard deviation of the composition, i.e. |${\sqrt{\sigma _X}}/{\overline{X}}$|⁠. In the bulk of the convective region, the deviation remains small (around 1 per cent), justifying the general 1D modelling approach used in stellar evolution. We note two locations however where this deviation becomes larger for all the chemical elements: both boundaries of the convective zone. The convective boundaries are not flat or regular, but distorted and turbulent, making important variations with respect to the average of the abundance of a given element. This behaviour is much more important at the top boundary, where the relative deviation can reach about 30 per cent for silicon.

Figure 17.

Standard deviation of a chemical element abundance, normalized by its abundance. The coloured regions have the same meaning as in Fig. 11, and the time average is performed over the same time window. The results come from the 512³ resolution simulation, with a boosting factor of 10.

We can investigate the origins of the variance using the RA-ILES chemical variance equation:

$$\begin{eqnarray} \overline{\rho }\widetilde {{D}}_t\sigma _i = -\nabla _x f_i^\sigma -2f_i\partial _x\widetilde {X}_i + 2\overline{X_i^{\prime \prime }\rho \dot{X}_i^\text{nuc}}. \end{eqnarray}$$

(7)

|$\sigma _i = \widetilde {X_i^{\prime \prime } X_i^{\prime \prime }}$| is the composition variance of the element i, |$f_i^\sigma = \overline{\rho X_i^{\prime \prime } X_i^{\prime \prime }u_x^{\prime \prime }}$| is the turbulent flux of the composition variance, |$2f_i\partial _x\widetilde {X}_i$| is a source term linked to the flux of the chemical element i, and |$2\overline{X_i^{\prime \prime }\rho \dot{X}_i^\text{nuc}}$| is a source term linked to nuclear reactions.

We discuss below the terms of equation (7), which are presented in Fig. 18.

Figure 18.

Left panel: RA-ILES decomposition of the mean ²⁰Ne abundance variance equation. The meaning of each curve is specified in the legend. The coloured regions have the same meaning as in Fig. 11, and the time average is performed over the same time window. The results come from the 512³ resolution simulation, with a boosting factor of 10. The thin red line shows a modelling of the residuals, done as in Mocák et al. (2018). Right panel: Radial integration of each term shown on the left-hand side of the figure.

4.3.1 Time evolution

|$\bar{\rho }\widetilde {D_t}\widetilde {\sigma }_{i}$| (solid black line) and |$\partial _t(\bar{\rho }\widetilde {\sigma }_{i})$| (dashed black line) represent the Lagrangian and Eulerian time derivatives of the composition variance. The two curves are almost identical throughout the domain, and are close to 0 everywhere, except near the top convective boundary. At this location, entrainment of ²⁰Ne inside the convective region is a source of composition variance.

4.3.2 Transport terms

∇_x(f_i^σ) is the turbulent composition variance flux, and is responsible for the redistribution of the composition variance inside the domain. As shown in Fig. 18 (see the blue curve), the variance produced near the top convective boundary is transported in both directions, outwards inside the stable radiative layer and inside the convective region. It progressively builds the variance peak visible in Fig. 17.

4.3.3 Turbulent production

|$2 f_i \partial _x(\widetilde {X_i})$| term is shown in magenta in Fig. 18, and is the dominant term affecting the creation of composition variance. As this term is effective only in regions where a composition gradient and a turbulent composition flux exist at the same time, it is mostly present near the top convective boundary (and to a much smaller extent near the bottom one). Looking at the integrated values (right panel), one clearly sees that composition variance inside the domain is mostly built by this interaction between the composition flux and the composition gradient.

4.3.4 Nuclear burning

The rate at which composition variance changes due to nuclear burning is traced by |$2\overline{X_i^{\prime \prime }\rho \dot{X}_{i}}$|⁠. In our case, it only occurs at a quite low level near the bottom boundary of the convective zone (see the green curve in the zoom-in on this region), where nuclear burning takes place. This term is negative, indicating that nuclear burning decreases the composition variance.

4.3.5 Numerical residual

|$\mathcal {N}_{\sigma _i}$| represents the residual in our simulation, i.e. the sum of all the other terms (with the same sign as shown in the legend of Fig. 18). As the resolution of our simulation is finite and is therefore limited in spatial resolution, it does not reproduce the exact behaviour of the flow at very small scales, where turbulence is dissipated. The residual thus represents the dissipation at the grid level in our simulations. Despite its numerical origin, the behaviour of this dissipation term is appropriate: we added on the left panel of Fig. 18 a theoretical modelling of this dissipation (thin solid red line), which fits well the residual obtained in our simulations. We adopted the same model as in Mocák et al. (2018), which assumes that the dissipation time-scale is the same as the Kolmogorov damping time-scale.

Here also, the radial integration over the whole domain of the previous terms provides insight into the main mechanisms leading to the change of the variance (here the ²⁰Ne one). We see that the turbulent transport is building the chemical variance inside the domain, while the dissipation at the grid level compensates for this net production.

4.4 Characteristic time-scales

Another way to study the interplay and relative importance of turbulence and nuclear processes is by using characteristic time-scales. The following characteristic time-scales are relevant in this context.

• The convective turnover time-scale τ_conv, which is defined as the typical time needed for a fluid element to cross twice the convective zone (a ‘convective loop’). We thus have

  $$\begin{eqnarray} \tau _\text{conv} = \frac{2d_\text{conv}}{v_\text{rms}}, \end{eqnarray}$$

  (8)

  where d_conv is the radial extend of the convective zone and v_rms is the average root-mean-square velocity inside the convective zone.
• The local nuclear-burning time-scale τ_nuc, which is the typical time-scale for nuclear burning to take place at a given location inside the star:

  $$\begin{eqnarray} \tau _{\text{nuc},i} = \frac{X_i}{\dot{X_i}}, \end{eqnarray}$$

  (9)

  whereas above X_i is the mass fraction of a given chemical element and |$\dot{X}_i$| is the rate of creation/annihilation of the element i due to nuclear reactions.
• The local transport time-scale τ_trans, which is the characteristic time for transport to remove/bring a chemical element at a given place:

  $$\begin{eqnarray} \tau _{\text{trans},i} = \frac{\bar{\rho }\widetilde {X_i}}{\partial _x f_i}, \end{eqnarray}$$

  (10)

  where f_i is the flux of the chemical species.

Based on these characteristic time-scales, one can define the so-called Damköhler (Da) number as the ratio between turbulent transport and nuclear time-scales for a given chemical element:

$$\begin{eqnarray} \text{Da}_i = \frac{\tau _{\text{trans},i}}{\tau _{\text{nucl},i}}, \end{eqnarray}$$

(11)

with the time-scales defined above. The Da number discriminates between two different regimes inside the convective zone: when Da is smaller than 1, the transport time-scale is smaller than the nuclear-burning time-scale. Chemical species are thus well mixed inside the convective region and inhomogeneities (illustrated by the variance) remain small inside the convective zone (except near the boundaries, as explained above). This is how the convective zones look like in 1D stellar evolution. However, it is possible in some cases that the Da number becomes close or even larger than 1. In this case, nuclear burning is faster than transport, and the mixing inside the convective zone is less efficient. This is the ‘convective–reactive regime’ discussed in Herwig et al. (2011) (see also Mocák et al. 2018).

The characteristic time-scales and the Da number are shown in Fig. 19 for ²⁰Ne for our Ex10 simulation. The fastest of the time-scales is the convective time-scale (solid black line), which is slightly less than |$100\, \text{s}$|⁠. The transport time-scale (magenta dot–dashed line) is roughly constant over a significant fraction of the convective zone (about |$10^4\, \text{s}$|⁠), with a peak at the location where the gradient of flux vanishes (see Fig. 15), and has a smaller value near the top boundary (less than |$10^3\, \text{s}$|⁠). Finally, the nuclear-burning time-scale (blue dashed line) is shortest where nuclear burning is strongest, at the bottom of the convective zone (about |$10^4\, \text{s}$|⁠, and increases outwards, to reach |$10^6\, \text{s}$| near the top boundary). The Da number (green dotted line, scale on the right) is thus significantly smaller than 1 everywhere inside the convective region. This means that the mixing of ²⁰Ne is faster than nuclear burning, and that the chemical composition of the convective region (at least concerning ²⁰Ne, but it is true also for the other chemical species followed in this work) is well homogenized. The same is true for the nominal luminosity case, which means that 1D stellar models are sufficient to reproduce standard neon burning (the same is true for the Ex10 and Ex100).

Figure 19.

Characteristic times for the ²⁰Ne. The convective turnover time-scale is shown by the solid black line. The nuclear time-scale is shown in dashed blue. The transport time-scale is shown by the dot–dashed purple line. The Damkhöler number is shown by the dotted green line (value to be read on the right axis). The coloured regions have the same meaning as in Fig. 11, and the time average is performed over the same time window. The results come from the 512³ resolution simulation, with a boosting factor of 10.

The Ex1000 simulation, on the other hand, behaves very differently. This extreme case is illustrated in Fig. 20 for ²⁰Ne (right panel). Variances (left panel) are about 10 times larger than in the Ex10 case shown in Fig. 17 in the whole convective region. In the inner part of the convective zone, nuclear and transport time-scales are about the same, with Da number equal or larger than 1. Looking at Fig. 14, we see that the composition profiles at the end of the simulation (solid lines) for the Ex1000 case are not perfectly flat, but show a gradient in the inner part of the convective zone, as expected from the high Da number. In such cases, nuclear burning and convective mixing can still be treated correctly in 1D models if, for example, these two processes are solved in a coupled manner as done e. g. in the mesa code (Paxton et al. 2011). Non-negligible variances, however, cannot be reproduced in 1D models and such cases need 3D simulations to simulate all the details. While the Ex1000 simulation is not representative of neon burning inside massive stars, there are several cases where we expect Da ∼ 1. In addition to the two cases mentioned above (Herwig et al. 2011; Mocák et al. 2018), other ‘convective–reactive’ environments are expected in various H–He shell interaction events (Hirschi 2007; Clarkson & Herwig 2021) and shell merger situations already found in 1D models (Rauscher et al. 2002; Tur, Heger & Austin 2007; Ritter et al. 2018; Côté et al. 2020), and in the 3D pre-supernova simulation of Yadav et al. (2020).

Figure 20.

Left panel: Relative variance for the four elements followed in this study, in the case of the Ex1000 simulation. Right panel: Characteristic time-scales and Da number for the ²⁰Ne, in the case of the Ex1000 simulation. For both panels, the quantities are averaged between 4 and |$29\, \text{s}$|⁠.

5 DISCUSSION AND CONCLUSIONS

This paper focuses on the 3D hydrodynamics study of a neon-burning shell during the advanced stages of the life of a massive star. We performed a large set of simulations with the prompi code at different resolutions (from 128³ to 1024³ cells), and with different boosting factors for the nuclear energy generation rate (from nominal to 1000 times boosted energy generation rate). The simulations were followed for a long enough time to ensure the flow has significant statistical properties (generally more than five turnover time-scales).

Our simulation at a nominal luminosity allows us to draw conclusions about the flow without any extrapolations. In the initial phases, we find that the convective velocities found in 3D are ≈2 times larger than those expected from MLT approximations. Due to CBM at the convective boundaries, over the course of the simulation (about seven convective turnovers), the convective shell grows noticeably, and the extra fuel entrained leads to an increase in the convective velocity.

A comparison of different resolutions shows that even at the lowest resolution (128³), the large-scale convective flows are reasonably resolved. They attain comparable values of kinetic energy and similar turbulent cascade in the velocity spectra at large wavelengths to their high-resolution counterparts. Using our RA-ILES analysis framework, we also performed a thorough comparison of the mean-field equations at the different available resolutions. As expected from previous works (Viallet et al. 2013; Cristini et al. 2017), an analysis of the mean TKE equation confirms that the statistical properties of the flow are already well captured by the simulation at a relatively low resolution. Higher resolutions, however, are necessary to resolve the convective boundaries and the corresponding mixing across them.

We analyse the mean field composition variance equation, which allows us to determine the degree of deviation of the 3D simulations from 1D averages. We find that the deviation is small (around 1 per cent) in most of the convection zone. At the convective boundary, matter is entrained from the stable region inside the convective region. The boundaries are less homogeneous than the rest of the convective zone, and the variance (in terms of mass fraction) can become as high as about 30 per cent (for silicon). Our mean-field analyses clearly show a net production of variance near the boundaries. This result justifies the general 1D stellar modelling approach for convection but reiterates the need for stellar modellers to use a CBM prescription at convective boundaries.

The various boosting factors we used allowed us to determine how different quantities (flow velocity, buoyancy, and abundance variance) scale with the luminosity of the model. This is similar to the work done by Cristini et al. (2019) and Rizzuti et al. (2022, 2023) for carbon- and neon-burning shells, respectively, however, for the first time, we span a large range of boosting comparing the nominal luminosity runs to those boosted by 1000. As it is not computationally feasible to run very low Mach number flows without boosting the luminosity (e.g. Cristini et al. 2017; Horst et al. 2021), this comparison is important as a first step to determine the feasibility and correctness of very high boosting.

In general, increasing the boosting gave results that are in line with our previous findings (Cristini et al. 2019; Rizzuti et al. 2022, 2023). In particular, the convective velocities, rate of CBM, and the corresponding growth of the convective shell all scale with the boosting applied. The present results [as well as those published in Rizzuti et al. (2022) using the same simulations as the ones presented in this paper] show that these quantities for the non-boosted (nominal luminosity) case can be extrapolated from simulations with a higher boosting factor (typically 10× or 100×), which are cheaper to perform in terms of computing resources. The turbulent spectra in each case all show similar shapes, but with different energies, as expected. Inside the convective zone, we find a very efficient mixing, making the chemical composition in the bulk of the convective region extremely homogeneous for reasonable boosting factors.

For the highest boosting factor we used (Ex1000 simulation), the nuclear burning becomes comparable to the transport time-scale (Da ≈ 1) in a significant fraction of the convection zone. This is qualitatively different from all previous boosting cases, where nuclear burning is only faster than the transport time-scale in a very thin region at the base of the convective shell. In the Ex1000 case, the mixing is not efficient enough to homogenize the convective zone, and a chemical gradient becomes visible in the simulations. This leads to the composition variance showing significant deviations (⁠|${\gt} 10 {{\, \rm per\, cent}}$|⁠) throughout the convection zone, breaking down the approximations that are normally used for 1D stellar evolution calculations. While in our study, this Da ≈ 1 event occurs due to artificial luminosity boosting, they are expected to occur naturally in dynamical processes in stars, such as ingestion events (Herwig et al. 2011; Mocák et al. 2018), and possibly during shell mergers.

Our results show that a correct treatment of convection in 1D stellar evolution codes is required to accurately follow the CBM and reproduce the behaviour appearing in our 3D simulations. The assumption that the mixing is very efficient and almost instantaneous in convective regions appears to be valid, except in cases where the nuclear burning is very fast (Da ∼ 1). Finally, we want to stress that the results presented here are valid for deep convection during advanced stages of stellar evolution, and we need to remain cautious when using them in other phases of the evolution, particularly for convection during the early stage (main sequence), or in stellar envelopes, where radiative effects play a key role.

ACKNOWLEDGEMENTS

RH acknowledges support from the World Premier International Research Centre Initiative (WPI Initiative), MEXT, Japan, and the IReNA AccelNet Network of Networks (National Science Foundation, grant no. OISE-1927130). CG has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Framework Programme (grant no. 833925). WDA acknowledges support from the Theoretical Astrophysics Program (TAP) at the University of Arizona and Steward Observatory. CG, RH, and CM acknowledge ISSI, Bern, for its support in organizing collaboration. This paper is based upon work from the ChETEC COST Action (CA16117) and the European Union’s Horizon 2020 Framework Programme (ChETEC-INFRA, grant no. 101008324). The authors acknowledge PRACE for awarding access to the resource MareNostrum 4 at Barcelona Supercomputing Center, Spain, and the STFC DiRAC HPC Facility at Durham University, UK (grants ST/P002293/1, ST/R002371/1, ST/R000832/1, ST/K00042X/1, ST/H008519/1, ST/K00087X/1, and ST/K003267/1). The University of Edinburgh is a charitable body, registered in Scotland, with Registration No. SC005336. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) Licence to any Author Accepted Manuscript version arising.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

1

For a general introduction to the RA-ILES method, see Chassaing et al. (2002).

2

In this work, y and z are both horizontal directions, while x is directed outwards (opposite to the gravity field).