© F. Ahlborn et al. 2022

1. Introduction

Internal transport processes substantially shape the structure of stars. This concerns the transport of energy as well as the transport of chemical elements and angular momentum. In intermediate and high mass main-sequence stars, the energy released by nuclear fusion in the centre is transported by convection, which is the transport of energy by means of fluid motion. In deep convective zones, convection determines the temperature gradient and dominates the chemical mixing processes due to its efficiency. Hence, convection is a crucial aspect of stellar structure and evolution models.

Notwithstanding its importance, convection remains one of the major uncertainties in stellar structure and evolution modelling. The most commonly used description of convection for stellar models is the so-called mixing-length theory (MLT; Biermann 1932; Böhm-Vitense 1958). Despite its simplicity – MLT is a local and time-independent theory – it has been used very successfully over many years to model stellar interiors. With ever-improving observational facilities and methods, however, more and more deficits of MLT have become apparent. One of the main problems of MLT concerns the treatment of convective boundaries. When ignoring compositional effects, the acceleration drops to zero at the so-called Schwarzschild boundary, while the velocity generally does not. From a theoretical point of view, it has been shown that convective motions should pass the Schwarzschild boundary and penetrate into the stable layers (Roxburgh 1978, 1992; Zahn 1991). In stable layers, convective motions are braked and the material carried with the flow mixes with the surroundings. In the literature, different terms are used to describe the processes at convective boundaries. Zahn (1991) differentiated between thermally efficient and inefficient convection. Thermally efficient convection is able to modify the model temperature gradient and is therefore referred to as sub-adiabatic penetration. Thermally inefficient convection leaves the temperature gradient unchanged while still mixing chemical elements. Zahn (1991) referred to this as overshoot mixing. Finally, the notion of convective entrainment refers to the continuous ingestion of material at convective boundaries into a convective region (Turner 1986). This effect is found to represent the convective boundary mixing in 3D simulations of an oxygen-burning shell (Meakin & Arnett 2007).

In the MLT picture, the convective velocities drop to zero at the formal Schwarzschild boundary, preventing convective mixing beyond this point. In a physical configuration, however, only the acceleration of fluid elements disappears, while the velocity generally remains finite. To include the effects of convective overshooting in stellar models, additional descriptions have to be applied. Early attempts were proposed by, for example, Saslaw & Schwarzschild (1965) and Shaviv & Salpeter (1973). For a critical review of these theories, we refer the reader to Renzini (1987). The main effect of convective overshooting on the stellar structure can be mimicked by introducing additional mixing at convective boundaries during stellar evolution, which we refer to as ad hoc overshooting. To account for the required mixing, Freytag et al. (1996) introduced an additional diffusion constant, which decreases exponentially as a function of the distance to the Schwarzschild boundary. This is also known as ‘diffusive overshooting’. Although this approach was originally based on 2D simulations of envelopes in A-type main-sequence stars and DA-type white dwarfs that have thin convective zones subject to strong radiative losses, it is commonly applied to all convective boundaries in stellar evolution models.

When assuming instantaneous mixing, ad hoc overshooting is commonly introduced by extending the chemical composition of the convective region by a certain fraction of the pressure scale-height into the formally stable region. This approach is also known as ‘step overshooting’. Neither of these descriptions provides constraints on the temperature gradient in the overshooting region or on the extent of this region.

Due to the high Reynolds numbers, convection is a highly turbulent process. The dynamic timescales of the involved flows are many orders of magnitude shorter than the nuclear timescales of stellar evolution in most evolutionary phases. This poses serious problems for numerical descriptions of convection. Full 3D hydrodynamic simulations are computationally expensive and can cover physical time spans of the order of years only. In contrast, stellar evolution calculations usually need to be computed over at least a couple of million years. At the same time, the computational costs need to be low. Therefore, the direct inclusion of 3D hydrodynamics into stellar evolution calculations is not feasible on present-day computers. This means that the effects of convection need to be included in stellar evolution models by some other means. One way to include the results of 3D simulations in 1D stellar evolution codes is via Reynolds-averaged Navier-Stokes analysis, as outlined, for example, by Viallet et al. (2013) and Arnett et al. (2015).

However, since 3D simulations suffer from the aforementioned computing limitations, in practice this method is essentially a variant of the Reynolds stress approach (Keller & Friedmann 1925; Chou 1945; André et al. 1976; Canuto 1992) and, more generally speaking, of turbulent convection models (TCMs). The main idea of TCMs is to construct higher order moment equations from the hydrodynamic equations and reduce the dimensionality of the problem by averaging over two spatial directions. These equations describe the dynamics of convection. One main difference of TCMs compared to MLT is the occurrence of transport terms. These terms describe the transport of physical quantities, for example the turbulent kinetic energy (TKE), by means of convective flows and naturally lead to the emergence of phenomena such as convective overshooting without any ad hoc description. As these transport terms connect different layers, they are often termed ‘non-local’. Furthermore, TCMs provide the convective flux that allows the temperature gradient to also be computed in the overshooting region. To date, a number of TCMs have been developed for stellar astrophysics (e.g., Xiong 1978, 1986; Stellingwerf 1982; Kuhfuß 1986, 1987; Canuto 1992, 1993, 1997, 2011; Canuto & Dubovikov 1998; Li & Yang 2007). Due to the non-linearity of the Navier-Stokes equation, higher order terms appear in the equations of the TCMs that cannot be computed consistently within the set of equations. To compute the higher order moments, so-called closure relations need to be applied. These closure relations are one of the major sources of uncertainty for any TCM. Ultimately, these closures could be supplied by 3D simulations (Chan & Sofia 1989, 1996; Kupka et al. 1999b, 2007a,b,c, 2008; Kupka 2007; Viallet et al. 2013; Arnett et al. 2015).

In this work we present the results of a TCM applied in a stellar evolution code. We implemented the TCM by Kuhfuß (1987) into the GARching STellar Evolution Code (GARSTEC; Weiss & Schlattl 2008). The additional partial differential equations of the TCM are solved simultaneously with the stellar structure equations via the commonly used Henyey method (Flaskamp 2003). Using this implementation, we computed stellar evolution models of low and intermediate mass main-sequence stars from 1.5 to 8 M_⊙ and studied their convective cores. For the present paper, surface convection zones are computed with conventional MLT, but models using the Kuhfuß TCM for envelope convection as well are already underway. We demonstrate that the original description of the Kuhfußconvection model leads to erroneous convective properties. In Kupka et al. (2022, hereafter Paper I), we show that the dissipation rate of the kinetic energy plays an important role in the description of turbulent convection. We implemented the description of Paper I, and here we further show that it leads to physically reasonable properties of convection in the framework of the Kuhfuß convection model, which we briefly introduce and review in Sect. 2. Its application in stellar models is the subject of the following section. In Sect. 4 we compute stellar models in a mass range from 1.5 M_⊙ to 8 M_⊙ with the new model and compare its results with previously used models. Our discussion in Sect. 5 analyses the origin of the structural differences between the new 3-equation model and the standard 1-equation model of Kuhfuß (1987) and reviews constraints on core sizes obtained by other available methods. In our conclusions, in Sect. 6, we provide a summary of the advantages and limitations of the new model and an outlook on further developments.

2. Implementation of the Kuhfuß (1987) model

We implemented the Kuhfuß (1987) convection model as described in Appendix A of Paper I into GARSTEC. This includes the three partial differential equations for the TKE, ω, the convective flux, Π, and the entropy fluctuations, Φ, as well as the increased dissipation rate in the overshooting zones. We refer to this model as the 3-equation model. For completeness, we here repeat the final model equations:

$$\begin{array}{cc}& {\text{d}}_t\omega =\frac{{\mathrm{\nabla }}_{\mathrm{ad}}T}{H_{\mathrm{p}}}\mathrm{\Pi }-\frac{C_D}{\mathrm{\Lambda }}{\omega }^{3/2}-{\mathcal{F}}_{\omega }\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\text{d}}_t\mathrm{\Pi }=\frac{2{\mathrm{\nabla }}_{\mathrm{ad}}T}{H_{\mathrm{p}}}\mathrm{\Phi }+\frac{2c_{\mathrm{p}}}{3H_{\mathrm{p}}}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\mathrm{ad}})\omega -{\mathcal{F}}_{\mathrm{\Pi }}-\frac{1}{{\tau }_{\text{rad}}}\mathrm{\Pi }\hfill \end{array}$$(2)

$$\begin{array}{cc}& {\text{d}}_t\mathrm{\Phi }=\frac{c_{\mathrm{p}}}{H_{\mathrm{p}}}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\mathrm{ad}})\mathrm{\Pi }-{\mathcal{F}}_{\mathrm{\Phi }}-\frac{2}{{\tau }_{\text{rad}}}\mathrm{\Phi },\hfill \end{array}$$(3)

where ∇ and ∇_ad refer to the actual and adiabatic temperature gradient, respectively. The substantial derivative is defined as ${\mathrm{d}}_t={\partial }_t+\overline{\mathit{v}}·\mathrm{\nabla }$. The radiative dissipation timescale is given as

$$\begin{array}{c}\hfill {\tau }_{\text{rad}}=\frac{c_{\mathrm{p}}\kappa {\rho }^2{\mathrm{\Lambda }}^2}{4\sigma T^3{\gamma }_R^2}.\end{array}$$

Non-local fluxes ℱ_a are modelled as

$$\begin{array}{cc}\hfill {\mathcal{F}}_a& =\frac{1}{\overline{\rho }}\text{div}(-{\alpha }_a\overline{\rho }\mathrm{\Lambda }\sqrt{\omega }\mathrm{\nabla }\overline{a})\hfill \end{array}$$

for a = ω, Π, Φ. The symbols C_D, γ_R, α_ω, α_Π and α_Φ are parameters. Kuhfuß (1987) suggested a value of $C_D=8/3·\sqrt{2/3}$ and ${\gamma }_R=2\sqrt{3}$ to be compatible with MLT in the flavour of Böhm-Vitense (1958) in the local limit of the model. The length scale of TKE dissipation is denoted by Λ and is discussed in detail below. Furthermore, c_p refers to the specific heat capacity at constant pressure, κ to Rosseland opacity and σ to the Stefan-Boltzmann-constant. The variables T, ρ, and H_p are temperature, density, and pressure scale height, as usual in stellar structure models. For the details of the derivation and the definitions of the symbols, we refer to the appendix A in Paper I and to Kuhfuß (1987) and Flaskamp (2003).

