3D stellar evolution: hydrodynamic simulations of a complete burning phase in a massive star

ABSTRACT

Our knowledge of stellar evolution is driven by one-dimensional (1D) simulations. 1D models, however, are severely limited by uncertainties on the exact behaviour of many multidimensional phenomena occurring inside stars, affecting their structure and evolution. Recent advances in computing resources have allowed small sections of a star to be reproduced with multi-D hydrodynamic models, with an unprecedented degree of detail and realism. In this work, we present a set of 3D simulations of a convective neon-burning shell in a 20 M_⊙ star run for the first time continuously from its early development through to complete fuel exhaustion, using unaltered input conditions from a 321D-guided 1D stellar model. These simulations help answer some open questions in stellar physics. In particular, they show that convective regions do not grow indefinitely due to entrainment of fresh material, but fuel consumption prevails over entrainment, so when fuel is exhausted convection also starts decaying. Our results show convergence between the multi-D simulations and the new 321D-guided 1D model, concerning the amount of convective boundary mixing to include in stellar models. The size of the convective zones in a star strongly affects its structure and evolution; thus, revising their modelling in 1D will have important implications for the life and fate of stars. This will thus affect theoretical predictions related to nucleosynthesis, supernova explosions, and compact remnants.

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

1 INTRODUCTION

To have good understanding and reliable predictions of the evolution and properties of stars, researchers have developed stellar evolution models that can follow the life cycle of an entire star from birth to death (e.g. Heger & Woosley 2002; Paxton et al. 2011; Ekström et al. 2012). However, there are many physical processes to take into account that can affect the evolution and fate of a star. To make the computation possible and affordable, stellar models rely on simplifications and assumptions, one of the most important of which is spherical symmetry. An equally important though less easily envisages simplification is the treatment of convection in stars with mixing length theory (MLT; Böhm-Vitense 1958).

To improve and refine these prescriptions, sections of a star may be studied in more detail employing multidimensional hydrodynamic models, which can realistically reproduce complex physical processes without using any of the prescriptions of the one-dimensional (1D) case (e.g. José & Hernanz 1998; Asplund et al. 2009; Lentz et al. 2015; Janka, Melson & Summa 2016; Burrows & Vartanyan 2021). Unfortunately, these simulations are computationally very expensive, requiring millions of core-hours to simulate hours of stellar evolution, so reproducing the evolutionary time-scales has so far been inaccessible. Hydrodynamic simulations of stellar interiors mostly focus on reproducing turbulent convection and mixing between layers (known as ‘convective boundary mixing’, CBM), resulting in the growth over time of convective zones due to entrainment of material from the stable regions into the convective ones. This is particularly impactful for the late phases of massive stars, where convective speeds are high and entrainment is vigorous, affecting the structure of the star and the end of its life. The very last phases before the onset of iron-core collapse are often studied with multi-D hydrodynamic supernova (SN) simulations (Couch et al. 2015), which recently started to take into account also the effects of rotation (Yoshida et al. 2021; McNeill & Müller 2022) and magnetic fields (Varma & Müller 2021) on the convective motions.

CBM is also important for the convective cores of main-sequence stars, with implications for the stellar lifetime and type of remnant (Kaiser et al. 2020), but here the convective velocities are much smaller and it is difficult for hydrodynamic models to simulate significant time ranges without boosting the rate of energy release to accelerate the simulations (Meakin & Arnett 2007; see also the study in Baraffe et al. 2023), as is also often done for later phases (Cristini et al. 2017; Horst et al. 2021).

To validate results, it is crucial to investigate stellar convection in hydrodynamic simulations that do not modify the energy release. Previous studies of 3D stellar simulations do alter the initial conditions introducing some energy boosting, with some exceptions such as the very late oxygen and silicon convective shells run for short time-scales (Couch & Ott 2015; Müller et al. 2017; Yoshida et al. 2019). However, Rizzuti et al. (2022), with their analysis of CBM in a nominal-luminosity neon shell, revealed the presence of strong entrainment in the convective zone, confirming it is not an effect of the boosting in luminosity but a natural result of the convectively unstable structure predicted by the 1D stellar models. This shows the importance of correctly implementing entrainment also in 1D models, to ensure consistency and agreement between results. We call this approach ‘321D-guided’, referring to 1D stellar models that have been improved using the results of 3D hydrodynamic simulations.

It is worth noting though that the nominal-luminosity simulation of Rizzuti et al. (2022) has been run only up to 3000 s, which while long enough for studying CBM with good statistics is too short to give any information on the evolution of the burning shell towards fuel exhaustion. Until now, multi-D studies have not investigated the interplay between entrainment and fuel exhaustion at the end of a burning phase, leaving many questions about the fate of convection unanswered. Specifically, does convection stop when the fuel is exhausted or does the mixing of fresh fuel from entrained material extend the convective growth indefinitely? This issue, on which 1D models cannot give any information, has a pivotal role for correctly predicting the later evolution of the star, in particular the pre-SN structure and the type of remnant.

We organize the paper as follows: In Section 2, we present the initial conditions and general set-up of the hydrodynamic simulations. In Section 3, we analyse the results from the sets of simulations. Finally, we discuss results and draw conclusions in Section 4.

2 METHODS

2.1 Initial conditions from a 1D stellar model

