© The Authors 2024

1. Introduction

Numerous observations strongly suggest that the convective, hydrogen-burning cores of intermediate-mass and massive stars are larger than those predicted based on the linear-stability criteria of Schwarzschild and Ledoux (see e.g., Kippenhahn et al. 2012). The evidence has traditionally been based on colour–magnitude diagrams of open clusters (e.g., Maeder & Mermilliod 1981; Demarque et al. 1994) and observations of eclipsing binary systems (Claret & Torres 2016, and references therein). More recently, core sizes have also been measured using asteroseismology (Aerts 2013; Anders & Pedersen 2023, and references in the latter), confirming the large core radii. The size of the mixed core on the main sequence influences stellar lifetimes as well as stellar structure and evolution in later evolutionary stages. The umbrella term ‘convective overshooting’ has traditionally been used to describe all physical processes that contribute to extending stellar convection zones, although the terms ‘convective penetration’ and ‘convective boundary mixing’ are also used, often interchangeably. We reserve the term convective penetration for a process that is fast enough to change the thermal stratification beyond the Schwarzschild boundary (see also Anders & Pedersen 2023). Classical one-dimensional (1D) stellar evolution theory resorts to simple parametric prescriptions for such processes, which severely limits the predictive power of current stellar-evolution models. Two prescriptions are particularly popular (see e.g., Kippenhahn et al. 2012): ‘step overshoot’ and ‘exponentially decaying diffusion’. The first represents an approximate model of convective penetration. It assumes that mixing that is approximately as fast as that in the formally convective layer extends some distance beyond the Schwarzschild boundary, making the stratification chemically homogeneous and adiabatic with a discontinuous transition to the radiative stratification. The second prescription describes the mixing of chemical species using a diffusion coefficient decreasing exponentially with distance from the convective boundary while assuming that the thermal stratification remains radiative¹. Both prescriptions involve a free parameter that describes the radial extent of the mixing.

Core convection is characterised by low Mach numbers (of order 10⁻⁴–10⁻³) in main sequence stars. The overlying stratification is so stable that it stops the slow convective flows within a tiny fraction of the pressure scale height (Roxburgh 1965; Saslaw & Schwarzschild 1965). This kind of dynamical ‘overshoot’ cannot reconcile stellar-evolution models with the observations mentioned above. Three-dimensional (3D), hydrodynamic simulations show that convective flows turning around at the convective boundary can entrain relatively small mass fractions of material originally located beyond the formal boundary of convective instability. This process of mass entrainment occurs even at very stiff convective boundaries (e.g., Meakin & Arnett 2007; Woodward et al. 2015; Horst et al. 2021). However, simulations of core convection on the main sequence predict unrealistically high mass-entrainment rates (Meakin & Arnett 2007; Gilet et al. 2013; Herwig et al. 2023), suggesting that the observed fast entrainment is just a transient phenomenon. Indeed, most simulations of this kind do not include radiative diffusion, making it impossible to sustain the expected stellar structure with a convective core and a radiative envelope on long timescales. Particularly relevant is the thermal timescale on which heat transport processes, including convective penetration, set the thermal structure of the star.

Processes occurring on the thermal timescale can be qualitatively described in terms of how they affect the radial profile of entropy. Both nuclear burning and hydrodynamic entrainment of high-entropy material from the radiative envelope increase the mean entropy of the mixed core. On the other hand, radiative diffusion carries energy outwards, decreasing the core entropy. It is reasonable to assume that the two processes reach equilibrium and the core stops changing its size on the thermal timescale (ignoring slow changes in chemical composition occurring on the nuclear timescale). Somewhat surprisingly, analytical constraints can be placed on how large the convective core can be in the equilibrium state. Roxburgh (1978, 1989) averages the 3D equations of fluid motion and radiative heat transport in space and time and, assuming statistically stationary convection, derives a simple integral criterion for the size of the convective core:

$$\begin{array}{c}\hfill {\displaystyle {\int }_V(\overline{F}-\overline{\mathrm{\Gamma }})\frac{1}{T_0^2}\frac{\mathrm{d}{T}_0}{\mathrm{d}r}\mathrm{d}V={\int }_V\frac{{\mathrm{\Phi }}_0}{T_0}\mathrm{d}V,}\end{array}$$(1)

where $\overline{F}$ and $\overline{\mathrm{\Gamma }}$ are the mean radial components of the radiative flux and of the energy flux from nuclear sources, respectively, T is the temperature, and Φ the dissipation rate of kinetic energy. The subscript ‘0’ refers to the mean state and the integration is performed over volume V of the core. The radial profile of the dissipation rate Φ₀ depends on details of the turbulent convective flow. Neglecting this term, Roxburgh obtains the maximum possible core mass, concluding that it can be substantially larger than the mass inside the formal Schwarzschild boundary. Roxburgh (1978, 1989) as well as Zahn (1991) argue that the whole core must be nearly adiabatic, including its extension beyond the Schwarzschild boundary now known as the ‘convective penetration layer’. However, the thickness of this layer cannot be further constrained without constraining the turbulent dissipation rate.

Multi-dimensional hydrodynamic simulations can encompass all physical processes needed to compute detailed models of penetrative convection. However, the thermal timescale associated with typical stellar convective layers is orders of magnitude longer than the convective turnover timescale and the only way to to obtain thermally relaxed simulations is to boost the heat flux by several orders of magnitude. Using this technique, a few research groups recently managed to approach the thermal timescale and to observe the formation of a convective penetration layer in various stellar environments (Hotta 2017; Käpylä 2019; Baraffe et al. 2023; Blouin et al. 2023; Mao et al. 2023). Baraffe et al. (2023) ran 2D simulations of a range of main sequence stars with convective cores but only one of their simulations (of a 3 M_⊙ star) uses sufficient luminosity boosting for the penetration layer to start approaching thermal equilibrium. The growth of the convective core of a 25 M_⊙ main sequence star slows down on the thermal timescale in the 3D simulations of Mao et al. (2023) but the core does not stop growing. This might be a consequence of continued chemical mixing in the simulation. In the simpler case of a chemically homogeneous stratification, Anders et al. (2022) show that the penetration layer stops growing on a sufficiently long timescale, reaching a statistically stationary state. However, their simulations use a simplified stratification rather than a realistic stellar model and their use of the Boussinesq approximation eliminates any effects of compressibility, which may be important in the strongly stratified stellar interior.

Our study focuses on thermal aspects of the stellar convective-penetration problem. We eliminate compositional effects by using a realistic model of a 15 M_⊙ zero-age main sequence (ZAMS) star. Employing the standard approach of luminosity boosting, we show that a penetration layer forms at the core boundary and stops growing on the thermal timescale. We measure the core size in the thermally relaxed state and extrapolate towards the star’s nominal luminosity. We run both 2.5D (a sphere with assumed axial symmetry) and 3D simulations using our fully compressible, low-Mach number SEVEN-LEAGUE HYDRO (SLH) code (Miczek 2013; Edelmann 2014; Edelmann et al. 2021) to explore both the limit of low turbulent dissipation, as expected in 2D convection and assumed by Roxburgh, and the more realistic case of 3D turbulent convection.

We describe our 1D initial stellar model in Sects. 2.1 and 2.2. The numerical setup of the 2.5D and 3D SLH simulations is detailed in Sect. 2.3. We use a simple numerical experiment to illustrate the importance of radiative diffusion in Sect. 3.1. Relevant properties of our 2.5D and 3D simulations and their numerical convergence are described in Sects. 3.2–3.4. We then extrapolate the thickness of the penetration layer to the actual luminosity of the stellar model in Sect. 3.5. Our results are summarised and discussed in Sect. 4.

2. Methods

2.1. One-dimensional stellar model

We use version 15140 of the stellar-evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) to compute a model of a 15 M_⊙ star with a metallicity of Z = 0.02 at the ZAMS. The evolution model is stopped at the age of 7.96 × 10⁴ yr, when the convective core has fully developed and the star has reached thermal equilibrium. At this stage, the central hydrogen mass fraction has decreased by 0.156% from its initial value of 0.700. The star’s radius and luminosity are R_⋆ = 3.47 × 10¹¹ cm (4.98 R_⊙) and L_⋆ = 7.49 × 10³⁷ erg s⁻¹ (1.95 × 10⁴ L_⊙), respectively. The model does not include any convective ‘overshoot’ or penetration. The latter evolves in our multi-dimensional models in a self-consistent way. The Schwarzschild boundary of the core is located at 8.74 × 10¹⁰ cm (0.252 R_⋆). The MESA inlist for the model as well as the resulting profile file are available on Zenodo².

2.2. Simplification of the initial 1D model

