© ESO 2021

1. Introduction

In recent decades the study of global acoustic modes has revealed precious information on the internal properties of the Sun (see, e.g., Thompson et al. 2003; Basu & Antia 2008; Kosovichev 2011; Buldgen 2019; Christensen-Dalsgaard 2021, for interesting reviews). Among the noteworthy results, the sound speed and rotation rate could be probed with a high precision down to 0.08 R_⊙ and 0.2 R_⊙, respectively. Nevertheless, the study of acoustic modes hardly gives access to the properties of deeper layers, where most of the solar luminosity is produced. In contrast, gravity modes, which propagate in the central radiative zone, have the potential to probe the micro- and macrophysics inside the solar core and test further stellar models (e.g., Appourchaux et al. 2010). For instance, the detection of gravity modes, combined with neutrinos flux measurements, can be expected to improve our representation of nuclear reaction rates and electron screening (e.g., Mussack & Däppen 2011; Lopes & Turck-Chièze 2014; Bonventre & Orebi Gann 2018). The measurement of their rotational splittings can also give us the possibility to reveal the deep rotation profile and put stringent constraints on the angular momentum history of the Sun (e.g., Eggenberger et al. 2019a). From a wider perspective, all the information brought by solar gravity modes can permit, on the one hand, to calibrate stellar models and help improve our understanding of the whole stellar evolution. On the other hand, potential discrepancies between theory and observations in the earlier or later evolutionary phases can also give evidence of a change of regimes in the internal processes at work throughout the star lifetime (e.g., Gehan et al. 2018; Eggenberger et al. 2019b).

The quest for solar gravity modes began more than 40 years ago. However, gravity modes are evanescent in the convective envelope and have very small amplitudes at the solar surface; they are thus very difficult to observe. Detections have been claimed by several studies (e.g., Brookes et al. 1976; Severnyi et al. 1976; Delache & Scherrer 1983; Thomson et al. 1995; Turck-Chièze et al. 2004; García et al. 2007), but none of these results has been independently confirmed (e.g., Appourchaux & Pallé 2013). The most recent claim of their detection was made by Fossat et al. (2017) and Fossat & Schmider (2018). Using data of the GOLF instrument (e.g., Gabriel et al. 1995), the authors indirectly found evidence of the signature of low-frequency gravity modes in the time variations of the large frequency separation of high-frequency acoustic modes. Their analysis led to the identification of hundreds of gravity modes with angular degrees between ℓ = 1 and ℓ = 4, from which they could infer the asymptotic period-spacing and the mean core rotation rate. Although this work was subsequently reproduced by other studies, the robustness of the analysis was put into question (Schunker et al. 2018; Appourchaux & Corbard 2019). The result was shown to be very sensitive to the parameters of the time series considered in the analysis (e.g., the start time), and appears to be an artifact of the methodology. Theoretical studies of the coupling between acoustic and gravity modes reinforced the fragility of this result (Scherrer & Gough 2019; Böning et al. 2019). At the time of writing no robust detection of the solar gravity modes has therefore been unequivocally confirmed.

Within this framework, theoretical estimates of the amplitude of solar gravity modes are needful to help design future observational missions and guide future seismic studies. The amplitude of gravity modes results from a balance between the driving by convective motions and damping processes (e.g., Belkacem et al. 2011). Previous theoretical estimates mainly considered the Reynolds stress of turbulent convective eddies as the source mechanism, or stochastic excitation (e.g., Gough 1985; Kumar et al. 1996; Belkacem et al. 2009). Guided by global 3D numerical simulations of the solar convective zone, and using reasonable values for their model parameters, Belkacem et al. (2009) determined that the amplitudes of asymptotic gravity modes, that is, with high radial orders and oscillation frequencies between 10 µHz ≲ ν ≲  100 µHz, are likely to lie slightly below the current GOLF detection threshold. Because of the non-detection of gravity modes, these predictions set an upper constraint on the Sun’s convective velocity in the excitation region (around 0.8 R_⊙). However, all these previous estimates did not account for the contribution from the penetration of convective plumes at the base of the convective region (or penetrative convection) to the mode driving.

Convective plumes are strong downdrafts originating from diving cool granules at the solar surface. They develop by turbulent entrainment of matter as coherent structures when crossing the convective region (e.g., Turner 1986; Rieutord & Zahn 1995). As the plumes reach the bottom of the convective bulk, they can penetrate into the underlying stably stratified radiative layers; there the plumes are braked by buoyancy and can transfer a part of their kinetic energy into gravity waves. While this excitation mechanism is ubiquitous in numerical simulations of extended convective envelopes overlying radiative zones (e.g., Andersen 1996; Dintrans et al. 2005; Kiraga et al. 2005; Rogers et al. 2006, 2013; Alvan et al. 2014; Edelmann et al. 2019), the covered values of the dimensionless control parameters are far from stellar regimes. Quantitative estimates by means of semi-analytical excitation models are thus required and complementary (e.g., Rempel 2004). Motivated by the issue of the redistribution of angular momentum in stellar interiors, Pinçon et al. (2016) modeled this driving process to predict the amplitude of very low-frequency progressive gravity waves propagating in the radiative zone of the Sun. They demonstrated that this process generates low-frequency gravity waves more efficiently than turbulent pressure and can have an important impact on the angular momentum evolution of low-mass stars (Pinçon et al. 2017). However, no application regarding the amplitude of the solar gravity modes with much higher frequencies has been undertaken so far.

In this work we estimate the amplitude of the solar gravity modes generated by penetrative convection following the model of Pinçon et al. (2016). We limit the study to asymptotic gravity modes in the frequency range between 10 µHz and 100 µHz (i.e., with a number of radial nodes much higher than unity). The damping of such modes is dominated by radiative losses, and is analytically tractable in the considered frequency range under the quasi-adiabatic limit (e.g., Dziembowski 1977b; Belkacem et al. 2009). In contrast, at higher frequencies the computation of the mode damping requires accounting for the interaction between oscillations and convection, which is much more complex and beyond the scope of this paper (e.g., Belkacem et al. 2011).

The paper is organized as follows. In Sect. 2, an analytical expression for the mode energy is derived based on the model of Pinçon et al. (2016). In Sect. 3 the apparent surface radial velocity of gravity modes is estimated from this expression for a solar model and is compared to the GOLF data detection threshold. The results are discussed in Sect. 4. The conclusions are presented in Sect. 5.

2. Excitation model by penetrative convection

In this section we derive an analytical expression for the energy of gravity modes excited by penetrative convection. The asymptotic (i.e., short-wavelength) approximation is used and the quasi-adiabatic limit is considered (i.e., non-adiabatic effects are globally considered as small perturbations). Both approximations are justified for the Sun in the considered frequency range (between 10 µHz and 100 µHz). The detailed derivation steps and technical issues are described in Appendix A.

2.1. Oscillation equations forced by penetrative plumes

Following Pinçon et al. (2016), the velocity field in the stellar frame is decomposed into a component associated with the convective plumes and a perturbation associated with the gravity modes. The source term in the linearized momentum equation is assumed to be the ram pressure exerted by the ensemble of convective plumes at the top of the radiative region. The feedback from the oscillations on the plume structure and dynamics is neglected. In other words, we assume that the wave energy is much lower than the plume kinetic energy at the base of the convection zone, which is checked a posteriori. The effect of the Coriolis force on both the oscillations and the plumes is not considered. This is justified for the gravity modes as the solar rotation period is much shorter than the modal periods. For the plumes, the effect of the buoyancy work is predominant so that the plume Rossby number is expected to be very low in the case of slow rotators as the Sun.

Within this framework, the forced linear non-adiabatic oscillation equation reads (see Appendix A.1 for details)

$$\begin{array}{c}\hfill {\partial }_t^2\mathit{\xi }+{\mathcal{L}}^{\mathrm{ad}}\left(\mathit{\xi }\right)+{\mathcal{L}}^{\mathrm{nad}}\left(H\frac{\delta S}{c_p}\right)=-\frac{1}{\rho }\mathbf{\nabla }·(\rho {\mathcal{V}}_{\mathrm{p}}\otimes {\mathcal{V}}_{\mathrm{p}}),\end{array}$$(1)

where ξ(r, t) is the mode displacement field, H is the temperature scale height, δS is the Lagrangian perturbation of specific entropy, c_p is the specific heat capacity at constant pressure, ρ is the equilibrium density, 𝒱_p(r, t) is the velocity field associated with the ensemble of plumes, ℒ^ad is the adiabatic differential operator provided in Eq. (A.7), ℒ^nad is a differential operator given in Eq. (A.6) and resulting from non-adiabatic effects, ∂_t denotes the partial time derivative, ∇ is the gradient operator, and (⊗) is the outer product.

Monitoring the evolution of δS requires us to consider the perturbed heat equation. Accounting only for radiative losses in the case of asymptotic gravity modes, the perturbed heat equation can be expressed within the diffusion approximation as (see Appendix A.2.1)

$$\begin{array}{c}\hfill {\partial }_t\left(\frac{\delta S}{c_p}\right)=\frac{1}{t_{\mathrm{R}}}[{\mathcal{L}}^{\mathrm{nad}1}\left(\frac{\mathit{\xi }}{H}\right)+{\mathcal{L}}^{\mathrm{nad}2}\left(\frac{\delta S}{c_p}\right)],\end{array}$$(2)

where t_R = ρc_pTH/F_R is the local radiative thermal timescale, which describes the exchange rate of energy between the modes and the radiation, with F_R and T the equilibrium radiative flux and temperature, respectively. We note that ℒ^nad1 and ℒ^nad2 in this equation are dimensionless linear differential operators with respect to radius and represent the perturbation of the divergence of the radiative flux. We also note that t_R is equal to the local thermal timescale in the radiative zone and in a very thin near-surface layer where the energy flux is mostly carried by radiation. This is not the case in deeper convective layers where the energy flux is mostly carried by the convective flux, but whose influence on the considered gravity modes can be neglected according to the numerical computations of Belkacem et al. (2009).

2.2. Mode amplitude in the quasi-adiabatic limit

2.2.1. Decomposition on the adiabatic eigenfunction basis

The eigenfunctions ξ_nℓm of the ℒ^ad operator, with radial orders n, angular degrees ℓ, and azimuthal numbers m, were shown to form a complete basis of the oscillation displacement field under the zero-boundary conditions (e.g., Chandrasekhar 1964; Unno et al. 1989). More precisely, the demonstration relied on the assumption that the oscillations are adiabatic (i.e., δS = 0). For our purpose we show in Appendix B that they also form a complete basis of the oscillation displacement field in the non-adiabatic case (i.e., δS ≠ 0). It is therefore possible to project the non-adiabatic mode displacement field onto the basis of the eigenfunctions of the ℒ^ad operator. This gives

$$\begin{array}{c}\hfill \mathit{\xi }(\mathit{r},t)={\displaystyle \sum _{n=-\infty }^{+\infty }\sum _{\ell =0}^{+\infty }\sum _{m=-\ell }^{+\ell }{a}_{n\ell m}(t){\mathit{\xi }}_{n\ell m}(\mathit{r}),}\end{array}$$(3)

where we have introduced the instantaneous amplitude a_nℓm. The eigenfunctions satisfy the eigenvalue relation

$$\begin{array}{c}\hfill {\mathcal{L}}^{\mathrm{ad}}\left({\mathit{\xi }}_{n\ell m}\right)={\omega }_{n\ell m}^2{\mathit{\xi }}_{n\ell m},\end{array}$$(4)

with

$$\begin{array}{c}\hfill {\mathit{\xi }}_{n\ell m}(\mathit{r})={\xi }_{n\ell m}^r(r)Y_{\ell }^m(\theta ,\phi ){\mathit{e}}_r+{\xi }_{n\ell m}^h(r)r\mathbf{\nabla }{Y}_l^m(\theta ,\phi ),\end{array}$$(5)

where ω_nℓm is the angular eigenfrequency; (r, θ, φ) are the spherical coordinates in the stellar frame; ${\xi }_{n\mathcal{l}m}^r$ and ${\xi }_{n\mathcal{l}m}^h$ are the (real) radial and poloidal components of the eigenfunctions, respectively, which are normalized such as ${\xi }_{n\mathcal{l}m}^r$ = 1 at the photosphere; and $Y_{\mathcal{l}}^m$ are the orthonormal spherical harmonics. The eigenfunctions are orthogonal with respect to the density-weighted inner product and their mode mass is defined as