To run a stellar hydrodynamic simulation, initial conditions need to be assumed from a 1D stellar evolution model that simulates the entire lifetime of the star, so that the realism of the simulated environment can be ensured. For this purpose, we used the mesa stellar evolution code (Paxton et al. 2011, 2013, 2018, 2019) to model the evolution of a 20 M_⊙ star of solar metallicity (Z = 0.014 using the relative abundances of Asplund et al. 2009) from the pre-main sequence until core collapse. Mass-loss rates were taken from the ‘Dutch’ options. This includes several prescriptions: for O-stars, the mass-loss rates from Vink, de Koter & Lamers (2000, 2001) are used; if the star enters the Wolf–Rayet stage, i.e. when the surface hydrogen mass fraction drops below 0.4, the mass-loss rate switches to the scheme from Nugis & Lamers (2000); if T_eff < 10⁴ K, the empirical mass-loss rate from de Jager, Nieuwenhuijzen & van der Hucht (1988) is used. The MLT (Böhm-Vitense 1958) of convection describes the treatment of convection in our model (using the ‘Henyey’ and ‘MLT++’ options), where we applied an efficiency of α_MLT = 1.67 (Arnett et al. 2018). The Schwarzschild criterion defines the convective boundaries in our model and as such we did not need to implement semiconvective mixing. For CBM, we included the exponential decaying diffusive model (Freytag, Ludwig & Steffen 1996; Herwig 2000) with f_ov = 0.05 for the top of convective cores and shells and with f_ov = 0.01 for the bottom of convective shells (with f₀ = f for both cases). We furthermore used the decline of the diffusion coefficient near the boundary (Jones et al. 2017). The value of f_ov = 0.05 for the top boundaries is larger than the majority of published large grids of stellar models (e.g. using α_ov = 0.1 in Ekström et al. 2012, α_ov = 0.335 in Brott et al. 2011). The value of 0.05 is motivated by the study of Scott et al. (2021), where values for f_ov up to at least 0.05 for 20 M_⊙ and above (see caption of their fig. 9 and references therein for the relation between α_ov, f_ov, and entrainment) best reproduce the observed width of the main sequence in the spectroscopic Hertzsprung-Russell diagram (Castro et al. 2014). For the bottom boundary, a CBM value of 1/5 the value of the top boundary is based on 3D hydrodynamic simulations (Cristini et al. 2017, 2019; Rizzuti et al. 2022), finding that CBM is slower at the bottom boundary due to it being stiffer and therefore harder to penetrate. Scott et al. (2021) also show that the amount of CBM should increase with initial mass since more massive stars are much more luminous (L ∼ M³ between 1 and 20 M_⊙) and thus our chosen value of 0.05 for 20 M_⊙ is consistent with the range of values inferred from asteroseismology for less massive stars (f_ov = 0.02–0.04; Bowman 2020). As we describe below, the inclusion and the extent of CBM, guided by both 3D hydrodynamic simulations and observations (main-sequence width in Hertzsprung-Russell diagram and asteroseismology), is a key aspect for the novelty of the results.

2.2 The 3D hydrodynamic model set-up

The 1D simulation has been run from the pre-main-sequence phase until core collapse (see the structure evolution diagram in Fig. 1, left). The hydrodynamic simulations have been started with initial conditions taken from the first neon-burning shell that develops in the 1D simulation (see Fig. 1, right). We employ here the stellar hydrodynamical code prompi (Meakin & Arnett 2007), successfully used over the years to simulate and study convection and CBM in advanced phases of massive stars (e.g. Arnett, Meakin & Young 2009; Cristini et al. 2017, 2019; Mocák et al. 2018; Rizzuti et al. 2022). prompi has been also recently compared to other stellar hydrodynamical codes in Andrassy et al. (2022), who showed prompi to be fully consistent with the other codes.

Figure 1.

Left: Structure evolution diagram of the 20 M_⊙ 1D mesa input stellar model as a function of the time left until the predicted collapse of the star (in years, log scale). In blue the convective zones, and in green the CBM zones. The red arrow indicates the neon-burning shell the 3D simulations were started from. Top right corner is a zoom-in on the neon-burning shell. Right: Zoom-in on the model used as initial conditions for the 3D simulations. The horizontal axis is the time in seconds relative to the start of the 3D simulations. The vertical axis is the radius in units of 10⁸ cm. In colour scale, the squared convective velocity. Isomass contours are shown as black lines. The lines show that the shell undergoes significant expansion during the Ne-burning phase. The vertical red bar indicates the start time and radial extent of the hydrodynamic simulations.

The radial variables density, pressure, temperature, entropy, mass, and chemical composition have been remapped on a tri-dimensional grid adding relatively small perturbations (∼10⁻⁷) to density and temperature as seeds for convective instabilities. The grid we used is in spherical coordinates with a radial extension from 3.6 to 8.5 × 10⁸ cm and an angular size of about 26° in both angular dimensions. Since the radial size of the grid is roughly twice the size of the other dimensions, we used a resolution with double the cells in radius. Making use of the spherical coordinates, gravity is recomputed at each time-step by integrating the mass inside each shell as a function of the radius. This allows for the contraction or expansion of the layers. A resolution of 256 × 128² grid points in r, θ, and ϕ, respectively, was initially used to explore the evolutionary time-scale of the simulations, and then for this study we analysed a more detailed set of simulations with resolution 512 × 256², and finally results have been validated with our most detailed simulations of 1024 × 512² and 2048 × 1024².

Convection has been fuelled with nuclear energy. Nuclear burning and nucleosynthesis are not always explicitly tracked in hydrodynamic simulations, because of the high computing cost: instead, fixed heating profiles are often used to drive convection. Rizzuti et al. (2022) included a simple energy generation routine of five isotopes (⁴He, ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si) to reproduce neon burning. We have now extended this network to include 12 isotopes: n, p, ⁴He, ¹²C, ¹⁶O, ²⁰Ne, ²³Na, ²⁴Mg, ²⁸Si, ³¹P, ³²S, and ⁵⁶Ni. We employed the most recent nuclear rates from the JINA REACLIB data base (Cyburt et al. 2010). While this extension does not have important effects on the energy release, which was already accurate with the five-isotopes burning routine, it allows us to study the nucleosynthesis and transport of other species, paving the way towards an extended multi-D nucleosynthesis in stellar models.

Finally, following the same approach as Rizzuti et al. (2022), in this study we included simulations with and without modified energy generation. In the nominal-luminosity simulations, the luminosity remains the same as in the 1D model. In this way, we make sure to validate our results ruling out the possibility that the boosting in luminosity introduces additional differences from the nominal case. When a boosting in luminosity was included, the nuclear rates for neon burning ²⁰Ne(γ, α)¹⁶O and ²⁰Ne(α, γ)²⁴Mg have been multiplied by a boosting factor. Since these are the reactions that dominate the energy release, it does not make any difference that the other reactions have not been boosted. The boosting factors used here are 1 (nominal-luminosity), 5, 10, and 50.

2.3 Entrainment law computation

prompi has a long history of modelling entrainment in stellar environments, having investigated a massive star O-shell (Meakin & Arnett 2007), C-shell (Cristini et al. 2019), and Ne-shell (Rizzuti et al. 2022). We can express the entrainment rate for a convective boundary, i.e. entrainment velocity v_e divided by convective velocity v_rms, according to Meakin & Arnett (2007):

$$\begin{eqnarray} E = \dfrac{v_\text{e}}{v_\text{rms}} = A\ Ri_\text{B}^{-n} \end{eqnarray}$$

