© ESO 2018

1 Introduction

Convective overshoot is a long-standing and thorny problem in stellar physics. Shaviv & Salpeter (1973) were the first to try to compute the overshooting above a convective core. They computed the slowdown of a rising convective bubble above the boundary, as given by the Schwarzschild criterion, assuming that the temperature gradient is slightly sub-adiabatic. However, they did not discuss the motion of the downward material, which remained unexplained. Their work was followed by many others, which are critically discussed in Renzini (1987). Due to the large uncertainties of the theories available at that time, the overshooting distance was quickly considered as a free parameter, taken proportional to the pressure scale height at the classical surface of the core (i.e. as given by the classical Schwarzschild criterion). In addition, there is still no agreement on the value of the temperature gradient to be adopted in the overshooting layers. Some codes assume an adiabatic gradient while others prefer the radiative one, so that the extent obtained through fitting of evolutionary tracks using observations depends on this choice.

From the late seventies, the group led by Xiong (see Xiong & Deng 2001 and references therein) developed a theory of convection based on a Reynolds stress approach and were able to incorporate it into an evolution code (see also the more recent work by Kupka & Montgomery 2002). This type of approach is quite powerful and widely used in geophysics but it provides results that are very sensitive to the adopted closure model (and associated free parameters) as well as to the domain of applicability of the closure model. Consequently, only a handful of evolutions have been computed using it. A very complete review giving the present state of the problem and including a clear historic of the modelling of convection is given by Kupka & Muthsam (2017). Recently, based on 3D hydrodynamical simulations of convection in massive stars (e.g. Meakin & Arnett 2007b; Arnett et al. 2009; Arnett & Meakin 2011; Viallet et al. 2013); Arnett et al. (2015) proposed a model based on a Reynolds-averaged Navier-Stokes equations (RANS) analysis to propose a set of equations, which they call the 321D algorithm, to include in 1D evolution codes. This promising approach presents the advantage of solving the problem and of dealing with convective boundary layers without making an assumption concerning the flow geometry.

On the other hand, Schmitt et al. (1984) introduced the theory of turbulent plumes in astrophysics. It was applied to stellar envelopes by Rieutord & Zahn (1995) and adapted to convective cores by Lo & Schatzman (1997). In these works, the spherically-averaged equations are ignored with the exception of the continuity equation, which is used to modify the advection term given by Taylor’s hypothesis. This assumption, also known as the turbulent entrainment hypothesis, states that the radial inflow of matter is proportional to the central vertical velocity inside a plume (see Turner 1986, for details). This approach has several consequences, the main one being that the coefficient of proportionality in Taylor’s hypothesis is to be specified. It is commonly considered as a universal value but has only been constrained in geophysical flows or by laboratory experiments. Therefore, its value in the stellar context remains subject to uncertainties. One also has to assume that the temperature gradient in the overshooting layers is everywhere close to the adiabatic value. Finally, one is left with a free parameter that is the number of plumes and one has to assume that all the energy, which is not carried out by radiation, is carried out in the plumes. However, this last hypothesis overlooks the fact that the energy generation rate is, with a very good approximation, the same at all points of a spherical surface, as the convective fluctuation of the temperature and the density are very small and one needs to explain how the energy flux would be concentrated into the plumes. Moreover, this assumption cannot be maintained as soon as the equations averaged over a spherical surface are taken into account. In addition, it is also fundamental to consider the kinetic energy flux because, as will be shown without any assumption in this article, if it is neglected, overshoot is impossible. We must also take into account that this flux has different signs in the up- and down-flows (while the convective flux has the same sign in both regions) and that it is non-negligible at every point on a spherical surface. This holds even in the local mixing-length theory (LMLT) for which only its average value vanishes. Therefore, an extra equation, which enables us to compute the super-adiabatic gradient at every point of the convective region, must be considered. It allows us to remove the adiabatic hypothesis in the overshooting layers.

Consequently, we will take into account the equations of continuity, mechanical energy conservation, and thermal energy conservationintegrated over the surface of a sphere. This set of equations will be supplemented by the same equations but integrated over the “horizontal” section of a plume. Using the above-mentioned equations enables the system to be closed and it is no longer necessary to use Taylor’s hypothesis. We have, however, to make an approximation, present in all the theories of convection, because we do not know the turbulent energy spectrum: we have to give a formula for the kinetic energy dissipation rate, which is a phenomenon taking place at the end of the turbulent spectrum (i.e. at high wave numbers) using variables characteristic of the low wave numbers. Therefore, we must choose an expression for the kinetic energy dissipation rate that is only an order of magnitude and that introduces a parameter similar to the mixing length. As a consequence, the theory we propose is essentially a fully coherent non-local mixing-length theory, which offers a self-consistent method to compute the structure of the overshooting layers. Unfortunately, the price to pay is that we are left with an eight-order system (plus essentially two algebraic equations). We also note that recent results obtained by 3D simulations of deep stellar convection (Arnett et al. 2009) were able to justify the formula giving the kinetic energy dissipation rate that has been used for many years for the above mentioned reason. However, they show that the mixing length is not the pressure scale height, as usually assumed, but the size of the larger eddies, which is close to the size of the convection zone (see also Sect. 2 in Arnett et al. 2015).

The theory we will develop can in principle be applied to any convective region, however in the outermost layers of convective envelopes, where the temperature gradient is significantly non-adiabatic and where the convective cells observed at the surface are formed, the plume theory is no longer applicable. Therefore to be able to apply the plume theory to stellar envelopes, it will be necessary to use numerical simulations as initial conditions on top of the quasi-adiabatic region. Alternatively, it can be applied if the solutions become quickly independent of the upperboundary, as suggested by Rieutord & Zahn (1995) for a somewhat idealized problem. On the other hand, nothing seems to be in contradiction with the application of the theory in convective cores. This motivated our choice to develop the theory for a convective core with gas moving upwards in the plumes. Nevertheless, we will indicate how easy it is to modify the equations for applications to envelopes.