2.1. The 1-equation model

In addition to the 3-equation model, Kuhfuß (1987) also suggested a simplified version of his TCM (see also Kuhfuß 1986). The number of equations is reduced to one by introducing the following approximation for the convective flux:

$$\begin{array}{c}\hfill \mathrm{\Pi }=-{\alpha }_{\mathrm{s}}\mathrm{\Lambda }\sqrt{\omega }\frac{c_{\mathrm{p}}}{H_{\mathrm{p}}}(\mathrm{\nabla }-{\mathrm{\nabla }}_{\mathrm{ad}}).\end{array}$$(4)

The parameter ${\alpha }_{\mathrm{s}}=(1/2)\sqrt{2/3}$ has been calibrated to MLT and Λ is again the length scale of TKE dissipation. As for the dissipation parameter C_D, the parameter value of α_s is obtained by calibrating convective velocities and fluxes of the local 1-equation model to MLT analytically. This approximation of the convective flux allows the convective flux from the ω equation to be eliminated such that it is only necessary to solve a single equation. We refer to this model as the 1-equation model. When applying the 1-equation model, the length scale for the dissipation of the TKE is defined as Λ = αH_p – with a freely adjustable parameter α – throughout the rest of the paper, equivalent to the usual mixing length. Additionally, we applied the reduction of the mixing length towards the centre following Wuchterl (1995) to counteract the divergence of the pressure scale height in this region (see Eq. (6) with β_s = 1). For a value of order unity for α results are indeed comparable to MLT results.

2.2. Convection equations

As most modern stellar evolution codes GARSTEC makes use of the implicit Henyey scheme to solve the four stellar structure equations (Henyey et al. 1964, 1965; Kippenhahn et al. 1967). The equations describing convection by MLT are solved algebraically outside the four stellar structure equations. To incorporate the three equations describing the convection model, Flaskamp (2003) extended the Henyey-scheme of GARSTEC to solve for in total seven variables (four stellar structure variables plus three convection variables). A solution for both the stellar structure and the convective variables is found by iterating over all variables simultaneously. The coefficients of the convection model depend on the stellar structure variables, such that the behaviour of the convection model is strongly coupled to the stellar structure. On the other hand, the stellar structure is coupled to the convective variables through the temperature gradient and the chemical composition. As described in Paper I the temperature gradient of the stellar model is computed self-consistently from the convective flux in each iteration (Eq. (A.7) in Paper I). The chemical mixing in convective zones is computed in the framework of a diffusion equation, alongside the composition changes due to nuclear burning after the structure equations have been solved. The diffusion constant is computed from the TKE determined by the convection model. Following Langer et al. (1985) the diffusion coefficient is computed as

$$\begin{array}{c}\hfill D={\alpha }_{\mathrm{s}}\mathrm{\Lambda }\sqrt{\omega ,}\end{array}$$(5)

where we chose the same parameter α_s as for the diffusion coefficient of entropy in Eq. (4).

The equations can also describe time-dependent effects. This was demonstrated for example by Flaskamp (2003), when computing models through the core helium flash at the tip of the red-giant branch, by Wuchterl & Feuchtinger (1998) in an application to protostars and non-linear pulsations, and by Feuchtinger (1999) for RR Lyrae stars. In this work, however, we focus on main-sequence stars that evolve on the nuclear timescale of hydrogen burning. This means that structural changes are sufficiently slow to neglect time-dependent terms and immediately solve for the stationary solution of the convection equations (left-hand sides of the TCM equations Eqs. (1)–(3) are set to zero). By iterating for the stationary solution, the code searches for the converged stellar structure and convection variables for a given chemical composition. When non-local effects are included in the convection model, the stationary solution describes the overshooting zone self-consistently. Its extent and temperature gradient are only constrained by the convection model, without any external descriptions.

As described in Paper I and at the beginning of this section, the Kuhfuß (1987) convection model contains a number of parameters. The values for these parameters need to be set. As described above, C_D and γ_R are obtained by calibrating a local model to MLT. The parameter values for the non-local terms α_ω, α_Π and α_Φ cannot be calibrated to MLT as they describe intrinsically non-local effects. Kuhfuß (1987) suggested a default value of α_ω = 0.25 by comparing kinetic energy and dissipation in a ballistic picture. No default values were given for the parameters α_Π and α_Φ in the non-local case. Although both the 1- and 3-equation models still contain a number of parameters, they are advantageous compared to for example MLT because the parameters describe physically more fundamental properties of the theory. For example, the parameter α_ω describes the impact of the non-local flux of the TKE, which is responsible for the extent of the overshooting region. However, compared to diffusive overshooting or step overshooting α_ω does not set the actual length scale of the overshooting. The extent of the overshooting is determined self-consistently from the solution of the TCM equations. We further investigate the impact of these other parameters below.

2.3. Dissipation rate

In Paper I we show that the original description of the dissipation rate proposed by Kuhfuß (1987) leads to an excessively large overshooting region. Therefore, in Paper I the dissipation rate was increased by taking into account buoyancy waves as a sink for the TKE. The increase in the dissipation rate was realised through a modification of its associated dissipation length scale. The dissipation rate is inversely proportional to this length scale, ϵ = c_ϵω^3/2/Λ, such that a decrease in the latter leads to an increase in the dissipation rate. This modification of the TKE dissipation length scale was implemented through a harmonic sum:

$$\begin{array}{c}\hfill \frac{1}{\mathrm{\Lambda }}=\frac{1}{\alpha H_{\mathrm{p}}}+\frac{1}{{\beta }_{\mathrm{s}}r,}\end{array}$$(6)

where the newly introduced parameter β_s is defined as

$$\begin{array}{c}\hfill {\beta }_{\mathrm{s}}={(1+{\lambda }_{\mathrm{s}}\stackrel{\sim }{N})}^{-1}\end{array}$$

and

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{s}}=c_4\mathrm{\Lambda }{\omega }^{-1/2},\end{array}$$(7)

where c₄ = c₃/(c₂c_ϵ) (Eq. (21) in Paper I). The parameters c₂ = 1.92 and c₃ = 0.3 are model parameters from Canuto & Dubovikov (1998) and c_ϵ is the dissipation parameter of the convection model. The buoyancy frequency $\stackrel{\sim }{N}$ is computed according to

$$\begin{array}{c}\hfill {\stackrel{\sim }{N}}^2=\frac{g^2\rho }{p}({\mathrm{\nabla }}_{\text{ad}}-\mathrm{\nabla }+{\mathrm{\nabla }}_{\mu }),\end{array}$$

assuming an ideal gas law. Here, ∇_μ indicates the dimensionless mean molecular weight gradient and g refers to the gravitational acceleration. Close to the stellar centre the pressure scale height and in turn the dissipation length scale, Λ, if defined through the pressure scale height, diverge. To avoid this divergence, Flaskamp (2003) introduced the modification by Wuchterl (1995) in the Kuhfuß model. This extension of the model we also implemented equivalently to Eq. (6), but with a constant parameter β_s = 1. In unstably stratified regions, the new correction factor is not applied and c₄ = 0 as c₃ drops to 0 (Canuto & Dubovikov 1998). Hence, the harmonic sum Eq. (6) recovers the Wuchterl (1995) expression automatically. We note that at this point where ∇ − ∇_ad → 0, also $\stackrel{\sim }{N}\to 0$ and β_s smoothly transitions to 1, such that setting c₃ = 0 does not introduce any discontinuity. The harmonic sum Eq. (6) can be converted into an equation for the dissipation length Λ. Rewriting Eq. (6) yields

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\underset{<1}{\underbrace{\frac{r}{r+\alpha H_{\mathrm{p}}\frac{1}{{\beta }_{\mathrm{s}}}}}}\alpha H_{\mathrm{p}},\end{array}$$(8)

which immediately shows that the derived expression is in essence a reduction factor for the mixing length. Plugging in the definitions for β_s and λ_s one finds the following expression:

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{r}{r+\alpha H_{\mathrm{p}}(1+c_4\mathrm{\Lambda }{\omega }^{-1/2}\stackrel{\sim }{N})}\alpha H_{\mathrm{p}}.\end{array}$$

This is a quadratic equation in Λ that can be solved to obtain the reduced length scale. Here, we expressed the dissipation rate timescale τ = 2ω/ϵ in terms of ω and Λ by noting that ϵ = c_ϵω^3/2/Λ and τ = 2Λ/(c_ϵω^1/2), which allowed us to rewrite an expression proportional to the ratio of TKE to buoyancy timescales, τ/τ_b, into one proportional to $\mathrm{\Lambda }{\omega }^{-1/2}\stackrel{\sim }{N}$ (we refer to Paper I for details). The final model of Λ in the convection zone reads

$$\begin{array}{c}\hfill \mathrm{\Lambda }(r)=-\frac{r+\alpha H_{\mathrm{p}}}{2c_4\stackrel{\sim }{N}{\omega }^{-1/2}\alpha H_{\mathrm{p}}}+\sqrt{{\left[\frac{r+\alpha H_{\mathrm{p}}}{2c_4\stackrel{\sim }{N}{\omega }^{-1/2}\alpha H_{\mathrm{p}}}\right]}^2+\frac{r}{c_4\stackrel{\sim }{N}{\omega }^{-1/2}}}\end{array}$$(9)

for ∇ < ∇_ad, and

$$\begin{array}{c}\hfill \mathrm{\Lambda }(r)={(\frac{1}{\alpha H_{\mathrm{p}}}+\frac{1}{{\beta }_{\mathrm{c}}r})}^{-1},\end{array}$$(10)

for ∇ > ∇_ad, where β_c = 1. To obtain a physically reasonable, positive dissipation length scale, Λ, the plus sign in front of the square root in the solution of the quadratic equation has to be chosen.

We note that the parameter c_ϵ takes different values in the Canuto and Kuhfuß convection models. Here, we take c_ϵ = C_D, which is the dissipation parameter in the Kuhfuß model. In unstably stratified regions with c₄ = 0 Eq. (6) solves explicitly for Λ. In stably stratified regions, the parameter takes a value of c₄ ≈ 0.072 using the parameters c₂ and c₃ from Canuto & Dubovikov (1998) and c_ϵ = C_D.