The MESA model contains a 0.093% jump in the mean molecular weight at the core boundary because some hydrogen burning occurred during the initial thermal adjustment of the model. We remove this step by applying the central value of the mean molecular weight μ = 0.6174 to the whole star in order to eliminate any compositional effects. To obtain an initial condition for multidimensional simulations, we re-integrate the hydrostatic stratification assuming that the difference between the actual temperature gradient ∇ and the adiabatic gradient ∇_ad is given by the MESA model everywhere but in the convective core, where it is set to zero. The equation of state used in this integration process as well as in the hydrodynamic simulations described below accounts for a mixture of fully ionised, monatomic gas and radiation. The radial profile of gravitational acceleration is interpolated from the MESA model without any modification. The resulting profiles of the pressure p as well as specific entropy s closely match those of the original MESA model, see Fig. 1.

Fig. 1. Stratification of the pressure p, energy generation rate per unit volume ρεnuc, opacity κ, and specific entropy s in the MESA model (thick grey lines) and in the SLH model at t = 0 (thin coloured lines). The profiles of p and ρεnuc are normalised such that their maximum values are unity. The entropy of the convective core is set to 0. The jump in ρεnuc at the boundary of the fully-mixed core is caused by a slight change in chemical composition, which is neglected in the SLH model.

Figure 1 also shows the profiles of opacity κ and energy generation rate per unit volume ρε_nuc. The opacity is simply interpolated without any modification. The profile of the energy generation rate is for reasons of computational efficiency modelled using the fitting function

$$\begin{array}{c}\hfill \rho {\epsilon }_{\mathrm{nuc}}=\{\begin{array}{c}4.135\times {10}^{5}exp[-\frac{r^2}{{(3.15\times {10}^{10}\mathrm{cm})}^2}]{\mathrm{erg}\mathrm{cm}}^{-3}\mathrm{for}r<r_0,\hfill \\ 0\mathrm{otherwise},\hfill \end{array}\end{array}$$(2)

where r₀ = 8.736 × 10¹⁰ cm is the radius of the initial core boundary. This profile is almost indistinguishable from that provided by MESA for r < r₀, see Fig. 1. There is a small-amplitude spike in the energy generation rate at r = r₀ in the MESA model due to a slight change in composition. We do not model this spike. Additionally, our 2.5D and 3D simulations have a central cut-out with a radius of 8.7 × 10⁹ cm (0.025 R_⋆, see Sect. 2.3 for details). Taking into account all of these effects, the total luminosity of the modified model is 5.5% lower than that of the original MESA model.

2.3. 2.5D and 3D simulations

As mentioned in Sect. 1, the time and length scales involved in the process of convective penetration in stars make the problem prohibitively costly for multidimensional hydrodynamic simulations. We illustrate this in Table 1, which is based on results of our simulations (see Sect. 3 for details). The last row of the table is an extrapolation to our 15 M_⊙ star’s actual parameters. The thermal timescale of the model is 5.9 × 10⁴ yr, which corresponds to 6.5 × 10⁵ convective turnover timescales. This is well beyond what can currently be achieved using multidimensional simulations. Additionally, we must consider the need to resolve the length scale of radiative diffusion on the turnover timescale. Table 1 shows that this length scale evaluated at the core boundary is only 0.4 cell widths on a grid with 512 radial grid cells. To marginally resolve it by 4 cells, we would need a grid with ≈5000 radial cells, making the simulations extremely expensive.

We solve these problems by increasing the energy-generation rate by a boost factor b ranging from 10³ to 10⁶, which shortens both the thermal and convective turnover timescales. We decrease the opacity by the same factor, which increases the thermal diffusion length scale. This choice also has the advantage that the radiative temperature gradient ∇_rad, which is proportional to the product of the luminosity and opacity, becomes independent of b at a given temperature and density. The stratification of the radiative layer remains largely the same, although slight changes (≲25% in ∇/∇_ad) still occur near the outer boundary of the simulation box (r/R_⋆ ≳ 0.7) as the hydrostatic stratification adjusts in response to (1) an initial growth of the convective core (see Sect. 3.3 for details) and (2) the decreased luminosity as compared with the MESA model (by 5.5%, see Sect. 2.2 for details).

We use the SLH code to simulate the convective penetration process using the simplified 1D model as an initial condition. The code solves the following set of compressible, inviscid Euler equations with gravity and diffusive heat transport:

$$\begin{array}{cc}\hfill {\partial }_t\rho +\mathbf{\nabla }·(\rho \mathit{u})& =0,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\partial }_t(\rho \mathit{u})+\mathbf{\nabla }·(\rho \mathit{u}\otimes \mathit{u}+p\mathbb{I})& =\rho \mathit{g},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\partial }_t(\rho e_{\mathrm{t}})+\mathbf{\nabla }·[(\rho e_{\mathrm{t}}+p)\mathit{u}-\chi \mathbf{\nabla }T]& =\rho {\epsilon }_{\mathrm{nuc}},\hfill \end{array}$$(5)

where ρ, u, p, 𝕀, g, χ, T, ε_nuc denote the density, velocity, pressure, unit tensor, gravity, thermal conductivity, temperature, and energy generation rate by nuclear reactions, respectively. The specific total energy e_t = e_i + e_k + Ψ includes internal energy e_i(ρ, T), kinetic energy $e_{\mathrm{k}}=\frac{1}{2}{|\mathit{u}|}^2$, and a time-independent gravitational potential Ψ.

The equations are solved on two types of grids: 2.5D and 3D. Both of them match the spherical geometry of stars to suppress discretisation errors. The 3D grid is a spherical grid with uniform spacing in the radius r, polar angle (colatitude) ϑ, and azimuthal angle (longitude) φ. The 2.5D grid is a polar grid with uniform spacing in r and ϑ that describes a 3D spherical system with all variables constant in φ (rotational symmetry around the polar axis). This results in a geometric source term in the numerical scheme, which is taken into account. The grids cover the radial range from 8.7 × 10⁹ cm (0.025 R_⋆) to 2.8 × 10¹¹ cm (0.808 R_⋆). We limit the range of polar angles to 30° ≤ϑ ≤ 150° because the cell width in the azimuthal direction drops close to the polar axis, severely limiting the time-step length. In the 3D case, we include azimuthal angles in the range −60° ≤φ ≤ 60°.

To judge numerical convergence of our results, we run simulations with each boost factor b on a range of grids. The number of radial grid cells is N_r ∈ (128, 256, 512, 1024) in 2.5D and N_r ∈ (128, 192, 256) in 3D. The only exception is that we do not include a 2.5D simulation with b = 10³ and N_r = 1024, which would be too costly. The number of cells in each of the two angular directions is always N_r/2.

We impose reflective boundary conditions (see e.g., LeVeque 2002) at all boundaries of the computational domain. This is implemented via ghost cells filled such that the component of momentum normal to the boundary becomes an odd function and all other variables become even functions of the distance from the boundary³. The only exception is the outer radial boundary, where the temperature in the first ghost cell is set to the value given by the initial 1D model when and only when the radiative-diffusion term is computed⁴. The reflective boundaries in the angular directions eliminate horizontal shear flows, which often become strong with periodic boundaries.

The SLH code is based on a finite-volume discretisation and offers both explicit and implicit time steppers via the method of lines. In this work, we use the ESDIRK23 implicit stepper (Hosea & Shampine 1996). The bulk Mach numbers reported in Sect. 3.2 may suggest that the flows in some of our simulations are fast enough to be better suited for explicit time steppers. However, this impression is misleading. Even without considering radiative diffusion, the length of explicit time steps would be limited by the inner edge of the computational grid, where the sound speed is fastest and grid cells are narrowest in the angular directions. For this reason, implicit time steps can be 50–80 times longer than explicit ones even with the highest boost factor b = 10⁶. This ratio is of similar order as the cost-per-step ratio of the implicit to a comparable explicit stepper, making both approaches comparably expensive with b = 10⁶. However, fast radiative diffusion caused by low density at the outer edge of the simulation domain imposes much stricter time-step constraints. The ratio of the lengths of implicit to diffusion-limited explicit time steps in the 2.5D simulation with b = 10⁶ and computed on the 1024 × 512 grid is ≈5000. Diffusion is much less constraining in simulations with b = 10³ but convection is so slow in those that than longest possible explicit time steps would be ≈500 times shorter than implicit steps. All in all, the implicit approach saves a significant amount of computing time⁵.

We use unlimited parabolic reconstruction in space. Because the stratification spans five orders of magnitude in pressure (see Fig. 1), hydrostatic equilibrium requires a special treatment to suppress spurious flows caused by discretisation errors, see Edelmann et al. (2021) for details. Specifically, we use the well-balancing method of Berberich et al. (2021), later referred to as the ‘deviation method’ by Edelmann et al. (2021). We further increase our resolving power by employing the low-Mach-number numerical flux function AUSM+-up of Liou (2006), which is much less dissipative than standard Riemann solvers in slow flows. We slightly modify the flux function as described by Edelmann et al. (2021). The modified flux uses the parameters f_a = 10⁻¹⁰ and $f_a^p={10}^{-1}$ in the notation of the latter paper.

