Issue 
A&A
Volume 644, December 2020



Article Number  A9  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201935858  
Published online  24 November 2020 
Transport of angular momentum by stochastically excited waves as an explanation for the outburst of the rapidly rotating Be star HD49330^{⋆}
^{1}
LESIA, Paris Observatory, PSL University, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France
email: coralie.neiner@obspm.fr
^{2}
Astronomical Institute, Graduate School of Science, Tohoku University, Sendai 9808578, Japan
^{3}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 GifsurYvette Cedex, France
^{4}
Physics Department, Mount Allison University, Sackville, NB E4L 1E6, Canada
Received:
9
May
2019
Accepted:
7
July
2020
Context. HD 49330 is an early Be star that underwent an outburst during its fivemonth observation with the CoRoT satellite. An analysis of its light curve revealed several independent p and g pulsation modes, in addition to showing that the amplitude of the modes is directly correlated with the outburst.
Aims. We modelled the results obtained with CoRoT to understand the link between pulsational parameters and the outburst of this Be star.
Methods. We modelled the flattening of the structure of the star due to rapid rotation in two ways: ChandrasekharMilne’s expansion and 2D structure computed with ROTORC. We then modelled κdriven pulsations. We also adapted the formalism of the excitation and amplitude of stochastically excited gravitoinertial modes to rapidly rotating stars, and we modelled those pulsations as well.
Results. We find that while pulsation p modes are indeed excited by the κ mechanism, the observed g modes are, rather, a result of stochastic excitation. In contrast, g and r waves are stochastically excited in the convective core and transport angular momentum to the surface, increasing its rotation rate. This destabilises the external layers of the star, which then emits transient stochastically excited g waves. These transient waves produce most of the lowfrequency signal detected in the CoRoT data and ignite the outburst. During this unstable phase, p modes disappear at the surface because their cavity is broken. Following the outburst and ejection of the surface layer, relaxation occurs, making the transient g waves disappear and p modes reappear.
Conclusions. This work includes the first coherent model of stochastically excited gravitoinertial pulsation modes in a rapidly rotating Be star. It provides an explanation for the correlation between the variation in the amplitude of frequencies detected in the CoRoT data and the occurrence of an outburst. This scenario could apply to other pulsating Be stars, providing an explanation to the longstanding questions surrounding Be outbursts and disks.
Key words: stars: emissionline, Be / stars: oscillations / stars: rotation / stars: individual: HD49330
© C. Neiner et al. 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Be stars are lateO to earlyA stars surrounded by a decretion disk fed by discrete massloss 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 gravitoinertial 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 groundbased 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 at the equator, T_{eff} = 28 000 ± 1500 K, log g = 3.9 ± 0.2, v sin i = 280 ± 40 km s^{−1}, 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 is commonly used in observational studies, seismic models rather consider , which is the rotation frequency normalised by , where R is the mean stellar radius. For HD 49330, corresponds to .
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 stochasticallydriven 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:

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.

