Issue 
A&A
Volume 588, April 2016



Article Number  A122  
Number of page(s)  21  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201527663  
Published online  30 March 2016 
Generation of internal gravity waves by penetrative convection
LESIA, Observatoire de Paris, PSL Research University, CNRS, Université Pierre et Marie Curie, Université Paris Diderot, 92195 Meudon, France
email: charly.pincon@obspm.fr
Received: 29 October 2015
Accepted: 13 December 2015
Context. The rich harvest of seismic observations over the past decade provides evidence of angular momentum redistribution in stellar interiors that is not reproduced by current evolution codes. In this context, transport by internal gravity waves can play a role and could explain discrepancies between theory and observations.
Aims. The efficiency of the transport of angular momentum by waves depends on their driving mechanism. While excitation by turbulence throughout the convective zone has already been investigated, we know that penetrative convection into the stably stratified radiative zone can also generate internal gravity waves. Therefore, we aim at developing a semianalytical model to estimate the generation of IGW by penetrative plumes below an upper convective envelope. The formalism is developed with the purpose of being implemented in 1D stellar evolution codes.
Methods. We derive the wave amplitude considering the pressure exerted by an ensemble of plumes on the interface between the radiative and convective zones as source term in the equation of momentum. We consider the effect of a thermal transition from a convective gradient to a radiative one on the transmission of the wave into the radiative zone. The plumeinduced wave energy flux at the top of the radiative zone is computed for a solar model and is compared to the turbulenceinduced one.
Results. We show that, for the solar case, penetrative convection generates waves more efficiently than turbulence and that plumeinduced waves can modify the internal rotation rate on shorter time scales. The result is solid since it holds despite a wide range of values considered for the parameters of the model. We also show that a smooth thermal transition significatively enhances the wave transmission compared to the case of a steep transition.
Conclusions. Driving by penetrative convection must be taken into account as much as turbulenceinduced waves for the transport of internal angular momentum. We propose a simple prescription that has the advantage of being easily implementable into 1D stellar evolution codes. We expect this mechanism to work in evolved stars, which will be subject to future investigations.
Key words: hydrodynamics / stars: interiors / stars: rotation / waves / convection
© ESO, 2016
1. Introduction
Rotation is a fundamental ingredient of stellar evolution. It induces transport of angular momentum and of chemical elements that modifies the internal structure of stars (e.g., Maeder 2009). For instance, rotationally induced mixing is able to refresh nuclear burning cores with hydrogen and thus to increase the star lifetime on the main sequence. Among other consequences, stellar agedating is significantly affected by rotation (e.g., Lebreton et al. 2014) such that the rotational history of stars has to be fully understood to obtain a relevant and complete picture of their evolution.
This requirement is supported by an increasing number of observational facts. First, indirect constraints, such as anomalies in surface chemical abundances observed in stellar clusters, highlighted the importance of rotationally induced mixing (see Charbonnel & Talon 2008, for a review). Second, the development of asteroseismology made direct measurements of the internal rotation profile possible. In the Sun, seismic measurements showed that its internal radiative zone rotates almost as a solid body (e.g., Brown et al. 1989; García et al. 2007). More recently, seismic data provided by the spaceborne missions CoRoT (Baglin et al. 2006a,b) and Kepler (Borucki et al. 2010) enabled us to extend the study from the main sequence up to the red giant evolutionary phase. Thanks to the detection of mixedmodes (e.g., Bedding et al. 2010; Mosser et al. 2012b), it has been possible to estimate the core rotation rate of several subgiant stars observed by Kepler (Deheuvels et al. 2012, 2014) as well as the spinning down of the red giants core during their ascent of the vertical branch (Mosser et al. 2012a).
All these observational constraints emphasized the need for including transport of angular momentum as well as rotational mixing in stellar models. However, stellar evolution codes that take meridional circulation and shearinduced turbulence into account are unable to reproduce the observations by several orders of magnitude (Eggenberger et al. 2012; Marques et al. 2013; Ceillier et al. 2013). This suggests that other efficient transport processes must be included. Several mechanisms have already been investigated such as the effects of magnetic fields (Spruit 2002; Heger et al. 2005; Cantiello et al. 2014; Rüdiger et al. 2015) or mixed modes (Belkacem et al. 2015a,b). Internal gravity waves (hereafter IGW) can also play a role in the radiative zone of stars. They have been shown to be able to explain the nearly flat internal rotation profile observed in the Sun (Zahn 1997; Talon et al. 2002) and to be responsible for the cold side of the lithium dip observed in lowmass stars (Talon & Charbonnel 1998a,b, 2003, 2004, 2005).
The efficiency of the transport by IGW crucially depends on the excitation mechanism. IGW are preferentially generated at the boundary between the convective and radiative regions by convective motions. Press (1981) considered turbulent pressure at the base of the convective zone as driving process. He found that the wave energy flux could be expressed as the product of the mechanical convective energy flux at the base of the convection region with the ratio of the wave impedances in both regions (i.e., in the convective and the radiative regions). Garcia Lopez & Spruit (1991) followed a similar approach but considered a distribution of convective eddies. They assumed the turbulence to follow a Kolmogorov spectrum and took into account the incoherent behavior of the eddies. Later, Zahn (1997) generalized this model to a continuous wave spectrum. Finally, Kumar et al. (1999) adapted the work of Goldreich et al. (1994) to the case of the excitation of IGW. In this model, IGW are generated by Reynold’s stress throughout the convective zone and the approach has the advantage of clearly taking into account the spatial and temporal correlations between the turbulent stochastic source and the waves. Most of the estimates of the transport by IGW used the latter formulation to include the effects of IGW in stellar radiative zones (e.g., Talon & Charbonnel 2005; Charbonnel et al. 2013; Mathis et al. 2013; Fuller et al. 2014).
However, penetration of convective plumes into stably stratified layers can also generate IGW. Indeed, in the penetration region, plumes are decelerated by buoyancy braking and can release a part of their kinetic energy into waves. In the geophysical context, this mechanism has already been observed in laboratory experiments and investigated theoretically for atmospheric flows (Townsend 1966; Stull 1976). In astrophysics, excitation of internal waves by penetrative convection has been investigated by means of 2D numerical simulations (Hurlburt et al. 1986; Andersen 1994; Dintrans et al. 2005; Rogers & Glatzmaier 2006; Rogers et al. 2013) and more recently, in 3D, spherical numerical simulations for the Sun (Brun et al. 2011; Alvan et al. 2014, 2015). Nevertheless, it is difficult today to extrapolate numerical results to stellar regimes. For instance, numerical constraints require higher thermal diffusivities than in stellar interiors, resulting in much more important convective energy flux (Rempel 2004; Dintrans et al. 2005) and so an expected overestimated wave energy flux.
In this work, we propose a complementary approach and elaborate a semianalytical model of excitation of IGW by penetrative convection. We consider the impact of the plumes on the stably stratified layers as the driving mechanism. The effect of the thermal transition on the transmission of the wave into the radiative zone is taken into account at the interface between the radiative and convective regions (hereafter, radiative/convective interface). Our goal is to estimate how efficient is this mechanism compared to the excitation by turbulence as proposed by Kumar et al. (1999). If it is shown to be as efficient or more than the latter, such a simplified description will enable us to account for the transport of angular momentum by plumeinduced waves in 1D stellar evolution codes.
The article is organized as follows. In Sect. 2, we introduce the theoretical formalism. Section 3 presents the characteristics of the plumes and describes buoyancy braking process in the penetration zone. In Sect. 4, we summarize the method used to derive the wave energy flux and the main results are given. The wave energy flux is computed for a solar model in Sect. 5. The results are discussed and compared to those obtained with the formalism of Kumar et al. (1999) in Sect. 6. Conclusions are formulated in Sect. 7.
2. Theoretical formalism
Our objective is to estimate the amount of plume kinetic energy transferred into wave energy in the penetration region. For this purpose, we introduce the theoretical formalism and general concepts regarding IGW and spectral analysis. In the following, we will deal with a stationary random process of excitation by an ensemble of plumes and will use a statistical approach.
2.1. Wave equation and source term
The determination of the wave energy flux requires the calculation of the wave amplitude, which depends on the excitation mechanism. As convective plumes penetrate into the stably stratified medium, they perturbate the equilibrium state and generate waves by exerting pressure on the radiative/convective interface. We assume that the total velocity field can be split into two components : the wave velocity field v(r,t), and the plume velocity field V_{p}(r,t). This actually corresponds to neglecting the dynamical effect of the turbulence inside the plumes. Moreover, we assume that the plumes evolve independently of the wave velocity and that there is no feedback from the waves on the plumes. This also suggests that the plumes and the waves follow their own continuity equation independently of each other. This is actually a good approximation if the wave velocity field is much smaller than the plume velocity field  v  ≪  V_{p} , which will be verified a posteriori.
Under all this set of approximations, we will then focus our attention on the generation of waves by the stress exerted by the convective plumes, represented by ∇·(ρV_{p} ⊗ V_{p}) as the source term in the momentum equation, with ρ the density at the equilibrium, ∇ the nabla operator and (⊗) the tensorial product. For sake’s of simplicity, we neglect any effects of rotation and adopt the Cowling approximation. Therefore, the linearized equations of momentum and of continuity with respect to the equilibrium state for the wave perturbations read where p′ and ρ′ denotes Eulerian perturbations of pressure and density, δρ is the Lagrangian perturbation of density and g the gravitational acceleration at equilibrium.
2.2. Spectral density of the wave specific energy
In the following, we want to determine the spectral distribution of the wave energy flux at each level in the star, generated by a population of penetrative plumes.
2.2.1. Stationarity and ergodicity
Fig. 1 Schematic view of the star. The convective plumes occur in the upper layers of the star and go deeper into the convective region. They grow by turbulent entrainment of matter at their edges and reach the top of the radiative zone. There, each of them, characterized by their angular position (θ_{i},φ_{i}), releases a part of their kinetic energy and generate internal waves which can propagate towards the center. In our work, we suppose that there are no reflected waves propagating from the center to the radiative/convective interface. 

