© C. Neiner et al. 2020

1. Introduction

Be stars are late-O to early-A stars surrounded by a decretion disk fed by discrete mass-loss events (see the review by Rivinius et al. 2013). The ejection of material from the star into the disk requires enough angular momentum to reach the breakup velocity at the stellar surface. Most of the angular momentum is provided by the rapid rotation of the star (on average, 90% of the critical velocity in our galaxy, see Frémat et al. 2005). The physical process leading to the transport of additional angular momentum to the surface has yet to be identified. However, transport by pulsation modes has been put forwards as the most probable explanation up to now (Lee & Saio 1993; Pantillon et al. 2007; Mathis et al. 2008). Indeed, Be stars are known to host pulsations excited by the κ mechanism, similar to those observed in β Cephei and SPB stars (e.g. Neiner et al. 2009). More recently, pulsations that are excited stochastically have also been observed in Be stars (Neiner et al. 2012a) and in O stars (Aerts & Rogers 2015). Gravity waves are also found to be excited in numerical simulations, to efficiently transport angular momentum in hot stars (Rogers et al. 2013), and to provide predictions for rotation profiles that match the seismic observations (Rogers 2015). Moreover, theoretical calculations have shown that rotation is a crucial ingredient for exciting gravito-inertial waves (Mathis et al. 2014).

HD 49330 is a rapidly rotating B0.5IVe star observed with the CoRoT satellite (Auvergne et al. 2009) from October 18, 2007 to March 3, 2008, that is over ∼137 consecutive days with a gap of 3.5 days. During the CoRoT observation, an outburst of several hundredths of magnitude occurred. Outbursts, and their subsequent dimming, are known to occur frequently in this star. This was, however, the first time that an outburst could be observed with an extremely high photometric accuracy and cadence.

The analysis of the CoRoT light curve of HD 49330 (Huat et al. 2009) and the associated ground-based spectroscopic data (Floquet et al. 2009) led to the detection of 30 independent groups of frequencies of variation composed of more than 300 frequencies in total. These frequencies were interpreted in terms of p and g pulsation modes. The amplitude of the various modes are directly correlated with the phases of the outburst: p modes have a higher amplitude during the quiescent and relaxation phases, while g modes increase in amplitude and become prominent just before (precursor phase) and during the outburst; see Huat et al. (2009) for more details. Although the correlation between the changes in the pulsation parameters and the outburst is clear, the mechanism linking the two phenomena still needs to be investigated. This can be done thanks to seismic modelling and this investigation is the purpose of this paper.

To create a seismic model, we need to know the pulsational parameters as well as the fundamental parameters of the star. The pulsational parameters are provided in Huat et al. (2009). Floquet et al. (2009; see their Table 2) determined that for $\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_{\mathrm{c}}}=0.9$ at the equator, T_eff = 28 000 ± 1500 K, log g = 3.9 ± 0.2, v sin i = 280 ± 40 km s⁻¹, i = 43 ± 11 deg, M = 14.4 ± 0.3 M_⊙, R_eq = 8.7 ± 2.5 R_⊙, and log L = 4.5 ± 0.2 L_⊙. While $\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_c}$ is commonly used in observational studies, seismic models rather consider $\overline{\mathrm{\Omega }}$, which is the rotation frequency normalised by $\sqrt{GM/R^3}$, where R is the mean stellar radius. For HD 49330, $\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_c}=0.9$ corresponds to $\overline{\mathrm{\Omega }}=0.7833$.

In Sect. 2, we explain the models we use for the structure of the rapidly rotating star and present our models of κ-driven pulsation modes. We then propose in Sect. 3 an alternative explanation for the CoRoT observations based on stochastically-driven g waves. We discuss the excitation of such waves and present seismic modelling results, including stochastic excitation. We discuss the angular momentum transport by those waves and the ignition of the outburst. Finally, we draw conclusions about the treatment and interpretation of seismic observations of early Be stars and about the impact of stochastic g waves on the ignition of Be outbursts in Sect. 4.

2. Seismic modelling of HD 49330 with the classical κ mechanism

The modelling of the pulsations of a star first requires a calculation of the structure of the star and then a pulsation stability analysis.

2.1. Modelling the rapidly rotating stellar structure

For the stellar structure we use two different approaches:

1. A 2D structure model calculated with the ROTORC code (Deupree 1990, 1995). ROTORC solves the conservation equations for mass, momentum, energy and composition, as well as Poisson’s equation on a 2D spherical grid. The model rotates uniformly on the zero age main sequence and the interior angular momentum is conserved locally as the model evolves. The surface is assumed to be an equipotential. The independent variables are the fractional surface equatorial radius and the colatitude. This allows us to obtain the rotationally flattened structure of the star.

2. A deformation following Chandrasekhar-Milne’s expansion, in which the effects of rotation on the shape of the star are taken into account with a method similar to that of Lee & Baraffe (1995) but slightly modified with the 2D model as input. In this case an isobaric surface is expressed as r = a[1 + α(a) + β(a)P₂(cos θ)] where a is the mean radius, P₂(cos θ) is the Legendre function, and the functions α and β represent the horizontal averaged effects of the centrifugal force and the deformation of the equilibrium state, respectively. Here we obtain the value of β(a) for a given rotation frequency by interpolating the coefficients determined from the 2D ROTORC model. This modification is necessary because Chandrasekhar-Milne’s expansion is only valid when the difference between the rotating and non-rotating structures can be approximated by Ω² order corrections, which is not the case when the rotation velocity is very high.

Details about the above two approaches are described in Sect. 3 of Neiner et al. (2012b).

Figure 1 shows the coefficient β as a function of the fractional radius for a 10 M_⊙ star rotating at various rates. The figure shows that for relatively slow rotation rates, the Chandrasekhar-Milne expansion and the 2D structure model are similar. For $\overline{\mathrm{\Omega }}>0.5$, however, the deformations from the 2D models are considerably larger than those from Chandrasekhar-Milne’s expansion, in particular in the outer layers.

Fig. 1. Coefficient β of the Legendre function P2(cos θ) as a function of the fractional radius for a 10 M⊙ star rotating at various rates. Solid lines represent the deformation from the Chandrasekhar-Milne’s expansion and dotted lines the deformation from 2D structure models. $\overline{\mathrm{\Omega }}$ is the rotation frequency normalised by $\sqrt{GM/R^3}$ where R is the mean radius. $\overline{\mathrm{\Omega }}$ increases from top to bottom in the figure. $\overline{\mathrm{\Omega }}=0.7833$ corresponds to $\frac{\mathrm{\Omega }}{{\mathrm{\Omega }}_{\mathrm{c}}}=0.9$ at the equator as observed for HD 49330.

2.2. Tohoku models of κ-driven pulsations

To model the pulsations, we use the Tohoku non-adiabatic pulsation code. This code includes the effects of the centrifugal and of the Coriolis accelerations. With this code, we perform a non-radial pulsation stability analysis for HD 49330 models. This analysis allows us to predict the modes that will be excited by the κ-mechanism as well as their frequencies and growth rates. More details of this oscillation code are described in Neiner et al. (2012b).

We first perform a non-adiabatic pulsation analysis for rapidly rotating star models whose evolutionary track pass close to the position of HD 49330 in the HR diagram. The evolutionary models are computed in the same manner as in Walker et al. (2005) and the stability analysis is based on the method by Lee & Baraffe (1995). We use the two stellar structures presented above as an input to the models.

We adopt the convention that a negative or positive m order represents a prograde or retrograde mode, respectively. Modes with even or odd parity, which corresponds to even or odd values of l + |m|, respectively possess either a symmetry or an antisymmetry about the equator of the star.

In a rapidly rotating star, a mode cannot be represented by a single l value for a given m order because of the latitudinal couplings between the various spherical harmonics owing to the Coriolis and the centrifugal accelerations. Thus, each mode must be expanded into a series of spherical harmonics. Consequently, for each mode we choose the component of this series with the highest amplitude at the stellar surface and use its l value to compare with the l value determined from the observations.

It is not our goal with these models to reproduce each individual observed frequency in HD 49330 but, rather, to reproduce the ranges of frequencies excited in the star.

We compute non-adiabatic models using a stellar structure deformation based either on Chandrasekhar-Milne’s expansion or on a 2D stellar structure, as presented in Sect. 2.1.

Using Chandrasekhar-Milnes’s expansion, our models reproduce g modes in the observed frequency range only for models with masses significantly lower than the mass determined from spectroscopic observations. See examples in Fig. 2 for a 10 and 13 M_⊙ star: the bottom panel shows the frequencies observed in the CoRoT data during or outside of the outburst detected by Huat et al. (2009) and which have a signal-to-noise ratio (S/N) above 4, the middle panel shows the growth rate of the modelled modes, the top panel shows the azimuthal order m with colors indicating the even and odd parity modes. The typical growth rates of excited modes are 10⁻⁴, meaning that it would take 10⁴ times the oscillation period for the amplitude to grow. g modes are excited in the observed frequency range for the lower mass model but not for M = 13 M_⊙. On the contrary, p modes are excited for the higher mass model but not at lower masses. It is thus difficult to excite both p and g modes at the same time whatever the set of stellar parameters with the Chandrasekhar-Milnes’s expansion.

Fig. 2. Statistically significant observed frequencies of HD 49330 (bottom panel) during and outside of the outburst phase compared to the modelled frequencies in a nearly critically rotating B star (middle panel) of 10 M⊙ (left) and 13 M⊙ (right) calculated using rotational deformation based on Chandrasekhar-Milne’s expansion. Top panel: azimuthal order m of each frequency peak. The dashed line represents the approximate separation between p and g modes. Red full lines correspond to even modes and dashed blue lines to odd modes.

