Convection Theory and Relevant Problems in Stellar Structure, Evolution, and Pulsational Stability I Convection Theory and Structure of Convection Zone and Stellar Evolution

A non-local and time-dependent theory of convection was briefly described. This theory was used to calculate the structure of solar convection zones, the evolution of massive stars, lithium depletion in the atmosphere of the Sun and late-type dwarfs, and stellar oscillations (in Part Ⅱ). The results show that: 1) the theoretical turbulent velocity and temperature fields in the atmosphere and the thermal structure of the convective envelope of the Sun agree with the observations and inferences from helioseismic inversion very well. 2) The so-called semi-convection contradiction in the evolutionary calculations of massive stars was removed automatically, as predicted by us. The theoretical evolution tracks of massive stars run at higher luminosity and the main sequence band becomes noticeably wider in comparison with those calculated using the local mixing-length theory (MLT). This means that the evolutionary mass for a given luminosity was overestimated and the width of the main sequence band was underestimated by the local MLT, which may be part of the reason for the contradiction between the evolutionary and pulsational masses of Cepheid variables and the contradiction between theoretical and observed distributions of luminous stars in the H-R diagram. 3) The predicted lithium depletion, in general, agrees well with the observation of the Sun and Galactic open clusters of different ages. 4) Our theoretical results for non-adiabatic oscillations are in good agreement with the observed mode instability from classic variables of high-luminosity red giants. Almost all the instability strips of the classical pulsating variables (including the Cepheid, δ Scuti, γ Doradus, βCephei, and SPB strips) were reproduced (Part Ⅱ).