Open with DEXTER 
We assume that at each time, an ensemble of several uniformly distributed incoherent plumes generate waves by penetrative convection (see Fig. 1, for a schematic representation). The process of excitation is then supposed to be random, stationary and ergodic, and a statistical approach is valid since the number of plumes is high enough. Hence, we can define the mean specific wave energy at a point r as (3)where v_{T}(r,t) is equal to the truncated part of the wave velocity field in the interval [− T/ 2,T/ 2] and is null outside. The ParsevalPlancherel’s theorem allows us to write (4)where ω is the temporal radian frequency and the symbol ( ) denotes the time Fourier transform^{1}. Then, using Eq. (4)and since the integral and the limit commute, Eq. (3)becomes (5)where we have defined the spectral density of the wave specific energy (6)
2.2.2. Ensemble of plumes
At each time, we assume that, on average, identical plumes are penetrating into the stably stratified zone. We note (θ_{i},φ_{i}) the angular position of the center of the plume i, that penetrates at the time t_{i} with a velocity field V_{p,i}(r,t). Since the plumes are spatially and temporally well separated with each other, the source term in Eq. (1)can be rewritten as (7)Each plume i generates a wave velocity field v_{i}(r,t;θ_{i},φ_{i},Δt = t_{i}) where Δt denotes the time shift of the excitation event compared to the instant t = 0. Given the linearity of the wave equation, the total wave velocity field is the superposition of all the contributions. In a similar way, at a fixed position r, the truncated wave velocity field, v_{T}, in the interval [− T/ 2,T/ 2], must be the superposition of the wave velocity fields generated by a finite number of plumes at different times^{2}(8)When writing Eq. (8), we assume that the plumeinduced wave packet v_{i} at any point r has a finite lifetime (otherwise, this would lead to an accumulation of wave energy in the radiative zone under the adiabatic hypothesis). Indeed, as suggested by the momentum equation Eq. (1), it must be close to the characteristic plume lifetime; this is the consequence of the temporal correlation with the plume and the fact that no reflected wave is considered, but only propagative ones towards the center of the star.
By the linearity and the time shift properties of Fourier transform, we then obtain^{3}(9)where the bar denotes the complex conjugate. The statistical average covers the spatial and temporal distribution of the plumes. The plumes are incoherent with each other and so the phase lag between each component is randomly distributed. Note that the same assumptions were first used by Garcia Lopez & Spruit (1991). As a consequence, the second term in the righthand side of Eq. (9)vanishes. Moreover, the plumes being uniformly distributed and independent with each other, the probability for the plumes to be located in the solid angles , respectively, is equal to . Thus, using Eq. (9), Eq. (6)becomes (10)Since we suppose that the plumes are identical and differ from each other only by their angular position, Eq. (10)becomes (11)where is the Fourier transform of the wave velocity field generated by one single plume with the angular position (θ_{0},φ_{0}) at t_{0} = 0.
To go further, the plume destruction rate must equal the plume emerging rate in order to ensure a constant number of penetrating plumes over time. If τ_{p} denotes the characteristic plume lifetime and if it is supposed to be the same for all of them, this rate is equal to . In other words, over a time T ≫ τ_{p}, a total number of plumes occur and contribute to the wave energy in Eq. (11). Finally, the spectral density of specific wave energy converges on average to (12)where ν_{p} = 1 /τ_{p}.
We first notice that the spectral density of wave specific energy is proportional to the number of plumes, which is the consequence of their incoherent behavior (in the case of coherent plumes, the wave specific energy would be proportional to ). Second, we see that the use of a statistical approach simplifies the calculation since the case of an ensemble of several plumes is reduced to the study of the excitation by one single plume.
2.3. Mean radial wave energy flux per unit of frequency
For convenience, the Eulerian wave velocity field v_{0} is expanded onto the spherical harmonics (13)where (r,θ,φ) are the spherical coordinates, l is the angular degree, m the azimuthal number, e_{r} the radial unit vector, and where v_{r,l,m}, v_{P,l,m} and v_{T,l,m} are the radial, poloidal and toroidal components of the wave velocity field, respectively. Given the decomposition in Eq. (13), Eq. (12)averaged over the solid angle represents the angular distribution of energy and reads (14)where we have defined (15)with . The mean radial wave energy flux by unit of frequency for an angular degree l and an azimuthal number m is thus obtained by multiplying the (l,m) contribution to the spectral density of the specific wave energy with the radial group velocity V_{gr} in a similar way to Zahn (1997)(16)Indeed, IGW are dispersive waves and the energy carried in the radial direction by such a wave packet travels at first order with the radial group velocity (see Lighthill 1978) given by (17)where k_{r} is the local radial component of the wave vector (Press 1981) (18)with the horizontal component of the wave vector and N the BruntVäisälä frequency. The total mean radial wave energy flux per unit of frequency due to the penetration of several plumes is then obtained by summing all the (l,m) components given by Eq. (16).
3. A simple description of convective plumes
The computation of the wave energy flux requires the knowledge of the wave velocity field whose the amplitude directly depends on the excitation mechanism. We then need to model the driving process expressed by the plumerelated term in Eq. (1). For this purpose, we present here a simple physical description of convective plumes and their velocity profile in the penetration zone.
3.1. Plume velocity field
The disturbance of the radiative/convective interface due to a single plume is localized in space and in time. Therefore, a plume is described by a characteristic radius b and a characteristic lifetime τ_{p} in the penetration region. We assume that the plume velocity field follows the Gaussian form^{4} proposed by Townsend (1966) and we also assume that the horizontal profile is maintained in the penetration region (19)with (20)where (θ_{0},φ_{0}) are the angular coordinates of the center of the plume and S_{h} corresponds to the distance on the sphere from the center of the plume.
3.2. Plume radius
The plume radius at the bottom of the convective zone is estimated from the model of turbulent plumes by Rieutord & Zahn (1995) who derived the expression (21)where z_{0} is the thickness of the convective zone, Γ_{1} the adiabatic coefficient and α_{E} the entrainment coefficient whose value is usually taken as 0.083 (Turner 1986). In the solar case and for a monoatomic gas (Γ_{1} = 5/3), we find b ≈ 10^{4} km, which is about five times smaller than the size of the biggest turbulent eddies at this radius as given by the MLT.
3.3. Vertical velocity profile in the penetration zone
Once the plume has penetrated into the stably stratified layers, it is less dense than the surrounding medium and is slowed down by buoyancy braking. We adopt the description of Zahn (1991) who considers a stationary flow in balance with buoyancy. Such an approach is justified by the high Péclet number at the base of the convective zone. Zahn (1991) could estimate the penetration length L_{p} and showed that it depends on a filling factor, the radiative diffusivity scale height, the total energy flux and the plume kinetic energy at the base of the convective envelope. In the solar case, he found L_{p} ~ 0.5H_{p} ~ 10^{4} km. Nevertheless, using helioseismology, Basu (1997) set an upper limit for the overshoot at 0.05H_{p}, with H_{p} the pressure scale height. In the framework of the model proposed by Zahn (1991), Basu (1997) used stellar models with a simple adiabatic extent of the convective zone above a discontinuous thermal transition. In contrast, ChristensenDalsgaard et al. (2011) demonstrated that smoother transitions were also possible to explain seismic observations of the Sun. Such a kind of transition is supported by more realistic models of penetrative convection taking into account interaction with the upflows and a distribution for the velocities of the plumes (e.g., Rempel 2004). Therefore, given the lack of knowledge and constraints about the penetration process, we will consider the penetration length as a parameter of the model whose value will be chosen around 0.1H_{p}.
Following the work of Zahn (1991), the vertical plume velocity V_{0} in the penetration region is given by^{5}(22)where z = r_{b}−r with r_{b} the radius of the base of the convective zone as given by the Schwarzschild’s criterion and V_{b} the initial vertical plume velocity in the penetration region, at the radius r_{b}. To go further, we use the model of Rieutord & Zahn (1995), so that (23)where F_{1p} represents the total luminosity (kinetic + enthalpic) carried by one single plume and ρ_{b} is the mean density at the base of the convection zone. To estimate F_{1p}, we assume that each plume carries an energy that is equal, on average, to the one carried by turbulent eddies in the interplume medium such that (24)where L_{∗} is the star luminosity and is the filling factor related to the fraction of the area occupied by the downdrafts at the base of the convective zone. Then, Eq. (23)becomes (25)which leads to V_{b} ≈ 185 m s in the case of the Sun, i.e about 7 times higher than the convective velocity v_{c} ~ 25 m s as given by the MLT at the base of the convective zone.
3.4. Plume lifetime
The plume lifetime is certainly the more difficult plume characteristic to estimate. After having reached the base of the convective zone, the plume penetrates into the stably stratified zone. There, it is slowed down and destroyed after a characteristic time τ_{p}. For sake’s of simplicity, we prefer here to tackle the issue by using orders of magnitude. We identify three potential processes working on the destruction of the plume:

1.
Radiative thermalization: in the penetration zone, the plumeloses its identity due to radiative thermal diffusion on a time scale,t_{rad}. It can be estimated by (26)where K_{b} is the radiative conductivity at the top of the stably stratified layer. In the Sun, we find t_{rad,⊙} ~ 10^{10} s, so it is much longer than the dynamical time scale t_{d, ⊙} ~ L_{p}/V_{b} ~ 10 s. As already mentioned, advection of adiabatic convective matter is faster than thermal exchange because of a high Péclet number at the base of the convective zone. Plume thermalization by radiative diffusion is thus inefficient, except in the transition region where the plume is slowed down enough for the dynamical time scale to be large enough.

2.
Turbulence inside the plume: although the plume is described like a coherent flow with welldefined radius and vertical profile, it is in fact turbulent at small scales as explained in Montalbán & Schatzman (2000) who suppose that turbulence is generated by shear flow in the penetration zone. In this case, the plume lifetime should be equal to a few turnover time scales of the biggest turbulent eddies whose the velocity and the size are about the ones of the plume, V_{b} and b, respectively. The turbulent time scale, t_{turb}, can be then approximated by (27)For the Sun, we find t_{turb,⊙} ~ 10^{5} s, which is much lower than the radiative time scale, but still one order of magnitude higher than dynamical time scale. Note that this value is close to the convective turnover time scale as predicted by the MLT, i.e., t_{c} ~ 10^{6} s at the base of the solar convective envelope.

3.
Restratification by lateral baroclinic eddies: finally, it is worth considering the restratification phenomenon observed in the terrestrial oceans (Jones & Marshall 1997) as a potential process operating in the stellar penetration region. In most locations in the oceans, the surface layers are stably stratified. However, convection driven by cooling at the seasurface can occur in some particular region. Then, convective plumes develop by turbulent entrainment of matter while going deep in the sea. When convection stops, the phase during which the plume is destroyed and loses its identity is called restratification. At this point, the density gradient between the plume and the surrounding stable medium generates a thermal wind. This latter is subject to a baroclinic instability creating lateral eddies able to transfer density between both mediums and thus able to homogenize the region. By adopting the model of Jones & Marshall (1997) to the solar case (see Appendix A for details), the restratification time scale is about t_{res} ~ 10^{6}−10^{7}s, which is in the order of magnitude of the convective eddy turnover.
The plume lifetime is likely to be the consequence of all these potential processes. The simple abovementioned estimates give reasonable values in the range of 10^{5}−10^{7} s for the Sun, i.e around the convective turnover time scale as given by the MLT. Thus, the plume lifetime will be considered as a parameter of the model whose the value will be chosen to be of the order of magnitude of the convective turnover time scale.
4. Generation of IGW by penetrative convection
4.1. Modeling the penetration zone
Fig. 2 Modeling of the interface between the radiative zone an the convective zone. r_{b} corresponds to the radius at the base of the convective envelope as given by the Schwarzschild’s criterion, L_{p} is the penetration length over which the plume is decelerated and generates waves, d denotes the extent of the transition region (where the plume velocity is negligible compared to V_{b}, its value at the entry of the penetration region) and r_{d} is the radius at the top of the radiative zone. 