(1)

as a function of the ‘bulk Richardson number’ Ri_B, which we can see as a measure of boundary resistance against convective fluid penetration, and defined as

$$\begin{eqnarray} Ri_\mathrm{ B}=\dfrac{\ell \Delta b}{v_\text{rms}^2}\ ; \quad \Delta b=\int \limits _{r_\text{b}-\ell /2}^{r_\text{b}+\ell /2} N^2 \mathrm{ d}r \end{eqnarray}$$

(2)

with Δb the buoyancy jump across the convective boundary, N the Brunt–Väisälä frequency, ℓ the length-scale of turbulent motions, and r_b the convective boundary location.

For the computation of the quantities just defined, we refer to Cristini et al. (2019) and Rizzuti et al. (2022) since we used here their same approach and definitions. It is worth mentioning here that r_b has been obtained using the neon abundance gradient, that v_e is the time derivative of r_b, that v_rms is equal to |$\left(v_r^2+v_\theta ^2+v_\phi ^2\right)^{1/2}$| inside the convective zone, that ℓ has been set equal to 1/12 of the local pressure scale height to enclose the peak in N², and that all quantities have been averaged for the entire time the simulations spent in the entrainment regime.

3 RESULTS

In this article, we present a set of 3D hydrodynamic simulations run continuously from the beginning to the end of a neon-shell burning phase. The level of realism of these new simulations is guaranteed by the updated 1D model for initial conditions, which includes a stronger CBM than usually implemented, the complex burning routine including 12 isotopes for energy generation, the presence of unaltered nominal-luminosity runs, spherical geometry, and some high-resolution runs.

In Table 1, we summarize the most important properties of the hydrodynamic simulations we have run for this study, each of them used for some part of the analysis, contributing to the final results we present here. The code name for each simulation summarizes its radial resolution (r) and the boosting factor (e).

Table 1.

Properties of the 3D hydrodynamic simulations presented in this study: model name; resolution N_rθϕ; boosting factor of the driving luminosity ϵ; starting t_start and ending t_end time of the simulation; time spent in the entrainment regime τ_entr; convective turnover time τ_conv; number of convective turnovers simulated in the entrainment regime n_conv; root-mean-square convective velocity v_rms; sonic Mach number Ma; and number of CPU core-hours required to run the simulation.

Name	Nrθϕ	ϵ	tstart (s)	tend (s)	τentr (s)	τconv (s)	nconv	vrms (cm s−1)	Ma (10−2)	Mcore-h
r256e1	256 × 1282	1	0	60 000	15 000	155	96	3.29 × 106	0.83	2.08
r256e5	256 × 1282	5	0	29 000	1500	59	25	6.55 × 106	1.76	0.89
r256e10	256 × 1282	10	0	19 000	800	50	16	8.06 × 106	2.15	0.60
r256e50	256 × 1282	50	0	30 000	150	30	5	1.31 × 107	3.48	0.96
r512e1	512 × 2562	1	16 000	19 000	3000	136	22	3.83 × 106	0.99	1.66
r512e5	512 × 2562	5	0	2000	1500	59	25	6.65 × 106	1.80	0.80
r512e10	512 × 2562	10	0	1000	800	49	16	8.28 × 106	2.23	0.50
r512e50	512 × 2562	50	0	490	150	30	5	1.34 × 107	3.61	0.20
r1024e1	|$1024\times 512^2$|	1	10 000	10 400	400	127	3	3.26 × 106	0.84	2.88
r2048e1	2048 × 10242	1	10 010	10 030	20	113	0	3.85 × 106	0.99	2.02

3.1 Analysis of the fluid dynamics

To provide a visual representation of our simulations, we show in Fig. 2 a cross-section of the neon mass fraction from our high-resolution nominal-luminosity r2048e1 simulation. In addition to showing the fine detail that hydrodynamic simulations can reveal with modern computing resources, this image helps us understand the importance of CBM for stellar evolution. Indeed, it is possible to see the mixing of neon-rich and neon-poor material inside the convective zone, as well as the entrainment of neon-rich material from the upper stable region, due to shear mixing at the interface of the two layers. The entrainment of fresh fuel extends the nuclear burning time-scale and therefore the convective shell lifetime, showing the importance of including CBM in all stellar models.

Figure 2.

Cross-section of the neon mass fraction (values in colour scale) from the r2048e1 simulation. The frame shows entrainment of some neon-rich material from the upper stable region into the convective zone. Videos of the evolution of the r1024e1 simulation are available online as Supplementary material, showing fluid speed and mass fractions in colour scale.

The effects of the mixing on the convective boundaries can also be seen in Fig. 3, where we show the difference in the radial profiles of mean atomic mass and entropy between the 1D initial conditions and a 3D r512e1 test simulation. This 3D r512e1 test simulation was run for five convective turnovers to highlight similarities and differences between the 1D and 3D profiles near the start of the simulations (the r512e1 simulation listed in Table 1 and discussed in the rest of the paper was restarted from the r256e1 simulation after 16 000 s to make the most of our computing budget). We added in Fig. 3 the grid points on the 3D simulation profiles, to indicate the resolution of the simulation across the convective boundary. The reason why the values for atomic mass and entropy in the convective plateaus at r < 5.4 × 10⁸ cm do not match perfectly (∼10⁻⁴ absolute difference) is the remapping and the small mixing that takes place during the initial transient. Importantly, what we can see from this plot is that, as the boundary moves outwards due to CBM, both chemical composition and entropy are consequently mixed. This point is also underlined by Fig. 4, which shows the temperature gradients for the same models. Since both 1D and 3D stellar models deviate from the adiabatic temperature gradient outside the convective region (although the 1D curve is affected by numerical noise), it is clear that CBM in the 3D simulation contributes to altering the temperature gradient in the ‘overshooting’ region, which becomes adiabatic due to the fast entropy mixing. These plots show that in the 1D mesa code the composition mixing is accurately modelled, since the shape of |$\bar{A}$| is similar to the 3D, but the entropy profile is not compatible, so it will be necessary in the future to also add the mixing of entropy in the 1D models.

Figure 3.

Radial profiles of the mean atomic mass |$\bar{A}$| (black) and entropy (blue) for the 1D input model (dashed lines) and at the end of a 3D r512e1 test simulation (angularly averaged quantities plotted as solid lines) run for five convective turnovers. The dots on the 3D curves indicate the location of the simulation cells.