Using rotational deformation based on 2D structure models, both p and g modes can be excited at the same time. Indeed, g modes seem rather insensitive to the way we treat stellar rotational deformation. The excitation of high-frequency p modes, however, seems to occur only when the rotation is nearly critical and thus the treatment of the rotational deformation is very important. See the left panel of Fig. 3 for an example of a star of 10 M_⊙ with both p and g modes excited, which is to be compared with Fig. 2. In this model, f = 11.9 c.d⁻¹ can be identified with a retrograde p mode with m = 1 or 2, while the frequency groups around 1.5 and 3 c d⁻¹ are identified with prograde g modes of m = −1 and −2 respectively. However, the stellar parameters used here do not correspond to the ones derived from spectroscopy (Floquet et al. 2009).

Fig. 3. Same as Fig. 2 but the modelled frequencies are calculated using a 2D structure model.

With stellar parameters compatible with the error boxes derived from spectroscopic observations (see Sect. 1), and with rotational deformation extracted from 2D structure models, p modes are easily excited but g modes are less easy to excite. The right panel of Fig. 3 shows that already for a 13 M_⊙ (and higher), the group of observed g modes around 1.5 c d⁻¹ is not present anymore in the model. In addition, the model predicts several strong modes between 6 and 10 c d⁻¹ which are not observed in the CoRoT data.

Using our modelling results, we can plot colour maps of pulsation amplitudes of modes observed in HD 49330, assuming that they are driven by the κ mechanism. Plotted colour quantities are the dimensionless pressure fluctuations normalised by gravity ($\frac{p\prime }{\sqrt{\rho }rg}$). The y and x axes respectively correspond to the distance from the centre of the star along the rotation axis and perpendicular to it that have been normalised by the mean radius of the star.

We find that g modes are trapped at the equator of the star (see Fig. 4) while p modes are more evenly distributed (see Fig. 5). The trapping of κ-driven g modes is due to rapid rotation.

Fig. 4. Colour maps of pulsation amplitudes for the g modes with l = 2, m = −2 and f = 2.9 c.d−1 (left) and l = 1, m = −1 and f = 1.5 c.d−1 (right), the two main frequencies observed in the CoRoT data during the outburst of HD 49330 (Huat et al. 2009).

Fig. 5. Colour maps of pulsation amplitudes for the p modes with l = 2, m = 2 and f = 11.7 c.d−1 (left) and l = 2, m = 2 and f = 16.9 c.d−1 (right), the two frequencies observed both in the CoRoT data in the quiet phase of HD 49330 (Huat et al. 2009) and in spectroscopic data (Floquet et al. 2009).

2.3. Issues with the κ mechanism

From the models presented above we conclude that (1) non-adiabatic treatment is needed, in particular to study the modes below 2Ω and to know whether the modes are excited; and (2) a 2D structure model is required to be able to reproduce p modes for which the rotational deformation has a strong impact. Nevertheless, whatever model of stellar structure and pulsations we use, we cannot excite p and g modes at the same time with the spectroscopically-derived stellar parameters of HD 49330 if we consider only pulsation modes excited by the κ mechanism.

Huat et al. (2009) and Floquet et al. (2009) showed that p modes are visible in the CoRoT data of HD 49330 when the star is in a quiet or relaxation phase, in other words outside the observed outburst. This is when the material at the surface of the star is close to but not at its critical velocity. During the precursor and outburst phase, g modes appear in the CoRoT data. This is when the velocity of the material at the surface of the star at its equator becomes critical and material gets ejected from the star.

Therefore several possibilities can be considered to try to answer the problem of excitation of both p and g modes:

First, the modelling of the observed g modes, assuming they are excited by the κ mechanism, requires a set of parameters with a much lower mass than the spectroscopically determined one. During the outburst, the star ejects material and its mass thus decreases. However, the mass ejected from the star during the outburst is much lower than the difference between the mass determined from observations and required in the model. Consequently, the mass loss during the outburst cannot explain the excitation of g modes.

Second, a change in the radius of the star would influence the exact frequency of the pulsations. Huat et al. (2009), however, did not observe changes in the frequencies along the CoRoT light curve. They only witnessed changes in amplitudes of the frequency peaks. Therefore, the change in radius during the outburst cannot be the cause of the excitation of g modes.

Finally, the outburst lasted 73 days (including the precursor phase), which is much shorter than the typical growth rate of κ-driven g modes of 10⁴ times the oscillation period (see e.g. the middle panel of Fig. 3). This indicates either that the observed g modes were excited already before the precursor phase but not visible at the surface of the star, or that these modes should not be considered as global oscillations of the star.

3. Alternative scenario: κ-driven p modes and stochastically-excited g waves

3.1. Proposed scenario

The results of our modelling show that p and g modes cannot be excited by the κ mechanism at the same time in HD 49330 using the spectroscopically-determined stellar parameters. Moreover, the power spectrum of Huat et al. (2009) shows broad frequency groups rather than sharp frequency peaks around 1–2 c d⁻¹ (see e.g. the bottom panel of Fig. 3). Considering the above arguments, we propose that what Huat et al. (2009) and Floquet et al. (2009) attributed to κ-driven g modes are rather stochastic g waves.

During the quiet phase, stochastic gravity waves can be excited in the convective core (Browning et al. 2004; Rogers et al. 2013; Augustson et al. 2016). These waves transport angular momentum from the core to the surface (Lee et al. 2014; Rogers 2015). When enough angular momentum has accumulated in the outer layers of the star, these layers become unstable and emit transient gravity waves. The surface layers reach the critical velocity, in particular at the equator where the rotation was already the closest to critical. The destabilisation of the surface layers thus ignites the outburst and breaks the cavity where the p modes were propagating. This explains both the disappearance of the p modes during the precursor and outburst phase of HD 49330 and the ejection of material from the surface into the disk, meaning the occurrence of the outburst, as observed by Huat et al. (2009) with CoRoT. Relaxation then occurs, recreating the cavity and letting the p modes reappear while the surface transient g waves disappear.

This outburst phenomenon will reoccur whenever the g modes stochastically excited in the core have transported a sufficient amount of angular momentum to the outer layers of the star. Both the stochastically excited g modes from the core and the transient g waves from the outer layers could be observed if their amplitude is high enough. A sketch of this scenario is presented in Fig. 6.

Fig. 6. Sketch of the Be outburst of HD 49330. In the quiescent phase, p modes are present and g modes stochastically excited at the interface between the convective core and the radiative envelope transport momentum towards the surface. Once enough momentum has accumulated, transient g modes get excited in the surface layer while p modes are destabilised and thus disappear. An ejection of material finally occurs. The star then relaxes back to a quiescent phase and p modes reappear.

Below we test this scenario by modelling stochastically excited gravito-inertial modes with the non-adiabatic Tohoku code modified for this purpose.

3.2. Stochastic excitation of g waves

3.2.1. Theoretical background

Convective regions, such as convective cores and subsurface convection zones in massive stars, are able to stochastically excite oscillation modes (Cantiello et al. 2009; Belkacem et al. 2010) and particularly gravity modes (Samadi et al. 2010). The latter become gravito-inertial modes in fast rotators such as HD 49330 because of the action of the Coriolis acceleration (see e.g. Lee & Saio 1997; Dintrans & Rieutord 2000; Mathis 2009; Ballot et al. 2010). This excitation is now observed both in realistic numerical simulations of convective cores surrounded by a stably stratified radiative envelope (see Browning et al. 2004; Rogers et al. 2013; Augustson et al. 2016; Edelmann et al. 2019) and in asteroseismic data (Neiner et al. 2012a; Aerts & Rogers 2015). Gravito-inertial waves are excited through their couplings with the volumetric turbulence in the bulk of convective regions (where purely gravity modes in a slowly rotating star are evanescent while gravito-inertial modes in a rapidly rotating star become inertial) and by the impact of structured turbulent plumes at the interfaces between convective and radiative regions. As a first step, we choose to focus on the volumetric stochastic excitation.

In this context, Belkacem et al. (2009) derived the general formalism to treat the stochastic excitation of non-radial modes in rotating stars taking into account the Coriolis acceleration. Their method is based on the results first obtained by Samadi & Goupil (2001) for radial p modes, which have been extended to non-radial modes in a non-rotating star by Belkacem et al. (2008a). However, this formalism has been applied only to slowly rotating stars. Therefore we revisit it in the case of rapidly rotating stars, still ignoring centrifugal acceleration since Ballot et al. (2010) demonstrated that it can be neglected for g modes given the relatively weak distortion of the excitation region. The details of our new equations are provided in Appendix A. The main differences with Belkacem et al. (2009) are related to the strong action of the Coriolis acceleration. First, the expansion of the wave displacement now involves several spherical harmonics (see Eq. (A.7)). Next, the Coriolis acceleration modifies the orthogonality relationships for eigenmodes (Schenk et al. 2002). The normalisation of the amplitude becomes therefore more complex (see Eq. (A.10)). In this work, we still consider a simplified model for the convective turbulent source term, where the turbulence is assumed to be statistically stationary, incompressible, homogeneous, and isotropic, as in Samadi et al. (2010). These assumptions are assumed as a first step because of the lack of available robust prescriptions for convective rotating turbulence. A better model of turbulence would require specific developments which are beyond the scope of this work (we refer the reader to the Appendix A.3 and Mathis et al. 2014 for a detailed discussion).

3.2.2. Tohoku models with stochastically excited g modes

We calculated seismic models of a M = 15 M_⊙ mass star with X = 0.7 and Z = 0.02 using the formalism described above and in Appendix A for $\overline{\mathrm{\Omega }}=0.1$, 0.3, 0.4, 0.5, 0.6, and 0.8. We use series expansion for the perturbations of rotating stars. For the stochastic excitation of the low frequency modes, we compute ∇ − ∇_ad for the core using the usual mixing length theory with physical quantities of a ZAMS model. It is the first time that such models are computed for a rapidly rotating hot star.

We find that the stochastic excitation of g and r modes occurs in this 15 M_⊙ star model and thus the mass of HD 49330 (M = 14.4 M_⊙) is not a concern anymore to produce the observed frequencies with this excitation mechanism (contrary to κ-driven modes).

