A two-dimensional perspective of the rotational evolution of rapidly rotating intermediate-mass stars Implications for the formation of single Be stars

Context. Recently, the ﬁrst successful attempt at computing stellar structure and evolution models in two dimensions was presented with models that include centrifugal deformation and self-consistently compute the velocity ﬁeld. Aims. The aim of the present study is to explore the rotational evolution of two-dimensional models of stars rotating at a signiﬁcant fraction of their critical angular velocity. From the predictions of these models, we aim to improve our understanding of the formation of single Be stars. Methods. Using the ESTER code, which solves the stellar structure of a rotating star in two dimensions with time evolution, we computed evolution tracks of stars of between 4 and 10 M (cid:12) for initial rotation rates ranging between 60 and 90% of the critical rotation rate. Furthermore, we compute models for both a Galactic metallicity and an SMC metallicity. Results. A minimum initial rotation rate at the start of the main sequence is needed to spin up the star to critical rotation within its main sequence lifetime. This threshold depends on the stellar mass, and increases with increasing mass. The models do not predict any stars above 8 M (cid:12) to reach (near-)critical rotation during the main sequence. Furthermore, we ﬁnd the minimum threshold of initial angular velocity is lower for SMC metallicity than for Galactic metallicity, which is in agreement with the increased fraction of observed Be stars in lower metallicity environments. The strong di ﬀ erence in the rotational evolution between di ﬀ erent masses is not predicted by any one-dimensional stellar evolution models. Conclusions. Our self-consistent two-dimensional stellar evolution models provide further insight into the rotational evolution of intermediate-mass stars, and our predictions are consistent with observations of velocity distributions and the fraction of Be stars amongst B-type stars. We ﬁnd that stars with a mass above 8 M (cid:12) do not increase their fraction of critical rotation during the main sequence. As a fraction of stars above 8 M (cid:12) have been observed to display the Be phenomenon, other processes or formation channels must be at play, or, alternatively, critical rotation is not required for the Be phenomenon above this mass.