Although our numerical solver is deterministic, turbulent convection is well known to be highly sensitive to initial conditions and perturbations. A pair of simulations with initial conditions differing by small amounts or one with identical initial conditions but a small change to the implicit time stepper (e.g., allowing one more iteration per time step) are likely to produce completely different flow patterns after a few convective timescales. Time-averaged quantities become well defined when averaging over many convective timescales, assuming that there is a statistically stationary state. However, the statistical nature of turbulence must still be taken into account when comparing averages from different simulations. We take great care to quantify statistical variation our simulation data. We express statistical-variation ranges using ±1σ intervals, where σ is the uncorrected standard deviation of the statistical sample, unless mentioned otherwise.

3. Results

3.1. Importance of radiative diffusion

We start by demonstrating that the presence of radiative diffusion is essential for our simulations to reach a statistically stationary state. In Fig. 2, we compare two simulations performed with a boost factor of 10⁴ on a 256 × 128 grid. The convective core keeps growing in the simulation that does not include radiative diffusion until, ultimately, the whole model becomes convective. This effect can be understood even without considering hydrodynamic mass entrainment. The heat source in the core keeps increasing the mean entropy of the mixed core. Whenever that entropy becomes greater than the entropy of a stably stratified layer atop the core, that layer gets mixed into the core. Ultimately, the entropy of the whole model becomes approximately the same and the star becomes fully convective. This process is further accelerated by the presence of strong internal gravity waves, which seem to induce some mixing and flatten the entropy profile close to the outer boundary of the computational domain (i.e., where the waves reach their maximum amplitude). On the other hand, changes to the entropy profile are much more subtle in the simulation with radiative diffusion, in which the initial outward propagation of the convective boundary stops early on and the simulation reaches a statistically stationary state. This is possible because radiative diffusion is sufficiently fast to transport entropy generated in the core through the radiative envelope. The heat then leaves the box owing to the constant-temperature outer boundary condition. This statistically stationary state can be maintained indefinitely because we do not model changes to the chemical composition of the core due to nuclear burning. We only discuss simulations with radiative diffusion in the rest of this section.

Fig. 2. Time evolution of the radial profiles of the specific entropy s in 2.5D simulations without and with radiative diffusion. Both simulations have b = 104 and Nr = 256. The profiles are shifted such that the entropy of the convective core is 0 at t = 0. Simulation time is given as a fraction of the thermal timescale τth on the colour bar. Each curve, except for that showing the initial condition, corresponds to a time average over 0.05τth.

3.2. Velocity field

The velocity field depends primarily on the boost factor b. This factor sets both convective velocity and efficiency. Because we scale κ ∝ b⁻¹, the rate of radiative diffusion increases in proportion to b and fluid parcels exchange increasing amounts of heat with their surroundings as they rise and sink in the convective core. This effect makes convection less efficient at transporting heat. Figure 3 shows the velocity field in 2.5D simulations with boost factors ranging from 10⁶ down to 10³. The velocity field in the core is dominated by large-scale vortices as typical of 2D convection. Small-scale motions are further suppressed by radiative damping. The structure of the convective flow is completely different in 3D, see Fig. 4. Instead of large vortices, we obtain a turbulent cascade from large to small scales as expected for 3D convection (without strong rotation or magnetic fields).

Fig. 3. Velocity fields in 2.5D simulations with boost factors b ∈ (106, 105, 104, 103) performed on the 512 × 256 grid. The velocity is normalised using the typical convective velocity uconv as given by the scaling law shown in Fig. 6. The dashed and solid lines give the time-averaged radii of the Schwarzschild and convective boundaries, respectively. The simulations are shown at t = 0.5τth(b).

Fig. 4. As in Fig. 3, but comparing 2.5D and 3D simulations with a boost factor of b = 105 performed on grids of 256 × 128 and 256 × 1282 cells, respectively. In the 3D case, a slice with the spherical angle φ = 0 is shown.

The convective core generates internal gravity waves (IGWs), which propagate through the radiative envelope. When we decrease the boost factor b, the convection becomes slower and it generates waves at lower temporal frequencies. Basic theory of linear IGWs (see e.g., Lighthill 2001; Sutherland 2010) predicts that a decrease in frequency corresponds to a decrease in the radial wavelength for the same horizontal wavelength. This effect is evident in Fig. 3 and it is further amplified by the dependence on b of the radiative diffusivity, which filters out the shortest wavelengths.

The simulations in Fig. 3 are shown after a statistically stationary state has been reached. The figure also shows the radius of the convective boundary r_cb (solid lines), the exact definition of which we defer to the next section. For now, we only use it to measure convective velocity u_conv. We define this quantity to be the mass-weighted, root-mean-square (rms) velocity averaged between the inner boundary of the computational domain and the radius r_cb of the convective boundary. Figure 5 shows u_conv as a function of time, boost factor, grid resolution, and the dimensionality of the simulation. The convective velocity rapidly reaches a statistically stationary state. For each boost factor b, the velocity is not only statistically constant in time but it is also the same for 2.5D and 3D simulations and all grid resolutions. This is rather surprising since the structure of the convective flows differs substantially between 2.5D and 3D (Fig. 4). Although the same total energy flux has to be transported in both cases, the velocity is only one of several variables that contribute to the fluxes of enthalpy and kinetic energy. Indeed, numerical studies of stellar convection usually show faster convective velocities in 2D than in 3D with the same initial condition (Muthsam et al. 1995; Meakin & Arnett 2007; Pratt et al. 2020; Horst et al. 2021). It is possible that the growth of convective instability is in our simulations limited by the same phenomenon in both geometries, which could be radiative damping or buoyancy braking in the penetration layer.

Fig. 5. Time evolution of the convective velocity uconv in all of our simulations with radiative diffusion. The time axis shows simulation time as a fraction of the thermal timescale τth(b). Each data point corresponds to a time average over 0.05τth. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed cref = 7.96 × 107 cm s−1.

To shed more light on this issue, we show in Appendix A a set of 2.5D simulations with b = 10⁴ that include radiative diffusion but eliminate mass entrainment from the stably stratified envelope. These simulations develop smooth, large-scale convective cells much faster than those in the analogous simulations of convective penetration. Moreover, the velocity in the large cells does not converge upon grid refinement, which is understandable considering the fact that numerical viscosity is the only mechanism dissipating kinetic energy in these simulations. On the other hand, the convective velocity is numerically converged in the simulations of convective penetration, suggesting that the growth of kinetic energy in those simulations is limited by a resolved, physical effect. We suspect that the buoyancy of the high-entropy material entrained from the radiative envelope into the core plays an important role in determining the convective velocity.

Convective velocity scales with b^1/3 in simulations of adiabatic convection (Jones et al. 2017; Cristini et al. 2019; Edelmann et al. 2019; Horst et al. 2021; Herwig et al. 2023). Our best-fit scaling is u_conv ∝ b^{0.285 ± 0.002} for the four 2.5D simulations with N_r = 512, see Fig. 6. The mean exponent of the scaling differs by −26σ from the adiabatic value of 1/3. To obtain a measure of statistical variation for each data point, we consider the convective velocities averaged over each of the n_bins = 14 time bins after the initial transient (t > 0.3τ_th) as statistically independent measurements. The standard deviation of the average of all of the bins is then estimated to be $\sigma ={\sigma }_{\mathrm{series}}/n_{\mathrm{bins}}^{1/2}$, where σ_series is the standard deviation of the time series. We then compute 10⁵ statistical realisations of the scaling law assuming a normal distribution of statistical fluctuations, fit a power law to each of them, and compute the mean scaling and the statistical spread around the mean as a function of b. All 2.5D and 3D data points shown in Fig. 6 are statistically consistent with the scaling law⁶. Our best estimate of convective velocity at the star’s nominal luminosity (i.e., b = 1) is (7.49 ± 0.13)×10⁴ cm s⁻¹. Taking the mass-weighted average sound speed c_ref = 7.96 × 10⁷ cm s⁻¹ inside the initial Schwarzschild radius as a reference, we obtain a nominal Mach number of (9.41 ± 0.16)×10⁻⁴. If we only fit the two lowest-luminosity 2.5D data points the resulting scaling exponent is 0.291 ± 0.004, which is statistically consistent with the exponent based on the full 2.5D data set. Finally, if we only fit the two 3D data points we obtain a scaling exponent of 0.295 ± 0.006, which differs by −6.3σ from the adiabatic value of 1/3 and is statistically consistent with the exponent based on the full 2.5D data set. The 3D-based scaling gives a convective velocity of (6.65 ± 0.50)×10⁴ cm s⁻¹ at b = 1, which is consistent with the 2.5D-based prediction. All of the 2.5D data points are also statistically consistent with the 3D-based scaling⁷. It is possible that radiative damping makes the scaling shallower, although it is unclear why its exponent does not approach 1/3 at low boost factors, i.e., when radiative diffusion becomes relatively slow in the core. One might also suspect numerical effects. However, Edelmann et al. (2021) demonstrate that the SLH code produces the expected 1/3 velocity scaling for adiabatic convection down to Mach numbers of 2 × 10⁻⁴ in a star-like environment. Simulations of 3D turbulent convection with mass entrainment performed using the SLH code have also been shown to closely agree with those obtained using another four major hydrodynamic codes used in the field, see Andrassy et al. (2022). Finally, we provide a set of test simulations in Appendix B, which show that the results of SLH simulations converge upon grid refinement.