$$\begin{array}{c}\hfill {\mathcal{M}}_{n\ell m}\equiv {\displaystyle {\int }_V\rho {\mathit{\xi }}_{n\ell m}·{\mathit{\xi }}_{n\ell m}^{\star }\text{d}V,}\end{array}$$(6)

where V is the stellar volume beyond which the stellar density vanishes and (^⋆) denotes the complex conjugate.

2.2.2. Globally quasi-adiabatic gravity modes

Asymptotic gravity modes are incompressible (e.g., Dintrans & Rieutord 2001), so that ℒ^nad1 and ℒ^nad2 in Eq. (2) are dominated by second-order derivatives with respect to radius, and their local norms scale as $H^2/{\lambda }_{n\mathcal{l}m}^2$ when applied on a harmonic (n, ℓ, m), where λ_nℓm(r) is the local radial wavelength (see Appendix A.2.2). Consequently, the time evolution of the amplitude a_nℓm is simultaneously governed by the local damping timescale given by $t_{\text{damp}}\sim {({\lambda }_{n\mathcal{l}m}^2/H)}^2{t}_{\text{R}}$ in Eq. (2) and the dynamical timescale t_dyn ∼ 1/ω_nℓm in Eq. (1). Almost everywhere in the star the quasi-adiabatic limit is supposed to be met, that is,

$$\begin{array}{c}\hfill t_{\mathrm{dyn}}\ll t_{\mathrm{damp}}(r),\end{array}$$(7)

and the non-adiabatic effects represented by the quantity δS/c_p can be locally treated as a small perturbation (e.g., Dziembowski 1977b). This is not met in a very thin near-surface layer where the density vanishes and t_R, which is equivalent to the thermal timescale in this region, becomes much smaller than dynamical timescale (see, e.g., Fig. 2b in Berthomieu & Provost 1990). However, this region is so thin and the part of the total mode energy contained inside is so small that its impact on the global mode damping is expected to be negligible. This expectation is again supported by the numerical computations of Belkacem et al. (2009), who demonstrated that the work performed by the radiative flux variations on the oscillations mainly originates from the radiative cavity for the considered frequency range.

In order to reason in a global way we therefore define a global damping timescale, denoted T_damp, which aims to measure the impact of the non-adiabatic effects on the global behavior of a harmonic (n, ℓ, m). Owing to the place of t_damp in Eq. (2), T_damp is taken equal to the inverse of the harmonic mean of 1/t_damp(r) weighted by the local mode energy. Considering only the contribution from the radiative zone to the mode damping, T_damp reduces to Eq. (A.54) in the asymptotic frequency range, that is,

$$\begin{array}{c}\hfill T_{\mathrm{damp}}^{-1}\approx {\displaystyle \frac{{\int }_0^{r_{\mathrm{b}}}{t}_{\mathrm{damp}}^{-1}\text{d}r/{\lambda }_{n\ell m}}{{\int }_0^{r_{\mathrm{b}}}\text{d}r/{\lambda }_{n\ell m}}},\end{array}$$(8)

where r_b denotes the radius of the base of the convective zone. In the following, we thus assume that the oscillations are globally quasi-adiabatic in the considered frequency range:

$$\begin{array}{c}\hfill t_{\mathrm{dyn}}\ll T_{\mathrm{damp}}.\end{array}$$(9)

The definition in Eq. (8) appears to be relevant and eases the computation of the mode amplitude.

2.2.3. Forced amplitude

Projecting Eq. (1) on ξ_nℓm and expressing δS by means of Eq. (2), a set of coupled third-order linear differential equations for the dependent variables a_nℓm(t) can then be obtained, each of them accounting for the first-order non-adiabatic perturbations (see Appendix A.2.3). These equations involve a fast timescale, which is t_dyn, and a slow timescale, which is actually T_damp defined in Eq. (8). It thus appears judicious to solve the differential system using a two-timing method (e.g., Kevorkian 1961). This method can provide us with a uniformly valid solution up to timescales on the order of T_damp. A two-timing analysis is sufficient for the present purpose since most of the mode energy is dissipated on a timescale on the order of T_damp, so that the error made is expected to be small compared to other uncertainties related, for instance, to the modeling of the generation process.

Expressing explicitly the partial derivatives with respect to the fast and slow timescales in the amplitude equations, and grouping together the terms of the same order, a two-timing approach ultimately leads to the mode amplitude (see Appendix A.2.4)

$$\begin{array}{c}\hfill a_{n\ell m}(t)\approx A_{n\ell m}(t)e^{+i{\omega }_{n\ell m}t-{\eta }_{n\ell m}t}+A_{n\ell m}^{\star }(t)e^{-i{\omega }_{n\ell m}t-{\eta }_{n\ell m}t},\end{array}$$(10)

with

$$\begin{array}{c}\hfill A_{n\ell m}(t)\approx \frac{1}{2i}{\displaystyle {\int }_{-\infty }^t{\omega }_{n\ell m}{\stackrel{\sim }{\mathcal{F}}}_{n\ell m}(t^{\prime })e^{-i{\omega }_{n\ell m}{t}^{\prime }+{\eta }_{n\ell m}{t}^{\prime }}\text{d}{t}^{\prime },}\end{array}$$(11)

where ${\stackrel{\sim }{\mathcal{F}}}_{n\ell m}$ is provided by Eq. (A.27). This term results from the scalar projection of the plume driving term in Eq. (1) on ξ_nℓm normalized by the mode mass. In addition, the expression of the damping rate η_nℓm is compatible, within the asymptotic limit, with the usual quasi-adiabatic expression provided, for example, by Dziembowski et al. (2001) or Godart et al. (2009). Its expression is provided in Eqs. (A.59) and (A.60) from which we deduce that η_nℓm ∼ 1/T_damp. This confirms a posteriori our choice for the definition of the global damping timescale in Eq. (8).

At this point we note that Eq. (10) can be retrieved if we assume that the solutions of the homogeneous equations for each amplitude a_nℓm take the form of an exponentially damped harmonic oscillator, as usually done in quasi-adiabatic analyses (e.g., Dziembowski 1977b; Unno et al. 1989), and then use this ansatz to find the particular solution of the forced equations. Nevertheless, the present analysis also clearly highlights two essential points. First, Eq. (10) is uniformly valid up to a timescale on the order of T_damp only. Second, it appears that Eq. (10) holds true if the coupling induced by the ℒ^nad operator in Eq. (1) between the different (n, ℓ, m) components remains negligible. Assuming that modes with adjacent radial orders n and n + 1 have comparable amplitudes for given values of ℓ and m, this is met if we have (see Appendix A.2.4)

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{t}_{\mathrm{core}}}{T_{\mathrm{damp}}}\approx n\frac{t_{\mathrm{dyn}}}{T_{\mathrm{damp}}}\ll 1,\end{array}$$(12)

where Δt_core is the time spent by a wave energy ray of frequency ω_nℓm to cross the radiative core. In this case the radiative losses accumulated during a one-travel path of a wave through the radiative core are negligible; as a result, the resonant modes keep at leading order the same global structure as the eigenfunctions of ℒ^ad, except for the small exponential decay of their amplitudes with time (see discussion in Appendix A.2.4). We show in Sect. 3 that this is met for a solar model in the considered frequency range. The relative error made when writing Eq. (10) thus turns out to be on the order of Δt_core/T_damp at most.

To express further Eq. (10) we need to model the ${\stackrel{\sim }{\mathcal{F}}}_{n\ell m}$ term, and thus to specify the plume velocity field.

2.3. Modeling the plumes and the penetration zone

While the mean ensemble behavior of convective plumes has been widely studied in numerical simulations, for example, how they structure together, how they can merge, or how efficiently they can transport energy (e.g., von Hardenberg et al. 2008; Pieri et al. 2016; Pratt et al. 2017), to the authors’ knowledge there is no quantitative description of their detailed structures. Moreover, though useful to guide the investigations, the outcomes of the current numerical simulations of extended convective envelopes still have to be taken with caution (e.g., Rempel 2004).

In the present work we thus choose to follow an analytical description based on the model of Pinçon et al. (2016). The penetration of convective plumes is supposed to be a random, stationary, and ergodic process. We assume that, on average, 𝒩 identical plumes (i.e., with the same velocity and shape) are penetrating at each time and can transfer their energy into gravity modes. The plumes are supposed to be incoherent with each other and uniformly distributed on the sphere. At the base of the convective region, the Péclet number, which represents the ratio of the efficiency of the advection of heat by the plumes to the efficiency of the radiative diffusion, is expected to be much larger than unity in the Sun (see Sect. 3.2). As a result, the convective plumes are braked over a very small penetration length (i.e., relative to the local pressure scale height) below the radius where the Schwarzschild criterion is met (Zahn 1991; Dintrans et al. 2005). The plume ram pressure gradient is thus supposed to be maximum inside the penetration region, and to vanish outside. The analytical form of the velocity field inside the penetration region that is associated with a plume penetrating at t = 0 and whose center has for latitudinal and azimuthal coordinates (θ₀, φ₀) reads (Pinçon et al. 2016, cf. Eqs. (19) and (20))

$$\begin{array}{c}\hfill {\mathcal{V}}_{\mathrm{p},0}(\mathit{r},t)=f\left(\frac{t}{{\tau }_{\mathrm{p}}}\right){\mathcal{V}}_r(r)e^{-S_h^2/2b^2}{\mathit{e}}_r,\end{array}$$(13)

where the f function represents the time evolution of the plume velocity, with τ_p the plume lifetime, 𝒱_r is the radial profile, b is the plume radius, and S_h(r; θ₀, φ₀) corresponds to the distance on a concentric sphere from the center of the plume.

Owing to the large Péclet number, an adiabatic temperature gradient is imposed in the penetration zone; 𝒱_r can thus follow the model of penetrative convection of Zahn (1991). The value of the plume velocity at the entry of the penetration region, denoted V_b, and that of the radius b are adapted from the model of turbulent plumes developed by Rieutord & Zahn (1995). Below this region, the transition from an adiabatic to a radiative temperature gradient is supposed to be very sharp, so that the Brunt-Väisälä frequency discontinuously changes from about zero in the adiabatic region to N_t at the top of the radiative zone. We note that this hypothesis is supported by the recent seismic inversions of the Brunt-Väisälä frequency in the Sun performed by Buldgen et al. (2020). Their observations show that the transition from an adiabatic to a radiative gradient at the base of the convective zone occurs over a distance d_trans representing about 0.5% of the solar radius and that the value of N_t is equal to about 550 µHz. As a result, we find that d_trans is much smaller than the local wavelength in the considered ranges of frequencies and degrees. This justifies the fact that, from the point of view of the gravity modes, the profile of the Brunt-Väisälä frequency can be supposed to be discontinuous at the top of the radiative zone.

Finally, the time evolution of the plumes in the penetration region is badly understood. Pinçon et al. (2016) identified two probable plume destruction processes: baroclinic instabilities and turbulence. Using orders of magnitude, they estimated that the plume lifetime is most likely in the range around the convective turnover timescale of turbulent eddies at the base of the convective region, τ_conv. Regarding the profile, they assumed that the penetration of a plume is a short-lived event and follows a Gaussian law in time, as in the initial work of Townsend (1966). Nevertheless, given the related uncertainties, it is relevant to test other prescriptions. Within the framework of the stochastic excitation of gravity modes, the spectral density of turbulent kinetic energy of convective eddies was also usually assumed to be Gaussian until Belkacem et al. (2009), based on 3D numerical simulations, showed that a Lorentzian form is more appropriate. Gaussian and Lorentzian laws represent two limiting cases in the time Fourier domain, since they stand for very rapidly and very slowly decreasing functions of frequency, respectively. By analogy, the spectral density of kinetic energy of the plumes will be also supposed to be Lorentzian, which, in the time domain, is equivalent to assuming an exponential time profile. Therefore, the f(t) function in Eq. (13) will be equal in the two limiting cases to either

$$\begin{array}{c}\hfill f_{\mathrm{G}}\left(\frac{t}{{\tau }_{\mathrm{p}}}\right)\equiv e^{-t^2/{\tau }_{\mathrm{p}}^2}\mathrm{or}{f}_{\mathrm{E}}\left(\frac{t}{{\tau }_{\mathrm{p}}}\right)\equiv e^{-|t|/{\tau }_{\mathrm{p}}}.\end{array}$$(14)

