Issue 
A&A
Volume 683, March 2024



Article Number  A97  
Number of page(s)  16  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202347407  
Published online  12 March 2024 
Towards a selfconsistent model of the convective core boundary in upper main sequence stars
I. 2.5D and 3D simulations
^{1}
Heidelberger Institut für Theoretische Studien, SchlossWolfsbrunnenweg 35, 69118 Heidelberg, Germany
email: robert.andrassy@hits.org
^{2}
HighPerformance Computing Center Stuttgart, Nobelstraße 19, 70569 Stuttgart, Germany
^{3}
Computer, Computational and Statistical Sciences (CCS) Division and Center for Theoretical Astrophysics (CTA), Los Alamos National Laboratory, Los Alamos, NM 87545, USA
^{4}
Zentrum für Astronomie der Universität Heidelberg, Astronomisches RechenInstitut, Mönchhofstr. 1214, 69120 Heidelberg, Germany
^{5}
Zentrum für Astronomie der Universität Heidelberg, Institut für Theoretische Astrophysik, Philosophenweg 12, 69120 Heidelberg, Germany
Received:
8
July
2023
Accepted:
14
December
2023
There is strong observational evidence that the convective cores of intermediatemass and massive main sequence stars are substantially larger than those predicted by standard stellarevolution models. However, it is unclear what physical processes cause this phenomenon or how to predict the extent and stratification of stellar convective boundary layers. Convective penetration is a thermaltimescale process that is likely to be particularly relevant during the slow evolution on the main sequence. We use our lowMachnumber SEVENLEAGUE HYDRO code to study this process in 2.5D and 3D geometries. Starting with a chemically homogeneous model of a 15 M_{⊙} zeroage main sequence star, we construct a series of simulations with the luminosity increased and opacity decreased by the same factor, ranging from 10^{3} to 10^{6}. After reaching thermal equilibrium, all of our models show a clear penetration layer; its thickness becomes statistically constant in time and it is shown to converge upon grid refinement. The penetration layer becomes nearly adiabatic with a steep transition to a radiative stratification in simulations at the lower end of our luminosity range. This structure corresponds to the adiabatic ‘step overshoot’ model often employed in stellarevolution calculations. The simulations with the highest and lowest luminosity differ by less than a factor of two in the penetration distance. The high computational cost of 3D simulations makes our current 3D data set rather sparse. Depending on how we extrapolate the 3D data to the actual luminosity of the initial stellar model, we obtain penetration distances ranging from 0.09 to 0.44 pressure scale heights, which is broadly compatible with observations.
Key words: convection / hydrodynamics / turbulence / methods: numerical / stars: interiors / stars: massive
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Numerous observations strongly suggest that the convective, hydrogenburning cores of intermediatemass and massive stars are larger than those predicted based on the linearstability 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 onedimensional (1D) stellar evolution theory resorts to simple parametric prescriptions for such processes, which severely limits the predictive power of current stellarevolution 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^{1}. 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^{−4}–10^{−3}) 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 stellarevolution models with the observations mentioned above. Threedimensional (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 massentrainment 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 highentropy 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:
where and 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 Φ_{0} 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.
Multidimensional 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 convectivepenetration problem. We eliminate compositional effects by using a realistic model of a 15 M_{⊙} zeroage 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, lowMach number SEVENLEAGUE 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. Onedimensional stellar model
We use version 15140 of the stellarevolution 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^{4} 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^{11} cm (4.98 R_{⊙}) and L_{⋆} = 7.49 × 10^{37} erg s^{−1} (1.95 × 10^{4} L_{⊙}), respectively. The model does not include any convective ‘overshoot’ or penetration. The latter evolves in our multidimensional models in a selfconsistent way. The Schwarzschild boundary of the core is located at 8.74 × 10^{10} cm (0.252 R_{⋆}). The MESA inlist for the model as well as the resulting profile file are available on Zenodo^{2}.
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 reintegrate 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 fullymixed 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
where r_{0} = 8.736 × 10^{10} cm is the radius of the initial core boundary. This profile is almost indistinguishable from that provided by MESA for r < r_{0}, see Fig. 1. There is a smallamplitude spike in the energy generation rate at r = r_{0} 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 cutout with a radius of 8.7 × 10^{9} 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^{4} yr, which corresponds to 6.5 × 10^{5} 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.
Time and length scales for the 15 M_{⊙} ZAMS star.
We solve these problems by increasing the energygeneration rate by a boost factor b ranging from 10^{3} to 10^{6}, 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:
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 , and a timeindependent 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^{9} cm (0.025 R_{⋆}) to 2.8 × 10^{11} 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 timestep 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^{3} 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^{3}. 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 radiativediffusion term is computed^{4}. 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 finitevolume 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^{6}. This ratio is of similar order as the costperstep ratio of the implicit to a comparable explicit stepper, making both approaches comparably expensive with b = 10^{6}. However, fast radiative diffusion caused by low density at the outer edge of the simulation domain imposes much stricter timestep constraints. The ratio of the lengths of implicit to diffusionlimited explicit time steps in the 2.5D simulation with b = 10^{6} and computed on the 1024 × 512 grid is ≈5000. Diffusion is much less constraining in simulations with b = 10^{3} 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^{5}.
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 wellbalancing 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 lowMachnumber 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^{−10} and 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. Timeaveraged 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 statisticalvariation 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^{4} 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 constanttemperature 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 = 10^{4} and N_{r} = 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^{−1}, 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^{6} down to 10^{3}. The velocity field in the core is dominated by largescale vortices as typical of 2D convection. Smallscale 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 ∈ (10^{6}, 10^{5}, 10^{4}, 10^{3}) performed on the 512 × 256 grid. The velocity is normalised using the typical convective velocity u_{conv} as given by the scaling law shown in Fig. 6. The dashed and solid lines give the timeaveraged 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 = 10^{5} performed on grids of 256 × 128 and 256 × 128^{2} 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 massweighted, rootmeansquare (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 u_{conv} 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 c_{ref} = 7.96 × 10^{7} cm s^{−1}. 
To shed more light on this issue, we show in Appendix A a set of 2.5D simulations with b = 10^{4} that include radiative diffusion but eliminate mass entrainment from the stably stratified envelope. These simulations develop smooth, largescale 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 highentropy 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 bestfit 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 , where σ_{series} is the standard deviation of the time series. We then compute 10^{5} 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^{6}. Our best estimate of convective velocity at the star’s nominal luminosity (i.e., b = 1) is (7.49 ± 0.13)×10^{4} cm s^{−1}. Taking the massweighted average sound speed c_{ref} = 7.96 × 10^{7} cm s^{−1} inside the initial Schwarzschild radius as a reference, we obtain a nominal Mach number of (9.41 ± 0.16)×10^{−4}. If we only fit the two lowestluminosity 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 3Dbased scaling gives a convective velocity of (6.65 ± 0.50)×10^{4} cm s^{−1} at b = 1, which is consistent with the 2.5Dbased prediction. All of the 2.5D data points are also statistically consistent with the 3Dbased scaling^{7}. 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^{−4} in a starlike 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 u_{conv} on the boost factor b. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed c_{ref} = 7.96 × 10^{7} cm s^{−1}. Statisticalvariation 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 MonteCarlo procedure described in Sect. 3.2. The grey shading shows the 3σ statisticalvariation interval around the mean fit. Parameters of the scaling law are given in the legend. The usual adiabatic scaling u_{conv} ∝ b^{1/3} is also shown for comparison. 
We define the convective turnover timescale as
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^{6} to 5.0 d at b = 10^{3}. 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^{4} 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 nearlyadiabatic core is much more pronounced. A convectivepenetration layer develops above the Schwarzschild boundary. The temperature gradient in the bulk of this layer is slightly subadiabatic but significantly superradiative, 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 latetime 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 = 10^{4} and N_{r} = 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 convectivepenetration 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 timeaveraged data. The time evolution of α is shown in Fig. 8 for four example simulations. The statistical variation in the simulations with b = 10^{6} 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 α = (r_{cb} − r_{Sb})/H_{p, Sb}, where r_{cb} is the radius of the convective boundary, r_{Sb} the radius of the Schwarzschild boundary, and H_{p, Sb} the pressure scale height at r_{Sb}. The four example simulations are 2.5D with N_{r} = 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 subadiabatic 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 N_{r} = 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 = 10^{3} 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 convectivepenetration 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^{8}. 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 stepovershoot 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^{6}) 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 secondorder code, we would asymptotically expect the dependence , where α_{∞} is the unknown value of α for N_{r} → ∞, and the constant C sets the magnitude of the error term. We observe firstorder convergence, , 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, , so that firstorder convergence corresponds to a straight line. The error bars indicate 3σ intervals of statistical variation. The grey shading shows σ, 2σ, and 3σ statisticalvariation 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^{3}, 10^{4}), possibly explaining the slowerthanexpected 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 secondorder accurate. However, our implementation of the constanttemperature boundary condition turns out to be only firstorder 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 oddeven 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 oddevendecoupling 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 firstorder convergence law to the data. This way, we obtain a welldefined, 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 statisticalvariation range in addition to the mean fit. We exclude the simulations with N_{r} = 128 () for the two lowest boost factors because these simulations seem to be severely underresolved 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.
Dimensionless penetration distance extrapolated to infinite grid resolution.
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 singlepowerlaw model.
Fig. 11. Penetration distance α extrapolated to an infinitely fine grid (see Fig. 10) as a function of the boost factor b. Assumed powerlaw dependencies α(b) are fitted to the 2.5D and 3D data points separately. The grey shading shows σ, 2σ, and 3σ statisticalvariation intervals of the fits. 
It is also possible that there is a powerlaw dependence α(b) but only for sufficiently low boost factors b, e.g., due to compressibility effects^{9} Indeed, the 2.5D data point with b = 10^{6} deviates from the scaling law discussed above and shown in Fig. 11 by 3.2σ. Therefore, we now exclude the b = 10^{6} 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^{5}. 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^{6} 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.5Dbased power law to make it pass through the 3Dbased data point at b = 10^{5}.
Fig. 12. As in Fig. 11, but excluding the simulations with b = 10^{6}, 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 = 10^{5}. 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 intermediatemass 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^{3} to 10^{6}. This scaling procedure preserves the value of the radiative temperature gradient. We run both 2.5D and 3D simulations using our timeimplicit, lowMachnumber 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 due to nuclear sources scales in proportion to b; as does the radiative flux because we scale the opacity in inverse proportion to b. The turbulent dissipation rate Φ_{0} is expected to scale with , 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^{3} ∝ 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 nonadiabatic 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 terminalage 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.
This prescription (see e.g., Paxton et al. 2011; Kippenhahn et al. 2012; or Anders & Pedersen 2023 for details) is physically inconsistent close to the convection zone, where the chemical diffusion coefficient is often assumed to be much larger than the coefficient of radiative diffusion but the exchange of heat associated with the rapid exchange of species between fluid elements is ignored. Some authors also assume the temperature gradient to be radiative when using the step overshoot prescription (see Anders & Pedersen 2023 for an overview).
The alternative of keeping the outgoing radiative flux constant becomes unstable when the opacity is assumed constant as is the case in our simulations. A largeenough negative temperature fluctuation at the boundary strongly reduces thermal conductivity χ ∝ T^{3}, blocking the flux from deeper layers. The constantflux boundary condition makes the temperature drop further, closing a positive feedback loop, which ultimately reduces the temperature to zero close to the boundary.
The ‘step overshoot’ prescription that can be activated in MESA via an inlist file only mixes composition in the ‘overshoot’ layer and the temperature gradient ∇ is kept radiative. However, it is possible to enforce ∇ = ∇_{ad} and model convective penetration by adding a simple subroutine to run_star_extras.f90 as we did in this example. Our MESA setup is available on Zenodo (https://doi.org/10.5281/zenodo.8127093).
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 LAUR2326456.
References
 Aerts, C. 2013, in EAS Publications Series, eds. K. Pavlovski, A. Tkachenko, & G. Torres, 64, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anders, E. H., & Pedersen, M. G. 2023, Galaxies, 11, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Anders, E. H., Jermyn, A. S., Lecoanet, D., & Brown, B. P. 2022, ApJ, 926, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Andrassy, R., Higl, J., Mao, H., et al. 2022, A&A, 659, A193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baraffe, I., Clarke, J., Morison, A., et al. 2023, MNRAS, 519, 5333 [NASA ADS] [CrossRef] [Google Scholar]
 Berberich, J. P., Chandrashekar, P., & Klingenberg, C. 2021, Comput. Fluids, 219, 104858 [Google Scholar]
 Blouin, S., Mao, H., Herwig, F., et al. 2023, MNRAS, 522, 1706 [CrossRef] [Google Scholar]
 Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A., & Torres, G. 2016, A&A, 592, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cristini, A., Hirschi, R., Meakin, C., et al. 2019, MNRAS, 484, 4645 [NASA ADS] [CrossRef] [Google Scholar]
 Demarque, P., Sarajedini, A., & Guo, X. J. 1994, ApJ, 426, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F. 2014, Ph.D. Thesis, Technische Universität München, Germany [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Edelmann, P. V. F., Horst, L., Berberich, J. P., et al. 2021, A&A, 652, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gilet, C., Almgren, A. S., Bell, J. B., et al. 2013, ApJ, 773, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Herwig, F., Woodward, P. R., Mao, H., et al. 2023, MNRAS, 525, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Horst, L., Hirschi, R., Edelmann, P. V. F., Andrássy, R., & Röpke, F. K. 2021, A&A, 653, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hosea, M., & Shampine, L. 1996, Appl. Numer. Math., 20, 21 [CrossRef] [Google Scholar]
 Hotta, H. 2017, ApJ, 843, 52 [Google Scholar]
 Jermyn, A. S., Anders, E. H., Lecoanet, D., & Cantiello, M. 2022, ApJ, 929, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, S., Andrassy, R., Sandalski, S., et al. 2017, MNRAS, 465, 2991 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J. 2019, A&A, 631, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (Berlin Heidelberg: SpringerVerlag) [Google Scholar]
 LeVeque, R. J. 2002, Finite Volume Methods for Hyperbolic Problems (Cambridge: Cambridge University Press), 31 [CrossRef] [Google Scholar]
 Lighthill, J. 2001, Waves in Fluids (Cambridge: Cambridge University Press) [Google Scholar]
 Liou, M.S. 2006, J. Comput. Phys., 214, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A., & Mermilliod, J. C. 1981, A&A, 93, 136 [NASA ADS] [Google Scholar]
 Mao, H., Woodward, P., Herwig, F., et al. 2023, ApJ, submitted [arXiv:2304.10470] [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Miczek, F. 2013, Dissertation, Technische Universität München, Germany [Google Scholar]
 Muthsam, H. J., Goeb, W., Kupka, F., Liebich, W., & Zoechling, J. 1995, A&A, 293, 127 [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Pratt, J., Baraffe, I., Goffrey, T., et al. 2020, A&A, 638, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roxburgh, I. W. 1965, MNRAS, 130, 223 [NASA ADS] [Google Scholar]
 Roxburgh, I. W. 1978, A&A, 65, 281 [NASA ADS] [Google Scholar]
 Roxburgh, I. W. 1989, A&A, 211, 361 [NASA ADS] [Google Scholar]
 Saslaw, W. C., & Schwarzschild, M. 1965, ApJ, 142, 1468 [Google Scholar]
 Sutherland, B. 2010, Internal Gravity Waves (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Woodward, P. R., Herwig, F., & Lin, P.H. 2015, ApJ, 798, 49 [Google Scholar]
 Zahn, J. P. 1991, A&A, 252, 179 [Google Scholar]
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 convectivepenetration 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^{4}. We keep the temperature at the outer radial boundary fixed at its initial value. The outer boundary is at r = 8.5003125 × 10^{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 coreonly 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^{4} 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 coreonly 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 smallscale 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 coreonly 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 highestresolution simulation when it was stopped and it also showed much less variability on short timescales (not visible in Fig. A.1) than the lowerresolution simulations. The velocities in all of the coreonly simulations are much larger than those in the analogous simulations of convective penetration.
Fig. A.1. Time evolution of the convective velocity u_{conv} 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 = 10^{4}. 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 c_{ref} = 7.96 × 10^{7} cm s^{−1}. Grid resolution is given in terms of the dimensionless radial grid spacing Δr such that a grid with N_{r} = 1024 in the ‘core + envelope’ case has Δr = 1. 
These effects can be understood as follows. Because the coreonly simulations lack a stable layer on top of the convective core, they cannot develop smallscale shear instabilities at that interface, leading to much more mundane convective flows. There is also no mass entrainment of highentropy 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 lowviscosity flows with a welldefined velocity in two spatial dimensions.
Appendix B: Accuracy of the SLH code
Motivated by the firstorder convergence of the penetration distance α upon grid refinement (see Sect. 3.4 and Fig. 10), we reevaluated 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_{1} errors in the distribution of an arbitrary quantity q to be
where q_{i, j} is the solution on a given grid and is the reference solution binned onto the same grid of N_{x} × N_{y} cells.^{10} The L_{1} 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 KelvinHelmoltz 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 secondorder accurate in this case. In the second test problem, we eliminate hydrodynamics and solve a diffusion problem with constanttemperature boundary conditions (Sect. B.2). This problem reveals a subtle effect that makes our original implementation of constanttemperature boundaries firstorder 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 secondorder 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 constanttemperature boundary condition is secondorder accurate also when combined with hydrodynamics. However, we find that errors related to the infamous oddevendecoupling effect limit the convergence rate to first order once relative errors have dropped below ≈10^{−2}. As discussed in Sect. B.3, this effect is intrinsic to collocated numerical schemes and it can only be suppressed, not eliminated.
B.1. KelvinHelmholtz instability
The initial condition for this test problem is
where ρ is the density, u and v are the components of the velocity vector, p is the pressure, and the function
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^{−5} ℛ, where ℛ is the gas constant per unit of mass. We stop the simulations at t = 80 when the instability has evolved well into its nonlinear 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 () 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 KelvinHelmholtz problem in terms of the specific internal energy e_{i} 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_{1} 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 secondorder is obtained on grids with at least 256 cells along the x axis. Our use of parabolic reconstruction, which is formally thirdorder accurate, is the most likely reason for the fast convergence.
Fig. B.2. Numerical convergence of the solutions to the KelvinHelmholtz problem. Relative L_{1} errors in the specific internal energy e_{i} are shown as a function of the grid resolution N_{x} along the x axis. A secondorder scaling is shown for reference. 
B.2. Heat diffusion with constanttemperature boundaries
To test our constanttemperature 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^{10} cm to 8 × 10^{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^{11} (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. Angleaveraged 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_{1} errors in the specific internal entropy s are shown in Fig. B.4. The constanttemperature 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. Reevaluating 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^{3}/ρ factor that enters the thermal conductivity at the cell interface. These approximations are secondorder 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 secondorder 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 firstorder error term and explaining the firstorder convergence observed in Fig. B.4. Because our simulations of convective penetration are long enough to form a statistically stationary state, firstorder 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 constanttemperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L_{1} errors in the specific entropy s are shown as a function of the radial grid resolution N_{r}. First and secondorder scaling laws are shown for reference. 
We improved our implementation of the constanttemperature boundary condition such that the state at the interface is obtained by extrapolating the state from the physical domain using secondorderaccurate expressions. We also compute the temperature gradient using onesided, secondorderaccurate 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 constanttemperature boundary conditions. The 2.5D simulation domain spans the radial range from 2 × 10^{10} cm to 8 × 10^{10} cm and polar angles 60° ≤ϑ ≤ 120°. The ‘hot bubble’ is modelled as an entropy perturbation
where s_{0} is the specific entropy of the background stratification, Δs the amplitude of the perturbation, r_{b} = 10^{10} cm the radius of the bubble, and ζ the distance from the centre of the bubble located at r = 3.5 × 10^{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^{3}, 10^{4}, 10^{5}). The perturbation amplitude Δs is 10^{5}, 4 × 10^{5}, and 1.6 × 10^{6} erg g^{−1} K^{−1} and the simulations are stopped after 1.6 × 10^{5}, 8 × 10^{4}, and 4 × 10^{4} s with boost factors of 10^{3}, 10^{4}, and 10^{5}, respectively. We keep the temperatures at both radial boundaries at their initial values using both implementations of the constanttemperature 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^{2} to 2048^{2} 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 hotbubble 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 L_{1} errors shown in Fig. B.6, are marked using the white dashes lines. 
We compute L_{1} errors in the specific entropy s in two subdomains, which we call ‘bubble’ and ‘boundary’ as depicted in Fig. B.5. The L_{1} errors for the two subdomains, three boost factors b, and two implementations of the constanttemperature 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_{1} errors for the ‘bubble’ subdomain are essentially the same with both kinds of boundary conditions and secondorder 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^{5}, which may be caused by the influence of the boundary effect discussed below. Indeed, the improved constanttemperature 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 hotbubble problem with the constanttemperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L_{1} errors in the specific entropy s are shown as a function of the radial grid resolution N_{r} for three boost factors b. First and secondorder 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^{−2}, a threshold almost independent of the boost factor. The reason for this behaviour is the infamous effect of oddeven 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 gridscale 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^{4} 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 hotbubble problem with b = 10^{4} 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 oddevendecoupling 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 pressurediffusion 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 fullymixed 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 = 10^{4} and N_{r} = 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 ∈ (10^{6}, 10^{5}, 10^{4}, 10^{3}) performed on the 512 × 256 grid. The velocity is normalised using the typical convective velocity u_{conv} as given by the scaling law shown in Fig. 6. The dashed and solid lines give the timeaveraged 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 = 10^{5} performed on grids of 256 × 128 and 256 × 128^{2} 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 u_{conv} 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 c_{ref} = 7.96 × 10^{7} cm s^{−1}. 

In the text 
Fig. 6. Dependence of the convective velocity u_{conv} on the boost factor b. The right vertical axis shows the corresponding Mach numbers computed using the reference sound speed c_{ref} = 7.96 × 10^{7} cm s^{−1}. Statisticalvariation 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 MonteCarlo procedure described in Sect. 3.2. The grey shading shows the 3σ statisticalvariation interval around the mean fit. Parameters of the scaling law are given in the legend. The usual adiabatic scaling u_{conv} ∝ b^{1/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 = 10^{4} and N_{r} = 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 convectivepenetration layer with dots marking the location of the steepest gradient in ∇/∇_{ad}. 

In the text 
Fig. 8. Time evolution of the dimensionless penetration distance α = (r_{cb} − r_{Sb})/H_{p, Sb}, where r_{cb} is the radius of the convective boundary, r_{Sb} the radius of the Schwarzschild boundary, and H_{p, Sb} the pressure scale height at r_{Sb}. The four example simulations are 2.5D with N_{r} = 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 N_{r} = 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 = 10^{3} 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 convectivepenetration 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, , so that firstorder convergence corresponds to a straight line. The error bars indicate 3σ intervals of statistical variation. The grey shading shows σ, 2σ, and 3σ statisticalvariation 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 powerlaw dependencies α(b) are fitted to the 2.5D and 3D data points separately. The grey shading shows σ, 2σ, and 3σ statisticalvariation intervals of the fits. 

In the text 
Fig. 12. As in Fig. 11, but excluding the simulations with b = 10^{6}, 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 = 10^{5}. 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 u_{conv} 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 = 10^{4}. 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 c_{ref} = 7.96 × 10^{7} cm s^{−1}. Grid resolution is given in terms of the dimensionless radial grid spacing Δr such that a grid with N_{r} = 1024 in the ‘core + envelope’ case has Δr = 1. 

In the text 
Fig. B.1. Reference solution to the 2D KelvinHelmholtz problem in terms of the specific internal energy e_{i} 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 KelvinHelmholtz problem. Relative L_{1} errors in the specific internal energy e_{i} are shown as a function of the grid resolution N_{x} along the x axis. A secondorder scaling is shown for reference. 

In the text 
Fig. B.3. Angleaveraged 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 constanttemperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L_{1} errors in the specific entropy s are shown as a function of the radial grid resolution N_{r}. First and secondorder scaling laws are shown for reference. 

In the text 
Fig. B.5. Reference solutions to the hotbubble 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 L_{1} 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 hotbubble problem with the constanttemperature boundary condition used in the simulations of convective penetration (old BC) and with an improved version of the same boundary condition (new BC). Relative L_{1} errors in the specific entropy s are shown as a function of the radial grid resolution N_{r} for three boost factors b. First and secondorder 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 hotbubble problem with b = 10^{4} 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 oddevendecoupling effect. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.