Introduction
The rotational evolution of stars born with a convective core and radiative envelope is still a great unknown in the theory of stellar structure and evolution (SSE).Since a star's rotation profile (or shear profile rather) influences the efficiency of rotationallyinduced chemical mixing (e.g.Zahn 1992;Chaboyer & Zahn 1992), uncertainties in the predicted rotation profile in turn introduces uncertainties on the evolutionary pathway of the star in the Hertzsprung-Russell diagram (HRD).Studies using rotating stellar models are limited to one (spatial) dimension and the validity of such models for stars that are spinning at a significant fraction of their critical rotation rate is questionable.In Mombarg et al. (2023), we have successfully ran the first two-dimensional SSE models of a 12 M ⊙ star with the ESTER code (Espinosa Lara & Rieutord 2013;Rieutord et al. 2016).We have calibrated our models with a β Cephei pulsator that has been well-characterized by Burssens et al. (2023) and we have shown that all observationally derived quantities such as luminosity, effective temperature, core mass, and core-and surface rotation, could be reproduced by the ESTER models.
Rotation also plays an important role in the formation of classical Be stars (Struve 1931), which we define here as stars of spectral type B that are rotating sufficiently fast to form a decretion disk.The formation of Be stars is not fully understood, and formation channels have been proposed including binary interaction (e.g.Kriz & Harmanec 1975;Pols et al. 1991) or single star evolution (Bodenheimer 1995;Ekström et al. 2008;Granada et al. 2013;Hastings et al. 2020).It is generally believed that the formation of a disk requires the star to spin at large fraction of its critical angular velocity, but the exact limit is unknown and studies have argued that (near) critical rotation is not a necessity for the Be phenomenon based on measured subcritical rotation rates.While Townsend et al. (2004) argue that the observed subcritical rotation rates are underestimated because gravity darkening (von Zeipel 1924) has not been properly taken into account, later studies accounting for this effect still find subcritical rotation rates (e.g.Cranmer 2005;Frémat et al. 2005).Moreover, Espinosa Lara & Rieutord (2011) have shown that the von Zeipel law with a 1/4-exponent, as used by Townsend et al. (2004), overestimates the difference between the effective temperature at the pole and the equator, still reducing the impact of gravity darkening on Vsin i measurements.
The present paper follows up the work of Mombarg et al. (2023).In this previous work, we discovered that a 12 M ⊙ model never reaches a quasi-steady state, which would occur if nuclear evolution were extremely slow.It turns out that nuclear evolution of such a model is too fast, especially near the end of the main sequence, for the baroclinic waves to be damped.Indeed, as shown by Busse (1981), the initial distribution of angular momentum of a star relaxes to a steady state when all baroclinic waves are damped out.This damping occurs on a time scale1 t relax ∼ (R2 /K)(N 2 /Ω 2 ), namely on the thermal relaxation time R 2 /K (where R is the radius of the star and K the heat diffusivity) times the square ratio of the Brunt-Väisälä frequency N and the angular velocity Ω.If the nuclear evolution is much longer than this relaxation time then the evolution of the star can be computed as a series of steady states with the same total angular momentum (e.g.Gagnier et al. 2019b).In such a case, it has been shown that provided the initial angular velocity is fast enough, and mass-loss is negligible, the star inevitably reaches the critical rotation during the main sequence.The result of Mombarg et al. (2023) nuances this simple picture showing that nuclear evolution can be fast enough to prevent an evolution in a series of relaxed quasi-steady states.The interesting conclusion is that even if initial angular velocity is high and mass-loss unimportant, the star may not reach its critical angular velocity.This is a crucial point to investigate if we want to understand the Be phenomenon.Hence, we shall here extend the investigations of Mombarg et al. (2023) to the broad mass range of intermediate-mass stars, namely 4 to 10 M ⊙ , where nuclear evolution is slower and possibly longer than the relaxation time of baroclinic waves.The comparison of time scales is however not easy since parameters like K, N or Ω vary by orders of magnitude inside a star.It turned out that the computation of the full evolution of 2D-ESTER models was the most direct way to get a broad view of the rotational evolution of these early-type stars, including the effects of metallicity.
Hence, we organised the paper as follows: in Section 2, we briefly describe the physics of the ESTER models that we use.In Section 3, we present the predicted evolution of the rotation velocity at the stellar surface and compare these predictions with the measured projected surface velocities of B-type stars by Huang et al. (2010).We discuss the effect of metallicity on the rotational evolution in Section 4, and finally present our conclusions in Section 5.