2.4. Mean mode energy

Owing to the random properties of the generation process, a statistical approach has to be considered. As the time-averaged mode amplitude vanishes, we need to estimate either the mean mode square amplitude or the mean mode energy. We find it more appropriate to reason in a first step on the mean mode energy since it does not depend on the chosen normalization of the eigenfunction basis, whereas the amplitude does. Because of the incoherence between the convective plumes, we first note that the wave velocity fields generated by different plumes cancel each other out. Assuming in addition that excitation is a stationary process and the plumes are all identical and uniformly distributed over the sphere, we show in Appendix A.3.1 that the total mean mode energy is merely equal to the oscillation energy generated by one single plume penetrating at t = 0 and averaged over its angular position (θ₀, φ₀), multiplied by the instantaneous number of penetrating plumes 𝒩. Therefore, using Eq. (13) in Eq. (10), it is possible to compute the mean oscillation energy associated with each orthogonal (n, ℓ, m) harmonic. In the asymptotic limit, the WKB form of the eigenfunctions can be used to formulate the result analytically (see Appendix A.3.2). In the case of a large Péclet number at the base of the convective region, that is, for a very small penetration region, we find that the mean energy of the (n, ℓ, m) harmonic can be finally expressed as (see Appendix A.3.5)

$$\begin{array}{c}\hfill ⟨E_{n\ell m}⟩\approx {\displaystyle \frac{\left[({\omega }_{n\ell m}\mathrm{\Delta }{\mathrm{\Pi }}_{\ell }/{\pi }^2)\overline{L_{\mathrm{p}}}{F}_{\mathrm{d},\ell }{e}^{-\ell (\ell +1)b^2/2r_{\mathrm{b}}^2}{\mathcal{C}}_{n\ell m}\right]}{2{\eta }_{n\ell m}}},\end{array}$$(15)

with

$$\begin{array}{c}\hfill \overline{L_{\mathrm{p}}}=\mathcal{A}{\mathcal{S}}_{\mathrm{p}}\frac{{\rho }_{\mathrm{b}}{V}_{\mathrm{b}}^3}{2},\end{array}$$(16)

where ΔΠ_ℓ is the asymptotic period spacing between two consecutive adiabatic gravity modes of degree ℓ given in Eq. (A.39), $\overline{L_{\mathrm{p}}}$ is the mean plume kinetic luminosity through the shell of radius r_b at the base of the convective zone, $\mathcal{A}=\mathcal{N}{b}^2/4r_{\text{b}}^2$ is the plume filling factor, 𝒮_p = πb² is the area occupied by a single plume, ρ_b and V_b represent the density and the plume velocity at r_b, and F_d, ℓ = V_bk_h, b/N_t is the Froude number at the top of the radiative zone, with $k_{h,\mathrm{b}}=\sqrt{\ell (\ell +1)}/r_{\mathrm{b}}$ the horizontal wavenumber of the mode.

We note that the numerator of Eq. (15) represents the amount of power injected into the mode per unit of time. The term inside the brackets refers to the mode mass. The Gaussian term represents the horizontal correlation between the plumes and the mode, while 𝒞_nℓm measures the temporal correlation. The general analytical expression of 𝒞_nℓm is provided by either Eq. (A.56) in the case of a Gaussian plume time evolution (i.e., f = f_G) or Eq. (A.55) in the case of an exponential plume time evolution (i.e., f = f_E). In the considered frequency range we usually have η_nℓm ≪ ν_p ≪ ω_nℓm, where ν_p = 1/τ_p. As a result, 𝒞_nℓm is merely equal to (see Appendix A.3.3)

$$\begin{array}{c}\hfill {\mathcal{C}}_{n\ell m}\approx 4\sqrt{\pi }\frac{{\eta }_{n\ell m}}{{\nu }_{\mathrm{p}}}\frac{{\nu }_{\mathrm{p}}^3}{{\omega }_{n\ell m}^3}\mathrm{if}f=f_{\mathrm{G}},\end{array}$$(17)

and

$$\begin{array}{c}\hfill {\mathcal{C}}_{n\ell m}\approx 16\frac{{\nu }_{\mathrm{p}}^3}{{\omega }_{n\ell m}^3}\mathrm{if}f=f_{\mathrm{E}}.\end{array}$$(18)

In the considered frequency range, the temporal correlation is thus expected to be much smaller in the Gaussian case than in the exponential case; as a result, the mean mode energy is also expected to be smaller. This will be discussed in greater detail in Sect. 3.

Moreover, within the quasi-adiabatic and asymptotic limits ${\eta }_{n\mathcal{l}m}\propto \mathcal{l}(\mathcal{l}+1)/{\omega }_{n\mathcal{l}m}^2$, as shown by Eq. (A.60) and illustrated in Fig. 1 (e.g., Godart et al. 2009). In the considered frequency range Eqs. (15)–(18) thus demonstrate that ⟨E_nℓm⟩ is independent of the frequency when f = f_E, whereas ⟨E_nℓm⟩ is inversely proportional to the squared frequency when f = f_G. In both cases the mean mode energy decreases as ℓ increases. As in the case of low-frequency progressive internal gravity waves, we see according to Eq. (15) that the excitation efficiency is mainly measured by the Froude number at the top of the radiative region, which is expected to be much smaller than unity in stars (see Sect. 3.1). As the correlation and mode mass terms are also smaller than unity, the mean mode energy turns out to be much smaller than the mean plume kinetic energy at the base of the convective zone, and we check a posteriori that the feedback from the modes to the plume dynamics inside the penetration region is negligible.

Fig. 1. Radiative damping rate as a function of the oscillation frequency νnℓm and the angular degree ℓ. Each circle corresponds to a given eigenmode.

2.5. Apparent mode radial velocity

Owing to a high signal-to-noise ratio, the search for the solar gravity modes is usually performed through radial velocity measurements integrated over the full solar disk (e.g., Appourchaux et al. 2010). In order to be compared to the observed measurements, the theoretical mean mode energy has to be converted into mean mode radial velocity accounting for the line-of-sight projection and the limb-darkening effects (Dziembowski 1977a; Berthomieu & Provost 1990). Among other observational effects that need to be corrected, they are expected to be the predominant ones. Following the computation of Belkacem et al. (2009), we show in Appendix C that the mean apparent radial mode velocity for a (n, ℓ, m) harmonic is equal, in the slow rotator limit, to

$$\begin{array}{c}\hfill v_{n\ell m}^{\mathrm{app}}=\sqrt{\frac{⟨E_{n\ell m}⟩}{{\mathcal{M}}_{n\ell m}}}{\stackrel{\sim }{V}}_{n\ell m},\end{array}$$(19)

where ${\stackrel{\sim }{V}}_{n\ell m}$ is identified as

$$\begin{array}{c}\hfill {\stackrel{\sim }{V}}_{n\ell m}=|{\alpha }_{\ell }^m{\xi }_{n\ell m}^r(r_{\mathrm{abs}})+{\beta }_{\ell }^m{\xi }_{n\ell m}^h(r_{\mathrm{abs}})|,\end{array}$$(20)

with r_abs the radius in the atmosphere where the absorption line considered to measure the radial velocity is formed. The visibility coefficients ${\alpha }_{\mathcal{l}}^m$ and ${\beta }_{\mathcal{l}}^m$ are provided in Eqs. (C.6) and (C.7). They depend on the limb-darkening law and on the angle between the stellar rotation axis and the line of sight.

We note that, contrary to the mode energy, it is not possible to express Eq. (20) using the WKB form of the eigenfunctions as it becomes questionable in the convective region and the atmosphere of stars. The estimate of Eq. (20) therefore requires the numerical computation of the mode eigenfunctions and their mode masses.

3. Applications to gravity modes in the Sun

In this section the excitation model is used to predict the apparent surface radial velocity of gravity modes generated by penetrative convection for a solar model. As in Belkacem et al. (2009), we choose to compare our results to the detection threshold of the GOLF instrument on board the SoHO spacecraft (e.g., Gabriel et al. 1995).

3.1. Solar model

We consider a calibrated solar model computed with the stellar evolution code CESTAM (Marques et al. 2013). The chemical composition follows the solar mixture, as given in Asplund et al. (2009), with initial helium and metal abundances Y₀ = 0.25 and Z₀ = 0.013. The OPAL 2005 equation of states and opacity tables as well as the NACRE nuclear reaction rates were used to build the model. The atmosphere was constructed following an Eddington gray approximation, and convection was modeled using the mixing-length theory with a parameter α_MLT = 1.62. Microscopic diffusion, overshooting, and rotation were not taken into account. We note that the obtained solar structure, although imprecise from the point of view of seismic inversions, is sufficient for our purpose as the uncertainties related to the assumptions considered in the excitation model are dominant. For example, it is known that microscopic diffusion can slightly modify the stratification at the base of the convective zone (e.g., Christensen-Dalsgaard et al. 1993) and hence can affect the properties of the penetrating plumes. Nevertheless, this effect is expected to be much smaller than the uncertainties related, for instance, to the plume formation mechanism across the convection zone.

The adiabatic eigenfrequencies ν_nℓm = ω_nℓm/2π, eigenfunctions, and mode masses needed to compute Eq. (19) were obtained via the oscillation code ADIPLS (Christensen-Dalsgaard 2008). The eigenfunctions were normalized such that ${\xi }_{n\mathcal{l}m}^r$ = 1 at the photosphere. Owing to the large drop in the density and pressure at the stellar surface, the reflective zero-boundary conditions are considered. We check that the period spacing between consecutive ℓ = 1 eigenmodes is very close to the asymptotic value, that is, ΔΠ_ℓ = 1 ≈ 26 min. Moreover, the mode mass is nearly proportional to ℓ³ and ${\nu }_{n\ell m}^{-6}$ (see, e.g., Fig. 6 in Belkacem et al. 2009). To be complete, the radiative damping rate occurring in Eq. (15) and computed according to the quasi-adiabatic expression provided in Eq. (A.60) is plotted in Fig. 1. In this figure, we first see that the damping rate is very similar to that computed numerically by Belkacem et al. (2009, see Fig. 9) who also accounted for the influence of the upper layers. This confirms our assumption that the work performed in the inner radiative cavity is the main contributor to the damping of the asymptotic gravity modes. Second, we also see in this figure that the damping timescale T_damp ∼ 1/η_nℓm is more than six orders of magnitude higher than the oscillation timescale t_dyn ∼ 1/ω_nℓm. This justifies a posteriori the use of the global quasi-adiabatic approximation in this frequency range. Moreover, since n ≲ 10³ in the considered frequency range, it is obvious that Δt_core ≈ nt_dyn ≪ T_damp, so that the hypothesis made in Eq. (12) appears a posteriori to be valid as well (i.e., the coupling induced by non-adiabatic effects between adjacent radial orders is negligible up to a timescale on the order of T_damp).

3.2. GOLF apparent radial velocity with standard parameters

The internal structure of the solar model provides us with all the equilibrium quantities to compute the parameters of the excitation model and the mode apparent velocity in Eq. (19). The radius at the base of the convective zone and the local density are equal to r_b ≈ 5.2 10⁵ km and ρ_b ≈ 125 kg m⁻³. Based on the model of plumes of Rieutord & Zahn (1995), we find b ≈ 10⁴ km and V_b ≈ 190 m s⁻¹. The Péclet number at the base of the convective region is estimated to be about equal to V_bH(r_b)/K_rad(r_b)∼10⁷, where H and K_rad are the temperature scale height and the radiative diffusivity, respectively. The large Péclet number assumption used in Sect. 2.3 is thus justified for the Sun. Following Buldgen et al. (2020), the Brunt-Väisälä frequency just below r_b is taken equal to N_t = 550 µHz. The Froude number at the top of the radiative zone is thus equal to F_d, 1 ≈ 10⁻³. Moreover, we use a reasonable value 𝒩 ≈ 1000, as previously estimated by Rieutord & Zahn (1995), which corresponds to 𝒜 ≈ 0.08, in qualitative agreement with previous numerical simulations (e.g., Stein & Nordlund 1998; Brummell et al. 2002). The convective turnover timescale at the base of the convective region is equal to about τ_conv ≈ 10 days according to the convection mixing length theory, which leads to ν_p ≈ 1/τ_conv ≈ 1 µHz. At this point, the considered values of the plume parameters are referred to as the standard values. We consider that the sodium NaD1 and NaD2 absorption lines used by the GOLF instrument to measure radial velocities at height h ≈ 300km above the photosphere (Bruls & Rutten 1992). We also use the limb-darkening law as formulated by Ulrich et al. (2000), and take for the angle between the stellar rotation axis and the line of sight Θ₀ ≈ 83°. Moreover, we systematically focus on the azimuthal numbers |m|=ℓ in what follows as we find they are less attenuated by the visibility effects and we want to study the more optimistic case.