Open with DEXTER 
Penetrative convection generates IGW at the radiative/convective interface and modifies the stellar equilibrium stratification. Therefore, we need to model this region to properly describe the excitation process. The scenario we consider in this work follows in part the description given by Zahn (1991). The situation is illustrated in Fig. 2. The plume reaches the base of the convective zone located at the radius r_{b} as prescribed by the Schwarzschild’s criterion. By inertia, the column of fluid penetrates into the stably stratified zone over a penetration length L_{p}. The density contrast between the plume and the surrounding medium causes buoyancy braking and the plume slows down, releasing a part of its kinetic energy into waves. Indeed, as explained by Zahn (1991), since the Péclet number at the bottom of the convective zone is very high (P_{e} ~ 10^{5}−10^{6}), advection of convective matter is more efficient than thermal exchanges so that the plume keeps its identity. Then, it imposes a quasiadiabatic stratification in the penetration zone, resulting in a nearby vanishing BruntVäisälä frequency N^{2} ≈ 0. Once its velocity is small enough (i.e., much smaller than its value at the entry of the penetration zone, V_{b}), the Péclet number decreases and the thermalization of the plume’s material occurs in a thermal transition region; there, the temperature gradient switches from quasiadiabatic one to radiative one over a distance d. In a first approach, we merely suppose that the BruntVäisälä frequency N^{2} increases linearly to , the value at the top of the radiative zone, as in the work of Lecoanet & Quataert (2013). By this way, we aim at investigating how the plumeinduced wave energy depends on the steepness of the transition.
4.2. Methodological approach
We now have to derive the time Fourier transform of the wave velocity field which leads to the spectral description of the wave energy flux through Eq. (6). For this purpose, we take the time Fourier transform of Eq. (1)and Eq. (2). We also assume that the oscillations are adiabatic so as to close the system. The firstorder system of two differential equations obtained can be rewritten in the following form (see Appendix B.1 and B.2 for a detailed derivation) (28)where we have defined where g is the gravitational acceleration, c is the sound speed, S_{l} is the Lamb frequency and d_{r} ≡ d/dr is the derivative with respect to r. The term F_{2} is defined by Eq. (B.12); it contains information on the spatial and temporal correlations between the wave and the plume, and we assume that it only exists in the penetration region where the wave driving is stronger. In this region, the inhomogeneous equations are solved using the method of variation of parameters and using the plume description presented in Sect. 3. In the convective and radiative zones, we use the WKB asymptotic solution for the homogeneous differential system, which is valid in these regions contrary to the transition region where Airy functions are considered.
To go further, we impose the continuity of the radial and horizontal displacements at the top of the radiative zone in r_{d} and at the beginning of the transition region in r_{d} + d. For boundary conditions, we first consider a pure regressive wave propagating towards the center in the radiative zone. That means we do not consider the possible reflection of the wave in the center of the star and so do not allow for standing waves to establish. Finally, in the convective zone, we consider a pure evanescent wave which is generated in the penetration zone and is damped towards the surface.
4.3. Wave energy flux
In the following, we summarize the main results of the detailed calculation given in Appendix B. The mean radial wave energy flux per unit of frequency generated by plumes for an angular degree l and an azimuthal number m reads (31)with where , dm = 4πρr^{2}dr, D ≈ 3.7 a numerical factor, s_{b} = S_{h}(r_{b},θ,φ;θ_{0},φ_{0}) following Eq. (20)and d is the length of the transition region (see Fig. 2) . We can first notice that the spectral density of wave luminosity, ℒ = 4πr^{2}ℱ_{r}, is conserved at each level in the star, which is the consequence of the adiabatic hypothesis for the waves in a nonrotating star. Now, let us discuss the different terms:

1.
ℋ_{l} corresponds to the wave driving term and it is representative of the instantaneous power injected into the wave in the whole penetration zone^{6}. It decreases with the angular degree l since the decay length of an evanescent wave in the penetration zone scales as Λ^{1}. In our model, it is also the only term which depends on the penetration length through the domain of integration and the plume velocity. If L_{p} decreases, the buoyancy braking is stronger and the excitation is more efficient since the plume deposits its energy where the evanescent wave has a higher amplitude. Inversely, if L_{p} increases, the plume energy is transferred into wave on a longer distance where the wave amplitude decreases rapidly.

2.
represents the temporal correlation term. The width of the spectral envelope is around ν_{p}, which means that the transit time at a point r of a wave packet generated by one single plume is around 1 /ν_{p} = τ_{p}, so that the assumptions made in Sect. 2.2.2 are verified. An increase of the characteristic plume lifetime causes an increase of the spectrum amplitude but a decrease of the spectrum width. In other words, the plume energy is transferred into higher frequencies at the expense of the lower ones, but the total wave energy is conserved since the integration of the flux over ω does not depend on ν_{p}. This is the consequence of the hypothesis of stationarity for the driving mechanism.

3.
expresses the horizontal spatial correlation between the wave and the plume.

4.
ℬ_{l} corresponds to the horizontal correlation term averaged over the angular position of the plume (θ_{0},φ_{0}). Since the plumes are supposed to be horizontally, uniformly distributed at each time, the spherical symmetry is conserved and ℬ_{l} does not depend on the azimuthal order m. This term decreases with l since the spherical harmonics presents more and more zeros with higher degrees. Moreover, a change of behavior is expected when the horizontal wavelength, λ_{h} = 2πr_{b}/ Λ, becomes smaller than plume radius, i.e., for a degree l_{crit} ~ 2πr_{b}/b. In the case of small degrees l ≪ l_{crit}, an approximated analytical expression in agreement with the numerical integration is given by (37)Since the wave energy flux per unit of frequency is proportional to Λℬ_{l} (neglecting the dependence of ℋ_{l} with l), Eq. (37)shows that the energy flux should be maximal for a degree l_{max} ~ r_{b}/b.

5.
f(γ_{d}) is a transmission function. It discriminates between the limiting case of a very sharp, almost discontinuous transition for γ_{d} ≪ 1 and the one of a smoother transition for γ_{d} ≫ 1. It physically corresponds to the ratio of the transmission coefficients in both cases. A very sharp transition has no effect on the wave dynamics. In this case, the flux is inversely proportional to the value of the BruntVäisälä frequency at the top of the radiative zone, N_{0}. That means that the higher is the step, the higher is the wave impedance in the radiative zone (Press 1981) and the smaller is the wave energy transferred. The smoother case occurs when the distance of the transition, d, is in the same order of magnitude or larger than the radial wavelength, so that the transition has an impact on the waves and improves its transmission into the propagative zone. In the special case of a linear transition profile, the transmission of the wave energy flux is enhanced by a multiplying factor (k_{r}d)^{1/3} compared to the discontinuous case, in agreement with Lecoanet & Quataert (2013). We highlight here the effect of the steepness of the transition on the flux; this latter depends on the slope of the BruntVäisälä frequency in the transition zone which plays the role of a pseudoimpedance (38)so that the smoother the transition, the higher the transmission of the wave energy.
4.4. Simplified expression for the wave energy flux
Finally, it is possible to derive an approximate expression for Eq. (31)in the whole radiative zone. This gives the advantage of expressing the result in terms of the characteristic physical quantities of the problem and it permits to avoid the numerical calculations of integrals. By using Eq. (37)for l ≪ l_{crit} and the assumption (which is valid for a penetration length much smaller than the characteristic decay length of the evanescent wave, i.e., L_{p} ≪ r_{b}/ Λ ), we find (39)where is the area occupied by a single plume, represents the plume kinetic energy flux and F_{R,l} = V_{b}k_{h,b}/N_{0} is the Froude number at the top of the base of the convective zone, with the horizontal wavenumber. This latter represents the ratio of the advection by the mean plume flow to the buoyancy force at the top of the radiative zone when the plume penetrates into. In the Sun, we find F_{R,1} ~ 10^{3} for l = 1. In the following Sect. 5, we will demonstrate the validity of this expression (for a penetration length small enough) since the major part of the wave energy flux is distributed over angular degrees below the critical angular degree l_{crit}. Therefore, the implementation of angular momentum transport by plumeinduced waves in a stellar evolution code will be easier using this simplified expression.
5. The solar case
In this section, we use a solar model to compute the wave energy flux given by Eq. (31). We analyze the effects of the plume characteristics on the excitation process and of the thermal transition length, d, on the wave transmission into the radiative zone (see Fig. 2).
5.1. Input physics
The calibrated solar model used here was obtained using the stellar evolution code CESTAM (Marques et al. 2013), with a mixing length parameter α_{MLT} = 1.69 and initial abundances Y_{0} = 0.25 and Z_{0} = 0.013. It provides the radial profile of all the quantities at the equilibrium that we need to compute the wave energy flux. The parameters computed from the model are the plume radius b, the plume velocity V_{b} and the BruntVäisälä frequency N_{0} at the bottom of the transition region (r = r_{d}). Thus, using Eq. (21)and Eq. (25)with ρ_{b} ≈ 150 kg m^{3} the density at the base of the convective zone and r_{b} = 5,1 × 10^{5} km the radius at the base of the convective zone, we find b = 10^{4} km and V_{b} ≈ 185m s^{1}. In most of stellar codes, standard mixinglength treatment of the convection is local. It produces a very sharp thermal transition at the radiative/convective interface and so a very sharp slope for the BruntVäisälä frequency which remains difficult to estimate in this region. For this reason, we consider the average of the BruntVäisälä frequency below the base of the convective zone (as prescribed by the Schwarzschild’s criterion) over a distance chosen as equal to 0.1H_{p} in the radiative zone (see Sect. 6.4 for a discussion). Such a procedure leads to the value N_{0} = 282 μH_{z}. To finish, there is still to choose three additional parameters, the plume filling factor , the penetration length L_{p} and the characteristic plume turnover frequency ν_{p}. For the filling factor, we use a reasonable value which corresponds to the values observed in numerical simulations of penetrative convection (e.g., Brummell et al. 2002) and of the uppermost layers of the convective zone (Stein & Nordlund 1998). This leads to a number of plumes that is equal to . We also take L_{p} = 0.05H_{p}, the upper limits given by Basu (1997) for the Sun from seismic observations, and ν_{p} = 1 μH_{z}, a frequency around the convective turnover frequency (see discussion in 3.4).
5.2. The case of a discontinuous transition (d = 0)
Fig. 3 Mean radial wave energy flux per unit of frequency at the top of the radiative zone as a function of the radian frequency, for any azimuthal number m and angular degrees l = 1,2 and 3, in the case of the plumeinduced waves (solid lines) and the one of turbulenceinduced waves from the formalism of Kumar et al. (1999) (dashed lines). 

Open with DEXTER 
In a first approach, we set d = 0. The result of the numerical calculation of Eq. (31)at the top of the radiative zone is plotted in Fig. 3 (solid lines). The Gaussian profile of the spectrum comes from the term ). For ω ~ 10^{6} rad s^{1} and l = 1, we read ℱ_{r} ≈ 3 × 10^{6} W m^{2} rad^{1} s. For comparison, using Eq. (39)and the value of the physical quantities given in Sect. 5.1, we find ℱ_{r} ≈ 3.2 × 10^{6} W m^{2} rad^{1} s. As a check, we calculate the total wave energy flux integrated over a large frequency range, from 0 to 25 × 10^{6} rad s^{1}, for all the (l,m) contributions, with l ∈ [1,200] and −l ≤ m ≤ l. We find a low value amounting to about 1% of the solar flux while the mean plume kinetic energy flux at the base of the convective zone, , represents about 40% of the solar flux. These results justify a posteriori the assumption of no feedback from the waves on the plumes and verify the conservation of energy.
Hereafter, we specifically comment the effects of the different terms in Eq. (31)and the dependence of the spectrum on the plume characteristics ν_{p}, b and the penetration length L_{p}.
Fig. 4 Mean radial wave energy flux per unit of frequency at the top of the radiative zone as a function of the radian frequency, for any azimuthal number m, an angular degree l = 1 and varying plume frequencies. 

Open with DEXTER 