ESTER models
The study of Gagnier et al. (2019b) presents the rotational evolution of intermediate-mass stars using 2D ESTER 2 models, where the chemical evolution is emulated with a succession of stationary models with the hydrogen mass fraction in the core given by, Here, m p is the proton mass, Q the energy released per reaction, and ⟨ϵ⟩ m the mass averaged nuclear energy production.The models of Gagnier et al. (2019b) do not account for the diffusion of material and thus the composition remains constant throughout  2021) reveals that, for most stars, SSE models with rotational mixing provide a better solution to the observed gravity-mode pulsation periods, than those with mixing by internal gravity waves.Moreover, a correlation was observed between the efficiency of chemical mixing and the rotation rate.
In ESTER, rotational mixing is based on the work of Zahn (1992) where the vertical diffusion coefficient scales with the shear, n • ∇Ω (projection perpendicular to the isobars), and the heat diffusivity, K, namely, (2) is the volume-averaged value of the Brunt-Väisälä frequency in the radiative envelope and η is a free constant that we set to a value such that η N 2 0 −1 V = 10 6 s 2 .This value prevents numerical issues due to steep chemical gradients building up during evolution.We fix the value of the horizontal chemical diffusion coefficient to 10 5 cm 2 s −1 .
Furthermore, in Espinosa Lara & Rieutord (2013) and Gagnier et al. (2019b) viscous effects are taken into account in the surface Ekman layers only, while in the new ESTER models which we presently use, viscous stresses are included throughout the star.This means we include a constant kinematic viscosity of ν = 10 7 cm 2 s −1 in both the horizontal and vertical direction, when we solve the angular momentum equation, Here, s is the distance to the rotation axis, v is the meridional velocity field, and ρ is the density.We note that this choice for the kinematic viscosity was adequate to reproduce the measured rotation profile of the β Cephei pulsator studied by Burssens et al. (2023) at the correct age, as demonstrated in Mombarg et al. (2023).As we limit ourselves to masses ≤ 10 M ⊙ , we assume mass loss can be neglected.
The ESTER models are discretized in 12 domains, which are spheroidal shells separated by isobars (Rieutord et al. 2016).Figure 1 shows the distribution of angular momentum per domain for a 10 M ⊙ star.From this figure, we can conclude that most of the star's angular momentum is located in the deep interior and therefore any angular momentum extracted from the stellar surface via winds is extremely small in the mass regime we study here.As we do not assume solid-body rotation, the angular momentum loss at the surface is less than predictions from 1D models like the ones used by Hastings et al. (2020).In Appendix A we show an example of a 2D rotation map at the zero-age main sequence (ZAMS), as well as for a model at the terminal-age main sequence (TAMS), and a model rotating at the critical rotation rate.The Keplerian critical angular velocity that we use in this paper as the point of criticality is defined by where R eq is the radius at the equator.Gagnier et al. (2019a) have shown that the reduction of Ω c due to the radiative acceleration at the stellar surface can be neglected in the mass regime we study here.We investigate the rotational evolution for both a metallicity typical for stars in the Milky Way (MW, Z MW = 0.02), and for a metallicity typical for stars in the Small Magellanic Cloud (SMC, Z SMC = 0.003).Lastly, we start our models at the ZAMS with a certain fraction of the critical angular velocity.These initial models are in a steady state like those of Espinosa Lara & Rieutord (2013).Throughout this paper, the mentioned ages concern the time spent on the MS.Pre-MS lifetimes range from 4.4 Myr for 4 M ⊙ to 0.34 Myr for 10 M ⊙ .

Evolution of rotation
We have computed evolution tracks for masses 4, 5, 6, 8, and 10 M ⊙ , all starting with an initial rotation rate of 75 per cent of the critical rotation frequency3 .
As we do not take into account any physics for the formation of a circumstellar disk, we stop the models once critical rotation is reached.Table 1 lists the corresponding initial velocities at the equator.

Model predictions
The evolution tracks in the HRD are shown in Fig. 2 for the effective temperature at the pole and at the equator.It should be noted that for stars close to the critical rotation limit, the difference between the effective temperature at the pole and equator can differ by as much as a factor two.The hooks seen at the end of some evolution tracks are the result of core contraction due to the steep drop in hydrogen in the convective core.Some models were stopped before the hook due to numerical problems.The hook seen in the equatorial-T eff track for the 6 M ⊙ model in the right panel of Fig. 2 is the result of the star reaching near critical rotation for a brief moment, but then decreasing the fraction of critical rotation again.For the Z SMC -case, we also show the 1D rotating evolution tracks from Georgy et al. (2013) (not available for Z MW = 0.02 as we assume in this paper).In general, the T eff predicted from the 1D rotating models is the average of the polar and equatorial T eff predicted by our ESTER models.The 1D models remain below Ω/Ω c = 0.8, while our 2D models predict stars below 8 M ⊙ to reach critical rotation during the main sequence.
Figure 3 shows the evolution of the fraction of critical rotation Ω eq /Ω c as a function of time spent on the main sequence and as a function of the hydrogen mass fraction in the core (X c ).We note several interesting findings.We find that for (Ω eq /Ω c ) i = 0.75, the stars of 4, 5 and 6 M ⊙ increase their fraction of critical rotation during the MS.While the 6 M ⊙ star spins down4 after about 45 Myr when it reaches a maximum of Ω eq /Ω c = 0.84, the 4 and 5 M ⊙ models continue to spin up and reach critical rotation around 77 Myr and 64 Myr, respectively.Our models for 8 and 10 M ⊙ predict a continuous decrease in the fraction of critical rotation during the MS.
Additionally, we have run models starting at 90 per cent of the initial critical rotation frequency.At this very high initial rotation rate, the 5 M ⊙ star reaches critical rotation after roughly 22 Myr, and the 6 M ⊙ star now also spins up to critical rotation only a few Myr after.In case of a 8 M ⊙ star, even an initial rotation rate as high as (Ω eq /Ω c ) i = 0.9 is not sufficient to reach critical rotation during the MS, and the star spins down and reaches Ω eq /Ω c = 0.7 near the end.Finally, we have also tested whether an initial rotation rate of 60 per cent of the critical rotation is still enough to spin up a 5 M ⊙ star to criticality.As can been seen in Fig. 3 (red dotted line), the star reaches a maximum of roughly Ω eq /Ω c = 0.7 and afterwards spins down.In summary, the threshold on the minimum initial rotation rate needed to reach critical rotation during the MS decreases with decreasing mass.