Let us recall that a convective zone boundary is the point where both the radial component of the velocity is cancelled out (V _r = 0) and radiation carries out all the energy (L = L_R, where L and L_R are the total and radiative luminosities, respectively). We now still have to define precisely what the overshooting region is. We would like to define it as the layers of the convective core laying above the core boundary predicted by the LMLT. The boundary given by the LMLT satisfies the two conditions V _r = 0 and L = L_R. It is worth recalling that it is mathematically not allowed to find that boundary using an interpolation of either ∇_R − ∇_a (Schwarzschild criterion, where ∇_R and ∇_a are the radiative and adiabatic gradients, respectively) or ∇_R −∇_L (this is the Ledoux criterion, see Gabriel et al. 2014, their Eq. (10), for a definition of ∇_L ) when that function is discontinuous at the surface of the core (see Gabriel et al. 2014). When a better theory is used, we see that the conditions L = L_R and ∇_a = ∇_R are not fulfilled at the same point and therefore these two conditions are not equivalent to define the top of the convective core without overshooting. We will thus take the bottom of the overshooting layers at the point where L = L_R and V _r ≠ 0, as this is the condition that is physically meaningful while the other one, ∇_a = ∇_R, is equivalent only in the frame of the LMLT. Meakin & Arnett (2007b) give another definition of the convection zone boundaries. They are given by the points where the bulk Richardson number becomes larger than 1∕4. Even if there is some uncertainty over what expression of that number should be used (e.g. Meakin & Arnett 2007b; Arnett et al. 2015; Cristini et al. 2017), when the convective core is receding our definition and theirs are practically equivalent. The rise of the potential wall produced by the molecular weight gradient, which gives rise to a fast increase of the Brunt-Väisälä, will quickly stop raising matter already strongly slowed down in the layers with a sub-adiabatic temperature gradient, which startsfrom below the overshooting region as we define it. In contrast, the definition used by Meakin & Arnett (2007b) will be very useful when the convection core expends, all the more so as they give an expression for the entrainment velocity. We will come back to this point in Sect. 8 devoted to expanding cores.

The paper is organized as follows: Sect. 2 introduces the basic equations of the problem while Sects. 3 and 4 discuss the mean-field and plume-averaged equations, respectively. The system of equations to be integrated is summarized in Sect. 5 while a discussion on the treatment of parameters is provided in Sect. 6 and the boundary conditions in Sect. 7. In Sect. 8, the case of an expanding convective core is discussed. Section 9 is dedicated to the conclusions.

2 Equations of the problem

In the following, we will make use of two usual hypotheses for convection, namely: stationarity and vanishing pressure fluctuations as justified in Massaguer & Zahn (1980). These two simplifying assumptions are indeed open to criticism. First, the hypothesis of vanishingpressure fluctuations is not sustained by numerical simulations. Second, as could be expected, the hypothesis of stationarity contradicts a basic characteristic of convection. At its surface, waves are irregularly produced with time in the stable region causing extra mixing. Arnett & Meakin (2011) discovered that during late phases of stellar evolution numerical simulations present spatial and time variations in the velocity field, which lead to oscillations similar to those observed in irregular and semi irregular variables. These oscillations caused by luminosity fluctuations directly associated with turbulent convective cells are inherently nonlinear, that is to say they are not detectable through the usual linear studies, and are produced by what they call the τ-mechanism (see also Smith & Arnett 2014). Also, our theory predicts one kind of plume, which represents the average of all plumes, but strong plumes will lead to more mixing and this mixing is irreversible.

Moreover, we will consider chemically homogeneous convective cores. This is justified for a shrinking convective core as encountered in massive main-sequence stars. Indeed, for those stars the core mass grows quickly enough so that when it reaches its maximum value, the star is still chemically homogeneous with a good approximation. In that case, the total mass of the core (including the overshoot region) decreases during the evolution as does the core mass obtained from the LMLT. The case of an expanding convective core with a discontinuity of chemical composition at its surface will be briefly considered in Sect. 8.

The equations as given in this section are taken from Ledoux & Walraven (1958) in which any variable y is split into two parts (1)

where is the mean value of y, averaged over a spherical surface. We will consider the variables y = {p, T, ρ, ρV _i}, where p is the pressure, T the temperature, ρ the density, and ρV _i the mass flux (V being the convective velocity) with . Since we consider a hypothetic stationary convection, we have here to consider spatial averages only and we do not, as in the more realistic numerical simulations, have to distinguish between spatial and time averages.

2.1 Equation of state

As already mentioned, the pressure fluctuations are neglected (Δp = 0) and the convective zone is considered to be chemically homogeneous. Thus, one has (2)

or equivalently (3)

Rigorously, instead of Eq. (2) one should write (4)

where X_i is the fraction by weight of nuclei of type i. One should also add to the kinetic equations for nuclear reactions a turbulent diffusion term (which is not well known). This would make the problem much more complex while introducing a non-negligible effect only in a very tiny layer at the surface of the core where the lifetime of the convective elements becomes of the same order as the thermal diffusion time. Therefore, we will ignore it and consider only Eq. (2).

Let us also recall the well-known formula for the adiabatic gradient (see for instance Landau & Lifchitz 1967, Eq. 16.3) (5)

where C_p is the heat capacity at constant pressure.

2.2 Continuity equation

For stationary convection, if we assume as is usually done that V _ϕ is independent of ϕ, the mass conservation equation can be written (6)

where V _r, V _θ, V _ϕ are the three components of the convective velocity (V).

From Eq. (6), it follows that there is no net mass flux over a spherical surface, that is, that . Therefore, one has (7)

As , we prefer to use ρV as a fundamental variable rather than V.

2.3 Mechanical energy conservation

The equation of motion, assuming Δp = 0, is (8)

The coefficient α as dissipation factor was first introduced by Prandl (see references in Arnett et al. 2015) but the same formula was also introduced independently by Cowling (1935). The last term, which comes from the dissipation by viscosity of kinetic energy at the low end of the energy cascade, is defined by using 1∕τ = α|V _r|∕l_LMLT for the characteristic time (τ) of kinetic energy dissipation into heat. Most often, l_LMLT is approximated by l_LMLT ≈ H_p, as has been the case since Böhm-Vitense. In contrast, the recent work by Arnett et al. (2015) suggests the use of 1∕τ =|V |∕l_d with l_d ≈ 0.8l_CZ (see their Table 1), where l_CZ is the thickness of the convection zone. This is an important change as it makes convection much more efficient and reduces the departure from adiabatic gradient. We momentarily keep the old definition of “the mixing length” to allow for the following discussion, but we will adopt Arnett et al.’s definition from Eq. (10) on.

We also assumed that the static equilibrium equation is given by (9)

with p the gas pressure. The latter equation is verified to a good accuracy in the deep interior of stars where V ≪c (c being the sound speed).

The last term of Eq. (8) represents the divergence of the viscous stress tensor where l_MLT is the mixing length. The coefficient α is introduced in order to recover the equations of the LMLT for the local solution of the linearized equations for convection (see Gabriel et al. 1974, 1975; Gabriel 1996, for details). For instance, to obtain the equation of Henyey et al. (1965), we must use the value α = 8∕3. However, there is no reason to also adopt this value in a convective core. Nuclear reactions are going on in the core and also the energy losses are not necessarily isotropic as assumed in the outer layers.

However, to simplify the notations, we have chosen here to set α equal to unity. It is equivalent to define a new mixing length l = l_MLT∕α and to define the characteric dissipation time of kinetic energy by τ = l∕|V _r|. This means that here the mixing length must not be interpreted as the mean travel distance of a turbulent eddy but as the length that allows us to define the correct characteristic time τ. It was already so in the formalism adopted by Gabriel et al. (1974, 1975) and Gabriel (1996).

To go further, we multiply Eq. (8) by V to get the equation of kinetic energy conservation (10)

where we introduced our above-mentioned definition of l. The parameter l used in the above equation is related to that defined by Arnett et al. (2009, 2015) by (11)