Time correlation (Fig. 4): we consider fourvalues in range of one order of magnitude around 1μH_{z} for the plume turnover frequency, ν_{p} = 0.5,1,2.5 and 5μH_{z}. The result is plotted in Fig. 4 as a function of the cyclic frequency. As said above, an increase of ν_{p} involves a redistribution of the wave energy over higher frequencies at the expense of the lower ones since the total wave energy does not depend on the plume lifetime (under the assumption of stationarity). Its influence on the shape of the wave energy flux is significative: an increase by a factor five transforms a bellshaped spectrum into an nearly flat one in the considered frequency range.
Fig. 5 Mean radial wave energy flux at the top of the radiative zone and at the frequency ω = 10^{6} rad s^{1}, as a function of the angular degree and for varying plume radii, b. The numerical results from Eq. (33)(solid lines) are compared to the results using Eq. (37)(dashed lines). The filling factor is supposed to be constant so that a change in the plume radius results in a change of the number of plumes.
Open with DEXTER 
2.
Horizontal correlation (Fig. 5): we assume a constant filling factor, , so that a change in the plume radius involves a change in the number of plumes. We consider four values in a range of one order of magnitude around 10^{4} km for the plume radius, b = 0.5,1,2.5 and 5 × 10^{4} km, leading to a number of plumes and 40, respectively. Although the last value seems unrealistic for the Sun, we use it to study the general behavior of the wave energy spectrum with respect to the plume radius. The wave energy flux is maximal for a degree that corresponds to the horizontal extension of the plume, l_{max} = r_{b}/b, independently of the frequency. We find l_{max} = 102,51,20 and 10, respectively. that agrees with the numerical results plotted in Fig. 5 as a function of the angular degree at the fixed frequency of 10^{6} rad s^{1}. As mentioned in Sect. 4.3, we can verify on the same figure that the use of Eq. (37)(dashed lines) instead of Eq. (33)(solid lines) for the term ℬ_{l} is valid below a critical degree l_{crit} = 2πr_{b}/b ≈ 125 and 63 for b = 2.5 and 5 × 10^{4} km, respectively. Furthermore, a decrease in the plume radius induces an increase of the wave energy flux at low degrees at the expense of higher degrees. The total wave energy flux integrated over the frequency range ω = 0−25 × 10^{6} rad s^{1} and for the degrees l = 1−200 is equal to 1.3%,1%,0.5% and 0.25% of the solar flux for b = 0.5,1,2.5 and 5 × 10^{4} km, respectively. Thus, the total wave energy flux decreases with larger plume radii, particularly because the number of penetrating plumes decreases in order to keep a constant filling factor.

3.
Driving term (Fig. 6): we choose three values for the penetration length within a range of two orders of magnitude, L_{p} = 0.01,0.1 and 1H_{p}. The wave energy spectrum depends on this parameter only through the term H_{l} in Eq. (32). The result is plotted in Fig. 6 as a function of the angular degree at the fixed frequency of 10^{6} rad s^{1}. An increase of the penetration length induces a decrease of the transmission of power from the plume into the waves. The higher the angular degree, the higher the decrease. Although this effect is moderate for low degrees (decrease by a factor about two for L_{p} = H_{p} and l< 50), it becomes critical for higher degrees (decrease of about one order of magnitude for L_{p} = H_{p} and l> 150). Moreover, this effect results in a shift of the maximum of the wave energy spectrum to lower angular degrees with an increasing penetration length. It goes from l_{max} ≈ 50 for L_{p} = 0.01H_{p} to l_{max} ≈ 30 for L_{p} = 1H_{p}. The total wave energy flux integrated over the frequency range ω = 0−25 × 10^{6} rad s^{1} for l = 1−200 represents about 1.3%,0.8% and 0.2% of the solar flux for L_{p} = 0.01,0.1 and 1H_{p}, respectively. Finally, we check that the assumption made on ℋ_{l} to derive Eq. (39)(black dashed line) is valid if L_{p} ≪ r_{b}/ Λ, i.e., l ≪ r_{b}/L_{p} = 940,94 and 9 for the three considered cases.
In summary, Eq. (39)has been shown to be a good approximation of Eq. (31)provided that l ≪ l_{crit} and L_{p} ≪ r_{b}/ Λ. We have seen that ν_{p}, b and L_{p} have a significative effect on the shape of the wave energy spectrum and the total energy transferred into the waves. This could have serious consequences especially for the issue of transport by IGW. Last, the set of values used hereabove for ν_{p}, b and L_{p} will enable us to compare penetrative convection to driving by turbulent convection considering a wide domain for these parameters. This will be discussed in Sect. 6.
Fig. 6 Mean radial wave energy flux at the frequency ω = 10^{6} rad s^{1}, as a function of the angular degree l for different penetration lengths. The black dashed line represents the result obtained using Eq. (39). 

Open with DEXTER 
5.3. Effect of a thermal transition layer (d ≠ 0)
We now consider the case of a thermal transition at the radiative/convective interface (d ≠ 0 in Fig. 2). We choose a transition lentgh d = 500 km, i.e., d ≈ 0.01H_{p} = L_{p}/ 5 (see Sect. 6.4 for a discussion).
In the derivation of Eq. (36), we have considered the limiting cases of very sharp transition and of a very smooth one; the first one, corresponding to γ_{d}(ω,l) ≪ 1 in Eq. (35), allows for using a firstorder Taylor expansion of Airy functions at 0 in the transition layer, while the second one, corresponding to γ_{d}(ω,l) ≫ 1, enables us to use their asymptotic expansion at −∞ (see Appendix B.6 for details). To avoid a discontinuity between the two regimes at γ_{d}(ω,l) = 1 (as in Eq. (36), see the dashed line in Fig. 7) and to be more realistic, we assume that the switchover from the first case to the second one is realized on the range γ_{d}(ω,l) ∈ [0.80,1.20]. Instead of using Eq. (36), we define the transmission function, F, as a weighted average of the two limiting cases, (40)where with g the weighting function. The result is plotted in Fig. 7 as a function of the cyclic frequency for l = 1,2,3 and 10. Physically, it represents the ratio of the wave energy flux obtained in the case of a linear transition (d = 500 km), to the flux in the case of a discontinuous transiton (d = 0). We see that the wave energy flux can be enhanced by factor two to five, depending on the ratio k_{h}dN_{0}/ω. Indeed, the higher is this ratio, the smaller is the radial wavelength (relatively to d), and the smoother is the transition from the point of view of this component (ω,l). This result agrees with the work of Lecoanet & Quataert (2013) who took into account a smooth thermal transition in the case of wave driving by turbulent pressure.
Fig. 7 Transmission function defined in Eq. (40) for d = 0.01H_{p} (solid line), as a function of the radian frequency for l = 1,2,3 and 10. It is compared to the result obtained using Eq. (36)for l = 10 (dashed line). 

Open with DEXTER 
The total wave energy flux at the top of the radiative zone integrated over the frequency range ω = 2−25 × 10^{6} rad s^{1} for l = 1−200 represents about 1.2% of the solar flux, which is about 10 times higher than in the case of a discontinuous transition. For ω = 0−25 × 10^{6} rad s^{1} and l = 1−200, integration gives about 11% of the solar flux, i.e., about one fourth of the mean plume kinetic energy flux at the base of the convective zone. However, the latter result is likely to be overestimated because of unmodeled physical processes at very low frequencies (see Sect. 6.5 for a discussion).
Hence, a smooth thermal transition at the radiative/convective interface is able to significatively enhance the transmission of the wave energy flux into the radiative zone. It also has an impact on the shape of the wave energy spectrum at the top of the radiative zone since it does not affect in the same way the different wave components. Therefore, the length of the thermal adjustment layer, d, is a key parameter in the generation process of IGW. Unfortunately, numerical simulations of penetrative convective are not currently able to quantify the extent of this thermal adjustment layer in the solar parameters regime (e.g., Brummell et al. 2002). Nevertheless, several semianalytical models of penetrative convection can give us some hints to its value (see Sect. 6.4).
6. Discussion
In this section, we compare our results with those based on the excitation by turbulent pressure through the whole convective zone as formulated by Kumar et al. (1999). We also compare them with the studies of Press (1981) and Lecoanet & Quataert (2013) that focus on the excitation by convective eddies at the radiative/convective boundary. Lastly, we investigate the sensitivity of the present excitation model to the input physics and discuss its limits.
6.1. Comparison with turbulenceinduced wave energy flux (Kumar et al. 1999)
IGW can also be generated by turbulent pressure in the convective zone. Kumar et al. (1999) proposed a model of excitation by Reynold’s stress due to turbulent motions in the convective zone. The energy cascade is supposed to follow the Kolmogorov’s spectrum. The wave energy flux per unit of frequency at the base of the convective zone for an angular degree l and an azimuthal number m is given by (43)where ξ_{r} and ξ_{h} are the radial and horizontal wave displacements normalized to unit IGW energy flux just below the convective zone, R_{s} is the outer radius of the convective zone, v is the velocity of the biggest convective eddies with a size L and a turnover time τ_{L} = L/v, and h_{ω} = Lmin [1,2(ωτ_{L})^{− 3/2}] is the size of the convective eddies with frequency ω at radius r. The lower limit frequency for these waves is the convective turnover frequency at the base of the convective zone, ω_{c}, such as . To compute this expression, L and v are deduced from the MLT and the wave displacement is approximated by the normalized WKB wave functions as in Eq. (25) of Kumar et al. (1999). In our solar model, we find a convective turnover frequency ω_{c} ≈ 0.8 × 10^{6} rad s^{1} and a convective velocity v_{c} ≈ 25 m s at the base of the convective zone, which corresponds to a convective flux that is equal to about 5% of the solar flux. The result for the solar model is plotted in Fig. 3 (dashed lines) as a function of the cyclic frequency for l = 1,2 and 3, and in Fig. 8 as a function of the angular degree at ω = 10^{6} rad s^{1}, since we consider that most of the wave energy is distributed around this frequency in both cases.
The total turbulenceinduced wave energy flux integrated over ω = 0.8−25 × 10^{6} rad s^{1} for angular degrees l = 1−200 is equal to 0.1% of the solar flux. With 0.6% of the solar flux, the excitation by penetrative convection is 6 times more efficient in this frequency range. By integrating over the frequency range 1−25 × 10^{6} rad s^{1}, the total turbulenceinduced wave energy flux falls to 0.02% of the solar flux, mainly because most of the energy is contained at low frequencies and because the spectrum decreases rapidly around 10^{6} rad s^{1}. In turn, the plumeinduced wave energy flux remains at 0.5% of the solar flux. Indeed, the turbulenceinduced wave energy flux is about one order of magnitude smaller than the plumeinduced one in the range 2−5 × 10^{6} rad s^{1} (Fig. 3), but it dominates otherwise. In Fig. 8, we compare the wave energy flux for both driving processes as a function of the horizontal degree l. In the case of the turbulenceinduced waves, the wave energy flux is maximal for l ~ 5 which corresponds to the size of the biggest convective eddies as given by the MLT, L = α_{MLT}H_{p} ~ 10^{5} km, while the plumeinduced wave energy flux peaks at l ~ 50 because of the smaller horizontal extent of the plumes.
In Fig. 3, the plumeinduced wave energy flux drops drastically below the turbulenceinduced one for frequencies higher than 5 × 10^{6} rad s^{1}. In a similar way, in Fig. 8, the turbulenceinduced wave energy flux overtakes the plumeinduced one for degrees higher than a few characteristic widths of the plumeinduced wave spectrum, r_{b}/b. In reality, the flow of the plume should become turbulent on higher temporal (and spatial) frequencies and should generate waves by Reynold’s stress due to shear flow, as proposed by Montalbán & Schatzman (2000). Indeed, we assumed a coherent plume flow, with a spatial and a temporal characteristic extent, and neglect the turbulence inside the plume. As a consequence, the present model does not take into account plumerelated perturbations in the high frequency domain. The power transmitted from the plume into the waves on temporal (spatial) frequencies higher than few ν_{p} (few r_{b}/b) vanishes rapidly and the excitation by convective turbulence then can become predominant.
Finally, wave driving by penetrative convection is five to ten times more efficient than the excitation by turbulence as formulated by Kumar et al. (1999). While turbulence excites preferentially at very low frequencies, penetrative convection efficiently produces waves on a larger frequency range (it nevertheless depends on the value of ν_{p} relatively to ω_{c}). The plumeinduced wave spectrum peaks at higher angular degrees than the turbulenceinduced one because the plume radius is smaller than the size of the energybearing convective eddies. We emphasize that the same conclusions are reached with the wide set of values considered for ν_{p}, b and L_{p} in Sect. 5.2. All the differences observed between wave energy flux spectrum produced by both mechanisms will have consequences on the transport of angular momentum by IGW which are discussed in the following section.
Fig. 8 Mean radial wave energy flux per unit of frequency, as a function of the angular degree l at the frequency ω = 10^{6} rad s^{1}, generated by penetrative convection (in blue) and by turbulent pressure in the convective zone (in green) from the formalism of Kumar et al. (1999). 