Figure 7 shows the variation in stellar surface magnitude produced by these waves as a function of their frequency for $\overline{\mathrm{\Omega }}=0.1$, 0.4, 0.6, and 0.8. Since the star is considered spherical here, $\overline{\mathrm{\Omega }}=0.1$, 0.4, 0.6, and 0.8 would correspond to 0.12, 0.46, 0.69, and 0.92 at the equator if the star was deformed as in the 2D ROTORC model shown above. The amplitude is computed from the maximum value of the luminosity variation at the surface induced by the stochastically excited waves. We refer the reader to Appendix A of Lee & Saio (1993) for the equations solved for non-adiabatic oscillations that allow the computation of the luminosity variation in the code.

Fig. 7. Amplitude of the stochastically excited g and r modes for various rotation rates: $\overline{\mathrm{\Omega }}=0.1$ (left), $\overline{\mathrm{\Omega }}=0.4$ (middle left), $\overline{\mathrm{\Omega }}=0.6$ (middle right), and $\overline{\mathrm{\Omega }}=0.8$ (right). The black lines represent l = |m| = 1, red and green lines represent l = 2 with |m| = 2 and |m| = 1 respectively, and blue lines represent l = 3, |m| = 2 modes. The solid and dotted lines are respectively for prograde and retrograde modes. The left panel has a different y-axis range to better see the low amplitudes.

For the figure, since for rapidly rotating stars mode identification becomes unclear particularly for retrograde modes, we consider non-adiabatic pulsations for which the kinetic energy l = |m| even parity and l = |m|+1 odd parity modes is above 50% of the total kinetic energy of the sum of all components constituting an eigenmode for a given m. We include g and r modes only for |m| = 1 and 2 with l = |m| and |m|+1. We find that the amplitude of the even parity, prograde modes is the strongest.

The computed amplitude for $\overline{\mathrm{\Omega }}\ge 0.4$ is similar to the amplitudes detected with CoRoT, which confirms that we likely see stochastically excited gravito-inertial modes in HD 49330. In fact, the modelled amplitude is slightly higher than the amplitude observed for HD 49330 (see Fig. 8, which shows in the bottom panel only the detected frequencies considered statistically significant, that is with a S/N above 4, while many more lower amplitude frequencies are present in the data). This may be because we considered a M = 15 M_⊙ model, while a M = 14.4 M_⊙ star such as HD 49330 would have a lower luminosity. In addition, our model is not complete since only the action of Reynolds stresses associated to convective small scales has been considered while large scale convective structures such as turbulent plumes likely also contribute to the stochastic excitation of gravito-inertial waves (Schatzman 1993; Rogers et al. 2013; Alvan et al. 2014; Pinçon et al. 2016), and the assumed turbulence source term is very simple (see Appendix A.3).

Fig. 8. Top: same as Fig. 7 for $\overline{\mathrm{\Omega }}=0.8$. Bottom: statistically significant observed frequencies of HD 49330 during and outside of the outburst phase as detected by Huat et al. (2009).

Nevertheless, the models show that an increasing rotation rate modifies the amplitude, allowing certain frequencies to be excited with higher amplitudes for rapid rotators. This is the result of a better coupling between gravito-inertial waves and the convective forcing, in particular in the sub-inertial regime in which gravito-inertial waves become propagative inertial waves in convective regions. Such a behaviour has been observed in numerical simulations (Rogers et al. 2013) and predicted theoretically (Mathis et al. 2014).

However, the variation in the amplitude of some of the modes in our models is not a monotonic function of the rotation rate. For example, the amplitude peaks for $\overline{\mathrm{\Omega }}=0.4$ and 0.8 for l = |m| = 1 modes, while it increases with rotation rates for l = |m| = 2 modes. This is because both the damping rate of modes and the injection of energy are not monotonic. Moreover, above $\overline{\mathrm{\Omega }}=0.4$, the fact that we ignored the deformation from a spherical star owing to centrifugal acceleration in the computation of waves (Ballot et al. 2010; Mathis & Prat 2019) as well as the increase in the size of the core observed in rapidly rotating Be stars (Neiner et al. 2012b) may have more impact.

Finally, as expected, we find that more frequencies are excited when the star rotates rapidly (Rieutord 2009), and most g modes are in the sub-inertial regime. This forest of stochastically excited g modes with higher amplitudes on average could help to transport angular momentum in rapidly rotating stars.

3.3. Transport of angular momentum

The aim of this subsection is to test the correlation between the transport of angular momentum by gravito-inertial waves, which are stochastically excited inside the star, and Be outbursts.

First, we recall that the idea that non-radial oscillations excited in massive stars are able to efficiently transport angular momentum through their damping by thermal diffusion below the stellar surface and that they may impel the surface of a Be star to reach its critical velocity has already been proposed by Ando (1986), Lee & Saio (1993), and Lee (2006). However, in these previous works, gravito-inertial waves (and r modes) were excited by the κ-mechanism. In the case of HD 49330, we propose that gravito-inertial waves are stochastically excited in turbulent convective regions. Pantillon et al. (2007) and Lee et al. (2014) demonstrated that such waves are also able to transport angular momentum.

Following the method of Lee (2012; see also Grimshaw 1984; Lee 2006; Mathis & de Brye 2012), we consequently study the transport of angular momentum induced by stochastically excited gravito-inertial waves. We recall that here we ignore the centrifugal acceleration (i.e. the star is considered spherical) both for the stellar structure model and the pulsations computations (contrary to the κ-driven modes computed in Sect. 2). Following this method, we calculate the acceleration of the mean azimuthal flow on an isobar. We consider the horizontally averaged equation describing the interaction between the wave and the mean flow, given by:

$$\begin{array}{c}\hfill \overline{\rho }\frac{\mathrm{d}\ell }{\mathrm{d}t}=-{\displaystyle \sum _{\alpha }\{\frac{1}{r^2}{\partial }_r\left(r^2\overline{F_{\alpha }^{\mathrm{AM};\mathrm{V}}}\right)+\overline{\langle {\rho }_{\alpha }\prime {\partial }_{\phi }{\mathrm{\Phi }}_{\alpha }\prime \rangle }\},}\end{array}$$(1)

where the subscript α identifies the excited mode, $\langle \dots \rangle =\frac{1}{2\pi }{\int }_0^{2\pi }\dots \mathrm{d}\phi$ and $\overline{\dots }=\frac{1}{\pi }{\int }_0^{\pi }\dots sin\theta \mathrm{d}\theta$. Here, $\ell =\overline{\langle r^2{sin}^2\theta \mathrm{\Omega }+rsin\theta {\overline{v}}_{\phi }\rangle }$ is the specific angular momentum around the rotation axis, where Ω is the considered uniform rotation, and $\overline{v_{\phi }}/(rsin\theta )$ is the induced differential rotation. Each scalar field (X) has been expanded as $X=\overline{X}+X\prime$, where $\overline{X}$ is the mean value on an isobar and X′ is the modal fluctuation. The vertical flux of angular momentum is given by

$$\begin{array}{c}\hfill F_{\alpha }^{\mathrm{AM};\mathrm{V}}=\overline{\rho }rsin\theta \langle v_r\prime (v_{\phi }\prime +2\mathrm{\Omega }cos\theta {\xi }_{\theta })\rangle ,\end{array}$$(2)

where we recall that, for modes calculated taking only the given uniform rotation into account,v = i(ω₀+mΩ)ξ.

We define the associated timescale for the transport of angular momentum τ as:

$$\begin{array}{c}\hfill \frac{1}{\tau }=\frac{\mathrm{d}ln\ell }{\mathrm{d}t}·\end{array}$$(3)

Figure 9 shows the inverse of the timescale of stochastically excited prograde and retrograde g and r modes for slowly and rapidly rotating stars with M = 15 M_⊙. For a sufficiently rapid rotation rate, angular momentum is transported just below the surface. As expected, retrograde modes extract angular momentum, whereas prograde modes deposit angular momentum at the surface, leading respectively to a deceleration and acceleration of the surface rotation. Below the surface (at about r = 0.975 R), the modes stochastically excited in the core experience an additional κ excitation at the position of the opacity bump, which creates the reversal of sign in the angular momentum flux (see Fig. 9). In the case of a slowly rotating star, the net transport (sum of deposition and extraction) is very low. In the case of a rapidly rotating star, however, most of the retrograde modes are r modes, for which the pressure and density perturbations are small so that their angular momentum transfer is inefficient. Therefore, the prograde modes deposit much more momentum than is extracted by the retrograde modes at the surface. Thus, the net angular momentum transport by stochastically excited modes is large and clearly increases with the rotation rate at the surface. These modelling results confirm the prediction by Mathis et al. (2014) and Lee et al. (2014).

Fig. 9. Inverse of the timescale of retrograde (red line), prograde (blue line) and the sum of the two (black line), g and r modes for a M = 15 M⊙ star with $\overline{\mathrm{\Omega }}=0.1$ (left), 0.4 (middle left), 0.6 (middle right) and 0.8 (right), where τ is normalised to $\sqrt{GM/R^3}$. The value of 1/τ is proportional to the angular momentum transport. Note that the right panel (fastest rotator) has a different y-axis range to accommodate the larger 1/τ values and shows an extra y-axis on the right side in years−1.

In addition, the minimum timescale τ is attained around r = 0.975 R and is of the order of 8.3 × 10⁵ (1/τ = 1.2 × 10⁻⁶) for the sum of prograde and retrograde modes in a rapidly rotating star. Since Ω_c is about 1.5 × 10⁻⁴ s⁻¹ in this model, the acceleration timescale is of the order of 125 years. This timescale would correspond to the time needed for the first outburst to occur, that is for the required angular momentum to be transported from the core to the surface and accumulated there.

A characteristic recurrence timescale of outbursts defined as

$$\begin{array}{c}\hfill \frac{1}{\tau \prime }=\frac{\mathrm{d}ln\overline{\langle r^2{sin}^2\theta {\mathrm{\Omega }}_{\mathrm{c}}-(r^2{sin}^2\theta \mathrm{\Omega }+rsin\theta {\overline{v}}_{\phi })\rangle }}{\mathrm{d}t}\end{array}$$(4)