Using all these physical ingredients, we can compute the apparent radial velocity of the asymptotic gravity modes, $v_{n\ell \ell }^{\mathrm{app}}$, between 10 µHz and 100 µHz. The result is plotted in Fig. 2 as a function of the mode eigenfrequencies and for typical degrees ℓ = 1−5. Both limiting cases of a Gaussian and exponential plume time evolution are considered and compared.

Fig. 2. Apparent surface radial velocity of solar gravity modes as a function of the oscillation frequency νnℓm and the angular degree ℓ, for azimuthal numbers m = ℓ and for standard plume parameters. The results obtained with a Gaussian and an exponential plume time evolution (see Sect. 2.3) are represented by the diamonds and the plus signs, respectively, for each considered eigenmode. The thick dashed gray line corresponds to the 22-year GOLF detection threshold (see Sect. 3.3). The blue dash-dotted line represents the upper limit for the amplitude of the dipolar modes when accounting for the largest plausible variations in the plume parameters (i.e., with a factor of two increase in the inverse plume lifetime and radius; see Sect. 4.3).

First, Fig. 2 shows in both cases that the apparent surface velocity of the gravity modes generated by penetrative convection is maximum for ℓ = 1 and sharply drops as a function of ℓ. At a given frequency, the difference of apparent velocity between the ℓ = 1 and ℓ = 5 modes is larger than about two orders of magnitude. In addition, $v_{n\ell \ell }^{\mathrm{app}}$ appears to depend linearly on the frequency in the exponential case, while it is almost independent of frequency in the Gaussian case. In order to disentangle the reasons for such trends from the different terms occurring in Eq. (19), we propose to rewrite Eq. (20) in a more simple way. Owing to the reflective boundary conditions, the Lagrangian perturbation of pressure vanishes at the solar surface, which leads to (e.g., Unno et al. 1989)

$$\begin{array}{c}\hfill {\xi }_{n\ell m}^h(R_{\odot })=\frac{G{M}_{\odot }}{R_{\odot }^3{\omega }_{n\ell m}^2}{\xi }_{n\ell m}^r(R_{\odot })\approx {\left(\frac{100\mu \mathrm{Hz}}{{\nu }_{n\ell m}}\right)}^2,\end{array}$$(21)

where G is the gravitation constant, M_⊙ and R_⊙ are the solar mass and radius, and where we considered the normalization ${\xi }_{n\mathcal{l}m}^r(R_{\odot })=1$. Therefore, over the considered frequency range, we have ${\xi }_{n\mathcal{l}m}^h(R_{\odot })\ge {\xi }_{n\mathcal{l}m}^r(R_{\odot })$. Moreover, we can assume that ${\xi }_{n\mathcal{l}m}^h$ does not significantly vary in the solar atmosphere, that is, ${\xi }_{n\mathcal{l}m}^h(r_{\text{abs}})\approx {\xi }_{n\mathcal{l}m}^h(R_{\odot })$. We checked the validity of this approximation using the numerical eigenfunctions. Since $\left|{\alpha }_{\mathcal{l}}^{\mathcal{l}}\right|$ is approximately equal to or smaller than $\left|{\beta }_{\mathcal{l}}^{\mathcal{l}}\right|$ (see Table 1), Eq. (19) can thus be expressed at first approximation as

$$\begin{array}{c}\hfill v_{n\ell \ell }^{\mathrm{app}}\approx \sqrt{\frac{⟨E_{n\ell \ell }⟩}{{\mathcal{M}}_{n\ell \ell }}}\left|{\beta }_{\ell }^{\ell }\right|{\left(\frac{100\mu \mathrm{Hz}}{{\nu }_{n\ell \ell }(\mu \mathrm{Hz})}\right)}^2.\end{array}$$(22)

For ℓ = 1−5, the Gaussian term in Eq. (15) remains close to unity since ℓ ≪ r_b/b ∼ 50. As a consequence, since ${\eta }_{n\mathcal{l}m}\propto \mathcal{l}(\mathcal{l}+1)/{\nu }_{n\mathcal{l}m}^2$, we find according to Sect. 2.4 that $⟨E_{n\ell m}⟩\propto {\nu }_{n\ell m}^{-2}$ in the Gaussian case and ⟨E_nℓm⟩∝ℓ⁻² in the exponential case. In addition, as the mode mass depends approximately on ν_nℓm and ℓ as ${\mathcal{M}}_{n\ell m}\propto {\ell }^3{\nu }_{n\ell m}^{-6}$, Eq. (22) leads in the Gaussian case to

$$\begin{array}{c}\hfill v_{n\ell \ell }^{\mathrm{app}}\propto \left|{\beta }_{\ell }^{\ell }\right|{\ell }^{-3/2},\end{array}$$(23)

and in the exponential case to

$$\begin{array}{c}\hfill v_{n\ell \ell }^{\mathrm{app}}\propto \left|{\beta }_{\ell }^{\ell }\right|{\ell }^{-5/2}{\nu }_{n\ell \ell }.\end{array}$$(24)

As shown in Table 1, the value of $\left|{\beta }_{\mathcal{l}}^{\mathcal{l}}\right|$ smoothly varies with ℓ, and the trends predicted in Eqs. (23) and (24) are in qualitative agreement with Fig. 2. It thus appears that the simultaneous influence of the mode mass and the damping rate counterbalances the decrease in the mode driving with frequency (via the temporal correlation term). This explains the trends observed in Fig. 2.

Moreover, we clearly show in Fig. 2 that the value of $v_{n\ell \ell }^{\mathrm{app}}$ predicted assuming a Gaussian plume time evolution at a given ℓ is at least three orders of magnitude smaller than that assuming an exponential plume time evolution. This huge difference results from the much smaller temporal correlation between the plumes and the considered eigenmodes in case of a Gaussian time evolution, which is represented by the 𝒞_nℓm term. In this case the spectral density of the plume kinetic energy is mostly carried by the frequencies lower than ν_p. As ν_nℓm ≫ ν_p, the coupling between the plumes and the eigenmodes is very small, and the excitation is very inefficient. In contrast, in the case of an exponential time evolution, the spectral density of the plume kinetic energy is distributed over higher frequencies, and the coupling between the plumes and the eigenmodes is maximum for frequencies around ν_nℓm (see the computation in Appendix A.3.3). As a result, the energy transfer from the plumes to the modes is much more efficient. Considering Eqs. (17) and (18) in Eqs. (15) and (19), the ratio of the mode apparent velocity in the Gaussian case to that in the exponential case scales as $\sqrt{{\eta }_{n\ell m}/{\nu }_{\mathrm{p}}}\ll 1$. Since ${\eta }_{n\ell m}\propto 1/{\nu }_{n\ell m}^2$ in the quasi-adiabatic limit, this ratio is expected to decrease as 1/ν_nℓm, in agreement with Fig. 2. Our estimate therefore demonstrates that the amplitude of the asymptotic solar gravity modes generated by penetrative convection critically depends on the plume time evolution at the base of the convective zone. Our lack of knowledge on the plume dynamics in this region thus prevents us from quantitatively predicting the efficiency of this excitation process. In turn, if detected, solar gravity modes can be expected to put important constraints on penetrative convection at the top of the radiative zone (see discussion in Sect. 4.4).

Finally, we point out that over the considered frequency range the GOLF apparent radial velocity of the gravity modes computed with the present model and standard plume parameters is about one order of magnitude smaller than that predicted considering turbulent pressure as the driving mechanism. For the ℓ = 1 modes in the exponential case, for which the excitation is the most efficient, our predictions lie between $v_{n\ell \ell }^{\mathrm{app}}\approx 0.01-0.05$ cm s⁻¹ from ν_nℓm = 10 µHz to ν_nℓm = 100 µHz, while those of Belkacem et al. (2009) lie between $v_{n\ell \ell }^{\mathrm{app}}\approx 0.1-0.5$ m s⁻¹ (see Sect. 4.2 for a more detailed comparison). These results therefore suggest that penetrative convection is less efficient than turbulent convection to generate asymptotic gravity modes in the Sun. Nevertheless, owing to a significant sensitivity of the predictions to the plume parameters (i.e., velocity, radius, lifetime), whose values are currently affected by uncertainties, we note that reasonable variations in these parameters are likely to reduce the gap between the two excitation mechanisms (see blue dash-dotted line in Fig. 2). This point is discussed further in Sect. 4.3.

3.3. Comparison with the GOLF detection threshold

Using a statistical approach, Appourchaux et al. (2010) could analytically express the threshold signal-to-noise ratio above which a peak in the power spectral density (PSD) of an observed time series can be considered as a relevant oscillating signal over a frequency interval Δf, with a false alarm probability p_th that the measurement is due to pure noise. As shown in Appendix D, this threshold can be converted into a detection threshold for the root mean square oscillation velocity, denoted with v_th, above which a signal would be detected with a certain confidence level. It reads

$$\begin{array}{c}\hfill v_{\mathrm{th}}\approx \sqrt{ln\left(\frac{T\mathrm{\Delta }f}{p_{\mathrm{th}}}\right)\frac{\stackrel{\sim }{s}}{T}},\end{array}$$(25)

where T is the observation time and $\stackrel{\sim }{s}$ is the mean noise level in the PSD over the range Δf.

We consider a GOLF observation time equal to T = 22 yr. The mean noise level $\stackrel{\sim }{s}$ is estimated from García et al. (2007, see Appendix D). We choose Δf = 10 µHz and p_th = 0.01. Instead of a common rejection level of 10%, we adopt a rejection level of 1% since it allows the posterior probability that a peak is due to noise to reach lower values close to 10%, as already mentioned in Appourchaux et al. (2010). The resulting 22-year GOLF detection threshold is represented by a thick gray dashed line in Fig. 2. We see that the threshold velocity smoothly decreases from about 1 cm s⁻¹ at ν_nℓm = 10 µHz to about 0.5 cm s⁻¹ at ν_nℓm = 100 µHz. This trend results from the decrease in the solar granulation noise with frequency in the GOLF data. Its value is about five orders of magnitude higher than the apparent mode velocity throughout the considered frequency range when assuming a Gaussian plume time evolution. According to Eq. (25) the amplitude of gravity modes will be lower even if the observation time is equal to ten billion years. In the case of an exponential plume time evolution, the amplitude predicted with standard plume parameters is also at least one order of magnitude lower than the 22-year GOLF detection threshold over the considered frequency range. Based on our simple estimate of the GOLF mean noise level (see Eq. (D.2)), Eq. (25) requires an unreachable GOLF observation time T ∼ 2000 yr to detect the dipolar n = 6 gravity mode at ν_nℓm ≈ 90 µHz. In conclusion, using standard values for the plume parameters, our model indicates that the solar asymptotic gravity modes generated by penetrative convection remains far from being detected in an acceptable observation time by radial velocity measurements such as performed by the GOLF instrument. It is worth mentioning that the uncertainties on the predicted amplitude, which are inherent to uncertainties on the plume parameters, may nevertheless considerably reduce the observation time required for a reliable detection. This potential is also evaluated in Sect. 4.3.

4. Discussion

4.1. Comparison with the estimate of Andersen (1996)