Fig. 6. Dependence of the convective velocity uconv on the boost factor b. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed cref = 7.96 × 107 cm s−1. Statistical-variation ranges for the individual data points are smaller than the markers and not shown. A straight line is fitted to the 2.5D data points using the Monte-Carlo procedure described in Sect. 3.2. The grey shading shows the 3σ statistical-variation interval around the mean fit. Parameters of the scaling law are given in the legend. The usual adiabatic scaling uconv ∝ b1/3 is also shown for comparison.

We define the convective turnover timescale as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{conv}}=\frac{2r_{\mathrm{cb}}}{u_{\mathrm{conv}}}.\end{array}$$(6)

For simplicity, we set r_cb = r_Sb, 0 + α(b)H_p, Sb, 0 in this estimate, where the r_Sb, 0 and H_p, Sb, 0 are the radius of and pressure scale height at the Schwarzschild boundary in the initial MESA model, and we take the scaling law α(b) derived from our 2.5D simulations in Sect. 3.5 and shown in Fig. 12. The resulting values of τ_conv, summarised in Table 1, range from 0.75 d at b = 10⁶ to 5.0 d at b = 10³. Extrapolating to the nominal luminosity (b = 1), we obtain τ_conv = 33 d.

3.3. Thermal evolution

In all of our simulations, we initially observe rapid growth in the size of the convective core. To illustrate this, we plot in Fig. 7 the time evolution of the actual and radiative temperature gradients ∇ and ∇_rad in the 2.5D simulation with b = 10⁴ and N_r = 1024. The gradients are normalised using the adiabatic temperature gradient ∇_ad, so the formal Schwarzschild boundary is located at the radius r_Sb where ∇_rad/∇_ad = 1. The Schwarzschild boundary moves slightly inwards over the thermal timescale. In contrast, the growth in the radial extent of the nearly-adiabatic core is much more pronounced. A convective-penetration layer develops above the Schwarzschild boundary. The temperature gradient in the bulk of this layer is slightly sub-adiabatic but significantly super-radiative, i.e., the total flux could in principle be transported by radiative diffusion but vigorous convection is still present. The outward propagation rate of the convective boundary drops rapidly after ≈0.1τ_th and a statistically stationary state seems to have been reached by ≈(0.2 − 0.3)τ_th. There is some slight late-time adjustment of the temperature gradient near the outer boundary of the simulation box (r/R_⋆ ≳ 0.7). However, the layer affected by it is sufficiently far out to not influence the penetration distance.

Fig. 7. Time evolution of the radial profiles of the actual and radiative temperature gradients ∇ and ∇rad, respectively, in the 2.5D simulation with b = 104 and Nr = 1024. The gradients are normalised using the adiabatic temperature gradient ∇ad. Simulation time is given as a fraction of the thermal timescale τth on the colour bar. Each curve, except for that showing the initial condition, corresponds to a time average over 0.05τth. The inset zooms in onto the convective-penetration layer with dots marking the location of the steepest gradient in ∇/∇ad.

To better quantify when the penetration distance stops changing, we compute the radii r_Sb and r_cb of the Schwarzschild and convective boundaries, respectively, using simulation data averaged over the spherical angles and over time bins 0.05 τ_th long. We define r_cb as the radius, where the drop in ∇(r)/∇_ad(r) is the steepest, see the inset in Fig. 7. The dimensionless penetration distance is then α = (r_cb − r_Sb)/H_p, Sb, where the pressure scale height H_p, Sb at the Schwarzschild boundary is derived from the same space- and time-averaged data. The time evolution of α is shown in Fig. 8 for four example simulations. The statistical variation in the simulations with b = 10⁶ is so large that all 20 data points seem to be consistent with a constant state. However, simulations with smaller boost factors cover many more convective timescales (see Table 1) and the statistical variation is largely suppressed in the time averages. These simulations show a clear initial increase in α followed by a statistically constant state. Based on Fig. 8, we define the first 0.3τ_th to be an initial transient for the purpose of measuring α. This transient is excluded from further analysis.

Fig. 8. Time evolution of the dimensionless penetration distance α = (rcb − rSb)/Hp, Sb, where rcb is the radius of the convective boundary, rSb the radius of the Schwarzschild boundary, and Hp, Sb the pressure scale height at rSb. The four example simulations are 2.5D with Nr = 512. The time axis shows the time as a fraction of the thermal timescale τth(b). The dotted vertical line at t = 0.3τth marks the end of the initial transient. Each data point corresponds to a time average over 0.05τth.

We show the radial profiles of the temperature gradient ∇ in 2.5D simulations with four different boost factors in Fig. 9. The stratification in the core becomes increasingly closer to being adiabatic as the boost factor is decreased. However, the stratification is at least slightly sub-adiabatic at r ≳ 0.17 R_⋆, which is well inside the Schwarzschild boundary at all four boost factors. The drop in the temperature gradient at the convective boundary (the outer boundary of the penetration layer) becomes increasingly steep with decreasing boost factor.

Fig. 9. Radial profiles of the temperature gradient ∇ in 2.5D simulations with Nr = 512 and the boost factors given in the legend. The gradients are normalised using the adiabatic temperature gradient ∇ad. The curves are derived from simulation data averaged over the time interval 0.3τth ≤ t ≤ τth except for the b = 103 simulation, for which we use the interval 0.3τth ≤ t ≤ 0.5τth. MESA models with α = 0.0 (our initial condition) and α = 0.6 are shown for comparison. The inset zooms in onto the convective-penetration layer with dots marking the location of the steepest gradient in ∇/∇ad.

Our simulations suggest that a simplified model, often referred to as ‘step overshoot’, which assumes that the penetration layer is fully mixed, perfectly adiabatic, and ends with a discontinuous jump to the radiative temperature gradient at the convective boundary is a good approximation to the thermal structure to be expected at the nominal luminosity (b = 1). A MESA model of a 15 M_⊙ ZAMS star computed using a parametric description of the penetration process is shown in Fig. 9 for comparison⁸. We use α = 0.6 in this model, which is approximately the value we obtain by extrapolating the results of our 2.5D simulation to the actual stellar luminosity (b = 1) in Sect. 3.5. In comparison, the popular alternative approach to model convective overshooting – exponentially decaying diffusion of chemical elements – assumes the temperature gradient to be radiative beyond the convective boundary. We do not show such a model in Fig. 9 because it would be indistinguishable from the model without any extra mixing at the convective boundary. It is possible if not likely that some slower form of mixing exists beyond the penetration layer. One could use the diffusive model beyond the layer covered by the adiabatic step-overshoot model, obtaining a combined model that Anders & Pedersen (2023) refer to as ‘extended convective penetration’.

3.4. Numerical convergence of the penetration distance

For each boost factor and grid dimensionality, we run simulations on a range of grids in order to quantify and suppress discretisation effects. Figure 10 shows that 2.5D simulations with the highest boost factor (b = 10⁶) are consistent with having the same value of α for all of our 2.5D grids (N_r = 128 to N_r = 1024), i.e., the simulations resolve all relevant scales well enough and α is numerically converged. The remaining sets of 2.5D and 3D simulations show some resolution dependence. For a second-order code, we would asymptotically expect the dependence $\alpha ={\alpha }_{\infty }+C{N}_{\mathrm{r}}^{-2}$, where α_∞ is the unknown value of α for N_r → ∞, and the constant C sets the magnitude of the error term. We observe first-order convergence, $\alpha ={\alpha }_{\infty }+C{N}_{\mathrm{r}}^{-1}$, instead, see Fig. 10.

Fig. 10. Numerical convergence of the penetration distance α for 2.5D and 3D simulations with different boost factors b. The radial grid spacing is expressed as the inverse number of radial grid cells, $N_{\mathrm{r}}^{-1}$, so that first-order convergence corresponds to a straight line. The error bars indicate 3σ intervals of statistical variation. The grey shading shows σ, 2σ, and 3σ statistical-variation intervals of the fits. The dashed lines show the mean fit in each case and the dotted lines extrapolate the fits towards data points excluded from the fitting procedure.