They also show that this expression for the kinetic energy dissipation is obtained because that term has the Kolmogorov form.

We also introduce the kinetic energy flux of turbulence F _K = ρV ²V∕2 and the expression adopted for the kinetic energy dissipation rate into heat ρϵ₂ = ρV ²∕τ = ρV ²|V _r|∕l . It provides for the linearized equations of convection a local stationary solution identical to the LMLT ones (see Gabriel et al. 1974, 1975; Gabriel 1996).

Equation (10) averaged over a spherical surface then gives (12)

2.4 Thermal energy conservation

Let us now consider the thermal energy equation (13)

where F_R is the radiative flux, ϵ is the energy generation rate, and F_H is the enthalpy flux defined as (14)

with S being the entropy. Using Eq. (6), the divergence of the enthalpy flux can be written as (15)

where the convective flux, F_C, is defined by (16)

The isentropic assumption has not been made so that Eq. (15) can be considered as a generalization of the energy equation used in Rieutord & Zahn (1995).

Using Eqs. (10) and (15), Eq. (13) finally becomes (17)

Averaging over a spherical surface, we recover the usual equation of thermal energy conservation (18)

or equivalently (19)

where are the average of the kinetic, convective, radiative luminosities, and the total luminosity. By taking the pressure fluctuations into account and by introducing a new term that implies that energy generation also supplies the work needed to maintain the convective flow, Arnett et al. (2009) (Sect. 3) and Meakin & Arnett (2007b) (Appendix A) obtained a different expression for Eq. (17). Also if the driving of turbulence may be local, turbulent flows redistribute turbulence more uniformly, implying non-zero turbulent flux.

At this point, it is worth mentioning that convection is more efficient in rising plumes rather than in the downward flow because the convective flux is larger (except if the surface area occupied by up- and down-flows is the same) and also because of the positive kinetic energy flux (while it is negative in the down-moving gas). Therefore, rising plumes carry out more energy than produced inside the plumes since ϵ is, with a good approximation, the same everywhere on a spherical surface. For instance, if in a thin layer the contribution of the radial component of the fluxes to ∇⋅(F_K + F_C) averaged over a plume is larger than its value averaged over a spherical surface, the difference must be counterbalanced by the contribution of the horizontal components directed inwards to the plume to prevent the gas of that layer from cooling down. This horizontal flux is necessary to allow a stationary state. This remark holds also for the LMLT even if, in that case, the flux of kinetic energy from the down to the up-flow allows a kinetic energy flow averaged over a spherical surface equals to zero.

Subtracting Eq. (18) from Eq. (17) we further get (20)

In the framework of the LMLT, the first term of the latter equation reads (21)

and the kinetic energy flux (F_K) is ignored. However, F_K is everywhere different from zero and only its average over a spherical surface vanishes. Therefore, it is not justified to assume that ΔF _K = 0 and we have to add an extra term. To this end, we will assume that (22)

The computation of the fluctuation of the convective and kinetic energy fluxes would normally require knowledge of higher-order fluctuations of the temperature and of the impulsion. Therefore, Eq. (22) may be considered as the second closure hypothesis of the problem, the first one being the expression chosen for ϵ₂ . Assuming that the turbulent cascade follows the Kolmogorov law proportional to the third power of the velocity, Arnett et al. (2009) have justified this expression for ϵ₂. Also numerical simulations do not need any closure equation as departures from averages are obtained through the simulations and all average quantities required, for instance, in Eqs. (20) to (24) of Arnett et al (2009), can be obtained through time or spatial averages.

Within the LMLT, the first term of the right hand side of Eq. (22) is introduced because, after having travelled over an average mixing length, a rising element releases its excess thermal energy into the average medium. Here, it will be interpreted as the dissipation rate of thermal energy. We add a second term in Eq. (22) because the kinetic energy must also be released into the average medium. Here, this term gives the dissipation rate of convective kinetic energy. When convection at one point is compared with the average over a spherical surface, the last term of Eq. (22) must be present. In contrast, in the LMLT, this term is absent. This can be explained as follows; the LMLT has been developed for outer layers, for which the two right-hand-side terms of Eq. (22) are of the same order of magnitude. These terms are only orders of magnitudes so that they can be gathered together into the first oneand it may be considered that the LMLT takes both terms into account.

Finally, using Eq. (20) together with Eq. (22), one gets (23)

with (24)

where we used the definitions and .

3 Mean equations

3.1 Horizontal dependence of the convective variables

To go further, we now have to specify the horizontal dependence of the convective variables. More precisely, we will separate the variables in the following way: (25)

where and are averages over a spherical surface as defined by Eq. (1).

For the mass flux, as already mentioned, its average vanishes so that (26)

where is the value of ρV at the point where h is chosen to be equal to unity. This holds also for Δρ₀ and ΔT₀.

The normalization of the horizontal function, h, gives (27)

where dΩ = sinθ dθ dϕ is the elementary solid angle, and S the surface of a sphere. Equation (27) ensures that the convective flux of matter averaged over a spherical surface is zero.

For the convective velocity, we must separate the variables with another horizontal function because the average of the velocity does not vanish. Therefore, we write (28)

The functions h and h₁ can be provided by numerical simulations. Moreover, using Reynold-averaged Navier-Stokes equations (RANS), there is no need to provide such a decomposition and subsequently to make a choice concerning their formulation (see below). Using a Taylor expansion at first order, one can also write for the convective velocity (29)

In convective cores, the Mach number (ratio between the convective velocity and the sound speed) is very small. Therefore, , because it depends on the squared Mach number. Consequently, equating Eq. (28) and Eq. (29 ), one has (30)

It is now possible to express the averaged convective quantities. Let us first consider the averaged convective flux. Using Eq. (26), it reads (31)

where we have defined (32)

For the convective velocity, using Eq. (7) together with Eq. (26), we have (33)

with (34)

It is now worth considering the averaged kinetic energy flux and more precisely its radial component. One can write (35)

To obtain Eq. (35) we have introduced a hypothesis concerning the isotropy of convection and we have written (36)

with k =3 for isotropic convection and k = 1 for purely radial convection. This is a significant advantage of using the equation of kinetic energy conservation rather than the radial component of the equation of motion as the latter takes a much more complex form when convection is anisotropic so that the previous works implicitly assume isotropy (e.g. Rieutord & Zahn 1995; Lo & Schatzman 1997). Of course, we introduce a new free parameter but at least it will be possible to check the influence of the choice of its value on the results while previous works were confined to k=3. Arnett et al. (2009) discussed the anisotropy far away from the boundaries (see their Eq. (10) in Sect. 2.4) and found that 1.5 < k < 2.28.

Using Eq. (35), the averaged kinetic energy flux immediately reads (37)

where we have introduced the notation (38)

and (39)

The averaged luminosity associated with the radial component of the averaged kinetic energy flux can then be written as .

Finally, one can express the kinetic energy dissipation rate into heat using the isotropic parameter (Eq. (36)), such as (40)

