© ESO, 2012

1. Introduction

Since the beginning of the theoretical modelling of stellar structure and evolution, many hydro- and magnetohydro-dynamical processes have been studied. In particular, the impact of the internal rotation and the associated shear has been modelled (e.g. Zahn 1992; Mathis & Zahn 2004) and evaluated (e.g. Meynet & Maeder 2000, for massive stars and Talon & Charbonnel 2005, for low-mass stars). For massive stars, especially for fast rotators close to their break-up velocity such as Be stars, the modification of stellar structure and evolution caused by the rotation has to be taken into account. Indeed, rotation impacts many aspects of stellar evolution such as the stellar life-time, internal structure (for both the thermic and chemical stratification), chemical yields, etc. However, such an hydrodynamical modelling always implies prescription, in particular for turbulence both in convective and radiative regions, and their limits. Therefore, it is absolutely necessary to obtain constraints from observations of these physical processes and their consequences, i.e. we have to probe stellar interiors to acquire their direct and indirect signatures.

In this context, high-precision asteroseismology from space is the best way to study the internal structure in the Hertzsprung-Russel (HR) diagram using oscillation mode frequencies, which give the properties of cavities where waves propagate. This is the framework of the present paper, in which we study the two fast rotating late Be stars HD 181231 and HD 175869. They have been observed with the CoRoT photometric satellite (for a description of CoRoT, see Baglin et al. 2006, 2009) during about five consecutive months and 27 consecutive days, respectively. An analysis of their CoRoT light curves showed that several independent periods of variations are present in these stars, which we interpreted in terms of pulsation g-modes (Neiner et al. 2009; Gutiérrez-Soto et al. 2009). Late Be stars indeed exhibit pulsations of g-modes similar to those of slowly pulsating B (SPB) stars. Therefore they have recently been baptised SPBe stars (Walker et al. 2005). Apart from a slightly different effective temperature, HD 181231 and HD 175869 have similar fundamental parameters. In addition, the amplitude of variations detected in HD 175869 is lower than that in HD 181231, as expected because of its later spectral type. Because the pulsations of HD 181231 and HD 175869 have been characterized with great precision thanks to CoRoT by Neiner et al. (2009) and Gutiérrez-Soto et al. (2009) and their fundamental parameters were determined from spectroscopy by the same authors, these late Be stars are very good candidates for seismic modelling. Furthermore, because we know that these Be stars rotate with an angular velocity close to their break-up velocity, they are relevant targets for probing and obtaining tight constraints on the physical processes related to their internal rotation.

We therefore studied HD 181231 and HD 175869 using state-of-the-art modelling tools to unravel rotation processes using seismology. First, we modelled stellar oscillations using two numerical codes (Lee & Baraffe 1995; Clement 1998) that take into account the Coriolis and the centrifugal accelerations, which play a crucial role in pulsation mode dynamics in fast rotators. This allows us to calculate mode frequencies and amplitudes for each given stellar model. We chose to use classical stellar structure models for the target structures, adding a perturbative Chandrasekhar expansion (Lee & Baraffe 1995) to model stars flattened by the centrifugal force, as well as 2D stellar models (Deupree 1990, 1995) that take into account the flattening in a non-perturbative way.

In both cases of models, the transport mechanisms and the associated “non-standard” mixing are not treated in a coherent way but are parametrised using the classical overshoot parameter α_ov that gives the extension of the induced convective core compared to the “standard” model for a single-star evolution without rotation. Thanks to the confrontation of mode computations to the CoRoT observational results, we were able to determine the best value of α_ov for each target. This motivated us to study transport processes that may modify the observed mixing. To this aim, four mechanisms can be considered: the hydrodynamical processes, in particular the meridional circulation and shear-induced turbulence caused by the radiative envelope differential rotation, the possible magnetic field, the penetrative convection at the top of the convective core, and the transport by internal waves. For a complete review on these processes, we refer the reader to Talon (2008) and Mathis (2010).

In Sect. 2 we summarize the pulsational characteristics and fundamental parameters of the two Be stars. We then present the codes and assumptions we used to model the stars (Sect. 3) and the results we obtained for the two stars (Sect. 4). We interpret (Sect. 5) and discuss (Sect. 6) these results in terms of core overshooting and rotational mixing by comparing them with theoretical predictions and models. Finally, we conclude and discuss perspectives of this study.

2. The stars HD 181231 and HD 175869

2.1. HD 181231

HD 181231 is a B5IVe star. Its stellar parameters, assuming a ratio between the angular velocity and its critical value , are summarized in Table 1. See also Table 3 of Neiner et al. (2009).

HD 181231 is not known as a binary in the literature. Inspection of its spectra revealed no sign of a companion. In addition, the spectral modelling performed by Neiner et al. (2009) using solar abundances fitted the observed spectra well and provided no sign of anomalous abundances and therefore no sign of former accretion.

Neiner et al. (2009) detected at least ten independent significant frequencies of variations among the 54 frequencies detected in the CoRoT data of HD 181231, which they interpreted in terms of non-radial pulsation g-modes and rotation. Two longer-term variations were also detected: one at about 14 days resulting from a beating effect between the two main frequencies of short-term variations, the other at about 116 days caused either by a beating of frequencies or by a zonal pulsation mode.

In addition, Neiner et al. (2009) observed HD 181231 with high-precision spectropolarimetry: none of the measurements exibits a magnetic signature, and all longitudinal field measurements are compatible with 0. From the error bars on these measurements, an uppper limit of the possible non-detected longitudinal magnetic field has been set to about 650 G.

2.2. HD 175869

HD 175869 is a B8IIIe star. Its stellar parameters, assuming , are summarized in Table 1 and are very similar to those of HD 181231 except for the cooler temperature. See also Table 3 of Gutiérrez-Soto et al. (2009).

