Issue 
A&A
Volume 667, November 2022



Article Number  A97  
Number of page(s)  17  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202243126  
Published online  11 November 2022 
Stellar evolution models with overshooting based on 3equation nonlocal theories
II. Mainsequence models of A and Btype stars
^{1}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStraße 1, 85748 Garching, Germany
email: fahlborn@mpagarching.mpg.de
^{2}
Dept. Applied Mathematics and Physics, Univ. of Applied Sciences, Technikum Wien, Höchstädtplatz 6, 1200 Wien, Austria
^{3}
WolfgangPauliInstitute c/o Faculty of Mathematics, University of Vienna, OskarMorgensternPlatz 1, 1090 Wien, Austria
^{4}
MaxPlanckInstitut für Sonnensystemforschung, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
Received:
15
January
2022
Accepted:
18
July
2022
Context. Convective overshoot mixing is a critical ingredient of stellar structure models but is treated in most cases by ad hoc extensions of the mixinglength theory for convection. Advanced theories that are both more physical and numerically treatable are needed.
Aims. Convective flows in stellar interiors are highly turbulent. This poses a number of numerical challenges for the modelling of convection in stellar interiors. We included an effective turbulence model in a 1D stellar evolution code in order to treat nonlocal effects within the same theory.
Methods. We used a turbulent convection model that relies on the solution of second order moment equations. We implemented this into a stateoftheart 1D stellar evolution code. To overcome a deficit in the original form of the model, we took the dissipation due to buoyancy waves in the overshooting zone into account.
Results. We compute stellar models of intermediate mass mainsequence stars of between 1.5 and 8 M_{⊙}. Overshoot mixing from the convective core and modified temperature gradients within and above it emerge naturally as a solution of the turbulent convection model equations.
Conclusions. For a given set of model parameters, the overshooting extent determined from the turbulent convection model is comparable to other overshooting descriptions, the free parameters of which had been adjusted to match observations. The relative size of the mixed cores decreases with decreasing stellar mass without additional adjustments. We find that the dissipation by buoyancy waves constitutes a necessary and relevant extension of the turbulent convection model in use.
Key words: convection / turbulence / stars: evolution / stars: interiors
© F. Ahlborn et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model.
This Open access funding provided by Max Planck Society.
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 mainsequence 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 socalled mixinglength theory (MLT; Biermann 1932; BöhmVitense 1958). Despite its simplicity – MLT is a local and timeindependent theory – it has been used very successfully over many years to model stellar interiors. With everimproving 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 socalled 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 subadiabatic 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 oxygenburning 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 Atype mainsequence stars and DAtype 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 scaleheight 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 presentday 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 Reynoldsaveraged NavierStokes 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 ‘nonlocal’. 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 nonlinearity of the NavierStokes 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, socalled 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 mainsequence 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 3equation model and the standard 1equation 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 3equation model. For completeness, we here repeat the final model equations:
where ∇ and ∇_{ad} refer to the actual and adiabatic temperature gradient, respectively. The substantial derivative is defined as . The radiative dissipation timescale is given as
Nonlocal fluxes ℱ_{a} are modelled as
for a = ω, Π, Φ. The symbols C_{D}, γ_{R}, α_{ω}, α_{Π} and α_{Φ} are parameters. Kuhfuß (1987) suggested a value of and to be compatible with MLT in the flavour of BöhmVitense (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 StefanBoltzmannconstant. 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 1equation model
In addition to the 3equation 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:
The parameter 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 1equation 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 1equation model. When applying the 1equation 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 Henyeyscheme 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 selfconsistently 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
where we chose the same parameter α_{s} as for the diffusion coefficient of entropy in Eq. (4).
The equations can also describe timedependent effects. This was demonstrated for example by Flaskamp (2003), when computing models through the core helium flash at the tip of the redgiant branch, by Wuchterl & Feuchtinger (1998) in an application to protostars and nonlinear pulsations, and by Feuchtinger (1999) for RR Lyrae stars. In this work, however, we focus on mainsequence stars that evolve on the nuclear timescale of hydrogen burning. This means that structural changes are sufficiently slow to neglect timedependent terms and immediately solve for the stationary solution of the convection equations (lefthand 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 nonlocal effects are included in the convection model, the stationary solution describes the overshooting zone selfconsistently. 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 nonlocal terms α_{ω}, α_{Π} and α_{Φ} cannot be calibrated to MLT as they describe intrinsically nonlocal 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 nonlocal case. Although both the 1 and 3equation 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 nonlocal 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 selfconsistently 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:
where the newly introduced parameter β_{s} is defined as
and
where c_{4} = c_{3}/(c_{2}c_{ϵ}) (Eq. (21) in Paper I). The parameters c_{2} = 1.92 and c_{3} = 0.3 are model parameters from Canuto & Dubovikov (1998) and c_{ϵ} is the dissipation parameter of the convection model. The buoyancy frequency is computed according to
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_{4} = 0 as c_{3} 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 and β_{s} smoothly transitions to 1, such that setting c_{3} = 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
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:
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 (we refer to Paper I for details). The final model of Λ in the convection zone reads
for ∇ < ∇_{ad}, and
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_{4} = 0 Eq. (6) solves explicitly for Λ. In stably stratified regions, the parameter takes a value of c_{4} ≈ 0.072 using the parameters c_{2} and c_{3} 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_{4} 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 hydrogenburning 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 mainsequence phase. Then we first switched to the 1 and subsequently to the 3equation model to generate a starting model for the computation with the 3equation model. As we are interested in the effects of overshoot mixing, all models shown in the following were computed, including the nonlocal terms in the 1 and 3equation 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 colourmagnitude 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 scaleheight 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 nonlocal 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 3equation model
As a representative example, we first discuss the evolution and the internal structure of a 5 M_{⊙} model applying the 3equation, nonlocal convection theory. Figure 1 shows the profiles of the TKE variable ω (panel a), the convective flux variable Π (panel b), the superadiabatic 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 nonlocal 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_{⊙} mainsequence model, calculated with the 3equation model. The dashed black line indicates the Schwarzschild boundary. The different panels show: (a) the TKE, (b) the convective flux variable, (c) the superadiabatic temperature gradient, and (d) the diffusion coefficient according to Eq. (5). The selected stellar model has a central hydrogen abundance of X_{c} = 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 3equation theory (blue lines) and the 1equation model (red lines). Panel a: compares the hydrogen profiles at an early stage on the main sequence, when X_{c} = 0.6 (solid lines), and at the end of it (X_{c} ≈ 0; dashed lines). Panel b: convective fluxes of an MLT model with diffusive overshooting (dotted black line), a 3equation model (solid blue line), and a 1equation model (solid red line). These models have been selected to have the same central hydrogen abundance of X_{c} = 0.6 and the same chemically homogeneous core size. 
Finally, panel c shows the superadiabatic 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 superadiabatic 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^{−3}) 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 subadiabatic. 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ß 1equation model, the convective flux is directly coupled to the superadiabatic gradient and the convective velocities (Eq. (4) for the 1equation model). This of course inhibits any region in which the convective flux and the superadiabatic 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 3equation model lifts the strong coupling of convective flux and superadiabatic 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 nonlocal term in the equation for Φ that supports the positive heat flux for subadiabatic 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 superadiabatic temperature gradient) is present. This allows for the outward directed transport of entropy fluctuations, which is a positive convective flux even in subadiabatic 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 3equation 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_{⊙} mainsequence model as in Fig. 1. At the formal Schwarzschild boundary the model temperature gradient has already a slightly subadiabatic 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 subadiabatic to radiative values in a rather narrow mass range. As a consequence, the model temperature gradient takes slightly superradiative 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 superradiative region and the small deviation from the radiative temperature gradient, the overshooting zone in the 3equation nonlocal 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 counterbalanced 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_{⊙} mainsequence 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 3equation nonlocal convection model. The dashed black line indicates the Schwarzschild boundary, while the dotted black line indicates the boundary of the wellmixed 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 X_{c} = 0.6. 
The temperature gradient of the 3equation 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 subadiabaticity 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 3equation nonlocal 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 1equation nonlocal 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 Kolmogorovtype 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 superadiabatic 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_{⊙} mainsequence 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_{4}, 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_{4} 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 nonlocal 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 nonlocal 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 superadiabatic gradient increases, while the extent of this region stays the same. Decreasing the value of the parameter α_{Φ} is increasing the size of the superadiabatic region, therefore reducing the size of the Deardorff layer, while the magnitude stays the same. This is expected because the nonlocal term in the Φ equation is the one that drives convection in the Deardorff layer.
3.2. Comparison to MLT
Mixinglength theory is still the most commonly used theory to describe convection in stars. Therefore, we now compare the results of the 3equation 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 3equation 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 3equation 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 3equation model results from the solution of the model equations. The lower panel of Fig. 2 shows a comparison of the convective fluxes in the 3equation 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 3equation 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 nonlocal model can be identified. This region is shown enlarged in the inset.
Finally, in Fig. 5 we compare the evolutionary tracks of the 3equation 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 3equation 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 nonlocal model shows a higher luminosity throughout the main sequence, as expected for the larger convectively mixed core.
Fig. 5. Evolutionary tracks in the HertzsprungRussell diagram of a 5 M_{⊙} model computed with MLT (yellow line), ad hoc overshooting (black line), the 3equation, nonlocal model (blue line) and the 1equation, nonlocal model (red line). The black dot marks the model selected at a central hydrogen abundance of X_{c} = 0.6. 
3.3. Comparison to the 1equation nonlocal models
In addition to the 3equation nonlocal models, we also computed stellar models in which convection was described by the 1equation model (Sect. 2.1). As before, we included the nonlocal terms. Figure 6 shows a comparison of the TKE profiles for the 1 and 3equation, nonlocal 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 1equation nonlocal model is very comparable to the 3equation nonlocal model for the same choice of the parameter α_{ω}. The overall behaviour of the 1equation nonlocal model and the new 3equation nonlocal 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 3equation model has much smaller TKE than the 1equation model. Nevertheless, the energies are still high enough to fully mix the overshooting region in the 3equation model.
Fig. 6. Comparison of the TKE in the 1equation nonlocal model and the 3equation nonlocal model with an enhanced dissipation rate of TKE in a 5 M_{⊙} mainsequence model on a logarithmic scale. The models have been selected to have the same central hydrogen abundance as the model in Fig. 1. The 1equation model computes Λ according to Eq. (10). 
In the upper panel of Fig. 2 we compare the hydrogen profiles of the 1equation model with the results from the 3equation 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 1equation model extends past the local MLT model and looks very similar to the 3equation model. Towards the end of the main sequence, the 1equation model has a smaller core than the 3equation model, leading to a slightly different slope in the hydrogen profiles. Likewise, the evolutionary track shown in Fig. 5 of the 1equation model looks very similar to the 3equation 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 1equation model to the 3equation model and an MLT model including ad hoc overshooting. As for the 3equation 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 3equation 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 1equation 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 1equation model to have a nearly adiabatic overshooting zone, as compared to the nearly radiative overshooting zone in the 3equation model.
4. Nonlocal convection for varying initial masses
We computed stellar models in a mass range of 1.5–8 M_{⊙}, using the 3equation nonlocal 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ß 1equation model to compare the results of the 3equation 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 3equation model is larger than the convective core from an MLT model, as expected. This shows that when applying the 3equation model, an overshooting zone emerges across the whole mass range investigated. Comparing the 3equation model to the 1equation model and the ad hoc overshoot model, the mixed core sizes show good qualitative agreement. We repeat that this is achieved without finetuning 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 3equation Kuhfuß models. For low stellar masses, the results show larger discrepancies. Stellar models applying the 3equation nonlocal model have the smallest cores, and the core size decreases faster with decreasing stellar mass than in the 1equation 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 3equation 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 X_{c} = 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:
When using MLT, the convective velocities were computed as
(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 nonlocal 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 3equation 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, v_{c,iso}, is plotted. The dotted lines indicate a linear fit to the logarithmic data. We note that the data of the 1 and 3equation 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 3equation models
As discussed in Sect. 2.1 the 1equation model is a simplification of the 3equation model, for which Kuhfuß (1987) assumed that the convective flux is proportional to the superadiabatic temperature gradient and the squareroot 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 3equation 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 3equation 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 nonlocal 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 3equation models. As soon as the temperature gradient becomes subadiabatic, which happens already well within the Schwarzschild boundary in the 3equation 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 3equation model is probably the temperature stratification that results from the solution of the model equations. Due to the coupling of the convective flux to superadiabatic gradient and convective velocities in the 1equation 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 superadiabatic 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 subadiabatic zone in the convective region for the case of the 3equation 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
where u and l are a typical velocity and length scale of the convective flow. The radiative diffusivity is defined as
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
The ratio of the Peclet numbers obtained for the 1 and 3equation models for a 5 M_{⊙} mainsequence model is shown in Fig. 9. In the bulk of the convection zone, the 1 and 3equation models have very similar Peclet numbers, which means the transport of energy by convection behaves very comparably. In the overshooting zone, however, the 1equation model has a Peclet number that is up to seven times higher than that of the 3equation model. This indicates that convection as described by the 1equation model is much more efficient in the overshooting zone than when described by the 3equation model.
Fig. 9. Ratio of the Peclet numbers for the 1 and 3equation models for a 5 M_{⊙} mainsequence model. 
This change in efficiency also impacts the temperature gradient. Following the approximation of the convective flux in the 1equation model, the convective flux scales as . Using the ratio of the Peclet numbers, one can therefore write a Pecletscaled convective flux,
to mimic the convective flux in the 3equation 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 3equation model are visualised by a green dotted and solid line, respectively.
Fig. 10. Comparison of convective fluxes and temperature gradients obtained from the Pecletscaling and the 1 and 3equation models. Upper panel: convective flux as a function of fractional mass. The red and blue lines show the results obtained from the 1 and 3equation 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 3equation models. The dashed green line shows the temperature gradient of the 1equation 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 3equation model is due to the different penetration depths of TKE and convective flux in the 3equation 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 3equation model. Following the terminology proposed by Zahn (1991), in the 1equation model the overshooting zone is best described by subadiabatic penetration while the more inefficient convection in the 3equation model concerns overshooting of chemical element distributions only.
Considering the chemical mixing, both models result in a more steplike 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 3equation model, the application of the 1equation model seems to be sufficient to obtain the chemical structure from a nonlocal convection model. The parameter α_{ω} can be tuned to obtain the correct size of the convective core. However, the stratification obtained from the 1equation 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 HertzsprungRusselldiagram of massive stars, and Maeder & Mermilliod (1981) in relation to the HertzsprungRusselldiagram of open clusters. In the latter case, isochrones derived from stellar models that include core overshooting match the morphology of the turnoff 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 mainsequence 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 gmode 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 pmode 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_{⊙} mainsequence star, which is in agreement with the temperature gradient obtained from the 3equation model but disagrees with the 1equation model. However, they point out that this result is only obtained for this single Btype star and might not be generalisable for all Btype stars.
With the increase in computational resources in recent years, more and more multidimensional hydrosimulations 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 3equation model predicts mixed core masses similar to the GARSTEC models that include this geometric cutoff (Fig. C.1), this indicates that our 3equation 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 nonlocal 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 wellmixed region increases with increasing stellar mass. Both results are in agreement with our findings employing the 3equation 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 1equation 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 3equation 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 mainsequence models of intermediate mass mainsequence 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 3equation 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ß 3equation 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_{⊙} mainsequence model. We find that the TKE extends beyond the formal Schwarzschild boundary of convective neutrality. This is the result of the nonlocal 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ß 3equation 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 3equation model to the simplified 1equation 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 subadiabatic 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 subadiabatic to a radiative value, exhibiting a small region with a superradiative 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 3equation model, the overshooting zone of the 1equation 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 1equation model (Eq. (4)). The approximation Eq. (4) does not allow for a Deardorff layer, as the convective flux is proportional to the superadiabatic gradient. From a theoretical point of view, the existence of a Deardorff layer can be attributed to the nonlocal 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 3equation model, to capture the temperature structure in the overshooting zone more accurately. The comparison of the Peclet numbers of the 1 and 3equation models further shows that the mostly radiative temperature gradient in the overshooting zone of the 3equation model can be explained by the reduced TKE or velocities, respectively, compared to the 1equation 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 1equation 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 3equation models over a broader mass range show qualitative agreement with other overshoot descriptions. For given values of α_{ω} and f_{OV}, this is achieved without finetuning 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 1equation and the 3equation 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 3equation model shows a reduced amount of overshooting compared to the 1equation model, while the 3equation 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 3equation 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 finetuned 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 3equation model from Kuhfuß (1987) shows that the resulting stellar structure depends sensitively on the details of the convection model. Although the original 1 and 3equation 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 3equation model compares very well with the 1equation 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 1equation model describes the convective core similarly well as the 3equation 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 3equation model need to be used. A selfconsistent 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 3equation 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 hydrosimulations 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 3equation 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 P29172N and P33140N and support from European Research Council (ERC) Synergy Grant WHOLESUN #810218
References
 André, J. C., de Moor, G., Lacarrère, P., & Du Vachat, R. 1976, J. Atmos. Sci., 33, 476 [CrossRef] [Google Scholar]
 Angelou, G. C., Bellinger, E. P., Hekker, S., et al. 2020, MNRAS, 493, 4987 [NASA ADS] [CrossRef] [Google Scholar]
 Arnett, W. D., Meakin, C., Viallet, M., et al. 2015, ApJ, 809, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Biermann, L. 1932, ZAp, 5, 117 [NASA ADS] [Google Scholar]
 BöhmVitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
 Bressan, A. G., Chiosi, C., & Bertelli, G. 1981, A&A, 102, 25 [Google Scholar]
 Canuto, V. M. 1992, ApJ, 392, 218 [Google Scholar]
 Canuto, V. M. 1993, ApJ, 416, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V. M. 1997, ApJ, 482, 827 [Google Scholar]
 Canuto, V. M. 2011, A&A, 528, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Canuto, V. M., & Dubovikov, M. 1998, ApJ, 493, 834 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. L., & Gigas, D. 1992, ApJ, 389, L87 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. L., & Sofia, S. 1989, ApJ, 336, 1022 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, K. L., & Sofia, S. 1996, ApJ, 466, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Chou, P. Y. 1945, Q. Appl. Math., 3, 38 [CrossRef] [Google Scholar]
 Claret, A., & Torres, G. 2019, ApJ, 876, 134 [Google Scholar]
 Constantino, T., & Baraffe, I. 2018, A&A, 618, A177 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deardorff, J. W. 1966, J. Atmos. Sci., 23, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Brandão, I., Silva Aguirre, V., et al. 2016, A&A, 589, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [Google Scholar]
 Feuchtinger, M. U. 1999, A&A, 351, 103 [NASA ADS] [Google Scholar]
 Flaskamp, M. 2003, PhD Thesis, MaxPlanckInstitut fürAstrophysik, Technische Universität München [Google Scholar]
 Freytag, B., Ludwig, H. G., & Steffen, M. 1996, A&A, 313, 497 [NASA ADS] [Google Scholar]
 Gilet, C., Almgren, A. S., Bell, J. B., et al. 2013, ApJ, 773, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Grevesse, N., & Noels, A. 1993, in Origin and Evolution of the Elements, eds. N. Prantzos, E. VangioniFlam, & M. Casse, 15 [Google Scholar]
 Henyey, L. G., Forbes, J. E., & Gould, N. L. 1964, ApJ, 139, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841 [Google Scholar]
 Higl, J., Siess, L., Weiss, A., & Ritter, H. 2018, A&A, 617, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Higl, J., Müller, E., & Weiss, A. 2021, A&A, 646, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Rheinhardt, M., Brandenburg, A., et al. 2017, ApJ, 845, L23 [CrossRef] [Google Scholar]
 Keller, L. V., & Friedmann, A. A. 1925, in Proceedings of the First International Congress of Applied Mechanics, Delft, ed. J. B. C. Biezeno, 395 [Google Scholar]
 Kippenhahn, R., Weigert, A., & Hofmeister, E. 1967, Methods. Comput. Phys., 7, 129 [NASA ADS] [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution, Astronomy and Astrophysics Library (Berlin Heidelberg: SpringerVerlag) [CrossRef] [Google Scholar]
 Kuhfuß, R. 1986, A&A, 160, 116 [NASA ADS] [Google Scholar]
 Kuhfuß, R. 1987, PhD Thesis, MaxPlanckInstitut für Astrophysik, Technische Universität München [Google Scholar]
 Kupka, F. 1999a, ApJ, 526, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Kupka, F. 1999b, in Stellar Structure: Theory and Test of Connective Energy Transport, eds. A. Gimenez, E. F. Guinan, & B. Montesinos, 173, 157 [NASA ADS] [Google Scholar]
 Kupka, F. 2007, in Convection in Astrophysics, eds. F. Kupka, I. Roxburgh, & K. L. Chan, 239, 92 [NASA ADS] [Google Scholar]
 Kupka, F., & Montgomery, M. H. 2002, MNRAS, 330, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2007a, in Convection in Astrophysics, eds. F. Kupka, I. Roxburgh, & K. L. Chan, 239, 80 [NASA ADS] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2007b, in Convection in Astrophysics, eds. F. Kupka, I. Roxburgh, & K. L. Chan, 239, 83 [NASA ADS] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2007c, in Convection in Astrophysics, eds. F. Kupka, I. Roxburgh, & K. L. Chan, 239, 86 [NASA ADS] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2008, in The Art of Modeling Stars in the 21st Century, eds. L. Deng, & K. L. Chan, 252, 463 [NASA ADS] [Google Scholar]
 Kupka, F., & Muthsam, H. J. 2017, Liv. Rev. Comput. Astrophys., 3, 1 [Google Scholar]
 Kupka, F., Zaussinger, F., & Montgomery, M. H. 2018, MNRAS, 474, 4660 [NASA ADS] [CrossRef] [Google Scholar]
 Kupka, F., Ahlborn, F., & Weiss, A. 2022, A&A, 667, A96 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Langer, N., El Eid, M. F., & Fricke, K. J. 1985, A&A, 145, 179 [NASA ADS] [Google Scholar]
 Li, Y. 2012, ApJ, 756, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Y. 2017, ApJ, 841, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Y., & Yang, J. Y. 2007, MNRAS, 375, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A., & Mermilliod, J. C. 1981, A&A, 93, 136 [NASA ADS] [Google Scholar]
 Magic, Z., Serenelli, A., Weiss, A., & Chaboyer, B. 2010, ApJ, 718, 1378 [NASA ADS] [CrossRef] [Google Scholar]
 Meakin, C. A., & Arnett, D. 2007, ApJ, 667, 448 [NASA ADS] [CrossRef] [Google Scholar]
 Michielsen, M., Pedersen, M. G., Augustson, K. C., Mathis, S., & Aerts, C. 2019, A&A, 628, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michielsen, M., Aerts, C., & Bowman, D. M. 2021, A&A, 650, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mombarg, J. S. G., Van Reeth, T., Pedersen, M. G., et al. 2019, MNRAS, 485, 3248 [NASA ADS] [CrossRef] [Google Scholar]
 Montgomery, M. H., & Kupka, F. 2004, MNRAS, 350, 267 [NASA ADS] [CrossRef] [Google Scholar]
 Muthsam, H. J., Goeb, W., Kupka, F., Liebich, W., & Zoechling, J. 1995, A&A, 293, 127 [Google Scholar]
 Muthsam, H. J., Göb, W., Kupka, F., & Liebich, W. 1999, New A, 4, 405 [NASA ADS] [CrossRef] [Google Scholar]
 Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nat. Astron., 5, 715 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrinferni, A., Cassisi, S., Salaris, M., & Castelli, F. 2004, ApJ, 612, 168 [Google Scholar]
 Renzini, A. 1987, A&A, 188, 49 [NASA ADS] [Google Scholar]
 Roxburgh, I. W. 1978, A&A, 65, 281 [NASA ADS] [Google Scholar]
 Roxburgh, I. W. 1992, A&A, 266, 291 [NASA ADS] [Google Scholar]
 Roxburgh, I. W., & Kupka, F. 2007, in Convection in Astrophysics, eds. F. Kupka, I. Roxburgh, & K. L. Chan, 239, 98 [NASA ADS] [Google Scholar]
 Saslaw, W. C., & Schwarzschild, M. 1965, ApJ, 142, 1468 [Google Scholar]
 Shaviv, G., & Salpeter, E. E. 1973, ApJ, 184, 191 [Google Scholar]
 Stellingwerf, R. F. 1982, ApJ, 262, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblay, P. E., Ludwig, H. G., Freytag, B., et al. 2015, ApJ, 799, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, J. S. 1986, J. Fluid Mech., 173, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Meakin, C., Arnett, D., & Mocák, M. 2013, ApJ, 769, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Viallet, M., Meakin, C., Prat, V., & Arnett, D. 2015, A&A, 580, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weiss, A., & Schlattl, H. 2008, Ap&SS, 316, 99 [CrossRef] [Google Scholar]
 Wuchterl, G. 1995, Comput. Phys. Commun., 89, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Wuchterl, G., & Feuchtinger, M. U. 1998, A&A, 340, 419 [NASA ADS] [Google Scholar]
 Xiong, D.R. 1978, Chin. Astron., 2, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Xiong, D.R. 1986, A&A, 167, 239 [NASA ADS] [Google Scholar]
 Xiong, D. R., & Deng, L. 2001, MNRAS, 327, 1137 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J. P. 1991, A&A, 252, 179 [Google Scholar]
 Zhang, Q. S. 2016, ApJ, 818, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, Q. S., & Li, Y. 2012, ApJ, 746, 50 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Comparison of 1 and 3equation models
Figure A.1 shows the TKE in the 1 and 3equation model on a linear scale. Even though the TKE of the 3equation model has approximately the same extent as the 1equation 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 3equation model are very comparable.
Fig. A.1. Comparison of the TKE of the 1equation and 3equation nonlocal models in a 5 M_{⊙} mainsequence model with limited dissipation length scale, Λ, for the 3equation model. 
Appendix B: Parameter dependence
The model for the modification of the dissipation length scale comes with a number of parameters (c_{2}, c_{3}, 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_{4}. Given the dissipation parameter of the 3equation model, we find c_{4} ≈ 0.072 while for the dissipation parameter of the Canuto & Dubovikov (1998) model c_{4} ≈ 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_{4} to some extent cancels out in the calculation of Λ, as c_{ϵ} appears in the denominator of c_{4} = c_{3}/(c_{2}c_{ϵ}) and in the numerator of Λ = c_{ϵ}ω^{3/2}/ϵ. Hence, if c_{4} is adjusted according to c_{4} = c_{3}/(c_{2}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_{4} 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_{4} 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 c_{4} in a 3equation nonlocal 5 M_{⊙} mainsequence 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 1equation nonlocal 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 3equation nonlocal 5 M_{⊙} mainsequence 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_{4} will change, while in Fig. B.1 only c_{4} is changed. The effects of changing c_{4} 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 C_{D} in a 3equation nonlocal 5 M_{⊙} mainsequence 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 3equation nonlocal 5 M_{⊙} mainsequence 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 scaleheight at the Schwarzschild boundary to define the overshooting distance. In the description of Freytag et al. (1996) the diffusion coefficient is computed according to
where f_{OV} is an adjustable parameter and z is the radial distance to the Schwarzschild boundary. The pressure scaleheight 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:
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:
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:
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 3equation 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 3equation 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_{⊙} mainsequence model, calculated with the 3equation model. The dashed black line indicates the Schwarzschild boundary. The different panels show: (a) the TKE, (b) the convective flux variable, (c) the superadiabatic temperature gradient, and (d) the diffusion coefficient according to Eq. (5). The selected stellar model has a central hydrogen abundance of X_{c} = 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 3equation theory (blue lines) and the 1equation model (red lines). Panel a: compares the hydrogen profiles at an early stage on the main sequence, when X_{c} = 0.6 (solid lines), and at the end of it (X_{c} ≈ 0; dashed lines). Panel b: convective fluxes of an MLT model with diffusive overshooting (dotted black line), a 3equation model (solid blue line), and a 1equation model (solid red line). These models have been selected to have the same central hydrogen abundance of X_{c} = 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_{⊙} mainsequence 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 3equation nonlocal convection model. The dashed black line indicates the Schwarzschild boundary, while the dotted black line indicates the boundary of the wellmixed 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 X_{c} = 0.6. 

In the text 
Fig. 4. Dissipation length scale Λ (a) and dissipation rate (b) in a 5 M_{⊙} mainsequence 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 HertzsprungRussell diagram of a 5 M_{⊙} model computed with MLT (yellow line), ad hoc overshooting (black line), the 3equation, nonlocal model (blue line) and the 1equation, nonlocal model (red line). The black dot marks the model selected at a central hydrogen abundance of X_{c} = 0.6. 

In the text 
Fig. 6. Comparison of the TKE in the 1equation nonlocal model and the 3equation nonlocal model with an enhanced dissipation rate of TKE in a 5 M_{⊙} mainsequence model on a logarithmic scale. The models have been selected to have the same central hydrogen abundance as the model in Fig. 1. The 1equation 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 X_{c} = 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, v_{c,iso}, is plotted. The dotted lines indicate a linear fit to the logarithmic data. We note that the data of the 1 and 3equation 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 3equation models for a 5 M_{⊙} mainsequence model. 

In the text 
Fig. 10. Comparison of convective fluxes and temperature gradients obtained from the Pecletscaling and the 1 and 3equation models. Upper panel: convective flux as a function of fractional mass. The red and blue lines show the results obtained from the 1 and 3equation 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 3equation models. The dashed green line shows the temperature gradient of the 1equation model computed from the scaled convective flux shown in the upper panel. 

In the text 
Fig. A.1. Comparison of the TKE of the 1equation and 3equation nonlocal models in a 5 M_{⊙} mainsequence model with limited dissipation length scale, Λ, for the 3equation 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 c_{4} in a 3equation nonlocal 5 M_{⊙} mainsequence 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 3equation nonlocal 5 M_{⊙} mainsequence 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 C_{D} in a 3equation nonlocal 5 M_{⊙} mainsequence 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 3equation nonlocal 5 M_{⊙} mainsequence 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.