The diffusion length scales listed in Table 1 suggest that the smallest features one would expect to see in the entropy distribution around the core boundary are only marginally resolved on the finest of our grids in simulations with b ∈ (10³, 10⁴), possibly explaining the slower-than-expected numerical convergence. However, there are other numerical effects that affect the mean thermal stratification in our simulations and possibly also the penetration distance. We include an a posteriori study of these effects in Appendix B. In Appendix B.1, we show that the core numerical solver of SLH is second-order accurate. However, our implementation of the constant-temperature boundary condition turns out to be only first-order accurate when there is a density gradient inside the computational domain, see Appendix B.2 for details. Furthermore, we show in Appendix B.3 that the effect of odd-even decoupling starts to dominate the remaining numerical errors when the physical solution becomes sufficiently well resolved. This effect, which is always present in numerical schemes on collocated grids, also reduces the convergence rate in our test problem. Both the inaccurate boundary condition and the odd-even-decoupling effect may have influenced the penetration distance via their influence on the mean thermal stratification, although it is difficult to judge which of them is more important.

Nevertheless, Fig. 10 shows that the penetration distance does converge and we can fit the first-order convergence law to the data. This way, we obtain a well-defined, extrapolated penetration distance α_∞ that would correspond to a hypothetical, infinitely fine grid and is not affected by the aforementioned problems restricting the convergence rate to first order. We employ the Monte Carlo procedure described in Sect. 3.2, so we obtain a statistical-variation range in addition to the mean fit. We exclude the simulations with N_r = 128 ($N_{\mathrm{r}}^{-1}=7.8\times {10}^{-3}$) for the two lowest boost factors because these simulations seem to be severely under-resolved and they deviate significantly from the fits based on finer grids. None of the points included in the fits deviate from the fitting lines by more than 2σ. The resulting values of α_∞ are summarised in Table 2.

3.5. Luminosity dependence of the penetration distance

Having suppressed discretisation effects, we are left with a single estimate of the penetration distance and its range of statistical variation for each boost factor and grid dimensionality. We show these data points in Fig. 11. The 2.5D data points suggest a weak, although highly statistically significant, dependence of α on the boost factor b. We first blindly fit power laws to the 2.5D and 3D data sets individually using the Monte Carlo procedure described in Sect. 3.2. In the 2.5D case, the scaling exponent is only 0.043 ± 0.001 and the penetration distance extrapolated to the nominal luminosity (b = 1) is α_2.5D = 0.580 ± 0.005. If we assume that the α(b) dependence for 3D simulations can also be described by a single power law, we obtain a much steeper scaling with an exponent of 0.181 ± 0.011. After extrapolating over five orders of magnitude, the power law would then imply a penetration distance of only α_3D = 0.089 ± 0.012. However, our current set of 3D simulations on its own does not contain enough evidence for the single-power-law model.

Fig. 11. Penetration distance α extrapolated to an infinitely fine grid (see Fig. 10) as a function of the boost factor b. Assumed power-law dependencies α(b) are fitted to the 2.5D and 3D data points separately. The grey shading shows σ, 2σ, and 3σ statistical-variation intervals of the fits.

It is also possible that there is a power-law dependence α(b) but only for sufficiently low boost factors b, e.g., due to compressibility effects⁹ Indeed, the 2.5D data point with b = 10⁶ deviates from the scaling law discussed above and shown in Fig. 11 by 3.2σ. Therefore, we now exclude the b = 10⁶ data points and explore the alternative assumption that the power law’s slope is the same for 2.5D and 3D simulations with b ≤ 10⁵. Figure 12 shows that the scaling exponent based on the three remaining 2.5D data points remains essentially unchanged (0.042 ± 0.001) because the statistical weight of the b = 10⁶ data point considered in the fit in Fig. 11 is low. The largest deviation from the power law is only 0.5σ now. The penetration distance implied by the 2.5D simulations at the nominal luminosity does not change significantly either (α_2.5D = 0.585 ± 0.005). However, the much shallower scaling implies a much larger extrapolated penetration distance of α_3D = 0.443 ± 0.006 for the 3D simulations if we fix the scaling exponent and simply shift the 2.5D-based power law to make it pass through the 3D-based data point at b = 10⁵.

Fig. 12. As in Fig. 11, but excluding the simulations with b = 106, fitting a power law to the remaining 2.5D data points and, assuming that the scaling exponent is the same in 3D, shifting the power law to pass through the data point based on 3D simulations with b = 105. The statistical variation associated with that data point is added to that associated with the scaling law.

4. Summary and discussion

Observational evidence strongly suggests that the convective cores of intermediate-mass and massive, main sequence stars are substantially larger that what is predicted by linear stability analysis. The process of convective penetration, first described by Roxburgh (1978), seems to be able to enlarge the convective cores of main sequence stars enough to reconcile them with observations. However, the penetration distance depends on the complex physics of 3D turbulent convection in the core and such effects cannot be consistently captured in analytic treatments.

We performed a large set of simulations of convective penetration in a chemically homogeneous model of a 15 M_⊙ ZAMS star. We make it possible to cover the thermal timescale by increasing the luminosity and decreasing the opacity by the same large boost factor b ranging from 10³ to 10⁶. This scaling procedure preserves the value of the radiative temperature gradient. We run both 2.5D and 3D simulations using our time-implicit, low-Mach-number code SLH and quantify their numerical convergence using a range of computational grids.

We argue that mass entrainment into the convective core cannot stop in the adiabatic approximation and show an adiabatic simulation that demonstrates this. Indeed, the star becomes fully convective in the adiabatic simulation. On the other hand, simulations with radiative diffusion reach a statistically stationary state with a convective core and radiative envelope on the thermal timescale.

All of our simulations with radiative diffusion show a clear penetration layer, the thickness of which becomes statistically constant in time once thermal equilibrium has been reached. As expected, simulations performed with large luminosity boost factors show large deviations from adiabacity in the core. However, the core, including the bulk of the penetration layer, becomes nearly adiabatic as the luminosity is decreased. This result qualitatively agrees with the simulations of penetrative convection in similar stellar environments performed by Baraffe et al. (2023), Mao et al. (2023), and Blouin et al. (2023), as well as with the much more idealised simulations of Anders et al. (2022).

There is a steep jump from the adiabatic to the radiative temperature gradient at the outer boundary of the penetration layer. This means that the popular simplified model, often referred to as ‘step overshoot’, which assumes that the penetration layer is fully mixed, perfectly adiabatic, and ends with a discontinuous jump to the radiative temperature gradient at the convective boundary, is a good approximation of the thermal structure to be expected at the nominal luminosity (b = 1).

We find strong evidence for a weak dependence of the dimensionless penetration distance α on the boost factor b. Weak if any dependence of the core size on b is expected also from the theoretical point of view. In the Roxburgh criterion (Eq. (1)), the energy flux $\overline{\mathrm{\Gamma }}$ due to nuclear sources scales in proportion to b; as does the radiative flux $\overline{F}$ because we scale the opacity in inverse proportion to b. The turbulent dissipation rate Φ₀ is expected to scale with $u_{\text{conv}}^3/L_{\text{turb}}$, where L_turb is an integral length scale of the turbulence. If L_turb is fixed and we consider adiabatic convection, we have the usual scaling u³ ∝ b, and so both sides of Eq. (1) scale in proportion to b and the size of the core cannot depend on b. However, we show in Sect. 3.2 and Fig. 6 that the velocity scaling in our non-adiabatic simulations is slightly shallower than that of adiabatic convection. This changes the relative importance of the dissipation term as b is varied over several orders of magnitude. The slight change in the core radius r_cb itself may affect the dissipation rate given that it is reasonable to assume that L_turb ∝ r_cb. We must also keep in mind that Roxburgh (1989) assumes the velocity vector to vanish at the convective boundary and that contributions from terms quadratic in fluctuations can be neglected. Some of those terms may not be negligible in our simulations, especially because some energy flux is needed to drive the waves we clearly see in the radiative envelope.

Extrapolating the results of our 2.5D simulations to b = 1, we derive a rather large penetration distance of α = 0.585 ± 0.005. Our set of 3D simulations covers only two relatively large boost factors at the moment. For this reason, we show two possible extrapolations to b = 1 in Sect. 3.5, which give penetration distances α = 0.089 ± 0.012 and α = 0.443 ± 0.006. These two values differ by a factor of five, but both are broadly consistent with observations. For stars of similar masses, the compilation of (Anders & Pedersen 2023, their Fig. 12) of observational constraints on the penetration distance from asteroseismology shows measurements and upper limits that roughly span our two α estimates based on 3D simulations. Brott et al. (2011) discuss a drop in rotational velocity observed at a certain value of surface gravity in a sample of stars in the Large Magellanic Cloud. Interpreting the drop as the terminal-age main sequence, their calibration gives α = 0.34 ± 0.1 for stellar models in the mass range from approximately 10–20 M_⊙. The calibration of Claret & Torres (2016) based on eclipsing binaries shows that α increases from zero at 1.2 M_⊙ to ≈ 0.2 at 2.1 M_⊙, where it stops changing and remains approximately the same up to 4.4 M_⊙. The fitting formula given by Eqs. (11)–(16) of Jermyn et al. (2022), which is based on the simulations and theory of Anders et al. (2022), gives α = 0.19 for a 15 M_⊙ star.