can also be computed. In the case of a shellular rotation that depends only on the radial coordinate, it becomes

$$\begin{array}{c}\hfill \tau \prime =\tau |{\overline{\mathrm{\Omega }}}^{-1}-1|.\end{array}$$(5)

It evaluates how much time is needed to reach again critical velocity at the surface after a first outburst, assuming that the transport triggered by waves occurs continuously; if $\overline{\mathrm{\Omega }}\to 1$ then τ′→0, meaning if the star rotates at critical velocity it looses mass all the time, while if $\overline{\mathrm{\Omega }}\to 0$ then τ′→∞, meaning if a star does not rotate no outburst occurs. $\overline{\mathrm{\Omega }}=0.8$ means $\overline{\mathrm{\Omega }}=0.92$ in the non-spherical case. We then predict that an outburst should occur every ∼11 years in HD 49330. Photometric observations by ASAS-3 have shown that an outburst occurred around December 2002, while the outburst observed by CoRoT occurred in February 2008 (Huat et al. 2009). The time difference between these two outbursts is about half but still similar to the recurrence timescale computed here. However, a more regular monitoring of the lightcurve of HD 49330 is needed to confirm the observed recurrence of its outbursts.

3.4. Outburst

The accumulation of angular momentum near the surface destabilises the outer layers of the star. This situation is similar to local forces and heating or cooling, such as a turbulent layer, wind, or convective structure in the terrestrial atmosphere, which generate gravity waves (e.g. Vadas & Fritts 2001; Townsend 1965). Realistic models for such excitation have been developed in atmospheric physics but never applied in the context of hot stars. Nevertheless, by analogy, we expect the perturbation excited in the outer layers of HD 49330 to produce local mechanical and thermal forcing that generate transient g waves. Moreover, once the amount of angular momentum deposited in the surface layers of the star impels it to the breakup velocity, the star starts loosing mass, meaning the outburst occurs.

Both the g modes from the core and the g waves from the outer layers could be observable. In the CoRoT observations of HD 49330, g waves are mostly visible during the precursor and outburst phases, therefore they must be the transient g waves from the outer layers. Nevertheless, some peaks around 3 and 3.6 c d⁻¹ seem to also be present during quieter phases (see Fig. 4 in Huat et al. 2009) and could originate in the core.

4. Conclusions

4.1. Summary of the results

Here, we summarise the main findings of this work:

1. The treatment of rotational deformation is very important for high-frequency p modes, because the excitation of these modes seems to occur only when the rotation is nearly critical. Consequently, Chandrasekhar-Milne’s expansion is not a good approximation to model p modes of rapidly rotating early B stars, such as early Be stars or rapidly rotating β Cephei stars. Using a 2D structure model is essential in this case.

2. Contrary to what is usually assumed, Be stars seem to host not only κ-driven pulsation modes, but also stochastically excited g waves produced continuously in the convective core as well as in a transient way in the outer layers of the star during an outburst. A direct detection of stochastically excited gravito-inertial modes has been obtained by Neiner et al. (2012a) in the Be star HD 51452. This should be carefully taken into account for future interpretations of seismic data of Be stars. Indeed, while for early Be stars, it is clear that κ-driven g modes cannot be excited and thus that any observed g mode is stochastically excited (Neiner et al. 2012a), for late Be stars g modes or waves can be excited both by the κ mechanism and stochastically.

3. The amplitude of stochastically excited modes increases with rotation rate on average. However, we find that the amplitude of some modes does not vary monotonically with $\overline{\mathrm{\Omega }}$. This is likely related to the competition between the injection of energy and damping rates, which are not monotonic functions of the rotation rate. The l = 2, m = −2 modes, however, have a clear increase in magnitude with rotation rate and contribute significantly to the transport of angular momentum. Those modes are the ones most commonly observed in Be stars (Rivinius et al. 2003).

4. The transport of angular momentum by stochastically excited modes clearly increases with rotation rate and becomes very strong for stars rotating close to critical. This is because prograde modes dominate retrograde modes. In particular, the l = 2, m = −2 modes have a strong amplitude when the star rotates fast.

5. The recurrence timescale of outbursts produced by this angular momentum transport mechanism is of the order of 11 years for HD 49330, which is similar to the outbursts occurrence observed in this star.

6. The outburst in HD 49330 is due, in addition to the almost critical velocity of the star, to the transport of angular momentum to the surface layers by stochastic gravity waves generated in the convective core. The accumulation of angular momentum in the surface layers is at the origin of the destabilisation of the surface layers, of the production of transient g waves in these layers, and of the outburst. This explains the correlation between the variation in the amplitude of the pulsation frequencies observed in the CoRoT data and the occurrence of the outburst. A similar scenario could occur in other pulsating Be stars, thus providing an explanation for the long-standing question of the origin of Be outbursts and disks, at least for some Be stars.

4.2. Limitations and future work

This work is a first attempt at coherently modelling stochastically excited gravito-inertial pulsation modes, and the transport of angular momentum they induce, in a rapidly rotating Be star. The formalism we developed here for Be stars can be applied to other rapidly rotating g-mode pulsators, such as γ Dor or SPB stars.

However, this first attempt can be improved in the future. In particular, the model presented here for stochastically excited gravito-inertial modes does not include centrifugal acceleration, neither for the stellar structure model nor for the pulsation computation, contrary to the models of κ-driven modes that we present. Rotation is nevertheless partly taken into account via the Coriolis acceleration for both the stellar structure and oscillations. In the future, it would be interesting to use a 2D stellar structure to see how this may change the amplitude of the modes. Moreover, Neiner et al. (2012b) showed that mixing induced by the penetrative movements at the bottom of the radiative envelope and by the secular hydrodynamical transport processes induced by the rotation in the envelope extend the size of the core of Be stars, which may increase stochastic excitation.

In addition, we used a standard description of turbulent convection in the core, which does not account for rapid rotation. Indeed, for a rapidly rotating star, convection becomes strongly anisotropic and the transfer of energy may be modified by rapid rotation. However, no robust prescription for these effects are available as of this writing. Furthermore, the formalism derived here for stochastic excitation focuses on the action of Reynolds stresses associated to convective small scales while large scale convective structures such as turbulent plumes likely also contribute to the stochastic excitation of gravito-inertial waves.

Finally, we considered κ-driven p and g modes as well as stochastically excited gravito-inertial (g and r) modes. Saio et al. (2018) and Saio (2018) showed that r modes could also be excited mechanically in most rotating stars, including Be stars. In this case, even r modes of order m appear as a group of frequencies just below m times the rotation frequency in the Fourier spectrum of those stars. Such mechanically excited r modes can be a transient phenomenon at the surface and could therefore be an alternative explanation to the low frequencies observed during the outburst of HD 49330, which will be studied in a separate paper (Saio et al., in prep.).

Acknowledgments

The authors thank the referee for their comments, which have allowed us to improve the article. The CoRoT space mission, launched on December 27th 2006, has been developed and was operated by CNES, with the contribution of Austria, Belgium, Brazil, ESA (RSSD and Science Program), Germany and Spain. CN acknowledges fundings from the SIROCO ANR project, CNES, and PNPS. SM and KCA acknowledge support from the ERC SPIRE 647383 grant and PLATO CNES grant at CEA/DAp-AIM. This research has made use of the SIMBAD database operated at CDS, Strasbourg (France), and of NASA’s Astrophysics Data System (ADS).

References

Appendix A: Formalism for stochastically excited g modes modified by rapid rotation

In this appendix, we describe the equations derived and implemented in the Tohoku code for the treatment of stochastically excited gravity modes modified by rapid rotation.

A.1. Stochastic excitation

The inhomogeneous wave equation taking the Coriolis acceleration into account and assuming a solid body rotation can be written

$$\begin{array}{c}\hfill (\frac{{\partial }^2}{\partial t^2}-{\mathit{L}}_{\mathrm{\Omega }}){\mathit{v}}_{\mathrm{osc}}+{\mathcal{C}}_{\mathrm{osc}}={\mathcal{S}}_t.\end{array}$$(A.1)

All physical quantities are split into an equilibrium term and a perturbation. In the following equations, the subscripts 1 and 0 denote Eulerian perturbations and equilibrium quantities respectively, except for the velocity where the subscript 1 has been dropped to ease the notation. The latter is split into two contributions, namely the oscillation velocity (v_osc) and the turbulent velocity field (u_t), such that u = v_osc + u_t. L_Ω is the linear oscillation operator which, in presence of uniform rotation, becomes

$$\begin{array}{c}\hfill {\mathit{L}}_{\mathrm{\Omega }}=\mathbf{\nabla }[{\alpha }_s{\mathit{v}}_{\mathrm{osc}}·\mathbf{\nabla }{s}_0+c_s^2\mathbf{\nabla }·\left({\rho }_0{\mathit{v}}_{\mathrm{osc}}\right)]-{\mathit{g}}_{\mathrm{eff}}\mathbf{\nabla }·\left({\rho }_0{\mathit{v}}_{\mathrm{osc}}\right)\\ \hfill & -{\rho }_0\mathrm{\Omega }\frac{{\partial }^2{\mathit{v}}_{\mathrm{osc}}}{\partial t\partial \phi }-2{\rho }_0\mathbf{\Omega }\times \frac{\partial {\mathit{v}}_{\mathrm{osc}}}{\partial t},\hfill \end{array}$$(A.2)

where c_s is the sound velocity and α_s = (∂p/∂s)_ρ. Next, the operator

$$\begin{array}{cc}\hfill {\mathcal{C}}_{\mathrm{osc}}& =\frac{\partial }{\partial t}[\frac{\partial \left({\rho }_1{\mathit{v}}_{\mathrm{osc}}\right)}{\partial t}+2\mathbf{\nabla }:\left({\rho }_0{\mathit{v}}_{\mathrm{osc}}{\mathit{u}}_t\right)\hfill \\ \hfill & +{\rho }_1\mathrm{\Omega }\frac{\partial {\mathit{v}}_{\mathrm{osc}}}{\partial \phi }+2{\rho }_1\mathbf{\Omega }\times {\mathit{v}}_{\mathrm{osc}}\hfill \\ \hfill & +\mathbf{\nabla }({\alpha }_s{\mathit{v}}_{\mathrm{osc}}·\mathbf{\nabla }{s}_1+c_s^2\mathbf{\nabla }·({\rho }_1{\mathit{v}}_{\mathrm{osc}}))]\hfill \end{array}$$(A.3)