As mentioned already in Sect. 3.6 in Paper I, the effect of changing c_ϵ on changing c₄ to some extent cancels out. We discuss this further in Appendix B.

3. Stellar models

We used GARSTEC to compute stellar models in a mass range of 1.5–8 M_⊙ in the core hydrogen-burning phase. We used the OPAL equation of state, OPAL opacities (Iglesias & Rogers 1996), extended by low temperature opacities by J. Ferguson (private communication and Ferguson et al. 2005), both for the Grevesse & Noels (1993, GN93) mixture of heavy elements. For the initial mass fractions we chose X = 0.7 and Z = 0.02 for all models. Convective chemical mixing was described in a diffusive way. We used MLT plus diffusive overshooting as described by Freytag et al. (1996) to evolve the models through an initial equilibration to the beginning of the main-sequence phase. Then we first switched to the 1- and subsequently to the 3-equation model to generate a starting model for the computation with the 3-equation model. As we are interested in the effects of overshoot mixing, all models shown in the following were computed, including the non-local terms in the 1- and 3-equation version of the Kuhfuß theory.

For the diffusive overshooting, we used the default GARSTEC parameter value of f_OV = 0.02, which was calibrated by fitting GARSTEC isochrones to the colour-magnitude diagrams of open clusters (Magic et al. 2010). The diffusion coefficient is computed according to Eq. (C.1). For small convective cores, excessively large overshooting zones can occur when applying the ad hoc overshooting schemes due to the diverging pressure scale height in the centre. To avoid such unfavourable conditions, a reduced overshooting parameter value is determined in GARSTEC, by applying a geometrical cutoff depending on the comparison between the radial extent of the convective region and the scale-height at its border. For a brief discussion of different geometric cutoff descriptions, we refer to Appendix C. For the results in this section, the ‘tanh’ cutoff according to Eq. (C.3) was used. For the parameter α_ω we chose a value of 0.3 as this value results in a similar convective core size as the ad hoc overshooting models for f_OV = 0.02 in the 5 M_⊙ model, which is in the middle of our mass range. For the parameters α_Π and α_Φ we chose the same value, assuming that the non-local transport behaves similar for all convective variables. We would like to point out that α_ω is an adjustable parameter and an external calibration will be necessary. This is discussed further below.

3.1. The 3-equation model

As a representative example, we first discuss the evolution and the internal structure of a 5 M_⊙ model applying the 3-equation, non-local convection theory. Figure 1 shows the profiles of the TKE variable ω (panel a), the convective flux variable Π (panel b), the super-adiabatic gradient ∇ − ∇_ad (panel c), where ∇ is the temperature gradient resulting from the convection model, and the diffusion coefficient according to Eq. (5) on a logarithmic scale (panel d). The TKE clearly extends beyond the Schwarzschild boundary, which we computed as usual as the point where ∇_ad = ∇_rad. Such a behaviour cannot be observed in MLT models, as it is a direct result of the non-local terms in the Kuhfuß convection model. The associated diffusion coefficient shows a rather high value beyond the Schwarzschild boundary and throughout the overshooting zone. This increases the size of the mixed convective core and therefore naturally creates an overshooting zone. Given the profile of the diffusion coefficient, the chemical mixing resembles the step overshooting rather than the diffusive overshooting scheme. The convective flux variable shows a region of negative flux beyond the Schwarzschild boundary. In the lower panel of Fig. 2 the region of negative convective flux is shown enlarged in the inset (likewise for the convective flux variable Π in the middle panel of Fig. 1). This is due to the braking of the convective motions in the stable, radiative stratification. The extent and magnitude of the negative convective flux are very comparable to early results from Xiong (1986), who found that the convective flux penetrates less deeply into the stable layers than for example the kinetic energy and has a nearly negligible magnitude compared to the total flux.

Fig. 1. Summary of the interior structure of the convective core of a 5 M⊙ main-sequence model, calculated with the 3-equation model. The dashed black line indicates the Schwarzschild boundary. The different panels show: (a) the TKE, (b) the convective flux variable, (c) the super-adiabatic temperature gradient, and (d) the diffusion coefficient according to Eq. (5). The selected stellar model has a central hydrogen abundance of Xc = 0.6. The inset in panel b shows the region of negative convective flux just beyond the Schwarzschild boundary.

Fig. 2. Comparison of models (yellow lines) using MLT without overshooting to those computed with the 3-equation theory (blue lines) and the 1-equation model (red lines). Panel a: compares the hydrogen profiles at an early stage on the main sequence, when Xc = 0.6 (solid lines), and at the end of it (Xc ≈ 0; dashed lines). Panel b: convective fluxes of an MLT model with diffusive overshooting (dotted black line), a 3-equation model (solid blue line), and a 1-equation model (solid red line). These models have been selected to have the same central hydrogen abundance of Xc = 0.6 and the same chemically homogeneous core size.

Finally, panel c shows the super-adiabatic temperature gradient of this stellar model (we note the vertical scale). In the inner part of the convection zone the model shows a very small super-adiabatic gradient as it is expected for regions with convective driving. At a fractional mass of about 0.05 the temperature gradient drops below the adiabatic value. In contrast to local models, this point does not coincide with the formal Schwarzschild boundary; however, the sign change happens substantially before the formal boundary. At the formal Schwarzschild boundary, the temperature gradient dropped to about 𝒪(10⁻³) below the adiabatic value. The comparison of panels b and c shows that there exists an extended region in the model in which the convective flux is positive while the temperature gradient is already sub-adiabatic. This region is also known as a Deardorff layer (Deardorff 1966) and has been observed in simulations of stellar convection (Chan & Gigas 1992; Muthsam et al. 1995, 1999; Tremblay et al. 2015; Käpylä et al. 2017; Kupka et al. 2018) or other Reynolds stress models (Kupka 1999a; Xiong & Deng 2001; Kupka & Montgomery 2002; Montgomery & Kupka 2004; Zhang & Li 2012).

Such a layer cannot exist in convection models, which do not have enough degrees of freedom. In MLT and the Kuhfuß 1-equation model, the convective flux is directly coupled to the super-adiabatic gradient and the convective velocities (Eq. (4) for the 1-equation model). This of course inhibits any region in which the convective flux and the super-adiabatic gradient have a different sign like in the Deardorff layer and forces the convective flux and the TKE to have the same penetration depth. The 3-equation model lifts the strong coupling of convective flux and super-adiabatic gradient by directly solving for two more variables (Π and Φ) and therefore allows for the existence of such a layer. Deardorff (1966), in the case of atmospheric conditions, argued that it is mainly the non-local term in the equation for Φ that supports the positive heat flux for sub-adiabatic temperature gradients (we also refer to Sect. 13 of Canuto 1992, and references therein). The diffusion term ℱ_Φ in the equation for the entropy fluctuations Eq. (3) acts as a source of entropy fluctuations even though no local source (a super-adiabatic temperature gradient) is present. This allows for the outward directed transport of entropy fluctuations, which is a positive convective flux even in sub-adiabatic layers. This has little impact on the stellar structure and evolution, as the temperature gradient remains nearly adiabatic, and the whole convection zone is chemically well mixed. The existence of such a layer is therefore expected, and confirms the physical relevance of the 3-equation Kuhfuß model. The extent of this layer is difficult to determine from a priori arguments, and asteroseismic analyses might provide observational constraints in the future.

In Fig. 3 we show the temperature gradients in the overshooting zone of the same 5 M_⊙ main-sequence model as in Fig. 1. At the formal Schwarzschild boundary the model temperature gradient has already a slightly sub-adiabatic value as discussed previously. Beyond the Schwarzschild boundary, the model temperature gradient does not drop to the radiative gradient immediately. Instead, it gradually transitions from slightly sub-adiabatic to radiative values in a rather narrow mass range. As a consequence, the model temperature gradient takes slightly super-radiative values in this transition region. However, the temperature gradient reaches a radiative value well before the boundary of the mixed region, indicated with the black dotted line in Fig. 3. Considering the small extent of the super-radiative region and the small deviation from the radiative temperature gradient, the overshooting zone in the 3-equation non-local model is mostly radiative. We would like to point out that the shape of the temperature gradient is not subject to assumptions about the thermal stratification (e.g., adiabatic, radiative or any gradual transition between both) in the overshooting zone but instead is a result of the convection model. In the transition region, the convective flux is negative due to the buoyancy braking, which effectively means that energy transport by convection is directed inwards instead of outwards. This effect is counter-balanced by increasing the energy transport by radiation, through an increased model temperature gradient (e.g., Chan & Sofia 1996).

Fig. 3. Temperature gradients in the overshooting zone of the same 5 M⊙ main-sequence model as in Fig. 1. The blue and orange lines indicate the radiative (∇rad) and adiabatic (∇ad) temperature gradients, respectively. The dashed green line indicates the model temperature gradient, ∇, obtained from the 3-equation non-local convection model. The dashed black line indicates the Schwarzschild boundary, while the dotted black line indicates the boundary of the well-mixed overshooting region. The inset shows the three temperature gradients from the centre to the surface of the stellar model. The selected stellar model has a central hydrogen abundance of Xc = 0.6.

The temperature gradient of the 3-equation model is comparable to results of different TCM approaches (Xiong & Deng 2001; Li & Yang 2007) for the base of the solar convective envelope. Both Zhang & Li (2012; their Figs. 6 and 7) and Xiong & Deng (2001; their Fig. 8) find a temperature gradient that transitions gradually from the adiabatic to the radiative value. They also find a Deardorff layer with a degree of sub-adiabaticity at the formal Schwarzschild boundary comparable to our findings. From the convective flux as presented in Xiong (1986) one also would expect a similar temperature gradient in the overshooting zone. Furthermore, the shape of the model temperature gradient is also in qualitative agreement with the discussion in Viallet et al. (2015). They argue that under the physical conditions in convective cores, in regions of overshooting efficient chemical mixing and a gradually transitioning temperature gradient are expected. In the 3-equation non-local model, the extent of the nearly adiabatic overshooting zone is controlled by the shape of the negative convective flux in the overshooting zone. For smaller (more negative) values of the convective flux (i.e. more efficient buoyancy braking) the temperature gradient is expected to be closer to the adiabatic value, while for larger (less negative) values it will be closer to the radiative temperature gradient. In Eq. (1) the negative convective flux and the dissipation term act as sink terms in the overshooting zone. Hence, the behaviour of the dissipation term will impact also on the convective flux and in turn on the value of the temperature gradient in the overshooting zone. In computations with the 1-equation non-local version of the theory, the negative convective flux is the dominant sink term for the TKE and the actual dissipation term is negligible (Fig. 8 of Paper I). This leads to more negative values of the convective flux and thus to a mostly adiabatic temperature gradient in the overshooting zone. We discuss this in more detail in Sect. 5.1 (we also refer to Fig. 10).