INTRODUCTION
Convection occurs within most stars. Convection causes the transfer of energy and momentum, and the mixing of matter in stellar interiors. Therefore, it strongly influences the internal structure, evolution, and pulsational stability of stars. Stellar convection is always turbulent due to its large scale. It is impossible to expect to have a perfect convection theory until a better understanding of turbulence has been achieved. Up to now the most popularly applied convection theory is still the MLT developed by Bohm-Vitense (1958), and its many versions. The most obvious advantages of the MLT are its straightforwardness in physical picture and simplicity of use. However, the MLT is not a dynamic theory following the hydrodynamic equations and turbulence theory, but is a phenomenological theory based on a simple analogy of turbulence with the kinetic theory of gas molecules. In fact, turbulence is more complex than the kinetic theory of gas molecules. The fundamental shortcoming of the MLT is that it cannot give a correct description for the dynamic behaviors of turbulent convection. When dealing with the dynamic problems of turbulent convection, such as time-dependent and non-local convection, this shortcoming becomes prominent. A possible approach to turbulent convection is direct hydrodynamic numerical simulation. Taking into account the fact that both the time and space scales of turbulence are more than ten orders of magnitude less than the stellar ones, a direct hydrodynamic numerical simulation of stellar convection would be impossible for the calculation of stellar structure and evolution in the foreseeable future. Therefore, it is a reasonable choice to develop a theoretical approach exact enough and simple enough for treating turbulent convection in the calculation of stellar evolution and oscillations. Helio-and astero-seismology have made significant progress in the past two decades thanks to the GONG, SOHO, OGLE, MACHO, 2MASS, CoRoT, and Kepler ground-based and space projects. They offer the best opportunity to test stellar evolution and convection theories. In the present paper we pay close attention to convection theory and the relevant problems in stellar structure, evolution, and pulsational stability. Based on the reflections above, we decided to abandon the phenomenological MLT, and to develop a non-local and time-dependent theory of convection based on hydrodynamic equations and turbulence theory (Xiong, 1978;Xiong, 1980;Xiong, 1981;Xiong, 1989;Xiong et al., 1997;Deng et al., 2006). In contrast to the MLT, it is a dynamic theory of correlation functions of turbulent velocity and temperature (and, in addition, of the element abundance for chemically inhomogeneous stars) following the hydrodynamic equations and turbulence theory. So, it can be expected that our theory has a more solid hydrodynamic foundation and can give a more exact description of the hydrodynamic behavior of turbulent convection than the MLT does. Our purpose is to improve the treatment of overshooting mixing in calculations of stellar evolution and the treatment of dynamic and thermodynamic coupling between convection and oscillations in calculations of stellar oscillations. Canuto (1993Canuto ( , 1997Canuto ( , 1999 and Li of the Yunnan group of China (Zhang, 2012a;Zhang and Li, 2012a;Zhang, 2012b;Zhang and Li, 2012b;Zhang, 2013) proposed a similar theory. In the present paper, we do not attempt to make a comprehensive review and comparison for convection theories and relevant problems, which would be a large and difficult amount of work. In A Non-Local and Time-Dependent Theory of Convection section, a brief description of our non-local and time-dependent convection theory is given. The following sections review the progress of its applications in theoretical calculations of the structure of the solar convection zone (Structure of the Solar Convection Zone section), of stellar evolution and lithium depletion in the atmosphere of the Sun and late-type dwarfs (Overshooting Mixing and Stellar Evolution section), and of stellar oscillations (in PartⅡ). These applications are not only the primary motivations of our research on the theory of stellar convection, but also an important means for testing convection theory. A summary and discussion follow in section 7 of PartⅡ, where we try to evaluate the successes and failures of our theory, analyze their reasons and identify direction for improvement. It is known that this paper will contain strong personal bias. Our viewpoint on the excitation mechanism for red giants may not gain widespread approval. Nowadays, our understanding of many problems is still inconclusive, and will need more time to verify. Disputes over different academic viewpoints will contribute to the development of science. Our theory is far from perfect, however it has shown obvious improvements in comparison with MLT.
This article is divided into two parts. Convection theory and its applications in theoretical calculations of the structure of the solar convection zone and for stellar evolution are contained in PartⅠ; the applications of the theory in calculations of stellar oscillations, and summary and discussions are contained in PartⅡ.

A NON-LOCAL AND TIME-DEPENDENT THEORY OF CONVECTION
In order to overcome the defects of the MLT, we developed a nonlocal and time-dependent theory of convection based on the hydrodynamic equations and turbulence theory (Xiong, 1978;Xiong, 1980;Xiong, 1981;Xiong, 1989;Xiong et al., 1997;Deng et al., 2006). In our theory, the classical Reynolds decomposition is used. Each of the physical variables X is expressed as the sum of its average value X and (Eulerian) turbulent fluctuation X ′ , substituting Eq. 1 into the hydrodynamic equations and the radiation transfer equation, and making a Taylor expansion for X ′ , retaining only its first-order terms, ignoring all the second and higher-order terms, then averaging for each of the equations, we can obtain the following average hydrodynamic equations and the radiation transfer equation, where D Dt is the comoving differential, and T, P, ρ are, respectively, the temperature, pressure, and density of gas, C P is the specific heat at constant pressure, α −(zln ρ/zln T) P the expansion coefficient, ε N the effective nuclear energy-generation rate per gram, σ ik the viscosity tensor, Φ the gravitational potential, κ the radiative opacity, and F i r the radiation flux vector. It can be seen that, when convection sets in, an extra term ρu ′ i u ′ k in the average equation of momentum conservation and another extra term C p ρu ′ i T ′ in the average equation of energy conservation will emerge. They are, respectively, the well-known Reynolds stress and convective enthalpy flux. In the present paper, all the equations are expressed in the tensor format. g ik is the metric tensor for the current framework. The implicit summation rule of tensors is used, i.e., a summation should be performed for an index k running from 1 to 3 if both the subscript and superscript k appear in a term. Subtracting the average equations, Eqs. 3, 4, from the corresponding original equations of momentum and energy conservation, after a long derivation, it is not difficult to obtain the following dynamic equations of turbulent velocity and temperature: where w ′ i ρu ′ i /ρ is the density-weighted turbulent velocity, ∇ ad αP/ρC p T (Γ 2 − 1)/Γ 2 is the adiabatic temperature gradient, C P,T (zln C p /zln T) p and C P,P (zln C p /zln P) T are, respectively, the partial derivative of Cp with respect to T and P, and Γ 2 (along with Γ 1 and Γ 3 below) are the adiabatic exponents introduced by Chandrasekhar (1939). Starting from Eqs. 7, 8, it is not difficult to constitute the dynamic equations of auto-and crosscorrelations for turbulent velocity and temperature. For example, we have the equation of velocity correlations as, The terms in square brackets on the left-hand side of Eq. 9 are the energy fluxes of pressure and viscous stress. They are negligible in comparison with the next two terms when it is assumed that the size scale of turbulent elements is far smaller than the characteristic length of the mean fluid fields. This implies that the turbulence is near quasi-isotropic. The terms on the right-hand side of Eq. 9 are the turbulent dissipation due to molecular viscosity. Based on isotropic turbulence theory, this dissipation can be expressed as (Hinze, 1975) where η e is the Heisenberg eddy coupling constant, l e is the characteristic length of energy-containing eddies of turbulence, and x is the rms turbulent velocity, In the dynamic equations for second-order correlations, the third-order correlations must appear due to nonlinearity of the hydrodynamic equations. The fourth-order correlations will appear if we constitute the dynamic equation for the third-order correlations. The equations of the correlation functions are never closed, as is well-known in turbulence theory (Hinze, 1975). In order to close the equation of second-order correlations, we use the gradient-type diffusion approximation for treatment of the thirdorder correlations (Xiong, 1980;Xiong, 1981;Xiong, 1989) u′ k w′ i w′ j −w′ k w′ d τ c ∇ α w′ i w′ j , where is the lifetime of turbulence and ∧ is the diffusion length of turbulence. We then further assume that both l e and ∧ are proportional to the local pressure scale height H p , The third-order correlations represent the non-local transport of turbulent convection. ρu ′ k w ′ i w ′ j is the energy flux of turbulent stress. The gradient-type diffusion approximation means that turbulence diffuses from one location where it is strong to another location where it is weak. This is a very simple but reasonable physical assumption. Another closure scheme is to construct dynamic equations of third-order correlations (Canuto, 1993;Xiong et al., 1997). The fourth-order correlations will appear; they are expressed as three products of two secondorder correlations of turbulent fluctuations in terms of the fourth-order correlation (Orszag, 1977;Lesieur, 1987). This standard normal approximation seems to have a more solid basis in stochastic theory. However, the physically positive quantities, such as the auto-correlations of turbulent velocity and temperature (namely x 2 and Z in the present paper), sometimes become negative. The numerical calculations will eventually fail to converge. Grossman (1996) compares in detail these two closure schemes (namely, our gradient-type diffusion approximation for third-order correlations and the quasi-normal approximations for fourth-order correlations). He concluded ". . . the Xiong solution predicts second moments better, with third-moment agreement not as good. The Full solution [namely with the quasi-normal approximation] predicts third moments better, but the second moments [namely x 2 , Z, and V in present paper] show inferior agreement. . . . Since the second moments are most important for constructing stellar models, we conclude that the Xiong closures perform impressively well." Therefore, we abandoned the closure schemes of the quasinormal approximations, and adopted the gradient-type diffusion approximation for the closure of the dynamic equations of turbulent correlations, because it is a simpler and more practical scheme than the quasi-normal approximation for calculations of stellar evolution and oscillations (Xiong et al., 1997).
Convection results from the thermal instability of a fluid medium. In a gravitationally stratified fluid medium, a perturbed fluid element will be accelerated by buoyant forces along the direction of gravity, when the local temperature gradient exceeds the adiabatic one. The original direction of convective motion is along with gravity, i.e., along the radial direction of the star. Therefore, convection is highly anisotropic in the low wave number range of the turbulent spectrum. Due to the continuity and nonlinearity of hydrodynamics, a part of the kinetic energy of convective elements will be converted into horizontal motion. Turbulence becomes more and more isotropic in the high wave-number range. Rotta (1951) pointed out that the correlation of the pressure and velocity gradients tends to make turbulence isotropic. Therefore, we assume (Deng et al., 2006) where x 2 and χ ij are, respectively, the isotropic and anisotropic component of the velocity correlation.
Substituting Eq. 17 into Eq. 9 and contracting with respect to indices i and j and noting Eqs. 10-16, we have 3 2 Subtracting the product of g ik and Eq. 18 from Eq. 9 and noting Eq. 16, we can obtain the dynamic equation of the anisotropic component of the velocity correlation χ ij , In the same way, we can derive the dynamic equations of the temperature auto-correlation and cross-correlation of velocity and temperature from Eqs. 7, 8: P e x/x c is the effective Peclet number of turbulent convection. Eqs. 3, 4 can be rewritten as where F k r , F k c , F k t are, respectively, the radiative, convective enthalpy, and turbulent kinetic energy fluxes, The contribution of pressure fluctuations to the enthalpy flux Eq. 26 has been neglected, because it is much less than the contribution of the temperature fluctuation. The combination of the average hydrodynamic equations Eqs. 2, 39, 49, 5, and the dynamic equations of auto-and cross-correlations for turbulent velocity and temperature Eqs. 18-21 form a complete and closed system of radiation-hydrodynamic equations for the calculation of stellar structure and oscillations.
The details of the derivation and main simplifications can be found in our original works (Xiong, 1978;Xiong, 1980;Xiong, 1981;Xiong, 1989;Xiong et al., 1997;Deng et al., 2006). The main simplifications and assumptions can be summarized as follows: 1) Turbulence is near quasi-isotropic; the characteristic linear size of turbulent elements is far less than the characteristic length of the average fields.
2) The relative fluctuations of temperature T′/T and density ρ ′ /ρ are much less than 1, so our theory is applicable only for subsonic convection.
3) The anelastic approximation (Gough, 1969;Xiong et al., 1997) is adopted. We assume that the only effect of pressure fluctuations is to make turbulence more isotropic (Eq. 16), and all of the other dynamic effects of pressure fluctuations were neglected. These assumptions are slightly weaker than the Boussinesq approximation. In our theory, the density-weighted turbulent velocity is used, so the compressibility is partially taken into account. However, noise sound waves are filtered due to the assumption that D(ρ ′ /ρ)/Dt 0 4) The turbulent fluctuations of the gravitational potential Φ′ have been neglected because they cancel out between each other, so gravity waves are also filtered. 5) The gradient-type diffusion approximation for the treatment of the third-order correlations is adopted (Xiong, 1980;Xiong, 1981;Xiong, 1989).
The turbulent dissipation, diffusion, and anisotropy are carefully taken into account. Based on the turbulence theory, they are expressed, respectively, by Eqs. 10, 12, 16, 19. There are three convective parameters (C 1 , C 2 , C 3 ) in our convection theory. They are associated, respectively, with turbulent dissipation, diffusion, and anisotropy. They can be calibrated by using the comparison between the observed and theoretical depth and T-P structure of the solar convection zone, the turbulence velocity-temperature fields in the solar atmosphere, the lithium depletion in the atmosphere of the Sun and late-type dwarfs, and 3D simulations (Deng et al., 2006). The calibration of the convective parameters (C 1 , C 2 , C 3 ) is somewhat complex. First, the convective parameters (C 1 , C 2 , C 3 ) are not a set of constants, but they vary slowly as a function of stellar parameters such as mass M, luminosity L, effective temperature T e (Ludwig et al., 1999); even in the interior of the same star they vary with radius r. Fortunately, in the deep interior of stars, apart from the surface layer, convective transport of energy is very effective (P e x/x c ≫ 1). Therefore, the temperature gradient is very close to the adiabatic one, independent of the choice of convective parameters. The depth of the convection zone is almost defined by the structure of the surface super-adiabatic convection zone, where convective energy transport is inefficient (P e x/x c ≤ 1). Our research shows that, if only the pulsational stability is of concern, by using a fixed set of (C 1 , C 2 , C 3 ) calibrated to the Sun, all of the pulsational instability strips in the H-R diagram, such as the RR Lyrae (Xiong et al., 1998b), δ Scuti − c Doradus (Xiong et al., 2016), Mira (Xiong et al., 1998a), and LPV (Xiong et al., 2018) instability strips, can be defined. The effects of (C 1 , C 2 , C 3 ) on convection are not completely independent and the available observations for calibration are too rare, so it is difficult to calibrate (C 1 , C 2 , C 3 ) independently. After theoretical considerations, these entanglements for calibration of (C 1 , C 2 , C 3 ) can be removed, at least partially. C 3 , unlike C 1 and C 2 , is a rather independent convective parameter. It is related to the anisotropy of turbulent convection. In the deep interior of a convectively unstable zone, the ratio of the squared radial component of turbulent velocity to the horizontal one w 2 r / w 2 h ≈ (3 + C 3 )/ 2 C 3 . Anisotropy increases and convective overshooting decreases with decreasing C 3. C 3 3 is a reasonable value, which agrees well with the observations of turbulent velocity in the solar atmosphere, the lithium depletion of the Sun and late-type main-sequence stars, and hydrodynamic simulations (Deng et al., 2006;Deng and Xiong, 2008). C 1 and C 2 are two convective parameters related to turbulent dissipation and diffusion. It can be seen from Eqs. 14, 15 that l e and ∧ are, respectively, the dissipation and diffusion length of turbulence. So, one can expect that C 2 / C 1 is approximately a constant independent of the stellar parameters M, L, and T e . C 2 /C 1 1/2 is a good choice. Once C 2 /C 1 and C 3 are identified, C 1 becomes the only adjustable convective parameter. Our research shows that a non-local and anisotropic convective model of the Sun with (C 1 , C 2 , C 3 ) (0.64, 0.50, 3.0) can reproduce almost all of the observed characteristics of the solar convection zone very well (Deng et al., 2006;Xiong and Deng, 2001).