Andersen (1996) first investigated the generation of solar gravity modes by penetrative convection using both numerical simulations and simple energy considerations. Including an ad hoc forcing term at the top of the radiative region that mimics the influence of penetrative plumes in a solar model, he numerically solved the wave equations and computed the transmission of the generated wave energy throughout the convective envelope, the structure of which was taken from the turbulent numerical simulations of Andersen (1994). He found that the wave transmission up to the surface is on the order of 0.05%. Using these results, he could estimate that the horizontal velocity of gravity modes near the solar surface ranges between 0.01−1 mm s⁻¹ for ℓ = 6 and ν_nℓm = 50−100 µHz. When accounting for an appropriate GOLF visibility factor, that is $\left|{\beta }_6^6\right|\approx 0.03$, this upper value turns out to be much larger than our predictions in the case of a Gaussian plume time evolution, but on the same order of magnitude in the exponential case. However, we note that this qualitative agreement is more likely to be fortuitous. In order to estimate the amplitude of gravity modes, Andersen (1996) assumed that the total mode energy is equal to the part of the convective kinetic energy transferred on average to progressive gravity waves in the numerical simulation of Andersen (1994), the upper and lower boundary conditions of which are open. As a consequence, this estimate actually did not properly account for the quasi-stationarity of the modes or for the damping processes. Therefore, although qualitatively comparable for the ℓ = 6 gravity modes, the physical origins of our predictions radically differ from those of Andersen (1996), and drawing conclusions from their comparison ultimately appears irrelevant.

4.2. Comparison with turbulence-induced gravity modes

Belkacem et al. (2009) modeled the generation of asymptotic gravity modes by considering the turbulent Reynolds stress as the driving mechanism. As mentioned at the end of Sect. 3.2, the surface mode velocity they predicted is one order of magnitude higher than our estimate in the exponential case when considering standard values for the plume parameters based on semi-analytical models and orders of magnitude. Nevertheless, for both excitation mechanisms, the dependence of the apparent radial velocities on ℓ and ν_nℓm turns out to be comparable (e.g., Belkacem et al. 2009, see Fig. 11). This similar dependence results on the one hand from the mode mass and the damping rate, which play the same role in both excitation processes, and on the other hand from a comparable temporal correlation between the driving source and the modes. Belkacem et al. (2009) assumed that the time coherence of the convective eddies is exponential; as a result, the temporal correlation between the modes and the eddies has a similar form to that between the modes and the plumes when assuming an exponential plume time profile (e.g., a dependence on about ${\nu }_{n\ell m}^{-3}$; see Eqs. (1)–(2) and Fig. 3 in Belkacem et al. 2009). Regarding the magnitude, as the velocity of the turbulent eddies that drive gravity modes is about equal to v_MLT ≈ 60 m s⁻¹ (according to the mixing length theory), which is about three times lower than the plume velocity at the base of the convective zone, penetrative convection could have been at first sight expected to generate gravity modes with higher amplitudes. However, while the excitation by penetrative convection occurs in a very thin shell at the base of the convective zone, convective eddies can drive gravity modes in a larger volume of the envelope. Moreover, as the eddy turnover frequencies are higher than ν_p in the upper layers of the convective region, the time correlation is larger for the excitation by turbulent convection than by penetrative plumes. The lower velocities of the convective eddies are thus compensated by a more extended excitation region and a better temporal correlation with the modes, leading to higher amplitudes in the exponential case.

In contrast, when assuming that the eddy-time coherence is Gaussian, as in Kumar et al. (1996), Belkacem et al. (2009) found the apparent radial velocities of gravity modes generated by turbulent eddies is on the order of 10⁻³ cm s⁻¹. Although their predictions also depend on ν_nℓm and ℓ in a similar way to our predictions in the Gaussian case, their magnitude is much larger than ours. While the same arguments as in the previous paragraph hold true to explain the similarities in the dependence on ν_nℓm and ℓ, the temporal correlation with the modes in the Gaussian case is very sensitive to the timescales associated with the driving source. As the values of the turnover frequencies of the turbulent eddies driving the modes are larger than ν_p, the temporal correlation is expected to be much larger than for the excitation by plumes. This certainly explains the huge difference between the excitation by turbulent convection and penetrative convection in the Gaussian case.

4.3. Sensitivity of the mode amplitude to the plume parameters

As shown by Pinçon et al. (2016), the efficiency of the wave excitation by penetrative convection depends in a significant way on the plume parameters b, ν_p, and V_b, whose values are subject to uncertainties. For example, the plume model of Rieutord & Zahn (1995) is expected to provide an upper limit on V_b since it does not take into account the upward counterflow that can exchange momentum with the downward plumes and hence slow them down (e.g., Rempel 2004). In contrast, the plume radius b is likely to be underestimated since the model of Rieutord & Zahn (1995) does not account for the possible clustering of plumes, as observed for instance in numerical simulations (von Hardenberg et al. 2008). In addition, as already discussed in Pinçon et al. (2016), turbulence inside the penetration region could lead to values of ν_p that are substantially larger than the convective eddy turnover frequency predicted by the MLT.

To investigate the influence of these uncertainties on the mode amplitudes, we arbitrarily consider in this section the effect of a decrease in V_b of 30%, and an increase in ν_p and b by a factor of two, while keeping the filling factor 𝒜 constant (i.e., 𝒩 ≈ 250 when the value of b is two times larger). Neglecting the variations of the horizontal correlation term in Eq. (15) for ℓ = 1−5, Eqs. (15)–(19) show that, at given frequency and degree, $v_{n\ell \ell }^{\mathrm{app}}\propto b{V}_{\mathrm{b}}^2{\nu }_{\mathrm{p}}^{3/2}$ in the exponential case and $v_{n\ell \ell }^{\mathrm{app}}\propto b{V}_{\mathrm{b}}^2{\nu }_{\mathrm{p}}$ in the Gaussian case. We see that a decrease in V_b of 30% leads to a decrease by a factor of two in $v_{n\ell \ell }^{\mathrm{app}}$, while an increase of b or ν_p by a factor of two results in an increase by a factor of between two or three. As the mode amplitudes are insignificantly low in the Gaussian case, such variations in the parameters, and even much larger variations, do not modify the large gap with the exponential case and the GOLF detection threshold; gravity modes thus remain undetectable in this case. In the exponential case, such variations in the parameters b and ν_p can in contrast affect the predictions about the detectability of the modes. Assuming an increase in ν_p by a factor of two, which results in an increase by a factor of about three in the amplitudes, the GOLF observation time required to detect a plume-induced dipolar gravity modes at ν_nℓm = 100 µHz is decreased by about one order of magnitude (i.e., to T ∼ 200 yr). Considering the largest plausible variations with a simultaneous increase in b by a factor of two, the observation time required to detect such a mode is reduced to T ∼ 50 yr. In this most favorable case considered here, the predicted amplitudes of the dipolar plume-induced gravity modes are plotted in Fig. 2 (blue dash-dotted line). This upper limit turns out to be close to the amplitudes of the turbulent-induced gravity modes predicted by Belkacem et al. (2009).

In conclusions, the time evolution of the plumes at the base of the convective zone represents the major source of uncertainties in our model. The error related to the uncertainties on the plumes parameters is insignificant compared to the effect of the assumption on the plume time evolution. In addition, while the plume-induced gravity modes remain undetectable in the Gaussian assumption whatever the values of the plume parameters, the uncertainties on these parameters can significantly modify the predictions of the detectability of few modes in the exponential assumption. In the most favorable plausible case, the GOLF observation time required to detect the plume-induced ℓ = 1 gravity modes around ν_nℓm ≈ 100 µHz is reduced to T ∼ 50 yr, with amplitudes close to those predicted considering turbulent pressure as the driving mechanism.

4.4. Can the amplitudes of gravity modes bring constraints on penetrative convection?

We see in Sect. 4.3 that our lack of knowledge on penetrative convection affects our predictions. In turn, we can wonder to what extent the measurement of the amplitudes of gravity modes in the Sun can provide us with constraints on penetrative convection and the properties at the base of the convective zone.

Our model shows either that the contribution from penetrative convection to the mode amplitude is undetectable for a Gaussian-like plume time evolution or that it behaves as a function of the frequency and the degree in a similar way to the contribution from turbulent convection for an exponential plume time profile. As both excitation mechanisms are subject to uncertainties (Belkacem et al. 2009), it can thus appears difficult at first sight to disentangle the two contributions and interpret the observed amplitudes in terms of structure without any independent information. Assuming that the plume time evolution is exponential, the measurement of the mode amplitude can nevertheless easily put upper limits on each excitation process, by ensuring that the predictions in each case remain lower than the observations. As the solar gravity modes have not been detected yet, we can similarly proceed by ensuring that the predictions remain lower than the current detection threshold. Imposing that $v_{\mathrm{th}}\gtrsim v_{n\ell \ell }^{\mathrm{app}}$ for ν_nℓm = 100 µHz and ℓ = 1 in the exponential case, we find, using Eqs. (15), (18), (19) and (25), that for 22-year GOLF observations

$$\begin{array}{c}\hfill F_{\mathrm{d},1}\left(\frac{\overline{L_{\mathrm{p}}}}{L_{\odot }}\right){\left(\frac{{\nu }_{\mathrm{p}}}{1\mu \mathrm{Hz}}\right)}^3\lesssim {10}^{-6},\end{array}$$(26)

where L_⊙ is the solar luminosity.

Although the potential of the gravity mode amplitudes to probe penetrative convection is highlighted here, the present work suggests that the current theoretical uncertainties on the modeling of penetrative convection, and in particular on the plume time evolution, would limit their physical interpretation. The issue of disentangling the contributions of different sources to the excitation is especially brought to the fore. Theoretical efforts, by means of numerical simulations and semi-analytical models, and possible complementary observational constraints will be needed in the future to improve our modeling of this phenomenon and to elaborate relevant and robust seismic diagnoses to interpret the amplitude of gravity modes.

5. Conclusions

In this work our aim was to estimate the amplitude of the asymptotic solar gravity modes generated by penetrative convection in the frequency range between 10 µHz and 100 µHz. Following Pinçon et al. (2016), we considered the ram pressure of an ensemble of incoherent, uniformly distributed convective plumes penetrating into the top layers of the radiative zone as the driving mechanism. The forced oscillation equation was solved in the global quasi-adiabatic approximation using a two-timing method. As a result, we obtain an analytical expression of the mean mode energy, which is converted into apparent radial mode velocity through appropriate visibility factors. The standard plume modeling (i.e., plume radius, velocity, and lifetime) follows semi-analytical models, and their time evolution in the penetration region is assumed to be either Gaussian or exponential, by analogy with stochastic excitation by turbulent eddies. The apparent mode radial velocity is computed for a solar model in both cases and the result is compared to the 22-year GOLF detection threshold.

We find that the mean mode energy drastically depends on the assumption about the time evolution of the plumes inside the penetration region. On the one hand, in the limiting case of a Gaussian time evolution, asymptotic gravity modes turn out to be undetectable by means of radial velocity measurements such as those performed by the GOLF instrument. This is the consequence of a plume lifetime that is too long compared to the oscillation period. This result holds true despite a wide range of values considered for the parameters of the model. In the other limiting case of an exponential time evolution we find that penetrative convection can generate gravity modes in a much more efficient way than in the Gaussian case. In this case, the lower the angular degree or the higher the frequency, the larger the apparent mode radial velocity. Using standard values for the plume parameters, the apparent radial mode velocity appears to reach about 0.05 cm s⁻¹ for ℓ = 1 and ν_nℓm ≈ 100 µHz. These predictions are one order of magnitude smaller than those predicted considering turbulent pressure as the driving mechanism, and remain well below the current 22-year GOLF detection threshold. Nevertheless, accounting for uncertainties in the plume parameters, we find in the most favorable plausible case that the predicted apparent mode radial velocity can be increased by a factor of six, that is, around 0.3 cm s⁻¹ for ℓ = 1 and ν_nℓm ≈ 100 µHz. The observation time required to detect such a mode is reduced to about 50 yr with the GOLF instrument. These variations mainly result from an important sensitivity of the mode amplitude to the plume parameters, and, in contrast, from a limited sensitivity of the detection threshold to the observation time. Our findings thus indicate that, in the most favorable plausible case, penetrative convection can drive asymptotic gravity modes as efficiently as turbulent convection and with amplitudes close to the detection limit. We highlight that, if detected, the measurement of the gravity mode amplitudes is expected to bring constraints on penetrative convection, but that the current uncertainties on the modeling of penetrative convection, and in particular their temporal evolution, will certainly limit their physical interpretation.