Figure 4.

Same as Fig. 3, but radial profiles of the actual temperature gradient (red) and the adiabatic temperature gradient (yellow), defined as ∇ = dln T/dln P, for the 1D input model (dashed lines) and the 3D r512e1 test simulation (angularly averaged quantities plotted as solid lines).

Another way of studying the evolution of the different simulations is to analyse how the kinetic energy of the fluid builds up and changes with time. In Fig. 5, we show the time evolution of the specific turbulent kinetic energy, defined as |$\dfrac{1}{2}\left(\overline{v_r^2}-\overline{v_r}^2+\overline{v_\theta ^2}-\overline{v_\theta }^2+\overline{v_\phi ^2}-\overline{v_\phi }^2\right)$| where the single quantities have been averaged in r, θ, and ϕ inside the convective zone. The figure shows the most complete set of simulations we have run, which consists of r256e1, r256e5, r256e10, and r256e50. These simulations have the same resolution but different luminosity boosting, and this explains the different evolution of the tracks in the figure. After an initial transient, during which convection develops in the 3D grid and the kinetic energy builds up, all simulations evolve in a similar way, with a more or less gradual increase in kinetic energy during the nuclear burning phase, and a slow decrease after neon is consumed and nuclear burning does not sustain the kinetic energy any longer. The differences in magnitude and time-scale can be all traced back to the boosting in luminosity, which strongly affects both the amount of kinetic energy produced during the simulation, and how rapidly neon is consumed, therefore the time-scale for nuclear burning. The peak in kinetic energy is approximately at 20 000, 2200, 960, and 180 s for r256e1, r256e5, r256e10, and r256e50, respectively.

Figure 5.

Time evolution of the specific turbulent kinetic energy |$\dfrac{1}{2}\left(\overline{v_r^2}-\overline{v_r}^2+\overline{v_\theta ^2}-\overline{v_\theta }^2+\overline{v_\phi ^2}-\overline{v_\phi }^2\right)$| for the four simulations with different luminosity boosting, r256e1, r256e5, r256e10, and r256e50.

It is also interesting to notice that towards the end all simulations seem to converge to the same constant value of turbulent kinetic energy, regardless of the luminosity boosting. This comes from the fact that after neon is exhausted the reactions ²⁰Ne(γ, α)¹⁶O and ²⁰Ne(α, γ)²⁴Mg stop occurring, so the only nuclear reaction that can proceed is the secondary ²⁴Mg(α, γ)²⁸Si, which was not boosted and has the same rate for all simulations. This is confirmed by a slight decrease in magnesium and increase in silicon abundances during this late phase. However, this reaction is not energetic enough; therefore, turbulence decays and the shell growth halts.

The time evolution of the turbulent kinetic energy is also shown in Fig. 6 for the same set of simulations, but as a colour map with angularly averaged values in colour scale. First, it is interesting to note the strong effect that the luminosity boosting has on the evolutionary time-scale. Apart from this, the shell is evolving in the same way in all simulations, with the first period being dominated by nuclear burning and entrainment, and demonstrating a linear growth of the shell. An additional weak burning front is visible around r ∼ 8 × 10⁸ cm, produced by the impact of gravity waves on a carbon shell above the neon one, but its energy is so small compared to the neon-burning shell (three orders of magnitude lower) that it has no impact on convection and entrainment. Once fuel is exhausted in the neon-burning shell (neon abundance is 6 per cent at the peak of kinetic energy) convection slowly diminishes, as evident from the drop in kinetic energy, and therefore entrainment also stops, halting the shell growth. This result has important implications for stellar evolution. In contrast to previously suppositions (e.g. Cristini et al. 2019; Horst et al. 2021; Rizzuti et al. 2022), CBM does not lead to the convective engulfment of the entire star, but our simulations show that the shell naturally stops growing when its fuel is exhausted. Indeed, making use of network calculations with a one-zone model to exclude convective mixing, we estimated the nuclear burning time-scale (⁠|$X_\text{Ne}/\dot{X}_\text{Ne}$|⁠) to be approximately 4000 s in this environment, which is much shorter than the time-scale for mass entrainment (⁠|$M_\text{e}/\dot{M}_\text{e}$|⁠), here around 30 000 s; therefore, entrainment cannot sustain convection on its own. This finding puts a limit on the size of convective zones, which can have a strong effect on the final structure of massive stars.

Figure 6.

Time evolution of the angularly averaged specific turbulent kinetic energy (in colour scale) for four simulations with different luminosity boosting factors: from top left to bottom right, r256e1, r256e5, r256e10, and r256e50. Overlaid in white, the isomass contours.

There are other important points that can be drawn from Fig. 6. Comparing the nominal-luminosity run (first panel) to the corresponding evolution in the 1D model (Fig. 1, right), it is evident that the convective shell in the 3D model evolves ∼5 times faster than that in its 1D equivalent. To understand this difference, in Fig. 6 we overlay isomass contours in white that can be used to track expansion. The contours show some expansion in the bulk of the convective zone, but no expansion is allowed closer to the upper and lower domains, because mass flow is not allowed through those closed boundaries. If compared again to Fig. 1 (right), it is clear that an expansion is present in the 1D layers that are not limited by domain boundaries, but which are free to contract or expand. It is this difference in convective zone size between 1D and 3D that explains the difference in burning time-scales. Indeed, the neon-burning energy release is strongly dependent on the temperature, sensitive to a power of ∼T⁵⁰ (Woosley, Heger & Weaver 2002) due to the temperature dependence of nuclear reaction rates and the α-particle mass fraction. Furthermore, the temperature of a gas is dependent on its volume, according to the equation of state. It is possible to show with a simple calculation that the volume difference between the final states of our convective zone in 1D versus 3D can explain the shorter time-scale of the 3D simulation. If we assume for simplicity that the two states are separated by an adiabatic expansion (it is reasonable to assume no heat exchange with the surroundings), then

$$\begin{eqnarray} \dfrac{T_\text{3D}}{T_\text{1D}} =\left(\dfrac{V_\text{1D}}{V_\text{3D}} \right)^{\gamma -1} =\left(\dfrac{R^3_\text{1D}-r^3_\text{1D}}{R^3_\text{3D}-r^3_\text{3D}} \right)^{\gamma -1} , \end{eqnarray}$$

(3)