As discussed above and in Paper I, we implemented the increase in the dissipation rate in the overshooting zone through a decrease in the dissipation length scale, Λ. In the original version, Kuhfuß (1987) modelled the dissipation of the TKE by a Kolmogorov-type term (ϵ = C_Dω^3/2/Λ with Λ(r) as in Eq. (10)). The dissipation length scale, Λ, describes the scale over which the kinetic energy is dissipated such that in the Kolmogorov model, at fixed TKE, a shorter length scale results in an increased dissipation rate. In Fig. 4 we show the length scale, Λ (panel a) and the dissipation rate (panel b) computed in the same stellar model as in Fig. 1. At a fractional mass of about 0.05 the profile of the dissipation length scale shows a slight kink in this model, where the transition starts. At about 0.28 in fractional mass, the dissipation length scale, Λ, drops to zero, which means that convective motions stop at that point. This coincides with the dying of the TKE beyond the Schwarzschild boundary, as seen in Fig. 1, panel a. The onset of the decrease in the dissipation length scale, Λ, also coincides with the sign change of the super-adiabatic gradient (Fig. 1, panel c). Towards the centre of the model, the dissipation length scale, Λ, drops to zero as well. This is a result of the Wuchterl (1995) correction, also used in our implementation (cf. Paper I). We have also investigated the prescription of Roxburgh & Kupka (2007) as an alternative to that one of Wuchterl (1995) but found little difference with respect to the overshooting region (we also refer to Paper I).

Fig. 4. Dissipation length scale Λ (a) and dissipation rate (b) in a 5 M⊙ main-sequence model. The stellar model is the same as in Fig. 1. The vertical dashed line indicates the Schwarzschild boundary of the model.

Parameter sensitivity. We explore some of the parameter dependences in Appendix B. We studied the impact of the new parameter c₄, appearing in Eq. (7) (and Eq. (21) in Paper I), on the structure of the convective core. A comparison of TKE profiles for different values of c₄ is shown in Fig. B.1. The comparison demonstrates that although the parameter has some impact on the overshooting extent, the main properties of the model are not changed substantially. The parameter α_ω, appearing in Eq. (1), controls the non-local flux of the TKE. Because this flux is mainly responsible for the extension of the kinetic energy beyond the Schwarzschild boundary, one expects that this parameter impacts on the overshooting distance. Figure B.2 shows a comparison of different hydrogen profiles of a 5 M_⊙ star computed with different values of α_ω. For a higher value of the parameter, the size of the convective core is larger throughout the evolution of the model. Smaller parameter values reduce the convective core size. Although Kuhfuß (1987) provided a default value for α_ω, this value is not known a priori from the theory. Hence, a calibration will be necessary, for example from observations or 3D hydrodynamic simulations. The parameters of the non-local terms in the Π and Φ equations α_Π and α_Φ, appearing in Eqs. (2) and (3), respectively, have a negligible impact on the overshooting extent (Fig. B.4). Instead, they have a larger impact on the temperature gradient as the variables Π and Φ are more closely related to the temperature structure. By increasing the parameter α_Π, the magnitude of the super-adiabatic gradient increases, while the extent of this region stays the same. Decreasing the value of the parameter α_Φ is increasing the size of the super-adiabatic region, therefore reducing the size of the Deardorff layer, while the magnitude stays the same. This is expected because the non-local term in the Φ equation is the one that drives convection in the Deardorff layer.

3.2. Comparison to MLT

Mixing-length theory is still the most commonly used theory to describe convection in stars. Therefore, we now compare the results of the 3-equation model with results obtained from standard MLT, that is, without and with an additional treatment of overshooting. The fractional hydrogen abundances early on and at the end of the main sequence computed with the 3-equation model and an MLT model without overshooting are shown in the upper panel of Fig. 2. At the same central hydrogen abundance, the fully mixed region of the 3-equation model always extends past that of the MLT model. This is comparable to models that include ad hoc overshooting beyond the formal Schwarzschild boundary. In contrast to the ad hoc overshooting models, the overshooting in the 3-equation model results from the solution of the model equations. The lower panel of Fig. 2 shows a comparison of the convective fluxes in the 3-equation model and an MLT model. Both models were selected to have the same central hydrogen abundance. Ad hoc overshooting was included in this MLT model, and tuned in such a way as to obtain a model with the same core size as obtained from the 3-equation model. In the bulk of the convection zone, both fluxes show very close agreement. Beyond the Schwarzschild boundary, the narrow region with negative convective flux in the non-local model can be identified. This region is shown enlarged in the inset.

Finally, in Fig. 5 we compare the evolutionary tracks of the 3-equation model, indicated by the solid blue line, with an MLT model without ad hoc overshooting (yellow line) and an MLT model including ad hoc overshooting (black line). The black dot indicates the position of the stellar model with X_c = 0.6 discussed in Sect. 3.1 in Figs. 1–2. The computation of the 3-equation model starts at the beginning of the main sequence from an MLT model including diffusive overshooting as described above and then evolves through core hydrogen burning up until core hydrogen exhaustion. Compared to the MLT model, the non-local model shows a higher luminosity throughout the main sequence, as expected for the larger convectively mixed core.

Fig. 5. Evolutionary tracks in the Hertzsprung-Russell diagram of a 5 M⊙ model computed with MLT (yellow line), ad hoc overshooting (black line), the 3-equation, non-local model (blue line) and the 1-equation, non-local model (red line). The black dot marks the model selected at a central hydrogen abundance of Xc = 0.6.

3.3. Comparison to the 1-equation non-local models

In addition to the 3-equation non-local models, we also computed stellar models in which convection was described by the 1-equation model (Sect. 2.1). As before, we included the non-local terms. Figure 6 shows a comparison of the TKE profiles for the 1- and 3-equation, non-local models on a logarithmic scale. The TKE on a linear scale can be found in Fig. A.1. The models were selected at the same central hydrogen abundance to ensure that they are in the same evolutionary stage. This shows that the overshooting extent in the 1-equation non-local model is very comparable to the 3-equation non-local model for the same choice of the parameter α_ω. The overall behaviour of the 1-equation non-local model and the new 3-equation non-local model looks very similar. The overshooting extent is clearly limited, and there is a steep drop in the TKE at the overshooting boundary. Also, the absolute values of the TKE look comparable in the bulk of the convection zone (Fig. A.1). We analyse the absolute value of the TKE in more detail below. We note that this result is obtained without tuning the parameters of the models. For the parameters that both models have in common, the same values were chosen. In the overshooting zone, the 3-equation model has much smaller TKE than the 1-equation model. Nevertheless, the energies are still high enough to fully mix the overshooting region in the 3-equation model.

Fig. 6. Comparison of the TKE in the 1-equation non-local model and the 3-equation non-local model with an enhanced dissipation rate of TKE in a 5 M⊙ main-sequence model on a logarithmic scale. The models have been selected to have the same central hydrogen abundance as the model in Fig. 1. The 1-equation model computes Λ according to Eq. (10).

In the upper panel of Fig. 2 we compare the hydrogen profiles of the 1-equation model with the results from the 3-equation model and an MLT model at the beginning and at the end of the main sequence. As expected from the similar TKE profiles shown in Fig. 6 the hydrogen profile of the 1-equation model extends past the local MLT model and looks very similar to the 3-equation model. Towards the end of the main sequence, the 1-equation model has a smaller core than the 3-equation model, leading to a slightly different slope in the hydrogen profiles. Likewise, the evolutionary track shown in Fig. 5 of the 1-equation model looks very similar to the 3-equation model. Towards the end of the main sequence, the luminosity is slightly lower owing to the smaller convective core. In the lower panel of Fig. 2 we compare the convective flux of the 1-equation model to the 3-equation model and an MLT model including ad hoc overshooting. As for the 3-equation model, the convective flux shows close agreement with the other two models in the bulk of the convection zone. Only in the overshooting zone, where the convective flux becomes negative, discrepancies become apparent. Compared to the 3-equation model, the zone of negative convective flux is more extended and the absolute value is larger, that is, the convective flux is more negative. This can be attributed to the parametrisation of the convective flux in the 1-equation model according to Eq. (4). As discussed in Sect. 4 of Paper I (their Fig. 8) the buoyancy term, proportional to the convective flux, acts as the main sink term in the overshooting zone, requiring larger absolute values in the overshooting zone. We would like to note that this strongly negative convective flux in the overshooting zone will also cause the 1-equation model to have a nearly adiabatic overshooting zone, as compared to the nearly radiative overshooting zone in the 3-equation model.

4. Non-local convection for varying initial masses

We computed stellar models in a mass range of 1.5–8 M_⊙, using the 3-equation non-local model. The models were constructed in the same way and using the same parameters as for the 5 M_⊙ model presented so far. For comparison, we computed three other sets of stellar models with different convection descriptions: (i) with the Kuhfuß 1-equation model to compare the results of the 3-equation model to a simpler TCM; (ii) with MLT plus diffusive overshooting as described by Freytag et al. (1996) to compare to one of the standard ad hoc descriptions of convective overshooting with the same parameter value f_OV = 0.02 as discussed above. Finally, (iii) we computed MLT models without overshooting to compare the results to a local convection theory. At least in terms of core size and temperature structure, models using the local Kuhfuß theories would be equivalent to MLT models. For the Kuhfuß theory, we used the same value of 0.3 for the parameter α_ω, as before. To allow for a comparison across the mass range and the different convection descriptions, the models were selected at the same central hydrogen abundance of X_c = 0.6.