STRUCTURE OF THE SOLAR CONVECTION ZONE
Astronomy is a science based on the observations of astronomical objects. Stellar convection theory originated from the requirement for a treatment of convective Frontiers in Astronomy and Space Sciences | www.frontiersin.org May 2021 | Volume 7 | Article 438864 transport of energy and momentum, and of convective mixing of materials in the stellar interior. In the previous section, we described a non-local and time-dependent theory of convection. In this section we will present its application in studies of the solar convection zone. The Sun is the closest star to us, and is the only star that can provide high spatial resolution observations. Therefore, the Sun is an ideal natural laboratory for testing convection theory. By using our non-local and anisotropic convection theory described in the previous section, we calculated a model of the solar convective envelope (Unno, et al., 1985;Xiong and Cheng, 1992;Xiong and Deng, 2001). The MHD equation of state Hummer and Mihalas, 1988; and OPAL opacity (Rogers and lglesias, 1992) supplemented by low-temperature opacities (Alexander and Ferguson, 1994) were used. The chemical abundance is X 0.70 and Z 0.02. The surface boundary was located at optical depth τ 10 −3 . By choosing a first trial set of T and r at the surface boundary T T 0 , r R 0 , and P aT 4 0 , by using a relaxation procedure, it is not difficult to obtain that T T e , r R ⊙ , and L 4πR 2 ⊙ T 4 e , at τ 2/3. The original fundamental equations are our radiation-hydrodynamic equations for the calculation of stellar structure and oscillations, Eqs. 2, 39, 49, 5, 18-21. The Henyey method (Henyey et al., 1964) for the integration of differential equations was used. The details of the working equations and boundary conditions can be found by referring to the original works of the author cited above.