Uncertainties from chemical mixing
As mentioned before, the feedback of a chemical gradient on the efficiency of chemical diffusion is not taken into account in the ESTER models and Eq. ( 2) has a free parameter.In Mombarg et al. (2023), we have shown that this description of chemical mixing does reproduce the measured core mass by Burssens et al. (2023).The choice for the value of η will influence the exact age at which the star reaches critical rotation.More efficient mixing will result in a larger convective core and a slightly larger radiative envelope, thereby increasing the time scale on which baroclinic modes are damped and angular momentum is redistributed (Eq.(3) in Mombarg et al. 2023).We have recomputed the evolution track of the (5 M ⊙ , (Ω eq /Ω c ) i = 0.75, Z MW ) model, where we have changed η = 1.85 to 1.Such a change in the global scaling of the vertical diffusion coefficient increases the time to reach criticality to 67 Myr, that is, 3 Myr later (see Appendix B).Therefore, our conclusions will remain the same for changes to the value of η that are of order unity.

Comparision with observations
The work of Huang et al. (2010) presents the V sin i/V c distribution 5 of a large sample of cluster and field B-type stars according to the spectroscopic surface gravity at the rotational pole, log g pole (Be stars excluded).These authors find that stars between 2 and 4 M ⊙ on average have larger values for V sin i/V c than stars between 7 and 13 M ⊙ .In Fig. 4, we show for each mass V eq /V c averaged over the MS life time (or the time until critical rotation is reached) computed from the 2D ESTER models for Z MW .The points are normalised to the initial value at the ZAMS, (V eq /V c ) i .We observe a decrease in V eq /V c /(V eq /V c ) i with respect to the stellar mass that has a similar slope compared to upper limits of V eq /V c derived by Huang et al. (2010) for nonemission B-type stars (grey points in Fig. 4).In other words, observations show that the upper limit on V eq /V c decreases with mass, while our models show the same trend once V eq /V c is rescaled by its initial value at the ZAMS.This latter rescaling could be avoided by running a large set of models representing the evolution of a star cluster and computing the maximum angular velocity per bin of mass, but such a simulation is presently too demanding in computing time.