Figure 7 shows the sizes of mixed cores in models derived from these four descriptions of convection and convective overshooting over a range of initial masses. For all masses under consideration, the mixed core from the 3-equation model is larger than the convective core from an MLT model, as expected. This shows that when applying the 3-equation model, an overshooting zone emerges across the whole mass range investigated. Comparing the 3-equation model to the 1-equation model and the ad hoc overshoot model, the mixed core sizes show good qualitative agreement. We repeat that this is achieved without fine-tuning any of the involved parameters. The relative size of the mixed cores decreases with decreasing stellar mass for all four descriptions, but differences in the details between the different convection descriptions are evident. For higher masses, the derived values for the mixed core sizes are almost identical among the ad hoc overshooting and the 1- and 3-equation Kuhfuß models. For low stellar masses, the results show larger discrepancies. Stellar models applying the 3-equation non-local model have the smallest cores, and the core size decreases faster with decreasing stellar mass than in the 1-equation models and the ad hoc overshoot models. As the ad hoc overshoot model was calibrated to observations, this allows at least for an indirect comparison of the 1- and 3-equation model with observations.

Fig. 7. Comparison of mixed core sizes of stellar models in units of stellar mass, M*, over a range of initial stellar masses computed with different convection models. The models with ad hoc overshoot include a geometric cutoff to limit the size of small convective cores in lower mass stars. The models were selected to have a central hydrogen abundance of Xc = 0.6. The MLT models are computed without overshooting.

In addition to the core sizes, we also analysed the absolute values of the convective velocities. The Kuhfuß model does not solve for the convective velocity itself, but rather for the TKE ω. We approximated the mean convective velocity from the TKE by assuming full isotropy:

$$\begin{array}{c}\hfill v_{\mathrm{c},\mathrm{iso}}=\sqrt{\frac{2}{3}\omega }.\end{array}$$(11)

When using MLT, the convective velocities were computed as

$$\begin{array}{c}\hfill v_{\mathrm{c}}=\frac{1}{\sqrt{8}}{\left(\frac{g}{H_{\mathrm{p}}}\right)}^{1/2}\mathrm{\Lambda }{(\mathrm{\nabla }-{\mathrm{\nabla }}_{\mathrm{ad}})}^{1/2},\end{array}$$

(e.g., Kippenhahn et al. 2012). The computation of convective velocities is the same in MLT models with and without diffusive overshoot, as the inclusion of ad hoc overshooting does not impact on the description of convection. To compare MLT and the Kuhfuß theories we evaluate the maximum convective velocity in the convection zone. The maximum is reached well within the Schwarzschild boundaries for all models, and therefore allows for a consistent comparison.

We show this comparison of the maximum convective velocities in the core as a function of stellar luminosity in Fig. 8. The convection descriptions are the same as in Fig. 7. For all cases, the scaling relation of the convective velocities has the same slope. A linear fit (dotted lines) to the data results in a value of ∼0.3 for all of them. However, the absolute values from the Kuhfuß and MLT models differ by a constant factor of about two, indicated by the offset between the two pairs of lines. This difference in the absolute value is the result of two different effects. The change of the mixing length has the largest impact on the velocity. A reduced mixing length leads to an increased dissipation rate and smaller velocities as a consequence. In the Kuhfuß models we used a smaller value of α = 1 as obtained by a solar calibration instead of α = 1.6 for the MLT models. In addition, the mixing length is reduced towards the centre according to the Wuchterl (1995) formulation (Eq. (6)). The convective velocity is reduced further compared to the local MLT models by taking the non-local terms into account, which act as a sink term in the bulk of the convection zone. The slope of 0.3 is very close to the v_c ∝ L^1/3 scaling relation expected from MLT. This comparison also demonstrates that the absolute values of the TKE are very similar between the 1- and the 3-equation Kuhfuß theories over the full mass range. For the 5 M_⊙ model, this was already apparent when comparing the TKE profiles in Fig. 6.

Fig. 8. Comparison of the maximum convective velocities for different convection models as a function of stellar luminosity. For the Kuhfuß model, the isotropic convective velocity, vc,iso, is plotted. The dotted lines indicate a linear fit to the logarithmic data. We note that the data of the 1- and 3-equation Kuhfuß models (red and blue points, respectively) are largely overlapping, as are the black and yellow dots.

5. Discussion

5.1. Relation between 1- and 3-equation models

As discussed in Sect. 2.1 the 1-equation model is a simplification of the 3-equation model, for which Kuhfuß (1987) assumed that the convective flux is proportional to the super-adiabatic temperature gradient and the square-root of the TKE. This allows the removal of two of the three equations, namely for the entropy fluctuations, Φ, and the velocity entropy correlations, Π. The equation for the TKE ω remains unchanged. The approximation for the convective flux allows the temperature gradient to be expressed as a function of the TKE. This couples the thermal structure and the TKE very closely. In the 3-equation model, the convective flux is evolved with an additional equation, which reduces the coupling of the thermal structure and the TKE. Despite the increased model complexity, the behaviour of the TKE in the 1- and 3-equation model is quite similar, as seen in Fig. 6.

In the bulk of the convection zone, both models result in the same absolute value of the TKE. This can be also seen in Fig. 8 by the agreement between red and blue points for a wider mass range. It can be attributed to the similarity of the equations. The additional dissipation is not or only weakly operating in the bulk of the convection zone. Hence, both the dissipation and the non-local term have the same functional form. The convective flux is adjusted such that a nearly adiabatic stratification is achieved in both models; therefore, the buoyancy term is also very comparable in the 1- and 3-equation models. As soon as the temperature gradient becomes sub-adiabatic, which happens already well within the Schwarzschild boundary in the 3-equation model, the dissipation by buoyancy waves is taken into account. This means that the functional form of the dissipation term changes. Due to the difference in the thermal structure in this region also the convective flux looks different. This leads to a different solution for the TKE in the overshooting zone, as is evident from Fig. 6 towards the edge of the convection zone.

The main difference between the 1- and 3-equation model is probably the temperature stratification that results from the solution of the model equations. Due to the coupling of the convective flux to super-adiabatic gradient and convective velocities in the 1-equation model, the convective flux is forced to have the same penetration depth into the stable layers as the TKE and at the same time to have the same sign as the super-adiabatic gradient. As seen in our models and pointed out by Xiong & Deng (2001) this strong coupling leads to a nearly adiabatic temperature gradient in the overshooting zone and prevents the existence of a Deardorff layer (we also refer to Paper I). Also, with respect to other models of convection (e.g., MLT), the existence of a large sub-adiabatic zone in the convective region for the case of the 3-equation model is striking. To better understand the behaviour of the temperature gradient in the overshooting zone, we analyse it in terms of the Peclet number. The Peclet number is the ratio of the timescales of radiative and advective transport and can be interpreted as an indicator for convective efficiency. A common definition of the Peclet number is

$$\begin{array}{c}\hfill \text{Pe}=\frac{u·l}{\chi },\end{array}$$

where u and l are a typical velocity and length scale of the convective flow. The radiative diffusivity is defined as

$$\begin{array}{c}\hfill \chi =\frac{16\sigma T^3}{3\kappa {\rho }^2{c}_{\mathrm{p}}}.\end{array}$$

As a typical convective velocity we again used the isotropic velocity, Eq. (11). Due to the usage of typical scales for velocities and length scales that are not rigorously defined, the interpretation of absolute values of the Peclet number remains difficult. Hence, we only looked at ratios of the Peclet number to compare different models. We also assumed that the length scale and the radiative diffusivity are the same, when comparing different models. Under these assumptions, it is easy to see that

$$\begin{array}{c}\hfill \frac{{\text{Pe}}_1}{{\text{Pe}}_3}\propto \sqrt{\frac{{\omega }_1}{{\omega }_3}.}\end{array}$$

The ratio of the Peclet numbers obtained for the 1- and 3-equation models for a 5 M_⊙ main-sequence model is shown in Fig. 9. In the bulk of the convection zone, the 1- and 3-equation models have very similar Peclet numbers, which means the transport of energy by convection behaves very comparably. In the overshooting zone, however, the 1-equation model has a Peclet number that is up to seven times higher than that of the 3-equation model. This indicates that convection as described by the 1-equation model is much more efficient in the overshooting zone than when described by the 3-equation model.

Fig. 9. Ratio of the Peclet numbers for the 1- and 3-equation models for a 5 M⊙ main-sequence model.

This change in efficiency also impacts the temperature gradient. Following the approximation of the convective flux in the 1-equation model, the convective flux scales as $\mathrm{\Pi }\propto \sqrt{\omega }$. Using the ratio of the Peclet numbers, one can therefore write a Peclet-scaled convective flux,

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{3,\text{Pe}}={\mathrm{\Pi }}_1·\frac{{\text{Pe}}_3}{{\text{Pe}}_1},\end{array}$$(12)

to mimic the convective flux in the 3-equation model. The resulting convective flux is shown with a blue dashed line in the upper panel of Fig. 10. Using this scaled convective flux, a scaled temperature gradient can be computed, as illustrated by the green dashed line in the lower panel of Fig. 10. For comparison, the model temperature gradients of the 1- and 3-equation model are visualised by a green dotted and solid line, respectively.

Fig. 10. Comparison of convective fluxes and temperature gradients obtained from the Peclet-scaling and the 1- and 3-equation models. Upper panel: convective flux as a function of fractional mass. The red and blue lines show the results obtained from the 1- and 3-equation model, respectively. The dashed blue line indicates the convective flux scaled with the Peclet number according to Eq. (12). Lower panel: temperature gradients as a function of fractional mass. The blue and orange lines indicate the radiative and adiabatic temperature gradient, respectively. The dotted and solid green lines indicate the temperature gradient as obtained by the 1- and 3-equation models. The dashed green line shows the temperature gradient of the 1-equation model computed from the scaled convective flux shown in the upper panel.

