Issue 
A&A
Volume 553, May 2013



Article Number  A25  
Number of page(s)  14  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/201220559  
Published online  25 April 2013 
Populations of rotating stars
II. Rapid rotators and their link to Betype stars
^{1}
Geneva Observatory, University of Geneva,
Maillettes 51,
1290
Sauverny,
Switzerland
email:
anahi.granada@unige.ch
^{2}
Astrophysics group, EPSAM, Keele University,
LennardJones Labs,
Keele, ST5 5BG, UK
^{3}
Centre de Recherche Astrophysique, École Normale Supérieure de
Lyon, 46 allée d’Italie, 69384
Lyon Cedex 07,
France
^{4}
Masaryk University, Kotlářská 2, 61137
Brno, Czech
Republic
^{5}
Bartol Research Institute, University of Delaware,
Newark,
DE
19716,
USA
Received:
15
October
2012
Accepted:
20
February
2013
Context. Even though it is broadly accepted that single Be stars are rapidly rotating stars surrounded by a flat rotating circumstellar disk, there is still a debate about how fast these stars rotate and also about the mechanisms involved in the angularmomentum and mass input in the disk.
Aims. We study the properties of stars that rotate near their criticalrotation rate and investigate the properties of the disks formed by equatorial mass ejections.
Methods. We used the most recent Geneva stellar evolutionary tracks for rapidly rotating stars that reach the critical limit and used a simple model for the disk structure.
Results. We obtain that for a 9 M_{⊙} star at solar metallicity, the minimum average velocity during the mainsequence (MS) phase to reach the critical velocity is around 330 km s^{1}, whereas it would be 390 km s^{1} at the metallicity of the Small Magellanic Cloud (SMC). Red giants or supergiants originating from very rapid rotators rotate six times faster and show N/C ratios three times higher than those originating from slowly rotating stars. This difference becomes stronger at lower metallicity. It might therefore be very interesting to study the red giants in clusters that show a large number of Be stars on the MS band. On the basis of our singlestar models, we show that the observed Bestar fraction with cluster age is compatible with the existence of a temperaturedependent lower limit in the velocity rate required for a star to become a Be star. The mass, extension, and diffusion time of the disks produced when the star is losing mass at the critical velocity, obtained from simple parametrized expressions, are estimated to be between 9.4 × 10^{12} and 1.4 × 10^{7} M_{⊙} (3 × 10^{6} to 4.7 × 10^{2} times the mass of the Earth), 2000 and 6500 R_{⊙}, and 10 and 30 yr. These values are not too far from those estimated for disks around Betype stars. At a given metallicity, the mass and the extension of the disk increase with the initial mass and with age on the MS phase. Denser disks are expected in lowmetallicity regions.
Key words: stars: general / stars: evolution / stars: rotation / stars: emissionline, Be / stars: massloss
© ESO, 2013
1. Introduction
There is general consensus that classical Be stars are intrinsically rapidly rotating Btype stars, with analogs among O and early A types (although they are much less frequent, see e.g. Zorec et al. 2005), surrounded by a gaseous decretion viscous disk. Such a disk, geometrically thin and in nearly Keplerian rotation, can explain many of the characteristics that define Be stars (see for example the reviews by Porter & Rivinius 2003; and Carciofi 2011).
However, there are still many open questions regarding Be stars: do they reach the critical limit, or do they rotate at 70–80% of their critical limit? Does binarity play a role in the appearance of the Be phenomenon? Which are the mechanisms involved in the angularmomentum and mass input in the disk? As noted by Stee & Meilland (2009), it seems that even if we consider the single stars that present the Be phenomenon, they do not all share a common massejection mechanism.
From an observational point of view, some authors (Frémat et al. 2005; Cranmer 2005, 2009; Huang et al. 2010) found that Be stars, in particular earlytype Be stars, seem to be rotating below their critical limit, though in some cases the determination of rotational velocities could be underestimated (Townsend et al. 2004). Also, Rivinius et al. (2001) found for some earlytype Be stars such as μ Cen that the beating of multiperiodic nonradial pulsations seems to be linked with the mass outbursts. Nevertheless, the lack of observational evidence for largeamplitude pulsations or other strong mechanisms that could help to drive angular momentum and mass into a circumstellar disk in many of the stars, favours the idea of weaker processes that would be able to overcome the effective gravity and let the matter escape from the photosphere into the disk: thus rotation close to the critical limit is expected (e.g. Howarth 2007). This idea is also supported by many interferometric observations, which have allowed to determine the oblateness of some Be stars (Domiciano de Souza et al. 2003; Zhao et al. 2011; Touhami et al. 2011; Meilland et al. 2012; Kraus et al. 2012), suggesting that some of them are close to their critical rotation.
Interestingly, theoretical studies on stellar evolution show that along the mainsequence (MS) evolution, the single intermediatemass stars with a sufficiently high rotational velocity on the zeroage main sequence (ZAMS) can reach equatorial velocities near the critical value^{1}, thanks to the transfer of angular momentum from the inner contracting parts to the outer regions (Meynet et al. 2007; Ekström et al. 2008).
Ekström et al. (2008) studied how stars with different masses, different metallicities, and different initial rotational velocities evolve towards the critical limit. In that work, the computation stopped when the model reached the critical velocity. In the study we present here, the following changes and improvements with respect to Ekström et al. (2008) were made:

The numerical method for computing the surface velocity asresulting from the various transport processes and mass loss wasimproved (as explained in detail in Georgyet al. 2013, hereafter Paper I). In contrast to the computation by Ekströmet al. (2008), in which the deformation of the star and the angular momentum contained in theenvelope^{2} wereneglected, we account here for both effects. The envelope contains one tenthousandth of the total stellar mass during the MSphase. We assumed that this region has the same angular velocityas the first layer of the interior^{3}. This development allows for a much better tracing of theevolution of the total angularmomentum content of the star. Asa result, our models now reach the critical limit less easily than inEkström et al. (2008), as discussedin Sect. 3. In addition, the effects of anisotropicwinds were accounted for in the computations following Georgyet al. (2011).

The computations were now pursued all along the MS phase until the helium flash for stars with M_{ini} ≤ 2 M_{⊙}, the early asymptotic giant branch for 2.5 M_{⊙} ≤ M_{ini} ≲ 9 M_{⊙} or the end of coreC burning when M_{ini} ≳ 9 M_{⊙}, even for the models that reach the critical limit. This allows us to explore the consequences of fast rotation beyond the MS phase.

In the present work, the physical ingredients are the same as in the paper by Ekström et al. (2012), and differ from those used by Ekström et al. (2008), as discussed in Sect. 3.
Using this new version of the Geneva stellar evolution code, our group has produced grids of stellarevolution tracks at solar metallicity (Z = 0.014) in the mass range of 0.8 to 120 M_{⊙}, both rotating and nonrotating (Ekström et al. 2012), and showed that the models rotating with V_{ini}/V_{crit} = 0.4 provide a reasonable fit to important observed characteristics of stars. In Paper I we presented a grid to explore in detail the mass domain between 1.7 to 15 M_{⊙}, which corresponds to earlyA, B, and lateO spectral types, at different metallicities and with initial angular velocity between 0 < Ω_{ini}/Ω_{crit} < 0.95. In this last paper, the effects of rotation on the tracks in the HR diagram, the lifetimes, the surface velocities, and abundances were discussed.
Here we specifically focus on a subset of these models: those whose surface velocity at the equator reaches the critical limit during the MS phase. With the study of these models, we aim to address the following questions: what can be said about the mechanical massloss rates when the star reaches the critical limit? What would be the extension and lifetime of a stationary circumstellar viscous decretion disk formed by matter expelled from the star? What is the massloss rate expected from the stardisk system? How do the properties of the disk vary as a function of the mass, metallicity and age of the star? To find the answer, we used the analytical development in Krtička et al. (2011) together with the results presented in Paper I.
The paper is organised as follows: the main physical ingredients of the stellar models are briefly recalled in Sect. 2. In Sect. 3, we discuss the properties of models that reach the critical velocity during their MS phase: average values for the rotation velocity, properties of red giants and supergiants originating from very rapid rotators, and the frequency of rapid rotators in clusters of various ages. We discuss the mechanical massloss rates obtained when the models reach the critical limit in Sect. 4. Section 5 discusses the properties of the disks around models that rotate at the critical limit using the disk model described in Krtička et al. (2011). The possible relations between stars at the critical limit and classical Betype stars are discussed in Sect. 6. Finally, the main conclusions are presented in Sect. 7.
2. Stellar models
The stellar models analysed here have masses between 1.7 and 15 M_{⊙}, and metallicities of Z = 0.002 (Z_{SMC}), 0.006 (Z_{LMC}), and 0.014 (Z_{⊙}). They are a subset of the grid presented in Paper I, since we consider here only the models that reach the critical limit during the MS phase, i.e. models with ω_{ini} = Ω_{ini}/Ω_{crit} ≳ 0.80. A detailed description of the physics included in the models can be found in Paper I. We briefly recall a few relevant inputs:

Overshooting: the convectivecore baundaries are determined with the Schwarzschild criterion. The convective core is extended with an overshoot parameter d_{over}/H_{P} = 0.10 (after calibration, as explained in Ekström et al. 2012).

Radiative mass loss: the more massive models (M > 5 M_{⊙}) lose mass through linedriven winds, according to the massloss prescription by de Jager et al. (1988). Only stars with an initial mass greater than 12 M_{⊙} lose a nonnegligible amount of mass through radiative winds during the MS phase. For instance, our 15 M_{⊙} stellar model at Z = 0.014 loses about 10% of its initial mass during the MS phase through radiative winds. For the redgiant phase, the mass loss rates are obtained by the prescription of Reimers (1975, 1977) up to 12 M_{⊙}, whereas the de Jager et al. (1988) prescription is used for masses of 15 M_{⊙} or larger.

Mechanical mass loss: when overcritical velocities are reached at the equator of a star, some matter is released and presumably launched along Keplerian orbits in a disk around the star. In the following, we call such an event an overcritical episode^{4}. Because of numerical difficulties in computing of the stellar evolution when the model is at the critical limit, we need to define a value for the highest fraction of the critical angular velocity allowed, ω_{max}, above which the mechanical mass loss is switched on. In our computation, we set ω_{max} = 0.99. The quantity of mass ejected mechanically from the model is governed by the quantity of angular momentum that needs to be removed from the surface to maintain the model below the limit set by ω_{max}, as explained in Paper I. The mass injection into the disk allows the star to decrease its surface velocity below the critical value, until the secular evolution (and more precisely the transport of angular momentum by meridional circulation and shears) drives the surface to become overcritical again. Thus, once the critical rotation is reached, the surface velocity of the model will stay around the critical value until the end of the MS, i.e. until the evolutionary timescale becomes too short to allow the transport of angular momentum to the surface.

Rotation: the transport of angular momentum and chemical species inside a star is implemented following the prescription of Zahn (1992) for the horizontal diffusion coefficient, and that of Maeder (1997) for the shear diffusion coefficient.
3. Properties of stars at the critical limit
Fig. 1 Evolution of the Ω/Ω_{crit} ratio for 9 M_{⊙} at Z = 0.002, with initial values of Ω_{ini}/Ω_{crit} = 0.10, 0.30, 0.50, 0.70, 0.80, and 0.90 (bottom to top) during the MS phase. The solid lines show the present models, the dotted lines the models by Ekström et al. (2008). 
In Fig. 1, we show the evolution of the ω = Ω/Ω_{crit} ratio for the 9 M_{⊙} models at the SMC metallicity, together with the corresponding models of Ekström et al. (2008, dotted lines).
We recall that we start the computation of a model with solidbody rotation on the ZAMS, similarly to the work by Heger et al. (2000). The underlying assumption is that preMS stars are fully convective, and that convection drives rigidbody rotation. Once the solidbody constraint is released, the flat internal profile of rotation is modified very rapidly under the action of meridional currents, and the surface rotation decreases strongly, as can be clearly seen in Fig. 1. This means that, in the present scenario, the models begin their evolution with a slowdown of their surface rotational velocity, and may reach the critical limit only after some evolution on the MS. As a result, the present way of initiating the computation makes it impossible to obtain very rapidly rotating stars at a very young age, even when imposing a high value of ω_{ini}. The validity of the assumption made is now being explored with models of preMS evolution and will be the subject of a future publication. The first results show that the internal rotation profile depends to some extent on the way the star mass and angularmomentum content build up through accretion. At the moment, though, we take the present results at face value, keeping in mind that we might be missing some young critical rotators.
From Fig. 1, we also see that the critical velocity is reached at the end of the MS phase. This occured in the models by Ekström et al. (2008) and is confirmed by the present models, although significant differences between these two sets of models appear for the most rapidly rotating models. The change in the prescription for the diffusion coefficient explains the slight changes in the MS lifetimes, but has no impact on the evolution of the surface velocity. This is clearly shown in the work by Meynet et al. (2013), who compared the timeaveraged velocities during the MS phase obtained using different sets of diffusion coefficients, D_{shear} and D_{h}. These authors found that for a 15 M_{⊙} model at Z = 0.002 with V_{ini}/V_{crit} = 0.5 the differences caused by using the present prescription and that of Ekström et al. (2008) are very small, of the order of only 0.5%. The source of the difference in surface velocities for the most rapidly rotating models is the improved numerical procedure for computing the surface velocity. The present models maintain a lower surface velocity because we now account for the angular momentum in the envelope and for the deformation of the star, and we also ensure the angularmomentum conservation of the star+wind system (including radiative and mechanical mass loss). As a result, the present models evolve less easily towards the critical limit, while the models by Ekström et al. (2008) overestimated the surface velocity.
The main properties of the models reaching the critical velocity at their surface during the MS phase are summarised in Table 1.
From Table 1, we see that only models with ω_{ini} > 0.80−0.90 are able to reach the critical velocity during the MS phase. If we compare these limits with those given by Ekström et al. (2008), we see that for the metallicity and the mass range in common between the two works, namely for the 3 and 9 M_{⊙} at Z = 0.002, the present models require a minimum ω_{ini} = 0.95 , while the models by Ekström et al. (2008) require ω_{ini} = 0.90 and 0.80, respectively. This reflects the fact mentioned above that the present models approach the critical limit less easily. We also see that the lowermass models reach the critical rotation much more easily than the more massive ones and can spend up to 30% of the total MS time at the critical velocity. This trend was also obtained in Ekström et al. (2008) at Z = 0.002. This comes from the longer time available in lowermass models, for meridional circulation and shears to bring angular momentum at the surface. Moreover, in the higher initialmass models, mass loss by stellar winds may remove some angular momentum from the surface and thus keep the star away from the critical limit (although this effect is quite weak at low metallicity).
Data for the models that reach the critical velocity during the MS.
3.1. Surface velocities and abundances for very rapid rotators
For the 9 M_{⊙} models at Z_{⊙}, the minimum average velocity needed to reach the critical velocity is around 330 km s^{1}. At Z_{SMC} this minimum value is higher by about 20% (~390 km s^{1}). For the 1.7 M_{⊙} models, the minimum average velocity amounts to 240 and 270 km s^{1} at Z_{⊙} and Z_{SMC}.
At the end of the MS phase, quite large enhancement factors of the N/C abundance ratios are obtained at the surface of the initiallyrapid rotators. While at Z_{⊙} it amounts to a factor of 2 for the 1.7 M_{⊙} model, it amounts to a factor of 20 or more for the 9 M_{⊙} model. The corresponding enhancement factors at Z_{SMC} are between 20–40 for the 1.7 M_{⊙} model and around 500 for the 9 M_{⊙} model.
It is interesting to note that when these very rapid rotators evolve into the redgiant or redsupergiant stage, they keep some traces of their former very high rotation rate (see Tables 2–4 in Paper I). This is reflected in two observable quantities: 1) their rotation rate is higher, typically the model of 2.5 M_{⊙} at Z = 0.014 with ω_{ini} = 0.10 has a rotation period P_{rot} = 950 days at the end of the core Heburning phase, while the corresponding model with ω_{ini} = 0.95 has P_{rot} < 150 days at the same evolutionary stage, i.e. it rotates at least six times faster; 2) the N/C ratio is enhanced in formerly rapid rotators, typically the 2.5 M_{⊙} originating from the ω_{ini} = 0.10 model has a N/C ratio of 1.4 (in mass fraction) at the end of the core Heburning phase, while the corresponding model with ω_{ini} = 0.95 has N/C = 4.8, i.e. more than three times higher. These differences are stronger at lower metallicities. Typically, the rotation for the 2.5 M_{⊙} at Z = 0.002 at the end of the core heliumburning phase is enhanced by a factor of 7 between the ω_{ini} = 0.10 and 0.95 models. The N/C ratio is more than 12 times higher in the initially most rapidly rotating model than in the slower one.
From the above results, we see that it might be very interesting to study red giants in clusters that show a large number of Be stars on the MS band. From the present models, one expects that the red giants that evolved from the Bestar population present surface velocities and N/C ratios higher than normal.
Timestep dependence of the total mass lost mechanically and of the mean mechanical massloss rates (in the 250 000 yr interval).
Properties of the 7 M_{⊙} stellar models at Z = 0.014 computed with various values of ω_{max}.
Mean mass fluxes during the criticalrotation phase for stars with Ω_{ini}/Ω_{crit} = 0.95.
3.2. Expected number of rapidly rotating Btype stars in clusters as a function of age
The fraction of stars with surface velocities above a given limit varies in clusters of various ages, and it can be estimated from the computed grids. This is shown in Figs. 2 and 3 for solar and SMC metallicity. The populations were computed as follows. First, the mass range (1.7 M_{⊙}–15 M_{⊙}) and initial rotationrate range (0 ≤ ω_{ini} ≤ 1) were discretised in several intervals (1000 mass intervals and 100 velocity intervals), dividing the plane M_{ini} − ω_{ini} into 100 000 cells. For each cell, we extracted the expected number of stars (normalised to an initial population of 1), according to an initial mass function (Salpeter 1955), and to an initial rotationrate distribution (Huang et al. 2010). Each cell thus defines initial conditions in Z, M, and ω. We built an evolutionary track for these initial conditions interpolating between our computed tracks. At each time step, we checked the spectral type of the star in each of the cells, as well as its current rotation rate. If this rotation rate was higher than a given constant (0.70, 0.90, or 0.98), we added the number of stars in the cell to the total number of stars of a given spectral type that rotates more rapidly than the given constant. The population was normalised with respect to the total number of B stars (including Be stars) at the same time^{5}. We assumed that the distribution of initial velocities given by Huang et al. (2010) corresponds to the distribution of velocities on the ZAMS, i.e. before the slowing down due to the activity of the meridional currents described above in Sect. 3. We emphasise here once more that an appropriate account of the way the mass and angularmoment content of the star builds up during the accretion phase may change this picture. This will be studied in detail in a forthcomming paper (Haemmerlé et al., in prep.).
Fig. 2 Evolution of the fraction of stars that rotate faster than a given limit in clusters of various ages at Z = 0.014, taking the Huang et al. (2010) distribution of velocities as the initial velocity distribution. The observational points represent the observed Be/(B+Be) ratio (McSwain & Gies 2005). The squared symbols (solid error bars) are for confirmed Be stars, and the triangle symbols (dashed error bars) are for all Betype stars. The ticks at the top of the plot show the typical lifetime of stars with various initial masses. 
Fig. 3 Same as Fig. 2 for Z = 0.002. Observed Be/(B+Be) ratio extracted from Martayan et al. (2010), following the same procedure as in McSwain & Gies (2005). We consider the following age beams: 7 ≤ log (t) < 7.5, 7.5 ≤ log (t) < 8.0, 8.0 ≤ log (t) < 8.5, and log (t) ≥ 8.5. The error bars for the fraction of Betype stars are Poisson statistic error bars. 
Under the previous hypotheses we find that at Z_{⊙}, the fraction of stars with ω > 0.98 (very rapid rotators hereafter, blue solid line in Fig. 2) becomes greater than 1% after 100 Myr, when the initially rapidly rotating stars with masses between 4 and 5 M_{⊙} approach the critical limit. The fraction of very rapid rotators then increases with the age of the cluster until a maximum at around 4% is reached at an age of about 300 Myr.
The curve for ω > 0.90 (green dashed) presents a local minimum around 20 Myr. This reflects the fact that in the first part of the evolution of our lower initial mass model, the ω ratio decreases. Therefore this first behaviour is quite dependant on the initial conditions chosen. Then the curve follows a similar trend as the one for ω > 0.98, but shifted at higher values. In contrast, the line corresponding to ω > 0.70 (red dotted) is almost constant up to about 400 Myr.
Comparing the values obtained at Z = 0.014 (Fig. 2) and 0.002 (Fig. 3), a few interesting remarks can be made:

for most of the age interval between 30 and 400 Myr, the fraction of stars with ω > 0.98 at a given age is greater at Z_{⊙} than at Z_{SMC}. This occurs because 4 and 5 M_{⊙} models at Z_{SMC} do not attain values of ω close to 1 (see also Table 1): even though the rotational velocities do not differ significantly between SMC and solar metallicities for the mentioned masses, these models do not approach the critical limit because the stars are more compact and their critical velocities are higher at SMC metallicities.

Stars with ω > 0.98 are found in older clusters at Z_{SMC} than at Z_{⊙}. This is because at SMC metallicity, at ages older than 1 Gyr only the lower mass rapidly rotating Btype models are still in their MS. At this age and metallicity, all other models with masses higher than 1.7 M_{⊙} and those of the same mass but with medium or low rotation, will have ended their corehydrogen burning. The enhancement in the MS lifetime for 1.7 M_{⊙} rapidly rotating models occurs because their central temperature attains values for which the ppchain is almost entirely responsible for the energy production during a long part of the MS phase (see discussion and Fig. 5 in Paper I).

For stars with ω > 0.90, the differences between Z_{⊙} and Z_{SMC} are less marked, except that at the Z_{SMC} some older clusters show a significant fraction of ω > 0.90 stars and that at very high ages (above the billion years), a great fraction of the Btype stars may be rapid rotators. This is expected because only the most rapidly rotating models are still Btype stars on the MS band at these ages. They are some kind of blue stragglers^{6}.
Implications of these results concerning the question of the Be phenomenon are discussed in Sect. 6.
4. Mechanical mass loss at the critical limit
In this section, we explain how the mechanical mass losses (see Cols. 9 and 10 in Table 1) were computed. We want to stress that in the frame of the hypotheses described below, we obtain an outbursting behaviour for the mass loss. However, this outbursting behaviour results from our numerical procedure and does not come from ab initio physical principles. In this section we explain in detail how the mechanical mass losses were computed and study how the amount of mass lost depends on the time step and on ω_{max}.
The general principle is the following: owing to internal transport of angular momentum and to the stellar evolution itself, it may occur, depending on the initial rotation, that the stellar surface reaches (and even surpasses) the critical velocity (or rather ω_{max} in our computation, see Sect. 2) during the evolution. In this case, the excess angular momentum is removed from the outer layers in the form of an additional equatorial mass loss on a dynamical timescale, i.e. within one time step. This brings the surface velocity of the star to below 0.99 Ω_{crit}. Ideally, the stellar surface should rotate at 0.99 Ω_{crit}, but it does not, because of discretisation. The transport of angular momentum by meridional currents and shears will then accelerate the surface until the limit is crossed again and some more mass is shed. For illustration purposes, the evolution of the mass loss obtained this way is shown in Fig. 4 during a period spanning about 1000 years for our 7 M_{⊙} stellar model at solar metallicity. We see a very regular pattern where the quiescent phase between two massejection episodes lasts for about 40 years, and the mass ejections are of the order of 10^{8}M_{⊙} yr^{1}. For this model, the radiative winds are negligible and thus the only massloss mechanism taking place is mechanical mass loss, as described above.
In the following, we explore the dependency on the results on the time step used for the calculations and on the parameter ω_{max}.
Fig. 4 Evolution of the mechanical massloss rate as a function of time for the 7 M_{⊙} stellar model at solar metallicity with a time step of about nine years (case n = 105, see text). Each point represents the mass loss rate in one time step. The quiescent phase between two massejection episodes lasts for about 40 years. 
4.1. Discussion of the choice of the time step
First, we briefly examine how the choice of the time steps modifies the amplitude of the mechanical massloss rates. Depending on 1) the duration of the time steps; 2) the equatorial velocity of the star just before becoming overcritical; and 3) the efficiency of the angularmomentum transport, the surface velocity during the next time step may surpass the critical limit to a lesser or greater extent. In some cases, very weak mechanical massloss episodes will occur (typically when the amplitude of the evolution of V_{eq} beyond the critical limit will just graze the limit), while in some other cases, very strong massloss episodes will occur, when V_{eq} strongly surpasses the limit. Therefore our numerical method leads to some scatter in the massloss rates. In Fig. 4, the interval shown is too short for this scatter to appear.
To study the impact of the choice of the time step on the outputs of our models, we computed different 7 M_{⊙} models at Z = 0.014 with ω_{ini} = 0.95. The first column of Table 2 indicates the value of n, characterising the different time steps Δt = Δt_{standard}/n, where Δt_{standard} = 942 yr. The last two columns show the total mass lost and the mean massloss rate during the considered period, which in the present case corresponds to a period of 250 000 yr beginning when the age of the star is 5.21 × 10^{7} yr. When the time step becomes lower than about 100 years, these two quantities stabilise around 7.3 × 10^{4} M_{⊙} for the total mass lost mechanically during that period and towards 2.9−3.0 × 10^{9} M_{⊙} per year for the timeaveraged massloss rates. This shows that below a given value for the time step, the results are no longer varying. We also see that taking a time step greater by an order of magnitude does not change the results by an order of magnitude: the difference amounts at most to a factor of two in the range of the time steps studied here.
Fig. 5 Fraction of time spent with a given massloss rate for the 7 M_{⊙} models at Z = 0.014 with various time steps. 
Figure 5 shows the fraction of time spent near the critical limit with zero and nonzero massloss rates for our 7 M_{⊙} model at Z = 0.014. Depending on the time step used, the fraction of the time spent away from the critical limit, i.e. when the star is rotating fast but does not lose mass through mechanical mass loss, amounts between 50 and 80% of the critical period. Shorter time steps tend to increase the duration of the null mechanical massloss period and to increase the massloss rates in the period when mass ejections occur. The differences become small when time steps shorter than 100 yr are used.
As a conclusion, we see that the scatter brought to the results by a change in the time step is not dramatic. The adoption of the longest time steps does not change the results by more than a factor of two with respect those obtained with the shorter time steps, and also allows for a reasonable computational time for these extended grids.
To support these numerical results, we compared them with a simple analytical estimation of the mechanical mass loss rate. Asuming that the star rotates as a solid body, that the rotational velocity change is very small when the star rotates close to the critical limit, and also that the time variation of the stellar angular momentum of the star is due to the loss of matter through the equator, we obtain simply that . This shows that the mechanical mass loss rate when the star rotates at the critical limit depends basically on the evolution of the momentum of inertia mass and the equatorial radius, which both depend on the stellar mass. With this approximation, within the same interval of time as described above, we obtain a mass loss rate of 1.5 × 10^{9} for a 7 M_{⊙} model at Z = 0.014, which nicely agrees with the numerical results.
4.2. Choice of the value of ω_{max} at which mechanical mass loss occurs
As emphasised above, because of the numerical difficulties of computing the stellar evolution when the star is at the critical limit, we need to define a value for the highest fraction of the critical velocity allowed, which we call ω_{max} (see Sect. 2) above which the mechanical mass loss sets in.
To check to what extent the choice of ω_{max} affects the evolution of the star, we computed 7 M_{⊙} models at solar metallicity with ω_{ini} = 0.95, setting ω_{max} equal to 0.99, 0.95, and 0.90 just after the readjustment phase at the very beginning of the MS. The results are presented in Table 3. The first column gives ω_{max}, the second and third one the MS lifetime and the fraction of the MS lifetime spent near the critical limit; the total mass lost and the timeaveraged massloss rate over the critical period are given in Cols. 4 and 5.
Fig. 6 HR diagram for 7 M_{⊙} models at solar metallicity computed with different values of ω_{max} (blue: 0.99, green: 0.95, and red: 0.90). As a comparison 5 and 9 M_{⊙} evolutionary tracks are also plotted. 
We see that the lower ω_{max}, the longer are the MS lifetimes, and the higher it is the ratio τ_{crit}/τ_{MS}. This is a consequence of the surface braking by the mechanical mass loss, which occurs earlier for the models with the lower values of ω_{max}, producing a stronger internal mixing in the star. This mixing brings fresh hydrogen to the core, increasing its size, slightly increasing its luminosity, and also the MS lifetime. The total mass lost and the mechanical massloss rates are also higher for lower ω_{max}. This is expected, since a lower value of ω_{max} implies an earlier entrance into the critical regime, more mass removed when the limit is crossed, and a higher amount of total mass lost. Typically, passing from ω_{max} = 0.99 to 0.95 already increases the total mass lost by more than a factor of two. Because both the total mass lost and the duration of the critical period increase in a similar way when ω_{ini} decreases, the averaged mass loss rate remains nearly constant.
Figure 6 shows the HR diagram for the 7 M_{⊙} models computed with the different values of ω_{max}, and as a comparison the evolutionary tracks corresponding to the nonrotating case and the evolutionary tracks corresponding to ω_{max} = 0.99, for 5, 7 and 9 M_{⊙}. The differences in log (L) and log (T_{eff}) have different causes. As mentioned above, the luminosity shift is due to the larger core of the model with the lower ω_{max}. The T_{eff} shift has a geometrical origin. For a given luminosity, the surface of a star rotating with ω = 0.99 is 1.16 times larger than the surface of a star with ω = 0.90, because of the deformation of the shape of the surface (see Georgy et al. 2011). Because the mean T_{eff} is deduced according to the StefanBoltzmann law (with Σ the stellar surface), this produces a shift of 0.016 dex towards the red in the HR diagram.
In the rest of the paper, we discuss the results obtained with ω_{crit} = 0.99. To conclude this section we stress two points:

1.
Reality may be quite different from the situation described above. It might occur, for instance, that mechanisms other than critical rotation are involved in producing mass ejection. These additional mechanisms may allow the loss of matter well before the surface rotational velocity becomes critical.

2.
While the present approach is of course probably too schematic, it can at least provide a reasonable estimate of the mass lost mechanically when the following hypotheses are made: mechanical mass loss only occurs when the velocity is higher than a fixed limit and the quantity of mass lost is entirely determined by the necessity of the star to shed overcritically rotating material.
4.3. Radiative and mechanical mass loss
In Table 1, Col. 10, we show the massloss rates averaged over the period during which the star evolves near the critical limit for stars of different stellar masses and metallicities.
The total mass and angular momentum lost mechanically are quite small with respect to the total mass and angularmomentum content. Typically, for the 15 M_{⊙} model at Z = 0.002, the total mass lost mechanically corresponds to 0.005% of the initial mass. This means that this massloss mechanism has only very little influence on the evolution of the star.
For a given metallicity, the mean mechanical massloss rates are higher for higher stellar masses, As mentioned in Sect. 4.1, this dependence can also be deduced from simple analytical calculations: the massloss rates depend on the evolution of the stellar momentum of inertia and on the equatorial radius of the stars, and both quantities depend on the initial stellar mass. For a given mass, we find that the mechanical massloss rates are in general higher at lower Z. In the high mass star range, this can be explained easily: indeed, in the highmass range, the stronger radiative winds of higher Z remove mass and angular momentum at the surface, keeping the star more easily away from the critical limit and thus reducing the need for a mechanical mass loss. In the lowmass range, the behaviour is more complex, which is the result of counteracting effects. First, stars are more compact at low Z, so the extent of the region over which angular momentum must be transported is smaller; second, this compactness leads to weaker meridional currents (Maeder & Meynet 2001), which reduces the transport efficiency. The present results seem to imply that in the metallicity range presented here, the first effect overcomes the second. From the above, we conclude that our mechanical mass loss rates increase with mass, as is the case of mass loss rates via radiatively linedriven winds, but mechanical mass loss is higher at lower Z, in contrast to what happens for stellar winds.
An interesting point to be highlighted here is the impact of the mass loss on the angularmomentum content of the stars. The total initial angular momentum of models that rotate with different ω_{ini} remains within an order of magnitude for a given stellar mass. Along the MS evolution, the angular momentum of the star will be modified by the effect of stellar winds (for stars with masses larger than 7 M_{⊙}) and by the mechanical angularmomentum loss, for stars that reach the critical limit. The stellar winds dominate the highmass range, driving the loss of about 30% of the initial angular momentum at Z = 0.002. The mechanical mass loss reduces the angular momentum at most up to 10%, but is very effective in injecting angular momentum into the circumstellar environment, probably forming a circumstellar disk.
4.4. Polar and equatorial mass fluxes
If mass loss occurs only via radiatively driven winds, the equatorial to polar massflux ratio depends only on the value of ω (Maeder 2002; Georgy et al. 2011) and basically reflects the temperature contrast between the pole and equator described by the von Zeipel’s theorem. For a star rotating at ω = 0.99 the equatorial to polar massflux ratio is therefore predicted to be at most 0.15, if there is no mechanical mass loss (Georgy et al. 2011). We see in Table 4 that for the models with 7 and 9 M_{⊙}, the mean mass flux is higher at lower metallicities. This might be surprising at first sight since a stronger mass loss by radiative winds is expected at high metallicity because of the metallicity dependence of the linedriven winds (which scale as ~Z^{0.5}). However, polar winds occur only when the surface velocity is close to the critical limit and the evolutionary stage when this occurs depends on the metallicity, and at low metallicity the critical limit is reached closer to the end of the MS, where the radiative winds are stronger due to the higher luminosity of the star.
When mechanical mass loss occurs, the mass flux at the equator consequently becomes much greater than the mass flux at the pole (enhancement factors of up to 60) in most cases. The only exceptions occur for the highest initial masses in each metallicity. In that case, the radiative winds are stronger and thus remain of the same order of magnitude (or even stronger in one case) than the equatorial mechanical mass flux.
5. Decretion disks
5.1. Estimation of the physical properties
Struve (1931) first proposed the existence of a decretion disk formed from an equatorialmass ejection episode by a rapidly rotating star. Since then, the existence of Keplerian rotating disks has been increasingly supported (e.g. Porter & Rivinius 2003; Carciofi 2011).
In the previous section, we showed that the surface velocity remains at around the critical limit from the first moment it reaches it until the end of the MS phase, and that mass ejections in the frame of the hypotheses made here occur as intermittent outbursts. The ejected mass launched into an equatorial disk around the star explains the observed Balmer emission lines of Be stars.
We study here the physical characteristics of a stationary disk that can be deduced from the developments made by Krtička et al. (2011), who solved the hydrodynamic equations in cylindrical coordinates to obtain the disk structure. Krtička et al. (2011) provided parametrised expressions for the massloss rate as a function of the angularmomentumloss rate and for the outer radius of the disk. We used these simple parametrised expressions to study how the disk properties change along the criticalrotation phase.
There are three free parameters in this disk model: two of them, p and T_{0}, arise from the assumption that the temperature in the disk varyies as a powerlaw in radius, T = T_{0}(R_{eq}/r)^{p}, where R_{eq} is the equatorial radius of the star. The third free parameter is the viscosity term α of Shakura & Sunyaev (1973), which intervenes in the momentum equation (tangential component). This viscosity couples the inner parts of the disk to the central star and may be generated by the presence of small magnetic fields (Balbus & Hawley 1991).
The structure of the disk and our stellar model are connected through three variables: 1) , which is the quantity of angular momentum that the star must lose per unit of time to maintain the surface velocity at the critical limit (at ω_{max} = 0.99); 2) R_{eq}, the equatorial radius of the star; and 3) V_{K}(R_{eq}), the Keplerian velocity at the equator of the star (see more details below).
Part of the excess angular momentum that the star has to lose is lost through radiative stellar winds. To calculate the angularmomentumloss rate by this process, , we inserted the radiative massloss rates Ṁ given by our evolutionary models and assumed that this mass is lost isotropically at the mean stellar radius given by , which gives . By doing so, we neglected the wind anisotropies and thus slightly overestimate the loss of angular momentum due to this process. However, except for the 15 M_{⊙} stellar model at Z = 0.002, the stellar winds remove a negligible fraction of the mass and of the angular momentum even when the effects of the stellar winds are overestimated. The remaining angularmomentumloss rate, , is transferred to the disk ( hereafter for clarity).
In a stationary situation, the disk will lose angular momentum at the same rate. Therefore, knowing how and where the disk loses this angular momentum will set up the disk structure corresponding to the considered physical situation. The disk can lose angular momentum through two different channels:

1.
through disk winds (). The value for the maximum is obtained from expression (26) in Krtička et al. (2011), which corresponds to the highest diskwind angularmomentumloss rate^{7}. In all our models the highest is negligible with respect to (compare Col. 5 with Col. 3 in Table 5), therefore we neglect the effect of a disk wind hereafter;

2.
through disk mass loss (). Since this is the only active process in most cases, one has .
The quantity is given by our evolutionary calculations, while the value of is obtained by Krtička et al. (2011). These authors showed that at a radius R_{crit} the radial velocity of the matter in the disk becomes equal to the local sound velocity, , where μ is the mean molecular weight and m_{H} the mass of a hydrogen atom. Matter flowing beyond R_{crit} removes angular momentum from the stardisk system.
Accounting for the fact that at large distances the disk is not rotating at the Keplerian velocity, where Ṁ_{disk} is the mass lost by the disk, V_{K}(R_{crit}), the Keplerian velocity at R_{crit}, i.e. . The equality between and allows us to determine Ṁ_{disk}, provided R_{crit} is known.
Krtička et al. (2011) proposed an approximation of the critical radius, which can be considered as the outer radius of the disk: (1)where a(R_{eq}) is the sound speed evaluated at R_{eq}.
At this point, we have three interesting quantities that characterise the disk: its radial extension, R_{crit} (defined here as its extension up to the point when angular momentum is no longer conserved), the rate of angularmomentum loss, , and of mass loss, Ṁ_{disk}.
Main characteristics of the decretion disks for models with Ω_{ini}/Ω_{crit} = 0.95.
We easily obtain two more interesting quantities, the diffusion time and the mass of the disk. The characteristic timescale for the evolution of the disk is given by the viscous diffusion time, (2)where is the diffusivity of the disk, with the viscosity parameter, and H the disk scale height, H = a(R) R^{3/2}/(G M)^{1/2}, where G is the gravitational constant and M the mass of the star. Evaluating H at R_{crit} and assuming an isothermal disk leads to (3)From the diffusion timescale and massloss rate of the disk, we can estimate the mass of the decretion disk. It is obtained by multiplying the massloss rate of the disk by the diffusion timescale: (4)The mass of the viscous decretion disk does not depend on the outer radius R_{crit} but is set by the angularmomentuminput rate in the disk when the stellar surface reaches the critical limit.
The disk properties are computed assuming an isothermal disk (p = 0) of temperature T_{0} = 2/3 T_{eff} (which are the disk parameters frequently used in the literature to represent Bestar disks), and α = 1 as proposed by Carciofi et al. (2012).
5.2. Results from the stellar models
In Table 5, we give the properties of the disks that are obtained by applying the equations of the previous section to our stellar models with ω_{ini} = 0.95. The indicated quantities are obtained by computing their timeaveraged value throughout the critical period to remove the scatter due to the outbursting nature of the process.
We note the following points:

The massloss rates via the viscous decretion disk, ⟨ Ṁ_{disk} ⟩, are much lower than the massloss rates given in Table 1, Ṁ_{mech,mean}, corresponding to the mass lost at the equator of the star: in the disk, the matter is released at a farther distance, and thus a lower mass loss is needed to remove a given amount of angular momentum. The massloss rates from the disk are about a factor of 10 below the massloss rates at the surface of the star (Krtička et al. 2011).

The masses of the disks are quite small. Values between 9.4 × 10^{12} and 1.4 × 10^{7} M_{⊙} (3 × 10^{6} to 0.047 times the mass of the Earth) are obtained. These quantities are strongly dependent on the disk parameters p, T_{0}, and α.

Under the hypotheses considered, the extension of the disk amounts to a few thousands of solar radii (between about 10 and 30 Astronomical Units). Expressed in terms of stellar radii, the extension amounts to some hundreds of stellar radii.

The diffusion time is very short, only a few decades. This is much too short for any instability in the disk to set up and to form something like an asteroid belt!

At a given metallicity, the mass and the extension of the disk increase with the initial mass (see upper panels of Figs. 7 and 8). Typically at solar metallicity, the mass of the disk for the 9 M_{⊙} model is a factor of 1600 higher than that of the 1.7 M_{⊙} model, while the extension of the disk is increased by a factor of 2.2 only. Therefore, the disks are expected to be on average denser around more massive stars.
Fig. 7 Evolution of the radius and diffusion time of the disk as a function of time for critically rotating models with an initial mass between 4 and 9 M_{⊙} at solar metallicity.

For a given initial mass, the extension of the disk increases when the star evolves along the MS as shown in Figs. 7 and 8. The same is true for the disk diffusion time.

For a given initial mass, in general the mass of the disk slightly increases when the metallicity decreases and the extension of the disk slightly decreases. For instance, the disk of the 9 M_{⊙} model is 2.8 times more massive at Z = 0.002 than at Z = 0.014, while the extension of the disk, in contrast, is 16% smaller at Z = 0.002 than at Z = 0.014, i.e. some denser disks are expected in lowmetallicity regions. Overall, Table 5 shows that the differences due to the metallicity are small compared to those due to the initial mass.