where r and R are the inner and outer radii of the shells, respectively. Comparing the 1D and 3D states at the end of the nominal-luminosity neon burning (taken when the neon abundance in the bulk of the convective zone X_Ne ∼ 0.015), we have (in units of 10⁸ cm)

$$\begin{eqnarray} R_\text{1D}=7.34,\quad r_\text{1D}=4.63,\quad R_\text{3D}=6.93,\quad r_\text{3D}=4.44 \end{eqnarray}$$

and γ = 1.55, which gives T_3D/T_1D = 1.11. At the same point in the actual simulations, we find T_3D = 1.88 GK and T_1D = 1.78 GK at the temperature peak (located near the bottom of the convective shell) corresponding to a ratio of 1.06. Thus, we can conclude that the limited expansion in 3D due to the closed boundary conditions can account for the higher temperature found in the 3D simulations compared to the 1D stellar model.

Since the nuclear energy generation rate depends on a power of the temperature |$\dot{\epsilon }\sim T^\alpha$|⁠, the difference in nuclear burning time-scale |$\tau \sim 1/\dot{\epsilon }$| between 1D and 3D can be approximated as

$$\begin{eqnarray} \dfrac{\tau _\text{3D}}{\tau _\text{1D}} \sim \left(\dfrac{T_\text{1D}}{T_\text{3D}} \right)^\alpha. \end{eqnarray}$$

(4)

Using α = 50 from Woosley et al. (2002) results in a 3D time-scale 15 times shorter than that for the 1D, which is much faster than what we see in our simulations. Realistically, the energy generation rate is also dependent on the neon abundance, which decreases with time; therefore, α < 50 is expected, although the value changes in time. While both the 1D and 3D simulations are complex, this simple estimate of the effect on the temperature of the different expansion shows that this can explain the difference in the time-scales of neon burning between the 3D and 1D simulations.

3.2 Spectral analysis and turbulence theory

The well-established ‘Kolmogorov theory’ (Kolmogorov 1941) states that, for a quasi-steady isotropic regime of convection, the rate of energy dissipation is independent of the scale and of the type of dissipative process, and the kinetic energy is expected to behave according to

$$\begin{eqnarray} E_\text{K}\sim v^2_\text{rms}\sim k^{-5/3} , \end{eqnarray}$$

(5)

where k is the wavenumber associated with the fluid scale. The fluid is expected to follow this scaling throughout the so-called ‘inertial range’, while it deviates from it at the smallest (dissipation) scales due to dissipative effects, and at the largest (integral) scales because it stops being isotropic. These premises allowed hydrodynamic codes to employ the implicit large eddy simulation (ILES) scheme, which replaces the explicit physical viscosity with implicit numerical viscosity due to the finite grid resolution. An important advantage of ILES is that it overcomes the necessity of resolving the flow at the viscosity scale, which would be impossible for stellar simulations.

It is thus possible to determine the quality of our simulations in the inertial range by comparing the power spectrum of the turbulent velocity in the simulations to the scaling expected from the Kolmogorov theory. Since we are employing spherical geometry, the default method would be to compute the spherical harmonic decomposition of v_rms. However, our simulations cover only a small section of a spherical surface (0.2 sr, or 2 per cent of the total spherical surface); therefore, it is not trivial to do the decomposition, which would also be unable to show the lowest order modes due to the limited solid angle covered. One way of doing that would be to repeat the pattern periodically to fill a full spherical surface, as done in Horst et al. (2021), but this would introduce artefacts coming from the spherical mapping and the risk of underestimating the low-order modes. For these reasons, we prefer here to compute the spectrum using a 2D Fourier analysis, following a similar approach to e.g. Cristini et al. (2019) and Andrassy et al. (2022). Selecting a radius, r = 5 × 10⁸ cm, in the middle of the convective region, a 2D Fourier transform of the velocity magnitude as a function of the angular coordinates θ, ϕ was computed:

$$\begin{eqnarray} \hat{v}_\text{rms}(k_\theta ,k_\phi)=\dfrac{1}{N_\theta N_\phi }\sum _{n_\theta =0}^{N_\theta -1}\sum _{n_\phi =0}^{N_\phi -1}v_\text{rms}(\theta ,\phi){\rm e}^{-i2\pi \left(\dfrac{k_\theta n_\theta }{N_\theta }+\dfrac{k_\phi n_\phi }{N_\phi }\right)} , \end{eqnarray}$$

(6)

where N_θ, N_ϕ are the numerical resolution, n_θ, n_ϕ are the cell numbers, and k_θ, k_ϕ are the wavenumbers, which span the range:

$$\begin{eqnarray} \begin{aligned} & k_\theta = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}i, & {{\rm if}\quad } 0 \le i \lt N_\theta /2 \\ i-N_\theta , & {{\rm if}\quad } N_\theta /2 \le i \lt N_\theta \end{array}\right.\\ & k_\phi = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}j, & {{\rm if}\quad } 0 \le j \lt N_\phi /2 \\ j-N_\phi , & {{\rm if}\quad } N_\phi /2 \le j \lt N_\phi \end{array}\right. \end{aligned} . \end{eqnarray}$$

(7)

Finally, we plot in Figs 7 and 8 the power spectrum |$\dfrac{1}{2}|\hat{v}_\text{rms}|^{2}$|⁠, which can be interpreted as the specific kinetic energy, as a function of the wavenumber |$k=\sqrt{k_\theta^2+k_\phi^2}$|⁠. Since k_θ ∈ [−N_θ/2, N_θ/2] and likewise for k_ϕ, and the norm k draws a circle in the (k_θ, k_ϕ) space, we limit the plot to k ∈ [0, min{N_θ/2, N_ϕ/2}] to avoid the circle going beyond the domain and losing a fraction of the signal, resulting in a drop of the power spectrum. Spectra have been averaged over at least one convective turnover for all simulations, except r2048e1 that has been run for a very short time-scale given its very high resolution, so it has been averaged for the last 10 s.

Figure 7.

Spectra of the specific kinetic energy as a function of both the wavenumber k and the real space x, taken on a surface at fixed radius, for the simulations with different resolutions r256e1, r512e1, r1024e1, and r2048e1. The spectra were averaged over 200 s, except for r2048e1 that was averaged over 10 s. The dashed black line is the Kolmogorov scaling k^−5/3; the vertical dotted line at k ∼ 2 is the size of the convective region; and the vertical dotted lines at k ∼ 8–70 correspond to the size of 15 cells for each simulation.