The results presented here demonstrate only the most basic properties of the rich data set we have created. An analysis based on Reynolds averaging is already in progress and is expected to quantify differences between our 2.5D and 3D simulations and become a basis for constructing simplified models of the penetration process. We are also working on extending our simulations – especially those in 3D geometry – to lower luminosity boost factors. Although the simulations are based on a realistic stellar model, the model is chemically homogeneous and does not rotate or contain magnetic fields. Additionally, our computational grid currently covers a spherical wedge rather than the full sphere. All of these assumptions and simplifications should be gradually relaxed in order to generate predictions for a wide range of observed stars.

Acknowledgments

We are grateful to Raphael Hirschi and his group at Keele University as well as to Saskia Hekker at the Heidelberg Institute for Theoretical Studies for fruitful discussions and comments. We acknowledge support by the Klaus Tschira Foundation. This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 – 390900948 (the Heidelberg Structures Excellence Cluster). This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 945806). P.V.F.E. was supported by the U.S. Department of Energy through the Los Alamos National Laboratory (LANL). LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001). This work has been assigned a document release number LA-UR-23-26456.

References

Appendix A: 2.5D convection without mass entrainment

To judge the importance of mass entrainment on the convective velocity in 2.5D convection, we perform a simple numerical experiment: we run 2.5D simulations of core convection using the same setup as in the convective-penetration simulations except that we impose a reflective boundary condition in the outer parts of the initially isentropic core. This way, we eliminate any mass entrainment from the radiative envelope. The volume heating and opacity profile are the same as those in the penetration simulations with the boost factor b = 10⁴. We keep the temperature at the outer radial boundary fixed at its initial value. The outer boundary is at r = 8.5003125 × 10¹⁰ cm, so that we can achieve exactly the same grid spacing as in the simulations with the radiative envelope. We stop the set of core-only simulations after 0.3τ_th(b).

The elimination of the radiative envelope and of any mass entrainment leads to major qualitative and quantitative changes in the velocity field. The simulations with the radiative envelope and b = 10⁴ contain a large number of vortices of different sizes, the distribution and shapes of which change chaotically in time (see Fig. 3 for a representative snapshot). On the other hand, the core-only simulations develop a single, smooth convective cell filling the core early on, which later splits into two large cells. The relative sizes of those cells change on short timescales but the cells remain very smooth with essentially no small-scale structure.

Figure A.1 compares the time evolution of the convective velocity in the two sets of simulations. The velocity fluctuates around the same, nearly constant, function in all four simulations with the radiative envelope, independently of their grid resolution. On the other hand, we see a fast increase in the convective velocity in the core-only simulations. It ultimately reaches a constant value (after averaging out rapid changes as done in Fig. A.1), which increases upon grid refinement. Starting with the coarsest grid, the velocity on the next finer grid is larger by factors of 1.45, 1.51, and 1.21. The last ratio cannot be interpreted as a sign of numerical convergence because the convective velocity was still increasing in the highest-resolution simulation when it was stopped and it also showed much less variability on short timescales (not visible in Fig. A.1) than the lower-resolution simulations. The velocities in all of the core-only simulations are much larger than those in the analogous simulations of convective penetration.

Fig. A.1. Time evolution of the convective velocity uconv in two series of 2.5D simulations with radiative diffusion: the simulations of convective penetration with mass entrainment from the radiative envelope (‘core + envelope’) and similar simulations without the radiative envelope (‘core only’). All of the simulations have the boost factor b = 104. The time axis shows simulation time as a fraction of the thermal timescale τth(b). Each data point corresponds to a time average over 0.005τth. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed cref = 7.96 × 107 cm s−1. Grid resolution is given in terms of the dimensionless radial grid spacing Δr such that a grid with Nr = 1024 in the ‘core + envelope’ case has Δr = 1.

These effects can be understood as follows. Because the core-only simulations lack a stable layer on top of the convective core, they cannot develop small-scale shear instabilities at that interface, leading to much more mundane convective flows. There is also no mass entrainment of high-entropy material from the stable envelope and therefore no buoyancy braking associated with this phenomenon. The only remaining mechanism dissipating kinetic energy is numerical diffusion, the magnitude of which decreases upon grid refinement. Because the physical scale of the convective cells is always the same, the timescale of reaching a balance between the driving of convection and numerical dissipation becomes longer upon grid refinement and the flows accumulate more kinetic energy by the time a stationary state is reached. The fact that radiative diffusion decreases temperature and density fluctuations, hence also the buoyancy driving the convection, does not seem to be sufficient to obtain low-viscosity flows with a well-defined velocity in two spatial dimensions.

Appendix B: Accuracy of the SLH code

Motivated by the first-order convergence of the penetration distance α upon grid refinement (see Sect. 3.4 and Fig. 10), we re-evaluated the accuracy of the SLH code on a set of three test problems designed to isolate different aspects of the fluid dynamics simulations reported in this work. All three test problems are 2D and formulated such that one can obtain a ‘converged’ solution on a sufficiently fine grid and use it as a reference to quantify numerical errors affecting solutions on grids coarser than the reference grid. More specifically, we define L₁ errors in the distribution of an arbitrary quantity q to be

$$\begin{array}{c}\hfill L_1=\frac{1}{{\sigma }^{\mathrm{ref}}}\frac{{\sum }_{i=1}^{N_x}{\sum }_{j=1}^{N_y}|q_{i,j}-q_{i,j}^{\mathrm{ref}}|}{N_x{N}_y},\end{array}$$(B.1)

where q_i, j is the solution on a given grid and $q_{i,j}^{\mathrm{ref}}$ is the reference solution binned onto the same grid of N_x × N_y cells.¹⁰ The L₁ error is normalised using the standard deviation σ^ref of q in the reference solution because the fluctuations of interest are often much smaller than the mean value. The spatial and temporal discretisation in the same as in the simulations of convective penetration, see Sect. 2 for details.

In the first test problem, a Kelvin-Helmoltz instability with radiative diffusion (Sect. B.1), we exclude gravity and use periodic boundary conditions to eliminate any boundary effects. We demonstrate that the numerical scheme is second-order accurate in this case. In the second test problem, we eliminate hydrodynamics and solve a diffusion problem with constant-temperature boundary conditions (Sect. B.2). This problem reveals a subtle effect that makes our original implementation of constant-temperature boundaries first-order accurate when there is a density stratification inside the computational domain. We fix this issue and show that the new implementation, to be used in future simulations, is second-order accurate. Finally, we solve a problem involving the buoyant rise of a ‘hot bubble’ under the influence of radiative diffusion in the isentropic core of our 15 M_⊙ stellar model. This last problem demonstrates that the improved constant-temperature boundary condition is second-order accurate also when combined with hydrodynamics. However, we find that errors related to the infamous odd-even-decoupling effect limit the convergence rate to first order once relative errors have dropped below ≈10⁻². As discussed in Sect. B.3, this effect is intrinsic to collocated numerical schemes and it can only be suppressed, not eliminated.

B.1. Kelvin-Helmholtz instability

The initial condition for this test problem is

$$\begin{array}{cc}\hfill \rho & =\gamma [1.25-0.5\eta (y)],\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill u& =0.01[1-2\eta (y)],\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill v& =0.001sin(2\pi x),\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill p& =1,\hfill \end{array}$$(B.5)

where ρ is the density, u and v are the components of the velocity vector, p is the pressure, and the function

$$\begin{array}{c}\hfill \eta (y)=\{\begin{array}{cc}\frac{1}{2}\{1+sin[16\pi (y+0.25)]\},\hfill & \text{for}y>-\frac{9}{32}\text{and}y<-\frac{7}{32},\hfill \\ 1,\hfill & \text{for}y\ge -\frac{7}{32}\text{and}y\le \frac{7}{32},\hfill \\ \frac{1}{2}\{1-sin[16\pi (y-0.25)]\},\hfill & \text{for}y>\frac{7}{32}\text{and}y<\frac{9}{32},\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}\end{array}$$(B.6)