The results of this work call for further studies, either observational or theoretical. First, our estimate of the excitation by penetrative plumes and the previous estimates of the excitation by turbulent pressure clearly suggest that we are likely to be very close to the detection in the asymptotic frequency range, and encourage us to carry on our efforts in observations and data analyses. While we mainly focused on the 22-year GOLF data, we note that a myriad of other data is available, for instance the observations by the GONG and BiSON ground-based telescope networks, and are an important source of information to be analyzed. Second, it will be interesting in the future to extend the theoretical predictions to a higher frequency range. A simple extrapolation of the available predictions toward slightly higher frequencies suggests that the amplitudes of the gravity modes in this domain are also likely to be close to the current detection limit. As already pointed out by Belkacem et al. (2011), predicting the amplitude of such low radial order gravity modes will require us to account consistently for the interplay between oscillations and convection, which is challenging since it will call for combining a proper non-local time-dependent treatment of convection with a fully non-adiabatic treatment of pulsations. Furthermore, new developments, based both on numerical simulations and semi-analytical models, are needed to improve our understanding of the behavior of downward convective plumes at the interface with the radiative region. Although it is still a challenge to accurately simulate stellar regimes, our hope is that the promising combination of such theoretical advancements with future measurements of the gravity mode amplitudes will bring constraints on the dynamics and the mixing at work at the base of the convective zone.

Acknowledgments

We thank the anonymous referee for his careful reading and relevant comments that greatly helped improving the manuscript. We sincerely acknowledge K. Belkacem and M.-A. Dupret for the very interesting discussions on the present subject and their sensible comments. During this work, C. P. was funded by a postdoctoral fellowship of Chargé de Recherche from F. R. S.-FNRS (Belgium). G. B. acknowledges fundings from the SNF AMBIZIONE grant No 185805 (Seismic inversions and modelling of transport processes in stars). C. P. warmly thanks M. Huet, V. Huet, A. Leguillon and C. Houdmond for their sincere friendship and their encouragements during the writing of this paper.

References

Appendix A: Energy spectrum of plume-induced gravity modes

In the following we derive the energy spectrum of asymptotic gravity modes generated by penetrative convection. The mathematical derivation largely uses and sometimes extends the analysis of stellar oscillations presented in Unno et al. (1989), and also uses the model of Pinçon et al. (2016).

A.1. Forced oscillation equations

Following the model of Pinçon et al. (2016), the linearized momentum, continuity, Poisson’s equations, and the equation of state¹ read

$$\begin{array}{cc}& \rho {\partial }_t^2\mathit{\xi }+\mathbf{\nabla }{p}^{\prime }-{\rho }^{\prime }\mathit{g}+\rho \mathbf{\nabla }{\psi }^{\prime }=-\mathbf{\nabla }·(\rho {\mathcal{V}}_{\mathrm{p}}\otimes {\mathcal{V}}_{\mathrm{p}}),\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill & \delta \rho +\rho \mathbf{\nabla }·\mathit{\xi }=0,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill & {\mathrm{\nabla }}^2{\psi }^{\prime }=4\pi G{\rho }^{\prime },\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill & \frac{\delta \rho }{\rho }=\frac{1}{{\mathrm{\Gamma }}_1}\frac{\delta p}{p}-v_T\frac{\delta S}{c_p},\hfill \end{array}$$(A.4)

where p and g are respectively the equilibrium pressure and gravitational acceleration, ψ′ is the perturbation of the gravitational potential, Γ₁ is the first adiabatic index, G is the gravitation constant, v_T = −(∂lnρ/∂lnT)_p, and all the other quantities are introduced in Sect. 2.1. In these four equations, X′ and δX denotes the Eulerian and Lagrangian perturbations of the quantity X, respectively. Using the usual solution of Poisson’s equation

$$\begin{array}{c}\hfill {\psi }^{\prime }(\mathit{r},t)=-{\displaystyle {\int }_{V}G\frac{{\rho }^{\prime }(\mathit{x},t)}{|\mathit{r}-\mathit{x}|}\text{d}V,}\end{array}$$(A.5)

where V is the stellar volume beyond which the stellar density vanishes, as well as Eqs. (A.2) and (A.4), it is possible to rewrite Eq. (A.1) in the form of Eq. (1) in which

$$\begin{array}{c}\hfill {\mathcal{L}}^{\mathrm{nad}}\left(H\frac{\delta S}{c_p}\right)=\frac{1}{\rho }\mathbf{\nabla }\left({\mathrm{\Gamma }}_1{v}_{T}p\frac{\delta S}{c_p}\right),\end{array}$$(A.6)

with H = −(dr/dlnT) the temperature scale height and T the equilibrium temperature, and where ℒ^ad is the adiabatic differential linear operator defined as (e.g., Chandrasekhar 1964; Unno et al. 1989, see Chap. 14.3 for details)

$$\begin{array}{cc}& {\mathcal{L}}^{\mathrm{ad}}\left(\mathit{\xi }\right)=\frac{\mathbf{\nabla }p}{{\rho }^2}\mathbf{\nabla }·\left(\rho \mathit{\xi }\right)-\frac{1}{\rho }\mathbf{\nabla }(\mathit{\xi }·\mathbf{\nabla }p)-\frac{1}{\rho }\mathbf{\nabla }(\rho c^2\mathbf{\nabla }·\mathit{\xi })\hfill \\ \hfill & +\mathbf{\nabla }\left({\int }_{V}G\frac{{\mathbf{\nabla }}_{\mathit{x}}·[\rho (\mathit{x})\mathit{\xi }(\mathit{x})]}{|\mathit{r}-\mathit{x}|}\text{d}{\mathit{x}}^3\right),\hfill \end{array}$$(A.7)

with c² = Γ₁p/ρ the squared sound speed. We note that H has been introduced in Eq. (A.6) for ℒ^nad to represent the same physical quantity as ℒ^ad (i.e., a squared frequency). Equation (1) has to be completed by the energy equation that specifies the evolution of δS. Neglecting the contribution of the convective flux for asymptotic gravity modes (Belkacem et al. 2009) and the nuclear energy production rate in the considered layers, it reads

$$\begin{array}{c}\hfill T{\partial }_t\delta S=-\delta (\frac{1}{\rho }\mathbf{\nabla }·{\mathcal{F}}_{\mathrm{R}}),\end{array}$$(A.8)

where ℱ_R is the radiative flux.

A.2. Mode amplitude

A.2.1. Energy equation in the diffusion approximation

At this stage the oscillation equations are represented by Eqs. (1) and (A.8). Furthermore, we thus need to express the Lagrangian perturbation of the radiative flux. In the diffusion approximation the radiative flux is equal to

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{R}}=-\rho c_p{K}_{\mathrm{rad}}\mathbf{\nabla }T=-\frac{16{\sigma }_{\mathrm{SB}}{T}^3}{3\rho \kappa }\mathbf{\nabla }T,\end{array}$$(A.9)

where K_rad is the radiative diffusivity, σ_SB is the Stefan-Boltzmann constant, and κ is the Rosseland mean opacity. Perturbing this equation and replacing δρ/ρ by −∇ ⋅ ξ, Eq. (A.8) can be rewritten

$$\begin{array}{cc}& {\partial }_t\left(\frac{\delta S}{c_p}\right)=-\frac{1}{\rho c_{p}T}\mathbf{\nabla }·[-F_{\mathrm{R}}H\mathbf{\nabla }\left(\frac{\delta T}{T}\right)+(4\frac{\delta T}{T}+\mathbf{\nabla }·\mathit{\xi }-\frac{\delta \kappa }{\kappa }){\mathit{F}}_{\mathrm{R}}]\hfill \\ \hfill & +\frac{1}{\rho c_{p}T}\mathbf{\nabla }·[\mathbf{\nabla }({\mathit{F}}_{\mathrm{R}}·\mathit{\xi })-\mathit{\xi }(\mathbf{\nabla }·{\mathit{F}}_{\mathrm{R}})],\hfill \end{array}$$(A.10)

where F_R = F_Re_r is the radial equilibrium radiative flux. We note that the right-hand side of Eq. (A.10) is equivalent to the right-hand side of Eq. (21.14) of Unno et al. (1989) with δε_N = 0. Considering the equation of state T(S, ρ) and the opacity table κ(T, ρ) while using the continuity equation, it becomes

$$\begin{array}{cc}& \frac{\delta T}{T}=\gamma \frac{\delta S}{c_p}-{\mathrm{\Gamma }}_1{\mathrm{\nabla }}_{\mathrm{ad}}\mathbf{\nabla }·\mathit{\xi }\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill & \frac{\delta \kappa }{\kappa }=\gamma {\kappa }_T\frac{\delta S}{c_p}-({\kappa }_T{\mathrm{\Gamma }}_1{\mathrm{\nabla }}_{\mathrm{ad}}+{\kappa }_{\rho })\mathbf{\nabla }·\mathit{\xi },\hfill \end{array}$$(A.12)

where γ = c_p/c_V, with c_p and c_V the specific heat capacities at constant pressure and volume, respectively; ∇_ad is the adiabatic temperature gradient; κ_T = (∂lnκ/∂lnT)_ρ; and κ_ρ = (∂lnκ/∂lnρ)_T. Using Eqs. (A.11) and (A.12) in Eq. (A.10), the evolution of the Lagrangian perturbation of entropy is ruled by

$$\begin{array}{c}\hfill {\partial }_t\left(\frac{\delta S}{c_p}\right)=\frac{1}{t_{\mathrm{R}}}[{\mathcal{L}}^{\mathrm{nad}1}\left(\frac{\mathit{\xi }}{H}\right)+{\mathcal{L}}^{\mathrm{nad}2}\left(\frac{\delta S}{c_p}\right)],\end{array}$$(A.13)

where we have introduced the local radiative thermal timescale

$$\begin{array}{c}\hfill t_{\mathrm{R}}=\frac{\rho c_{p}TH}{F_{\mathrm{R}}},\end{array}$$(A.14)

and where ℒ^nad1 and ℒ^nad2 are two linear operators that involve derivatives with respect to the normalized variable r/H and whose expressions can be readily deduced from Eqs. (A.10)–(A.12).

A.2.2. Local scaling of the oscillation equations

Before going further, it is instructive to express locally the oscillation equations in a dimensionless form. In the following we focus on a mode with a characteristic angular frequency and a local wavelength denoted σ and λ(r, σ), respectively. The dynamical timescale is thus defined as t_dyn ≡ σ⁻¹. First, according to Eq. (4), the local norm of the ℒ^ad operator² is equal to σ² for such a mode. Owing to the place of ℒ^nad in Eq. (1), it seems reasonable to assume that its norm also scales as σ². This hypothesis is checked a posteriori in Sect. A.3.4. Equation (1) can thus be rewritten as

$$\begin{array}{c}\hfill {\partial }_{\tau }^2\mathit{\xi }+{\stackrel{\mathbf{\sim }}{\mathcal{L}}}^{\mathrm{ad}}\left(\mathit{\xi }\right)+{\stackrel{\mathbf{\sim }}{\mathcal{L}}}^{\mathrm{nad}}\left(H\frac{\delta S}{c_p}\right)=-\frac{1}{\rho {\sigma }^2}\mathbf{\nabla }·(\rho {\mathcal{V}}_{\mathrm{p}}\otimes {\mathcal{V}}_{\mathrm{p}}),\end{array}$$(A.16)

where we have defined τ = t/t_dyn and the differential operators ${\stackrel{\mathbf{\sim }}{\mathcal{L}}}^{\mathrm{ad}}={\sigma }^{-2}{\mathcal{L}}^{\mathrm{ad}}$ and ${\stackrel{\mathbf{\sim }}{\mathcal{L}}}^{\mathrm{nad}}={\sigma }^{-2}{\mathcal{L}}^{\mathrm{nad}}$ such that their norms remain on the order of unity. Second, regarding the energy equation, we first note that, owing to the incompressible character of the asymptotic low-frequency gravity modes (e.g., Dintrans & Rieutord 2001),

$$\begin{array}{c}\hfill \mathbf{\nabla }·\mathit{\xi }=\mathcal{O}(|\mathit{\xi }|/H),\end{array}$$(A.17)

where the big-𝒪 Bachmann-Landau notation is introduced. As ℒ^nad1 and ℒ^nad2 in Eq. (A.13) involve third-order and second-order derivatives at most with respect to r/H of ∇ ⋅ ξ and δS/c_p, respectively, their norms thus at most scale as (H²/λ²) since λ ≲ H and Eq. (A.17) are met within the asymptotic limit. Equation (A.13) can hence be locally expressed as