Structure of Turbulent Velocity and Temperature Fields
Convection stops suddenly at the boundary of the convection zone in the local convection theory; however, convection penetrates deeply into the convectively stable zone in the non-local convection model. Figure 1 shows the variation of x, Z, |V|, and χ 1 1 vs. log P in our non-local convection model. V and χ 1 1 change their sign passing through the boundary of the convection zone. In the convectively unstable zone, V and χ 1 1 are greater than zero. Convective flux F C > 0 and the ratio of the squared radial component of turbulent velocity to the horizontal one w 2 r / w 2 h (x 2 + χ 1 1 )/(x 2 − χ 1 1 ) ≈ (3 + C 3 )/2C 3 > 1/2, when C 3 < 3. Turbulent motions dominate in the radial direction. In the overshooting zones, V and χ 1 1 become negative, (V < 0, χ 1 1 < 0), convective flux Fc < 0 and w 2 r / w 2 h (x 2 + χ 1 1 )/(x 2 − χ 1 1 ) < 1/2, and turbulent motions become dominant in the horizontal direction. In the surface overshooting zone, the correlation coefficient of turbulent velocity with temperature R VT V/xZ 1 2 ≈ − 1. Passing through the boundary of the convection zone, R VT changes abruptly from -1 to +1. In most of the convection zone, far from the boundary, R VT is very near to 1, and decreases rapidly to zero toward the bottom boundary of the convection zone. This FIGURE 1 | Auto-and cross-correlations of turbulent velocity and temperature, x, χ ij (Eq. 17), Z (Eq. 22), and V (Eq. 23) versus log P for the local (dotted-lines) and non-local (dashed-lines) convective models of the Sun having the same convection zone depth.
FIGURE 2 | Fractional convective flux L c /L, turbulent kinetic energy flux L t /L, super-adiabatic temperature gradient ∇ − ∇ ad , and the ratio of turbulent pressure to gas pressure P t /P vs. log P for the local (dotted-line) and non-local (solid-line) convective models.
Frontiers in Astronomy and Space Sciences | www.frontiersin.org May 2021 | Volume 7 | Article 438864 asymmetry of R VT at the surface and bottom boundaries of the convection zone results from the fact that, in the surface zone, the effective Peclet number P e x e /x C ≪ 1, and convective energy transport is very inefficient; however, in the deep interior of stars, P e ≫ 1, and convective energy transport is very efficient. Convective overshooting is different for the different physical variables. It can be seen that turbulent velocity and temperature penetrate deeply into the convectively stable zone. x, Z, and V decrease exponentially with ln P in the overshooting zone, and the e-folding length of turbulent velocity is about 1.4 C 1 C 2 √ times the pressure scale height. However, overshooting of the convective flux is negligible because the convective energy transport is inefficient (x/ x C ≪ 1) in the surface overshooting zone and the correlation coefficient of turbulent velocity and temperature is far less than 1 (R VT V/x Z (1/2) in the bottom overshooting zone). The overshooting distance decreases with decreasing C 3 and x/ x C . It can be seen that x, Z, and V decline more rapidly in the surface overshooting zone, where x/x C ≤ 1, than that in the bottom overshooting zone, where x/x C ≫ 1.