Figure 8.

Same as Fig. 7, but for simulations with same resolution and different boosting factors: r256e1, r256e5, r256e10, and r256e50. All spectra were averaged over one convective turnover. The two vertical lines identify the inertial range, as in Fig. 7.

Fig. 7 shows the spectra for simulations with boosting factor equal to 1 (nominal luminosity) but different resolutions. In all simulations, the bulk of the spectrum follows the expected Kolmogorov scaling, which is a good confirmation of the presence of a large inertial range in our simulations. Also, as expected, the spectrum starts deviating from the k^−5/3 scaling both at the largest scales (around the vertical k ∼ 2 line) and at the smallest scales, because of the numerical dissipation near the grid scale. For this reason, as resolution increases in Fig. 7 the inertial range extends towards larger wavenumbers, because dissipation takes place on smaller and smaller spatial scales. The point where the spectrum slope starts deviating strongly from the k^−5/3 scaling corresponds approximately to 15 cells for all simulations, as we indicate in the plot with vertical dotted lines around k ∼ 8–70.

In Fig. 8, we present instead the spectra for a set of simulations with same resolution but different boosting factors. As expected, the specific kinetic energy of the spectra increases with the boosting factor, but the extent of the inertial range does not change, since the resolution is the same. This is a confirmation of the fact that introducing a boosting factor does not affect the general properties of the turbulent flow (apart from the magnitude of the kinetic energy). These findings are perfectly in line with previous simulations (e.g. Cristini et al. 2019) and especially with results from the code comparison study of Andrassy et al. (2022).

3.3 Entrainment analysis and parametrization

As mentioned earlier, entrainment of fresh fuel into the convective zone is one of the most important effects on stellar convection, profoundly affecting the stellar structure and its evolution. We remind that the new simulations we are presenting in this work have been started from initial conditions taken from a 1D model using strong CBM at all convective boundaries, including the late-phase shells. This last point is particularly important in order to understand how entrainment differs between 1D and 3D stellar models, since many studies underline strong discrepancies between the two (Staritsin 2013; Higl, Müller & Weiss 2021; Horst et al. 2021; Scott et al. 2021; Rizzuti et al. 2022). In particular, entrainment in 3D is always found to be much stronger than in 1D models, which often include little or no CBM at all. Starting from a 1D model with strong CBM, we can determine whether the corresponding entrainment in 3D is larger as usual or whether a convergence between 1D and multi-D stellar models can be achieved.

In Fig. 9, we present the entrainment rates estimated from the data presented in this work (blue), alongside one of the previous prompi simulations of a neon shell from Rizzuti et al. (2022) (red), and the 1D study of a convective hydrogen core in a 15 M_⊙ star from Scott et al. (2021) (black). To study entrainment, we used the data coming from the r512e1, r512e5, r512e10, and r512e50 simulations, which have been run for many convective turnovers with a high resolution. For each simulation, we computed the entrainment rate and bulk Richardson number for both the upper and lower convective boundaries, averaged through the entire entrainment phase; we list results in Table 2. As expected, a larger boosting factor in the simulations results in higher convective and entrainment velocities and smaller Ri_B due to the larger penetrability. In addition, the upper boundary has always larger entrainment rate and smaller Ri_B than the lower one. Error bars in the figure are standard deviations of the values at each time-step in our simulations, and since both their computation and the fitting have been done in real space, in some cases the log scale of the plot shows the bars going towards zero.

Figure 9.

Entrainment rate versus bulk Richardson number. Data from stellar simulations and respective linear regressions: ‘mesa’ Ne-shell from this study (blue), ‘GENEC’ Ne-shell from Rizzuti et al. (2022) (red), and H-core from 1D Scott et al. (2021) (black). Triangles are lower convective boundaries, circles are upper boundaries. The dashed vertical line indicates the bulk Richardson number in the convective H-core. Error bars are standard deviations. The legend shows parameter estimates for the entrainment law |$v_\text{e}/v_\text{rms} = A\ Ri_\text{B}^{-n}$|⁠.

Table 2.

Properties of the 3D hydrodynamic simulations used to study entrainment in this work: model name; root-mean-square convective velocity v_rms (cm s⁻¹); upper entrainment rate |$v_\text{e}^\text{up}/v_\text{rms}$|⁠; lower entrainment rate |$v_\text{e}^\text{low}/v_\text{rms}$|⁠; upper bulk Richardson number |$Ri_\text{B}^\text{up}$|⁠; and lower bulk Richardson number |$Ri_\text{B}^\text{low}$|⁠.

Name	vrms	|$v_\text{e}^\text{up}/v_\text{rms}$|	|$v_\text{e}^\text{low}/v_\text{rms}$|	|$Ri_\text{B}^\text{up}$|	|$Ri_\text{B}^\text{low}$|
r512e1	3.83 × 106	1.01 × 10−3	5.38 × 10−5	51.3	224
r512e5	6.65 × 106	5.03 × 10−3	3.69 × 10−4	13.8	64.7
r512e10	8.28 × 106	8.25 × 10−3	6.54 × 10−4	8.91	42.5
r512e50	1.34 × 107	2.72 × 10−2	1.84 × 10−3	2.63	15.3

It is evident from Fig. 9 that the value of Ri_B for convection in late phases (data points) is several orders of magnitude smaller than that during the main sequence (vertical line): the lack of entrainment data and the consequent need for extrapolation are the main reasons for the current disagreement on CBM between different stellar phases. From the comparison of the different trends in Fig. 9, several interesting conclusions may be drawn. The entrainment rates measured from our simulations are lower than all the previous multi-D studies done with prompi, as we show in Fig. 10 where we compare entrainment in all prompi simulations so far. Although our new rates are not as small as the ones predicted from studying the H-core in 1D (black line in Fig. 9), the larger steepness and lower dispersion of the new results imply much less entrainment than the previous studies, especially at larger Ri_B, reaching a surprisingly good agreement with predictions for the convective core in 1D models (dashed vertical line). This finding is an important step towards convergence of results between 1D and 3D stellar models.

Figure 10.

Same as Fig. 9, but comparison between prompi simulations: ‘mesa’ Ne-shell from this study (red, solid), ‘GENEC’ Ne-shell from Rizzuti et al. (2022) (blue, dotted), C-shell from Cristini et al. (2019) (green, dashed), and O-shell from Meakin & Arnett (2007) (in yellow, dot–dashed). In the legend, parameter estimates for the entrainment law |$v_\text{e}/v_\text{rms} = A\ Ri_\text{B}^{-n}$|⁠.