A deformation following ChandrasekharMilne’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_{2}(cos θ)] where a is the mean radius, P_{2}(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 ChandrasekharMilne’s expansion is only valid when the difference between the rotating and nonrotating structures can be approximated by Ω^{2} 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 ChandrasekharMilne expansion and the 2D structure model are similar. For , however, the deformations from the 2D models are considerably larger than those from ChandrasekharMilne’s expansion, in particular in the outer layers.
Fig. 1. Coefficient β of the Legendre function P_{2}(cos θ) as a function of the fractional radius for a 10 M_{⊙} star rotating at various rates. Solid lines represent the deformation from the ChandrasekharMilne’s expansion and dotted lines the deformation from 2D structure models. is the rotation frequency normalised by where R is the mean radius. increases from top to bottom in the figure. corresponds to at the equator as observed for HD 49330. 
2.2. Tohoku models of κdriven pulsations
To model the pulsations, we use the Tohoku nonadiabatic pulsation code. This code includes the effects of the centrifugal and of the Coriolis accelerations. With this code, we perform a nonradial 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 nonadiabatic 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 nonadiabatic models using a stellar structure deformation based either on ChandrasekharMilne’s expansion or on a 2D stellar structure, as presented in Sect. 2.1.
Using ChandrasekharMilnes’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 signaltonoise 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^{−4}, meaning that it would take 10^{4} 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 ChandrasekharMilnes’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 ChandrasekharMilne’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 highfrequency 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^{−1} can be identified with a retrograde p mode with m = 1 or 2, while the frequency groups around 1.5 and 3 c d^{−1} 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).
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^{−1} is not present anymore in the model. In addition, the model predicts several strong modes between 6 and 10 c d^{−1} 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 (). 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) nonadiabatic 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 spectroscopicallyderived 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^{4} 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 stochasticallyexcited 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 spectroscopicallydetermined 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^{−1} (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 gravitoinertial modes with the nonadiabatic 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 gravitoinertial 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). Gravitoinertial 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 gravitoinertial 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 nonradial 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 nonradial modes in a nonrotating 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 , 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 , 0.4, 0.6, and 0.8. Since the star is considered spherical here, , 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 nonadiabatic 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: (left), (middle left), (middle right), and (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 yaxis 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 nonadiabatic 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 is similar to the amplitudes detected with CoRoT, which confirms that we likely see stochastically excited gravitoinertial 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 gravitoinertial 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 . 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 gravitoinertial waves and the convective forcing, in particular in the subinertial regime in which gravitoinertial 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 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 , 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 subinertial 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 gravitoinertial waves, which are stochastically excited inside the star, and Be outbursts.
First, we recall that the idea that nonradial 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, gravitoinertial waves (and r modes) were excited by the κmechanism. In the case of HD 49330, we propose that gravitoinertial 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 gravitoinertial 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:
where the subscript α identifies the excited mode, and . Here, is the specific angular momentum around the rotation axis, where Ω is the considered uniform rotation, and is the induced differential rotation. Each scalar field (X) has been expanded as , where is the mean value on an isobar and X′ is the modal fluctuation. The vertical flux of angular momentum is given by
where we recall that, for modes calculated taking only the given uniform rotation into account,v = i(ω_{0}+mΩ)ξ.
We define the associated timescale for the transport of angular momentum τ as:
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 (left), 0.4 (middle left), 0.6 (middle right) and 0.8 (right), where τ is normalised to . The value of 1/τ is proportional to the angular momentum transport. Note that the right panel (fastest rotator) has a different yaxis range to accommodate the larger 1/τ values and shows an extra yaxis 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^{5} (1/τ = 1.2 × 10^{−6}) for the sum of prograde and retrograde modes in a rapidly rotating star. Since Ω_{c} is about 1.5 × 10^{−4} s^{−1} 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
can also be computed. In the case of a shellular rotation that depends only on the radial coordinate, it becomes
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 then τ′→0, meaning if the star rotates at critical velocity it looses mass all the time, while if then τ′→∞, meaning if a star does not rotate no outburst occurs. means in the nonspherical case. We then predict that an outburst should occur every ∼11 years in HD 49330. Photometric observations by ASAS3 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^{−1} 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:

The treatment of rotational deformation is very important for highfrequency p modes, because the excitation of these modes seems to occur only when the rotation is nearly critical. Consequently, ChandrasekharMilne’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.

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 gravitoinertial 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.

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 . 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).

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.

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.

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 longstanding 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 gravitoinertial 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 gmode 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 gravitoinertial 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 gravitoinertial waves.
Finally, we considered κdriven p and g modes as well as stochastically excited gravitoinertial (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/DApAIM. This research has made use of the SIMBAD database operated at CDS, Strasbourg (France), and of NASA’s Astrophysics Data System (ADS).
References
 Aerts, C., & Rogers, T. M. 2015, ApJ, 806, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Alvan, L., Brun, A. S., & Mathis, S. 2014, A&A, 565, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ando, H. 1986, A&A, 163, 97 [NASA ADS] [Google Scholar]
 Augustson, K. C., & Mathis, S. 2019, ApJ, 874, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Augustson, K. C., Brun, A. S., & Toomre, J. 2016, ApJ, 829, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Augustson, K. C., Mathis, S., & Astoul, A. 2020, ApJ, 903, 90 [CrossRef] [Google Scholar]
 Auvergne, M., Bodin, P., Boisnard, L., et al. 2009, A&A, 506, 411 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ballot, J., Lignières, F., Reese, D. R., & Rieutord, M. 2010, A&A, 518, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., Goupil, M.J., & Dupret, M.A. 2008a, A&A, 478, 163 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Samadi, R., & Goupil, M. J. 2008b, J. Phys. Conf. Ser., 118, 012028 [CrossRef] [Google Scholar]
 Belkacem, K., Mathis, S., Goupil, M. J., & Samadi, R. 2009, A&A, 508, 345 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Belkacem, K., Dupret, M. A., & Noels, A. 2010, A&A, 510, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BöhmVitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
 Browning, M. K., Brun, A. S., & Toomre, J. 2004, ApJ, 601, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Cantiello, M., Langer, N., Brott, I., et al. 2009, A&A, 499, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deupree, R. G. 1990, ApJ, 357, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Deupree, R. G. 1995, ApJ, 439, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Dintrans, B., & Rieutord, M. 2000, A&A, 354, 86 [NASA ADS] [Google Scholar]
 Edelmann, P. V. F., Ratnasingam, R. P., Pedersen, M. G., et al. 2019, ApJ, 876, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Floquet, M., Hubert, A.M., Huat, A.L., et al. 2009, A&A, 506, 103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frémat, Y., Zorec, J., Hubert, A.M., & Floquet, M. 2005, A&A, 440, 305 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gough, D. O. 1977, ApJ, 214, 196 [NASA ADS] [CrossRef] [Google Scholar]
 Grimshaw, R. 1984, Ann. Rev. Fluid Mech., 16, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Huat, A.L., Hubert, A.M., Baudin, F., et al. 2009, A&A, 506, 95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kolmogorov, A. 1941, Akademiia Nauk SSSR Doklady, 30, 301 [Google Scholar]
 Lee, U. 2006, MNRAS, 365, 677 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U. 2012, MNRAS, 420, 2387 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., & Baraffe, I. 1995, A&A, 301, 419 [NASA ADS] [Google Scholar]
 Lee, U., & Saio, H. 1993, MNRAS, 261, 415 [Google Scholar]
 Lee, U., & Saio, H. 1997, ApJ, 491, 839 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, U., Neiner, C., & Mathis, S. 2014, MNRAS, 443, 1515 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, S. 2009, A&A, 506, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & de Brye, N. 2012, A&A, 540, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., & Prat, V. 2019, A&A, 631, A26 [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, S., Talon, S., Pantillon, F., & Zahn, J. 2008, Sol. Phys., 251, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mininni, P. D., & Pouquet, A. 2010, Physics of Fluids, 22, 035105 [NASA ADS] [CrossRef] [Google Scholar]
 Neiner, C., GutiérrezSoto, J., Baudin, F., et al. 2009, A&A, 506, 143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neiner, C., Floquet, M., Samadi, R., et al. 2012a, A&A, 546, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neiner, C., Mathis, S., Saio, H., et al. 2012b, A&A, 539, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pantillon, F. P., Talon, S., & Charbonnel, C. 2007, A&A, 474, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pharasi, H. K., Kumar, K., & Bhattacharjee, J. K. 2014, Phys. Rev. E, 89, 023009 [CrossRef] [Google Scholar]
 Pinçon, C., Belkacem, K., & Goupil, M. J. 2016, A&A, 588, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M. 2009, The Rotation of Sun and Stars, eds. J.P. Rozelot, & C. Neiner, Lect. Notes Phys., 765, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Rivinius, T., Baade, D., & Štefl, S. 2003, A&A, 411, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&ARv., 21, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M. 2015, ApJ, 815, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Saio, H. 2018, ArXiv eprints [arXiv:1812.01253] [Google Scholar]
 Saio, H., Kurtz, D. W., Murphy, S. J., Antoci, V. L., & Lee, U. 2018, MNRAS, 474, 2774 [Google Scholar]
 Samadi, R., & Goupil, M.J. 2001, A&A, 370, 136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Nordlund, Å., Stein, R. F., Goupil, M. J., & Roxburgh, I. 2003, A&A, 404, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Samadi, R., Belkacem, K., Goupil, M. J., et al. 2010, Ap&SS, 328, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Schatzman, E. 1993, A&A, 279, 431 [NASA ADS] [Google Scholar]
 Schenk, A. K., Arras, P., Flanagan, É. É., Teukolsky, S. A., & Wasserman, I. 2002, Phys. Rev. D, 65, 024001 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F. 1967, Sol. Phys., 2, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, A. A. 1965, J. Fluid Mech., 22, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Vadas, S. L., & Fritts, D. C. 2001, J. Atmos. Sci., 58, 2249 [CrossRef] [Google Scholar]
 Walker, G. A. H., Kuschnig, R., Matthews, J. M., et al. 2005, ApJ, 635, L77 [NASA ADS] [CrossRef] [Google Scholar]
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
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
where c_{s} is the sound velocity and α_{s} = (∂p/∂s)_{ρ}. Next, the operator
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
with
where ℒ_{t} contains the linear terms source^{1}. 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 ℳ^{3} 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 ℳ^{2}. Note that in the case of “free” turbulence where the oscillations do not have any effect on the turbulent flow, we have ρ_{1} ≈ ρ_{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
where 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 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
where ω_{0} 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
where the operator ⟨…⟩ denotes a statistical average performed on an infinite number of independent realisations, η is the damping rate computed by the nonadiabatic 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):
where
with
where
and .
Finally, is the Reynolds stress contribution, is the one related to entropy fluctuations, contains the contributions of the Coriolis acceleration and the Doppler term, and C_{c} represents the crosssource 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:
where 𝒱_{CZ} is the volume of the convective regions, with
(we refer the reader to Eq. (A.8) in Belkacem et al. 2009). We have:
where
while c.c denotes the complex conjugate. This allows us to compute
with α = 1/3. Finally, we have to compute
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, is the solid angle with and ^{*} is the complex conjugate. We have:
where c.c. means the complex conjugate.
Upon integration over the surface, we obtain
To evaluate the last terms, we use the following procedure:
where we introduce the matrices (Λ_{lk}), (A_{lk}), and (B_{lk}) defined by
This leads to:
We finally obtain:
Furthermore,
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 quasinormal 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 eddytime 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 nonrotating stars (Belkacem et al. 2008a), that is
where ℋ is the anisotropy factor introduced in Samadi & Goupil (2001) which, in the current assumption (isotropic turbulence), is equal to 4/3. In addition,
where
and
where E_{s}(k) is the spatial entropy spectrum while χ_{k}(ω) is the same eddytime 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):
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. ∇ ⋅ (ρ_{0}u_{t}) = 0), is equivalent to a quasiisentropic 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 (ℳ^{3}) because they present a dependency to the perturbed mass flux in ρ_{1} u_{t}. Compared to the Reynolds contribution, which scales as ℳ^{2}, 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 gravitoinertial 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 HeisenbergKolmogorov spectrum, which is also the tack taken here. In contrast to the welltested HeisenbergKolmogorov theory, the spectral form of the Reynolds stress correlations that arise in buoyancydriven 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 , 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 firstorder perturbations on the average flow field consist of disorderly displacements of separate fluid volumes, one with respect to another, of a length scale ℓ_{0} (where ℓ_{0} is the Prandtl’s mixing path length); these firstorder 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 ℓ_{1} < ℓ_{0} and relative velocities v_{1} < v_{0}; 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_{n}v_{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 selfsimilar eddies that obey the scalebyscale 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 timeevolving 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_{0} and the entropy scalar variance defined as
and
where w is the estimate for the vertical convective velocity deduced from the mixing length theory (MLT, BöhmVitense 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ρ_{0}/∂lnT_{0})_{P}. Finally, Λ = α_{c}H_{p} is the usual mixing length, where H_{p} is the pressure scale height while the mixinglength parameter α_{c} = 1.5 in the computed stellar models which are used here to study the amplitude of stochasticallyexcited waves.
In addition, we define k_{0} as the wavenumber of the largest eddy in the inertial range. Following Samadi & Goupil (2001), we relate it to the mixing length by
where is a free parameter. Samadi & Goupil (2001) demonstrated that for Φ = 2, (we refer the reader to their Eq. (76)). However, Samadi et al. (2010) suggested choosing the size of the convective core for earlytype stars instead of 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);
where K = k/k_{0}, while
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
with the normalisation condition
We define the linewidth
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):
Since the lifetime of the largest eddies in the inertial range
cannot be longer than the characteristic time at which the convective energy dissipates into the turbulent cascade, we have
Using the definitions of k_{0} (Eq. (A.35)), u_{k} (Eq. (A.41)), and w (Eq. (A.33)), we get ; we here choose λ = 1.
All Figures
Fig. 1. Coefficient β of the Legendre function P_{2}(cos θ) as a function of the fractional radius for a 10 M_{⊙} star rotating at various rates. Solid lines represent the deformation from the ChandrasekharMilne’s expansion and dotted lines the deformation from 2D structure models. is the rotation frequency normalised by where R is the mean radius. increases from top to bottom in the figure. corresponds to at the equator as observed for HD 49330. 

In the text 
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 ChandrasekharMilne’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. 

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

In the text 
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). 

In the text 
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). 

In the text 
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. 

In the text 
Fig. 7. Amplitude of the stochastically excited g and r modes for various rotation rates: (left), (middle left), (middle right), and (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 yaxis range to better see the low amplitudes. 

In the text 
Fig. 8. Top: same as Fig. 7 for . Bottom: statistically significant observed frequencies of HD 49330 during and outside of the outburst phase as detected by Huat et al. (2009). 

In the text 
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 (left), 0.4 (middle left), 0.6 (middle right) and 0.8 (right), where τ is normalised to . The value of 1/τ is proportional to the angular momentum transport. Note that the right panel (fastest rotator) has a different yaxis range to accommodate the larger 1/τ values and shows an extra yaxis on the right side in years^{−1}. 

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