Remarks About MLT
The local MLT is still a good first approximation for the treatment of convection in the calculation of the thermal structure of stars, when the non-locality and time dependence of convection are not important for the problem concerned. Figure 3 illustrates log T vs. log P for the local and non-local models of the solar convective envelope, with the same convection zone depth. It can be seen that they agree with each other very well. Figure 4 illustrates the relative differences in the squared sound speed and density between the non-local and local convection models. Except for the regions closest to the surface, the relative differences are within approximately one percent. Figure 2 illustrates the fractional convective enthalpy and turbulent kinetic fluxes L C /L, L t /L (L C 4πr 2 F c , L t 4πr 2 F t and L L r + L C + L t ), the ratio of turbulent pressure to gas pressure P t / P, and the super-adiabatic temperature gradient ∇ − ∇ ad versus depth (log P). It can be seen that in the convectively unstable zone, away from the convection zone boundaries, these quantities are very nearly the same for the local and the non-local models, and the turbulent pressure and turbulent kinetic energy flux are negligible in comparison with the gas pressure and convective enthalpy flux. It is not difficult to prove that, neglecting turbulent pressure and the third-order correlation (i.e., turbulent kinetic energy flux), the dynamic equations of correlation functions Eqs. 18, 20, 21 will return to a form analogous to the local MLT formulae (Xiong, 1978;Xiong, 1989). It is not difficult to prove that ∇ ≈ ∇ ad + (∇ rad − ∇ ad )4x c /9x when x/x C ≫ 1, and ∇ ≈ ∇ rad − (∇ rad − ∇ ad )(x/x c ) 2 9/4 when x/x C ≪ 1. In the deep convection zone, where P e x/x C ≫ 1, convective energy transport is effective, and the temperature gradient is near the adiabatic one independent of the convection theory used. So, it is not difficult to understand why the T-P structures are so close to each other for the local and non-local convective envelope models in Figures 3, 4.

