© 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 $$y=\overline{y}+\text{UpDelta}y,$$(1)

where $\overline{y}$ 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 $\overline{\rho V_i}=0$. 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 $${(\frac{\partial \ln p}{\partial \ln T})}_{\rho }\frac{\text{UpDelta}T}{\overline{T}}+{(\frac{\partial \ln p}{\partial \ln \rho })}_T\frac{\text{UpDelta}\rho }{\overline{\rho }}=0,$$(2)

or equivalently $$\frac{\text{UpDelta}\rho }{\overline{\rho }}+Q\frac{\text{UpDelta}T}{\overline{T}}=0\text{with}Q=-{(\frac{\partial \ln \rho }{\partial \ln T})}_p.$$(3)

Rigorously, instead of Eq. (2) one should write $${(\frac{\partial \ln p}{\partial \ln T})}_{\rho }\frac{\text{UpDelta}T}{\overline{T}}+{(\frac{\partial \ln p}{\partial \ln \rho })}_T\frac{\text{UpDelta}\rho }{\overline{\rho }}+{\displaystyle \sum _i(}\frac{\partial \ln p}{\partial \ln X_i}{)}_{\rho ,T}\frac{\text{UpDelta}{X}_i}{\overline{X_i}}=0,$$(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) $${\nabla }_a=-\frac{\overline{p}}{C_p\overline{\rho }\overline{T}}{(\frac{\partial \ln \rho }{\partial \ln T})}_p=\frac{Q\overline{p}}{C_p\overline{\rho }\overline{T}},$$(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 $$\frac{1}{r^2}\frac{\partial }{\partial r}(r^2\rho V_r)+\frac{1}{r\sin \theta }\frac{\partial }{\partial \theta }(\sin \theta \rho V_{\theta })=0,$$(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 $\overline{\rho V}=0$. Therefore, one has $$\overline{\rho }\overline{V}+\overline{\text{UpDelta}\rho V}=0.$$(7)

As $\overline{\rho V}=0$, 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 $$\nabla \cdot (\rho VV)=\rho V\cdot (\nabla \cdot V)=\frac{\text{UpDelta}\rho }{\rho }\nabla \overline{p}-\alpha \frac{\rho V\left|V_r\right|}{l_{\text{MLT}}}.$$(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 $$\frac{\text{d}\overline{p}}{\text{d}r}=-\overline{\rho }g,$$(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 $$\nabla \cdot (\frac{1}{2}\rho V^{2}V)=\frac{\text{UpDelta}\rho }{\rho }V\cdot \nabla \overline{p}-\frac{\rho V^2\left|V_r\right|}{l},$$(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 $$l=l_d\frac{\left|V_r\right|}{\left|V\right|}.$$(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 $$\nabla \cdot \overline{F_K}=\overline{\frac{\text{UpDelta}\rho }{\rho }V}\cdot \nabla \overline{p}-\overline{\rho {ϵ}_2}.$$(12)

2.4 Thermal energy conservation

Let us now consider the thermal energy equation $$\nabla \cdot (F_H+F_R)-V\cdot \nabla p=\rho ϵ+\frac{\rho V^2\left|V_r\right|}{l},$$(13)

where F_R is the radiative flux, ϵ is the energy generation rate, and F_H is the enthalpy flux defined as $$F_H\equiv \rho UV+pV=(\overline{U}+\frac{\overline{p}}{\overline{\rho }}+\overline{T}\text{UpDelta}S)\rho V$$(14)

with S being the entropy. Using Eq. (6), the divergence of the enthalpy flux can be written as $$\nabla \cdot F_H=\rho V\cdot (\frac{\nabla \overline{p}}{\overline{\rho }}+\overline{T}\nabla \overline{S})+\nabla \cdot F_C,$$(15)

where the convective flux, F_C, is defined by $$F_C\equiv \overline{T}\text{UpDelta}S\rho V.$$(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 $$\nabla \cdot (F_K+F_C+F_R)+\rho V\cdot (\overline{T}\nabla \overline{S})=\rho ϵ.$$(17)

Averaging over a spherical surface, we recover the usual equation of thermal energy conservation $$\nabla \cdot (\overline{F_K}+\overline{F_C}+\overline{F_R})=\overline{\rho ϵ}$$(18)

or equivalently $$\overline{L_K}+\overline{L_C}+\overline{L_R}=\overline{L,}$$(19)

where ${\overline{L}}_K,{\overline{L}}_C,{\overline{L}}_R$ are the average of the kinetic, convective, radiative luminosities, and $\overline{L}$ 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 $$\nabla \cdot (\text{UpDelta}{F}_C+\text{UpDelta}{F}_K+\text{UpDelta}{F}_R)+\rho V\cdot (\overline{T}\nabla \overline{S})=\text{UpDelta}(\rho ϵ).$$(20)

In the framework of the LMLT, the first term of the latter equation reads $$\nabla \cdot (\text{UpDelta}{F}_C)=\frac{\rho T\text{UpDelta}S\left|V_r\right|}{l_{\text{MLT}}},$$(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 $$\nabla \cdot (\text{UpDelta}{F}_K+\text{UpDelta}{F}_C)=\frac{\rho T\text{UpDelta}S\left|V_r\right|}{\alpha l}+{\alpha }_2\frac{k}{2}\frac{\rho V_r^3}{l}.$$(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 $$\nabla \cdot (\text{UpDelta}{F}_R)+\rho V\cdot (\overline{T}\nabla \overline{S})=\text{UpDelta}(\rho ϵ)-\frac{\rho T\text{UpDelta}S\left|V_r\right|}{\alpha l}-{\alpha }_2\frac{k}{2}\frac{\rho V_r^3}{l},$$(23)

with $$\text{UpDelta}(\rho ϵ)=\overline{\rho ϵ}[\nu -(Q+1)\mu ]\frac{\text{UpDelta}T}{\overline{T}},$$(24)

where we used the definitions $\nu \equiv {(\frac{\partial \ln ϵ}{\partial \ln T})}_{\rho }$ and $\mu \equiv {(\frac{\partial \ln ϵ}{\partial \ln \rho })}_T$.

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: $$\begin{array}{l}\rho (r,\theta ,\varphi )=\overline{\rho (r)}+\text{UpDelta}\rho (r,\theta ,\varphi )=\overline{\rho (r)}+\text{UpDelta}{\rho }_0(r)h(\theta ,\varphi ),\\ T(r,\theta ,\varphi )=\overline{T(r)}+\text{UpDelta}T(r,\theta ,\varphi )=\overline{T(r)}+\text{UpDelta}{T}_0(r)h(\theta ,\varphi ),\end{array}$$(25)

where $\overline{\rho (r)}$ and $\overline{T(r)}$ are averages over a spherical surface as defined by Eq. (1).

For the mass flux, as already mentioned, its average vanishes so that $$\rho V(r,\theta ,\varphi )={(\rho V)}_0(r)h(\theta ,\varphi ),$$(26)

where ${(\rho V)}_0$ 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 $${\displaystyle {\int }_{S}h}(\theta ,\varphi )\text{d}\Omega =0,$$(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 $$V=V_0(r)h_1(\theta ,\varphi ).$$(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 $$V\simeq \frac{\rho V}{\overline{\rho }}(1-\frac{\text{UpDelta}\rho }{\overline{\rho }})=\frac{{(\rho V)}_0}{\overline{\rho }}h[1-(\frac{\text{UpDelta}{\rho }_0}{\overline{\rho }})h].$$(29)

In convective cores, the Mach number (ratio between the convective velocity and the sound speed) is very small. Therefore, $\text{UpDelta}\rho /\overline{\rho }\ll 1$, because it depends on the squared Mach number. Consequently, equating Eq. (28) and Eq. (29 ), one has $$h_1\simeq h\text{and}\overline{\rho }{V}_0\simeq {(\rho V)}_0.$$(30)

It is now possible to express the averaged convective quantities. Let us first consider the averaged convective flux. Using Eq. (26), it reads $$\overline{F_C}=C_p\overline{\text{UpDelta}T\rho V}=C_p\text{UpDelta}{T}_0{(\rho V)}_0{G}_1,$$(31)

where we have defined $$G_1=\frac{1}{S}{\displaystyle {\int }_S{h}^2}{r}^2\text{d}\Omega .$$(32)

For the convective velocity, using Eq. (7) together with Eq. (26), we have $$\overline{V}=-\frac{\text{UpDelta}{\rho }_0(r)}{\overline{\rho }}{V}_0(r)\frac{1}{S}{\displaystyle {\int }_{S}h}{h}_1{r}^2\text{d}\Omega =-\frac{\text{UpDelta}{\rho }_0(r)}{\overline{\rho }}{V}_0(r)G_2,$$(33)

with $$G_2=\frac{1}{S}{\displaystyle {\int }_{S}h}{h}_1{r}^2\text{d}\Omega \simeq \frac{1}{S}{\displaystyle {\int }_S{h}^2}{r}^2\text{d}\Omega =G_1.$$(34)

It is now worth considering the averaged kinetic energy flux and more precisely its radial component. One can write $${\overline{F}}_{K,r}=\overline{\frac{1}{2}\rho V^2{V}_r}=\frac{k}{2}\overline{\rho V_r^3}.$$(35)

To obtain Eq. (35) we have introduced a hypothesis concerning the isotropy of convection and we have written $$V^2=k{V}_r^2,$$(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 $$\begin{array}{ll}{\overline{F}}_{K,r}\hfill & =\overline{\frac{1}{2}\rho V^2{V}_r}=\frac{k}{2}\overline{\rho V_r^3}\hfill \\ \hfill & =\frac{k}{2}{(\rho V_r^3)}_0\frac{1}{S}{\displaystyle \int h}{h}_1^2{r}^2\text{d}\Omega =\frac{k}{2}{(\rho V_r^3)}_0{G}_3,\hfill \end{array}$$(37)

where we have introduced the notation $${(\rho V_r^3)}_0\equiv {(\rho V_r)}_0{V}_{0,r}^2$$(38)

and $$G_3=\frac{1}{S}{\displaystyle {\int }_{S}h}{h}_1^2{r}^2\text{d}\Omega \simeq \frac{1}{S}{\displaystyle {\int }_S{h}^3}{r}^2\text{d}\Omega .$$(39)

The averaged luminosity associated with the radial component of the averaged kinetic energy flux can then be written as ${\overline{L}}_K=4\pi r^2{\overline{F}}_{K,r}$.

Finally, one can express the kinetic energy dissipation rate into heat using the isotropic parameter (Eq. (36)), such as $$\overline{\rho {ϵ}_2}=\frac{\overline{\rho V^2\left|V_r\right|}}{l}=\frac{k\overline{\rho V_r^2\left|V_r\right|}}{l}.$$(40)

Therefore, one has $$\overline{\rho {ϵ}_2}=\frac{k{(\rho V_r^3)}_0}{l}\frac{1}{S}{\displaystyle {\int }_{S}h}{h}_1\left|h_1\right|r^2\text{d}\Omega =\frac{k{(\rho V_r^3)}_0}{l}{G}_4,$$(41)

with $$G_4=\frac{1}{S}{\displaystyle {\int }_{S}h}{h}_1\left|h_1\right|r^2\text{d}\Omega \simeq \frac{1}{S}{\displaystyle {\int }_S{h}^2}\left|h\right|r^2\text{d}\Omega .$$(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 $$h(\theta ,\varphi )\equiv \cos (\frac{\pi }{2}\frac{\theta }{{\theta }_M}).$$(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 $$\begin{array}{ll}\frac{\text{UpDelta}\rho }{\overline{\rho }}\hfill & =\frac{\text{UpDelta}{\rho }_0}{\overline{\rho }}\cos (\frac{\pi }{2}\frac{\theta }{{\theta }_M})\hfill \\ \frac{\text{UpDelta}T}{\overline{T}}\hfill & =\frac{\text{UpDelta}{T}_0}{\overline{T}}\cos (\frac{\pi }{2}\frac{\theta }{{\theta }_M}),\hfill \end{array}$$(44) (45)

and for the radial mass flux $$\rho V_r={(\rho V_r)}_0\cos (\frac{\pi }{2}\frac{\theta }{{\theta }_M}).$$(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,

$$h(\theta ,\varphi )=A_D\cos (\frac{\pi }{2}\frac{\theta }{{\theta }_{M,D}}),$$(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

$$h(\theta ,\varphi )=A_D\sin (\pi \frac{\theta -{\theta }_M}{{\theta }_{M,D}-{\theta }_M})$$(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 ($\overline{\rho V}=0$) 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 $$N{J}_{0,M}+A_D{N}_D{J}_{0,D}=0,$$(49)

where we have defined $$\begin{array}{ll}{J}_{0,M}\hfill & ={\displaystyle {\int }_0^{{\theta }_M}\cos }(\frac{\pi }{2}\frac{\theta }{{\theta }_M})\sin \theta \text{d}\theta \hfill \\ J_{0,D}\hfill & ={\displaystyle {\int }_0^{{\theta }_{M,D}}\cos }(\frac{\pi }{2}\frac{\theta }{{\theta }_{M,D}})\sin \theta \text{d}\theta .\hfill \end{array}$$(50)

Alternatively, if we use Eq. (48) instead of Eq. (47), one gets $$J_{0,M}+A_D{J}_{0,D}^{\prime }=0,$$(51)

where $$J_{0,D}^{\prime }={\displaystyle {\int }_{{\theta }_M}^{{\theta }_{M,D}}\sin }(\pi \frac{\theta -{\theta }_M}{{\theta }_{M,D}-{\theta }_M})\sin \theta \text{d}\theta .$$(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 $J_{i,D}^{\prime }$ 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 $$\begin{array}{ll}\frac{1}{r^2}\frac{\text{d}}{\text{d}r}\hfill & (\frac{k{r}^2}{2}{(\rho V_r^3)}_0{G}_3)\hfill \\ \hfill & =-(\frac{\text{UpDelta}{\rho }_0}{\overline{\rho }}){(\rho V_r)}_{0}g{G}_1-\frac{k{(\rho V_r^3)}_0}{l}{G}_4,\hfill \end{array}$$(53)

where, using Eq. (47), the angular integrals G_i are given by $$\begin{array}{ll}{G}_1\hfill & =G_2=\frac{N}{2}{J}_{1,M}+\frac{N_D{A}_D^2}{2}{J}_{1,D}\hfill \\ G_3\hfill & =\frac{N}{2}{J}_{4,M}+\frac{k_D}{k}\frac{N_D{A}_D^3}{2}{J}_{4,D}\hfill \\ G_4\hfill & =\frac{N}{2}{J}_{4,M}+\frac{l{k}_D}{l_{D}k}\frac{N_D\left|A_D^3\right|}{2}{J}_{4,D},\hfill \end{array}$$(54)

and if Eq. (48) is used instead of Eq. (47), the integrals J_i,D have to be replaced by $J_{i,D}^{\prime }$.

For the conservation of the total luminosity (Eq. (19)), using Eq. (37) and Eq. (31), we immediately get $$\frac{\overline{L}(r)-{\overline{L}}_R(r)}{4\pi r^2}=C_p\text{UpDelta}{T}_0{(\rho V_r)}_0{G}_1+\frac{k}{2}{(\rho V_r^3)}_0{G}_3.$$(55)

It is interesting to note that if ${\overline{F}}_{K,r}>0$ (i.e. if the last term of Eq. (55) is positive) in the overshooting layers, we are guaranteed to find a solution where ${\overline{F}}_{K,r}=0$ or V _r = 0 as shown in Appendix B. In the opposite case, i.e. ${\overline{F}}_{K,r}<0$ 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 (${\overline{F}}_{K,r}=0$), 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 $$\frac{\text{UpDelta}{\rho }_0}{\overline{\rho }}{(\rho V_r)}_{0}g{G}_1+\frac{k{(\rho V_r^3)}_0}{l}{G}_4=0,$$(56)

and $$\frac{\overline{L}(r)-{\overline{L}}_R(r)}{4\pi r^2}=C_p\text{UpDelta}{T}_0{(\rho V_r)}_0{G}_1.$$(57)

The second term of Eq. (56) is always positive and then $\text{UpDelta}{\rho }_0{(\rho V_r)}_0/\overline{\rho }$ must be negative. Consequently, in the second equation it is impossible to have $\overline{L}(r)-{\overline{L}}_R(r)<0$, 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: $$\nabla \cdot \overline{F_K}=\overline{\frac{\text{UpDelta}\rho }{\rho }V}\cdot \nabla \overline{p}-\overline{\rho {ϵ}_2}.$$(58)

We clearly see the consequence of the hypothesis stating that the averaged kinetic energy flux cancels everywhere. Since by definition $\overline{\rho {ϵ}_2}>0$, it is necessary that $\overline{\text{UpDelta}\rho V_r/\rho }<0$ to fulfil the equation of mechanical energy conservation. This means that $\overline{\text{UpDelta}T{V}_r/T}>0$, so does $\overline{L}(r)-{\overline{L}}_R(r)$. 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 $$\begin{array}{ll}{J}_{0,M}\frac{\text{d}}{\text{d}r}[r^2{(\rho V_r)}_0]\hfill & +J_{2,M}{r}^2{(\rho V_r)}_0\frac{\text{d}\ln {\theta }_M}{\text{d}r}\hfill \\ \hfill & +r\sin {\theta }_M{(\rho V_{\theta })}_{{\theta }_M}=0.\hfill \end{array}$$(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 $$\begin{array}{ll}{J}_{4,M}\hfill & \frac{\text{d}}{\text{d}r}[r^2\frac{k}{2}{(\rho V_r^3)}_0]+\frac{3}{2}k{r}^2{(\rho V_r^3)}_0\frac{\text{d}\ln {\theta }_M}{\text{d}r}{J}_{5,M}+r{(\sin \theta F_{K,\theta })}_{{\theta }_M}\hfill \\ \hfill & =-r^2[\frac{\text{UpDelta}{\rho }_0}{\overline{\rho }}{(\rho V_r)}_{0}g{J}_{1,M}+\frac{k{(\rho V_r^3)}_0}{l}{J}_{4,M}].\hfill \end{array}$$(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 $${(F_{K,\theta })}_{{\theta }_M}={(\rho V_{\theta }^3)}_{{\theta }_M}.$$(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 $$\text{begin{align}}\int \text{\_}0\text{^}{\theta }_M\text{nn\&}\{\frac{1}{r^2}\frac{\partial }{\partial r}[r^2\text{UpDelta}{F}_{R,r}]+\rho V_r\overline{T}\frac{\text{d}\overline{S}}{\text{d}r}-\text{UpDelta}(\rho ϵ)\}{r}^2\sin \theta \text{d}\theta \text{nn\&}+r^2\frac{C_p\overline{T}}{\alpha l}(\frac{\text{UpDelta}{T}_0}{\overline{T}}){(\rho V_r)}_0{J}_{1,M}+{\alpha }_2\frac{k}{2}\frac{{(\rho V_r^3)}_0}{l}{r}^2{J}_{4,M}\text{nn\&}+r{[\sin \theta (\text{UpDelta}{F}_{R,\theta })]}_{{\theta }_M}=0.\text{nnend{align}}$$(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 $${(\text{UpDelta}{F}_{R,\theta })}_{{\theta }_M}=\frac{4ac}{3}\overline{(\frac{T^4}{\kappa \rho })}(\frac{\text{UpDelta}{T}_0}{\overline{T}})\frac{\pi }{2}\frac{1}{r{\theta }_M},$$(63)

and for the perturbation of the radial component $$\begin{array}{ll}\text{UpDelta}{F}_{R,r}\hfill & ={\chi }_1{\overline{F}}_{R,r}\cos (\frac{\pi }{2}\frac{\theta }{{\theta }_M})+{\chi }_2\frac{\text{d}\overline{S}}{\text{d}r}[\ln (\frac{\text{UpDelta}{T}_0}{\overline{T}})]\cos (\frac{\pi }{2}\frac{\theta }{{\theta }_M})\hfill \\ \hfill & +{\chi }_2(\frac{\pi }{2}\frac{\theta }{{\theta }_M})\frac{\text{d}\ln {\theta }_M}{\text{d}r}\sin (\frac{\pi }{2}\frac{\theta }{{\theta }_M}),\hfill \end{array}$$(64)

where we defined $$\begin{array}{ll}{\chi }_1\hfill & =[(4-{\kappa }_T)+Q(1+{\kappa }_{\rho })](\frac{\text{UpDelta}{T}_0}{\overline{T}})\hfill \\ {\chi }_2\hfill & =-\frac{4ac}{3}\overline{(\frac{T^4}{\kappa \rho })}(\frac{\text{UpDelta}{T}_0}{\overline{T}}).\hfill \end{array}$$(65)

Finally, using Eq. (62) together with Eqs. (24), (63), and (64), we get $$\begin{array}{l}{C}_p\overline{T}{(}_{\rho }(\overline{\nabla }-{\nabla }_a)\frac{\text{d}\ln \overline{p}}{\text{d}r}=[\nu -(Q+1)\mu ]\overline{\rho ϵ}(\frac{\text{UpDelta}{T}_0}{\overline{T}})\\ -\frac{C_p\overline{T}}{\alpha l}(\frac{\text{UpDelta}{T}_0}{\overline{T}}){(}_{\rho }\frac{J_{1,M}}{J_{0,M}}-{\alpha }_2\frac{k}{2}\frac{{(\rho V_r^3)}_0}{l}\frac{J_{4,M}}{J_{0,M}}\\ -\frac{1}{r^2}\frac{\text{d}\overline{S}}{\text{d}r}\left\{r^2{\chi }_2\frac{\text{d}\ln (\text{UpDelta}{T}_0/\overline{T})}{\text{d}r}\right\}-\frac{1}{r^2}\frac{\text{d}\overline{S}}{\text{d}r}\left\{r^2{\chi }_2\frac{\text{d}\ln {\theta }_M}{\text{d}r}\right\}\frac{J_{2,M}}{J_{0,M}}\\ -\frac{1}{r^2}\frac{\text{d}\overline{S}}{\text{d}r}\left\{r^2{\chi }_1{\overline{F}}_{R,r}\right\}+{\chi }_1{\overline{F}}_{R,r}\frac{\text{d}\ln {\theta }_M}{\text{d}r}\frac{J_{2,M}}{J_{0,M}}\\ +{\chi }_2\{\frac{J_{2,M}+J_{3,M}}{J_{0,M}}{(}^{\frac{\text{d}\ln {\theta }_M}{\text{d}r}}-\frac{J_{2,M}}{J_{0,M}}\frac{\text{d}\ln {\theta }_M}{\text{d}r}+\frac{2}{\pi }\frac{\sin {\theta }_M}{{\theta }_M{r}^2}\frac{1}{J_{0,M}}\},\end{array}$$66

where, for the entropy gradient, we used the relation $$\frac{\text{d}\overline{S}}{\text{d}r}=C_p(\overline{\nabla }-{\nabla }_a)\frac{\text{d}\ln \overline{p}}{\text{d}r}.$$(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) $$(\overline{\nabla }-{\nabla }_a)=\frac{H_p}{\alpha l}(\frac{\text{UpDelta}{T}_0}{\overline{T}})\frac{J_{1,M}}{J_{0,M}}+{\alpha }_2\frac{k}{2}\frac{V_{r,0}^2{H}_p}{C_p\overline{T}l}\frac{J_{4,M}}{J_{0,M}}.$$(68)

In the overshooting layers $\text{UpDelta}{T}_0/\overline{T}<0$ 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 $\text{UpDelta}{T}_0/\overline{T}$ 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 $\text{UpDelta}{T}_0/\overline{T}$ 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 $$-\nabla \cdot \text{UpDelta}{F}_R\simeq \frac{4ac}{3}\overline{(\frac{T^4}{\kappa \rho })}\nabla \cdot \nabla (\frac{\text{UpDelta}{T}_0}{\overline{T}})\simeq \frac{4ac}{3{\mathcal{L}}^2}\overline{(\frac{T^4}{\kappa \rho })}.$$(69)

To recover Eq. (42) in Henyey et al. (1965), we must take $\mathcal{L}=\sqrt{2}\pi l_{\text{MLT}}$. Therefore, we obtain $$V_r(\overline{\nabla }-{\nabla }_{\text{ad}})\frac{\text{d}\ln \overline{p}}{\text{d}r}=-[\frac{1}{\alpha \tau }+\omega {\omega }_R](\frac{\text{UpDelta}{T}_0}{\overline{T}})$$(70)

with $${\omega }_R=\frac{4ac}{3}\overline{\frac{T^3}{C_p\kappa \rho {\mathcal{L}}^2}}.$$(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 ${(\rho V_r)}_0$. To this end, we start by considering Eq. (22), and simply noting that $$\text{UpDelta}{F}_{K,r}+\text{UpDelta}{F}_{C,r}=F_{K,r}+F_{C,r}-({\overline{F}}_{K,r}+{\overline{F}}_{C,r}),$$(72)

then, using Eq. (55) and performing integration gives $$\begin{array}{l}{J}_{4,M}\frac{\text{d}\overline{S}}{\text{d}r}[r^2\frac{k}{2}{(}_{\rho }]+3J_{5,M}{r}^2\frac{k}{2}{(}_{\rho }\frac{\text{d}\ln {\theta }_M}{\text{d}r}\\ +\frac{\text{d}\overline{S}}{\text{d}r}(r^2{F}_{C,r,0})J_{1,M}+r^2{F}_{C,r,0}\frac{\text{d}\ln {\theta }_M}{\text{d}r}{J}_{6,M}\\ -(1-\cos {\theta }_M)\frac{\text{d}\overline{S}}{\text{d}r}\{r^2[C_p\overline{T}(\frac{\text{UpDelta}{T}_0}{\overline{T}}){(}_{\rho }{G}_1+\frac{k}{2}{(}_{\rho }{G}_3]\}\\ +r\sin {\theta }_M{(}_{F_{K,\theta }}=r^2\frac{\rho C_{p}T}{\alpha l}(\frac{\text{UpDelta}{T}_0}{\overline{T}})J_{1,M}+{\alpha }_2\frac{k}{2}\frac{{(\rho V_r^3)}_0}{l}{r}^2{J}_{4,M},\end{array}$$(73)

where, usingEq. (16), we used the notation $F_{C,r,0}\equiv \overline{T}\text{UpDelta}{S}_0{(\rho V_r)}_0$.

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 $$\begin{array}{ll}\frac{\text{d}m}{\text{d}r}\hfill & =4\pi \rho r^2\hfill \\ \frac{\text{d}p}{\text{d}r}\hfill & =-G\frac{m}{r^2}\rho \hfill \\ \frac{\text{d}L}{\text{d}r}\hfill & =4\pi \rho r^2[{ϵ}_N-T\frac{\text{d}S}{\text{d}t}]=4\pi \rho r^2ϵ.\hfill \end{array}$$(74) (75) (76)

In radiative layers, we also have to add $$\frac{\text{d}T}{\text{d}r}=-\frac{4ac}{3}\frac{\kappa \rho }{T^3}\frac{L(r)}{4\pi r^2}.$$(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))

$$N{J}_{0,M}+A_D{N}_D{J}_{0,D}\text{or}{J}_{0,M}+A_D{J}_{0,D}^{\prime }=0,$$(78)

• the equation of mechanical energy conservation (Eq. (53)) $$\begin{array}{ll}\frac{1}{r^2}\frac{\text{d}\overline{S}}{\text{d}r}(\frac{k{r}^2}{2}{(\rho V_r^3)}_0{G}_3)\hfill & =-(\frac{\text{UpDelta}{\rho }_0}{\overline{\rho }}){(\rho V_r)}_{0}g{G}_1\hfill \\ \hfill & -\frac{k{(\rho V_r^3)}_0}{l}{G}_4,\hfill \end{array}$$(79)

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

$$\frac{\overline{L}(r)-{\overline{L}}_R(r)}{4\pi r^2}=C_p\text{UpDelta}{T}_0{(\rho V_r)}_0{G}_1+\frac{k}{2}{(\rho V_r^3)}_0{G}_3,$$(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))

$$\begin{array}{ll}{J}_{0,M}\frac{\text{d}\overline{S}}{\text{d}r}[r^2{(\rho V_r)}_0]\hfill & +J_{2,M}{r}^2{(\rho V_r)}_0\frac{\text{d}\ln {\theta }_M}{\text{d}r}\hfill \\ \hfill & +r\sin {\theta }_M{(\rho V_{\theta })}_{{\theta }_M}=0.\hfill \end{array}$$(81)

• the equation of mechanical energy conservation (Eq. (82)) $$\begin{array}{l}{J}_{4,M}\frac{\text{d}\overline{S}}{\text{d}r}[r^2\frac{k}{2}{(}_{\rho }]+\frac{3}{2}k{r}^2{(}_{\rho }\frac{\text{d}\ln {\theta }_M}{\text{d}r}{J}_{5,M}\\ +r{(}_{\sin }=-r^2[\frac{\text{UpDelta}{\rho }_0}{\overline{\rho }}{(}_{\rho }g{J}_{1,M}+\frac{k{(\rho V_r^3)}_0}{l}{J}_{4,M}],\end{array}$$(82)

  where

$${(F_{K,\theta })}_{{\theta }_M}={(\rho V_{\theta }^3)}_{{\theta }_M}.$$(83)

• the super adiabatic gradient (Eq. (66)) $$\begin{array}{l}{C}_p\overline{T}{(}_{\rho }(\overline{\nabla }-{\nabla }_a)\frac{\text{d}\ln \overline{p}}{\text{d}r}=[\nu -(Q+1)\mu ]\overline{\rho ϵ}(\frac{\text{UpDelta}{T}_0}{\overline{T}})\\ -\frac{C_p\overline{T}}{\alpha l}(\frac{\text{UpDelta}{T}_0}{\overline{T}}){(}_{\rho }\frac{J_{1,M}}{J_{0,M}}-{\alpha }_2\frac{k}{2}\frac{{(\rho V_r^3)}_0}{l}\frac{J_{4,M}}{J_{0,M}}\\ -\frac{1}{r^2}\frac{\text{d}\overline{S}}{\text{d}r}\left\{r^2{\chi }_2\frac{\text{d}\ln (\text{UpDelta}{T}_0/\overline{T})}{\text{d}r}\right\}-\frac{1}{r^2}\frac{\text{d}\overline{S}}{\text{d}r}\left\{r^2{\chi }_2\frac{\text{d}\ln {\theta }_M}{\text{d}r}\right\}\frac{J_{2,M}}{J_{0,M}}\\ -\frac{1}{r^2}\frac{\text{d}\overline{S}}{\text{d}r}\left\{r^2{\chi }_1{\overline{F}}_{R,r}\right\}+{\chi }_1{\overline{F}}_{R,r}\frac{\text{d}\ln {\theta }_M}{\text{d}r}\frac{J_{2,M}}{J_{0,M}}\\ +{\chi }_2\{\frac{J_{2,M}+J_{3,M}}{J_{0,M}}{(}^{\frac{\text{d}\ln {\theta }_M}{\text{d}r}}-\frac{J_{2,M}}{J_{0,M}}\frac{\text{d}\ln {\theta }_M}{\text{d}r}+\frac{2}{\pi }\frac{\sin {\theta }_M}{{\theta }_M{r}^2}\frac{1}{J_{0,M}}\},\end{array}$$(84)

  where

$$\frac{\text{d}\overline{S}}{\text{d}r}=C_p(\overline{\nabla }-{\nabla }_a)\frac{\text{d}\ln \overline{p}}{\text{d}r}.$$(85)

• and finally the advection velocity (Eq. (73)) $$\begin{array}{l}{J}_{4,M}\frac{\text{d}\overline{S}}{\text{d}r}[r^2\frac{k}{2}{(}_{\rho }]+3J_{5,M}{r}^2\frac{k}{2}{(}_{\rho }\frac{\text{d}\ln {\theta }_M}{\text{d}r}\\ +\frac{\text{d}\overline{S}}{\text{d}r}(r^2{F}_{C,r,0})J_{1,M}+r^2{F}_{C,r,0}\frac{\text{d}\ln {\theta }_M}{\text{d}r}{J}_{6,M}\\ -(1-\cos {\theta }_M)\frac{\text{d}\overline{S}}{\text{d}r}\{r^2[C_p\overline{T}(\frac{\text{UpDelta}{T}_0}{\overline{T}}){(}_{\rho }{G}_1+\frac{k}{2}{(}_{\rho }{G}_3]\}\\ +r\sin {\theta }_M{(}_{F_{K,\theta }}=r^2\frac{\rho C_{p}T}{\alpha l}(\frac{\text{UpDelta}{T}_0}{\overline{T}})J_{1,M}+{\alpha }_2\frac{k}{2}\frac{{(\rho V_r^3)}_0}{l}{r}^2{J}_{4,M}.\end{array}$$(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. ($m,p,T,L,{(\rho V_r)}_0,\text{UpDelta}{\rho }_0/\overline{\rho },V_{{\theta }_M},{\theta }_M,N$)) 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), $$\nabla \cdot (\text{UpDelta}{F}_K+\text{UpDelta}{F}_C)=\frac{\rho T\text{UpDelta}S\left|V_r\right|}{\alpha l}+{\alpha }_2\frac{k}{2}\frac{\rho V_r^3}{l}.$$(87)

The average over a spherical surface of Eq. (87) being null, α and α₂ must have different values in up- and down-flows. This implies $$\begin{array}{l}\frac{C_p\rho T}{\alpha }(\frac{\text{UpDelta}{T}_0}{\overline{T}})V_{r,0}N{J}_{1,M}+{\alpha }_2\frac{k}{2}{(}_{\rho }N{J}_{4,M}\\ =\frac{C_p\rho Tl}{{\alpha }_D{l}_D}(\frac{\text{UpDelta}{T}_0}{\overline{T}})A_D^2{V}_{r,0}{N}_D{J}_{1,D}+{\alpha }_{2,D}\left|A_D^3\right|\frac{k_{D}l}{2l_D}{(}_{\rho }{N}_D{J}_{4,D},\end{array}$$(88)