Therefore, one has (41)

with (42)

To be able to compute the coefficients G_i and in the following section the equations integrated over one plume, we must give an explicit form to the angular function h(θ, ϕ). We will first assume that there are N identical plumes, which have a conical shape of opening θ_M(r). Moreover, within each plume, we assume (43)

This means that we take h = 1 on the axis of the cone associated with each plume. Equation (43) is different from the usual one, which assumes a Gaussian profile that fits well the laboratory observations (e.g. Schmitt et al. 1984; Rieutord & Zahn 1995; Lo & Schatzman 1997). We have preferred it because it clearly defines the boundary for a plume, which is necessary when the descending regions are also taken into account. If we want to consider descending plumes as encounteredin stellar envelopes, we just have to multiply the expressions given above for the convective variables by − 1 or take h = −1 on the axis of the cone.

Therefore, using Eq. (43), the convective fluctuations of density and temperature read (44) (45)

and for the radial mass flux (46)

Only the equation of mechanical energy conservation is used in this work instead of the radial component of the equation of motion as done in the above-mentioned works on plumes. Consequently, only the value of ρV _θ at the boundary of a plume will appear in the equations and we do not have to make an assumption on its angular variation.

We also have to introduce an expression for h(θ, ϕ) outside the ascending plumes, where matter moves downwards everywhere. To do so, we will consider two hypotheses;

• First, we assume that there are N_D descending regions that have the same geometry as the plumes. Therefore,

(47)

where the coefficient A_D will be deduced from the mass conservation.

• Second, we assume that each plume is surrounded by a hollow cone with an opening between θ_M and θ_M,D, where the flow is directed downwards. In this case, one has

(48)

and N_D = N.

We will have to assume that the sum of the fraction of a spherical surface occupied by plumes and by descending gas equals to unity (i.e., S_M + S_D = 1). Obviously, this hypothesis is an approximation as it is impossible to cover completely a spherical surface with spherical caps. Also S_M = N(1 − cosθ_M)∕2, while the expression for S_D will be different for the two hypotheses considered (see Appendix A).

3.2 Explicit form of the equations averaged over a spherical surface

In the following, we will consider two cases: that one can use either Eq. (47) or Eq. (48), depending on the treatment of the downward flow. As already discussed in Sect. 3.1, there is no net mass flux () so that the average of h over a spherical surface vanishes. Thus, using the normalization condition, Eq. (27), together with Eq. (43) and Eq. (47), allows us to write (49)

where we have defined (50)

Alternatively, if we use Eq. (48) instead of Eq. (47), one gets (51)

where (52)

Therefore, the conservation of the mass flux, through Eq. (49) and Eq. (51), provides the coefficient A_D . The integrals J_i,M, J_i,D, and are given in Appendix A.

To proceed further, we must consider that numerical simulations show, for convective envelopes, that the physical properties of convection are very different inside and outside the plumes (e.g. Stein & Nordlund 1998). To take this into account, we must consider that the anisotropy parameter (k) may have different values inside and outside the plumes. Therefore, we will consider the parameter k in the plumes and k_D in the downward flow. Indeed, one would expect that convection is more isotropic outside the plumes and therefore that k is closer to three than inside them. This is what numerical simulations for the solar envelope show (R. Samadi, priv. comm.). They also show that convection is nowhere isotropic; this shows the advantage of using the equation of mechanical energy conservation rather than the equation of motion. It is also possible that the mixing length has different values (l and l_D ) in the up- and down-moving flows. We will thus formally consider l_D ≠ l in the following. However, as the ratio l_D∕l is unknown, a pragmatic approach would be to take l_D = l in order to avoid multiplying the number of free parameters. Also Arnett et al. (2009) and Viallet et al. (2013) showed from their numerical simulations that we can take l = l_D ≈ l_CZ.

Using the same procedure, Eq. (12), together with Eq. (37) and Eq. (41), gives for the conservation of mechanical energy (53)

where, using Eq. (47), the angular integrals G_i are given by (54)

and if Eq. (48) is used instead of Eq. (47), the integrals J_i,D have to be replaced by .

For the conservation of the total luminosity (Eq. (19)), using Eq. (37) and Eq. (31), we immediately get (55)

It is interesting to note that if (i.e. if the last term of Eq. (55) is positive) in the overshooting layers, we are guaranteed to find a solution where or V _r = 0 as shown in Appendix B. In the opposite case, i.e. in the overshooting layers, it is much more problematic except if the overshooting is very small.

Although Arnett et al. (2009); Viallet et al. (2013), and Arnett et al. (2015) have shown that the averaged kinetic energy flux cannot be null everywhere (), it is interesting to discuss this extreme case because it conforms with the LMLT and corresponds to the case S_M = S_D = 1∕2 when θ_M is small. Then, Eq. (53) and Eq. (55) become (56)

and (57)

The second term of Eq. (56) is always positive and then must be negative. Consequently, in the second equation it is impossible to have , which means that overshooting is impossible above the core boundary predicted by the LMLT. One can reach the same conclusion using the fundamental equation, Eq. (12), which does not make any assumption concerning the dissipation rate of kinetic energy: (58)

We clearly see the consequence of the hypothesis stating that the averaged kinetic energy flux cancels everywhere. Since by definition , it is necessary that to fulfil the equation of mechanical energy conservation. This means that , so does . Therefore, a non-vanishing turbulent kinetic energy flux is fundamental for the existence of overshooting. This is in opposition to a lot of works that compute overshooting and neglect the kinetic energy flux. They get results because they neglect some basic equations of the problem such as the equation of mechanical energy conservation, and also because they do not consider the downwards moving matter (see Shaviv & Salpeter 1973, and the works discussed in Renzini 1987). Also, Eq. (58) shows more generally that, in the overshooting layers, the divergence of the kinetic energy flux must be negative and in absolute value larger than the kinetic energy dissipation rate.

4 Equation integrated on a spherical section of a plume

We now consider the equations describing plumes. The computation of the integrated expressions is somewhat long but requires only fundamental treatments of the equations. One has just to keep in mind that θ_M is a function of the radius and also that most terms in the equations are linear in the convective variables with the exception of some terms in Eqs. (7), (10), (16), (22), and (23). As long as the order of each term is maintained during the algebraic manipulations, the local value of any variable can be taken as equal to its average over a spherical surface and conversely.

4.1 Continuity equation

We start by integrating the equation of mass conservation (Eq. (6)); this gives (59)

When only the plumes are taken into account, the lack of a physical boundary condition at the edge of the cone forces us to adopt a hypothetical boundary condition and Taylor’s hypothesis is used. However, this closure has been criticized by Schmitt et al. (1984) and Bonin & Rieutord (1996). The latter pointed out the need to replace Taylor’s hypothesis by a proper closure using the mean-field equations and also to take into account the dissipation term in the equation of motion. Here we will not have to use Taylor’s hypothesis as we will use the mean-field equations (see Sect. 3).

4.2 Mechanical energy conservation