This comparison confirms that the reduced convective efficiency obtained from the Peclet numbers is sufficient to change the behaviour of the temperature gradient from nearly adiabatic to more radiative in the overshooting zone, and implies that the behaviour of the temperature gradient in the overshooting zone can at least qualitatively be predicted from the TKE alone without invoking the other convective equations for Π and Φ. The fact that the region of negative values in the scaled convective flux is more extended compared to the actual convective flux from the 3-equation model is due to the different penetration depths of TKE and convective flux in the 3-equation model. We conclude that this is an important indication for internal consistency of the model. It furthermore shows that the mostly radiative temperature gradient in the overshooting zone is in fact a result of reduced convective efficiency in the overshooting zone in the 3-equation model. Following the terminology proposed by Zahn (1991), in the 1-equation model the overshooting zone is best described by sub-adiabatic penetration while the more inefficient convection in the 3-equation model concerns overshooting of chemical element distributions only.

Considering the chemical mixing, both models result in a more step-like chemical mixing profile. The extent of the mixed region is mainly dependent on the choice of the parameter α_ω, as discussed in Appendix B. This parameter cannot be determined from first principles. As for the ad hoc descriptions of convective core overshooting, a calibration of this parameter is required, which is discussed below. Furthermore, the resulting convective flux is also very similar in the bulk of the convection zone and differences become only obvious in the overshooting zone. Given the relative freedom in choosing α_ω and the similarity in the mixing properties of both models, resulting stellar models are basically indistinguishable when comparing the chemical structure. Once the observations become sensitive enough to the thermal structure in the overshooting zone as, for example, with the help of asteroseismology (Michielsen et al. 2019), it will become possible to detect differences between both models. In view of the general agreement between 1- and 3-equation model, the application of the 1-equation model seems to be sufficient to obtain the chemical structure from a non-local convection model. The parameter α_ω can be tuned to obtain the correct size of the convective core. However, the stratification obtained from the 1-equation model is less realistic.

5.2. Other constraints on core overshooting

To date, a range of different approaches to determining the extent of convective cores has been followed. In this subsection, we discuss a few examples of convective core size determinations. The need for a larger mixed core was already recognised in the 1980s. For example, Bressan et al. (1981) discussed it in relation to the Hertzsprung-Russell-diagram of massive stars, and Maeder & Mermilliod (1981) in relation to the Hertzsprung-Russell-diagram of open clusters. In the latter case, isochrones derived from stellar models that include core overshooting match the morphology of the turn-off region better than models without overshooting (Pietrinferni et al. 2004; Magic et al. 2010). The comparison to observations showed further that the overshooting distance in terms of pressure scale heights needs to increase with mass in the range between 1.2 and ∼2 M_⊙. This mass dependence can be included explicitly in the computations by expressing the overshooting parameter as a function of the total stellar mass. Alternatively, this mass dependence can also be introduced in the stellar models by limiting the radial extent of the overshooting zone geometrically. In both cases, the parameter of the overshooting scheme effectively needs to increase with stellar mass. For the TCM, however, this is a natural outcome without imposing it.

Eclipsing binary systems offer an excellent opportunity to put constraints on stellar physics. Claret & Torres (2019) used a large sample of eclipsing binaries to determine the overshooting parameter as a function of mass. They find a clear increase in the overshooting parameter (extent) with mass, even though the statistical significance of this result has been debated (Constantino & Baraffe 2018). In a detailed analysis, Higl et al. (2018) addressed the evolution of the binary system TZ Fornacis with an evolved red giant primary and a main-sequence secondary star, both with masses of ∼2 M_⊙. They find that a basically unrestricted overshooting extent, using the standard value for the free parameter, is required to explain the evolution of this system. This puts further constraints on the mass dependence of the overshooting parameter.

The high precision photometric data obtained from space telescopes such as Kepler or CoRoT allowed further constraints to be set on the stellar evolution models. Using asteroseismology of g-mode pulsators, the convective core masses of intermediate mass stars has been determined for larger samples of stars (e.g., Pedersen et al. 2021; Mombarg et al. 2019). Similarly, the seismology of p-mode pulsators allows the required overshooting efficiency in lower mass stars to be determined (e.g., Deheuvels et al. 2016; Angelou et al. 2020). In agreement with the previously mentioned studies, they find that the relative mass of the mixed core needs to increase with stellar mass. A more detailed comparison of the TCM models discussed in this work with asteroseismic observations needs to be addressed in future work. Finally, asteroseismology allows the temperature gradient in the overshooting zone to be probed. Michielsen et al. (2021) infer a predominantly radiative overshooting zone in a ∼3.5 M_⊙ main-sequence star, which is in agreement with the temperature gradient obtained from the 3-equation model but disagrees with the 1-equation model. However, they point out that this result is only obtained for this single B-type star and might not be generalisable for all B-type stars.

With the increase in computational resources in recent years, more and more multi-dimensional hydro-simulations of stellar core convection have been carried out (e.g., Meakin & Arnett 2007; Gilet et al. 2013; Edelmann et al. 2019; Higl et al. 2021). These simulations confirmed, for example, the scaling of the stellar luminosity with the third power of the convective velocities (e.g., Edelmann et al. 2019; Higl et al. 2021, and references therein). In Higl et al. (2021) the authors calibrated the overshooting parameter f_OV in GARSTEC to 2D simulations of core convection in low and intermediate mass stars. By matching the size of the mixed convective core in the 1D GARSTEC models to the size of the mixed region in the 2D simulations, they find that the effective parameter of f_OV needs to decrease with stellar mass. To limit the size of the convective cores for small stellar mass in the 1D models, they used the geometric ‘square’ cutoff according to Eq. (C.2). As they used the same stellar evolution code as we do, this allows for a direct comparison of the results. Higl et al. (2021) find that the size of the convective cores resulting from the 2D simulations need to be larger than the GARSTEC models computed including the square cutoff (Eq. (C.2)) at a constant overshooting parameter. This indicates that the geometric square cutoff is too restrictive. As the 3-equation model predicts mixed core masses similar to the GARSTEC models that include this geometric cutoff (Fig. C.1), this indicates that our 3-equation model might be too restrictive at this lower mass range as well.

Finally, other TCMs have been used to compute stellar evolution models including the effects of non-local convection. Xiong (1986) computed stellar models in the mass range of 7–60 M_⊙. He found that by solving the convection equations, the TKE extends beyond the formal Schwarzschild boundary and that this also increases the size of the mixed region. The size of the well-mixed region increases with increasing stellar mass. Both results are in agreement with our findings employing the 3-equation model. Xiong (1986) also found that the convective flux penetrates much less deeply into the stable layers than the TKE. The magnitude of the convective flux in this region is negligible, which causes the temperature gradient to be mostly radiative in the overshooting region. This is in good agreement with the decoupling of the thermal and the chemical structure we discussed in detail in Sect. 3 and Paper I. Zhang (2016) applied the TCM by Li & Yang (2007) in a similar mass range as in this work. They develop a simplified model comparable to the 1-equation model and find very good agreement between the full model and the simplified version. Li (2017) applied the simplified TCM by Li (2012) to compute stellar models of a 5 M_⊙ star and find an overshooting distance of about 0.2H_p, which is comparable to the overshooting distance obtained with our 3-equation model.

6. Conclusions

In this work we present results of stellar structure and evolution calculations using the TCM proposed by Kuhfuß (1987). We implemented the Kuhfuß model into GARSTEC, which solves the four stellar structure equations and the three equations of the convection model simultaneously using the implicit Henyey method. We computed main-sequence models of intermediate mass main-sequence stars between 1.5 and 8 M_⊙ that consistently compute the structure and evolution of the TKE, convective flux, and entropy fluctuations. This naturally includes the effects of convective overshooting for the thermal and chemical structure. In Paper I we demonstrate that the original 3-equation model with a standard MLT prescription for the dissipation length of TKE leads to convection zones that essentially extend throughout the entire stellar interior. We therefore implemented the dissipation by gravity waves as discussed in Paper I in addition to the original Kuhfuß model. We show that the Kuhfuß 3-equation model with an increased dissipation rate results in models with physically reasonable overshooting distances. This indicates that the dissipation was actually underestimated by the original description, and that dissipation by gravity waves is a relevant effect in core overshooting zones.

In Fig. 1 we show a summary of the TKE and flux and the temperature gradient for a convective core of a 5 M_⊙ main-sequence model. We find that the TKE extends beyond the formal Schwarzschild boundary of convective neutrality. This is the result of the non-local terms in the Kuhfuß model. The convective flux shows a region of negative values beyond the Schwarzschild boundary, which is, however, penetrating less deeply into the stable layers than the TKE. As the convective motions are very efficient in mixing chemical elements, the extended convective core has essentially the same composition as the convective core. This can be seen in Fig. 2 (upper panel), in which the hydrogen profiles at the end of the main sequence of an MLT and a Kuhfuß 3-equation model are compared. In the Kuhfuß model, this extension beyond the Schwarzschild boundary is the outcome of the solution of the model equations and not due to the inclusion of any sort of ad hoc overshooting. We also compared the results of the full 3-equation model to the simplified 1-equation model and find qualitative and quantitative agreement of the TKE throughout a large part of the convection zone. This is a result of the similarity of the model equations and the chosen parameters, which are the same for both models. In addition to the convective velocity and the associated mixing, the temperature gradient is also part of the model solution. This is another important difference compared to ad hoc descriptions of convective overshooting, in which the temperature gradient needs to be assumed separately and independently.