defines a smooth transition between the shearing layers. The smoothness of the transition makes it possible to obtain numerically converged solutions. The equation of state is that of an ideal monatomic gas with the ratio of specific heats γ = 1.4 and mean molecular weight μ = 1. The initial Mach number is approximately 0.01. We include radiative diffusion with a constant conductivity coefficient of χ = 3.61 × 10⁻⁵ ℛ, where ℛ is the gas constant per unit of mass. We stop the simulations at t = 80 when the instability has evolved well into its non-linear phase, see Fig. B.1. Our choice of χ is such that the diffusion affects the instability in a significant way without suppressing it completely. In the absence of hydrodynamics, the radiative diffusion on its own would widen the smooth transitions between the shear layers approximately by a factor of three by t = 80. Depending on whether we take the wavelength of the initial perturbation (unity) or the initial width of the transitions between the shear layers ($\frac{1}{16}$) as a reference length scale, the Péclet number of the flow is ≈ 870 and ≈ 55, respectively.

Fig. B.1. Reference solution to the 2D Kelvin-Helmholtz problem in terms of the specific internal energy ei normalised to range from 0 to 1. The solution was computed on a 4096 × 2048 grid.

We compute solutions to this test problem on equidistant, Cartesian grids ranging from 32 × 16 to 4096 × 2048 cells. The latter serves as the reference solution. The L₁ errors in the specific internal energy e_i are shown in Fig. B.2. As expected, the numerical convergence is slow on very coarse grids, which do not resolve important features of the solution. Convergence slightly faster than second-order is obtained on grids with at least 256 cells along the x axis. Our use of parabolic reconstruction, which is formally third-order accurate, is the most likely reason for the fast convergence.

Fig. B.2. Numerical convergence of the solutions to the Kelvin-Helmholtz problem. Relative L1 errors in the specific internal energy ei are shown as a function of the grid resolution Nx along the x axis. A second-order scaling is shown for reference.

B.2. Heat diffusion with constant-temperature boundaries

To test our constant-temperature boundary conditions, we set up a problem that involves diffusion but no hydrodynamics. We use part of the initially isentropic core of our 15 M_⊙ stellar model (Sect. 2.2), namely the radial range from 2 × 10¹⁰ cm to 8 × 10¹⁰ cm, as an initial condition. We fix the temperatures at both radial boundaries to their initial values, turn off volume heating, and include radiative diffusion with the original stellar opacity profile. We use 100 long, implicit time steps to reach 1% of the thermal timescale¹¹ (589 yr) and we stop the simulations at this point. The resulting change in the temperature gradient is shown in Fig. B.3.

Fig. B.3. Angle-averaged reference solution to the diffusion problem computed on a 2.5D grid of 2048 × 512 cells and visualised in terms of the temperature gradient. The isentropic initial state has ∇/∇ad = 1.

Although the stratification is spherically symmetric, we use 2.5D grids covering polar angles in the range 82.5° ≤ϑ ≤ 97.5°. We solve this problem on grids ranging from 16 × 4 to 2048 × 512 cells. The latter serves as the reference solution. The L₁ errors in the specific internal entropy s are shown in Fig. B.4. The constant-temperature boundary condition as implemented for the simulations of convective penetration reported in this work is denoted ‘old BC’ in the figure. Clearly, the solutions converge at first order only. Re-evaluating the implementation, we identified the following subtle issue that causes the slow convergence. When computing the diffusive flux at a given cell face, we use the states in the two cells around the interface to compute estimates of the temperature gradient, opacity (which is fixed), and the T³/ρ factor that enters the thermal conductivity at the cell interface. These approximations are second-order accurate for smooth solutions. At the domain boundary, we fix the temperature in the first ghost cell and use reflective boundary conditions for the density and momentum as described in Sect. 2.3. The temperature varies smoothly across the domain boundary, so the estimate of the temperature at the interface is second-order accurate. However, if there is a density stratification inside the computational domain the reflective boundary condition forces the component of the density gradient normal to the boundary to suddenly vanish at the domain boundary, introducing a first-order error term and explaining the first-order convergence observed in Fig. B.4. Because our simulations of convective penetration are long enough to form a statistically stationary state, first-order errors would have propagated from the outer domain boundary throughout the simulation domain, possibly reducing the numerical convergence rate of the mean thermal stratification and of the penetration distance.

Fig. B.4. Numerical convergence of the solutions to the diffusion problem with the constant-temperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L1 errors in the specific entropy s are shown as a function of the radial grid resolution Nr. First- and second-order scaling laws are shown for reference.

We improved our implementation of the constant-temperature boundary condition such that the state at the interface is obtained by extrapolating the state from the physical domain using second-order-accurate expressions. We also compute the temperature gradient using one-sided, second-order-accurate finite differences, fixing the temperature exactly at the domain boundary (rather than in the first ghost cell). The same test problem converges at second order with the new implementation, see Fig. B.4 (the ‘new BC’ curve). The improved implementation will be used in future simulations.

B.3. Buoyant rise of a hot bubble

Our third and most challenging test problem involves the buoyant rise of a single ‘hot bubble’ in the initially isentropic core of our 15 M_⊙ stellar model (Sect. 2.2) under the influence of radiative diffusion with constant-temperature boundary conditions. The 2.5D simulation domain spans the radial range from 2 × 10¹⁰ cm to 8 × 10¹⁰ cm and polar angles 60° ≤ϑ ≤ 120°. The ‘hot bubble’ is modelled as an entropy perturbation