The disk diffusion time does not present very strong variations either as a function of the initial mass or as a function of the metallicity.
5.3. Sensitivity to the disk parameters
It is important to check to what extent the choice of the parameters that describe the thermal structure of the disk and the viscosity parameter affects our results. We use the data of the 7 M_{⊙} model at solar metallicity, two different values for the parameters T_{0}/T_{eff} (0.66 and 0.80), for p (0.0 and 0.1), and for α (0.1 and 1.0). The results are summarised in Table 6.
A scatter is produced in the mass, extension, diffusion time, and massloss rate of the disk. The highesttolowest value ratio varies by a factor of 33 for the mass, of 2.3 for the extension, of 20 for the diffusion time and of 1.4 for the massloss rate. This illustrates the absolute values obtained in Table 5 are sensitive to the physical ingredients of the disk models.
The mass and the diffusion time of the disks are the most sensitive quantities, and α, the viscosity, is the most critical quantity. A stronger viscosity implies smaller masses for the disks and shorter diffusion times.
Mean values of some disk’s characteristics for the 7 M_{⊙} model at solar metallicity, varying T_{0}/T_{eff}, p, and α.
6. Possible link between very rapid rotation and the classical Betype star phenomenon
6.1. Is critical rotation a necessary and sufficient condition for Betype stars?
There is no doubt that Be stars are rapid rotators, i.e. with ω > 0.70, however, it is not clear whether Be stars are rotating at the critical limit or not. In the sample studied by Levenhagen & Leister (2004), eight Betype stars with a mass estimate between 8 and 10 M_{⊙} show Vsini values between 160 and 335 km s^{1}, with a mean value of 230 km s^{1}. Corrected for a random orientation of the rotational axis, this value translates into a mean V_{eq} around 290 km s^{1}. The critical velocity for our 9 M_{⊙} stellar model at the end of the MS phase (when the most rapidly rotating models reach the critical limit) is 410 km s^{1}. Comparing these two values, we obtain a mean ratio of V/V_{crit} = 0.71 (i.e. ω = 0.88) for the observed stars, which means that Betype stars may not necessarily be at the critical limit. The same conclusion can be deduced from the work of Frémat et al. (2005), who found a mean ratio of ω = 0.88 from the fit of photospheric lines of Be stars, accounting for gravity darkening, which also agrees with the result of Levenhagen & Leister (2004), who did not account for it. Beshell stars are assumed to be Be stars seen equatoron. Thus, the determination of Vsini for these objects gives direct information on V_{eq}. Rivinius et al. (2006) showed that for the brightest shell stars of different spectral types, the mean fraction of the critical velocity is 81%, which corresponds to a mean ω higher than 0.92. Meilland et al. (2012) found a mean value of ω = 0.95 for a sample of eight Be stars, considering Vsini and V_{crit} from Frémat et al. (2005) but i measurements derived from interferometry. Delaa et al. (2011) deduced from interferometric observations of Be stars 48 Per and Ψ Per that both stars rotate at nearly critical rotation. Therefore, the question to which extent Be stars are really critically rotating stars is not settled yet. Two additional factors have to be accounted for in that discussion: 1) when the surface velocity approaches the critical limit, the method of velocity measurement based on the observed width of the absorption lines might underestimate the values (Townsend et al. 2004); 2) it may be that Be stars are not always critically rotating, and undergo periods in the course of their evolution during which they rotate below the critical value and periods when they are rotating critically.
In Figs. 2 and 3, we show the fraction of Be stars observed in clusters with ages between 10 and 300 Myr at solar metallicity (data from McSwain & Gies 2005) and between 10 Myr and 1 Gyr at the SMC metallicity (data extracted from Martayan et al. 2010), and compare them to the ratio given by our models. Considering only single stars with ω > 0.98, the numbers are much too low to account for the observed fraction of Betype stars, particularly for young clusters. This last conclusion agrees with the study by Ekström et al. (2008). Stars with ω ~ 0.80 cannot reproduce the observed trend with age either. Then, in the frame our hypothesis, our present models do not provide a reasonable fit to the observations.
A different conclusion is obtained when one considers that the lowest velocity that a star has to reach to present the Be phenomenon may vary as a function of the stellar mass and/or the temperature. This was suggested by Cranmer (2005) and Huang et al. (2010), who found that more massive (or hotter) stars become Betype stars at lower V/V_{crit} values than lowermass (or cooler) stars. For instance, if we consider the lowest rotation rate values (V_{eq}/V_{crit})_{min} in bins of T_{eff} proposed by Cranmer (2005) to obtain a Be star, then the fractions of stars with velocities higher than these values (which appear as Be stars) given by our models reproduce the observed fraction of Be stars quite well (see Fig. 9). Our models accordingly support the analysis by Cranmer (2005).
Fig. 9 Evolution of the Be/(B+Be) ratio in clusters of various ages at Z = 0.014 obtained from our models assuming a mass (or T_{eff}) dependence for the Be phenomenon as proposed by Cranmer (2005, using (V_{eq}/V_{crit})_{min}. The observational points are the same as in Fig. 2. 
In addition to the variation of the fraction of Be stars with age, another important feature is the way this fraction varies as a function of metallicity at a given age. It has been established (Maeder et al. 1999; Wisniewski & Bjorkman 2006, see also our Figs. 2 and 3) that there is a larger fraction of Be stars at low metallicity. Indeed, the observed fraction of Betype stars in clusters with ages between 10 and 25 Myr is higher by about a factor of 2 at SMC metallicity (Be/(Be+B) fraction of about 30%) than at solar metallicity (15%). Starting from identical initial velocity distribution, we did not obtain such an enhancement in our models. Solutions of this problem can come from different directions. For instance, it may be that the initial rotation distribution is shifted to higher velocities at low metallicity. Another possibility could be that since disks are denser at low metallicity (for a given mass and age), they could be more easily detected at low metallicity. According to Silaj et al. (2010), the minimum density for a disk to produce Hα emission is 5 × 10^{13} g cm^{3}. Finally, the analysis of Cranmer (2005), which sets a minimum velocity for becoming a Betype star that varies as a function of the effective temperature was performed at solar metallicity. At present we cannot say whether such an analysis would hold at lower metallicities and whether this would help in resolving the above question.
On the basis of our singlestar models, we therefore conclude that there are much too few stars that rotate strictly at the critical limit to account for the observed Betype star population. If the Be phenomenon appears at lower rotation rate, as suggested by many authors, and has a dependence on the mass and/or the T_{eff}, as proposed by Cranmer (2005) and Huang et al. (2010), the models reproduce the fraction of Betype stars much better. As a consequence, some other mechanism is required in addition to the centrifugal force to push the matter into a disk.
This conclusion of course depends on the physics included in the models. Here, we did not account for any internal dynamo, which would impose a solidbody rotation during the MS phase and in turn would generate higher surface rotations during the MS phase. Moreover, in addition to the mechanism explored here, where the surface is accelerated by internal processes, some stars may acquire rapid rotation by tidal interactions or through mass exchange in close binaries. Furthermore, as already discussed previously, we started the computation on the ZAMS with a flat Ω profile, which implies an internal readjustment period at the very beginning of the evolution. As a result, the effective ω_{ini} is lower than the one by which we label our models, and we do not produce critical rotators on the ZAMS.
6.2. Comparison of the massloss rates
Fig. 10 Mass loss – luminosity diagram. The colour points correspond to data by Waters et al. (1987), while the black crosses correspond to data from Stee (2003) and the black circle to the data from Carciofi et al. (2012) for 28 CMa. The thin dotted curve traces the upper limit for IR massloss rates obtained by Waters et al. (1987). The thick black curve corresponds to the mean mechanical massloss rates along the critical rotation phase, obtained as described in Sect. 4 and shown in Table 1, while the dashed black line corresponds to the mean value during the time in which the star is effectively losing mass. The grey region corresponds to the range of instantaneous massloss rates. 
Waters (1986) and Waters et al. (1987) studied the density structure and massloss rates of Be stars and interpreted the observed IR excesses in terms of a simple poleon disk model, which allowed them to estimate upper limits for the IR (equatorial) massloss rates between 7 × 10^{9} and 2 × 10^{8} M_{⊙} yr^{1}, and masses around 1−5 × 10^{8} M_{⊙}. The thin dotted line in Fig. 10 shows this upper limit for IR massloss rates obtained by Waters et al. (1987), and the colour points in this figure show some of their data for which reliable stellarmass determinations are available in the literature. Stee (2003) also obtained massloss rates of some Be stars with SIMECA, a latitudedependent radiative wind plus dense disk model that produces both spectroscopic and interferometric synthetic data, plotted in crosses in Fig. 10. Carciofi et al. (2012) obtained similar values for 28 CMa by modelling the decline rate of the Vband excess. Since the values by Waters et al. (1987) and Stee (2003) are strongly model dependent, they cannot be taken as strong constraints of the models. However, it is interesting to compare them with the mass loss rates we obtained in the present work. In Fig. 10 we also plot the mean equatorial massloss rates obtained from our criticallyrotating models as described in Sect. 4 and presented in Table 1 as a solid black line for stars with Ω/Ω_{crit} = 0.95 on the ZAMS. In this plot we also represented the averaged massloss rates for overcritical episodes (black dashed line). As mentioned in Sect. 5, the disk massloss rates are lower than these values.
We see that our massloss rate estimates run below the determinations made by Waters et al. (1987). The origin of this shift between the two sets of values is difficult to assess, since both may be model dependent. However, it is remarkable that the models predict an increase of the mechanical loss rate with the luminosity similar to that obtained by Waters et al. (1987).
In this context, it is interesting to note that per se, an overcritical velocity is not sufficient to drive a mass loss. When it becomes overcritical, matter will be launched into a Keplerian orbit and therefore will remain gravitationally locked to the star (Jupiter is orbiting around the Sun with a critical velocity!). This matter may even fall back on the star, depending on the forces acting on it. Some additional force (or forces?) have to come into play and provide the momentum needed to drive an effective mass loss from the star and/or an outwardsexpanding disk. Here we did not account for any additional forces. If they were present, they would probably increase the mass “mechanically” lost by a fastrotating star, and therefore our present estimate must be considered as a lower limit.
Looking at Fig. 10, we can see that the upper values obtained in our models for the instantaneous mass loss rates fall in the range of the observed values. However, while the timeaveraged mass loss rates are meaningful quantities in the frame of our hypothesis, this is not the case of the instantaneous values, which are quite model dependent. Therefore, we cannot deduce any firm conclusion from this coincidence.
An interesting point is that stars rotating at or near the critical velocity may or may not have an equatorial disk, but all these stars should present some degrees of anisotropy in the radiative winds. As mentioned in Sect. 4.4, the radiative mass flux ratio between the equator and the pole reaches 0.15 at most. When mechanical mass loss occurs at the equator, the ratio increases to 60. Some observations have been performed on α Ara (a 9−10 M_{⊙} star estimated to rotate at V_{eq}/V_{crit} = 0.91, i.e. ω = 0.98−0.99) and α Eri (Achernar, a 6−8 M_{⊙} star estimated to rotate at ω = 0.992). For α Ara, Chesneau et al. (2005) found a polar flux of 1.7 × 10^{9} M_{⊙} yr^{1} sr^{1}, with a ratio between the polar and equatorial flux Φ_{eq}/Φ_{pol} = 30. These numbers have been revised in a second paper (Meilland et al. 2007) to Φ_{pol} = 7 × 10^{9} M_{⊙} yr^{1} sr^{1} and Φ_{eq}/Φ_{pol} = 0.03. For α Eri, Kanaan et al. (2008) found a polar flux of 4.0 × 10^{12} M_{⊙} yr^{1} sr^{1} and a ratio Φ_{eq}/Φ_{pol} = 0.01−1. The observational range is wide, and for the time being, it seems to be premature to use these kind of observations to put constraints on the models.
6.3. Comparison of the disk mass, extension, and diffusion time
In our model, we obtained that the masses of disks increase significantly with the stellar mass. Interestingly enough, Silaj et al. (2010) modelled Hα emission line profiles of Be stars with their nonLTE radiative transfer code BEDISK by fitting disk parameters: their Table 3 shows that disks around B0 type stars are more massive than disks around B8 type stars.
Our models predict extensions for the disk’s outer radius of some thousands of solar radii. Radio studies provide information on the outer regions of disks surrounding Be stars. Using the Very Large Array, Taylor et al. (1990) carried out a radio survey of 21 Be stars selected for their excess IR emission detected by IRAS. They found that the minimum radius of Bestar disks are of several hundreds to thousands of solar radii. Identical results were obtained by Clark et al. (1998), who observed a sample of 13 Be stars of the southern sky with the Australian Telescope Compact Array. The values of a few thousand solar radii we obtain with the different stellar masses and different disk parameters seem to be compatible with the observed quantities. The extension of the disk may be dependent on the wavelength used to investigate it. In Hα, the extensions are much shorter, a few tens of solar radii, since only the portion which is nearby the star and therefore hot and dense enough would be observable, while at longer wavelengths, cooler and shallower parts of the disk can be seen. Our predictions are valuable for isolated stars. Disk truncation can happen if a close companion star is present.
The diffusion timescales we obtain with our models can be compared to longterm variability cycles observed in Be stars. Some of the observed longterm variations in Be stars, consisting of transitions between normal B and Be phases, or brightness variations (see e.g. Porter & Rivinius 2003; Carciofi 2011, and references therein), can be understood in terms of the formation and dissipation of a circumstellar disk. They typically last from several months or years to decades. The values of a few of decades we obtain are consistent with these timescales, though different choices of disk parameters can affect this quantity.
7. Conclusions
With an improved version of the Geneva stellarevolution code, we have studied the characteristics of initially very rapidlyrotating stars as well as some properties of the equatorial disks that are produced when the surface velocity becomes overcritical.
The main conclusions are:

The angularmomentum and mass losses (between 4 × 10^{13} and 4 × 10^{9} M_{⊙} yr^{1}, depending on the initial mass and metallicity) when stars are at the critical limit are very modest and have no strong impact on the way the star evolves.

Nevertheless, signatures of rapid initial rotation on the ZAMS remain present in more evolved stages such as the redgiant or supergiant stages, for instance the higher values for rotation and surface nitrogen enrichment compared to stars that evolved from slowly rotating progenitors.

For stars that reach the critical limit, the timeaveraged massloss rates we obtain are one order of magnitude lower than those currently estimated from observations of Betype stars. However, among the range of instantaneous massloss rates we obtained, the highest values agree relatively well with the estimates derived from observations.

When the star is losing mass at the critical velocity, the mass of the disk produced is estimated to range between 9.4 × 10^{12} and 1.4 × 10^{7} M_{⊙} (from 3 × 10^{6} to 47 × 10^{2} M_{⊕}), its extension between 2000 and 6500 R_{⊙}, and its typical evolution timescale between 10 and 30 yr. Some hints were given on how these values vary as a function of the initial mass, age, and metallicity. With our simple assumptions, even though the massloss rates are somewhat lower than those derived from observations, the masses, extensions, and diffusion times of the disk obtained here are consistent with the current estimates of these quantities for observed disks around Betype stars.

If we assume that there exists a T_{eff} dependence for the minimum ω for the Be phenomenon to occur (as suggested in the literature), the models agree reasonably well with the observed values of Be/(B+Be).
Our results were obtained in the frame of the shellular rotation model. If an internal magnetic field would more strongly couple the core to the envelope, and provided no substantial magnetic braking occurs at the surface, the star may reach (at an earlier evolutionary stages) the critical limit more easily. In a forthcoming paper we will explore the consequences of a strong interior coupling.
By critical velocity we mean here the velocity required for the centrifugal force to counterbalance the gravity at the equator of the star (see e.g. Maeder & Meynet 2000).
It is important to properly distinguish between V/V_{crit} and Ω/Ω_{crit}. Both quantities indicate how close to its breakup limit a star rotates, but they are not equivalent. Writing the angular velocity as Ω = V_{eq}/R_{eq}, where V_{eq} is the rotational velocity at the equator and R_{eq} is the equatorial radius of the star, it turns out that Ω/Ω_{crit} = (V_{eq}/V_{crit})·(R_{eq,crit}/R_{eq}). Because R_{eq} depends on the rotational velocity of the star, Ω/Ω_{crit} is equal to V_{eq}/V_{crit} only when Ω = 0 or in a star rotating at the critical limit R_{eq,crit} = R_{eq}. Throughout this work, every time we refer to rotation close to the critical limit, we implicitly assume that we talk about Ω/Ω_{crit}.
If this highest were be greater than the angular momentum carried by the disk outflow, , the extension of the disk would be set by the radius R_{out}, given by Krtička et al. (2011) such that . This never happens in our models, see text.
Acknowledgments
The authors thank the referee of this article for his/her constructive comments. C.G. acknowledges support from EUFP7ERC2012St Grant 306901. J.K. was supported by grant GA ČR 205/08/0003.
References
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Carciofi, A. C. 2011, in IAU Symp. 272, eds. C. Neiner, G. Wade, G. Meynet, & G. Peters, 325 [Google Scholar]
 Carciofi, A. C., Bjorkman, J. E., Otero, S. A., et al. 2012, ApJ, 744, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Chesneau, O., Meilland, A., Rivinius, T., et al. 2005, A&A, 435, 275 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Clark, J. S., Steele, I. A., & Fender, R. P. 1998, MNRAS, 299, 1119 [NASA ADS] [CrossRef] [Google Scholar]
 Cranmer, S. R. 2005, ApJ, 634, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Cranmer, S. R. 2009, ApJ, 701, 396 [NASA ADS] [CrossRef] [Google Scholar]
 de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259 [NASA ADS] [Google Scholar]
 Delaa, O., Stee, P., Meilland, A., et al. 2011, A&A, 529, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Domiciano de Souza, A., Kervella, P., Jankov, S., et al. 2003, A&A, 407, L47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ekström, S., Meynet, G., Maeder, A., & Barblan, F. 2008, A&A, 478, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [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]
 Georgy, C., Meynet, G., & Maeder, A. 2011, A&A, 527, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Georgy, C., Ekström, S., Granada, A., et al. 2013, A&A, 553, A25 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Howarth, I. D. 2007, in Active OBStars: Laboratories for Stellare and Circumstellar Physics, eds. A. T. Okazaki, S. P. Owocki, & S. Stefl, ASP Conf. Ser., 361, 15 [Google Scholar]
 Huang, W., Gies, D. R., & McSwain, M. V. 2010, ApJ, 722, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Kanaan, S., Meilland, A., Stee, P., et al. 2008, A&A, 486, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kraus, S., Monnier, J. D., Che, X., et al. 2012, ApJ, 744, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Krtička, J., Owocki, S. P., & Meynet, G. 2011, A&A, 527, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levenhagen, R. S., & Leister, N. V. 2004, AJ, 127, 1176 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 1997, A&A, 321, 134 [NASA ADS] [Google Scholar]
 Maeder, A. 2002, A&A, 392, 575 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maeder, A., & Meynet, G. 2000, A&A, 361, 159 [NASA ADS] [Google Scholar]
 Maeder, A., & Meynet, G. 2001, A&A, 373, 555 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maeder, A., Grebel, E. K., & Mermilliod, J.C. 1999, A&A, 346, 459 [NASA ADS] [Google Scholar]
 Martayan, C., Baade, D., & Fabregat, J. 2010, A&A, 509, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McSwain, M. V., & Gies, D. R. 2005, ApJS, 161, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Meilland, A., Stee, P., Vannier, M., et al. 2007, A&A, 464, 59 [Google Scholar]
 Meilland, A., Millour, F., Kanaan, S., et al. 2012, A&A, 538, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meynet, G., Ekström, S., Maeder, A., & Barblan, F. 2007, in Active OBStars: Laboratories for Stellare and Circumstellar Physics, eds. A. T. Okazaki, S. P. Owocki, & S. Stefl, ASP Conf. Ser., 361, 325 [Google Scholar]
 Meynet, G., Ekström, S., Maeder, A., et al. 2013, in Seismology for studies of stellar rotation and convection, Lect. Notes Phys. (Berlin: Springer Verlag) [Google Scholar]
 Porter, J. M., & Rivinius, T. 2003, PASP, 115, 1153 [NASA ADS] [CrossRef] [Google Scholar]
 Reimers, D. 1975, Mem. Soc. Roy. Sci. Liège, 8, 369 [Google Scholar]
 Reimers, D. 1977, A&A, 57, 395 [NASA ADS] [Google Scholar]
 Rivinius, T., Baade, D., Štefl, S., et al. 2001, A&A, 369, 1058 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rivinius, T., Štefl, S., & Baade, D. 2006, A&A, 459, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Silaj, J., Jones, C. E., Tycner, C., Sigut, T. A. A., & Smith, A. D. 2010, ApJS, 187, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Stee, P. 2003, A&A, 403, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stee, P., & Meilland, A. 2009, in The Rotation of Sun and Stars, Lect. Notes Phys. (Berlin: Springer Verlag), 765, 195 [Google Scholar]
 Struve, O. 1931, ApJ, 73, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, A. R., Dougherty, S. M., Waters, L. B. F. M., & Bjorkman, K. S. 1990, A&A, 231, 453 [NASA ADS] [Google Scholar]
 Touhami, Y., Gies, D. R., & Schaefer, G. H. 2011, ApJ, 729, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Townsend, R. H. D., Owocki, S. P., & Howarth, I. D. 2004, MNRAS, 350, 189 [Google Scholar]
 Waters, L. B. F. M. 1986, A&A, 162, 121 [NASA ADS] [Google Scholar]
 Waters, L. B. F. M., Cote, J., & Lamers, H. J. G. L. M. 1987, A&A, 185, 206 [NASA ADS] [Google Scholar]
 Wisniewski, J. P., & Bjorkman, K. S. 2006, ApJ, 652, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Zahn, J.P. 1992, A&A, 265, 115 [NASA ADS] [Google Scholar]
 Zhao, M., Monnier, J. D., & Che, X. 2011, in IAU Symp. 272, eds. C. Neiner, G. Wade, G. Meynet, & G. Peters, 44 [Google Scholar]
 Zorec, J., Frémat, Y., & Cidale, L. 2005, A&A, 441, 235 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Timestep dependence of the total mass lost mechanically and of the mean mechanical massloss rates (in the 250 000 yr interval).
Properties of the 7 M_{⊙} stellar models at Z = 0.014 computed with various values of ω_{max}.
Mean mass fluxes during the criticalrotation phase for stars with Ω_{ini}/Ω_{crit} = 0.95.
Main characteristics of the decretion disks for models with Ω_{ini}/Ω_{crit} = 0.95.
Mean values of some disk’s characteristics for the 7 M_{⊙} model at solar metallicity, varying T_{0}/T_{eff}, p, and α.
All Figures
Fig. 1 Evolution of the Ω/Ω_{crit} ratio for 9 M_{⊙} at Z = 0.002, with initial values of Ω_{ini}/Ω_{crit} = 0.10, 0.30, 0.50, 0.70, 0.80, and 0.90 (bottom to top) during the MS phase. The solid lines show the present models, the dotted lines the models by Ekström et al. (2008). 

In the text 
Fig. 2 Evolution of the fraction of stars that rotate faster than a given limit in clusters of various ages at Z = 0.014, taking the Huang et al. (2010) distribution of velocities as the initial velocity distribution. The observational points represent the observed Be/(B+Be) ratio (McSwain & Gies 2005). The squared symbols (solid error bars) are for confirmed Be stars, and the triangle symbols (dashed error bars) are for all Betype stars. The ticks at the top of the plot show the typical lifetime of stars with various initial masses. 

In the text 
Fig. 3 Same as Fig. 2 for Z = 0.002. Observed Be/(B+Be) ratio extracted from Martayan et al. (2010), following the same procedure as in McSwain & Gies (2005). We consider the following age beams: 7 ≤ log (t) < 7.5, 7.5 ≤ log (t) < 8.0, 8.0 ≤ log (t) < 8.5, and log (t) ≥ 8.5. The error bars for the fraction of Betype stars are Poisson statistic error bars. 

In the text 
Fig. 4 Evolution of the mechanical massloss rate as a function of time for the 7 M_{⊙} stellar model at solar metallicity with a time step of about nine years (case n = 105, see text). Each point represents the mass loss rate in one time step. The quiescent phase between two massejection episodes lasts for about 40 years. 

In the text 
Fig. 5 Fraction of time spent with a given massloss rate for the 7 M_{⊙} models at Z = 0.014 with various time steps. 

In the text 
Fig. 6 HR diagram for 7 M_{⊙} models at solar metallicity computed with different values of ω_{max} (blue: 0.99, green: 0.95, and red: 0.90). As a comparison 5 and 9 M_{⊙} evolutionary tracks are also plotted. 

In the text 
Fig. 7 Evolution of the radius and diffusion time of the disk as a function of time for critically rotating models with an initial mass between 4 and 9 M_{⊙} at solar metallicity. 

In the text 
Fig. 8 Same is in Fig. 7 for an initial mass between 1.7 and 3 M_{⊙}. 

In the text 
Fig. 9 Evolution of the Be/(B+Be) ratio in clusters of various ages at Z = 0.014 obtained from our models assuming a mass (or T_{eff}) dependence for the Be phenomenon as proposed by Cranmer (2005, using (V_{eq}/V_{crit})_{min}. The observational points are the same as in Fig. 2. 

In the text 
Fig. 10 Mass loss – luminosity diagram. The colour points correspond to data by Waters et al. (1987), while the black crosses correspond to data from Stee (2003) and the black circle to the data from Carciofi et al. (2012) for 28 CMa. The thin dotted curve traces the upper limit for IR massloss rates obtained by Waters et al. (1987). The thick black curve corresponds to the mean mechanical massloss rates along the critical rotation phase, obtained as described in Sect. 4 and shown in Table 1, while the dashed black line corresponds to the mean value during the time in which the star is effectively losing mass. The grey region corresponds to the range of instantaneous massloss rates. 

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.