Like HD 181231, HD 175869 is not known as a binary in the literature, its spectra revealed no sign of a companion, and the spectral modelling performed by Gutiérrez-Soto et al. (2009) using solar abundances provided no sign of anomalous abundances and therefore no sign of former accretion.

In HD 175869 Gutiérrez-Soto et al. (2009) detected a frequency of variation compatible with the rotational frequency, as well as six of its harmonics. Low-amplitude periodic variations were also detected. The results were interpreted both in terms of stellar activity just above the photosphere and in terms of pulsation g-modes.

Moreover, Gutiérrez-Soto et al. (2009) observed HD 175869 using high-precision spectropolarimetry: none of the measurements exibits a magnetic signature, and all longitudinal field measurements are compatible with 0. From the error bars on these measurements, an upper limit of the possible non-detected longitudinal magnetic field has been set to about 400 G.

2.3. Evolutionary scenarios

Two evolution scenarios are proposed to explain the existence of Be stars. First, Be stars can be born as single stars with a high initial rotation rate, which evolves as the star progresses on the main sequence (Martayan et al. 2007; Ekström et al. 2008). Second, Be stars can result from the evolution of a B star in a binary system (Pols et al. 1991; Dewi 2007; de Mink et al. 2011). In this case, the high rotation results from an accretion phase from the companion star during the main sequence. It is not straightforward to distinguish a fast rotator that was born as a rapidly rotating single star from a fast rotator that resulted from binary evolution. Indeed, rapidly rotating stars resulting from these interactions will often appear to be single because the companion tends to be a low-mass, low-luminosity star or compact object in a long period orbit. Nevertheless, the lack of obvious surface abundance anomalies in HD 181231 and HD 175869 favours the single-star scenario. Most binary evolution scenarios indeed lead to N surface abundance enrichments (see Langer et al. 2008). In the remain of this paper, we therefore assume that the two stars are single stars. A complete study of the binary scenario is beyond the scope of this paper but we provide an outline in Sect. 6.

3. Models of stellar structure, pulsation, and evolution

Modelling the pulsations of a star first requires the calculation of the star’s structure (see Sect. 3.1) and then a pulsation stability analysis (see Sect. 3.2). In addition we used the Geneva stellar evolution code (see Sect. 3.3) to constrain the impact of rotational mixing on stellar structure.

3.1. Stellar structure

The calculation of the stellar structure of Be stars should take into account the deformation of the star through their rapid rotation. For the present seismic study, we modelled the stellar structure of the rapidly rotating late Be stars with two methods: a 2D model and a Chandrasekhar-Milne-like expansion with a method similar to that of Lee & Baraffe (1995) but slightly modified with the 2D model as input. The non-standard mixing (and transport) processes, which are not treated, are parametrised by an overshoot region at the top of the convective core, whose length is defined as d_ov ≡ α_ovH_p.

3.1.1. 2D rotational deformation: ROTORC

The 2D structure models were calculated with the rotorc code (Deupree 1990, 1995). This code solves the conservation equations for mass, momentum, energy and composition, as well as Poisson’s equation on a 2D spherical grid using a Henyey method. We made the models rotate uniformly on the zero age main sequence (ZAMS), and the interior angular momentum is conserved locally as the model evolves. The independent variables are the fractional surface equatorial radius and the colatitude, and the surface was assumed to be an equipotential. In this way, we needed no assumptions about the shape of the stellar surface, and we obtained the rotationally flattened structure of the star. In this model we used the OPAL opacities (Iglesias & Rogers 1996).

3.1.2. Chandrasekhar-Milne expansion: the Tohoku code

In Lee & Baraffe (1995) an isobaric surface is expressed as r = a [1 + ϵ(a,θ)]  where ϵ(a,θ) is a function of Ω². Here solid-body rotation is assumed (Ω is uniform). Lee & Baraffe (1995) then define ϵ(a,θ) so that the equipotential surface and physical quantities in the equilibrium state depend only on the mean radius a. This allows a mapping between the structure of the rotating and non-rotating stars in terms of a, and ϵ(a,θ) is determined by the Chandrasekhar-Milne expansion (Chandrasekhar 1933a,b). This mapping is only valid when the difference between the rotating and non-rotating structures can be approximated by the Ω² order corrections involving the function ϵ(a,θ), i.e. when the rotation velocity is not too fast. In these conditions, ϵ(a,θ) can be rewritten as ϵ(a,θ) = α(a) + β(a)P₂(cosθ) where P₂ = (3cos²θ − 1)/2. Consequently, r = a [1 + α(a) + β(a)P₂(cosθ)] . The functions α and β represent the horizontal averaged effects of the centrifugal force and the deformation of the equilibrium state, respectively. See Lee & Baraffe (1995) for more details.

Since the late Be stars we are studying here rotate close to their critical velocity, the Chandrasekhar-Milne approximation described above may not be adequate. Therefore we chose to use the 2D structure described in the previous section as an input to our 1D models. This is a first step towards non-adiabatic seismic modelling of rapidly rotating stars (see Sect. 6 for the discussion of the associated drawbacks).

Therefore we assume the isobaric surface to be expressed as , because the unperturbed 1D models already include the effect of the spherical symmetric component of the centrifugal force. The value of for a given rotation frequency is obtained by interpolating the coefficients determined from the 2D ROTORC equilibrium structures. In this model we use OP opacity tables (Badnell et al. 2005). However there is almost no difference in the results when using OPAL opacities.

3.2. Oscillation codes

In the rapidly rotating case, 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 due to the Coriolis and the centrifugal acceleration. For this reason the eigenfunctions for a fixed m value are expanded into a series of terms proportional to spherical harmonics Y_l,m(θ,ϕ) with l ≥ |m|. Indeed, we express the spatial dependence of the displacement vector ξ and a scalar variable f′ as (1)and (2)where with I = 0 for even modes (in which scalar variables are symmetric with respect to the equator) and I = 1 for odd modes. The terms proportional to T^j represent a toroidal displacement needed because of the Coriolis force term in the momentum equation. Obviously a single latitudinal degree l of is not a good parameter anymore for a rotating star. We designate the type of the latitudinal dependence of a mode by using ℓ, which corresponds to the latitudinal degree l of the component with highest amplitude at the stellar surface. This can be best compared to the ℓ values derived from the CoRoT observations.