Structure of the Convective Overshooting Zone
Overshooting is a very natural phenomenon from the viewpoint of hydrodynamics. However, it was an outstanding problem for a FIGURE 3 | log T as a function of log P in the solar atmosphere for the non-local (dashed-line) and local convection models of the Sun (dotted-line). The inverse triangles are the Harvard-Smithsonian reference atmosphere (Gingerich et al., 1971). long time in the astronomical community. Up to now, the community still cannot cast off the influence of the MLT. The first non-local MLT is the generalized mixing-length theory by Spiegel (Spiegel, 1963). Ulrich (1970) proposed an analogous non-local MLT. Their theories were used to construct a model of the Sun (Ulrich, 1970;Travis and Matsushima, 1973). Observations show that the temperature gradient of their theoretical models were too gentle in the atmosphere of the Sun, because the efficiency of convective energy transport was overestimated by their theory. Figure 3 shows log T as a function of log P for our non-local and local convection models of the Sun. The inverse triangles mark the Harvard-Smithsonian reference atmosphere (Gingerich et al., 1971). The obvious feature of our non-local model is a larger temperature gradient at the top of the convection zone, which agrees with the Harvard-Smithsonian reference atmosphere very well. Helioseismology provides a direct method for probing solar interior structure. Figure 5 shows the results of the sound-speed inversion (Zhang et al., 2012). The open circles are relative differences in squared sound speed between the Sun and the standard solar model S (Christensen-Dalsgaard et al., 1996), which is a local MLT model. It can be seen that the relative difference of squared sound speed between the Sun and the reference model has a bump below the bottom boundary of the convection zone. This bump was understood as an indication of gravitational diffusion (Richard et al., 1996;Brun et al., 1999;Christensen-Dalsgaard, et al., 2007) for a long time. However, we already indicated long ago that this bump was not a result of gravitational diffusion, but was an indication of non-local convective overshooting (Xiong and Deng, 2001). We can see from Figure 2 that, in our non-local convection theory, the convective flux L C / L becomes negative and the fractional radiative flux L r / L ∇/∇ rad becomes slightly larger than 1, and the temperature gradient (∇) will be greater than the radiative one (∇ rad ) in the bottom overshooting zone. Therefore, the temperature in the overshooting zone will be higher than that predicted by the local theory of convection as shown in Figure 4. So, we predicted that this bump will be removed, or at least reduced to a great extent, if the reference model of inversion is replaced by our non-local model for the Sun (Xiong and Deng, 2001). This theoretical prediction has been confirmed by sound-speed inversion. In Figure 5 the filled circles are for our non-local model for the Sun as the reference model for inversion. The relative difference in squared sound speed is indeed reduced greatly in comparison with the standard solar model S (open circles).
Up to now the MLT is still a dominant idea in the astronomical community. All of the non-local MLT (Spiegel, 1963;Ulrich, 1970;Shaviv and Salpeter, 1973;Maeder, 1975;Bressan et al., 1981;Zahn, 1991) models still use a ballistic-type phenomenological theory. Figure 6A illustrates a sketch of the structure for the bottom overshooting zone of the Sun in the nonlocal MLT (Monteiro et al., 2000). In the convection zone, the temperature gradient is close to and slightly higher than the adiabatic gradient (∇ − ∇ ad ≥ 0). There is an overshooting zone under the convection zone, where the temperature gradient is near to and slightly lower than the adiabatic gradient. Passing through a very thin transition layer, the temperature gradient jumps from the adiabatic value to the radiative one. So, there is a near discontinuity of the second derivative of sound speed at the bottom boundary of the solar convection zone. This discontinuity will reflect incident acoustic waves and induce an oscillation component in the frequencies of p-modes as a function of radial order n (Gough, 1990). However, helioseismology observations show that the Sun has a mostly smooth stratification (Gough and Sekii, 1993), or the discontinuity is very small (Basu and Antia, 1994;Monteiro et al., 1994;Roxburgh and Vorontsov, 1994). Figure 6B shows the temperature gradient as a function of fractional radius around the bottom boundary of the convection zone predicted by our non-local convection theory (Xiong and Deng, 2001). There is not any discontinuity in the temperature gradient, and it transforms smoothly from the adiabatic one in the convectively unstable zone into the radiative one in the convectively stable zone. In our non-local convection theory, the temperature gradient has already become sub-adiabatic (∇ < ∇ ad ) in the convectively unstable zone before reaching the bottom boundary of the convection zone. It becomes sub-adiabatic and super-radiative (∇ rad < ∇ < ∇ ad ) in the overshooting zone below, until it approaches the radiative value (∇ ∇ rad ) in the deep radiative zone. It can be seen from Figures 6A,B that the structure of the transition zone around the boundary of the convection zone is very different between our non-local convection theory and MLT, which reflects the profound difference between two convection theories in the understanding and treatment of convection. Recently, Christensen-Dalsgaard et al. (2011) revisited this problem. They support our viewpoint of a smoothly stratified bridge from the region of the lower convection zone to the radiative interior, such as shown in Figure 6B, and believe that it will be in better agreement with helioseismic data than that of the standard solar model.
Eq. 18 is the dynamic equation for the isotropic component of the auto-correlation of turbulent velocity; it can be also understood as the conservation equation of turbulent kinetic energy. The first term is the growth rate of turbulent kinetic energy, which could be set equal to (taking into account minus signs) the sum of the rest of the terms in the equation. The second and third terms are the transformation rate between the average motion and the turbulent kinetic energy. The fourth term is the net gain from the non-local transport of turbulent kinetic energy, which is the key term distinguishing our non-local convection theory from the local theory. The fifth term is the boundary work W b , which is directly proportional to the convective flux F C ( ρC P TV) and represents the transformation rate between turbulent kinetic energy and thermal energy. The convectively unstable zone, where F C and W b > 0, is the driving region for turbulent convection, where the thermal energy is transformed into turbulent kinetic energy by buoyancy work. On the other hand, the convectively stable (including overshooting) zone is the damping zone for turbulent convection, where the turbulent kinetic energy is transformed into thermal energy due to buoyancy work. Therefore, it is reasonable to define F C 0 as the boundary of the convection zone (Xiong and Deng, 2001;Deng and Xiong, 2008).
During the past 20 years, rapid progress has been made in numerical simulations, and these have become an important means to study stellar convection. We have tried to use 2D or 3D simulations to test some basic assumptions of our non-local and time-dependent theory of convection and to calibrate the convection parameters (Kupka, 2002;Deng et al., 2006;Kupka and Muthsam, 2007;Kupka, 2007a;Kupka, 2007b;Kupka, 2007c;Cai, 2008;Tian et al., 2009;Chan et al., 2008;Chan et al., 2011;Cai and Chan, 2012;Cai, 2014;Kupka, 2017;Cai, 2018;Cai, 2020a;Cai, 2020b;Cai, 2020c). Some meaningful results were obtained. Because of limitations on the number of figures, these results cannot be described in detail. Interested readers can refer to the original texts, as well as to extensive reviews by Asplund et al. (2009) andNordlund et al. (2009).