For the mechanical energy conservation, we use Eq. (10), which after integration gives (60)

At the plume boundary V _r = 0. Therefore, using the hypothesis that the two horizontal components of the convective velocity are equal, we have (61)

It is also worthwhile noticing that, in a plume, the radial component of the kinetic energy flux is positive while it is negative outside. Moreover, the difference between the values of ∇⋅F_K in these regions is due to either the difference of Archimede’s force, or to the difference in the dissipation rate of kinetic energy (through the amplitude term A_D ), or to the advection of kinetic energy at the plume boundary.

4.3 Thermal energy conservation

In this section, we aim at deriving the equation governing the temperature gradient. To this end, we start by integrating Eq. (23) over a horizontal section of a plume. We then obtain (62)

To go further, we have to express the perturbation of the radiative flux. To do so, we place ourselves in the limit of the diffusion and we get for the horizontal component of the radiative flux (63)

and for the perturbation of the radial component (64)

where we defined (65)

Finally, using Eq. (62) together with Eqs. (24), (63), and (64), we get 66

where, for the entropy gradient, we used the relation (67)

If we consider Eq. (66) in the overshooting zone above a convective core, because of the low radiative conductibility in stellar interior, the second and third terms of the right hand side will be the leading ones and we have with a very good approximation (this also holds deeper in the convective core) (68)

In the overshooting layers and the first term is negative while the second one is of course positive. At the bottom of the overshooting region, we expect a slightly sub-adiabatic temperature gradient and therefore will become increasingly negative as the convective flux must then become increasingly negative there. On the other hand, the radial velocity V _r,0 can only decrease. Therefore, after remaining for a while close to the adiabatic value, the temperature gradient will become more and more sub-adiabatic. This will slow down the decrease of making it possible to have the temperature gradient varying progressively from the adiabatic value to the radiative one in a significant fraction of the overshooting zone. As pointed out by Zahn (1991), the neglected terms in Eq. (68) will be significant only in a thin layer on top of the overshooting zone.

It is also interesting to compare Eq. (66) with its counterpart in the LMLT framework. There, ϵ = 0 and α₂ = 0. The second term of the right hand side member is also present and the other ones correspond to − ∇⋅ ΔF_R , which is approximated by (69)

To recover Eq. (42) in Henyey et al. (1965), we must take . Therefore, we obtain (70)

with (71)

The latter equation is identical to the usual one (except for our definition of the mean free path). The local solution of the set of equations obtained in this section simplified with the approximations F_K = 0, θ_M constant, noadvection, and also the approximate Eq. (69), allows us to recover the Böhm-Vitense cubic equation as obtained from the mean-field equations by Gabriel et al. (1974) or from the Lorenz solution or also the vortex model as shown by Arnett & Meakin (2011) and Smith & Arnett (2014).

4.4 Advection velocity

Finally, to close the system of equations, we need an additional equation to obtain the advection velocity or equivalently the radial mass flux . To this end, we start by considering Eq. (22), and simply noting that (72)

then, using Eq. (55) and performing integration gives (73)

where, usingEq. (16), we used the notation .

5 System to be integrated

In this section, we summarize the equations to be integrated. Obviously, they have to be complemented by an equation of state, the equations for the kinetic of nuclear reactions, and the three usual equations, namely (74) (75) (76)

In radiative layers, we also have to add (77)

In the convective core we have first to take into account the equations integrated over a spherical surface. More precisely, they are;

• the equation of continuity (Eqs. (49) or (51))

(78)

• the equation of mechanical energy conservation (Eq. (53)) (79)

• the equation of thermal energy conservation (Eq. (55))

(80)

with Eq. (54) for the definition of the coefficients G_i.

The same equations integrated over the “horizontal” section of a plume, together with the one giving the advection velocity, are also to be taken into account. Namely;

• the equation of continuity (Eq. (59))

(81)

• the equation of mechanical energy conservation (Eq. (82)) (82)

  where

(83)

• the super adiabatic gradient (Eq. (66)) (84)

  where

(85)

• and finally the advection velocity (Eq. (73)) (86)

It is worth noticing that Eq. (78) and Eq. (80) are algebraic and not differential.

Consequently, we have a total of ten unknowns (i.e. ()) and A_D (we may exclude Eq. (78) of the number of equations and then not consider A_D as a variable). Also, we could have added the extra variable θ_M,D given either by N(1 − cosθ_M) + N_D(1 − cosθ_MD) = 2 or N(1 − cosθ_M,D) = 2 (see Eq. (A.19)). Also, in the first case we can consider that N_D ≠N and then N_D will be a free parameter. However, we have as in the LMLT several parameters, namely the mixing length l_MLT and other more or less arbitrary constants. Therefore, we have as many equations as unknowns and we do not need the usual Taylor’s closure hypothesis. Finally, those equations have to be complemented by the boundary conditions at the centre.

6 Remarks concerning α and α₂

Before discussing the boundary conditions, it is worth considering the treatment of α and α₂ . While most often considered as constant parameters and more precisely equal in up- and down-flows, we will show that this is more a convention rather than a physically-grounded assumption, and we thus propose an alternative way of considering them.

First, let us consider Eq. (22), (87)

The average over a spherical surface of Eq. (87) being null, α and α₂ must have different values in up- and down-flows. This implies (88)

where, in the latter equation, α, α₂ and α_D, α_2,D are associated with the up- and down-flows, respectively. We draw attention to the fact that we need the mean temperature gradients obtained from the up and down flows to be the same.

Except in very small regions close to each boundary of the convective core but much thinner than one layer there, because of the very low thermal conductivity deep inside stars and also because , Eq. (84 ) reduces with a very good accuracy to (89)

Therefore we must have for the down-flow (90)

thus (91)

Equations (88) and (91) form a system connecting the two sets of variable (i.e. α, α₂ and α_D , α_2,D). Obviously, if one decides that α and α₂ are constants then the other pair must be variable. However, there is no reason to treat α, α₂ as constants and α_D, α_2,D as variables (or the other way around). The only way to avoid this is if these variables are the solution of the following two systems: (92) (93)

and (94) (95)

where C₁(r) and C₂(r) are the new arbitrary coefficients, constant over a spherical surface.

The coefficients C₁(r) and C₂(r) still depend on the radius even within the LMLT framework. This can be fixed by introducing to other coefficients, namely (96) (97)

With this choice, within the LMLT, and no longer depend on the radius and are constants. Indeed, within the LMLT, one has (98)

since the mixing length is proportional to the pressure scale height. Moreover, one has A_D = −1, NJ_1,M = N_DJ_1,D, J_4,M = N_DJ_4,D, J_1,M∕J_0,M = J_1,D∕J_0,D, J_4,M∕J_0,M = J_4,D∕J_0,D, the latter being constants. We see that the two sets of equations are identical and the two sets of unknowns are equal. However, if α and α₂ are constants, then α_D and α_2D will vary only slightly as ∇_a . This justifies once more the absence of the term in α₂ in the LMLT. There is a slight difference with the usual LMLT, which results from a different choice of the “horizontal” variation of the variables. In the LMLT, the function h( θ, ϕ) is assumed to be constant in the up- (h(θ, ϕ) = 1) and down-flows (h( θ, ϕ) = −1). This finally suggests that we should take and as true constants independent of the distance to the centre.