The analysis of the temperature gradient shows the existence of a Deardorff layer (Deardorff 1966), in which the temperature gradient is sub-adiabatic and the convective flux is still positive. The existence of the Deardorff layer has been confirmed in different numerical simulations of stellar convection (Chan & Gigas 1992; Muthsam et al. 1995, 1999; Tremblay et al. 2015; Käpylä et al. 2017) and other Reynolds stress models (Kupka 1999a; Xiong & Deng 2001; Kupka & Montgomery 2002; Montgomery & Kupka 2004; Zhang & Li 2012). Beyond the Deardorff layer, the convective flux becomes negative as the result of the stable stratification (cf. also Muthsam et al. 1995 and references in Canuto 1992). In the overshooting region, the model temperature gradient gradually transitions from a slightly sub-adiabatic to a radiative value, exhibiting a small region with a super-radiative temperature gradient (Fig. 3). However, this transition region is rather narrow, such that the overshooting zone has a mostly radiative temperature gradient. This is in agreement with very recent results from asteroseismology (Michielsen et al. 2021). In contrast to the 3-equation model, the overshooting zone of the 1-equation model shows a mostly adiabatic temperature gradient. As pointed out by Xiong & Deng (2001), this can be attributed to the assumption of a full correlation between the convective flux and the convective velocities, as is done in the 1-equation model (Eq. (4)). The approximation Eq. (4) does not allow for a Deardorff layer, as the convective flux is proportional to the super-adiabatic gradient. From a theoretical point of view, the existence of a Deardorff layer can be attributed to the non-local term of the Φ equation, Eq. (3), Deardorff (1966). This shows that an independent equation for the convective flux is required (see also Kupka 2020, Priv. comm.) and highlights the necessity to consider more complex turbulence models, such as the 3-equation model, to capture the temperature structure in the overshooting zone more accurately. The comparison of the Peclet numbers of the 1- and 3-equation models further shows that the mostly radiative temperature gradient in the overshooting zone of the 3-equation model can be explained by the reduced TKE or velocities, respectively, compared to the 1-equation model. The narrow range of the transition region from an adiabatic to a radiative temperature gradient can be attributed to the shallower penetration of the convective flux into the stable layers compared to the 1-equation model.

A comparison with stellar models that use MLT shows a qualitative agreement of the TKE and the convective flux in the part of the convection zone that is unstable according to the Schwarzschild criterion (∇_ad = ∇_rad). For the convective flux, we even find a very good quantitative agreement between the Kuhfuß model and MLT in that region (Fig. 2, lower panel). The convective velocities found in the Kuhfuß model are smaller than in MLT by a factor of two. This qualitative agreement indicates that the stellar structure has the largest impact on the convective properties in the bulk of the convection zone, irrespective of the convection model in use. The nuclear energy released in the centre determines the convective flux – about 80% of the local flux in the centre – and the coefficients of the convection model determine the absolute values of the other convective variables. Differences appear in the overshooting zone, which is sensitive to subtler changes in the convection model. In the overshooting zone the stellar flux is mainly transported by radiation, and as such the convective structure of this region is less constrained by the stellar structure. It is, however, rather constrained by the convection model.

The results of the 1- and 3-equation models over a broader mass range show qualitative agreement with other overshoot descriptions. For given values of α_ω and f_OV, this is achieved without fine-tuning any of the other model parameters, as we used the closure parameters suggested by the authors of the turbulence model (see Canuto & Dubovikov 1998, and references therein). Tests of those parameters – for different physical scenarios – are published in the literature (Kuhfuß 1986, 1987; Wuchterl 1995; Canuto 1992; Canuto & Dubovikov 1998; Kupka et al. 2007c). In Fig. 7 we show a comparison of mixed hydrogen core sizes computed with diffusive overshooting, the 1-equation and the 3-equation model, and MLT. The comparison shows that the exact extent of the convective core depends on the details of the model. For the same parameter choice of the parameter α_ω at lower masses, the 3-equation model shows a reduced amount of overshooting compared to the 1-equation model, while the 3-equation models have larger cores at higher masses. The newly introduced parameters that control the reduction of the dissipation length scale, Λ, are shown to have a moderate impact on the overshooting extent. For its default overshoot parameter, the diffusive mixing model produces even larger mixed core sizes. Towards the lower end of the mass range, both Kuhfuß models show a decrease in the mixed core size as implemented by other methods to match observations from open clusters and binaries (Claret & Torres 2019; Pietrinferni et al. 2004; Magic et al. 2010). Higl et al. (2018) find that the overshooting parameter needs to increase more steeply with mass than predicted by GARSTEC, including the geometric square cutoff from the analysis of the TZ For binary system. In agreement with the conclusion from the 2D simulations, (Higl et al. 2021), this indicates that the convective core size predicted by the 1- and 3-equation model increases too shallowly with stellar mass for a constant parameter α_ω. To match the observations and the results from the 2D simulations, this would require modifying the parameter α_ω as a function of mass. Further constraints on the parameter α_ω can be obtained by comparing the Kuhfuß model to results from asteroseismology of intermediate mass stars.

However, such tuning of α_ω seems ill advised: the physical incompleteness of the ad hoc model of convective overshooting, as an example, is well demonstrated by the shrinking of the convective core size with mass for stars with M < 2 M_⊙; this not only requires an extra cutoff function (Eq. (C.2)) to limit the convective core size to values compatible with observations, but also demands a more fine-tuned function (Eq. (C.3)) to pass such a stringent test. While applying a similar procedure to α_ω as for f_OV appears to be a convenient ad hoc solution to match the observational data exactly, it provides no new insights and requires redoing similar procedures in related, but different, scenarios. Finding the physical reason for the remaining, now already much smaller, discrepancies with the data, on the other hand, might allow for a model that does not require such measures in other applications either (cf. also the discussion on such requirements in Kupka & Muthsam 2017).

Our study of stellar models that apply the 3-equation model from Kuhfuß (1987) shows that the resulting stellar structure depends sensitively on the details of the convection model. Although the original 1- and 3-equation theories are very closely related, their application results in very different structures of the convection zone (Paper I). We have demonstrated that modifying the dissipation term of the TKE can remedy this discrepancy. The improved 3-equation model compares very well with the 1-equation model, both in terms of the TKE and the mixing properties. For applications in which the temperature structure of the overshooting zone is not important, the 1-equation model describes the convective core similarly well as the 3-equation model. The parameter α_ω allows us to obtain the convective core size required to match with the observations. Only when the thermal structure of the overshooting zone is of major interest does a more complex convection model such as the 3-equation model need to be used. A self-consistent prediction of the detailed convective core structure appears to require a physically more complete (and thus more complex) model. As discussed above, at least three equations are required to allow for more complex phenomena, such as a Deardorff layer. In Canuto & Dubovikov (1998), further partial differential equations for the dissipation rate and the vertical component of the TKE are discussed. This potentially allows some of the simplifications still present in the 3-equation model to be removed. However, increasing the number of equations comes at the price of increasing numerical complexity. Including these equations in stellar structure and evolution models has to remain the task of future work. Finally, a comparison with more and more realistic 3D hydro-simulations of convective cores will be useful for placing further constraints on our TCM model, in particular to restrict closure conditions. This is presently underway (Ahlborn & Higl, in prep.). The application of the extended 3-equation Kuhfuß model to stellar envelopes will be presented in another forthcoming paper (Ahlborn et al., in prep.). In any case, TCMs offer a convincing and feasible improvement of the treatment of convection in 1D stellar models beyond the standard MLT.

Acknowledgments

F.K. is grateful to the Austrian Science Fund FWF for support through projects P29172-N and P33140-N and support from European Research Council (ERC) Synergy Grant WHOLESUN #810218

References

Appendix A: Comparison of 1- and 3-equation models

Figure A.1 shows the TKE in the 1- and 3-equation model on a linear scale. Even though the TKE of the 3-equation model has approximately the same extent as the 1-equation model as seen in Fig. 6 it drops off faster. Despite the smaller values of the TKE, the efficiency of the chemical mixing is still high enough to fully mix the overshooting region. Therefore, the chemical profiles obtained from the 1- and 3-equation model are very comparable.

Fig. A.1. Comparison of the TKE of the 1-equation and 3-equation non-local models in a 5 M⊙ main-sequence model with limited dissipation length scale, Λ, for the 3-equation model.

Appendix B: Parameter dependence

The model for the modification of the dissipation length scale comes with a number of parameters (c₂, c₃, c_ϵ). These parameters are not necessarily free, but were calibrated in the framework of the convection models developed by Canuto (1992). These parameters enter the equation of the reduction factor as a single parameter, which was previously defined as c₄. Given the dissipation parameter of the 3-equation model, we find c₄ ≈ 0.072 while for the dissipation parameter of the Canuto & Dubovikov (1998) model c₄ ≈ 0.2 as described in Paper I. As mentioned already in Sect. 3.6 in Paper I, the effect of changing c_ϵ on changing c₄ to some extent cancels out in the calculation of Λ, as c_ϵ appears in the denominator of c₄ = c₃/(c₂c_ϵ) and in the numerator of Λ = c_ϵω^3/2/ϵ. Hence, if c₄ is adjusted according to c₄ = c₃/(c₂c_ϵ), a change of c_ϵ first of all influences the TKE dissipation rate ϵ throughout the whole model and is not specifically changing ϵ only within the overshooting zone. To check what impact these parameters actually have on the result, we varied their values. As they enter the equations as one parameter, we only varied this effective parameter value c₄ by ±60%. All other parameters take their default values. The resulting profiles of the TKE are shown in Fig B.1. It can be seen that the variation of the parameter c₄ leads to some noticeable variation in the TKE profile. However, within these ranges, the models keep their property of a limited overshooting range. The direction of the variation can be explained by the theory as well. An increase in the parameter will lead to a decrease in the dissipation length scale, Λ. A decreased length scale leads in turn to an increased dissipation. This reduces the overshooting, which can be observed for the dark red line. For the case with a decreased parameter value, the same argument applies in the opposite direction. The expected behaviour can be similarly observed for the yellow line.

Fig. B.1. Comparison of the TKE as a function of fractional mass on a logarithmic scale for different values for the parameter c4 in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.

Apart from the new parameters for the reduction of the dissipation length scale, Λ, the model still contains the original parameters of the Kuhfuß model. Here, we focus on the parameter α_ω, which controls the flux of the TKE. In the 1-equation non-local model, this parameter controls the extent of the overshooting region. Kuhfuß suggested a default value of α_ω = 0.25. In Fig. B.2 we show three different hydrogen profiles for the 5 M_⊙ model at the same age. The parameter α_ω takes values of 0.1, 0.3 and 0.5. It can be clearly seen that the original property of this parameter of controlling the overshooting extent is still given.

Fig. B.2. Comparison of the hydrogen profiles as a function of fractional mass for different values for the parameter αω in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.