$$\begin{array}{c}\hfill s=s_0+\{\begin{array}{c}\mathrm{\Delta }scos\left(\frac{\pi }{2}\frac{\zeta }{r_{\mathrm{b}}}\right)\mathrm{for}\zeta <r_{\mathrm{b}},\hfill \\ 0\mathrm{otherwise},\hfill \end{array}\end{array}$$(B.4)

where s₀ is the specific entropy of the background stratification, Δs the amplitude of the perturbation, r_b = 10¹⁰ cm the radius of the bubble, and ζ the distance from the centre of the bubble located at r = 3.5 × 10¹⁰ cm and ϑ = 0°. We do not include volume heating in the test problem. Radiative diffusion is included with the opacity profile identical to that used in our simulations of convective penetration with the boost factors of b ∈ (10³, 10⁴, 10⁵). The perturbation amplitude Δs is 10⁵, 4 × 10⁵, and 1.6 × 10⁶ erg g⁻¹ K⁻¹ and the simulations are stopped after 1.6 × 10⁵, 8 × 10⁴, and 4 × 10⁴ s with boost factors of 10³, 10⁴, and 10⁵, respectively. We keep the temperatures at both radial boundaries at their initial values using both implementations of the constant-temperature boundary condition discussed in Sect. B.2.

Buoyancy makes the bubble accelerate upwards. The central part of the bubble, where the entropy perturbation is largest, gradually overtakes the outer parts and the bubble turns into a pair of vortices. Radiative diffusion decreases the entropy contrast between the bubble and its surroundings. It also generates a mean entropy gradient in the background stratification with prominent boundary layers at the radial domain boundaries. We stop the simulations at a stage when the bubble has become strongly deformed but while it is still possible to obtain a numerically converged solution. We run this test problem on grids ranging from 16² to 2048² cells. The latter serves as the reference solution. The distributions of specific entropy at the end of the reference simulations are shown in Fig. B.5. The figure also includes the corresponding maximum Mach numbers reached in the simulation domain, which (by our choise of Δs) closely match the convective velocity in our simulations of convective penetration at the same boost factors, see Fig. 6.

Fig. B.5. Reference solutions to the hot-bubble problem with three different boost factors b in terms of the specific entropy s normalised to range from 0 to 1. The solutions were computed a 2.5D grid of 2048 × 2048 cells. The ‘bubble’ and ‘boundary’ subdomains, in which we compute the L1 errors shown in Fig. B.6, are marked using the white dashes lines.

We compute L₁ errors in the specific entropy s in two subdomains, which we call ‘bubble’ and ‘boundary’ as depicted in Fig. B.5. The L₁ errors for the two subdomains, three boost factors b, and two implementations of the constant-temperature boundary condition are shown in Fig. B.6. As expected, the convergence is slow on coarse grids and low boost factors, i.e. when important features of the solutions are unresolved. We now focus on finer grids (N_r ≳ 64). The L₁ errors for the ‘bubble’ subdomain are essentially the same with both kinds of boundary conditions and second-order convergence is obtained on sufficiently fine grids. The only exception is a slight decrease in the convergence rate on the finest grids in simulations with b = 10⁵, which may be caused by the influence of the boundary effect discussed below. Indeed, the improved constant-temperature boundary condition, detailed in Sect. B.2, gives slightly smaller errors for the ‘bubble’ subdomain with N_r ≳ 512.

Fig. B.6. Numerical convergence of the solutions to the hot-bubble problem with the constant-temperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L1 errors in the specific entropy s are shown as a function of the radial grid resolution Nr for three boost factors b. First- and second-order scaling laws are shown for reference.

The new boundary condition greatly improves the convergence rate in the boundary layer, in agreement with the results of the diffusion test reported in Sect. B.2. However, the convergence rate drops to (approximately) first order once the relative errors have dropped below ≈10⁻², a threshold almost independent of the boost factor. The reason for this behaviour is the infamous effect of odd-even decoupling that is unavoidable in numerical schemes working with collocated grids. The decoupling occurs when the mean stratification changes and it takes the form of grid-scale oscillations in the radial direction. We visualise this effect in Fig. B.7 using the signed residuals in s along the radial ray at ϑ = 97.5°, which passes through one of the vortices, in the simulation with b = 10⁴ computed on a 1024 × 1024 grid. Whereas the residuals are smooth in the radial range occupied by the bubble (c.f., Fig. B.5), they tend to oscillate between positive and negative values close to the radial domain boundaries.

Fig. B.7. Relative residuals in the specific entropy s along the radial ray at ϑ = 97.5° in a simulation of the hot-bubble problem with b = 104 computed on a 1024 × 1024 grid. The inset shows the outer boundary layer on a linear scale, which makes it easier to judge the magnitude of the odd-even-decoupling effect.

Similar oscillations occur close to the outer radial boundary in our simulations of convective penetration, possibly influencing the numerical convergence rate of the mean thermal stratification and of the penetration distance. The decoupling effect can be suppressed by increasing numerical diffusivity (via the pressure-diffusion term in the AUSM+-up flux function, see Liou (2006)). However, doing so only trades decoupling errors for diffusive errors, reducing the effective resolving power of the numerical scheme.

All Tables

All Figures

Fig. 1. Stratification of the pressure p, energy generation rate per unit volume ρεnuc, opacity κ, and specific entropy s in the MESA model (thick grey lines) and in the SLH model at t = 0 (thin coloured lines). The profiles of p and ρεnuc are normalised such that their maximum values are unity. The entropy of the convective core is set to 0. The jump in ρεnuc at the boundary of the fully-mixed core is caused by a slight change in chemical composition, which is neglected in the SLH model.

In the text

Fig. 2. Time evolution of the radial profiles of the specific entropy s in 2.5D simulations without and with radiative diffusion. Both simulations have b = 104 and Nr = 256. The profiles are shifted such that the entropy of the convective core is 0 at t = 0. Simulation time is given as a fraction of the thermal timescale τth on the colour bar. Each curve, except for that showing the initial condition, corresponds to a time average over 0.05τth.

In the text

	Fig. 3. Velocity fields in 2.5D simulations with boost factors b ∈ (106, 105, 104, 103) performed on the 512 × 256 grid. The velocity is normalised using the typical convective velocity uconv as given by the scaling law shown in Fig. 6. The dashed and solid lines give the time-averaged radii of the Schwarzschild and convective boundaries, respectively. The simulations are shown at t = 0.5τth(b).
In the text

	Fig. 4. As in Fig. 3, but comparing 2.5D and 3D simulations with a boost factor of b = 105 performed on grids of 256 × 128 and 256 × 1282 cells, respectively. In the 3D case, a slice with the spherical angle φ = 0 is shown.
In the text

	Fig. 5. Time evolution of the convective velocity uconv in all of our simulations with radiative diffusion. The time axis shows simulation time as a fraction of the thermal timescale τth(b). Each data point corresponds to a time average over 0.05τth. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed cref = 7.96 × 107 cm s−1.
In the text

Fig. 6. Dependence of the convective velocity uconv on the boost factor b. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed cref = 7.96 × 107 cm s−1. Statistical-variation ranges for the individual data points are smaller than the markers and not shown. A straight line is fitted to the 2.5D data points using the Monte-Carlo procedure described in Sect. 3.2. The grey shading shows the 3σ statistical-variation interval around the mean fit. Parameters of the scaling law are given in the legend. The usual adiabatic scaling uconv ∝ b1/3 is also shown for comparison.

In the text

Fig. 7. Time evolution of the radial profiles of the actual and radiative temperature gradients ∇ and ∇rad, respectively, in the 2.5D simulation with b = 104 and Nr = 1024. The gradients are normalised using the adiabatic temperature gradient ∇ad. Simulation time is given as a fraction of the thermal timescale τth on the colour bar. Each curve, except for that showing the initial condition, corresponds to a time average over 0.05τth. The inset zooms in onto the convective-penetration layer with dots marking the location of the steepest gradient in ∇/∇ad.

In the text

Fig. 8. Time evolution of the dimensionless penetration distance α = (rcb − rSb)/Hp, Sb, where rcb is the radius of the convective boundary, rSb the radius of the Schwarzschild boundary, and Hp, Sb the pressure scale height at rSb. The four example simulations are 2.5D with Nr = 512. The time axis shows the time as a fraction of the thermal timescale τth(b). The dotted vertical line at t = 0.3τth marks the end of the initial transient. Each data point corresponds to a time average over 0.05τth.

In the text

Fig. 9. Radial profiles of the temperature gradient ∇ in 2.5D simulations with Nr = 512 and the boost factors given in the legend. The gradients are normalised using the adiabatic temperature gradient ∇ad. The curves are derived from simulation data averaged over the time interval 0.3τth ≤ t ≤ τth except for the b = 103 simulation, for which we use the interval 0.3τth ≤ t ≤ 0.5τth. MESA models with α = 0.0 (our initial condition) and α = 0.6 are shown for comparison. The inset zooms in onto the convective-penetration layer with dots marking the location of the steepest gradient in ∇/∇ad.

In the text

Fig. 10. Numerical convergence of the penetration distance α for 2.5D and 3D simulations with different boost factors b. The radial grid spacing is expressed as the inverse number of radial grid cells, $N_{\mathrm{r}}^{-1}$, so that first-order convergence corresponds to a straight line. The error bars indicate 3σ intervals of statistical variation. The grey shading shows σ, 2σ, and 3σ statistical-variation intervals of the fits. The dashed lines show the mean fit in each case and the dotted lines extrapolate the fits towards data points excluded from the fitting procedure.

In the text

	Fig. 11. Penetration distance α extrapolated to an infinitely fine grid (see Fig. 10) as a function of the boost factor b. Assumed power-law dependencies α(b) are fitted to the 2.5D and 3D data points separately. The grey shading shows σ, 2σ, and 3σ statistical-variation intervals of the fits.
In the text

	Fig. 12. As in Fig. 11, but excluding the simulations with b = 106, fitting a power law to the remaining 2.5D data points and, assuming that the scaling exponent is the same in 3D, shifting the power law to pass through the data point based on 3D simulations with b = 105. The statistical variation associated with that data point is added to that associated with the scaling law.
In the text

Fig. A.1. Time evolution of the convective velocity uconv in two series of 2.5D simulations with radiative diffusion: the simulations of convective penetration with mass entrainment from the radiative envelope (‘core + envelope’) and similar simulations without the radiative envelope (‘core only’). All of the simulations have the boost factor b = 104. The time axis shows simulation time as a fraction of the thermal timescale τth(b). Each data point corresponds to a time average over 0.005τth. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed cref = 7.96 × 107 cm s−1. Grid resolution is given in terms of the dimensionless radial grid spacing Δr such that a grid with Nr = 1024 in the ‘core + envelope’ case has Δr = 1.

In the text

	Fig. B.1. Reference solution to the 2D Kelvin-Helmholtz problem in terms of the specific internal energy ei normalised to range from 0 to 1. The solution was computed on a 4096 × 2048 grid.
In the text

	Fig. B.2. Numerical convergence of the solutions to the Kelvin-Helmholtz problem. Relative L1 errors in the specific internal energy ei are shown as a function of the grid resolution Nx along the x axis. A second-order scaling is shown for reference.
In the text

	Fig. B.3. Angle-averaged reference solution to the diffusion problem computed on a 2.5D grid of 2048 × 512 cells and visualised in terms of the temperature gradient. The isentropic initial state has ∇/∇ad = 1.
In the text

Fig. B.4. Numerical convergence of the solutions to the diffusion problem with the constant-temperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L1 errors in the specific entropy s are shown as a function of the radial grid resolution Nr. First- and second-order scaling laws are shown for reference.

In the text

	Fig. B.5. Reference solutions to the hot-bubble problem with three different boost factors b in terms of the specific entropy s normalised to range from 0 to 1. The solutions were computed a 2.5D grid of 2048 × 2048 cells. The ‘bubble’ and ‘boundary’ subdomains, in which we compute the L1 errors shown in Fig. B.6, are marked using the white dashes lines.
In the text

Fig. B.6. Numerical convergence of the solutions to the hot-bubble problem with the constant-temperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L1 errors in the specific entropy s are shown as a function of the radial grid resolution Nr for three boost factors b. First- and second-order scaling laws are shown for reference.

In the text

	Fig. B.7. Relative residuals in the specific entropy s along the radial ray at ϑ = 97.5° in a simulation of the hot-bubble problem with b = 104 computed on a 1024 × 1024 grid. The inset shows the outer boundary layer on a linear scale, which makes it easier to judge the magnitude of the odd-even-decoupling effect.
In the text