7 Boundary conditions at the centre

Boundary conditions at the centre of convective cores constitute a non-trivial issue even if it is a coordinate singularity rather than a physical one. Regarding 3D hydrodynamical simulations, it is difficult to obtain the temperature gradient since the 3D hydrodynamical simulations are not yet able to reach the thermal relaxation as its characteristic time is very long. Also, for the vertical component of the velocity, two recent numerical simulations by Gilet (2012) and Cai et al. (2011) adopt different boundary condition. The first one considers V _r (0) = 0 while the second one considers a non-vanishing value.

In this section, we provide boundary conditions for our model but we would like to stress that the discussion that follows is closely linked to our hypothesis that the plumes theory we decided to develop is valid down to the centre. One of the consequences of this hypothesis is that the equations have a singularity there and that the behaviour of the variables of the problem close to it is given by the requirement that the solution is regular. Of course our hypothesis can be put in question as the physics of the problem does not necessarily require a singularity at the centre. One could for instance prefer the kind of behaviour adopted by Cai et al. (2011). Nevertheless we may reasonably hope that this choice will not significantly modify the solution in the outer part of the convective core where the temperature gradient departs significantly from the adiabatic value. The same point of view justifies our decision to ignore the very small super adiabatic central core we will find in the discussion here below and to assume that the central regions of the core are adiabatic.

7.1 The temperature gradient at the centre

Let us now consider Eq. (66) (or Eq. (84)). It contains several terms that may lead to a singularity at the centre. The second to last of the right hand side member is the most singular one, so let us first compare it to the fourth one as these terms are the two components of the divergence of the radiative flux convective fluctuation. Clearly, close to the centre the non-radial component will dominate, because the “horizontal” extent of a plume becomes very small, so that the “horizontal” temperature gradient is very large. This term is still larger than the radial one far from the centre. Indeed, we can roughly write the radial term such that (99)

Therefore, the radial term is the leading one if (100)

Consequently, except for large values of θ_M approachingone, this condition will not be fulfilled in most of the radiative cores so that we have to take the non-radial component into account. However, because of the low conductivity of the gas deep in stars, the second term of the right hand side member in Eq. (66) (or Eq. (84)) is expected to become quickly the leading one.

Finally, we conclude that, because the radiative dissipation through the “horizontal” component will be important in the vicinity of the centre, our formalism leads to a small non-adiabatic core near the centre. At first, we will estimate the extent of the non-adiabatic convective core.

Close tothe centre, Eq. (66) reduces to (101)

which must be solved together with Eq. (79) and Eq. (80) where (103) (103)

where is the radiative gradient.

Thus, using Eqs. (102) and (103) into Eq. (101), Eq. (80) becomes (104)

To go further, it is necessary to take the Taylor expansion of V _r,0 and at least to the second order to ensure that the temperature gradient does not remain equal to the radiative one and becomes adiabatic at some point off-centre. To that end, we assume that (105) (106)

with n − m = 3, as demanded by Eq. (101).

Further assuming that G₁ and G₃ are constant (which will be found justified later on) and that l = β₁r (which is justified by the large radiative dissipations), we get the relations (107) (108) (109)

and (110)

together with the values m = 5, n = 8, and p = q = 6.

We have computed these quantities for a few main-sequence models with an initial composition X = 0.72 and Z = 0.015 either close to the ZAMS or close to the middle main sequence (MMS). The values of C_T and C_V depend on the assumed values for the G_i (see Table 1). They allow nevertheless a good accuracy on r₁ because it appears at the power 12 in Eq. (110). The results are given in Table 1.

Though it is not very accurate to compute the extent of the super-adiabatic convective core with just a second-order expansion, we may expect that r₁ gives its right order of magnitude. Since r₁ ≪ r(2) (r(2) being the radius of the first off-centre point of the models) in all cases considered in Table 1 and since (T(1) − T(2) being the temperature difference between the centre and the second mesh point of the models) is small, it will in practice make no difference (it will be smaller than the accuracy of the models computation) whether the temperature gradient is adiabatic down to the centre or whether there is a small super-adiabatic core. Moreover, it would be necessary to introduce many mesh points in that small domain in order to compute accurately the structure of this non-adiabatic core and particularly to follow the variation of the temperature gradient. Therefore we will ignore the latter and we will assume that until we are far enough from the centre. Then, Eq. (84) may be used to compute the departure from the adiabatic gradient. Nevertheless, this small super-adiabatic core as predicted by the plume theory as well as the nuclear reactions might help the formation ofup-rising plumes.

7.2 The boundary conditions

We will nowconsider the Taylor expansion of the equations Eq. (79) to Eq. (86), assuming that the temperature gradient is adiabatic at the centre justified in the previous section.

Let us first recast Eq. (79) and Eq. (80) as (111)

with (112)

The second term in the bracket is negligible close to the centre and we will take again l = β₁r. More generally, we suggest taking l = min[β₁r, β₂R_CC, β₃H_p], where R_CC is the radius of the convective core and the β_i are undetermined constants. This would be the usual way to define l but, following the recent findings by Arnett et al. (2009, 2015), one can instead suggest to consider l = min[β₁r;β₂R_CC]. The solution of Eq. (111) is given by (113)

which showsthat . As a result, the kinematic convective energy flux is negligible compared to the convective flux and is given by (114)

and it is different from zero at the centre. Not surprisingly, we recover the same results as Lo & Schatzman (1997).

To go further, we consider a Taylor expansion of V _r,0 and and we will keep the notations as in the previous section, (115) (116)

but now m=1 and n=0 as imposed by Eq. (113) and Eq. (114).

We now consider the equations for a plume. Let us first consider Eq. (59). The second term of the equation is negligible as it is at least proportional to r³ so that we obtain (117)

or (118)

We see that close to the centre the advection is directed towards the interior of the plume.

Let us now turn to Eq. (82). It leads to (119)

and to avoid the singularity of d lnθ_M∕dr, we must have (120)

This equation shows that θ_M is constant over the domain of the Taylor expansion. Also through the G_i appearing in the definition of C_T and C_V, the three variables N, N_D, and θ_MD appear in this equation. There is a relation between N, N_D, θ_M, and θ_MD as we have either N(1 − cosθ_M) + N_D(1 − cosθ_MD) = 2 or N(1 − cosθ_MD) = 2 (see Eq. (A.19)).

Finally,we consider Eq. (86). The leading terms of this equation include F_C,r,0 , which gives (121)

Proceeding as previously, one has (122)

This equation provides the second relation between θ_M and N. We see that the values of N and θ_M depend on two poorly known parameters α and β₁.

It could also be tempting to cancel the terms in so that (123)

Thus, (124)