Open with DEXTER 
6.2. Comparison with models of excitation by turbulent eddies at the radiative/convective boundary
The excitation model as proposed by Kumar et al. (1999) assumes that waves are generated by Reynold’s stress in each layer of the whole convective region and tunnel towards the radiative zone where they can propagate. While this model is used in most of the actual computation of angular momentum transport by internal gravity waves, it is also interesting to compare our results with the excitation models of Press (1981) and Lecoanet & Quataert (2013). Indeed, their approaches are similar to the one used in the present model in the sense that they focus on the driving of waves in the confined overshoot region.
Press (1981) considered the convective motions of the biggest eddies above the boundary as the driving mechanism. The eddies are characterized by a frequency ω_{c}, a size λ = 2π/k_{c} and a velocity v_{c} = ω_{c}/k_{c}. He assumed that the eddies generate waves with a frequency ω ~ ω_{c} and a horizontal wavenumber k_{h} ~ k_{c}, and he matched the turbulent pressure with the wave pressure perturbation, i.e., . Assuming a discontinuous thermal transition (d = 0 in Fig. 2) and using the continuity of the Eulerian perturbation of pressure at the boundary, the amplitude of the vertical wave displacement just below the radiative/convective interface is deduced from the relation ξ_{r} = p′/Z_{v}, where Z_{v} = ρN_{0}ω/k_{h} is the wave impedance at the top of the radiative zone. Using the Boussinesq’s approximation, the specific wave energy just below the boundary equals . The radial wave energy flux is then obtained by multiplying ℰ with the vertical group velocity given by Eq. (17)so that we obtain (44)Equation (44) have some similarities and differences with our model and its simplified expression Eq. (39):

1.
Equation (44) is equal to the product of themechanical convective flux with the Froude number (or theconvective Mach number) at the base of the convective zone,which also plays the role of a transmission factor. This is similar tothe product of the plume kinetic energy flux with F_{R,l} in Eq. (39)with the substitution v_{c} ← V_{b}.

2.
Equation (44) contains no equivalent of the exponential terms in Eq. (39)since Press (1981) assumed that waves are efficiently excited in a narrow frequency range around ω_{c}.

3.
Excitation by convective eddies occurs in each element of surface of the spherical shell at the radiative/convective boundary whereas convective penetration generates waves only in localized regions of the spherical shell. As a result, Eq. (39) differs from Eq. (44)by an additional geometrical factor, , which represents the fraction of the area occupied by the ensemble of plumes.
The excitation models proposed by Garcia Lopez & Spruit (1991) and Zahn (1997) are very similar to the one of Press (1981), except that they consider in addition a Kolmogorov’s distribution of convective eddies with an incoherent behavior.
More recently, Lecoanet & Quataert (2013) investigated also the excitation of internal gravity waves by turbulent convection. In their model, the wave amplitude is derived by considering the turbulent Reynold’s stress as the driving mechanism in a similar way to Kumar et al. (1999). However, their approach differs from the latter since they focused on the excitation of waves in the overshooting layer and took into acount the thermal transition layer at the radiative/convective interface (see Fig. 2). Moreover, they assumed that the injection length scale of the turbulent cascade at the boundary is of the same order of magnitude as the thickness of the overshooting region (i.e., smaller than the usual mixing length of the order of H_{p}) while conserving the same convective velocity as given by the MLT. In other words, the velocity of the eddies with a typical size h<H_{p} at the boundary is higher than in the case of the description of the turbulence by the MLT as used in Kumar et al. (1999). The spectral density of the wave energy flux in the case of a smooth linear thermal transition, Eq. (42) in Lecoanet & Quataert (2013), is given by (45)which we can rewrite as (46)with ω ≥ ω_{c} and k_{h} ≤ k_{h,max}(ω), where we have used v_{c} = ω_{c}/k_{c} and k_{c} ~ 1 /H_{p}, the typical convective scales as given by the MLT. We have also defined which depends on the characteristic length scales of the turbulence that is considered in the excitation zone ( in the case of the MLT). Lecoanet & Quataert (2013) showed that such considerations may enhance the wave flux by a factor two to five. The excitation model of Lecoanet & Quataert (2013) differs from the present one by several points:

4.
The points 1 and 3 mentionedhereabove are still valid in this case.

5.
Equation (46) depends on ω and k_{h} as power laws coming from the assumption of a Kolmogorov’s spectrum for the turbulence, whereas Eq. (39)depends on them as exponential laws coming from the Gaussian profile used to model the plume velocity in the penetration region.

6.
The last term in Eq. (46) is related to the effect of the thermal adjustment layer of length d on the wave transmission into the radiative zone. As a consequence, Eq. (46) is proportional to , which is consistent with our model.
Finally, the expressions derived by Press (1981) and Lecoanet & Quataert (2013) for the wave energy flux at the top of the radiative zone scale as whereas the plumeinduced wave energy flux depends on the plume velocity as . With a plume velocity one order of magnitude higher than the convective velocity in the solar case, we can conclude that the driving by penetrative plumes is more efficient than the excitation by turbulent eddies at the radiative/convective boundary as in the model proposed by Press (1981) and Lecoanet & Quataert (2013).
6.3. Efficiency of the angular momentum transport by plumeinduced IGW
IGW are able to carry angular momentum in the radiative zone and to modify the mean rotation rate. The wave flux of angular momentum is deduced from the wave energy flux through the relation (Zahn 1997) (47)In the adiabatic case, the wave luminosity of angular momentum is conserved through the whole star and no exchange is possible with the mean flow. However, as they propagate towards the center of the star, IGW are radiatively damped from the top of the radiative zone. The total wave flux of angular momentum is deduced at each depth from the one at the top of the radiative zone, modulated by a damping term (48)Press (1981) investigated the damping of IGW and found in the quasiadiabatic approximation (49)where K is the radiative diffusion coefficient and (50)is the Dopplershifted frequency in the corotating frame, with δΩ the differential rotation with respect to rotation rate at the base of the convective zone. Hence, IGW can deposit angular momentum in presence of differential rotation since prograde waves (m> 0) and retrograde waves (m< 0) are asymmetrically damped.
The transport of angular momentum follows an advectivediffusive equation (e.g., Maeder 2009). Solving numerically the problem is out of scope here but a preliminary study can give some clues on the efficiency of the transport by waves. The characteristic time scale, T_{L}, on which IGW can affect the rotation meanflow can be estimated through (e.g., Belkacem et al. 2015a) (51)where represents the divergence of the total mean radial wave flux of angular momentum (52)In order to get an order of magnitude and to simplify the calculation, we assume that the radiative zone rotates like a solid body, so that δΩ is constant in the whole radiative zone. We consider the case of a low differential rotation with respect to the base of the convective zone, δΩ = 0.15 × 10^{6} rad s^{1} (hereafter, case l), and the case of a stronger differential rotation, δΩ = 10^{6} rad s^{1} (hereafter, case s). The case l is almost similar to the initial rotation profile in Talon et al. (2002, Fig. 3). We assume also that the waves deposit all their angular momentum when they meet a critical layer (i.e., ) and that they only travel oneway towards the core of the star. We take into account only degrees from 1 to 50 since we consider it is representative of the wave spectrum in both processes of excitation. The frequency range of integration is taken between 1 and 6 × 10^{6} rad s^{1}, since most of the wave energy is concentrated in this range and since theory is uncertain below. The result is plotted in Fig. 9.
We find that transport by plumeinduced waves is able to considerably affect the internal rotation rate in the Sun in less than 0.1 Gyr in the case l and in less than 0.01 Myr in the case s. Just below the convective envelope, characteristic time scale is of the order of the year; this is the result of the rapid damping of very high degree and low frequency prograde waves; this can generate a shear layer oscillation as obtained, for instance, in Talon et al. (2002). In the radiative interior, plumeinduced waves can affect the rotation meanflow on time scales about one order of magnitude shorter than turbulenceinduced ones. Note that the core could be decelerated in less than 0.01 Gyr in the case l and in less than 1000 yrs in the case s whatever the excitation process is considered because the specific angular momentum becomes weaker and weaker in the innermost layers as it decreases as r^{2}. Moreover, the comparison between both cases l and s emphasizes the character of IGW which tends to work against strong rotation rate gradients. The stronger the differential rotation, the shorter the time scale on which IGW can affect the rotation profile.
Fig. 9 Characteristic time scale of evolution of the rotation rate in the case of a low differential rotation δΩ = 0.15 × 10^{6} rad s^{1} (dashed lines) and a stronger one δΩ = 10^{6} rad s^{1} (solid lines). We compare the results obtained with a plumeinduced wave spectrum (blue) and with a turbulenceinduced wave spectrum (red) from the formalism of Kumar et al. (1999). 