Furthermore, an excited non-radial mode, which has a frequency ν₀ in the co-rotating frame, should be observed at the frequency |ν₀ − mΩ|, where we adopt the convention that prograde modes have a negative azimuthal order m.

3.2.1. NRO adiabatic oscillation code

Clement (1998) developed an adiabatic code (Normal modes of Oscillation for Rotating stars, hereafter NRO) focusing on g-modes of SPB stars. This code is appropriate for late Be stars, which are located in the same instability strip as SPB stars.

We used a modified version of this NRO code that includes differential rotation terms. The variables used in the models are the radial and latitudinal components of the Lagrangian displacement, the Eulerian pressure, and the gravitational potential perturbations. These variables are evaluated on level surfaces and normalized to have nonzero values at the boundaries. The linearized adiabatic pulsations are solved on a finite difference grid along several radii, i.e. including several spherical harmonics.

Convergence cannot always be reached with the NRO code when the mode frequency is lower than twice the angular velocity of the star (i.e. the inertial frequency), and the modes calculated below 2Ω should therefore be considered with caution. This makes it impossible to compute high-order g-modes (n > 15) in the presence of fast rotation.

3.2.2. Tohoku non-adiabatic oscillation code

To be able to safely model frequencies below 2Ω, we used a non-adiabatic oscillation code, coupled to the Tohoku structure code (see Sect. 3.1.2), to model the non-radial pulsations of the two rapidly rotating Be stars. This code includes the effects of the centrifugal acceleration (which breaks the spherical symmetry of the star as described above) and of the Coriolis acceleration, which have to be taken into account in our seismic modelling since (see Lee 2008; Reese 2010).

With this non-adiabatic oscillation code, we performed a non-radial pulsation stability analysis for rapidly rotating star models, whose evolutionary track passes close to the position of the Be stars in the HR diagram. The evolutionary models were computed in the same way as in Walker et al. (2005) and the analysis allows us to predict the modes that will be excited by the κ-mechanism as well as their frequencies and growth rates (η). The effects of the Coriolis and the centrifugal acceleration are included, although the effect of the deformation is low for low-frequency modes. Moreover, we recall that the rotation is assumed to be uniform.

3.3. Geneva stellar evolution code including transport processes

In the Geneva stellar evolution code, the effects of rotation are included according to the shellular rotation prescription in radiation zones, which assumes a strongly anisotropic turbulence leading to an essentially constant angular velocity on the isobars in these regions (Zahn 1992). This assumes a closure turbulent model for the horizontal turbulence for which we must choose a prescription as for any hydrodynamical model that consistently treats large-scale motions. Here we used prescriptions by Zahn (1992). Moreover, the rotation is assumed to be uniform in convection zones.

With this hypothesis of shellular rotation, every quantity depends solely on pressure and can be split into a mean value and its latitudinal perturbation, i.e. , where P₂(cosθ) is the second Legendre polynomial. Then, the equations of stellar structure need to be modified to account for the hydrostatic effects of rotation because of the centrifugal acceleration (Kippenhahn & Thomas 1970). The usual spherical coordinates are then replaced by new coordinates characterizing the equipotentials. The Kippenhahn & Thomas (1970) method is only applicable in conservative cases, while the problem is non-conservative when the rotation law is shellular. In this case, however, the equations of stellar structure can still be written consistently in terms of a coordinate referring to the mass inside the isobaric surfaces (see Meynet & Maeder 1997). We used the Roche model in the computations. This approximation uses the gravitational potential emerging from the spherical distribution of mass for computing the gravitational force. Note that models including shellular rotation show an internal rotation profile far below the critical limit when they reach the critical velocity at the surface. Consequently, the conditions required for the Roche approximation to be valid are fulfilled in these models. This is not because the Roche approximation is imposed, but because the various physical processes affect the angular momentum distribution in the stellar interior (Meynet et al. 2010).

The transport and mixing processes related to (differential) rotation were treated in the framework of shellular rotation (Zahn 1992; Maeder & Zahn 1998). Then, a large-scale meridional circulation was generated by structural adjustments of the stars and the turbulent viscous transport owing to the secular shear instability. Next, because of the heat advection by this circulation, stellar radiation zones relax onto a new thermal state that generates a new differential rotation, as explained in Rieutord (2006) and Decressin et al. (2009). In this framework, the transport of angular momentum obeys an advection-diffusion equation, while the vertical transport of chemicals through the combined action of vertical advection and strong horizontal turbulence can be described as a pure diffusive process (Chaboyer & Zahn 1992). Because the description of modelling these rotational mechanisms has already been presented in detail in previous papers (see e.g. Meynet & Maeder 2000), we only recall here that the Geneva code includes rotational transport and mixing processes where the meridional circulation is treated as a truly advective process (see Eggenberger et al. 2008). The meridional circulation as well as thermal and chemical fluctuations induced by the 1D shellular rotation are described in 2D using projections on spherical harmonics (see e.g. Meynet & Maeder 2000; Decressin et al. 2009).

4. Seismic modelling of the two late Be stars

4.1. Evolutionary stage

Unperturbed 1D Tohoku models for the pulsation analysis were obtained from spherical symmetric evolutionary models including averaged centrifugal force. Figure 1 shows evolutionary tracks in the log T_eff − log g plane, where solid lines represent evolutionary tracks without convective-core overshooting and dotted lines are those including a core overshooting of 0.35H_p.