To better understand the reasons for this convergence, it is important to underline the differences between the previous and the new simulations. The two sets of hydrodynamic simulations we show in Fig. 9 are both of a Ne-burning shell, with a similar burning network and energy release, but with initial conditions taken from two different 1D models: one (Rizzuti et al. 2022, red) from a GENEC model with no CBM for this phase and the other (this study, blue) from a mesa model with strong CBM. The stellar mass is also different (15 M_⊙ for Rizzuti et al. 2022 and 20 M_⊙ for this study), but previous prompi simulations show that stellar mass does not seem to affect the entrainment law parameters (Cristini et al. 2019; Rizzuti et al. 2022). Furthermore, the present set of simulations has been run for much longer than before, covering the entire shell evolution until fuel exhaustion. What we can conclude from this is that a hydrodynamic simulation started from a 1D model already including CBM produces significantly less entrainment than simulations from a 1D model with no CBM. This is a clear sign of convergence in the old problem of comparing CBM between 1D and 3D models. Moreover, it is a significant confirmation of the potential of this novel approach towards developing 3D stellar evolutionary simulations.

3.4 Nucleosynthesis and evolution of the chemical composition

The simulations we present here have been produced employing a nuclear burning routine with an explicit list of isotopes to generate energy and drive convection. Making use of this routine, it is possible to study the time evolution of the chemical abundances and their distribution in the different layers of the simulations. We show in Fig. 11 the initial (dashed lines) and final (solid lines) mass fraction profiles from simulation r256e1 for the most important isotopes involved in neon burning. At the beginning of the simulation, the convective zone, identified by the central plateaus in the abundance profiles, is limited to the region between 4.5 and 5.8 × 10⁸ cm, while towards the end it has almost doubled in size, as we have already seen from Fig. 6. ²⁰Ne has been almost completely consumed in the convective zone via the reactions ²⁰Ne(γ, α)¹⁶O and ²⁰Ne(α, γ)²⁴Mg, while ¹⁶O and ²⁸Si have been produced as a result, as well as some ²⁴Mg that has been partially burnt to produce silicon according to ²⁴Mg(α, γ)²⁸Si.

Figure 11.

Angularly averaged mass fractions of ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si as a function of the stellar radius. Dashed: the abundances at the beginning of the r256e1 simulation. Solid: after ∼8 h, close to the end of the convective phase.

Another way of tracking the neon consumption is to look at the time evolution of the neon abundance inside the convective zone, as we show in Fig. 12. The plot shows the four simulations r256e1, r256e5, r256e10, and r256e50 with different boosting factors. During the initial transient phase (up to ∼100 s), some fluctuations in the neon abundance come from the initial propagation of plumes and eddies through the convection region, linked with the entrainment of some neon-rich material, and we can see that this trend is the same for all simulations except r256e50, where some neon is already being burnt due to the high energy boosting. After this phase, simulations consume neon on a different time-scale but all with a very similar trend: this is an additional confirmation of the fact that the boosting in luminosity affects mainly the simulation time-scale.

Figure 12.

Time evolution of the ²⁰Ne mass fraction in the bulk of the convective zone, for the four different simulations r256e1, r256e5, r256e10, and r256e50, each with a different boosting factor.

The chemical abundances can be also studied with a mean-field statistical analysis. We use here the Reynolds-averaged Navier–Stokes (RANS) framework developed for hydrodynamic simulations in spherical geometry by Mocák et al. (2014), making use of the dedicated open-source code ransx.¹ We refer to Mocák et al. (2014, 2018) for definitions and implementation. The RANS framework includes two types of averaging, a time averaging and an angular averaging. We will indicate here the Reynolds average (time and angular average) of a quantity q on a spherical shell at radius r as

$$\begin{eqnarray} \overline{q}(r)=\dfrac{1}{T\Delta \Omega }\int _0^T \int _{\Delta \Omega } q(r,\theta ,\phi ,t)\ \mathrm{d} t\ \mathrm{d}\Omega \end{eqnarray}$$

(8)

with dΩ = sin θdθdϕ the solid angle element, T the time window, and ΔΩ the solid angle of the shell. The Favre average (density-weighted average) is defined as |$\widetilde{q}=\overline{\rho q}/\overline{\rho }$|⁠; therefore, the field decomposition has been done as |$q=\overline{q}+q^{\prime }$| and |$q=\widetilde{q}+q^{\prime \prime }$|⁠, respectively, with |$\overline{q},\widetilde{q}$| the means and q′, q″ the fluctuations of the quantity q (see Mocák et al. 2014, 2018).

Fig. 13 shows the radial profiles of the mean turbulent flux in the r1024e1 simulation for ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si, defined as |$f_\text{i}=\overline{\rho }\ \widetilde{X^{\prime \prime }_\text{i} v^{\prime \prime }_\text{r}}$| for a species i. Positive and negative values in the flux represent upward and downward flows, respectively. It is evident that the flux is dominated by downward transport of ²⁰Ne, towards the bottom of the convective zone where the nuclear burning is taking place. On the other hand, ¹⁶O, ²⁴Mg, and ²⁸Si, the ashes of the burning, are all transported upwards in a similar way, through the entire convective zone. Additionally, a small non-zero flux is present immediately below the convective zone, slightly positive for silicon and negative for oxygen and magnesium. This is due to the mixing that takes place at the lower convective boundary, and the fact that silicon is more abundant below the boundary so it is brought inside the convective zone, while oxygen and magnesium are produced and more abundant above the boundary so they are transported downwards. Finally, the thin black line in Fig. 13 represents the sum of the flux profiles for all the 12 elements included in our nuclear network: it is always equal to zero, confirming that the sum of the mass fractions is conserved in our simulations.

Figure 13.

Turbulent flux profiles of ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si as a function of the stellar radius, from the r1024e1 simulation, averaged over the entrainment regime (three convective turnovers), and defined as |$f_\text{i}=\overline{\rho }\ \widetilde{X^{\prime \prime }_\text{i} v^{\prime \prime }_\text{r}}$|⁠. The thin black line is the sum of the flux profiles for all the 12 elements in the network.