Open with DEXTER 
Therefore, we conclude that penetrative convection can generate waves able to affect the rotation profile of solarlike stars over time scales much shorter than their lifetime on the main sequence (around 10 Gyr). Moreover, transport of angular momentum by plumeinduced waves is more efficient than by turbulenceinduced ones. A more quantitative study including all mechanisms of transport, as meridional circulation and shearinduced turbulence, have to be done to confirm this result.
6.4. Sensitivity on input physics
The proposed wave excitation model depends on several physical quantities, like the plume characteristics b, ν_{p}, V_{b}, and on other parameters such as N_{0}, , L_{p} and d. We discuss here the sensitivity of the model to these quantities.
The effect of b and ν_{p} has already been investigated in Sect. 5. We have shown that the total wave energy flux does not depend on the value of ν_{p} since we assume a stationary process of excitation and that it is almost inversely proportional with b. However, they have a real impact on the shape of the spectrum. From Eq. (39), the smaller b (ν_{p}), the larger (the smaller) the width of the spectrum in the spatial (temporal) frequency domain. Moreover, the amplitude of the spectrum depends on 1 /ν_{p} and b^{2}. Multiplication by a factor 2 for ν_{p} and b leads to a multiplication by a factor 0.5 and 4, respectively. Thus, the error made on the values of ν_{p} and b is not negligible in the estimate of the wave energy spectrum. The plume velocity V_{b} at the entry of the penetration zone is another important plume characteristic which highly influences the wave energy flux. Indeed, following Eq. (39), the wave energy flux is proportional to . This means that an increase of 10%, 20% and 30% in V_{b} leads to an increase of about 50%, 250% and 500% in the wave energy flux, respectively. Similarly, a decrease of 10%, 20% and 30% for V_{b} leads to a decline of about 35%, 60% and 75%. We think the model of turbulent plumes of (Rieutord & Zahn 1995) is reliable enough to give reasonable values for V_{b} and b which will not affect too much the wave energy flux, i.e., by a factor between 0.25 and 4. Note that the same issue arises in the excitation by turbulence since the wave power spectrum depends on the convective velocity to the power 3 and on the size of the energybearing to the power 4 in Eq. (43). Uncertainties linked to their determination by the MLT have also a serious effect on the estimate of the wave energy flux in their model.
Among the free parameters, the filling factor represents the fraction of the area occupied by the ensemble of plumes at the base of the convective zone. In our model, the wave energy flux is proportional to . Conservation of matter imposes since it exists upflows less dense and slower than the plume downflows. We use a reasonable value of about and expect it not to vary too much around this value given what is observed in numerical simulations of the uppermost layers of the convective zone (Stein & Nordlund 1998). Given Sect. 6.1, a filling factor lower than 0.01 is required for penetrative convection to become negligible compared to the excitation by turbulent convection.
The penetration length, L_{p} (see Fig. 2), has been theoretically shown to be of the order of magnitude of H_{p} by Zahn (1991). Basu (1997) observationally found a lower limit valued at 0.05H_{p} in the Sun, i.e., more than one order of magnitude below the theoretical prediction. The wide range of values for L_{p} considered in Sect. 5.2 reflects this discrepancy between theory and observations. From Fig. 6, we see that a change in L_{p} does not influence too much the wave energy spectrum for low angular degrees. A multiplication by 1000 of L_{p} leads to a division by about five of the wave energy flux for l< 50. However, the drop is drastically larger for higher degrees and can reach two orders of magnitude for l> 100. In the scope of the transport of angular momentum by IGW, L_{p} should highly influence high degrees and so the SLO below the convective zone (since the wave damping τ ∝ l^{3}), while low degrees that are damped deeper in the star would not be affected.
Another difficulty is the determination of the BruntVäisälä frequency at the top of the radiative zone, N_{0}, since the radiative/convective transition is very sharp in stellar evolution codes. In this work, we arbitrarily propose to take the mean value of N over 0.1H_{p}, i.e., 1% of the thickness of the solar radiative zone, just below the Schwarzschild’s criterion. Its value influences more the wave energy spectrum in the case of sharp transition (ℱ_{r} ∝ 1 /N_{0}) than in the case of a smooth transition (), but its influence remains significative. However, uncertainties on N_{0} does not change the comparison with the turbulenceinduced wave energy spectrum since it plays the same role in both models. Indeed, it plays a role only in the transmission of the wave into the radiative zone and not in the excitation process.
The above remark about N_{0} is also valid for the length d of the thermal transition zone (r_{d} ≤ r ≤ r_{b}−L_{p}, see Fig. 2). In addition, we have shown the larger d, the better the transmission of the waves, so that the wave energy flux at the top of the radiative zone obtained in the case of a discontinuous transition (d = 0) appears as a lower limit. Several estimates of the thickness of the thermal adjustment layer have theoretically been done in the solar case. Schmitt et al. (1984) found an upper limit at about 500 km, whereas Zahn (1991) evaluated d to be of the order of the kilometer. More recently, Rempel (2004), using a semianalytical model of penetrative convection taking into account the interaction with the upflows and a distribution of the plume velocity, estimated d ≈ 350 km. Given these previous works, we have chosen d = 500 km in Sect. 5.3 as an upper limit for the transition length. With such a value, the total plumeinduced wave energy flux at the top of the radiative zone have been shown to represent about 10% of the solar flux. Hence, in the case of a linear profile for the BruntVäisälä frequency in the transition region, the plumeinduced wave energy flux goes from 1% to 10% of the solar flux for the considered range d ∈ [0,500] km. Since d plays the same role in both excitation processes, the turbulenceinduced wave energy flux is expected to be between 0.1% to 1% of the solar flux for the same range of values for d. In the case of a linear smooth transition (γ_{d} ≫ 1), we have seen that the wave energy flux is proportional to d^{1/3}, so that a decrease of one order of magnitude for d induces only a decrease by about a factor two for the wave energy flux. However, the sensitivity of the wave transmission to the parameter d depends on the profile of the transition. For instance, Lecoanet & Quataert (2013) showed that the wave energy flux is actually proportional to d in the case of a smooth hyperbolic tangent profile for the BruntVäisälä frequency at the base of the convective zone. It is thus more strongly sensitive to the value of d than for a linear transition profile. At last, these results could explain the very large wave energy flux observed in numerical simulations. For example, Dintrans et al. (2005) found that IGW could carry up to 40% of the total kinetic energy despite a very low Péclet number at the base of the convective zone (imposed by numerical constraints). As explained in Dintrans et al. (2005), a low Péclet number means a weak buoyancy braking of the plume (and so an inefficient excitation by penetrative convection), but it also means a very smooth transition (e.g., Brummell et al. 2002) which can considerably enhance the transmission of IGW into the radiative zone.
Finally, despite the wide range of values considered in Sect. 5.2 and in this section for the parameters of the excitation model, the same result holds true: excitation by penetrative convection remains more efficient than by turbulent pressure in the Sun and generates a total wave energy flux larger than 0.1% of the solar flux.
6.5. Limits of the model
It remains to discuss the limits of the model and the improvements to be done in the future. First, in Sect. 6.4, we have mentioned the dependence of the plumeinduced wave energy spectrum on the parameters and the plume characteristics by assuming that they were independent. However, the reality is much more complex. For example, the penetration length has been shown to depend on the filling factor and the plume velocity at the base of the convective zone. Schmitt et al. (1984) and, later, Zahn (1991) demonstrated that L_{p} is in fact proportional to . Besides, we can intuitively understand that the plume radius and the plume lifetime must be correlated. The larger the plume, the longer its lifetime. For instance, the characteristic destruction time of the plume by baroclinic eddies is proportional to the plume radius from the model of Jones & Marshall (1997) in Eq. (A.15).
Second, the plume characteristics are not unique and must be distributed around mean values. This is contradictory with the assumption that all the plumes are identical. For example, a distribution of plume velocities at the base of the convective zone could affect the penetration length and the steepness of the thermal transition (Rempel 2004). A distribution of the plume parameters could also have an impact on the shape of the wave energy spectrum at the top of the radiative zone. Numerical simulations of convection will be essential to answer these questions. Exhaustive studies could be done by varying parameters such as the Reynolds or the Péclet number of the simulations in order to find scaling relations and to extrapolate to stellar regimes.
Lastly, we have found that the total wave energy flux at the top of the radiative zone represents one fourth of the mean plume kinetic energy flux in the case of a smooth thermal transition (d = 500 km). However, this result has to be taken with caution since it accounts for very low frequencies at which the present description is not valid. For example, physical processes like nonadiabatic effects cannot be neglected for frequencies that are typically lower than 10^{6} rad s^{1}. In addition, the derived wave transmission coefficient diverges when ω → 0 in the case of a very smooth transition (see Fig. 7). As a consequence, the amplitude of the transmitted lowfrequency waves rises steeply and the linear approximation is no more valid as soon as k_{r}ξ_{r} ≈ 1 and nonlinear processes like wave breaking can occur. Finally, while beyond the scope of this article, we note that the Coriolis acceleration is to be taken into account in the modeling of the driving processes for wave frequencies around (or lower) than twice the internal rotation rate. This is however a challenging task since it requires a full twodimensional description of both the wave field and the plume dynamics.
7. Conclusion and perspectives
In this paper, we propose a semianalytical model of excitation of IGW by penetrative plumes. We consider the pressure exerted by an ensemble of plumes at the radiative/convective interface as the driving mechanism that is supposed to be random and stationary. The wave equation with source term is solved by taking into account the thermal transition from a convective to a radiative gradient at the base of the convective zone. We then determine the amplitude of the waves and investigate the effect of the steepness of the thermal adjustment layer on the wave transmission into the radiative zone. The description of the plumes at the base of the upper convective region follows the model of Rieutord & Zahn (1995). We thus obtain the expression of the wave energy flux per unit of frequency in the whole radiative zone in Eq. (31)as well as its approximated version in Eq. (39). Finally, the model depends on four free parameters: the filling factor, the plume lifetime, the length of penetration and the one of the thermal transition layer.
Numerical calculations for a solar model with reasonable values for the parameters show that the excitation by penetrative convection is five to ten times more efficient than the driving by turbulence in the convective zone as formulated by Kumar et al. (1999). Moreover, we demonstrate that plumeinduced waves are able to modify the internal rotation profile of the Sun in about 0.1 Gyr, which is more than ten times faster than with turbulenceinduced waves. Plumeinduced waves are thus able to considerably affect the internal rotation rate of solarlike stars on times shorter than their lifetime on the mainsequence. Furthermore, we show that these results are conservative despite a wide range of values considered for the parameters of the model.
These preliminary results on the Sun will have to be numerically quantified including transport by meridional circulation and shearinduced turbulence with plumeinduced waves, and turbulenceinduced waves. The study will have to be extended to subgiant and red giant stars with appropriate values of the free parameters. The present excitation model enables an easy implementation in 1D stellar evolution codes, and could be used to investigate the redistribution of angular momentum all along the evolution of stars taking into account transport by plumeinduced internal gravity waves as well as turbulenceinduced ones, and other mechanisms.
At last, the effect of rotation on the wave dynamics will have to be included in order to get a more realistic view of the transport by internal waves in stars. Indeed, whatever the driving process is, low frequency waves are affected by the Coriolis force. Thus, it cannot be neglected and it will have an impact on the transport by gravitoinertial waves in the radiative zone (Mathis & Zahn 2004; Mathis 2009) and on the excitation process (e.g., Mathis et al. 2014). Such a development represents a theoretical challenge for the future since we will have to go beyond the traditional approximation for subinertial waves (i.e., ω ≤ 2Ω) and to tackle the problem in 2D, as a function of the radial and latitudinal variables.
We conserve the same velocity profile than in Zahn (1991), but we consider the penetration length L_{p} as a free parameter of the model.
Indeed, the integrand in Eq. (32) is the product of a first term related to the plumeinduced vertical force per unit of mass with a second term proportional to the radial wave velocity field, solution of the homogeneous wave equation in the WKB approximation.
References
 Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Alvan, L., Strugarek, A., Brun, A. S., Mathis, S., & Garcia, R. A. 2015, A&A, 581, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andersen, B. N. 1994, Sol. Phys., 152, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Baglin, A., Auvergne, M., Barge, P., et al. 2006a, in ESA SP 1306, eds. M. Fridlund, A. Baglin, J. Lochard, & L. Conroy, 33 [Google Scholar]
 Baglin, A., Auvergne, M., Boisnard, L., et al. 2006b, in COSPAR Meeting, 36th COSPAR Scientific Assembly, 36, 3749 [Google Scholar]
 Basu, S. 1997, MNRAS, 288, 572 [NASA ADS] [Google Scholar]
 Bedding, T. R., Huber, D., Stello, D., et al. 2010, ApJ, 713, L176 [NASA ADS] [CrossRef] [Google Scholar]
 Belkacem, K., Marques, J. P., Goupil, M. J., et al. 2015a, A&A, 579, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Marques, J. P., Goupil, M. J., et al. 2015b, A&A, 579, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Brown, T. M., ChristensenDalsgaard, J., Dziembowski, W. A., et al. 1989, ApJ, 343, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825 [NASA ADS] [CrossRef] [Google Scholar]
 Brun, A. S., Miesch, M. S., & Toomre, J. 2011, ApJ, 742, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Cantiello, M., Mankovich, C., Bildsten, L., ChristensenDalsgaard, J., & Paxton, B. 2014, ApJ, 788, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Ceillier, T., Eggenberger, P., García, R. A., & Mathis, S. 2013, A&A, 555, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Charbonnel, C., & Talon, S. 2008, in IAU Symp. 252, eds. L. Deng, & K. L. Chan, 163 [Google Scholar]
 Charbonnel, C., Decressin, T., Amard, L., Palacios, A., & Talon, S. 2013, A&A, 554, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 ChristensenDalsgaard, J., Monteiro, M. J. P. F. G., Rempel, M., & Thompson, M. J. 2011, MNRAS, 414, 1158 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dintrans, B., Brandenburg, A., Nordlund, Å., & Stein, R. F. 2005, A&A, 438, 365 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fuller, J., Lecoanet, D., Cantiello, M., & Brown, B. 2014, ApJ, 796, 17 [NASA ADS] [CrossRef] [Google Scholar]
 García, R. A., TurckChièze, S., JiménezReyes, S. J., et al. 2007, Science, 316, 1591 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Garcia Lopez, R. J., & Spruit, H. C. 1991, ApJ, 377, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., Murray, N., & Kumar, P. 1994, ApJ, 424, 466 [NASA ADS] [CrossRef] [Google Scholar]
 Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1986, ApJ, 311, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, H., & Marshall, J. 1997, J. Phys. Oceanogr., 27, 2276 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, P., Talon, S., & Zahn, J.P. 1999, ApJ, 520, 859 [NASA ADS] [CrossRef] [Google Scholar]
 Lebreton, Y., Goupil, M. J., & Montalbán, J. 2014, in EAS Pub. Ser., 65, 99 [Google Scholar]
 Lecoanet, D., & Quataert, E. 2013, MNRAS, 430, 2363 [NASA ADS] [CrossRef] [Google Scholar]
 Lighthill, J. 1978, Waves in fluids (CUP) [Google Scholar]
 Maeder, A. 2009, Physics, Formation and Evolution of Rotating Stars (Berlin, Heidelberg: Springer) [Google Scholar]
 Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S. 2009, A&A, 506, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Zahn, J.P. 2004, A&A, 425, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Decressin, T., Eggenberger, P., & Charbonnel, C. 2013, A&A, 558, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Montalbán, J., & Schatzman, E. 2000, A&A, 354, 943 [NASA ADS] [Google Scholar]
 Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012a, A&A, 548, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012b, A&A, 540, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Press, W. H. 1981, ApJ, 245, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Rempel, M. 2004, ApJ, 607, 1046 [NASA ADS] [CrossRef] [Google Scholar]
 Rieutord, M., & Zahn, J.P. 1995, A&A, 296, 127 [Google Scholar]
 Rogers, T. M., & Glatzmaier, G. A. 2006, ApJ, 653, 756 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G., Gellert, M., Spada, F., & Tereshin, I. 2015, A&A, 573, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schmitt, J. H. M. M., Rosner, R., & Bohn, H. U. 1984, ApJ, 282, 316 [NASA ADS] [CrossRef] [Google Scholar]
 Spruit, H. C. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stein, R. F., & Nordlund, A. 1998, ApJ, 499, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Stull, R. B. 1976, J. Atmos. Sci., 33, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Talon, S., & Charbonnel, C. 1998a, in Cool Stars, Stellar Systems, and the Sun, eds. R. A. Donahue, & J. A. Bookbinder, ASP Conf. Ser., 154, 973 [Google Scholar]
 Talon, S., & Charbonnel, C. 1998b, A&A, 335, 959 [NASA ADS] [Google Scholar]
 Talon, S., & Charbonnel, C. 2003, A&A, 405, 1025 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Charbonnel, C. 2004, A&A, 418, 1051 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., & Charbonnel, C. 2005, A&A, 440, 981 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Talon, S., Kumar, P., & Zahn, J.P. 2002, ApJ, 574, L175 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, A. A. 1966, J. Fluid Mech., 24, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, J. S. 1986, J. Fluid Mech., 173, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial oscillations of stars, 2nd edn. (Tokyo: University of Tokyo Press) [Google Scholar]
 Zahn, J.P. 1991, A&A, 252, 179 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1997, in SCORe’96: Solar Convection and Oscillations and their Relationship, eds. F. P. Pijpers, J. ChristensenDalsgaard, & C. S. Rosenthal, Astrophys. Space Sci. Libr., 225, 187 [Google Scholar]