$$\begin{array}{c}\hfill {\partial }_{\tau }\left(\frac{\delta S}{c_p}\right)={\epsilon }_{\sigma }[{\stackrel{\sim }{\mathcal{L}}}^{\mathrm{nad}1}\left(\frac{\mathit{\xi }}{H}\right)+{\stackrel{\sim }{\mathcal{L}}}^{\mathrm{nad}2}\left(\frac{\delta S}{c_p}\right)].\end{array}$$(A.18)

In this equation, we have defined

$$\begin{array}{c}\hfill {\epsilon }_{\sigma }(r)\equiv \frac{t_{\mathrm{dyn}}}{t_{\mathrm{damp}}(r,\sigma )}\mathrm{and}{\stackrel{\sim }{\mathcal{L}}}^{\mathrm{nad}k}\equiv \frac{{\lambda }^2(r,\sigma )}{H^2(r)}{\mathcal{L}}^{\mathrm{nad}k}\end{array}$$(A.19)

in such a way that the norms of ${\stackrel{\sim }{\mathcal{L}}}^{\mathrm{nad}k}$ for k = 1, 2 also remain on the order of unity, and where we can identify the local damping timescale as

$$\begin{array}{c}\hfill t_{\mathrm{damp}}(r,\sigma )\equiv \frac{{\lambda }^2(r,\sigma )}{H^2(r)}{t}_{\mathrm{R}}(r).\end{array}$$(A.20)

A.2.3. Globally quasi-adiabatic oscillations

At this point, we focus on a component (n, ℓ, m) with the amplitude a_nℓm(t) and the frequency ω_nℓm, and fix σ = ω_nℓm and ε_σ ≡ ε_nℓm. According to Eq. (A.18), δS locally scales as

$$\begin{array}{c}\hfill \frac{\delta S}{c_p}=\mathcal{O}\left({\epsilon }_{n\ell m}\frac{|{\mathit{\xi }}_{n\ell m}|}{H}\right).\end{array}$$(A.21)

In a given frequency range around ω_nℓm the damping timescale is expected to be much longer than the oscillation period almost everywhere in the star. As a result, ε_nℓm ≪ 1 and δS/c_p ≪ |ξ_nℓm|/H. In the limit of ε_nℓm → 0, δS → 0 and the problem locally tends toward the adiabatic case. Nevertheless, there obviously exists a very thin near-surface layer where ε_nℓm ≫ 1 and the quasi-adiabatic approximation locally fails, so that the ordering in Eq. (A.21) is not met. Therefore, in order to reason from a global point of view, we define the global damping timescale as the inverse of the harmonic mean of the local damping timescale throughout the star weighted by the local mode energy, that is,

$$\begin{array}{c}\hfill T_{\mathrm{damp}}^{-1}\equiv \frac{1}{{\mathcal{M}}_{n\ell m}}{\displaystyle {\int }_V{t}_{\mathrm{damp}}^{-1}\rho {\mathit{\xi }}_{n\ell m}·{\mathit{\xi }}_{n\ell m}^{\star }\text{d}V,}\end{array}$$(A.22)

where T_damp has to be related to the local wavelength of ξ_nℓm through Eqs. (A.20) and (A.22), and where ℳ_nℓm is the mode mass defined in Eq. (6). In the following we thus assume that the oscillations are globally quasi-adiabatic, that is,

$$\begin{array}{c}\hfill {\delta }_{n\ell m}\equiv t_{\mathrm{dyn}}/T_{\mathrm{damp}}\ll 1.\end{array}$$(A.23)

For convenience, we also define the quantity

$$\begin{array}{c}\hfill {\chi }_{n\ell m}={\mathcal{M}}_{n\ell m}{\delta }_{n\ell m}.\end{array}$$(A.24)

The global quasi-adiabatic approximation expressed by Eq. (A.23) will greatly ease the derivation of the mode amplitude, as we show here.

Taking the derivative of Eq. (A.16) with respect to τ, injecting Eqs. (3) and (A.18) into the obtained expression, and computing the inner product with ${\mathit{\xi }}_{n\mathcal{l}m}^{\star }$, we obtain a differential equation ruling the mode amplitude a_nℓm(τ). We point out that, owing to the spherical symmetry of the star, the inner product with ${\mathit{\xi }}_{n\mathcal{l}m}^{\star }$ selects the angular degree ℓ and the azimuthal number m in the expansion of the perturbations onto the orthonormal basis of the spherical harmonics (see, e.g., the definitions of the differential operators). Therefore, the subscript ℓm is often dropped to simplify the notation in what follows. The amplitude equation finally reads, at first order in δ_n,

$$\begin{array}{c}\hfill {\partial }_{\tau }^3{a}_n+{\partial }_{\tau }{a}_n-2{\delta }_n{\displaystyle ({\stackrel{\sim }{\eta }}_{\mathit{nn}}{a}_n+\sum _{n^{\prime }\ne n}{\stackrel{\sim }{\eta }}_{n{n}^{\prime }}{a}_{n^{\prime }})={\partial }_{\tau }{\stackrel{\sim }{\mathcal{F}}}_n+\mathcal{O}({\delta }_n^2),}\end{array}$$(A.25)

where

$$\begin{array}{cc}& {\stackrel{\sim }{\eta }}_{n{n}^{\prime }}=-\frac{1}{2{\chi }_{n\ell m}}{\displaystyle {\int }_V\rho {\mathit{\xi }}_{n\ell m}^{\star }·{\stackrel{\sim }{\mathcal{L}}}^{\mathrm{nad}}\left[H{\epsilon }_{n^{\prime }\ell m}{\stackrel{\sim }{\mathcal{L}}}^{\mathrm{nad}1}\left(\frac{{\mathit{\xi }}_{n^{\prime }\ell m}}{H}\right)\right]\text{d}V,}\hfill \end{array}$$(A.26)

$$\begin{array}{cc}\hfill & {\stackrel{\sim }{\mathcal{F}}}_n=-\frac{1}{{\omega }_{n\ell m}^2{\mathcal{M}}_{n\ell m}}{\displaystyle {\int }_V\mathbf{\nabla }·(\rho {\mathcal{V}}_{\mathrm{p}}\otimes {\mathcal{V}}_{\mathrm{p}})·{\mathit{\xi }}_{n\ell m}^{\star }\text{d}V.}\hfill \end{array}$$(A.27)

The second-order term in Eq. (A.25) is related to the Lagrangian perturbation of entropy and results from the scaling in Eq. (A.21) and the global quasi-adiabatic hypothesis in Eq. (A.23). In addition, the product ${\delta }_n{\stackrel{\sim }{\eta }}_{\mathit{nn}}$ stands for the damping rate (in terms of the τ variable). In contrast, the sum in the brackets on the left-hand side of Eq. (A.25) encapsulates coupling terms between the components n and n′ ≠ n. Such coupling terms result, by construction, from the expansion onto the eigenfunctions of the ℒ^ad operator as they are not natural solutions of the non-adiabatic problem. When writing Eq. (A.25), we have implicitly assumed ${\stackrel{\sim }{\eta }}_{n{n}^{\prime }}=\mathcal{O}(1)$. As shown in Eq. (A.26), the integral in the numerator is at most on the order of χ_nℓm since the norms of the operators are on the order of unity. Moreover, this integral corresponds to the inner product between an oscillating radial function, with n radial nodes and a characteristic wavelength λ_n, and another oscillating radial function, with about n′ radial nodes and a characteristic wavelength λ_n′. For n′ ∼ n we therefore have $|{\stackrel{\sim }{\eta }}_{n{n}^{\prime }}|\sim 1$, whereas for n′ ≫ n or n′≪n both oscillating functions are incoherent with each other and the inner product vanishes, that is, $|{\stackrel{\sim }{\eta }}_{n{n}^{\prime }}|\ll 1$. This justifies ${\stackrel{\sim }{\eta }}_{n{n}^{\prime }}=\mathcal{O}(1)$.

A.2.4. Two-timing analysis

As shown by Eq. (A.25), the evolution of the mode amplitude is ruled by a fast dynamical timescale, T₀ = τ = t/t_dyn, and a slow damping timescale, T₁ = δ_nτ = t/T_damp. Using a two-timing perturbation method therefore seems judicious to solve analytically this equation (e.g., Kevorkian 1961). This method can provide us with a uniformly valid solution up to a timescale on the order of 𝒪(1/δ_n) in terms of the independent variable τ, or equivalently 𝒪(T_damp) in terms of the dependent variable t.

Homogeneous equation. Taking advantage of the linearity of the problem, we search in a first step for a general solution of the homogeneous amplitude equation. For all n, we start from the two-timing perturbation ansatz, that is,

$$\begin{array}{c}\hfill a_n(\tau )=a_n^0(T_0,T_1)+{\delta }_n{a}_n^1(T_0,T_1)+\mathcal{O}({\delta }_n^2),\end{array}$$(A.28)

keeping in mind that τ, T₀, and T₁ also depend on n through the scaling by 1/ω_nℓm. Injecting Eq. (A.28) into Eq. (A.25) and using ∂_τ = ∂_T0 + δ_n∂_T1 via the chain rule, the homogeneous equation reads at leading order in δ_n

$$\begin{array}{c}\hfill {\partial }_{T_0}^3{a}_n^0+{\partial }_{T_0}{a}_n^0=0,\end{array}$$(A.29)

and at first order in δ_n

$$\begin{array}{cc}& {\partial }_{T_0}^3{a}_n^1+{\partial }_{T_0}{a}_n^1+\underset{\mathrm{Secular}\mathrm{terms}}{\underbrace{(3{\partial }_{T_0}^2{\partial }_{T_1}{a}_n^0+{\partial }_{T_1}{a}_n^0-2{\stackrel{\sim }{\eta }}_{\mathit{nn}}{a}_n^0)}}\hfill \\ \hfill & =2{\displaystyle \sum _{\begin{array}{c}{n}^{\prime }\ne n\end{array}}{\stackrel{\sim }{\eta }}_{n{n}^{\prime }}{a}_{n^{\prime }}^0.}\hfill \end{array}$$(A.30)

The general solution of Eq. (A.29) then reads

$$\begin{array}{c}\hfill a_n^0(T_1,T_0)={\alpha }_n(T_1)e^{i{T}_0}+{\beta }_n(T_1)e^{-i{T}_0},\end{array}$$(A.31)

where α_n(T₁) and β_n(T₁) are two functions to be determined with the first-order equation. Using Eq. (A.31) in Eq. (A.30), it is straightforward to show (e.g., via the method of variation of parameters) that, if non-null, the terms in brackets will induce terms proportional to τe^±iτ in the expression of $a_n^1$, up to a factor on the order of α_n and β_n. This will make $a_n^1$ on the order of $\mathcal{O}(a_n^0/{\delta }_n)$ up to a timescale on the order of 𝒪(1/δ_n) and, therefore, would break the perturbation ansatz in Eq. (A.28). Therefore, to obtain a uniformly valid expression up to such timescales, the terms in brackets have to vanish. This condition leads to

$$\begin{array}{c}\hfill {\alpha }_n(T_1)=A_n{e}^{-{\stackrel{\sim }{\eta }}_{\mathit{nn}}{T}_1}\mathrm{and}{\beta }_n(T_1)=B_n{e}^{-{\stackrel{\sim }{\eta }}_{\mathit{nn}}{T}_1},\end{array}$$(A.32)

with A_n and B_n two arbitrary constants. As a consequence, the amplitude $a_n^1$ is ruled only by the coupling term on the right-hand side of Eq. (A.30), which depends on the leading-order amplitudes $a_{n^{\prime }}^0$ with n′ ≠ n. For n′ ≠ n, it is obvious that the amplitudes follow the same zeroth-order and first-order equations as in Eqs. (A.29) and (A.30); however, because of the considered scaling, T_k for k = 0 or 1 in these equations are replaced by $T_k^{\prime }={\omega }_{n^{\prime }}t=\varpi T_k$ with ϖ = (ω_n′/ω_n). As a result, we get