Effect of metallicity
We have also computed evolution tracks for the same masses with a metallicity typical of the SMC.The results are shown in Fig. 5.We find that at lower metallicities, the stars spin up faster, and the 5 M ⊙ star reaches critical rotation around 40 Myr, whereas for Galactic metallicity, this would take 63 Myr.Similarly, the 5 M ⊙ star starting with (Ω eq /Ω c ) i = 0.9 reaches critical rotation after 13 Myr at SMC metallicity, whereas it takes 22 Myr at MW metallicity.Moreover, the lower metallicity is sufficient to also spin up the 6 M ⊙ star to (near) critical rotation around 45 Myr.Since the time scale on which baroclinic modes are damped scales with the square of the thickness of the radiative envelope (Busse 1981), we expect angular momentum to be faster redistributed in lower-metallicity stars as they are 5 The critical rotation velocity is estimated from the Roche model.more compact.On the other hand, the nuclear time scale will also decrease with decreasing metallicity since a lower metallicity results in higher central temperatures.Nonetheless, we find that for the 5 M ⊙ ZAMS model, the baroclinic time scale differs with 150 per cent between Z MW and Z SMC , while the nuclear time scale only differs by 25 per cent.Likewise, for the 10 M ⊙ ZAMS model, the difference in the baroclinic time scale is 50 per cent, compared to 14 per cent for the nuclear time scale.Again, our models predict no stars ≥ 8 M ⊙ that significantly increase the fraction of critical rotation compared to their initial value when (Ω eq /Ω c ) i equals 0.75.Also, even at SMC metallicity the spin up of the 5 M ⊙ star starting with (Ω eq /Ω c ) i = 0.6 stalls after 60 Myr having reached a maximum of Ω eq /Ω c = 0.7.
It is interesting to note that our 2D models suggest an increase in the fraction of Be stars with decreasing metallicity.Indeed, larger fractions of Be stars have been observed in lower metallicity environments (Martayan et al. 2010;Iqbal & Keller 2013), but this is not reproduced by the 1D models by Ekström et al. (2008Ekström et al. ( , 2012) ) and Granada et al. (2013).The predictions of Hastings et al. (2020, also based on 1D models) do predict this observed relation between the fraction of Be stars and the metallicity, but in their models, this is due the dependence of the mass-loss rate on the metallicity.Thus, this should only be of importance in more massive stars, namely above 10 M ⊙ , while we find a strong metallicity effect also at lower masses and we do not account for mass loss.

Conclusions
Following the successful calibration of the first 2D stellar structure and evolution models by Mombarg et al. (2023), we have presented here rotating stellar evolution tracks for intermediatemass stars computed with the 2D ESTER code.Solving the stellar structure in two dimensions instead of one has the added value that we account for the centrifugal deformation of the star and we compute the rotation profile self-consistently from the baroclinic torque (e.g.Eq. ( 26) in Espinosa Lara & Rieutord 2013).We have presented evolution models for stellar masses 4, 5, 6, 8 and 10 M ⊙ for initial rotation rates between 60 and 90 per cent of the critical angular velocity.Summarising, the most important findings from our 2D models are the following.
-An initial rotation rate of (Ω eq /Ω c ) i = 0.75 is not enough for stars with mass ≥ 6 M ⊙ and galactic metallicity to evolve towards critical rotation.Stars of 5 M ⊙ will reach critical rotation around 60 Myrs.
-At an initial rotation rate of (Ω eq /Ω c ) i = 0.9, stars of 5 and 6 M ⊙ do reach critical rotation around 20-25 Myrs.-Stars in the SMC, thus with a lower metallicity than Galactic stars, spin up faster and stars with masses 5 and 6 M ⊙ can reach critical rotation around 40 Myr.-Stars with mass ≥ 8 M ⊙ always move away from critical rotation during the main sequence, even for SMC metallicity or (Ω eq /Ω c ) i = 0.9.
If we assume that the emergence of the Be phenomenon requires the star to rotate close to criticality, we can make the following predictions for the fraction of single Be stars with respect to the total number of B-type stars.First, if we hypothesise that single Be stars are born as rapid rotators, then our models show Fig. 4. The average value of the equatorial velocity with respect to the critical one, V eq /V c , divided by the initial value at the ZAMS, as a function of stellar mass.Data points in grey show the observationally derived upper limits of V eq /V c , taken from the last column of Table 7 of Huang et al. (2010).that the minimum initial rotation rate needed to spin up stars to criticality decreases with decreasing metallicity.Second, our models predict that no stars with mass above 8 M ⊙ can stay rotating at or near critical rotation.We thus conclude that for massive B-type stars ≳ 8 M ⊙ the velocity threshold to expel matter from the surface must be lower than for the lower-mass ones.
In order to constrain the prevalence of the single-star formation channel compared to the binary formation channel, the upper limit on the fraction of critical rotation at which material can be expelled from the surface (for a given stellar mass) needs to be quantified.Additionally, while the predicted rotation profile from ESTER are in agreement with the observed one in case of the β Cephei pulsator HD192575 (Mombarg et al. 2023), future measurements of rotation profiles of massive stars with a precise mass and age will help us get a better picture of the accuracy of the current treatment of angular momentum transport.Even though we have limited ourselves to masses up to 10 M ⊙ , we expect the more massive B-type stars that weigh up to 16 M ⊙ to also decrease the fraction of critical rotation during the MS.Assuming a mass loss rate predicted by Björklund et al. (2021), a 16 M ⊙ has lost about 0.02 M ⊙ by the point of corehydrogen exhaustion (or about a factor 10 more when assuming Vink et al. (2001) mass loss rates).Hence, the decrease in the value of Ω c due to mass loss is negligible.Moreover, Gagnier et al. (2019a) have demonstrated that the reduction of the critical rotation velocity due to radiative acceleration is small for stars below 40 M ⊙ .However, the decrease of Ω eq /Ω c with respect to time can be accelerated as a result of angular momentum loss at the stellar surface (Gagnier et al. 2019b).While we expect this effect to be small, 2D evolution with mass loss need to be further studied to quantify this.
Observed fractions of Be stars compared to the total number of B-type stars in clusters show little increase over time (Mc-Swain & Gies 2005).In light of our results for single stars, this would imply that the majority of stars is born with a rotation rate ≲60 per cent of the critical angular velocity at the ZAMS.The higher fraction of Be stars in low-metallicity environments observed around an age of 40 Myr (Martayan et al. 2010) is in line with our predictions as we predict a lower threshold on the initial rotation for stars to reach critical rotation and a smaller time span to do so.Furthermore, the predicted evolution of the surface velocity by the ESTER models seems to be in good agreement with the measurements of the fastest rotating B-type stars in the sample of Huang et al. (2010).
In summary, the results of two-dimensional stellar structure and evolution models that include differential rotation and meridional currents show the limitations of one-dimensional models to predict the evolution of rapidly rotating stars.In addition to the (projected) velocity at the surface, the surface abundance of the nitrogen-14 isotope is also a useful observable to test the theory of rotational evolution in massive stars.In stars that are undergoing hydrogen burning via the CNO-cycle, an overabundance of nitrogen is created in the core and the enhancement of its abundance at the surface is indicative of the efficiency of chemical (rotational) mixing (e.g.Hunter et al. 2009;Brott et al. 2011).In a future paper, we will investigate the evolution of surface abundances predicted from rotational mixing with 2D evolution models.