Appendix A: Plume lifetime and restratification in the Sun
We adapt the model developed by Jones & Marshall (1997) in the penetration zone of the Sun where a lateral density contrast exists between the plume and the surroundings. We consider a fluid cylinder of height and radius equal to L_{p} and b, respectively (see Fig. A.1)
Fig. A.1 Schematic view of the plume (in blue) in the penetration region, surrounded by a stable stratitfied medium. μ represents the horizontal length scale of the density contrast between the two distinct parts. 

Open with DEXTER 
Appendix A.1: Geostrophic wind in the solar penetration zone
Jones & Marshall (1997) investigated the restratification of convective patches in the oceans for which the geostrophic approximation is valid. Under these conditions, the flows follow the geostrophic wind equations resulting in a balance between the pressure gradient and the Coriolis force. In order to adopt their model to the Sun, we must first verify the validity of the geostrophic approximation in the solar penetration zone.
 1.
Acceleration in the corotating frame with the base of the convection zone must be negligible compared to the Coriolis force, so that the Rossby number (A.1)is smaller than unity, with Ω the mean angular rotation rate, L the characteristic horizontal length scale and U the zonal velocity in the corotating frame. In the solar tachocline case, taking Ω ≈ 2.5 × 10^{6} rad s^{1}, L ≈ 2b ≈ 2 × 10^{4} km, we will verify a posteriori that R_{o} ≪ 1 provided that U ≪ 100 m s.
 2.
Viscous stress can be neglected compared to the Coriolis term in the equation of motion, the Ekman number (A.2)with ν the viscosity, is much smaller than unity. Even with a turbulent viscosity due to local shear expected to be around , the Ekman number stays very small.
 3.