In a post-main-sequence star, g-modes are damped by strong radiative damping in the dense radiative core where the spatial wavelengths of g-modes become very short. This means that for g-modes to be excited in a B-type star, the star must be in the main-sequence phase with a convective core. Because HD 181231 and HD 175869 show g-mode oscillations, they should be in the main-sequence phase. Their positions on the log T_eff − log g diagram in Fig. 1 indicate that considerable mixing should be occurring in the radiative layers around the convective-core boundary.

4.2. HD 181231

4.2.1. Tohoku models without mixing

We computed models of HD 181231 with the stellar structure and Tohoku oscillation codes described above. We first calculated models without mixing. We adopted a mass of 5   M_⊙, which is consistent with the spectroscopically estimated mass range (see Table 1) and marginally consistent with the position on the log T_eff − log g diagram shown in Fig. 1. Three values of rotation frequencies (f_rot =  0.432, 0.778 and 1.037 c d^-1) were examined. These values lie within the spectroscopically estimated range of f_rot = 0.76  ±  0.33 c d^-1.

Figure 2 compares frequencies of HD 181231 observed with CoRoT (bottom panel) with frequencies of low-degree (<2) non-radial modes excited in models with log T_eff ~ 4.15 (near the end of the main sequence). The azimuthal order m is also shown for one of the models. It shows that in the most slowly rotating model (f_rot = 0.432 c d^-1), frequencies of excited g-modes cover the observed frequency range but do not show the observed frequency groupings.

Frequency groupings become visible when most of the tesseral g-modes are damped, which occurs when , where R is the mean radius. The critical rotation corresponds to . Since f_rot = 0.432 c d^-1 corresponds to , no frequency groupings occur.

In more rapidly rotating models (f_rot = 0.778 and 1.037 c d^-1), frequency groupings indeed appear because most of the tesseral g modes (ℓ − |m| ≠ 0) are damped. The grouped frequencies, which consist of sectoral g-modes of m =  −1 and  − 2, are too high compared to the observed ones however.

Note that among the lowest frequency around 0.1 c d^-1, some frequencies are combinations of higher frequencies or multiples of the frequency resolution (see Neiner et al. 2009). The others can be modelled with r-modes with m = 1.

To obtain frequency groupings at f_rot ~ 0.4–0.5 c d^-1, we used models with a larger radius (or a lower surface gravity). Since the models adopted in Fig. 2 are close to the termination of the conventional main sequence, we adopted models beyond the main-sequence band by including extra mixing in the radiative layers. This is consistent with the positions of the two stars in the log T_eff − log g diagram (Fig. 1).

4.2.2. Tohoku models with mixing

We then introduced mixing in the models. We calculated models in the effective temperature range log T_eff = 4.155 ± 0.046. Figure 3 shows two models of 5.5 and 6.0 M_⊙, which roughly reproduce the frequency groups around 0.6 and 1.2 c d^-1 with prograde sectoral g-modes having m =  −1 and  − 2, respectively. Both models are near the end of the main-sequence phase with a core overshooting of 0.35 H_p. The positions on the log T_eff − log g plane are shown by open squares in Fig. 1.

The masses used in the two models lie in the upper part and slightly above the limit of the spectroscopically determined range M = 4.4 ± 1.3   M_⊙. While similar fits are obtained with slightly more massive models, no good fits are obtained with M < 5   M_⊙, i.e. our seismic modelling allows us to narrow the possible mass range of HD 181231 down to 5  < M/M_⊙ < 6.5.

The rotation frequencies are 0.475 c d^-1 in the 5.5 M_⊙ model and 0.518 c d^-1 in the 6.0 M_⊙ model, both of which are within the spectroscopically determined range f_rot = 0.76 ± 0.33 c d^-1 and were chosen to qualitatively fit the observed frequency groups. The normalized rotation frequency is 0.52 and 0.62 for the 5.5 and 6 M_⊙ models, respectively. These correspond to  ~75 and 90% of the critical rotation.

The models do not reproduce the comparatively high frequencies observed at 1.87 and 2.5 c d^-1. These observed frequencies are indeed aliases of the main frequency at 0.62 c d^-1, as noted in Neiner et al. (2009).

4.2.3. NRO models

We calculated three NRO models with a full 2D ROTORC stellar structure, with 4, 4.5 and 5 M_⊙ uniformly rotating on the ZAMS with velocities of 200, 250 and 300 km s^-1. The advantage of NRO compared to Tohoku models is that they are fully 2D. The drawback is that the NRO code is adiabatic. It is therefore useful to compare the results obtained with both methods.

Each NRO model was evolved with a convective core overshoot of 0, 0.15 and 0.3 H_p. The evolution sequence can be seen in Fig. 1 of Lovekin et al. (2010). For each model, we then computed a χ² value by comparing the modelled frequencies to the observed ones. We first fitted axisymmetric modes only and then included also non-axisymmetric modes. Details about these models are published in Lovekin et al. (2010).

With axisymmetric modes, the best model has M = 5   M_⊙, is rotating at 105 km s^-1, has an overshoot of 0.3 H_p, and was evolved to an age of 91 Myr. This model has log T_eff = 4.156,R = 5.33   R_⊙ and f_rot = 0.389 c d^-1. While T_eff and R fit the spectroscopically determined parameters well, and M = 5   M_⊙ is compatible with the results on mass obtained with the Tohoku models above, f_rot is lower than the observed f_rot = 0.76 ± 0.33 c d^-1.

Including non-axisymmetric modes, the best-fitting model as M = 5   M_⊙, is rotating at 110 km s^-1, has an overshoot of 0.15 H_p, and was evolved to an age of 84 Myr. This model has log T_eff = 4.154,R = 5.10   R_⊙ and f_rot = 0.426 c d^-1. In this case, f_rot becomes marginally compatible with the spectroscopically determined value.

Figure 4 shows the match between the modelled frequencies and the ten independent frequencies observed in HD 181231 (Neiner et al. 2009) and gives an indication of the ℓ values. From the CoRoT observations, only the ℓ value for the frequency 0.695 c d^-1 could be estimated: ℓ ~ 3 (Neiner et al. 2009). Here we obtain ℓ = 4 from the NRO models for a frequency at 0.708 c d^-1. Because the NRO code is adiabatic, the results show all possible modelled modes without an indication of their excitation or damping. In reality only some of the modelled frequencies would indeed be visible in the CoRoT observations.