This equation would enable the computation of α₂. However, we have to remember that the solutions for V _r,0 and therefore for are valid to the first order only and that if a solution taking higher-order terms had been computed, the lower order terms in would be mixed with higher order terms in F_C,r,0 so that Eq. (124 ) is not valid and α₂ remains a free parameter.

8 Case of an expanding core with a discontinuity of chemical composition at its surface

During some fraction of the main-sequence phase in stars slightly more massive than the Sun and in many stars during central helium burning, the expanding core creates a discontinuity of chemical composition at its surface. In that situation, rising convective matter will reach the discontinuity with a finite radial velocity and will try to overshoot above it, but the difference between the densities on both sides of the discontinuity will soon become very large compared to the convective density fluctuations. That is to say, if ρ₊ and ρ₋ are respectively the densities above and below the discontinuity, one has | ρ₊ − ρ₋| ≫|Δρ_0,−|. Therefore, ifwe consider the radial motion of matter above the discontinuity, it is given with a good approximation by (125)

as it is essentially slowed down by the buoyancy force, which is much larger that the dissipation term. This gives a penetration length, Δr, given by (126)

where V _r,0,− and c₋ are respectively the radial velocity and the sound speed just below the discontinuity. As using orders of magnitude given by the LMLT and (ρ₊ − ρ₋)∕ρ₋≃ 10⁻¹ − 10⁻², Δr∕H_p−≃ 10⁻⁴ and the overshooting distance will be very small. As a matter of fact, rather than to overshoot, the rising material will crush against the discontinuity as the ratio (ρ₊ − ρ₋)∕Δρ_0,− is enormous. This will generate internal gravity waves (see Meakin & Arnett 2007b,a; Arnett et al. 2009) propagating along the discontinuity and producing some mixing, which will be proportional to (X₊ − X₋) for main-sequence stars or to (Y ₊ − Y ₋) during central helium burning and to where S_M is the fraction of the spherical surface coinciding with the discontinuity occupied by rising material and is the average radial component of the impulsion of rising material at the position of the discontinuity of chemical composition.

The equations obtained in the preceding sections are still valid below this complex mixing layer coinciding just with the discontinuity, and therefore provide the value of Δρ_0,− and V _r,0,−. However, one must also determine the evolution of the mass of the convective core. To do so, we will consider stars during the main-sequence phase, but the transposition to the central helium burning phase is straightforward.

The variation of the mass of the convective core can be written as (127)

where F_X,r is the radial component of the hydrogen flux, which is given by (128)

where the subscript D corresponds to down-moving matter and X is the abundance in descending material after having interacted and exchanged some matter with the material above the discontinuity and it reads (129)

where ΔM∕M is the constant of proportionality. We recall that and stand for averaged values and therefore we have (130)

Finally, using Eq. (127) together with Eq. (128) and Eq. (130), one gets (131)

From Eq. (131), we can obtain the value of ΔM∕M, which allowsa growth rate of the core mass given by the location of the chemical composition discontinuity such that, at its lower side, that is to say in the convective core, L_R = L is verified. The value obtained is the value assumed implicitly in stellar evolution. It is also a minimal one (even though one could also consider evolutions with undershooting). This is what we consider to be an evolution without overshooting. If a larger value is chosen, it will lead to an evolution with overshooting and thus models with larger convective cores and also L_R > L at the location of the chemical discontinuity. However, decreases with increasing overshooting and this could quickly require unreasonable values of ΔM∕M to have the convective core mass still growing fast enough. This leads naturally to a maximum overshooting. We see that the mechanism that determines the extent of the overshooting is different in shrinking and in expanding convective cores. In shrinking cores, the extent of the overshooting layers is given by the rate of slowing down of the rising material computed in the frame of a coherent theory. In an expanding core, it is given by the ability of the rising material to erode the chemical discontinuity; the rate of this process fixes the extent of the overshooting region and there is no argument to choose it as proportional to the pressure scale height.

We notice from Eq. (55) that if the kinetic energy flux is positive (G₃ > 0) at the location of the chemical composition discontinuity where L ≤ L_R, the convective flux must be negative and there. If we take Spruit (2015)’s suggestion that the density fluctuation of the descending material must be positive, we can get a maximum value of ΔM∕M, assuming that since, if it was negative, the material would not move down spontaneously.

Since , if we expand this equation to the first order taking into account that the chemical composition of the down-moving material is different from that of the convective core, assuming that the rising material has the same composition as the average medium, we get (132)

and (133)

For a perfect gas, (134)

This requires a hypothesis concerning and it seems reasonable to assume that (135)

If we notice that X and X₋ in Eq. (129) correspond to X_D and X_M in Eq. (134) and since (X₊ − X₋) is given by the model, these two equations provide us with a value of ΔM∕M that according to Spruit (2015)’s suggestion is the maximum one. Also we see that the hypothesis made in Eq. (135) maximizes the value of ΔM∕M.

When the convective core grows, a parameter-free alternative approach is to use the entrainment velocity, obtained from numerical simulations, given by Meakin & Arnett (2007b). It provides the growth rate of the core boundary ν. It is given by with lnA = 0.027 ± 0.38 and n = 1.05 ± 0.21. We will define the Richardson number as in Cristini et al (2017) by (136)

where V _M is the mean velocity at the location of the discontinuity, Δr is the penetration depth in the stable region, and l is a length scale for the turbulent motions, which is often taken to be the horizontal integral scale of the turbulence near the interface, according Meakin & Arnett (2007b), so that l ≈ r_D sinθ_M(r_D). The lower boundary of the integral is not very important as the main contribution will come from the sub-adiabatic layers below the chemical discontinuity. When computing the integral, we should not forget that N² is a delta function at the location of the chemical composition discontinuity and then varies slowly in the stable region where N² > 0. As the penetration depth will be very small compared to the density scale height, we have (137)

where g is the gravitational acceleration and ρ₋ and ρ₊ are respectively the density below and above the chemical composition discontinuity so that (138)

If when the discontinuity coincides with the point where L = L_R (as assumedin stellar evolution without overshooting) the entrainment rate is larger than required by stellar evolution (Ri_B too small), the sub-adiabatic layer with ∇ < ∇_a (N²> 0) will grow until becomes sufficiently positive and V _M sufficiently small to have the required value of the entrainment rate (this will correspond to a Richardson number only slightly smaller than 1∕4). If without overshooting the entrainment rate is too small (Ri_B too large), then we have an evolution with undershooting and the discontinuity should be moved inwards in order to reduce the size of the sub-adiabatic region.

9 Conclusion

Up to now, there have been mainly three ways of handling the overshooting region at the upper boundary of convective zones. First, in most of the current existing codes, its extent is considered as a free parameter (generally a fraction of the local pressure scale height) and more importantly the temperature gradient is arbitrarily taken as either equal to the adiabatic one or to the radiative one. Second, higher-order methods based on Reynolds stress models naturally provide an estimate of the overshooting. However, such models are subject to many uncertainties and are quite difficult to implement in stellar evolutionary codes. Third, based on numerical simulations (see Meakin & Arnett 2007b; Arnett et al. 2009; Arnett & Meakin 2011; Viallet et al. 2013); Arnett et al. (2015; Sect. 3.8) propose an alternative model to account for the braking at the upper edge of convective cores and more generally to compute the structure of convective zones. They named it the 321D algorithm.