In Fig. B.3 we show hydrogen profiles of a 5 M_⊙ star at the end of the main sequence for different values of the dissipation parameter C_D = 0.79, 1, 2.18 and 3. We always assume C_D = c_ϵ, as both parameters have the same role in the Canuto & Dubovikov (1998) and in the Kuhfuß model. Here, the value of 0.79 refers to the default dissipation parameter from the Canuto & Dubovikov (1998) model, while 2.18 is the numerical default value in the Kuhfuß model. The extent of the hydrogen profile is largest for the model computed with the smallest dissipation parameter and smallest for the largest parameter. This behaviour is expected, as a decreased dissipation allows the TKE flux to extend farther out. Compared to the parameter of the TKE flux α_ω, the impact of the dissipation parameter on the overshooting extent is rather limited, as the variation is much smaller when compared to the results shown in Fig. B.2. The variation of the overshooting extent is also smaller, as one could have expected from the comparison shown in Fig. B.1. This is because by changing C_D also c₄ will change, while in Fig. B.1 only c₄ is changed. The effects of changing c₄ and C_D partially compensate for each other, resulting in a smaller net effect. Finally, one could find combinations of parameters α_ω and C_D that allow models with equal convective core sizes to be obtained. In Fig. B.4 we show the impact of the parameter α_Π and α_Φ in the upper and lower panel, respectively. They both have a negligible impact on the overshooting distance, which is yet smaller than the impact of the dissipation parameter C_D. This indicates again that the parameter α_ω, which determines the importance of the TKE flux, has the largest impact on the overshooting distance (we note that Fig. B.1 has a different scale for the fractional mass axis).

Fig. B.3. Comparison of the hydrogen profiles as a function of fractional mass for different values for the parameter CD in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.

Fig. B.4. Comparison of the hydrogen profiles as a function of fractional mass for different values for the parameters αΠ and αΦ in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.

Appendix C: Ad hoc overshooting for small convective cores

Parametrised descriptions of overshooting like the diffusive or step overshooting use the pressure scale-height at the Schwarzschild boundary to define the overshooting distance. In the description of Freytag et al. (1996) the diffusion coefficient is computed according to

$$\begin{array}{c}\hfill D_{\mathrm{adhoc}}=D_{0}exp\left(\frac{-2z}{f_{\mathrm{OV}}{H}_{\mathrm{p}}}\right),\end{array}$$(C.1)

where f_OV is an adjustable parameter and z is the radial distance to the Schwarzschild boundary. The pressure scale-height is, however, diverging towards the stellar centre. Hence, the inferred overshooting distance for a fixed overshooting parameter will increase for a decreasing convective core size as well as the size of the mixed region. The resulting mixed core sizes for a fixed parameter value have been shown to be too large when comparing to observations (Pietrinferni et al. 2004; Magic et al. 2010). To avoid the unphysical growth of the overshooting region, the overshooting parameter needs to be artificially restricted. In GARSTEC the unphysical growth of the overshooting zone is prevented by comparing the size of the Schwarzschild core to the pressure scale height and use the smaller one as the relevant length scale. As usual, there are different ways to implement this. Originally, Magic et al. (2010) suggested the following expression:

$$\begin{array}{c}\hfill \stackrel{\sim }{H_{\mathrm{p}}}=H_{\mathrm{p}}·\min (1,{\left(\frac{r_{\mathrm{CZ}}}{H_{\mathrm{p}}}\right)}^2),\end{array}$$

where r_CZ is the radius of the Schwarzschild boundary. A correction factor of the same type was also used in Higl et al. (2021) when comparing GARSTEC results to 2D simulations, while introducing a factor of two in the denominator:

$$\begin{array}{c}\hfill \stackrel{\sim }{H_{\mathrm{p}}}=H_{\mathrm{p}}·\min (1,{\left(\frac{r_{\mathrm{CZ}}}{2H_{\mathrm{p}}}\right)}^2).\end{array}$$(C.2)

As the size of the convective core is now compared to a length scale twice as large as the original expression, the size of the overshooting region is limited more strongly by Eq. (C.2). The comparison to 2D simulations (Higl et al. 2021) as well as the study of the eclipsing binary TZ For (Higl et al. 2018) reveal, however, that this expression finally limits the size of the convective core too strongly. This led to the introduction of a different functional form of the limitation:

$$\begin{array}{c}\hfill \stackrel{\sim }{H_{\mathrm{p}}}=H_{\mathrm{p}}·\min (1,\frac{1}{2}[\tanh \left(5(\frac{r_{\mathrm{CZ}}}{H_{\mathrm{p}}}-1)\right)+1]),\end{array}$$(C.3)

which less strongly limits the core sizes at 2 M_⊙ but very quickly limits the size of the mixed cores for smaller masses.

In Fig. C.1 we show a comparison of the mixed core sizes obtained without any cut, with the square cut according to Eq. (C.2) and with the tanh cut according to Eq. (C.3). For reference, results obtained with MLT and the 1- and 3-equation models are shown in the same figure. As discussed above, the square cut is more restrictive than the tanh cut at masses around and above 2 M_⊙. Only at a mass of 8 M_⊙ does the square cut no longer restrict the convective core size. In contrast, the tanh cut restricts the core size only marginally, already at 2 M_⊙. For lower masses below about 4 M_⊙, the results of the square cut are in good agreement with the 3-equation model.

Fig. C.1. Comparison of mixed core sizes obtained with different descriptions of convection and different geometric cuts to limit the overshooting (OV) for small cores in the ad hoc overshooting model.

All Figures

Fig. 1. Summary of the interior structure of the convective core of a 5 M⊙ main-sequence model, calculated with the 3-equation model. The dashed black line indicates the Schwarzschild boundary. The different panels show: (a) the TKE, (b) the convective flux variable, (c) the super-adiabatic temperature gradient, and (d) the diffusion coefficient according to Eq. (5). The selected stellar model has a central hydrogen abundance of Xc = 0.6. The inset in panel b shows the region of negative convective flux just beyond the Schwarzschild boundary.

In the text

Fig. 2. Comparison of models (yellow lines) using MLT without overshooting to those computed with the 3-equation theory (blue lines) and the 1-equation model (red lines). Panel a: compares the hydrogen profiles at an early stage on the main sequence, when Xc = 0.6 (solid lines), and at the end of it (Xc ≈ 0; dashed lines). Panel b: convective fluxes of an MLT model with diffusive overshooting (dotted black line), a 3-equation model (solid blue line), and a 1-equation model (solid red line). These models have been selected to have the same central hydrogen abundance of Xc = 0.6 and the same chemically homogeneous core size.

In the text

Fig. 3. Temperature gradients in the overshooting zone of the same 5 M⊙ main-sequence model as in Fig. 1. The blue and orange lines indicate the radiative (∇rad) and adiabatic (∇ad) temperature gradients, respectively. The dashed green line indicates the model temperature gradient, ∇, obtained from the 3-equation non-local convection model. The dashed black line indicates the Schwarzschild boundary, while the dotted black line indicates the boundary of the well-mixed overshooting region. The inset shows the three temperature gradients from the centre to the surface of the stellar model. The selected stellar model has a central hydrogen abundance of Xc = 0.6.

In the text

	Fig. 4. Dissipation length scale Λ (a) and dissipation rate (b) in a 5 M⊙ main-sequence model. The stellar model is the same as in Fig. 1. The vertical dashed line indicates the Schwarzschild boundary of the model.
In the text

	Fig. 5. Evolutionary tracks in the Hertzsprung-Russell diagram of a 5 M⊙ model computed with MLT (yellow line), ad hoc overshooting (black line), the 3-equation, non-local model (blue line) and the 1-equation, non-local model (red line). The black dot marks the model selected at a central hydrogen abundance of Xc = 0.6.
In the text

	Fig. 6. Comparison of the TKE in the 1-equation non-local model and the 3-equation non-local model with an enhanced dissipation rate of TKE in a 5 M⊙ main-sequence model on a logarithmic scale. The models have been selected to have the same central hydrogen abundance as the model in Fig. 1. The 1-equation model computes Λ according to Eq. (10).
In the text

Fig. 7. Comparison of mixed core sizes of stellar models in units of stellar mass, M*, over a range of initial stellar masses computed with different convection models. The models with ad hoc overshoot include a geometric cutoff to limit the size of small convective cores in lower mass stars. The models were selected to have a central hydrogen abundance of Xc = 0.6. The MLT models are computed without overshooting.

In the text

Fig. 8. Comparison of the maximum convective velocities for different convection models as a function of stellar luminosity. For the Kuhfuß model, the isotropic convective velocity, vc,iso, is plotted. The dotted lines indicate a linear fit to the logarithmic data. We note that the data of the 1- and 3-equation Kuhfuß models (red and blue points, respectively) are largely overlapping, as are the black and yellow dots.

In the text

	Fig. 9. Ratio of the Peclet numbers for the 1- and 3-equation models for a 5 M⊙ main-sequence model.
In the text

Fig. 10. Comparison of convective fluxes and temperature gradients obtained from the Peclet-scaling and the 1- and 3-equation models. Upper panel: convective flux as a function of fractional mass. The red and blue lines show the results obtained from the 1- and 3-equation model, respectively. The dashed blue line indicates the convective flux scaled with the Peclet number according to Eq. (12). Lower panel: temperature gradients as a function of fractional mass. The blue and orange lines indicate the radiative and adiabatic temperature gradient, respectively. The dotted and solid green lines indicate the temperature gradient as obtained by the 1- and 3-equation models. The dashed green line shows the temperature gradient of the 1-equation model computed from the scaled convective flux shown in the upper panel.

In the text

	Fig. A.1. Comparison of the TKE of the 1-equation and 3-equation non-local models in a 5 M⊙ main-sequence model with limited dissipation length scale, Λ, for the 3-equation model.
In the text

	Fig. B.1. Comparison of the TKE as a function of fractional mass on a logarithmic scale for different values for the parameter c4 in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.
In the text

	Fig. B.2. Comparison of the hydrogen profiles as a function of fractional mass for different values for the parameter αω in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.
In the text

	Fig. B.3. Comparison of the hydrogen profiles as a function of fractional mass for different values for the parameter CD in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.
In the text

	Fig. B.4. Comparison of the hydrogen profiles as a function of fractional mass for different values for the parameters αΠ and αΦ in a 3-equation non-local 5 M⊙ main-sequence model with limited dissipation length scale, Λ.
In the text

	Fig. C.1. Comparison of mixed core sizes obtained with different descriptions of convection and different geometric cuts to limit the overshooting (OV) for small cores in the ad hoc overshooting model.
In the text