4.3. HD 175869

4.3.1. Tohoku models with mixing

We modelled HD 175869, which is slightly cooler than HD 181231, in the same way with the Tohoku code. Reasonable models were searched for in the effective temperature range log T_eff = 4.079 ± 0.022. Figure 5 shows 4.5 and 5.0 M_⊙ models, in which sectoral prograde g-modes of m =  −1 and  − 2 are excited at frequencies around  ~0.55 and 1.3 c d^-1 with a core overshooting of 0.35 H_p. The rotation frequencies of the two models, f_rot = 0.475 and 0.518 c d^-1, correspond to about 70% and 80% of the critical rotation.

The mass range in which similar agreements can be obtained in the above effective temperature domain is 4.5 < M/M_⊙ < 5.5. This range is consistent with, and considerably narrower than, the spectroscopically determined range (4.4 ± 1.3   M_⊙).

In HD 175869, only four of the observed frequencies are independent ones: 0.639, 1.451, 2.0 and 2.475 c d^-1. Other observed frequencies are harmonics or combination of these four frequencies (see Gutiérrez-Soto et al. 2009). The models do reproduce the lower frequencies, but the highest frequencies at 2.0 and 2.475 c d^-1 cannot be explained by the models.

4.3.2. NRO models

The NRO models were calculated for HD 175869 in the same way as presented above for HD 181231. However, with only four independent frequencies and considering the high density of g-modes in these adiabatic models, many models with different parameters but with similar and reasonable χ² values are possible fits of the four frequencies. The NRO models are therefore not well constrained in this case. Nevertheless, we found that all models that fitted the observations with a low χ² require a high overshoot of 0.3 H_p.

4.4. Pulsation g- and r-modes

The theoretical pulsation models of the two late Be stars HD 181231 and HD 175869 that take the stellar deformation into account show that introducing a mixing of α_ov = 0.3–0.35 H_p is necessary to reproduce the observed pulsation frequencies with sectoral prograde g-modes. The value of overshooting α_ov we obtain is a signature of non-standard mixing processes, which are discussed below.

A few higher frequencies are observed in HD 175869 and are not reproduced by the Tohoku models. These modes could be either g-modes that are poorly reproduced by our models because we did not treat rotation perfectly, or g-modes excited by another mechanism than the κ mechanism, e.g. stochastic g-modes.

In addition, very low frequencies ( ≤ 0.2 c d^-1) are observed in both stars and especially in HD 181231. These frequencies could be attributed to m = 1 r modes of high radial order. As Townsend (2005) and Lee (2006) have shown, r modes are indeed excited by the κ mechanism at the Fe opacity bump in B-type stars (which also excites the g-modes). Therefore they appear in the models presented in Figs. 3 and 5.

Indeed, models show r-modes near the zero frequency (<0.15 c d^-1) of m = 1 with ℓ′ = ℓ − 1 = 1 (see Figs. 3 and 5), where ℓ′ is the latitudinal degree for the toroidal component. Frequencies of r-modes in the observers frame can be written as (3)This equation indicates that an r-mode of m = 1 with ℓ′ = 1 (i.e. ℓ = 2) has a low frequency of in the observer’s frame. These modes are odd modes (blue dashed lines in Figs. 3 and 5) because ℓ − m = 1 (the latitudinal degree for spheroidal component determines the symmetry of scalar variables such as temperature perturbations).

5. A laboratory for unravelling transport and mixing processes

As we pointed out in the introduction, four non-standard transport mechanisms may be modifying mixing in the studied stars:

• the hydrodynamical processes induced by rotation: i.e. the action of differential rotation and the associated shear-induced turbulence and large-scale meridional circulation; this combines with the centrifugal flattening due to the centrifugal acceleration;

• the action of a possible fossil magnetic field in the radiative envelope;

• the action of convective flows of the core, which penetrate the surrounding radiative envelope;

• the action of internal waves generated by this penetrative convection.

5.1. The rotational effects

To investigate the effects of (differential) rotation on the properties of the Be stars HD 181231 and HD 175869, and in particular on the size of their convective core, 5 M_⊙ models were computed with the Geneva stellar evolution code described in Sect. 3.3. These models were computed with a solar chemical composition as given by Grevesse & Noels (1993) and a solar calibrated value for the mixing-length parameter (α_⊙ = 1.59 with the input physics used for these computations).

Figure 6 shows the evolutionary tracks in the HR diagram for 5 M_⊙ models with and without rotation. The rotating model had an initial velocity on the ZAMS of 350 km s^-1 and was computed without overshooting from the convective core into the surrounding radiatively stable layers. A corresponding non-rotating model was computed with exactly the same initial parameters except for the inclusion of shellular rotation. The effects of an overshooting from the convective core into the surroundings layers (see Sect. 5.3.1) on a distance d_ov ≡ α_ovmin [H_p,r_core]  (Maeder & Meynet 1989) were also compared by computing the evolution of the same non-rotating model but with an overshoot parameter α_ov of 0.1.

As discussed in detail in Eggenberger et al. (2010) for 3 M_⊙ models, including rotation changes the properties of the star in various ways. First, the hydrostatic effects of rotation result in a shift of the track towards lower luminosities and effective temperatures (compare the dashed and continuous lines in Fig. 6). This is because the effective gravity of the star is reduced by the inclusion of the centrifugal acceleration term. As the evolution proceeds, rotational mixing also begins to play an important role through two main effects. First, this brings fresh hydrogen fuel into the convective core, which slows down its decrease in mass during the evolution on the main sequence. Second, rotational mixing transports helium and other H-burning products into the radiative zone. Including shellular rotation consequently increases the size of the convective core for rotating models compared to non-rotating models. This can be clearly seen in Fig. 6 since this increase results in a shift of the turn-off point at the end of the main-sequence to lower effective temperatures.