In this work, we have developed a theory of convection based on a plume model, which takes also into account the counter flow. We have also introduced a new equation that closes the system. The theory is comparable with a non-local mixing length theory but with a somewhat different interpretation of the mixing length. Its main advantage is that it allows the computation of the departure of the temperature gradient from the adiabatic value in the whole convective core, that is to say also in the overshooting layers. This finally provides eight differential equations (including the usual ones: the mass conservation, the hydrostatic equilibrium, and the thermal energy conservation) and essentially two algebraic relations for ten unknowns. It is worth noticing that the set of equations is complete if we assume that the number of regions with down- and up-moving flows is equal. If not, the number of down-moving regions is a free parameter. Interestingly enough, we do not have to make several hypotheses done in the classical plume theory, such as Taylor’s advection (or entrainment) hypothesis, which is replaced by a closure based on mean-field equations. Of course our theory still contains a few free parameters, some of which are the same as in the LMLT while others are related to the assumed geometry. In this respect, 3D hydrodynamical simulations are certainly valuable tools to help constrain those parameters. We also notice that the 321D approach as proposed by Arnett et al. (2015) is certainly more general since it does not require flow assumptions. Comparisons between our approach and the 321D model are thus desirable in the near future.

The equation giving the departure from the adiabatic gradient shows that at the bottom of the overshooting zone the temperaturegradient will be close to the adiabatic value, but that when the distance to the bottom becomes a significant fraction of the pressure scale height the temperature gradient becomes strongly sub-adiabatic and progressively reaches the radiative value. It is what is also found by numerical simulations (Brummell et al. 2002), with Xiong’s code (see Xiong & Deng 2001; Zhang & Li 2012, and the reference therein), or by use of semi-analytical models (Rempel 2004). Also the development of the theory allowed us to show rigorously, that is without any approximation, that it is necessary that the convective kinetic energy flux is not null everywhere in order to have overshooting at any boundary of a convective zone. This implies that all theories developed so far that neglect that term get results because they neglect at least one of the basic equations ofthe problem, very often the convective kinetic energy conservation. Consequently, their results must be considered with caution. Indeed, as discussed in Appendix B and shown in Arnett et al. (2015), the braking region must always be sub-adiabatic in order to have enough braking and this region does not exist in the LMLT.

Finally, we emphasize that our theory is not valid in the surface layers of stars where convective cells form (i.e. at the photosphere). However, it could be adapted to convective envelopes provided numerical simulations can provide the boundary conditions below the strongly non-adiabatic layers. We also point out that the bulk of our work focused on the case of shrinking convective cores. For an expanding core, the problem is much more difficult because one has to determine the mixing rate at the location of the chemical composition discontinuity. We did not solve that problem but we introduced an extra parameter the minimum value of which can be obtained from an evolution without overshooting. This will at least provide the order of magnitude of that parameter. This is certainly an issue that deserves more theoretical development. However, an alternative approach based on results of numerical simulations and which is parameter-free has been given by Meakin & Arnett (2007b).

Acknowledgement

We thank Arlette Noels for providing us with the models required to compute Table 1.

Appendix A: Computation of the coefficients J_i,M, J_i,D , and

In this appendix, we provide the coefficients J_i,M, J_i,D ,and that appear from angular integration on a spherical surface or over a plume horizontal section. The computation of those coefficients is simple even if somewhat lengthy. It requires only the use of simple fundamental trigonometric functions.

A.1:

The coefficients J_i,M, J_i,D

In the first hypothesis discussed at the end of Sect. 3.1, we have J_i,M = J_i,D. We thus start by defining the parameter ζ = π∕(2θ_M), so that one obtains

• for J_0,M

(A.1)

and, in the limit of small angle θ_M, one has (A.2)

• for J_1,M (A.3)

  and, in the limit of small angle θ_M, one has

(A.4)

• for J_2,M (A.5)

  and, in the limit of small angle θ_M, one has

(A.6)

• for J_3,M (A.7)

  with (A.8) (A.9)

  and, in the limit of small angle θ_M, one has

(A.10)

• for J_4,M (A.11)

  and, in the limit of small angle θ_M, one has

(A.12)

• for J_5,M (A.13)

  and, in the limit of small angle θ_M, one has

(A.14)

• for J_6,M (A.15)

  and, in the limit of small angle θ_M, one has

(A.16)

We also have (A.17)

When we assume that the geometry of the descending regions is the same as that of the plume, the coefficients J_i,D are given by the same formulae as the J_i,M (and θ_M becomes θ_D). However, when we assume that each plume is surrounded by a hollow cone with an opening between θ_M and θ_M,D where the gas is flowing downwards, the J_i,D are given by the different formulae that we write .

A.2:

The coefficients

In the following, we will consider only the coefficients that are required in the integration over a spherical surface.

First, we have for S_D (A.18)

and the condition S_M + S_D = 1 leads to (A.19)

The condition implies (see Eq. (51)) (A.20)

with (A.21)

Therefore, Eq. (A.20) gives A_D. If θ_M,D and θ_M are small, one has (A.22)

Consequently, using Eq. (A.2), one gets (A.23)

Now we consider the angular integrals G_i in Eq. (54 ). This gives

• for G₁ (A.24)

  where (A.25)

  which, for small values of θ_M,D and θ_M, reduces to

(A.26)

• for G₃ and G₄ (A.27) (A.28)

  where (A.29)

  which, for small values of θ_M,D and θ_M, reduces to

(A.30)

Appendix B: Influence of the sign of L_K on the extent of the overshooting

From Eq. (79) and Eq. (80) (B.1)

The latter equation can be recast such that (B.2)

where (B.3)

When G₃ > 0 (F_K,r > 0), then y > 0 and Eq. (B.2) shows that dy∕dr < 0 when L < L_R, that is to say above the point where L = L_R, as G₄ > 0. This implies that braking starts when the convective luminosity is negative, that is to say when the temperature fluctuation of rising material is negative, which is only possible if the temperature gradient is sub-adiabatic. This is also what is found in numerical simulations (see Arnett et al. 2015). This seems to exclude the second possibility discussed next. Then we are sure to obtain a solution y = 0. On the other hand, if G₃ < 0 (F_K,r < 0) then y < 0 and it is necessary that dy∕dr > 0 above the point where L = L_R to have a point where y = 0. However, the only term which contributes to make the derivative positive is ∇_a y∕H_p and as long as (L − L_R)∕4π − y < 0, the derivative is negative, y decreases as the radius grows, and it is impossible to have y = 0. Consequently, when G₃ < 0 the existence of a solution seems problematic or, if it exists, it will be found for a very small overshooting.

References

All Tables