$$\begin{array}{c}\hfill a_{n^{\prime }}^0(T_1,T_0)=A_{n^{\prime }}{e}^{-{\stackrel{\sim }{\eta }}_{n^{\prime }{n}^{\prime }}\varpi T_1}{e}^{i\varpi T_0}+B_{n^{\prime }}{e}^{-{\stackrel{\sim }{\eta }}_{n^{\prime }{n}^{\prime }}\varpi T_1}{e}^{-i\varpi T_0},\end{array}$$(A.33)

where A_n′ and B_n′ are two arbitrary constants. According to Eqs. (A.30)–(A.32) we thus find that $a_n^1$ is ruled at leading order by

$$\begin{array}{c}\hfill {\partial }_{T_0}^3{a}_n^1+{\partial }_{T_0}{a}_n^1=2{\displaystyle \sum _{n^{\prime }\ne n}{\stackrel{\sim }{\eta }}_{n{n}^{\prime }}{a}_{n^{\prime }}^0.}\end{array}$$(A.34)

Using again the method of variation of parameters, we can show that the amplitude $a_n^1$ is equal to a series of terms as a function of n′ that are $\mathcal{O}({\stackrel{\sim }{\eta }}_{nn\prime }{a}_{n\prime }^0/|\varpi -1|$). At this point, we assume that the magnitude of $a_{n^{\prime }}^0$ slowly varies with n′, so that it is $\mathcal{O}(a_n^0)$ for n′∼n; this will be checked a posteriori (see the resulting velocity amplitudes in Fig. 2). Provided ${\stackrel{\sim }{\eta }}_{nn\prime }=\mathcal{O}(1)$ for n ∼ n′ and ${\stackrel{\sim }{\eta }}_{nn\prime }\ll 1$ otherwise, as justified in Appendix A.2.3, the main contribution to $a_n^1$ results from the terms of the series associated with n′ ∼ n; as a result, $a_n^1=\mathcal{O}(a_n^0{\omega }_n/\Delta {\omega }_n)$ with Δω_n = |ω_n − ω_n − 1|. Within the asymptotic limit, ω_n/Δω_n ∼ n ≫ 1 and thus $a_n^1\gg a_n^0$ in order of magnitude. This highlights afterward that, instead of δ_n, it is more appropriate to scale the non-adiabatic perturbation of the mode amplitude by

$$\begin{array}{c}\hfill {\stackrel{\sim }{\delta }}_n\equiv ({\omega }_n/\mathrm{\Delta }{\omega }_n){\delta }_n,\end{array}$$(A.35)

in such a way that the perturbation ansatz becomes

$$\begin{array}{c}\hfill a_n(\tau )=a_n^0(T_0,T_1)+{\stackrel{\sim }{\delta }}_n{\stackrel{\sim }{a}}_n^1(T_0,T_1)+\mathcal{O}({\stackrel{\sim }{\delta }}_n^2),\end{array}$$(A.36)

and the new first-order perturbation ${\stackrel{\sim }{a}}_n^1=\mathcal{O}(a_n^0)$. As a result, the relative error when considering $a_n\approx a_n^0$ up to a timescale on the order of T_damp is on the order of ${\stackrel{\sim }{\delta }}_n$. We check a posteriori at the end of Sect. 3.1 that ${\stackrel{\sim }{\delta }}_n\ll 1$ in the considered asymptotic frequency range. Therefore, all the previous perturbation developments, fortunately, hold true and the coupling terms do not stand for secular terms. These terms thus have no effect on the two-timing solution provided by $a_n^0$, and are just responsible for the relative error on the order of ${\stackrel{\sim }{\delta }}_n$.

In summary, provided that the mode amplitude slowly varies with the radial order, the solution of the homogeneous equation as a function of the variable t, which is uniformly valid up to a timescale on the order of 𝒪(T_damp), takes the form, according to Eqs. (A.31) and (A.32), of the adiabatic solution modulated by an exponential damping term, that is,

$$\begin{array}{c}\hfill a_n(t)\approx A_n{e}^{+i{\omega }_{n}t-{\eta }_{n}t}+B_n{e}^{-i{\omega }_{n}t-{\eta }_{n}t},\end{array}$$(A.37)

where the relative error is on the order of ${\stackrel{\sim }{\delta }}_n$ and the damping rate is equal to

$$\begin{array}{c}\hfill {\eta }_n={\omega }_n{\delta }_n{\stackrel{\sim }{\eta }}_{\mathit{nn}},\end{array}$$(A.38)

with δ_n and ${\stackrel{\sim }{\eta }}_{\mathit{nn}}$ defined in Appendix A.2.3.

Physically speaking, we note that the parameter ${\stackrel{\sim }{\delta }}_n$ measures the coupling induced by the non-adiabatic effects between adjacent harmonics in the decomposition onto the eigenfunctions of the ℒ^ad operator. To show this, it is possible to express this parameter in more sensible ways. On the one hand, the eigenfrequencies of the adiabatic asymptotic gravity modes follow at leading order (e.g., Shibahashi 1979; Tassoul 1980)

$$\begin{array}{c}\hfill {\omega }_n\approx \frac{2\pi }{n\mathrm{\Delta }{\mathrm{\Pi }}_{\ell }},\mathrm{with}\mathrm{\Delta }{\mathrm{\Pi }}_{\ell }=\frac{2{\pi }^2}{\sqrt{\ell (\ell +1)}}{\displaystyle {\left({\int }_0^{r_{\mathrm{b}}}N\frac{\text{d}r}{r}\right)}^{-1},}\end{array}$$(A.39)

where ΔΠ_ℓ is the so-called period-spacing, r_b is the radius of the base of the convective zone, and N is the Brunt-Väisälä frequency defined as

$$\begin{array}{c}\hfill N^2=g(\frac{1}{{\mathrm{\Gamma }}_1}\frac{\text{d}lnp}{\text{d}r}-\frac{\text{d}ln\rho }{\text{d}r}).\end{array}$$(A.40)

Therefore, for n ≫ 1 in the asymptotic regime, we have

$$\begin{array}{c}\hfill \mathrm{\Delta }{\omega }_n\approx \pi {\displaystyle {\left({\int }_0^{r_{\mathrm{b}}}\frac{k_r}{{\omega }_n}\text{d}r\right)}^{-1}\mathrm{with}{k}_r\approx \frac{\sqrt{\ell (\ell +1)}}{r}{(\frac{N}{{\omega }_n}-1)}^{1/2},}\end{array}$$(A.41)

where k_r is the local radial wavenumber. As the radial group velocity of gravity waves is about equal to ω_n/k_r (e.g., Unno et al. 1989), Eq. (A.41) shows that Δω_n ∼ 1/Δt_core, where Δt_core is the time spent by a wave energy ray of frequency ω_n to cross the radiative core. We thus obtain ${\stackrel{\sim }{\delta }}_n\sim \mathrm{\Delta }{t}_{\mathrm{core}}/T_{\mathrm{damp}}$. On the other hand, using η_n ∼ 1/T_damp according to our scaling, we have ${\stackrel{\sim }{\delta }}_n\sim {\eta }_n/\mathrm{\Delta }{\omega }_n$. The hypothesis Eq. (12) is thus equivalent to assuming that the width of each oscillation peak around an adiabatic eigenfrequency is much smaller than the frequency spacing between two consecutive radial orders n and n + 1 in the power oscillation spectrum.

Forced amplitude. In a second step, we search for the solution of the forced amplitude equation that is uniformly valid up to timescales on the order of T_damp and satisfies the initial value condition a_nℓm(−∞) = 0. Using the method of variation of parameters, it is straightforward to show that this solution is

$$\begin{array}{c}\hfill a_n(t)\approx A_n(t)e^{+i{\omega }_{n}t-{\eta }_{n}t}+B_n(t)e^{-i{\omega }_{n}t-{\eta }_{n}t},\end{array}$$(A.42)

where A_n(t) is a function of time provided by Eq. (11) and B_n(t) = A_n(t)^⋆. As a check, Eq. (A.42) leads to

$$\begin{array}{cc}& {\partial }_{\tau }^3{a}_n={\partial }_{\tau }{\stackrel{\sim }{\mathcal{F}}}_n-2{\delta }_n{\stackrel{\sim }{\eta }}_{\mathit{nn}}{a}_n{\stackrel{\sim }{\mathcal{F}}}_n+\mathcal{O}({\delta }_n^2),\hfill \end{array}$$(A.43)

$$\begin{array}{cc}\hfill & {\partial }_{\tau }{a}_n=-2{\delta }_n{\stackrel{\sim }{\eta }}_{\mathit{nn}}{a}_n.\hfill \end{array}$$(A.44)

Injecting Eqs. (A.43) and (A.44) in Eq. (A.25), we see that the secular term represented by 2η_nna_n cancel out and that the relative residual is on the order of ${\stackrel{\sim }{\delta }}_n$ at most, as shown in the previous paragraph.

A.3. Mean mode energy within the asymptotic limit

Within the asymptotic limit, it is possible to use, to good approximation, the leading-order WKB analytical expressions of the eigenfunctions in order to express the mean mode energy in a simple but sensible way.

A.3.1. General expression for the mean mode energy

Owing to the equipartition of the specific kinetic and potential energies in the case of gravity waves (e.g., Lighthill 1978), we define the mean total oscillation energy as

$$\begin{array}{c}\hfill \langle E\rangle =\underset{T\to +\infty }{lim}\frac{1}{T}{\displaystyle {\int }_{-T/2}^{+T/2}\left({\int }_V\rho (r){\left|{\partial }_t\mathit{\xi }(\mathit{r},t)\right|}^2\text{d}V\right)\text{d}t,}\end{array}$$(A.45)

where (⟨ ⋅ ⟩) is the (ensemble or time) mean operator and V represents the stellar volume.

To express Eq. (A.45) in a simple way, the same reasoning as done in Sects. 2.2.1 and 2.2.2 of Pinçon et al. (2016) can be used. While Pinçon et al. (2016) considered progressive wave packets propagating toward the center of the star and never returning back upward, it is in contrast necessary in the present case to consider the modal structure of the oscillations. The only difference is that, instead of considering that the lifetime of the wave packet at a certain point of the radiative zone is on the order of τ_p, as in Pinçon et al. (2016), we have to consider in the present case that the lifetime of the mode oscillation is on the order of T_damp.

The reasoning is as follows. Owing to the linearity of the wave equations, the displacement field at any time can be represented as the superposition of all the waves induced by each individual plume. Each of these waves has a finite lifetime, which is on the order of the damping timescale T_damp. Over the time interval [ − T/2, T/2], it thus appears obvious that the set of waves constituting ξ is generated by a finite number of plumes, denoted 𝒩_T. In the limit of T ⋙ T_damp, as the plume emerging rate is constant and equal to 𝒩/τ_p, with τ_p the characteristic plume lifetime, the number of penetrating plumes that contributes to ξ in the time interval [ − T/2, T/2] is equal at leading order to 𝒩_T ∼ 𝒩ν_pT, with ν_p = 1/τ_p. Because of the incoherence of the convective plumes between each other, the wave velocity fields generated by different plumes negatively interfere. Assuming in addition that the excitation is a stationary process and that the plumes are all identical and uniformly distributed over the sphere, the time and volume integrals in Eq. (A.45) can be written as the mode energy associated with the displacement field ξ₀(r, t; θ₀, φ₀) that is generated by one single plume penetrating at t₀ = 0 in the solid angle dΩ₀ = sin θ₀dθ₀dφ₀ and with the velocity 𝒱_p, 0, multiplied by the number of penetrating plumes over [ − T/2, T/2] (i.e., 𝒩_T ∼ 𝒩ν_pT), and finally averaged over the plume angular position (θ₀, φ₀). Using the Parseval-Plancherel theorem in the time Fourier space, the mean mode energy can thus be expressed as

$$\begin{array}{c}\hfill \langle E\rangle =\mathcal{N}\frac{{\nu }_{\mathrm{p}}}{8{\pi }^2}{\displaystyle {\int }_V{\int }_{{\mathrm{\Omega }}_0}\left({\int }_{-\infty }^{+\infty }\rho (r){\left|\widehat{{\partial }_t{\mathit{\xi }}_0}(\mathit{r},\omega ;{\theta }_0,{\phi }_0)\right|}^2\text{d}\omega \right)\text{d}{\mathrm{\Omega }}_0\text{d}V,}\end{array}$$(A.46)