Moreover, Fig. 6 shows that including overshooting leads to a similar shift of the turn-off point to the red part of the HR diagram, which is directly related to the increase of the mass of the convective core. Interestingly, this effect is very similar for the rotating model computed without overshooting and for the non-rotating model computed with an overshoot parameter of α_ov = 0.1 (solid and dotted lines respectively in Fig. 6). This indicates that the increase of the mass of the convective core due to rotational mixing is roughly reproduced by a non-rotating model with an overshoot parameter of α_ov = 0.1.

For a 5 M_⊙ model computed with an inital velocity of 350 km s^-1 we find indeed that rotational mixing typically leads to an increase of the size of the convective core of about 0.10 to 0.13 H_p along the main sequence. The exact value of this increase within this range depends on the evolutionary stage of the model because the effects of rotational mixing become more and more visible as the star is close to the end of the main sequence.

Figure 7 shows the enrichment in nitrogen for this rotating model along the main sequence. The enrichment increases as the star evolves, but the final enrichment remains low and probably undetectable in the data mentioned in Sect. 2.

5.2. The contribution of the magnetic field

In the past decade, studies of the magnetic fields in massive stars (e.g. by the MiMeS collaboration, see Wade et al. 2011) have shown that a magnetic field exists in a fraction of these stars. The fields are probably fossil ones in the external radiative envelope of massive stars. The magnetic field is therefore another major mechanism that should be studied in HD 181231 and HD 175869, because it could strongly modify the transport of angular momentum and mixing in these stars (Maeder & Meynet 2004; Mathis & Zahn 2005).

Neiner et al. (2009) and Gutiérrez-Soto et al. (2009) showed that HD 181231 and HD 175869 do not display magnetic signatures in their Stokes V spectra. Moreover, each longitudinal field value determined from their many (18 per star) Stokes V measurements is compatible with 0. The authors estimated upper limits of a possible non-detected longitudinal magnetic field to the largest error bars of the longitudinal field measurements: B_l ~ 650 and  ~400 G for HD 181231 and HD 175869, respectively.

However, the relevant value for the present work is the upper limit of the polar magnetic field at the surface rather than the longitudinal field value, which depends on the magnetic configuration as seen by the observer. In addition, since the rotation period of the stars is not known with a very high precision, we do not know if the 18 Stokes V measurements have been obtained at a magnetic phase where the amplitude of the putative Zeeman signature is strong or flat. We therefore need to perform a statistical analysis.

Taking these facts into account, we performed a Monte Carlo simulation of the spectropolarimetic data to investigate the level of the magnetic field that could be hidden in the data noise.

First, we reanalysed the Narval spectropolarimetric data of Neiner et al. (2009) and Gutiérrez-Soto et al. (2009) using the least-squares deconvolution (LSD) technique (Donati et al. 1997). We used an LSD line mask including 80 and 227 photospheric lines for HD 181231 and HD 175869, respectively. The LSD method assumes that the intrinsic broadening of each line is similar, therefore hydrogen lines are never used in LSD line masks. For each spectral line, the mask contains the wavelength, depth, and Landé factor, to be used by the LSD programme. Landé factors were extracted from Kurucz models. Line depths were adjusted to the observed spectrum. For each spectrum and each star, LSD Stokes I and V profiles were extracted.

For each star and for various values of the polar magnetic field B_pol, we then calculated 1000 oblique dipole models of each of the 18 Stokes V profiles with random obliquity angle, random rotational phase, and noise corresponding to the signal-to-noise ratio of the Narval data.

To calculate these oblique dipole models, we used Gaussian local intensity profiles with a width calculated according to the resolving power of Narval and a typical macroturbulence. The depth of the intensity profile was determined by fitting the observed LSD I profiles. We then calculated local Stokes V profiles assuming the weak-field case and integrated over the visible hemisphere of the star. We obtained synthetic LSD Stokes V profiles, which we normalized to the intensity continuum. We used the mean Landé factor and wavelength obtained from the Narval observations. The model includes three parameters: the inclination angle i, the dipole obliquity angle β, and the polar field strength B_pol.

We then computed among the 1000 models of each 18 Stokes V profiles the chances of detection of the field in at least one of the 18 measurements. This translates into an upper limit for the possible non-detected polar field B_pol = 2300 G for HD 181231 and 700 G for HD 175869, for a 50% chance of detection at 3σ. Pushing the limit to 90% chance of detection at 5σ, a field with any oblique dipolar configuration above 5000 G for HD 181231 and 1500 G for HD 175869 would have been detected (see Fig. 8).

Local stronger magnetic fields could exist if the hypothetical magnetic field were more complex than dipolar. However, for massive (non-Bp) stars, detected fields usually are dipolar (see e.g. Neiner et al. 2003), a predicted behaviour for fossil fields (Braithwaite & Spruit 2004; Duez & Mathis 2010). A few non-dipolar examples are known, such as τ Sco (Donati et al. 2006), but they remain very rare.

The relative action between the magnetic field and (differential) rotation can be calculated by comparing the mean Alfvén pulsation (4)where is the average of the magnetic field amplitude and the mean density, and the characteristic surface angular velocity (5)where R and V_surf are the stellar radius and surface velocity, respectively (Moss 1992). Stellar parameters are similar for HD 181231 and HD 175869. The theory of dipolar fossil fields shows that the maximum internal field amplitude is only about 30 times stronger than the amplitude at the surface (see Fig. 8 in Braithwaite 2008; Braithwaite & Nordlund 2006). Therefore, using the conservative estimate  G, where B_pol = 5000 G is the highest value determined above for HD 181231 and a 90% chance of detection, we find that , i.e. that Ω² is more than six thousand times higher than . Consequently, rotational mechanisms would largely dominate magnetic effects, even if an undetected field were present in the studied stars. For the magnetic field to have an impact, the stars would need to have an averaged magnetic field amplitude  G. Such a mean amplitude is clearly rejected by the Narval observations (even with 100% chance of detection at 5σ) and is unrealistic in a main-sequence Be star. We conclude that in Be stars in general, magnetic field impact on internal mixing will be dominated by rotational processes.