involves both turbulent and pulsational velocities and contributes to the linear dynamical damping (see Samadi & Goupil 2001, for details). Finally, the 𝒮_t operator that contains the source terms of the inhomogeneous wave equation (Eq. (A.1)) is given by

$$\begin{array}{cc}\hfill {\mathcal{S}}_t& {\displaystyle =-\frac{\partial }{\partial t}\mathbf{\nabla }:({\rho }_0{\mathit{u}}_t{\mathit{u}}_t)+\mathbf{\nabla }({\alpha }_s{\mathit{u}}_t·\mathbf{\nabla }{s}_1)+{\mathcal{S}}_{\mathrm{\Omega }}+{\mathcal{S}}_M}\hfill \end{array}$$(A.4)

with

$$\begin{array}{cc}\hfill {\mathcal{S}}_{\mathrm{\Omega }}& =-\frac{\partial }{\partial t}\left[{\rho }_1(\mathrm{\Omega }\frac{\partial }{\partial \phi }{\mathit{u}}_t-2\mathbf{\Omega }\times {\mathit{u}}_t)\right],\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\mathcal{S}}_M& =\frac{\partial }{\partial t}({\rho }_1{\mathbf{g}}_1)+\mathbf{\nabla }[c_s^2\mathbf{\nabla }·({\rho }_1{\mathit{u}}_t)]-{\mathit{g}}_{\mathrm{eff}}\mathbf{\nabla }·({\rho }_1{\mathit{u}}_t)\hfill \\ \hfill & -\frac{{\partial }^2}{\partial t^2}({\rho }_1{\mathit{u}}_t)+{\mathcal{L}}_t,\hfill \end{array}$$(A.6)

where ℒ_t contains the linear terms source¹. The first two terms of Eq. (A.4) correspond to the Reynolds stress and entropy contributions, respectively. The rotational contributions are the two following terms in Eq. (A.4), that are respectively related to the Doppler effect and the Coriolis acceleration. Finally, the next four terms are neglected since they scale as ℳ³ where ℳ is the Mach number (see Samadi & Goupil 2001, for details). Therefore, those terms are negligible in comparison to the Reynolds stress contribution that scales as ℳ². Note that in the case of “free” turbulence where the oscillations do not have any effect on the turbulent flow, we have ρ₁ ≈ ρ_t.

A.2. Mean square amplitude

Using Eq. (A.1), we can determine the mean square amplitude of modes excited stochastically.

First, we recall that in a rapidly rotating star, a single l value does no longer represent a pulsation mode for a given m value because of the latitudinal couplings between the different spherical harmonics owing to the Coriolis and the centrifugal accelerations. For this reason the eigenfunctions for a fixed m are expanded into a series of terms proportional to spherical harmonics Y_l, m(θ,φ) with l ≥ |m|. We express the spatial dependence of the displacement vector ξ as

$$\begin{array}{c}\hfill {\displaystyle \mathit{\xi }={\displaystyle \sum _{j=1}^J\{{\xi }_{r;m}^{l_j}{Y}_{l_j}^m{\widehat{\mathbf{e}}}_r+{\xi }_{H;m}^{l_j}{\mathbf{\nabla }}_H{Y}_{l_j}^m+{\xi }_{T;m}^{l_{j\prime }}[{\mathbf{\nabla }}_H{Y}_{l_{j\prime }}^m\times {\widehat{\mathbf{e}}}_r]\},}}\end{array}$$(A.7)

where ${\displaystyle {\mathbf{\nabla }}_H=\frac{\partial }{\partial \theta }{\widehat{\mathbf{e}}}_{\theta }+\frac{1}{sin\theta }\frac{\partial }{\partial \phi }{\widehat{\mathbf{e}}}_{\phi }}$ and (r,θ,φ) are the usual spherical coordinates, l_j = |m|+2(j−1) + I and l_j′ = l_j + 1 − 2I with I = 0 for even modes and I = 1 for odd ones. The terms proportional to ${\xi }_{T;m}^{l_{j\prime }}$ represent a toroidal displacement needed because of the Coriolis acceleration term in the momentum equation.

In addition, the wave velocity field (v) is related to the displacement (ξ) by

$$\begin{array}{c}\hfill \mathit{v}=Ai({\omega }_0+m\mathrm{\Omega })\mathit{\xi }{e}^{i{\omega }_{0}t},\end{array}$$(A.8)

where ω₀ is the eigenfrequency and A(t) is the amplitude due to the turbulent forcing.

Following Belkacem et al. (2009) and generalising the results to the case of rapid rotation, one then obtains the mean square amplitude

$$\begin{array}{c}\hfill {⟨|A|}^2⟩=\frac{1}{8\eta {\omega }_0^2}{\left|\frac{1}{b}\right|}^2(C_R^2+C_S^2+C_{\mathrm{\Omega }}^2+C_c^2),\end{array}$$(A.9)

where the operator ⟨…⟩ denotes a statistical average performed on an infinite number of independent realisations, η is the damping rate computed by the non-adiabatic oscillation code in the radiative flux approximation (i.e. δ(∇⋅F_conv) = 0 where F_conv is the convective flux; see e.g. Lee 2006) as expected in massive stars (see the discussion in Belkacem et al. 2010; Samadi et al. 2010).

We introduce the modified mode inertia in the case of rapid rotation (Schenk et al. 2002):

$$\begin{array}{c}\hfill b={\displaystyle \sum _j{\int }_0^{R_{\mathrm{CZ}}}{I}_b^j\left(r\right){\rho }_0\left(r\right)r^2\mathrm{d}r,}\end{array}$$(A.10)

where

$$\begin{array}{cc}\hfill I_b^j\left(r\right)& ={\left({\xi }_{r;m}^{l_j}\right)}^{\ast }[{\xi }_{r;m}^{l_j}-\frac{i\mathrm{\Omega }}{{\omega }_0}(C_m^{l_j}+D_m^{l_j})]\hfill \\ \hfill & +{L_j}^2{\left({\xi }_{H;m}^{l_j}\right)}^{\ast }[{\xi }_{H;m}^{l_j}-\frac{i\mathrm{\Omega }}{{\omega }_0}(E_m^{l_j}+F_m^{l_j})]\hfill \\ \hfill & +{L_{j\prime }}^2{\left({\xi }_{T;m}^{l_j\prime }\right)}^{\ast }[{\xi }_{T;m}^{l_j\prime }-\frac{i\mathrm{\Omega }}{{\omega }_0}(G_m^{l_j\prime }+H_m^{l_j\prime })]\hfill \end{array}$$(A.11)

with

$$\begin{array}{cc}& C_m^{l_j}=\frac{1}{{L_j}^2}[l_j{B}_m^{l_j}{\xi }_{T;m}^{l_j-1}-(l_j+1)B_m^{l_j+1}{\xi }_{T;m}^{l_j+1}]\hfill \\ \hfill & D_m^{l_j}=-im{\xi }_{H;m}^{l_j}\hfill \\ \hfill & E_m^{l_j}=\frac{1}{{L_j}^2}[B_{l_j}^m{\xi }_{T;m}^{l_j-1}+B_m^{l_j+1}{\xi }_{T;m}^{l_j+1}]\hfill \\ \hfill & F_m^{l_j}=-\frac{\mathit{im}}{{L_j}^2}[{\xi }_{r;m}^{l_j}+{\xi }_{H;m}^{l_j}]\hfill \\ \hfill & G_m^{l_j\prime }=-\frac{\mathit{im}}{{L_{j\prime }}^2}{\xi }_{T;m}^{l_{j\prime }}\hfill \\ \hfill & H_m^{l_{j\prime }}=l_{j\prime }{A}_m^{l_{j\prime }}[{\xi }_{r;m}^{l_{j\prime }-1}-(l_{j\prime }-1){\xi }_{H;m}^{l_{j\prime }-1}]\hfill \\ \hfill & -(l_{j\prime }+1)A_m^{l_{j\prime }+1}[{\xi }_{r;m}^{l_{j\prime }+1}+(l_{j\prime }+2){\xi }_{H;m}^{l_{j\prime }+1}]\hfill \end{array}$$

where

$$\begin{array}{cc}& A_m^l=\frac{1}{l^2}\sqrt{\frac{l^2-m^2}{4l^2-1}}\mathrm{and}{B}_m^l=l^2(l^2-1)A_m^l\hfill \end{array}$$

and $L_l^2=l(l+1)$.

Finally, $C_R^2$ is the Reynolds stress contribution, $C_S^2$ is the one related to entropy fluctuations, $C_{\mathrm{\Omega }}^2$ contains the contributions of the Coriolis acceleration and the Doppler term, and C_c represents the cross-source terms, that is the interferences between the various excitation sources.

A.2.1. Reynolds stress contribution

We follow Belkacem et al. (2009) but have modified their R3 term to account for contributions that are important for rapid rotation. As a result, the turbulent Reynolds stress contribution is:

$$\begin{array}{cc}\hfill C_R^2& =4{\pi }^3{\displaystyle {\int }_{{\mathcal{V}}_{\mathrm{CZ}}}\mathrm{d}m\sum _j{R}_m^j(r)S_R({\omega }_0),}\hfill \end{array}$$(A.12)

where 𝒱_CZ is the volume of the convective regions, with

$$\begin{array}{c}\hfill R_m^j(r)=R_{1;m}^j(r)-R_{2;m}^j(r)+R_{3;m}^j(r)\end{array}$$(A.13)

(we refer the reader to Eq. (A.8) in Belkacem et al. 2009). We have:

$$\begin{array}{cc}\hfill R_{1;m}^j(r)& =2{\left|\frac{\mathrm{d}{\xi }_{r;m}^{l_j}}{\mathrm{d}r}\right|}^2+4{\left|\frac{{\xi }_{r;m}^{l_j}}{r}\right|}^2\hfill \\ \hfill & +L_{l_j}^2{\left|{\mathcal{A}}_{R;m}^{l_j}\right|}^2-2L_{l_j}^2[\frac{{\left({\xi }_{r;m}^{l_j}\right)}^{\ast }{\xi }_{H,m}^{l_j}}{r^2}+\mathrm{c}.\mathrm{c}.]\hfill \\ \hfill & +2L_{l_j}^2(L_{l_j}^2-1){\left|\frac{{\xi }_{H;m}^{l_j}}{r}\right|}^2\hfill \\ \hfill & +L_{l_j\prime }^2(L_{l_j\prime }^2-2){\left|\frac{{\xi }_{T;m}^{l_{j\prime }}}{r}\right|}^2\hfill \\ \hfill & +L_{l_j\prime }^2{\left|{\mathcal{B}}_{R;m}^{l_{j\prime }}\right|}^2,\hfill \end{array}$$(A.14)

where

$$\begin{array}{cc}\hfill {\mathcal{A}}_{R;m}^{l_j}& =\frac{\mathrm{d}{\xi }_{H;m}^{l_j}}{\mathrm{d}r}+\frac{1}{r}({\xi }_{r;m}^{l_j}-{\xi }_{H;m}^{l_j})\hfill \\ \hfill {\mathcal{B}}_{R;m}^{l_j\prime }& =\frac{\mathrm{d}{\xi }_{T;m}^{l_j\prime }}{\mathrm{d}r}-\frac{{\xi }_{T;m}^{l_j\prime }}{r},\hfill \end{array}$$(A.15)

while c.c denotes the complex conjugate. This allows us to compute

$$\begin{array}{c}\hfill R_{2;m}^j(r)=2\alpha R_{1;m}^j(r)\end{array}$$(A.16)

with α = 1/3. Finally, we have to compute

$$\begin{array}{c}\hfill R_{3;m}^j(r)=\beta (R_{1;m}^j(r)+2{\displaystyle {\int }_{4\pi }}\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Omega }}}{4\pi }{\displaystyle \sum _{i,j;i\ne j}}(B_{\mathit{ii}}^{\ast }{B}_{\mathit{jj}}+B_{\mathit{ii}}{B}_{\mathit{jj}}^{\ast })),\end{array}$$(A.17)

where we refer the reader to the formalism explained in Appendix A of Belkacem et al. (2009) for the definition of the B_ii coefficients (see their Eqs. (A.5) and (A.7)) while β = 1/5, $\stackrel{\sim }{\mathrm{\Omega }}$ is the solid angle with $\mathrm{d}\stackrel{\sim }{\mathrm{\Omega }}=sin\theta \mathrm{d}\theta \mathrm{d}\phi$ and ^* is the complex conjugate. We have:

$$\begin{array}{cc}& {\displaystyle \sum _{i,j;i\ne j}}(B_{\mathit{ii}}^{\ast }{B}_{\mathit{jj}}+B_{\mathit{ii}}{B}_{\mathit{jj}}^{\ast })\hfill \\ \hfill & =2{\left|\frac{{\xi }_r}{r}\right|}^2+2(\frac{\partial {\xi }_r^{\ast }}{\partial r}\frac{{\xi }_r}{r}+\frac{\partial {\xi }_r}{\partial r}\frac{{\xi }_r^{\ast }}{r})\hfill \\ \hfill & +[(\frac{\partial {\xi }_r^{\ast }}{\partial r}+\frac{{\xi }_r^{\ast }}{r}){\mathbf{\nabla }}_H·{\mathit{\xi }}_H+(\frac{\partial {\xi }_r}{\partial r}+\frac{{\xi }_r}{r}){\mathbf{\nabla }}_H·{\mathit{\xi }}_H^{\ast }]\hfill \\ \hfill & +\frac{1}{r}\frac{\partial {\xi }_{\theta }}{\partial \theta }{(\frac{{\xi }_{\theta }}{r}\frac{cos\theta }{sin\theta }+\frac{1}{rsin\theta }\frac{\partial {\xi }_{\varphi }}{\partial \varphi })}^{\ast }+\mathrm{c}.\mathrm{c}.,\hfill \end{array}$$

where c.c. means the complex conjugate.

Upon integration over the surface, we obtain

$$\begin{array}{cc}& {\displaystyle {\int }_{4\pi }\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Omega }}}{4\pi }\sum _{i,j;i\ne j}(B_{\mathit{ii}}^{\ast }{B}_{\mathit{jj}}+B_{\mathit{ii}}{B}_{\mathit{jj}}^{\ast })}\hfill \\ \hfill & ={\displaystyle \sum _j\{2{\left|\frac{{\xi }_{r;m}^{l_j}}{r}\right|}^2+2(\frac{\mathrm{d}{\left({\xi }_{r;m}^{l_j}\right)}^{\ast }}{\mathrm{d}r}\frac{{\xi }_{r;m}^{l_j}}{r}+\mathrm{c}.\mathrm{c}.)}\hfill \\ \hfill & -L_{l_j}^2[(\frac{\mathrm{d}{\left({\xi }_{r;m}^{l_j}\right)}^{\ast }}{\mathrm{d}r}+\frac{{\left({\xi }_{r;m}^{l_j}\right)}^{\ast }}{r})\frac{{\xi }_{H;m}^{l_j}}{r}+\mathrm{c}.\mathrm{c}.]\}\hfill \\ \hfill & +{\displaystyle {\int }_{4\pi }\frac{\mathrm{d}\stackrel{\sim }{\mathrm{\Omega }}}{4\pi }[\frac{1}{r}\frac{\partial {\xi }_{\theta }}{\partial \theta }{(\frac{{\xi }_{\theta }}{r}\frac{cos\theta }{sin\theta }+\frac{1}{rsin\theta }\frac{\partial {\xi }_{\phi }}{\partial \phi })}^{\ast }+\mathrm{c}.\mathrm{c}.].}\hfill \end{array}$$(A.18)

To evaluate the last terms, we use the following procedure:

$$\begin{array}{cc}\hfill \frac{1}{r}\frac{\partial {\xi }_{\theta }}{\partial \theta }& ={\displaystyle \sum _j\{\frac{{\xi }_{H;m}^{l_j}}{r}[-{\mathrm{\Lambda }}_{l_j}{Y}_{l_j}^m-\frac{1}{{sin}^2\theta }(cos\theta sin\theta \frac{\partial Y_{l_j}^m}{\partial \theta }-m^2{Y}_{l_j}^m)]}\hfill \\ \hfill & +m\left(i\frac{{\xi }_{T;m}^{l_j^{{}^{\prime }}}}{r}\right)\frac{1}{{sin}^2\theta }(sin\theta \frac{\partial Y_{l_j^{{}^{\prime }}}^m}{\partial \theta }-cos\theta Y_{l_j^{{}^{\prime }}}^m)\}\hfill \\ \hfill & ={\displaystyle \sum _{j,k}[\frac{{\xi }_{H;m}^{l_j}}{r}(-{\mathrm{\Lambda }}_{l_j{l}_k}+A_{l_j{l}_k})+m\left(\mathrm{i}\frac{{\xi }_{T;m}^{l_j^{{}^{\prime }}}}{r}\right)B_{l_j^{{}^{\prime }}{l}_k}]Y_{l_k}^m,}\hfill \end{array}$$(A.19)

$$\begin{array}{cc}\hfill \frac{1}{rsin\theta }\frac{\partial {\xi }_{\phi }}{\partial \phi }+\frac{cos\theta }{rsin\theta }{\xi }_{\theta }& ={\displaystyle \sum _j\{\frac{{\xi }_{H;m}^{l_j}}{r}\frac{1}{{sin}^2\theta }(cos\theta sin\theta \frac{\partial Y_{l_j}^m}{\partial \theta }-m^2{Y}_{l_j}^m)}\hfill \\ \hfill & +m\left(\mathrm{i}\frac{{\xi }_{T;m}^{l_j^{{}^{\prime }}}}{r}\right)\frac{1}{{sin}^2\theta }(cos\theta Y_{l_j^{{}^{\prime }}}^m-sin\theta \frac{\partial Y_{l_j^{{}^{\prime }}}^m}{\partial \theta })\}\hfill \\ \hfill & =-{\displaystyle \sum _{j,k}[\frac{{\xi }_{H;m}^{l_j}}{r}{A}_{l_j{l}_k}+m\left(\mathrm{i}\frac{{\xi }_{T;m}^{l_j^{{}^{\prime }}}}{r}\right)B_{l_j^{{}^{\prime }}{l}_k}]Y_{l_k}^m}\hfill \end{array}$$(A.20)

where we introduce the matrices (Λ_lk), (A_lk), and (B_lk) defined by

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{l_j{l}_k}={\delta }_{l_j{l}_k}{L}_{l_j}^2,\end{array}$$(A.21)

$$\begin{array}{c}\hfill A_{l_j{l}_k}=-{\displaystyle {\int }_{4\pi }\frac{d\stackrel{\sim }{\mathrm{\Omega }}}{4\pi }\frac{1}{{sin}^2\theta }(cos\theta sin\theta \frac{\partial Y_{l_j}^m}{\partial \theta }-m^2{Y}_{l_j}^m){(Y_{l_k}^m)}^{\ast },}\end{array}$$(A.22)

$$\begin{array}{c}\hfill B_{l_j\prime l_k}={\displaystyle {\int }_{4\pi }\frac{d\stackrel{\sim }{\mathrm{\Omega }}}{4\pi }\frac{1}{{sin}^2\theta }(sin\theta \frac{\partial Y_{l_j\prime }^m}{\partial \theta }-cos\theta Y_{l_j\prime }^m){(Y_{l_k}^m)}^{\ast }.}\end{array}$$(A.23)

This leads to:

$$\begin{array}{cc}& {\displaystyle {\int }_{4\pi }\frac{d\stackrel{\sim }{\mathrm{\Omega }}}{4\pi }[\frac{1}{r}\frac{\partial {\xi }_{\theta }}{\partial \theta }{(\frac{{\xi }_{\theta }}{r}\frac{cos\theta }{sin\theta }+\frac{1}{rsin\theta }\frac{\partial {\xi }_{\phi }}{\partial \phi })}^{\ast }+\mathrm{c}.\mathrm{c}.]}\hfill \\ \hfill & =-{\displaystyle \sum _{j,k}(\frac{{\xi }_{H;m}^{l_j}}{r}{\left(\frac{{\xi }_{H;m}^{l_k}}{r}\right)}^{\ast }+\mathrm{c}.\mathrm{c}.)(-A_{l_j{l}_k}{\mathrm{\Lambda }}_{l_k}+\sum _n{A}_{l_j{l}_n}{A}_{l_k{l}_n})}\hfill \\ \hfill & -{\displaystyle \sum _{j,k}[i\frac{{\xi }_{T;m}^{l_j\prime }}{r}{\left(\frac{{\xi }_{H;m}^{l_k}}{r}\right)}^{\ast }+\mathrm{c}.\mathrm{c}.](-B_{l_j\prime l_k}{\mathrm{\Lambda }}_{l_k}+2\sum _n{B}_{l_j\prime l_n}{A}_{l_k{l}_n})}\hfill \\ \hfill & -{\displaystyle \sum _{j,k}[i\frac{{\xi }_{T;m}^{l_j\prime }}{r}{\left(i\frac{{\xi }_{T;m}^{l_k\prime }}{r}\right)}^{\ast }+\mathrm{c}.\mathrm{c}.]\sum _n{B}_{l_j\prime l_n}{B}_{l_k\prime l_n}.}\hfill \end{array}$$(A.24)

We finally obtain:

$$\begin{array}{cc}\hfill R_m^j(r)& =\frac{16}{15}{\left|\frac{\mathrm{d}{\xi }_{r;m}^{l_j}}{\mathrm{d}r}\right|}^2+\frac{44}{15}{\left|\frac{{\xi }_{r;m}^{l_j}}{r}\right|}^2+\frac{4}{5}(\frac{{\xi }_{r;m}^{l_j}}{r}\frac{\mathrm{d}{\left({\xi }_{r;m}^{l_j}\right)}^{\ast }}{\mathrm{d}r}+\mathrm{c}.\mathrm{c}.)\hfill \\ \hfill & +\frac{8}{5}{L}_{l_j}^2{\left|{\mathcal{A}}_{R;m}^{l_j}\right|}^2-\frac{22}{15}{L}_{l_j}^2[\frac{{\left({\xi }_{r;m}^{l_j}\right)}^{\ast }{\xi }_{H,m}^{l_j}}{r^2}+\mathrm{c}.\mathrm{c}.]\hfill \\ \hfill & -\frac{2}{5}{L}_{l_j}^2[\left(\frac{\mathrm{d}{\left({\xi }_{r;m}^{l_j}\right)}^{\ast }}{\mathrm{d}r}\right)\frac{{\xi }_{H;m}^{l_j}}{r}+\mathrm{c}.\mathrm{c}.]\hfill \\ \hfill & +\frac{16}{15}{L}_{l_j}^2(L_{l_j}^2-1){\left|\frac{{\xi }_{H;m}^{l_j}}{r}\right|}^2\hfill \\ \hfill & -\frac{2}{5}{\displaystyle \sum _{j,k}(\frac{{\xi }_{H;m}^{l_j}}{r}{\left(\frac{{\xi }_{H;m}^{l_k}}{r}\right)}^{\ast }+\mathrm{c}.\mathrm{c}.)(-A_{l_j{l}_k}{\mathrm{\Lambda }}_{l_k}+\sum _n{A}_{l_j{l}_n}{A}_{l_k{l}_n})}\hfill \\ \hfill & +\frac{8}{5}{L}_{l_j\prime }^2(L_{l_j\prime }^2-2){\left|\frac{{\xi }_{T;m}^{l_{j\prime }}}{r}\right|}^2\hfill \\ \hfill & +\frac{8}{5}{L}_{l_j\prime }^2{\left|{\mathcal{B}}_{R;m}^{l_{j\prime }}\right|}^2\hfill \\ \hfill & -\frac{2}{5}{\displaystyle \sum _{j,k}[i\frac{{\xi }_{T;m}^{l_j\prime }}{r}{\left(i\frac{{\xi }_{T;m}^{l_k\prime }}{r}\right)}^{\ast }+\mathrm{c}.\mathrm{c}.]\sum _n{B}_{l_j\prime l_n}{B}_{l_k\prime l_n}}\hfill \\ \hfill & -\frac{2}{5}{\displaystyle \sum _{j,k}[i\frac{{\xi }_{T;m}^{l_j\prime }}{r}{\left(\frac{{\xi }_{H;m}^{l_k}}{r}\right)}^{\ast }+\mathrm{c}.\mathrm{c}.](-B_{l_j\prime l_k}{\mathrm{\Lambda }}_{l_k}+2\sum _n{B}_{l_j\prime l_n}{A}_{l_k{l}_n}).}\hfill \\ \hfill \end{array}$$(A.25)

Furthermore,

$$\begin{array}{c}\hfill S_R({\omega }_0)={\displaystyle {\int }_0^{\infty }\frac{\text{d}k}{k^2}{E}^2(k){\int }_{-\infty }^{\infty }\text{d}\omega {\chi }_k(\omega +{\omega }_0){\chi }_k(\omega ),}\end{array}$$(A.26)

where (k, ω) are the wavenumber and frequency of the turbulent eddies and E(k, ω) is the turbulent kinetic energy spectrum, which is expressed as the product E(k) χ_k(ω) for isotropic turbulence. Indeed, as it will be discussed in detail in Appendix A.2.3., here we follow Samadi et al. (2010) and Samadi & Goupil (2001) by making the quasi-normal approximation and the assumptions of stationary, incompressible, homogeneous, and isotropic turbulence. These hypotheses allow us to explicitly relate the Reynolds stress integrals to the energy spectra. We refer the reader to Eqs. (25)–(28) and (34) in Samadi & Goupil (2001) for more technical details. The same assumptions are made to express the entropy fluctuations contribution. E(k) is the spatial kinetic energy spectrum, while χ_k(ω) is the eddy-time correlation function. In Appendix A.3, a detailed discussion of how to describe those functions will be given.

A.2.2. Entropy fluctuations contribution

As shown by Samadi & Goupil (2001) and discussed briefly in Samadi et al. (2010), the Reynolds stress contribution is not the unique source of excitation. One also has to evaluate the contribution to the excitation by the entropy fluctuations advection. Following Belkacem et al. (2008a), this entropy source term depends on the mode compressibility. Consequently, since the divergence of the toroidal component, which is the rotational part of the spherical harmonic, vanishes, one can consider the poloidal contribution only. Therefore, the final expression for the contribution of entropy fluctuations remains the same for rapidly rotating stars as for non-rotating stars (Belkacem et al. 2008a), that is

$$\begin{array}{c}\hfill C_S^2=\frac{4{\pi }^3\mathcal{H}}{{\omega }_0^2}{\displaystyle {\int }_{{\mathcal{V}}_{\mathrm{CZ}}}{\text{d}}^3{x}_0{\alpha }_s^2\sum _j({\mathcal{A}}_{S;m}^{l_j}+{\mathcal{B}}_{S;m}^{l_j}){\mathcal{S}}_S({\omega }_0),}\end{array}$$(A.27)

where ℋ is the anisotropy factor introduced in Samadi & Goupil (2001) which, in the current assumption (isotropic turbulence), is equal to 4/3. In addition,

$$\begin{array}{cc}\hfill {\mathcal{A}}_{S;m}^{l_j}& \equiv \frac{1}{r^2}{\mathcal{D}}_m^{l_j}{|\frac{\mathrm{d}(ln\mid {\alpha }_s\mid )}{\mathrm{d}lnr}-\frac{\mathrm{d}{\mathcal{D}}_m^{l_j}}{\mathrm{d}lnr}|}^2,\hfill \end{array}$$(A.28)

$$\begin{array}{cc}\hfill {\mathcal{B}}_{S;m}^{l_j}& \equiv \frac{1}{r^2}{L_j}^2{\left|{\mathcal{D}}_m^{l_j}\right|}^2,\hfill \end{array}$$(A.29)

where

$$\begin{array}{c}\hfill {\mathcal{D}}_m^{l_j}(r,l_j)\equiv {\mathcal{D}}_r-\frac{{L_j}^2}{r}{\xi }_{H;m}^{l_j},{\mathcal{D}}_r\equiv \frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2{\xi }_{r;m}^{l_j}\right)\end{array}$$(A.30)

and

$$\begin{array}{c}\hfill {\mathcal{S}}_S({\omega }_0)\equiv {\displaystyle {\int }_0^{\infty }\frac{\text{d}k}{k^4}E(k)E_s(k){\int }_{-\infty }^{\infty }\text{d}\omega {\chi }_k({\omega }_0+\omega ){\chi }_k(\omega ),}\end{array}$$(A.31)

where E_s(k) is the spatial entropy spectrum while χ_k(ω) is the same eddy-time correlation function as in the Reynolds contribution.

A.2.3. Rotational contributions

To get clues about the contributions of the rotational terms, we express them in a suitable form following Belkacem et al. (2009):

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}(2{\rho }_1\mathbf{\Omega }\times {\mathit{u}}_t)-\frac{\partial }{\partial t}\left({\rho }_1\mathrm{\Omega }\frac{\partial {\mathit{u}}_t}{\partial \phi }\right)=-[2\mathbf{\Omega }\times {\mathit{u}}_t+\mathrm{\Omega }\frac{\partial {\mathit{u}}_t}{\partial \phi }]\mathbf{\nabla }·\left({\mathit{\rho }}_0{\mathit{u}}_t\right).\end{array}$$(A.32)