OVERSHOOTING MIXING AND STELLAR EVOLUTION
We know from the theory of stellar structure and evolution that the stellar structure should be defined once the mass and element abundances of a star are assigned. Therefore, stellar evolution in the H-R diagram is, in fact, a reflection of nuclear evolution in the stellar interior. Convection (including overshooting) is the most important means of element mixing in the stellar interior. Our non-local convection theory has been successfully used to treat non-local convective mixing of elements in the evolution of massive stars (Xiong, 1986) and to calculate lithium depletion in the atmosphere of the Sun (Xiong and Deng, 2002) and latetype main sequence stars ). Schwarzschild and Härm (1958) have shown that the hydrogenrich radiative envelope adjacent to the helium-rich convective core cannot be stable against convection for massive stars (M ≥ 9 M ⊙ ). This is the familiar so-called semi-convection contradiction. Schwarzchild and Härm continued the local treatment of convection and introduced a so-called semiconvective zone in order to overcome this contradiction. However, there is not a self-consistent method for the construction of this semi-convection zone. Various researchers constructed their respective models of the semi-convection zone using their own methods. A review of these methods was given by Stothers (1970). He shows that at least ten schemes are available in the literature to treat this unstable intermediate zone. The different treatments of semi-convection resulted in some discrepancies in the evolution of massive stars. The semiconvection contradiction, in our opinion, results from the local treatment of convection; it will be removed automatically if a non-local treatment of convection is applied. In order to treat non-local mixing of nuclear fuels in stellar evolution, we need to develop a non-local theory of convection for chemically inhomogeneous stars (Xiong, 1981). An advantage of our dynamic theory of correlation functions is its convenience for generalization. For a chemically inhomogeneous star, apart from the conservation of total mass, Eq. 2, the conservation equation for each type of nuclear fuel should be added,