Until now, we assumed that the magnetic field in the radiative envelope is fossil. For completion, below we discuss the possibility for dynamo mechanisms to occur in massive stars. First, Spruit (2002) suggested that a dynamo may develop in radiation zones because of the instability of a predominently axisymmetric toroidal field generated by the shearing of the poloidal field by the differential rotation: this is the so-called Tayler-Spruit dynamo. The only way to close the dynamo loop in this configuration is if the electromotive force from the coupling of the instability fluctuating velocity and magnetic field becomes sufficiently efficient to trigger the axisymmetric poloidal field against its Ohmic decay (Zahn et al. 2007). If this mechanism were to become active in radiative layers, it might modify the mixing (Maeder & Meynet 2004). However, Zahn et al. (2007) and Ruediger et al. (2011) showed how difficult it is to excite this dynamo in stellar interiors and we therefore exclude it as a probable scenario. Second, convective cores of intermediate-mass and massive stars host a dynamo action that generates a magnetic field because of the turbulence and the differential rotation driven by convection (Brun et al. 2005; Augustson et al. 2011) that can also interact with the fossil field at the border of the convective core (Featherstone et al. 2009). The net effect of these processes on the mixing is to modify the convective flows because of the feedback of the Lorentz force on the dynamics and therefore on the penetrative convection, which is discussed hereafter. Finally, the potential small-scale field generated by a dynamo occuring in the sub-surface convective layer (Cantiello & Braithwaite 2011) is not relevant for the large-scale mixing studied here.

5.3. Penetrative convection and internal waves

Here we discuss the extra-mixing rate, which could be induced by penetrative convection at the border of the convective core and the associated related phenomenon such as the generation of internal waves in the surrounding stably-stratified radiative envelope (we ignore here the action of a potential sub-surface convection zone).

5.3.1. Penetrative convection

Penetrative convection, commonly called overshooting in stellar evolution, occurs in stars when convective elements penetrate into adjacent radiative zones. This corresponds here to the turbulent convective flows of the core that penetrate into the radiative envelope due to their inertia (see the pedagogical review by Dintrans 2009). One therefore has to obtain the value of the associated extension of the convective core, which is added to its standard radius determined using the classical Schwarzschild or Ledoux criterion. This cannot be derived using the usual mixing-length phenomenology but requires a physical prescription using theoretical considerations or numerical simulations.

From the theoretical point of view, Roxburgh (1978, 1989, 1992) provided constraints on the size of the convective core and the extent of overshooting. He showed that an estimate of the total size of the convective core can be deduced by looking at the following criterion (see also Zahn 1991) (6)which is deduced from the entropy equation. L_nuc and L_rad are the total nuclear and the radiative luminosities, respectively; Φ is the viscous dissipation per unit volume, is the volume, T is the temperature, r is the radius, and r_c the radius of the convective core including overshooting.

In the nuclear core some of the nuclear energy must be carried by convection, i.e. L_nuc > L_rad. Consequently, if the viscous dissipation is neglected (Φ ~ 0), there must be an overshooting region where L_rad > L_nuc whose extent is estimated from Eq. (6). This criterion has been applied for example by Doom (1985), who claimed that evolutionary models that take into account the Roxburgh criterion represent the observational constraints much better than standard models (see also Meynet et al. 1993, for the impact of overshooting on isochrones). However, convective flows are highly turbulent and the effects of viscous friction are enhanced. Indeed, if the microscopic viscosity of stellar interiors (ν_m) is low, it has to be replaced by a turbulent viscosity ν_t ~ l_tv_t, where l_t and v_t are the characteristic sizes of the largest turbulent eddies. Then ν_t/ν_m ≫ 1 and the adiabatic version of Eq. (6) (Φ ~ 0) gives only an upper limit for the penetration (Zahn 1993).

However, the necessary characterisation of turbulent flows can only be tested using numerical simulations in the domain of parameters accessible today. Moreover, we have to take into account rapid rotation, which is one of the key parameters for our two rapid rotating Be stars.

This is the reason why it is interesting to discuss the most recent global highly non-linear numerical simulations of the dynamics of the rotating convective cores in intermediate-mass and massive stars (see e.g. Browning et al. 2004; Meakin & Arnett 2007; Augustson et al. 2011). If we focus on the detailed work by Browning et al. (2004), many interesting features can be isolated concerning convective cores, which are illustrated in their work devoted to the case of an A-star. First, for a given value of the entropy gradient at the convection/radiation border, the degree of turbulence (i.e. the Reynolds number) influences the penetration length, which decreases as the flow is more and more turbulent because of a smaller filling factor for the convective upflows (the plumes) (see e.g. the 3D non-linear simulations of penetrative compressible convection in localised planar domains, Brummell et al. 2002). Second, when the rotation rate is increased (for given stratification and diffusivities), the penetration becomes more intense because of the stabilisation of the flow by the Coriolis acceleration. This behaviour has also been predicted in planetary studies of, e.g., Takehiro & Lister (2001). Finally, Browning et al. (2004) showed that the 0.2 H_p extension needed here to complete the rotational transport contribution to match our seismic modelling is attainable.

Because of rapid rotation, penetrative convection therefore constitutes the natural candidate to explain the additional H_p extension that cannot be explained by the combined action of centrifugal acceleration and secular rotational transport computed in Sect. 5.1. However, to be able to conclude, one has to examine the last transport mechanism that may operate in the studied stars, namely the action of internal waves.

5.3.2. Internal waves