Finally, we show in Fig. 14 the standard deviation profiles of the mass fraction for ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si, defined as |$\sigma _\text{i}=(\widetilde{X^{\prime \prime }_\text{i} X^{\prime \prime }_\text{i}})^{1/2}$| for a species i, and in Fig. 15 the standard deviation normalized by the Reynolds-averaged mass fraction of the isotopes |$\sigma _\text{i}/\overline{X}_\text{i}$|⁠, presenting deviations as fractions of the mean values.² The standard deviation represents the dispersion of the chemical composition as a function of the radius; therefore, it can be seen as a way of measuring the departure from a perfect spherical symmetry, providing important information for the comparison between 1D and multi-D stellar models. The largest deviations of up to a few tens of per cents (see Fig. 15) are found at the convective boundary locations. These deviations are explained by the deformation of the boundary due to the convective flow plowing into it as well as entrainment and the interaction of different layers at the interface (see e. g. Fig. 2). Inside the convective zone, on the other hand, the mixing makes the composition more homogeneous and reduces the dispersion down to a per cent or less. It is interesting to note that deviations are still present below the convective region due to the fluctuations generated by entrainment and internal gravity waves.

Figure 14.

Same as Fig. 13, but standard deviation profiles of the mass fraction for ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si, defined as |$\sigma _\text{i}=(\widetilde{X^{\prime \prime }_\text{i} X^{\prime \prime }_\text{i}})^{1/2}$|⁠.

Figure 15.

Same as Fig. 14, but normalized standard deviation profiles for ¹⁶O, ²⁰Ne, ²⁴Mg, and ²⁸Si, defined as |$\sigma _\text{i}/\overline{X}_\text{i}$|⁠.

Overall, the magnitude of the standard deviation is quite small and we do not expect major deviation from spherical symmetry for nucleosynthesis in normal convective burning episodes. The situation, however, is expected to be different in more dynamic contexts, such as merging shells or in cases where fuel is ingested in an unusual burning region (see Mocák et al. 2018; Andrassy et al. 2020).

4 CONCLUSIONS

In this paper, we have presented a set of 3D hydrodynamic simulations of a complete stellar burning phase, a neon-burning shell in a 20 M_⊙ star. The accuracy of the simulations has been enhanced by improvements in geometry and resolution, nuclear network and burning routine, and initial conditions. We show that results from our simulations may be analysed in terms of nucleosynthesis, studying the abundance evolution and stratification of the isotopes included in the nuclear network, and of hydrodynamics. For the latter, we have analysed the convective motions and tracked the convective boundary evolution, allowing observation of the growth of the convective zone and its death when fuel is exhausted. Studying CBM is also an excellent way of comparing 1D and multi-D stellar simulations, where results are often in disagreement. CBM in 1D models is subject to uncertainties and needs calibration; hydrodynamic models can provide this calibration, but only if started from correct initial conditions. This shows how the two approaches are mutually dependent, and a convergence of results can only be achieved by improving one with the other.

In previous works, 1D models sometimes include little CBM in the convective core but totally ignore any CBM in later phases, when convection is even stronger and its effects are more important. In other works, multi-D models are started from initial conditions with little to no CBM, always leading to a very strong entrainment that is completely in disagreement with the initial 1D model. With this work, we make a step forward towards the convergence of 3D to 1D stellar models (321D approach). Our 3D simulations of a burning shell, run continuously from early development to fuel exhaustion, show that (a) entrainment in late-phase shells does not proceed indefinitely as previously supposed, engulfing the entire star, but it halts when fuel is exhausted and convection dies (see Fig. 6), and (b) starting from initial conditions already including strong CBM, the resulting entrainment in 3D is much more in agreement with the 1D model (see Fig. 9). In particular for this last point, our entrainment study produced a law that may be equally well applied to CBM in convective cores (large RiB, vertical line in Fig. 9) and to late-phase shells (small RiB, data points in Fig. 9). This law may finally close the gap between the 1D and 3D stellar models, usually in disagreement regarding the amount of CBM to be included. Our results show that significant CBM is required not only in the convective core, as the most recent 1D models are starting to include, but also in the late-phase convective shells. The presence of large CBM is also supported by asteroseismic observations (e.g. Bowman 2020; Pedersen et al. 2021).

The work presented in this paper introduces exciting prospects for stellar modelling. We have shown that simulating an entire burning phase in more than one dimension is now possible using the right tools and enough computing resources. The next few years will inevitably feature more simulations of significant fractions of the stellar lifetime in multi-D. Since covering the entire stellar evolution will probably never be possible in more than one dimension, the 1D stellar model will remain the main tool for predicting and explaining stellar evolution. However, it is the interplay between 1D and 3D models that really pushes forwards our knowledge of stellar evolution, and we show here that an agreement in results between the two is possible.

We recall here that the range of applications of stellar modelling to other branches of astrophysics is large and various. This includes the production of accurate progenitor models as initial conditions for SN explosion studies (Müller & Janka 2015; Yoshida et al. 2019; Burrows & Vartanyan 2021) with possible deviations from spherical symmetry, potentially solving the long-standing core-collapse SN engine problems, but also comparison to asteroseismic measurements (Aerts 2021; Pedersen et al. 2021), analysis and implementation of magnetic fields and dynamo effects in stars (Varma & Müller 2021; Leidi et al. 2022), predictions on the nature of the different remnants (white dwarfs, neutron stars, and black holes) with improvements to the final–initial mass relation (Kaiser et al. 2020; Scott et al. 2021), nucleosynthesis, and galactic chemical evolution. Furthermore, the synergy of theoretical models and observations will help tackle today’s open problems of stellar astrophysics, such as the red supergiant problem (Smartt 2009) and the black hole mass gap (Woosley & Heger 2021).

SUPPORTING INFORMATION

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

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). WDA acknowledges support from the Theoretical Astrophysics Program (TAP) at the University of Arizona and Steward Observatory. CG has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant No. 833925). The University of Edinburgh is a charitable body, registered in Scotland, with Registration No. SC005336. CG, RH, and CM acknowledge ISSI, Bern, for its support in organizing collaboration. This article is based on work from the ChETEC COST Action (CA16117) and the European Union’s Horizon 2020 research and innovation programme (ChETEC-INFRA, Grant No. 101008324). The authors acknowledge 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).

DATA AVAILABILITY

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

Footnotes

1

https://github.com/mmicromegas/ransX

2

Note that in previous prompi studies, such as Mocák et al. (2018), the variable σ_i is used to indicate the variance, rather than the standard deviation.