Finally, the vertical length scale must be negligible compared to horizontal one. In the solar case, taking L_{p} ≈ 0.1H_{p}, H_{p} ≈ 5 × 10^{4} km, the pressure scale height at the base of the convective zone, and b ≈ 10^{4} km, the ratio is equal to (A.3)
Appendix A.2: Plume destruction time scale by geostrophic eddies (Jones & Marshall 1997)
By assuming vertical hydrostatic balance (which is possible inside the plume if we consider an almost uniform and stationary flow), we can then use the thermal wind equation (A.4)where f = 2sinθΩ is the Coriolis parameter, ρ the density, p the pressure and k the unit vector perpendicular to the surface of the sphere. The density contrast δρ between the plume and the surrounding medium creates a vertical shear, which obeys to (A.5)where μ is the horizontal length scale on which density exchanges take place and δb = −gδρ/ρ represents buoyancy acceleration. To go further, we use a first order expansion at z ~ L_{p} of δb, (A.6)where we impose (δb)_{Lp} = 0 and we note (A.7)with the mean BruntVäisälä frequency in the penetration zone. By injecting Eq. (A.6)in Eq. (A.5), we obtain the zonal wind as a function of the altitude (A.8)where the constant of integration is chosen to satisfy . Jones & Marshall (1997) showed that the thermal wind undergoes a baroclinic instability when the radius of the fluid column is longer than the Rossby length scale of deformation (A.9)In the penetration zone, injection of convective material imposes a quasiadiabatic gradient before reaching a radiative regime where the BruntVäisälä frequency reaches N_{0} ≈ 300μH_{z}. Hence, we can estimate N_{t} ~ 1−10μH_{z}. With L_{p} ~ 0.1H_{p}, we find L_{D} ≈ 10^{2}−10^{3} km ≤ b and the fluid column in the Sun undergoes a baroclinic instability. Moreover, at z = 0 , we find (A.10)We then assume that the density exchanges take place on a more expanded length scale than the Rossby radius of deformation, so μ ≥ L_{D}. Therefore, u(z = 0) ≤ 10 m s^{1}, the Rossby number R_{o} ≪ 1 and we verify a posteriori the thermal wind approximation (cf. paragraph on Rossby number Eq. (A.1)).
By conservation of potential buoyancy, is conserved in the cylinder of fluid, where the bar denotes time average over a characteristic time scale. Let v be the eddiesinduced velocity field. The Reynolds transport theorem enables us to write (A.11)where n is the unit vector perpendicular to the lateral surface and where we have neglected the contribution of the buoyancy fluxes through the top and the bottom of the cylinder. To go further, we suppose (A.12)with c_{e} representing the efficiency coefficient of the transport by geostrophic eddies. The righthand side of Eq. (A.12)becomes after integration (A.13)The lefthand side of Eq. (A.12)gives (A.14)whence we can deduce the restratification time scale (A.15)From numerical experiments and observations, Jones & Marshall (1997) derived c_{e} = 0.027 in the case of the oceans. Given the lack of knowledge about this process in the Sun, we use this value to estimate the plume lifetime in our case. We find t_{res} ~ 10^{6}−10^{7} s with μ ~ L_{D}, i.e ν_{p} = 1 /τ_{p} ~ 0.1−1μH_{z}. Note that this value is of the same order of magnitude than the turnover convective frequency, ω_{c} = v_{c}/l_{c}, as given by the MLT in a standard model, with v_{c} the convective velocity and l_{c} the mixing length.
Appendix B: Derivation of the wave energy flux
In the following, we detail the derivation of the wave power spectrum in the radiative zone. We will focus on the pressure due to convective plumes as driving mechanism in the penetration zone, and we will take into account the effect of the thermal transition on the transmission of the generated waves, as illustrated in Fig. 2.
Appendix B.1: Adiabatic wave equation with source term
As explained in Sect. 2.1, the total velocity field is split into a wave component v and a plume component V_{p}. We neglect the nonlinear wave term ∇·(ρv ⊗ v) and the plumerelated coupling terms, like ∇·(ρV_{p} ⊗ v). We also assume that the time variation of the plume impulsion in the penetration zone, ∂_{t}(ρV_{p}), does not participate into the excitation of the wave, and that the plume destruction is mainly due to instabilities and turbulent motions in this region. The equations of dynamics are then given by Eq. (1)and Eq. (2).
To go further, we adapt the procedure used by Unno et al. (1989) in order to include the forcing term and the continuous frequency spectrum of the plumeinduced wave packet. We introduce the Lagrangian wave displacement vector ξ(r,t) related to the wave velocity field by v = ∂_{t}ξ. The Eulerian perturbations are decomposed onto the spherical harmonics, like the wave velocity field in Eq. (13)or, for example, the radial Lagrangian displacement (B.1)Deriving the spectral density of the wave specific energy Eq. (12)requires the computation of the time Fourier transform of the wave velocity field. By taking the time Fourier transform of Eqs. (1) and (2) and using the fact that spatial differential operators commute with time Fourier transform, we obtain The horizontal part of Eq. (B.2)enables to link the horizontal displacement vector to and the horizontal part of the forcing term (B.4)with ∇_{⊥} the horizontal nabla operator. To continue, we proceed in the same way than Unno et al. (1989, Eq. (13.40)). We compute the horizontal divergence of Eq. (B.2)and replace from Eq. (B.3)to find (B.5)Using the adiabatic closure hypothesis (in Fourier space), (B.6)and the decomposition onto spherical harmonics Eq. (B.1), Eq. (B.5)can be rewritten (B.7)Similarly, using and Eq. (B.6), the radial part of Eq. (B.2)becomes (B.8)with c the sound speed, S_{l} the Lamb frequency and N the BruntVäisälä frequency. After projecting onto , we obtain a firstorder linear system of two differential equations for and
The differential system is similar to Eq. (14.2) and (14.3) in Unno et al. (1989), with the additional righthand side source terms F_{1} and F_{2}, given by with the solid angle given by dΩ = sinθdθdφ.
Appendix B.2: Modeling of the source term in the wave equation
Considering the velocity profile in Eq. (19), the pressure gradient due to a single plume is given by (B.13)We assume that the radius at the interface, r_{b}, is much larger than the penetration length, L_{p} ≪ r_{b}. Hence, the operator 2 /r can be neglected compared to ∂_{r}, and the radius is taken constant in this region. Therefore, using Eq. (B.13), Eq. (B.12)can be written as (B.14)where we have defined with s_{b} = S_{h}(r_{b},θ,φ;θ_{0},φ_{0}) following Eq. (20). We also find that Eq. (B.11)cancels, F_{1} = 0, since Eq. (B.13) is collinear with the radial direction. For the same reason, using Eq. (B.13)in Eq. (B.4), we see that the toroidal part of the displacement vector is null, ξ_{T,l,m} = 0, and that the poloidal part, renamed , is linked to the perturbation of pressure through the expression (B.17)Finally, Eqs. (B.9) and (B.10) can be rewritten in the form of the following nonhomogeneous firstorder differential system (B.18)where z, b and A are given by Eq. (29)and Eq. (30).
Appendix B.3: Homogeneous equation
To solve the homogeneous differential system, i.e., Eq. (B.18)with b = 0, we use the change of variables (Unno et al. 1989) This leads to a more tractable secondorder differential linear equation for v_{l,m}(B.21)with k_{r} given by Eq. (18)and where we have neglected the term f(P) in Eq. (16.11) of Unno et al. (1989). Therefore, is obtained through Eq. (B.19)by solving Eq. (B.21), and is deduced from Eq. (B.17), Eq. (B.20)and the relation (B.22)which comes from Eq. (B.7). This approach will be used in the following paragraphs to find the two linearly independent solutions (subscript 1 or 2) of the homogeneous differential system for each region: the quasiadiabatic region, composed of the penetration zone and the convective zone, the transition region and the radiative interior (see Fig. 2).
B.3.1 Quasiadiabatic region (r_{b}−L_{p} ≤ r)
In the convective zone and the penetration zone, we assume that the temperature gradient is quasiadiabatic, N^{2} ~ 0 and the waves are evanescent. In this region, the WKB solutions of Eq. (B.21)are a good approximation. The wave functions, and , in the quasiadiabatic region (superscript a) are then estimated up to a constant by where the subscripts r and h refer to the radial and horizontal components of the vectors, r_{b} is the radius of the base of the convective zone (as prescribed by the Schwarzschild’s criterion), L_{p} is the penetration length and .
B.3.2 Radiative interior (r ≤ r_{d})
In the same way, in the radiative zone (superscript r), propagative waves are estimated with the WKB solutions of Eq. (B.21). The inward and outward wave functions, and , are respectively given up to a constant by
where r_{d} represents the radius at the top of the radiative zone r_{d} = r_{b}−L_{p}−d.
B.3.3 Thermal transition layer (r_{d} ≤ r ≤ r_{b}−L_{p})
Once the plume is slowed down enough, radiative diffusion can operates efficiently. In this socalled thermal adjustment layer, the temperature gradient undergoes a transition from a quasiadiabatic gradient to a radiative one. We follow here the work of Lecoanet & Quataert (2013) and assume that the BruntVäisälä frequency varies linearly in the layer (Fig. 2) (B.27)with z′ = r−(r_{b}−L_{p}) + d/ 2. By injecting Eq. (B.27)in the expression of the radial wavenumber Eq. (18), Eq. (B.21)becomes: (B.28)As in Lecoanet & Quataert (2013), we define with z_{t} the wave turning point (where N^{2} = ω^{2}) and k_{h} ≈ Λ /r_{b} the horizontal wavelength, supposed to be constant in the transition layer. Then, using the change of variable χ = K^{1/3}(z′−z_{t}), we find that Eq. (B.28)verifies the equation of Airy: (B.31)whose the solution is a linear combination of the Airy functions A_{i}(χ) and B_{i}(χ). To determine the horizontal part , we use Eq. (B.22)that can be written as a function of the variable χ(B.32)Therefore, the two linearly independent^{7} solutions of the homogeneous wave equation, and , in the transition layer (superscript t) read up to a constant
Appendix B.4: Particular solution of the nonhomogeneous system
We suppose that wave driving occurs exclusively in the penetration zone where buoyancy braking is the strongest (i.e., b ≠ 0 for r_{b}−L_{p}<r<r_{b}). We derive a particular solution of the nonhomogeneous system thanks to the method of variation of parameters. In the quasiadiabatic region (r_{b} ≤ r), we search a solution in the form of (B.35)By injecting Eq. (B.35)in Eq. (B.18), the functions μ and λ must verify the linear system (B.36)whose the solution is obtained by using the Cramer’s rule where b = iF_{2}/rω. From Eqs. (B.14), (B.23) and (B.24), integration leads to where α(ω) and are given by Eq. (B.15)and Eq. (B.16), and where we have arbitrarily chosen μ(r_{b}) = λ(r_{b}) = 0.
Appendix B.5: General solution and boundary conditions
The general wave function is obtained by ensuring the continuity of the radial and the horizontal displacements at the limits between each region. We have also to consider two boundary conditions which are taken at surface and at the center of the star. In the quasiadiabatic region (r_{b}−L_{p} ≤ r), the general solution z reads (B.41)We impose A_{1} = 0, since we assume that there is no evanescent inward component coming from the surface of the star. In the transition layer (r_{d} ≤ r ≤ r_{b}−L_{p}), the general solution reads (B.42)and in the radiative interior (r ≤ r_{d}), (B.43)We choose B_{2} = 0 because we only consider the inward component of the wave (travelling towards the center). In other words, we suppose that the waves will be damped before being reflected in the core, so that the emergence of modes is not allowed.
Appendix B.6: Continuity of the wave function
To derive the wave energy flux at the top of the radiative zone (r = r_{d}), we need to determine the modulus of the constant B_{1}. By imposing the continuity of the radial and the horizontal displacement (which is equivalent to continuity of the perturbation of pressure) at the points r_{p} = r_{b}−L_{p} and r_{d}, we get four linear equations linking B_{1} to A_{2}, C_{1} and C_{2}, so that we close the system.
At the point r = r_{p}, corresponding to z′ = d/ 2 in Eq. (B.27), the argument of the Airy functions is equal to (B.44)since k_{h}d ≤ 1 and ω/N_{0} ≪ 1. Therefore, we can use the firstorder Taylor expansion at 0 of the functions A_{i}(χ) et B_{i}(χ) and their derivatives^{8}where and are linked to C_{1} and C_{2} following the Taylor expansion of the Airy functions The continuity at the radius r_{p} between the adiabatic region to the transition layer yields to (B.49)where μ_{p} = μ(r_{p}) and λ_{p} = λ(r_{p}). The second line of the system Eq. (B.49)gives (B.50)Using Eq. (B.50)in the first line of Eq. (B.49), we obtain (B.51)At the point r = r_{d}(z′ = −d/ 2), the argument of the Airy functions is equal to (B.52)and the continuity condition between the top of the radiative zone and the transition layer at r_{d} yields (B.53)
B.6.1 Case of a sharp transition
In the case of a sharp transition, i.e  χ_{d}  ≪ 1, we can use again the first order Taylor expansion of the Airy functions. The second line of the system Eq. (B.53)gives directly as a function of B_{1}(B.54)Then, using Eqs. (B.54) and (B.50), we can deduce A_{2}(B.55)Now, Eq. (B.51)enables us to determine as a function of B_{1} by replacing A_{2} by Eq. (B.55)(B.56)To go further, and are replaced by Eqs. (B.56) and (B.54) in the first line of the system Eq. (B.53)to obtain B_{1}(B.57)In the case of a sharp transition, we have k_{h}d ≪ 1, and the modulus square of B_{1} is finally given by (B.58)
B.6.2 Case of a smooth transition
We focus now on the case of a smooth transition, i.e.,  χ_{d}  ≫ 1, which is equivalent to k_{r}d ≫ 1. We use the asymptotic expansion of the Airy functions and the one of their derivatives for x → + ∞Using Eq. (B.59), Eq. (B.60)and the relation , Eq. (B.53)can be rewritten as After some algebra, we find By computing the modulus of Eqs. (B.63) and (B.64), and dividing Eqs. (B.63) by (B.64), we obtain To go further, we inject Eqs. (B.51) and (B.50) into Eqs. (B.47) and (B.48), and we calculate^{9}Using Eq. (B.67), we get A_{2} as a function of C_{1} that we reinject into Eq. (B.68)to find (B.69)Finally, we can assume (B.70)and we compute the absolute square of Eq. (B.69)to obtain, using Eq. (B.65), (B.71)Thus, we find a similar expression for  B_{1}  ^{2} to what we obtained in the case of a sharp transition. They differ from each other according to a structural factor, corresponding to the ratio of the transmission coefficients in both cases, (B.72)
where D is a numerical factor (B.73)
Appendix B.7: Final expression for the wave energy flux
In the radiative zone, the wave function is equal to . Then, using Eq. (B.25), Eq. (15)becomes Given Eqs. (16) and (17), the mean radial wave energy flux reads (B.76)In the case of a sharp transition, using Eqs. (B.58), (B.39), (B.15) and (B.16), (B.76) becomes (B.77)where we have defined In the case of a smoother transition, the derivation is similar and simple if we use Eq. (B.72), so that we obtain (B.80)
All Figures
Fig. 1 Schematic view of the star. The convective plumes occur in the upper layers of the star and go deeper into the convective region. They grow by turbulent entrainment of matter at their edges and reach the top of the radiative zone. There, each of them, characterized by their angular position (θ_{i},φ_{i}), releases a part of their kinetic energy and generate internal waves which can propagate towards the center. In our work, we suppose that there are no reflected waves propagating from the center to the radiative/convective interface. 

Open with DEXTER  
In the text 
Fig. 2 Modeling of the interface between the radiative zone an the convective zone. r_{b} corresponds to the radius at the base of the convective envelope as given by the Schwarzschild’s criterion, L_{p} is the penetration length over which the plume is decelerated and generates waves, d denotes the extent of the transition region (where the plume velocity is negligible compared to V_{b}, its value at the entry of the penetration region) and r_{d} is the radius at the top of the radiative zone. 

Open with DEXTER  
In the text 
Fig. 3 Mean radial wave energy flux per unit of frequency at the top of the radiative zone as a function of the radian frequency, for any azimuthal number m and angular degrees l = 1,2 and 3, in the case of the plumeinduced waves (solid lines) and the one of turbulenceinduced waves from the formalism of Kumar et al. (1999) (dashed lines). 

Open with DEXTER  
In the text 
Fig. 4 Mean radial wave energy flux per unit of frequency at the top of the radiative zone as a function of the radian frequency, for any azimuthal number m, an angular degree l = 1 and varying plume frequencies. 

Open with DEXTER  
In the text 
Fig. 5 Mean radial wave energy flux at the top of the radiative zone and at the frequency ω = 10^{6} rad s^{1}, as a function of the angular degree and for varying plume radii, b. The numerical results from Eq. (33)(solid lines) are compared to the results using Eq. (37)(dashed lines). The filling factor is supposed to be constant so that a change in the plume radius results in a change of the number of plumes. 

Open with DEXTER  
In the text 
Fig. 6 Mean radial wave energy flux at the frequency ω = 10^{6} rad s^{1}, as a function of the angular degree l for different penetration lengths. The black dashed line represents the result obtained using Eq. (39). 

Open with DEXTER  
In the text 
Fig. 7 Transmission function defined in Eq. (40) for d = 0.01H_{p} (solid line), as a function of the radian frequency for l = 1,2,3 and 10. It is compared to the result obtained using Eq. (36)for l = 10 (dashed line). 

Open with DEXTER  
In the text 
Fig. 8 Mean radial wave energy flux per unit of frequency, as a function of the angular degree l at the frequency ω = 10^{6} rad s^{1}, generated by penetrative convection (in blue) and by turbulent pressure in the convective zone (in green) from the formalism of Kumar et al. (1999). 

Open with DEXTER  
In the text 
Fig. 9 Characteristic time scale of evolution of the rotation rate in the case of a low differential rotation δΩ = 0.15 × 10^{6} rad s^{1} (dashed lines) and a stronger one δΩ = 10^{6} rad s^{1} (solid lines). We compare the results obtained with a plumeinduced wave spectrum (blue) and with a turbulenceinduced wave spectrum (red) from the formalism of Kumar et al. (1999). 

Open with DEXTER  
In the text 
Fig. A.1 Schematic view of the plume (in blue) in the penetration region, surrounded by a stable stratitfied medium. μ represents the horizontal length scale of the density contrast between the two distinct parts. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.