Internal waves are generated by penetrative convective flows at the top of the convective core as observed in the previously cited simulations (Browning et al. 2004; Meakin & Arnett 2007). These internal waves can induce mixing in two ways:

• First, internal waves propagate in the radiative envelope and deposit their associated angular momentum at the location where they are damped by radiative damping. There, the differential rotation is modified because of Reynolds stresses, and mixing is modified in an indirect way because of the shear instability (Zahn et al. 1997; Talon & Charbonnel 2005).

• Second, internal waves act directly on the chemical distribution as a diffusion mechanism, as proposed by Press (1981) and Schatzman (1996), or because they propagate in their non-linear regime close to their excitation region, which makes them break (Mocák et al. 2011).

The Be stars we study here rotate close to their break-up velocity. The Coriolis acceleration (and the centrifugal one) therefore modifies the propagation of internal waves as soon as the 2Ω/ω_c ratio increases (ω_c is the frequency at which the internal wave is excited). The induced transport of angular momentum and mixing is modified in two ways:

• First, internal waves are trapped in an equatorial belt when ω_c < 2Ω because of the Coriolis acceleration action (Lee & Saio 1997; Mathis 2009; Ballot et al. 2010). The convective kinetic energy is then transmitted less efficiently to the internal waves in the radiative envelope. The associated transport is consequently decreased (see Pantillon et al. 2007; Mathis et al. 2008).

• Next, because of the Coriolis acceleration, internal waves are damped through thermal diffusion closer to the excitation region than in the slowly rotating case (i.e. close to the convective core) and just below the surface (see again Pantillon et al. 2007; Mathis et al. 2008).

Therefore it is reasonable to believe that transport and mixing processes by internal waves in extremely rapidly rotating stars, such as those studied here, may be confined close to the core or the surface and be less efficient than for slowly rotating stars.

6. Discussion

Our seismic models have shown that mixing of 0.3–0.35 H_p is necessary to reproduce the pulsation g-modes observed in the CoRoT data of the two late Be stars HD 181231 and HD 175869. We also reproduced r-modes at low frequencies.

Our calculation of stellar evolution models including rotation showed that 0.1 to 0.13 H_p of the non-standard convective core extension is due to the combined action of the centrifugal flattening and the rotational mixing mechanisms that lead to the increase of its size. Our present knowledge of penetrative convection and transport processes induced by internal waves suggests that the other 0.2 H_p non-standard core extension may be caused by the convective overshoot with a possible weak action of excited internal waves. Linearly combining the rotational mixing and overshoot contributions is allowed because of their different characteristic time-scales (secular and dynamical, respectively). Finally, we showed that a potential fossil magnetic field does not significantly contribute to mixing in the studied stars and in Be stars in general.

The drawbacks of our present analysis are the following:

• First, the current knowledge on internal waves generation and related transport processes in rapidly rotating stars is quite poor. Indeed, no secular evolution models, such as those applied to low-mass stars by Talon & Charbonnel (2005), are available in the literature for massive stars. Moreover, the study of the impact of the Coriolis and centrifugal accelerations on transport is only in its infancy, as is the prediction of the stochastic excitation of gravito-inertial waves.

• Second, we combined results obtained with several tools that have been computed in 1-, 2-, and 3D. To obtain a coherent picture in the near future, a more rigourous approach would be to combine completely 2D stellar evolution codes that include transport processes (Rieutord et al. 2005; Espinosa Lara & Rieutord 2007) with corresponding high-precision 2D stellar oscillation codes (Reese et al. 2009). However, this constitutes an extremely complex task, and the approach presented here, in which state-of-the-art models are used, constitutes a good first exploration of highly coupled dynamical processes induced by rapid rotation. Moreover, 3D simulations of penetrative convection that have been discussed in Sect. 5.3.1 do not take into account the stellar flattening by the centrifugal acceleration. However, this effect becomes important only for external layers and can be neglected for the convective core.

• Finally, we assumed a single-star evolution scenario. If the studied Be stars had followed a binary evolution, their structure and their mixing would have been modified by the accretion and tidal interactions (the dynamical tide, see e.g. Zahn 1975, 1977). These effects should then be coupled to previously studied transport processes, which are also efficient in binaries (Zahn 1994; Talon & Kumar 1998; de Mink et al. 2009). This constitutes one of the main challenge for the modelling of non-standard mechanisms in stellar evolution.

Therefore, even though we used several state-of-the-art tools that should be improved in the near future for extremely rapidly rotating stars, we conclude that the observed mixing is probably caused by the combined effect of strong rotational mixing owing to the rapid rotation of Be stars and convective overshooting with a possible contribution of excited internal waves. In this context, asteroseismology of critically rotating stars constitutes a challenge for the development of an improved modelling of stellar transport and mixing processes and gives us strong related constraints.

7. Conclusions

Using various state-of-the-art codes to model CoRoT observations of pulsations in two rapidly-rotating late Be stars, we have shown that these stars host prograde sectoral g-modes as well as possible r-modes.

With these models and theoretical calculations, we have also shown that strong mixing is present in these stars. Part of the mixing is due to penetrative convection at the top of the convective core with a possible contribution of excited internal waves. The rest of the mixing is due to rotational processes, in particular the meridional circulation and shear-induced turbulence due to the radiative envelope differential rotation and the centrifugal flattening.

Finally, we have shown that if a magnetic field is present in rapidly rotating stars, such as Be stars, its impact on mixing will be dominated by the impact of rotation.

Acknowledgments

The CoRoT space mission, launched on December 27, 2006, has been developed and is operated by CNES, with the contribution of Austria, Belgium, Brazil, ESA, Germany and Spain. We thank J.-P. Zahn for valuable discussions and the referee for his/her suggestions. C.N. thanks the Agence National pour la recherche (ANR) programme SIROCO, the Programme National de Physique Stellaire (PNPS) and CNES for their support. C.N. is grateful to E. Alecian for providing her Monte Carlo code.

References

All Tables

All Figures