Evolution of Massive Stars
where C n is the content per gram of the nth nuclear fuel, q n and J k n are, respectively, its destruction rate and the molecular diffusion flux for the nth nuclear fuel. By using the same method of Reynolds decomposition and after an averaging procedure, mentioned in A Non-Local and Time-Dependent Theory of Convection section, it is not difficult to obtain the average conservation equation of fractional mass, and the corresponding dynamic equations of auto-and crosscorrelations C ′2 n ,u ′ k C′n and C′nT′/T . Eqs. 2,39,49,18,19,20,21,29, and the dynamic equations for C ′2 n ,u ′ k C′n and C′nT′/T form a set of equations for the calculation of stellar evolution (Xiong, 1981;Xiong, 1986).
By using our non-local theory of convection for chemically inhomogeneous stars mentioned above (Xiong, 1981), we calculated the evolution of massive stars in the hydrogenburning stages (Xiong, 1986). The semi-convection contradiction was indeed removed automatically, as predicted by us. Figure 7 shows the outline of the hydrogen content at various evolution ages for a 60 M ⊙ star. It can be seen that non-local mixing of helium penetrates deeply into the convectively stable zone beyond the boundary of the convectively unstable core marked by the dashed line. A molecular-weight gradient region adjacent to the convective core formed automatically. The convective core size increases due to non-local overshooting mixing, so the evolutionary track runs at a higher luminosity for a star with fixed mass, and the main sequence band becomes noticeably wider and main sequence lifetimes become longer in comparison with those calculated using the local MLT (the dotted lines, Maeder, 1981). It can be seen from Figure 8 that these influences of non-local convective mixing on the structure and evolution of stars increase toward lower mass, because the relative increase of convective core size increases with the decreasing mass of the star. These theoretical predictions have been confirmed by subsequent research (Stothers, 1991;Chiosi, et al., 1992;Schootemeijer, et al., 2019;Higgins and Vink, 2020). This means that the evolution masses for a given luminosity were overestimated and the width of the main sequence band was underestimated when convective overshooting mixing was neglected. This might be one of the important causes for the Cepheid mass discrepancy (Christy, 1968;Stobie, 1969;Cogan, 1970;Rodgers, 1970) and the contradiction between the theoretical and observed distribution of luminous stars in the H-R diagram (Humphreys, 1979;Humphreys and Davidson, 1979). The MLT models can be made to agree with observations by adjusting core overshooting and semi-convection parameters, FIGURE 7 | Outlines of hydrogen content at various evolution ages for a 60 M ⊙ star. The dashed line is the boundary of the convectively unstable core, and the dotted line indicates the location of the nuclear energy generation rate ε 1000 erg/s (Xiong, 1986).
FIGURE 8 | Theoretical evolutionary tracks in the H-R diagram for 7-60 M ⊙ stars (Xiong, 1986). The solid and dashed lines are for our non-local convection theory with two different convective parameters, and the dotted lines are for the local MLT theory (Maeder, 1981 which makes such models less predictive, in contrast to our models. How one defines the boundary of convective zone is another contributor to the uncertainty in the treatment of non-local convection overshooting, and therefore to the uncertainty in the evolution. For a long time, a disputed question has been whether the criterion for convective instability should follow the Schwarzschild criteria or the Ledoux criteria. Our research shows that the Schwarzschild and Ledoux criteria are only applicable for local convective instability. In our non-local convection theory, the convective instability criterion follows neither Schwarzschild nor Ledoux. F C 0 is a more convenient and reasonable definition for the boundary of the convection zone, as described in Structure of the Convective Overshooting Zone section, and the details can be found in our earlier works (Xiong, 1986;Xiong and Deng, 2001;Deng and Xiong, 2008). The difficulty in determining the boundary of the convective core results from the local treatment of convection. This difficulty is also removed automatically. There is no longer any ambiguity in defining the convective core and constructing the overshooting zone in our non-local convection theory. Therefore, the uncertainties in massive star evolution due to ambiguous semi-convection and overshooting do not arise.  Chiosi et al. studied the mass discrepancy of the Cepheids. They found that "the mass discrepancy problem likely originates from the adoption of semi-convective models and insufficient accuracy in the determination of the mass by one of the two methods. When this is feasible, as in the ideal laboratory given by the young LMC clusters with Cepheids, the discrepancy no longer exists" (Chiosi, et al., 1992). By using the new and improved equation of state Hummer and Mihalas, 1988; and opacity (Roger and Iglesias, 1992), not only were the excitation mechanism of β Cephei and SPB stars explained (Cox et al., 1992;Moskalik and Dziembowski, 1992;Dziembowski and Pamyatnykh, 1993), but also the mass discrepancy between the evolutionary mass and the "bump" and "beat" masses of the Cepheids was reduced further (Carson and Stothers, 1988;Moskalik and Dziembowski, 1992;Petersen, 1992;Simon, 1995;Guzik et al., 2020).

Lithium Depletion in Late-Type Main Sequence Stars
Lithium abundance was a very important problem for nucleosynthesis during the Big Bang. In addition, lithium and beryllium are both very fragile elements, which are destroyed quickly due to nuclear reaction at T∼2.5 and 4.0 million K in the stellar interior, which makes them ideal trace elements for measuring the depth of the surface convection zone during stellar evolution. Since the original observational detection of the lithium abundance in solar-type stars by Herbig (1965), rich observational data have been accumulated. Observations of the lithium abundance of Galactic open clusters of various ages and metallicities provide very conclusive constraints on the depletion mechanism of lithium. Up to now, the depletion mechanism of lithium in the Sun and stars is still not fully understood. Various depletion mechanisms have been proposed, including mass loss (Weymann and Sears, 1965;Hobbs, et al., 1989;Schramm et al., 1990) and wavedriven mixing (Garcia Lopez and Spruit, 1991;Montalban and Schatzman, 1996), among which rotationally induced mixing (Pinsonneault et al., 1989;Charbonnel et al., 1992;Chaboyer et al., 1995) has been a dominant viewpoint. In our view, convective overshooting mixing and gravitational settling seems to be the most reasonable mechanism for lithium depletion.
By using our non-local theory of convection in chemically inhomogeneous stars described in A Non-Local and Time-Dependent Theory of Convection section and the previous subsection, we calculated the lithium depletion in the atmosphere of the Sun (Xiong and Deng, 2002) and late-type dwarfs . Apart from convective overshooting, the gravitational settling has been taken into account, because the timescale for gravitational settling in the atmosphere of warm stars becomes comparable to and even shorter than the evolutionary timescale. Referring to Chapman and Cowling (1970), the micro-diffusion flux J of lithium can be expressed as ) where D and K T are, respectively, the diffusion and thermal diffusion coefficients, μ o the mean molecular weight and φ the mean ionization degree per ion. The first, second, and third terms in Eq. 30 are, respectively, pressure, gravitational, and thermal diffusion. Figure 9A shows the evolution of surface lithium abundance with age for stars of different masses. We can see that lithium abundance tends to decrease exponentially with time. Figure 9B shows the e-folding time of lithium depletion as a function of stellar mass. The dotted lines are for the models in which only convective overshooting mixing is taken into account and the gravitational settling is neglected; the solid lines are for models with both convective overshooting mixing and gravitational settling. It can be seen that, for larger stars (M ≥ 1.1 M ⊙ ), lithium depletion mainly results from gravitational settling, because their surface convection zones are too shallow for the overshooting mixing to be efficient. The mass of surface convection zone and the timescale of gravitational settling increases with decreasing stellar mass. The effect of convective overshooting mixing on lithium depletion increases rapidly with decreasing stellar mass ( Figure 9B), and convective overshooting mixing becomes the main mechanism of lithium depletion for low-mass stars (M < 1.1 M ⊙ ) Figure 10 shows the lithium abundance versus depth (log ΔM r /M ⊙ ) at different ages (indicated on the curves) for 1.5 M ⊙ and 1.0 M ⊙ stars, clearly showing the physical picture for lithium depletion induced by overshooting mixing and gravitational settling. For moderately high-mass warm (T e ≥ 6200 K) stars, which have a shallow surface convection zone, surface lithium is drawn into the deeper radiative interior by gravitational settling and stays there, and the surface lithium abundance decreases with time ( Figure 10A). For low-mass stars (M ≤ 1.0 M ⊙ ), which have an extensive convective envelope, the surface lithium is brought into high-temperature interior regions by convective overshooting and destroyed by nuclear reactions there, and the surface lithium abundance decreases with time ( Figure 10B). Therefore, the mechanism of lithium depletion is very different for warm (M ≥ 1.1 M ⊙ ) and cool (M ≤ 1.0 M ⊙ ) main sequence stars. Lithium depletion is very sensitive to the treatment of convection, because lithium depletion increases exponentially with time. A small error will be enlarged exponentially with time. Therefore, observations of lithium depletion provide the most rigorous constraint for stellar convection theory. Figures 11A-D show the distribution of lithium abundance with effective temperatures for the member stars of several Galactic open clusters with different ages; the theoretical isochrones are also drawn on the figures. It can be seen that the theoretical predictions, in general, agree well with the observations for Galactic open clusters over a very wide range of ages. This means that our convection theory and the assumption of the convective overshooting mechanism for lithium depletion seem to be true.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.

FUNDING
This work was supported by National Natual Science Foundation of China (NSFC) through grants 11373069.