Fig. 1 .
Fig.1.Angular momentum integrated per domain as a function of the radial coordinate in the equatorial direction.Once for a steady-state model (no viscous stresses), and once for the ESTER models used in this work.Both models are computed at X core /X ini = 0.12 and have a rotation rate Ω eq /Ω c = 0.57.

Fig. 2 .
Fig. 2. Hertzsprung-Russell diagram showing the evolution of the ESTER models, colour-coded by the fraction of critical rotation.For each mass, a track is plotted once using the effective temperature at the pole and once at using the effective temperature the equator.The left (right) colour bar corresponds to the pole (equator).For comparison, non-rotating MESA (r23.05.1;Paxton et al. (2011, 2013, 2015, 2018, 2019); Jermyn et al. (2023)) tracks are shown with the grey lines.The MESA tracks were computed up to the base of the red giant branch, whereas the ESTER tracks are stopped close to the TAMS or when critical rotation is reached.In the right panel, also rotating 1D models (same initial rotation) from Georgy et al. (2013) are shown in blue.

Fig. 3 .
Fig. 3. Evolution of the fraction of critical rotation as a function of time spent on the MS (top panel) and as a function of the mass fraction of hydrogen in the core with respect to the initial mass fraction (middle panel).The bottom panel shows the velocity at the equator as a function of the surface gravity at the pole.The different line styles correspond to the different initial rotation rates at the start of the MS.The models shown here are for Z MW .

Table 1 .
The equatorial velocity at the ZAMS.The last column indicates the time spent on the MS before reaching critical rotation (no value means critical rotation is not reached).