Those terms vanish under the assumption that the turbulent field is anelastic. This approximation is verified because we are dealing with excitation of modes in convective zones. Indeed, the energy equation at the first order, together with the mass conservation equation (i.e. ∇ ⋅ (ρ₀u_t) = 0), is equivalent to a quasi-isentropic equilibrium state.

In addition, by using dimensional analysis (Samadi & Goupil 2001), it appears that all those terms scale as the Mach number to the third (ℳ³) because they present a dependency to the perturbed mass flux in ρ₁ u_t. Compared to the Reynolds contribution, which scales as ℳ², all rotational contributions thus become negligible in the subsonic regime (ℳ <  1).

A.3. Turbulence modelling

To model turbulence, we first have to specify the spatial kinetic energy spectrum (E(k)) and the spatial energy spectrum for entropy (E_s(k)).

While it is tempting to adopt relatively simple models for the nature of rotating convection, the stochastic excitation mechanisms that can drive gravito-inertial waves are subtle as they depend upon the precise nature of the turbulent convection in the driving region. The most widely used approach is to assume that the convective motions driving the waves follow a Heisenberg-Kolmogorov spectrum, which is also the tack taken here. In contrast to the well-tested Heisenberg-Kolmogorov theory, the spectral form of the Reynolds stress correlations that arise in buoyancy-driven turbulent rotating convection are highly model dependent and as yet do not have a closed form solution. This difficulty arises because of the broken symmetries of inhomogeneous and anisotropic turbulence that exist in these systems. Given that the Coriolis force depends upon the velocity, the random velocity distributions as defined by Kolmogorov depend not only upon the coordinates but also upon the initial velocity distributions. Therefore, the problem of defining simple and integrable spectral distributions of the Reynolds stresses that are linked to the kinetic energy spectrum is rendered completely inhomogeneous and therefore time dependent, breaking both of Kolmogorov’s hypotheses that permit the construction of structure functions.

Nevertheless, there are parameter regimes in which approximate turbulence models of rotating thermal convection can be constructed. Depending upon that regime, the spectral scaling of the Reynolds stresses has been found to vary as $\widehat{v^2}\propto k^{\eta (k)}$, where η(k) can range from being a constant (e.g. Mininni & Pouquet 2010) to a function of wavenumber to account for various integral scales such as the Bolgiano scale (e.g. Pharasi et al. 2014). Thus, for simplicity and as a first approximation, the turbulence model presented in Samadi & Goupil (2001), Samadi et al. (2010) and Belkacem et al. (2008b) is adopted here, that is the Kolmogorov one, with the testing of more sophisticated rotating convection models being reserved for future work.

Therefore, we assume that the spatial structure of the convection is taken to be that of isotropic, homogeneous, and incompressible turbulence, which has a power law in the spatial transform of the kinetic energy for scales far from the integral driving scale and away from the diffusive microscale (Eq. (A.36)). These assumptions that the motions are homogeneous and isotropic, have meaning only in a statistical sense. That is to say that the turbulence model is built using random variables whose correlations are constructed such that they obey the expectation value of the energy equation applied to them (we refer the reader to Eq. (14) in Kolmogorov 1941). To paraphrase Kolmogorov and Taylor, this turbulence model can be considered in the following fashion: for very large Reynolds numbers, the first-order perturbations on the average flow field consist of disorderly displacements of separate fluid volumes, one with respect to another, of a length scale ℓ₀ (where ℓ₀ is the Prandtl’s mixing path length); these first-order perturbations are themselves of such a high amplitude that they are in turn unsteady and on these structures are superposed fluctuations of the second order with a mixing length ℓ₁ <  ℓ₀ and relative velocities v₁ <  v₀; such a process of successive refinement of turbulent perturbations may be carried out until for some sufficiently high order n the Reynolds number R_n = (l_nv_n)/ν, where ν is the kinematic viscosity, becomes so small that the effect of viscosity on the perturbations of the order n finally prevents the formation of further fine structure at the order n + 1. To put that more succinctly, the spatial structure of the flow is a large scale flow with superposed self-similar eddies that obey the scale-by-scale averaged dissipation equation (see again Eq. (14) Kolmogorov 1941). Yet such motions are only captured in an ergodic statistical sense in the turbulence model of Kolmogorov used here.

Nevertheless, the precise and time-evolving spectrum of the waves will be coupled to the actual structure of the turbulent flows that excite them, which are modified by rotation and not fully captured in the Kolmogorov’s picture. However, our measurements are integrated in time and space (over the disk of the star), and therefore represent an averaging operation, which is somewhat akin to the averaging procedure used to construct the Kolmogorov structured turbulence and certainly averages over some of the timescales of the convective motions. So, as mentioned above, this is a first approximation for the modelling of the turbulence exciting the waves and additional theoretical efforts are ongoing to better address the interactions of waves with rotating turbulent structures from small to large scales (e.g. Mathis et al. 2014; Augustson & Mathis 2019; Augustson et al. 2020), but they are beyond the scope of this work.

Following Samadi & Goupil (2001), we introduce a velocity u₀ and the entropy scalar variance $\stackrel{\sim }{s}$ defined as

$$\begin{array}{c}\hfill \frac{3}{2}{u}_0^2=\frac{1}{2}\mathrm{\Phi }{w}^2={\displaystyle {\int }_0^{\infty }\mathrm{d}kE\left(k\right)}\end{array}$$(A.33)

and

$$\begin{array}{c}\hfill \frac{1}{2}{\stackrel{\sim }{s}}^2=\frac{1}{2}{\left(\frac{2\mathrm{\Phi }{c}_p}{g\mathrm{\Lambda }\delta }\right)}^2{w}^4={\displaystyle {\int }_0^{\infty }\mathrm{d}k{E}_s\left(k\right),}\end{array}$$(A.34)

where w is the estimate for the vertical convective velocity deduced from the mixing length theory (MLT, Böhm-Vitense 1958), while Φ is a factor introduced by Gough (1977) to take the anisotropy of the convective turbulence into account. Here, we follow Samadi & Goupil (2001) and we fix Φ = 2, which is the value that corresponds to the Bohm & Vintense’s MLT. c_p is the heat capacity at constant pressure, g is the gravity and δ = (∂lnρ₀/∂lnT₀)_P. Finally, Λ = α_cH_p is the usual mixing length, where H_p is the pressure scale height while the mixing-length parameter α_c = 1.5 in the computed stellar models which are used here to study the amplitude of stochastically-excited waves.

In addition, we define k₀ as the wavenumber of the largest eddy in the inertial range. Following Samadi & Goupil (2001), we relate it to the mixing length by

$$\begin{array}{c}\hfill k_0=\frac{2\pi }{\widehat{\mathrm{\Lambda }}}=\frac{2\pi }{\widehat{\beta }\mathrm{\Lambda }},\end{array}$$(A.35)

where $\widehat{\beta }$ is a free parameter. Samadi & Goupil (2001) demonstrated that for Φ = 2, $\widehat{\beta }\approx 1.9$ (we refer the reader to their Eq. (76)). However, Samadi et al. (2010) suggested choosing the size of the convective core for early-type stars instead of $\widehat{\beta }\mathrm{\Lambda }$ since this is the natural size of the possible largest eddy (see e.g. the 3D spherical numerical simulations by Browning et al. 2004; Augustson et al. 2016; Edelmann et al. 2019). We hereafter adopt this choice for the computation of the amplitude of the waves.

The Kolmogorov spectrum is then expressed following Samadi & Goupil (2001);

$$\begin{array}{c}\hfill E\left(K\right)=\frac{u_0^2}{k_0}{K}^{-5/3},\end{array}$$(A.36)

where K = k/k₀, while

$$\begin{array}{c}\hfill E_s\left(k\right)=\frac{1}{3}{\stackrel{\sim }{s}}^2\frac{E\left(k\right)}{u_0^2}·\end{array}$$(A.37)

Finally, we have to specify the time correlation function (χ_k). As shown by Samadi et al. (2003) for solar p modes and by Belkacem et al. (2009) for solar g modes, the choice of this function has a strong impact on the mode amplitudes. In the case of massive stars, this has been confirmed by Belkacem et al. (2010) and Samadi et al. (2010), who showed that the choice of a Lorentzian function is the most favorable. Therefore, we follow Belkacem et al. (2010) and we define

$$\begin{array}{c}\hfill {\chi }_k\left(\omega \right)=\frac{1}{\pi {\omega }_k}\frac{1}{1+{\left(\frac{\omega }{{\omega }_k}\right)}^2}\end{array}$$(A.38)

with the normalisation condition

$$\begin{array}{c}\hfill {\displaystyle {\int }_{-\infty }^{\infty }{\chi }_k\left(\omega \right)\mathrm{d}\omega =1.}\end{array}$$(A.39)

We define the linewidth

$$\begin{array}{c}\hfill {\omega }_k=\frac{k{u}_k}{\lambda },\end{array}$$(A.40)

where λ is a free parameter that accounts for the uncertainties in the definition of ω_k (see Samadi & Goupil 2001; Belkacem et al. 2010), while u_k is the velocity of the eddy with the wavenumber k, which was related to the kinetic energy spectrum by Stein (1967):

$$\begin{array}{c}\hfill u_k^2={\displaystyle {\int }_k^{2k}\mathrm{d}kE\left(k\right).}\end{array}$$(A.41)

Since the lifetime of the largest eddies in the inertial range

$$\begin{array}{c}\hfill {\tau }_{k_0}=\frac{\lambda }{\left(k_0{u}_{k_0}\right)}\end{array}$$(A.42)

cannot be longer than the characteristic time at which the convective energy dissipates into the turbulent cascade, we have

$$\begin{array}{c}\hfill {\tau }_{k_0}\le \frac{\mathrm{\Lambda }}{w}.\end{array}$$(A.43)

Using the definitions of k₀ (Eq. (A.35)), u_k (Eq. (A.41)), and w (Eq. (A.33)), we get $\lambda \le 2.7\sqrt{\mathrm{\Phi }}\approx 3.82$; we here choose λ = 1.