Astronomy & Astrophysics manuscript no.

Be_paper c ESO 2019

December 12, 2019

The Single Star Path to Be Stars

B. Hastings, Chen Wang, and N. Langer

Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
e-mail: bhastings@[Link]

December 12, 2019

ABSTRACT
arXiv:1912.05290v1 [[Link]] 11 Dec 2019

Context. Be stars are rapidly rotating B main sequence stars, which show line emission due to an outflowing disc. By studying the
evolution of rotating single star models, we can assess their contribution to the observed Be star populations.
Aims. We identify the main effects which are responsible for single stars to approach critical rotation as functions of initial mass and
metallicity, and predict the properties of populations of rotating single stars.
Methods. We perform population synthesis with single star models of initial masses ranging between 3 and 30 M , initial equatorial
rotation velocities between 0 and 600 km s−1 at compositions representing the Milky Way, Large and Small Magellanic Clouds. These
models include efficient core-envelope coupling mediated by internal magnetic fields and correspond to the maximum efficiency of
Be star production. We predict Be star fractions and the positions of fast rotating stars in the colour-magnitude diagram.
Results. We identify stellar wind mass-loss and the convective core mass fraction as the key parameters which determine the time
dependance of the stellar rotation rates. Using empirical distributions of initial rotational velocities, our single star models can repro-
duce the trends observed in Be star fractions with mass and metallicity. However, they fail to produce a significant number of stars
rotating very close to critical. We also find that rapidly rotating Be stars in the Magellanic Clouds should have significant surface
nitrogen enrichments, which may be in conflict with abundance determinations of Be stars.
Conclusions. Single star evolution may explain the high number of Be stars if 70 to 80% of critical rotation would be sufficient to
produce the Be phenomenon. However even in this case, the unexplained presence of many Be stars far below the cluster turn-off
indicates the importance of the binary channel for Be star production.
Key words. stars: massive – stars: Be – stars: rotation – stars: evolution

1. Introduction lar to coronal mass ejections as seen in the Sun, as suggested by

Balona & Ozuyar (2019) based on recent TESS results.
Ever since their discovery over 150 years ago (Secchi 1866), On the other hand, Townsend et al. (2004) have argued that
Be stars have offered a promising, although misted window into all Be stars in fact rotate very close (3rot /3crit > 0.95) to the criti-
massive star evolution and structure. It was proposed by Struve cal velocity, with those which have low measured rotation rates
(1931) that Be stars are fast rotators, whose emission lines origi- being strongly affected by gravity darkening. Following the Von
nate from a circumstellar decretion disc, a picture which is main- Zeipel law (von Zeipel 1924), gravity darkening in a fast rotat-
tained until today (Rivinius et al. 2013). Yet, it is still not clear ing star makes the stellar pole, which has a low rotational veloc-
how fast a B-type star must rotate in order to become a Be star. ity, more luminous than the equator which has a high rotational
velocity, resulting in the star appearing as though it is rotating
For a decretion disc to form, the equatorial rotation velocity slower than in reality.
3rot is expected to be a significant fraction of the critical rotation A further question surrounding Be stars is why the phe-
velocity, vcrit , defined as the rotation velocity at which material nomenon seems to be restricted mostly to B-type stars and why
at the equator becomes unbound from the star. Observational ev- Be stars are more common in certain spectral classes than oth-
idence suggests that the threshold rotation rate for the Be phe- ers. Observations in the Milky Way show that the fraction of Be
nomenon is mass dependant, and could be as low as 3rot /3crit = stars in a certain spectral class varies across spectral type with
0.6 for stars more massive than 8.6 M and as high as 3rot /3crit = the most Be-stars found at B1-B2 classification, where the Be
0.96 for stars with M < 4 M (Huang et al. 2010). Similarly fraction is 34% while in comparison the Be fraction for B9 stars
Zorec et al. (2016) find that the Be phenomenon is character- is 8% and the total fraction of Be stars to B stars was measured to
ized by a wide range of true velocity ratios (0.3 < 3rot /3crit < be 17% (Zorec & Briot 1997). Furthermore, Oe stars seem to be
0.95) and that the probability that Be stars are critical rotators rather rare, with less than 20 having been detected in the Milky
is small. In this case one must look for an additional mechanism Way (Li et al. 2018). It is not clear whether this is caused by pro-
to feed the Be disc. Pulsations seem promising as they can serve cesses within O stars themselves, the mechanisms responsible
to kick matter from the surface of a star, however it is found for forming a disc around a fast rotating star, or the conditions
that not all Be stars pulsate (Baade et al. 2002) and that among under which very massive stars form.
those that do there is a wide range of pulation frequencies and The existence of around 150 (Raguzova & Popov 2005) de-
types (Rivinius et al. 2013). Another possibility is that the disc tected binary systems consisting of a Be star and a compact ob-
is fed through outbursts of magnetically active starspots, simi- ject (so called Be/X-ray binaries) demonstrates that binary inter-
Article number, page 1 of 14
 A&A proofs: manuscript no. Be_paper

actions can spin up a star significantly (Kriz & Harmanec 1975;

Pols et al. 1991; Liu et al. 2006; Langer 2012). In order to build !0.43
a full model of the Be phenomenon, binary evolution must be as 1
Ṁ(Ω) = Ṁ(0) (1)
well understood as the channels for forming single Be stars. 1 − 33crit
rot

van Bever & Vanbeveren (1997) used binary star evolution

calculations to predict that at most 20% of the Be star popula- where
tion results from binary interactions. Despite this, observations r
of surface nitrogen abundances in Magellanic Cloud Be stars are GM κL
3crit = (1 − Γ) ; Γ = . (2)
in disagreement with fast rotating single star models (Dunstall R 4πcGM
et al. 2011). Thus questions such as, which is the dominant Be
star formation channel and what are the differences between the For a detailed description of the models see Brott et al.
single and binary Be population remain open. (2011).
Considerable efforts have been made to predict the relative
fractions of Be stars through single star modelling, most notably 2.2. Population Synthesis
by Ekström et al. (2008). These models include coupling be-
To predict properties of populations of rotating stars we use pop-
tween core and envelope by hydrodynamic viscosities through
ulation synthesis to model open star clusters (i.e. collections of
the shear instability. However, a stronger coupling may be pro-
coeval stars without any continuous star formation) at various
duced by internal magnetic fields known as the Tayler-Spruit dy-
ages. For a cluster age t, we select pairs of random values from
namo (Spruit 2002) which produces very efficient angular mo-
a Salpeter initial mass distribution (with exponent 2.35) and a
mentum transport throughout the star so that near solid-body
distribution of initial critical velocity fraction, Mi and 3rot /3crit i .
rotation occurs. The models analysed here include such effects,
Then we find the masses, M1 , M2 on the model grid that are
and therefore employ the maximum efficiency of spinning up the
straddling the chosen mass value, such that M1 < Mi < M2 . For
surface layers due to the core contraction during hydrogen burn-
M1 and M2 we interpolate the hydrogen burning lifetimes, t MS
ing, which is a key factor to produce single Be stars. At present
as a function of initial critical velocity fraction to obtain the hy-
the Tayler-Spruit dynamo is used to explain relatively slow rota-
drogen burning lifetimes at the chosen value, 3rot /3crit i . Next the
tion rates in white dwarfs (Suijs et al. 2008) and young pulsars
hydrogen burning lifetime, t MS ,i , of a model with mass Mi and
(Heger et al. 2005).
initial rotation 3rot /3crit i is found by interpolating between the
In Section 2 the models are introduced and our approaches hydrogen burning lifetimes of M1 and M2 . The fractional life-
are outlined. In Section 3, models from the grid are analysed un- time is then given by t/t MS ,i . If the fractional lifetime is greater
der conditions of fast and slow rotation for MW and SMC metal- than 1, the star will not be hydrogen burning anymore so the
licities, the factors governing the approach to the critical velocity process is abandoned and new samples are drawn. We then se-
are investigated and the expected surface nitrogen abundances of lect models with masses M1 and M2 at fractional hydrogen burn-
Be stars is investigated. In Section 4, population synthesis is per- ing times t/t MS ,i . An interpolation of the quantity of interest, Q
formed to calculate the expected fractions of single Be stars in across initial critical velocity fraction gives the values of Q for
clusters of differing ages and metallicities and predict the posi- masses M1 , M2 with initial rotation 3rot /3crit i and fractional hy-
tions of fast rotating stars in the colour-magnitude diagram. drogen burning time t/t MS ,i . One final interpolation between M1
and M2 gives the predicted quantity of the selected mass Mi at
the given cluster age. The quantities of interest are luminosity,
critical velocity fraction at the current time and effective temper-
2. Method ature.
The initial rotational velocity distribution used was taken
2.1. Stellar Models from VLT-FLAMES observations of early B stars in the 30 Do-
radus region of the LMC (Dufton et al. 2013) and is shown in
We analyse the single star evolutionary models of Brott et al. Fig. 1. The deconvolved distribution of equatorial rotational ve-
(2011) to predict the properties of rotating single stars through- locities was converted to a distribution in critical velocity frac-
out their main-sequence evolution. We consider masses from tion by applying a mapping between the two as determined from
3 M to 30 M at various initial rotation velocities ranging from the 15 M LMC models at ZAMS and then normalizing such that
zero to approximately 600 km s−1 . Because the model grid is the integral over the whole probability density function equals
spaced in initial rotational velocities, but the initial critical ro- unity. It is noted that for the heaviest masses on the grid, the dis-
tation velocity increases with mass, our grid does not contain tribution extends beyond its limits. When such a massive, very
models with initial values of 3rot /3crit greater than around 0.7 for fast rotating star is chosen from the distributions, instead the
initial masses greater than 25 M . fastest rotator in the grid is used. Observations in 30 Doradus
Three initial chemical compositions represent metallicities show that no O stars are observed to rotate with deconvolved
of the Milky Way (MW), Large Magellanic Could (LMC) and equatorial velocities much greater than 500 km s−1 (Ramírez-
Small Magellanic Cloud (SMC). The models include internal Agudelo et al. 2013; Dufton et al. 2011). For a 25 M star to
transport of angular momentum via the Taylor-Spruit dynamo rotate at a critical velocity fraction of 0.7, it would require an
(Spruit 2002) which has the effect of enforcing near solid-body equatorial rotation velocity of the order 700 kms−1 . It is thus
rotation throughout most of the main-sequence evolution. The safe to assume that O stars do not enter the ZAMS with initial
adopted mass-loss scheme is given by Vink et al. (2000). An en- critical velocity fractions much greater than 0.65, or that if they
hancement of the mass loss due to rotation is used as outlined do, they spin down very quickly.
in Yoon & Langer (2005), whereby the mass loss rates are in- We note that adopting the observed distribution of rotational
creased by a factor depending on the ratio of the rotation velocity velocities of Dufton et al. (2013) as the initial velocity distri-
to the critical velocity like bution for stars in our synthetic populations may introduce an
Article number, page 2 of 14
 Hastings, Wang, and Langer: The Single Star Path to Be Stars

vrot /kms −1 maintain hydrostatic and thermal equilibrium. Eq. 2 shows that
0 100 200 300 400 500 600 as the stellar radius increases, the critical velocity decreases, thus
2.5 during main-sequence evolution, the critical velocity will fall.
In the absence of internal angular momentum transport, as
core density increases, the local conservation of angular momen-
2.0 tum will demand that the angular velocity of the core increases.
Likewise as the envelope expands, the angular velocity of the
Prob. density

1.5 envelope will decrease. This results in an angular velocity gradi-

ent developing between the convective core and radiative enve-
lope. However when the core and envelope are coupled via angu-
1.0 lar momentum transport, angular momentum is transported from
the core to the envelope, decreasing the angular velocity gradi-
ent throughout the star. The physical processes responsible for
0.5 the angular momentum transport in the models studied here are
magnetic torques arising from the Tayler-Spruit dynamo (Spruit
0.0 2002), which leads to near solid body rotation. While the en-
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 velope is expanding and the star is rotating as a solid body, the
vrot /vcrit critical rotation velocity will decrease while the equatorial rota-
tion velocity drops only slowly or even increases (see Section
Fig. 1. The deconvolved rotation distribution of early B stars as ob- 3.3).
served by VLT-FLAMES Survey (Dufton et al., 2013). The distribution Fig. 2 gives examples of the evolution of critical velocity,
was converted to 3rot /3crit using 15 M LMC models at ZAMS. The up- equatorial velocity and fraction of critical velocity during main-
per scale shows how the critical velocity fraction values, 3rot /3crit match sequence evolution for MW and SMC models of initial masses
to the equatorial velocities, 3rot . 5,15 and 25 M . All models in the plot have an initial critical ve-
locity fraction of approximately 0.6. Although the less massive
inconsistency, since the the sample of Dufton et al. (2013) con- models have slower equatorial velocities, they also have lower
sists of field stars of all ages. However, as discussed in Dufton critical velocities (because of a relatively weak dependance on
et al. (2013) (see also Sect. 3.1 below), the rotational velocities radius with mass), making the critical velocity fraction nearly
of the considered single stars are expected to change very little the same for all models in the plot. It is seen that for the 5 M
during their main sequence evolution. If binary evolution affects models the equatorial rotation velocity remains nearly constant,
this distribution de Mink et al. (2013), then we would overesti- while the critical velocity decreases. For the 15 M models, the
mate the number of stars which are born rotating very rapidly. In equatorial rotation velocity increases due to the effects of angular
this case, the predicted number of Be stars from our models may momentum transport. As a result, the critical velocity fraction is
be considered as upper limits. generally increasing during hydrogen burning such that the stars
To compare the models with observations of the SMC open are evolving closer to critical rotation, unless angular momentum
cluster NCG 330 (Milone et al. 2018) in the colour-magnitude is drained at a high rate due to mass loss (see Section 3.2).
diagram, the bolometric luminosities and effective temperature
are converted to absolute magnitudes in the Hubble Space Tele- 3.2. The Effect of Mass-loss
scope Wide-Field Camera 3 filters F814W and F336W by inter-
polating tables based on synthetic stellar spectra (Girardi et al. The dominating effect of mass-loss through stellar winds is to
2002). The values of distance modulus and reddening adopted remove angular momentum, not mass. This is so because even
are (m − M) = 18.92 and E(B − V) = 0.06 respectively. The in the absence of magnetic fields, the rate of angular momen-
absorption coefficients used are AF814W = 2.04E(B − V) and tum loss relative to the total angular momentum is about 10
AF336W = 5.16E(B − V) (Milone et al. 2009). The effects of times larger than the rate of mass loss relative to the stellar mass
gravity darkening are included as described by Espinosa Lara (Langer 1998). Thus a star’s mass-loss may strongly affect the
& Rieutord (2011), whereby the effective temperature and lumi- spin evolution.
nosity of a star are multiplied by parameters that depend on the Panels a, c, e of Fig. 2 show that models at MW metallicity
inclination angle and fraction of angular critical velocity. Then experience a turn-over in the evolution of their critical velocity
using these corrected effective temperature and luminosity val- fraction. This is because in the late stages of hydrogen burn-
ues, we calculate the absolute magnitudes as described above. ing, the star’s mass-loss rate increases significantly which has
The inclination angles, i , in our synthetic population are chosen the effect of removing angular momentum from the surface at a
such that cos(i) is uniformly distributed between 0 and 1, mean- rate which cannot be compensated by internal angular momen-
ing that it is more likely to observe any given star equator-on tum transport mechanisms, meaning that solid body rotation is
than pole-on. Such a distribution describes a random orientation no longer a good approximation and the equatorial velocity de-
of the rotation axis. creases. This period of strong mass-loss is caused by the iron
opacity bistability in which partial recombination of Fe ions at
effective temperatures of around 22kK causes a sharp increase
3. Results in opacity and hence mass-loss (Lamers et al. 1995). This be-
3.1. Spin Evolution haviour is strongly dependant on metallicity and so weaker in
the SMC or LMC models. Comparing Panels c and d of Fig. 2, it
During the evolution of a slowly rotating star during core hydro- can be seen that despite both MW and SMC models starting with
gen burning, a strong chemical gradient develops between the approximately equal critical velocity fractions, the SMC model
convective core and the radiative envelope. The core density in- achieves a much larger critical velocity fraction at the end of the
creases, and as a reaction the envelope must expand in order to main-sequence. This comparison between MW and SMC mod-
Article number, page 3 of 14
 A&A proofs: manuscript no. Be_paper

900 900 900

1.0 1.0 1.0
800 a 800 c 800 e
700 0.8 700 0.8 700 0.8
600 600 600
v/kms −1

v/kms −1

v/kms −1
0.6 0.6 0.6

vrot /vcrit

vrot /vcrit

vrot /vcrit
500 500 500
400 0.4 400 0.4 400 0.4
300 300 300
5 M¯ 0.2 15 M¯ 0.2 25 M¯ 0.2
200 MW 200 MW 200 MW
100 0.0 100 0.0 100 0.0
900 900 900
1.0 1.0 1.0
800 b 800 d 800 f
700 0.8 700 0.8 700 0.8
600 600 600
v/kms −1

v/kms −1

v/kms −1
0.6 0.6 0.6
vrot /vcrit

vrot /vcrit

vrot /vcrit
500 500 500
400 0.4 400 0.4 400 0.4
vrot
300 300 300 vcrit
5 M¯ 0.2 15 M¯ 0.2 25 M¯ 0.2
200 SMC 200 SMC 200 SMC vrot /vcrit
100 0.0 0.2 0.4 0.6 0.8 1.0 0.0 100 0.0 0.2 0.4 0.6 0.8 1.0 0.0 100 0.0 0.2 0.4 0.6 0.8 1.0 0.0
t/tMS t/tMS t/tMS
Fig. 2. The evolution of equatorial rotational velocity 3rot (thick dot-dashed), critical velocity 3crit (dashed) and the critical velocity fraction
3rot /3crit (solid red) for 5 M (panels a, b), 15 M (panels c, d) and 25 M (panels e, f) models at MW metallicity (panels a, c, e) and SMC
metallicity (panels b, d, f). The initial critical velocity fraction values are all approximately 0.6. The X-axis indicates the fractional main sequence
lifetime, t/t MS .

els illustrates the effect of mass-loss on the approach to the crit- Another effect is rotationally enhanced mass loss. As a star
ical velocity. approaches the critical rotation velocity, material at the equator
becomes less tightly bound due to the centrifugal force, thus one
Fig. 3 shows the ratio of angular momentum at the end of
would expect angular momentum lost through winds to increase
hydrogen burning to the initial angular momentum for models
with rotation velocity. In the models this is governed by Eq. 1.
of varying mass and initial critical velocity fraction for MW and
With 3rot /3crit = 0.8, the mass loss rates are doubled, so rotation-
SMC models. It is clear that almost every SMC model loses
ally enhanced mass loss plays only a large role when very high
much less angular momentum than the corresponding model at
critical rotation fractions are achieved. As an example the 5 M
MW metallicity. The exceptions are fast rotating massive SMC
models at MW metallicity in Fig. 3 show that in the range of
models which undergo quasi-chemically homogeneous evolu-
initial 3rot /3crit from 0 to 0.4, the total angular momentum lost is
tion, and during so become very luminous which leads to in-
almost the same. However when approaching critical rotation,
creased mass loss rates. For SMC models one may judge that
the enhancement factor becomes divergent, so models rotating
mass-loss becomes irrelevant to the angular-momentum budget
near the critical velocity experience tremendous mass loss.
below masses of around 10 M , where most models (except very
fast initial rotators, say with 3rot /3crit > 0.7) retain more than
A further effect of rotationally induced mixing is to increase
90% of their angular momentum. For MW models we see that
the overall mean molecular weight, µ in the star, compared to
only slowly rotating models less massive than 5 M retain more
models with no rotationally induced mixing. Homologous mod-
than 90% of their angular momentum. As expected, the effect of
els suggest a strong dependance on luminosity with mean molec-
mass-loss is strongly metallicity dependant. It is also seen that
ular weight such that L ∝ µ4 (Kippenhahn & Weigert 1990). In
below 15. . . 20 M , angular momentum loss becoming less de-
turn, mass-loss rates are dependant on the luminosity, for the
pendant on mass.
wind prescription used in the models the dependance is approx-
Fig. 3 shows that for any given mass the fraction of angular imately Ṁ ∝ L2 (Vink et al. 2000). Thus rotationally induced
momentum lost is a strong function of the initial rotation. For ex- mixing leads to higher mass loss and angular momentum loss.
ample a 5 M model at MW metallicity can lose between 2 and For models which experience quasi-chemical homogeneous evo-
30 % of its initial angular momentum. One contributing factor lution, where the star can become a helium star, this effect be-
is the lifetime effect, whereby under the effects of rotation, rota- comes very apparent. Models which undergo quasi-chemical ho-
tional mixing causes hydrogen to be mixed into the centre of the mogeneous evolution for the duration of hydrogen burning (de-
star and so hydrogen burning can continue for a longer time. For fined by a monotonically increasing surface helium mass frac-
all of our models, the hydrogen burning lifetime enhancement tion) are marked with black circles in Fig. 3. Similarly models
between the non-rotating model and that with an initial rotation which undergo a phase of quasi-chemical homogeneous evolu-
velocity of 600 km s−1 never exceeds a factor of 1.5. Thus for tion (defined by having a slowly increasing difference between
models losing only small fractions of their total angular momen- surface and central helium mass fraction for longer than one third
tum as slow rotators (such as low mass MW metallicity models of the hydrogen burning lifetime) are marked by a red circle.
and all SMC metallicity models), the lifetime effect cannot ex- Although these models do not have high initial critical rotation
plain wholly the increase in angular momentum loss. fractions, they still lose large fractions of their angular momen-
Article number, page 4 of 14
 Hastings, Wang, and Langer: The Single Star Path to Be Stars

1.0 gular momentum budget. Thus to maintain a constant rotational

velocity during envelope expansion, a relatively low angular mo-
mentum transport rate is required. On the other hand, for an
0.8 envelope dominated star to rotate at a constant velocity while
the envelope is expanding, the angular momentum transport rate
L final /L initial

from core to envelope must be high. It is then likely that inter-

0.6 nal angular momentum transport mechanisms are unable to meet
this demand, and as a result the rotational velocity of the stellar
surface will decrease due to the effect of local angular momen-
0.4 1.0 tum conservation.
In our models, there is an inner region of the star from which
MW angular momentum is being transported and there is an outer re-
0.2 0.8 gion which the angular momentum is being transported to. Inbe-
5 10 15 20 25 30 tween these regions there must be a point which neither gains nor
M /M ¯

initial vrot /vcrit

loses any specific angular momentum. The location of this angu-
0.6 lar momentum "valve" will give an indication as to the strength
1.0 of the core mass effect as discussed above. Figs. 4. a, b show the
specific angular momentum profiles of a 5 and 15 M model at
0.4 one time early in their evolution and one time near the end of
0.8 hydrogen burning. Hydrogen mass fraction profiles are plotted
to show how advanced the evolution is. The models shown are
L final /L initial

the same as in Figs. 2b,d and have equal initial critical velocity
0.6 0.2 fractions of 0.6. It can be seen that there is a point for each mass
at which the specific angular momentum does not change.
Figure 4.c shows this more clearly, where the difference in
0.4 0.0 angular momentum at both times, divided by the angular mo-
mentum at the earlier time, ∆ j/ j1 , is plotted. Here, regions where
SMC angular momentum is gained have a positive value, whereas re-
0.2 gions were angular momentum is lost have a negative value. For
5 10 15 20 25 30 the more massive model, the point with a constant specific an-
M /M ¯ gular momentum is closer to the star’s edge than for the less
massive model. For the 15 M model, approximately 90% of the
Fig. 3. The fraction of angular momentum at the end of hydrogen burn- total mass is acting as a donor of angular momentum, while for
ing L f inal , to angular momentum at the start of hydrogen burning Linitial , the 5 M mode the figure is 80%. Using the arguments above,
as a function of mass for MW (top) and SMC (bottom) metallicities. therefore the 15 M model will approach the critical rotation ve-
The colour of the points represents the initial critical velocity fraction
locity more easily. Furthermore by inspecting the area under the
3rot /3crit . Models with a growing helium surface abundance throughout
the duration of their evolution are marked by a black circle. Models curves, in the region where ∆ j/ j1 is positive in Fig. 4.c one can
with a growing helium surface abundance for part of their evolution are determine how much relative angular momentum is gained. For
marked by a red circle. example, if a region from M m
= mi to M m
= 1 had doubled its
R1
total angular momentum, the integral m ∆ j/ j1 d( M m
) would be
i
equal to (2 − 1)(1 − mi ). Fig. 4.c shows that the relative angu-
tum. Quasi-homogeneous evolution occurs more readily in the
lar momentum gain of the matter in the envelope of the 5 M
lower metallicity models because mass loss being a strong func-
model is greater than that of the 15 M model. This tells us that
tion of metallicity, so the MW metallicity models slow down
to maintain near solid body rotation, a relatively smaller amount
relatively quickly, rotational mixing becomes less effective and
of angular momentum must be transported in the more massive
the homogeneous evolution stops (Yoon et al. 2006).
model.

3.3. The Effect of Convective Core Mass

3.4. The Effect of Efficient Rotational Mixing
During the main-sequence evolution of a massive star, the con-
vective core contracts while the radiative envelope expands. The When a star rotates initially at high velocities, quasi-chemically
conservation of angular momentum will therefore demand that homogeneous evolution can occur. During such evolution, ro-
in the absence of any internal angular momentum transport, the tational mixing is so efficient that any chemical gradient be-
convective core and radiative envelope increases and decreases tween core and envelope cannot develop, meaning that the ra-
respectively their overall angular velocity (i.e. the core "spins diative envelope does not expand and the star’s radius remains
up" while the envelope "spins down"). This tells us that to en- roughly constant (Maeder 1987; Yoon et al. 2006). However
force solid body rotation during core contraction and envelope because the luminosity of a quasi-chemically homogeneously
expansion, angular momentum must be transported from the evolving star approaches the Eddington limit, the critical veloc-
core to the envelope. This is achieved by magnetic interactions ity of such a star does decrease (through Eq. 2). Furthermore the
which transport angular momentum along the angular velocity increased luminosity causes a strong increase in the mass-loss
gradient within the star. rate, meaning that the equatorial velocities of stars undergoing
Let us now consider two extreme examples. In a rotating star quasi-chemically homogeneous evolution are likely to decrease
with a negligible envelope mass, the core will dominate the an- with time. Thus such stars will evolve with a decreasing criti-
Article number, page 5 of 14
 A&A proofs: manuscript no. Be_paper

19 0.6
a 0.7 b 0.7 c
19 0.4
0.6 0.6
18
0.5 18 0.5 0.2
log(j/cm2s 1)

log(j/cm2s 1)
17 0.4 0.4 0.0

j/j1
17

X
0.3 0.3 0.2
16
0.4
5 M 0.2 16 15 M 0.2 5M
15 SMC SMC
Vinit=370kms 1 0.1 15 Vinit=470kms 1 0.1 0.6 15M
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0
m/M m/M m/M
Fig. 4. a-b: The specific angular momentum (solid lines) and hydrogen mass fraction (dotted lines) profiles for two SMC models of masses 5 and
15 M and initial equatorial velocities of 370 and 470 kms−1 respectively. These are the same models plotted in panels b & d of Fig. 2 and both
models have initial critical velocity fractions of around 0.6. Profiles are plotted for models where central helium mass fraction is 0.45 (black) , and
0.91 (orange). c: For each model in panel a & b the fractional difference in specific angular momentum is plotted between the two times. The blue
line represents the 5 M model,the orange the 15 M model. The dotted red line gives a reference for no angular momentum transport. The X-axis
in all plots indicates the fractional mass co-ordinate.

cal velocity fraction while quasi-chemically homogeneous evo- From the bottom panels of Fig. 5, we expect that in the
lution occurs. Milky-way, nitrogen is never enhanced by much more than a
This behaviour is shown for 25 M models in Figs. 2 e,f. It factor of 10 for models that rotate near the critical velocity. On
can be seen that the highest critical velocity fractions occur dur- the other hand, the LMC and SMC models that attain near criti-
ing the early part of the stars’ lifetimes because the critical ve- cal rotation velocities show surface nitrogen enhancements of at
locities (dashed lines) decrease relatively slowly while the equa- least a factor 10 and up to approximately a factor 30.
torial velocities (dot-dashed lines) fall due to angular momentum We therefore judge that single Be stars in the LMC should
loss. In the MW model, quasi-chemically homogeneous evolu- have surface nitrogen abundances  = 12 + log(N/H) no smaller
tion is shutdown when the star reaches an age of around 80% than around 7.7, and in the SMC no smaller than around 7.4.
of the hydrogen burning lifetime and from this point on the star In the Milky-way, we no not expect the single Be stars to have
evolves normally and advances towards the critical velocity. The outstanding nitrogen surface abundances.
phase of quasi-chemically homogeneous evolution ends because
as the rotational velocity decreases, rotational mixing processes
4. Population Synthesis Results
become less efficient and a chemical gradient eventually devel-
ops in the star which provides a barrier to mixing through buoy- 4.1. Predicted Fractions of Be Stars
ancy forces (Heger et al. 2000) and effectively "turns off" quasi-
chemically homogeneous evolution. In this section, we discuss synthetic populations of coeval rotat-
ing single stars as described in Sec. 2.2. These results can then
For a star to evolve to a high critical velocity fraction, it must
be compared to the number of Be stars observed in young star
have a significant initial rotation velocity but also be rotating
clusters of various ages. From our models, we derive the frac-
slowly enough to avoid quasi-chemically homogeneous evolu-
tion of Be stars within one bolometric magnitude (assumed to
tion. As the minimum rotation rate required for quasi-chemically
be equal to one visual magnitude) of the turn-off. We consider a
homogeneous evolution decreases with increasing mass (Yoon
stellar model to correspond to a Be star when its rotational ve-
et al. 2006), very massive stars rotate at high critical velocity
locity exceeds a predefined fraction of critical rotation. Figure 6
fractions for only very short fractions of their lifetimes, thus ex-
shows the result as a function of age for various threshold critical
plaining the rarity of Oe stars.
velocity fractions and metallicities.
A striking feature of this plot is the maximum in Be frac-
3.5. Nitrogen Enrichment tion for all metallicities near 10 Myr. At t = 0 there are no Be
stars because, as discussed earlier, the initial rotation distribu-
Here we address the question of whether or not Be stars formed tion prevents O stars entering the ZAMS with 3rot /3crit fractions
through a single star evolving towards the critical velocity are greater than around 0.7. The 30 M models take approximately
expected to show significant surface nitrogen enrichment, where 5 Myr to evolve towards critical rotation, shortly after which
nitrogen is the product of hydrogen burning and is brought to point they leave the main-sequence. From 5 to 10 Myr the Be
the stellar surface through rotationally induced mixing. Fig. 5 fraction grows sharply as angular momentum loss from winds
shows the evolution of both surface nitrogen abundance and crit- diminishes. From 10 to 20 Myr the Be fraction falls because the
ical velocity fraction as a function of the fractional hydrogen- of the core-mass effect as discussed in Section 3.3. Comparing
burning lifetime. Displayed are models with initial masses 5, 15 the hydrogen burning lifetimes of non-rotating MW models, it is
and 25 M with SMC, LMC and MW metallicities. As discussed found that populations with ages from 10 to 20 Myr have a main
by Brott et al. (2011) the relative increase of the nitrogen abun- sequence turn off mass of around 17 M . These models are in
dance goes down with increasing metallicity, therefore we see a "Goldilocks" situation where they are massive enough to have
weaker nitrogen enrichment in the MW models than the LMC or an appreciable convective core but not so massive to lose large
SMC models. It is also clear that rotationally induced mixing is amounts of angular momentum.
more efficient in more massive stars, owing to the effects of in- Furthermore we see that the Be fraction increases at earlier
creased radiation pressure in more massive stars (Maeder 1987; times for the lower metallicity models and that there is a clear
Yoon et al. 2006). trend in metallicity which shows that single Be stars become
Article number, page 6 of 14
 Hastings, Wang, and Langer: The Single Star Path to Be Stars

45
a c e
40 5 M¯ 15 M¯ 25 M¯
35 SMC SMC SMC
N enrichment factor

30
25
20 ² = 7.8
15
10 ² = 7.5 1.0
5 ² = 7.0
45
b d f
40 5 M¯ 15 M¯ 25 M¯ 0.8
35 LMC LMC LMC
N enrichment factor

30
25
20 ² = 8.1 0.6

vrot /vcrit
15
10 ² = 7.8
5 ² = 7.3 0.4
45
g h i
40 5 M¯ 15 M¯ 25 M¯
35 MW MW MW
N enrichment factor

0.2
30
25
20 ² = 8.94
15 0.0
10 ² = 8.64
5 ² = 8.14
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
t/tMS t/tMS t/tMS
Fig. 5. The surface nitrogen enrichment factor, computed as the nitrogen mass fraction divided by the initial nitrogen mass fraction as a function
of the fractional hydrogen-burning lifetime, t/t MS for models with initial rotational velocities between 0 and 600 km s−1 and initial masses 5, 15,
25 M as marked in the plot. SMC, LMC and MW compositions are shown in the top, middle and bottom panels respectively. The colour scale
corresponds to the critical rotation fraction, 3rot /3crit . For each metallicity, various nitrogen abundances,  = 12 + log(N/H) are displayed as dotted
black lines with the value given in the left column plots.

more common with decreasing metallicity. Both of these fea- rotating stars are located very close to the turn-off, as can be ex-
tures are a result of the fact that the strength of angular momen- pected by the fact that our models only achieve high fractions of
tum loss is metallicity dependant as discussed in Section 3.2. critical rotation near core hydrogen exhaustion (see Figs. 2,A.1).
It is also found that the Be fraction is strongly depen- Fig. 7 also shows the effects of gravity darkening, with the slowly
dant on the chosen Be criterion. When critical velocity frac- rotating models almost confined to a single isochrone while the
tion, 3rot /3crit of 0.7 is chosen as the Be criterion, we predict fast rotators suffer strong gravity darkening and display a wider
Be fractions in the Magellanic Clouds in the range 15 to 35%. range of colours, due to a relatively large spread in effective tem-
Wheres when we restrict Be stars to being nearly critical rota- peratures. Appendix B shows the same colour-magnitude dia-
tors ( 3rot /3crit > 0.98), the Magellanic Cloud Be fraction lies in gram but ignoring the effects of gravity darkening for compari-
the range 0 to 10% and is almost 0 for population ages greater son.
than 50 Myr. In the following analysis we shall assume that the Hα emit-
ters in NGC 330 are Be stars. The observations in Fig. 7 show
4.2. Fast Rotators in the Colour-Magnitude Diagram that most of the observed Hα emitters are redder than the ordi-
nary main sequence stars. This segregation is not a unique fea-
Using the procedure outlined in Section 2.2 we build synthetic ture to NGC 330, with many LMC and SMC clusters exhibiting
colour-magnitude diagrams to indicate the expected positions of the same trait (Milone et al. 2018). Telting et al. (1998) have
fast rotators. Fig. 7 shows the colour-magnitude diagram posi- suggested through observations of Be stars whose spectra show
tions and critical velocity fractions of our SMC models with a rapid switching between containing emission lines and not, that
coeval age of 35 Myr. Over plotted on our theoretical predictions the decretion disc can contribute up to 40% of a Be star’s flux.
are Hubble Space Telescope observations of the SMC cluster Because our models do not include the contribution of a Be star’s
NGC 330 (Milone et al. 2009), We see that the nearly critically decretion disc to the observed fluxes, it is not meaningful to
Article number, page 7 of 14
 A&A proofs: manuscript no. Be_paper

M /M
turn − off ¯
27.0 17.3 10.7 8.6 6.8 5.5
0.35
vrot /vcrit > 0.7 vrot /vcrit > 0.8
0.30
0.25
0.20
frac.

0.15
0.10
0.05
0.00
0.35 MW
vrot /vcrit > 0.9 LMC vrot /vcrit > 0.98
0.30 SMC
0.25
0.20
frac.

0.15
0.10
0.05
0.00
2 10 50 100 2 10 50 100
t/Myr t/Myr
Fig. 6. The fractions of stars rotating faster than various values of the critical rotation fraction, 3rot /3crit , as given in the top left corner of each plot,
in our synthetic coeval single star populations as function of their age. Considered only are stars brighter than one bolometric magnitude below
the main-sequence turn-off. Metallicities are displayed as MW (black), LMC (blue) and SMC (green). For the top left panel the top scale gives the
main-sequence turn off mass for non-rotating MW models. Observations with error bars from Iqbal & Keller (2013) are shown as green and blue
stars for SMC and LMC observations respectively.

compare the colours of our fast rotating models with those of threshold for considering the models a Be star, or by adopting
observed Be stars. However, assuming that the error on the mag- larger initial rotation velocities. However, it is an intrinsic fea-
nitudes of our synthetic Be stars is at most 0.35 mag, we may ture of the rotating single star models to increase the ratio of
consider the relationship between the relative number of Be stars rotation to critical rotation velocity with time (Fig. 5; see als Ek-
and apparent magnitude. ström et al. 2008) . Thus, it appears quite unlikely that single star
Fig. 8 shows the Be fraction of our model predictions in evolution that the observed distribution of Be stars in NGC 330
F814W apparent magnitude bins for various threshold rotation can be explained solely by single star evolution.
rates for stellar models to be considered a Be star. We see that
our model Be stars are strongly biased to being located near the
main-sequence turn off, around mF814W = 15.4 mag, where de- 5. Discussion
pending on how fast we require a Be star to rotate, the Be fraction
5.1. Uncertainties
is between 30 and 90%. Fig. 8 compares our model predictions
to Be star counts in NGC 330 (Milone et al. 2018). The observa- It is important to keep in mind that models are simply that and
tions show that in NGC 330 the Be fraction is approximately 0.6, at some point they must fail to reflect the behaviour of real stars.
and it remains rather constant within a wide magnitude range Martins & Palacios (2013) have found that the hydrogen burning
of 17.0 to 15.6 mag. This result is quantitatively confirmed by lifetimes given by models studied here differ to those of MESA
a recent study of Bodensteiner et al. (2019), who used MUSE models by approximately 15%. As this is merely a discrepancy
spectroscopy to identify Be stars in the core of NGC 330. in the clock, but not in the physical behaviour of the models
Comparing our models with the observations, Fig. 8 reveals (Marchant 2017), this is concerning, but should not change the
that when assuming that Be stars are required to rotate only at main results presented here.
70% or more of the critical rotation velocity, our single star mod- Another issue is the treatment of mass-loss, which is a strong
els agree with the observed Be star count in NGC 330 at the turn factor in determining a model’s evolution towards the critical ve-
off region. However, they strikingly fail in two respects. Firstly, locity. The mass-loss prescription used (Vink et al. 2000) was
even adopting the least stringent threshold value for Be stars, our calibrated for models in the range 15-20 M , so it may be ques-
models fail to produce the large number of observed Be stars. tionable whether this scheme is accurate for models outside this
Secondly, our models can not reproduce the fact that the Be star range. Furthermore the correct treatment of a star rotating close
fraction in NGC 330 is constant over a range of 1.5 magnitudes. to the critical velocity is complex. Near the critical velocity,
Here, the second failure seems the worst. The total number two distinct winds are expected to form, a cold equatorial wind
of Be star can in principle be boosted by lowering the rotation (which carries away angular momentum) and a warm polar wind
Article number, page 8 of 14
 Hastings, Wang, and Langer: The Single Star Path to Be Stars

14.0 H α emitters
Stars 15.6
15.8 1.0
15.0
16.0 0.8
16.0 16.2
mF814W

mF814W
0.6

vrot /vcrit
16.4
17.0 0.4
16.6
16.8 0.2
18.0
17.0 0.0
-2.0 -1.0 0.0 -1.4 -1.2 -1.0
mF336W − mF814W mF336W − mF814W
Fig. 7. The synthetic colour-magnitude diagram of a 35Myr star cluster at SMC metallicity, where each dot represents one single star and the
colour gives the critical velocity fraction, 3rot /3crit as indicated by the colour bar. Gravity darkening is included assuming a random orientation of
the rotation axis. Over plotted are observations of SMC cluster NGC 330 (Milone et al. 2018), with H α emitters marked by orange squares and
normal stars as orange star symbols. The right panel shows the region indicated by the red box in the left panel. To convert the models to apparent
magnitudes a distance modulus of 18.92 mag and a reddening of 0.06 mag were used.

1.0 vrot /vcrit >0.7

ent which wind has the dominating effect, although the winds
of models presented here always carry away angular momen-
0.9 vrot /vcrit >0.8 tum. Eq. 1 demands that as a star approaches the critical rotation
0.8 vrot /vcrit >0.9 velocity, the mass loss rate becomes infinite. It can also be ques-
vrot /vcrit >0.95
0.7 tioned whether this is a correct treatment of the mass loss of a
Be/(Be + B)

critically rotating star.

0.6 In the SMC, Rímulo et al. (2018) matched light-curve mod-
0.5 els to observations of 54 Be stars and determined that the typical
0.4 mass-loss rate of a Be star in the mass range 10 to 20 M is of
the order 10−10 M yr−1 . Our models predict that only stars with
0.3 initial masses less than around 10 M have such mass-loss rates
0.2 on the main-sequence. This discrepancy highlights the fact that
0.1 the models studied here are 1-dimensional and hence may strug-
0.0 18.0 17.0 16.0 15.0
gle to represent accurately the mass-loss of a rotating star (which
is a 2-dimensional problem).
mF814W Our definition of a Be star is one that is rotating close to
the critical velocity. However it would appear that nature has a
Fig. 8. A histogram of the Be fraction of stars in SMC cluster NGC 330 slightly different definition, with pulsations perhaps playing a
(Milone et al. 2018) as a function of apparent magnitude in the Hubble role (Rivinius et al. 2001; Neiner & Mathis 2014). Pulsations
Space Telescope Wide-field Camera 3 filter F814W, mF814W plotted as could serve to kick matter off of the stellar surface and aid the
white columns. It is assumed that H α emitters are Be stars. Over plotted
formation of a circumstellar disc. Observations with the CoRoT
as coloured lines are our model predictions in the same magnitude bins
for various minimum rotation rates to be counted as a Be star, as given space telescope show that Be stars display pulsations that can
in the legend. transport angular momentum through the star (Huat et al. 2009),
thus affecting the evolution of rotation velocities. The interaction
of rotation and pulsations is out of the scope of this work and
the results of such endeavours are eagerly awaited. Furthermore
(which carries away less angular momentum). It is not appar- recent observations with the TESS space telescope imply that
Article number, page 9 of 14
 A&A proofs: manuscript no. Be_paper

the Be star disc could be fed by mass ejections from starspots tions increase with decreasing metallicity, in agreement with our
(Balona & Ozuyar 2019). Be stars are complex objects, with ro- models.
tation being a key ingredient of the Be phenomenon, but perhaps Both Maeder et al. (1999) and Martayan et al. (2010) found
not the only one. that Be stars are three to five times more frequent in the SMC
Our population synthesis results are dependant on the ini- than the galaxy, again in fairly good agreement with the model
tial rotation distribution that is assumed. Whereas Dufton et al. predictions. Martayan et al. (2010) reported that the distribution
(2013) used high quality data, a large and unbiased sample of of Be star frequency across spectral types does not depend on
star, and corrected for effects as macro turbulence, other rota- metallicity. As Fig. 6 shows similar trends for all metallicities,
tional velocity distributions are available (Martayan et al. 2006; this behaviour is confirmed by the models.
Hunter et al. 2008; Huang et al. 2010)While our quantitative re- Observations from Tarasov (2017) show that Be stars be-
sults might change using one of them, we do not expect a change come most common in clusters with ages of 12-20Myr, in rel-
in the qualitative behaviour of our results, which appears to be atively good agreement to the model predictions.
determined by the evolutionary factors discussed in Sect. 3 Golden-Marx et al. (2016) found that in the SMC the fre-
quency of Oe stars is strongly peaked around spectral types O9.
Furthermore the Oe to O star fraction was measured as 0.26,
5.2. Comparison with Previous Models compared to 0.03 for the MW. This measurement supports our
result that very few stars with high critical velocity fractions and
The frequencies of single Be stars have been predicted from
ages less than 10Myr should be found in the MW, but are found
models by Ekström et al. (2008), Granada et al. (2013) and
at lower metallicities (see lower Panels of Fig. 6).
Granada et al. (2016). A major difference between these sets of
Owing to the fact that the results of our population synthe-
models and the models studied here is the inclusion of the effects
sis rely strongly on the adopted initial rotation distribution, one
of an internal magnetic field which strongly couples the core and
may question whether it is appropriate to assume that stars in the
envelope, which may increase the predicted numbers of critically
MW and Magellanic Clouds have equivalent initial critical rota-
rotating stars.
tion fractions. Whether the observed trends in Be fraction with
Ekström et al. (2008) predict at solar metallicity and an age metallicity are due to stellar evolution or a metallicity dependant
of 25 Myr the fraction of stars with a brightness of up to two rotation distribution is not clear. If lower metallicity stars were to
magnitudes below the turn-off and rotating at the critical velocity rotate significantly faster, rotationally enhanced mass-loss would
is 5%, compared to 0% found here. However when one looks at hinder the formation of Be stars at lower metallicities, therefore
the fraction of stars rotating faster than 3rot /3crit =0.7 at 20 Myr, there is a limit to how much faster stars at lower metallicities
Ekström et al. (2008) finds a Be fraction of around 15%, com- can rotate. Even though Penny et al. (2004) found no signifi-
pared to 8% found here. Both sets of models agree that at ages cant difference between rotational velocities of O-type stars in
greater than around 40 Myr, almost no stars rotate at the critical the galaxy and Magellanic Clouds, Keller (2004) concludes that
velocity. Ekström et al. (2008) used a Gaussian-like initial rota- LMC stars are more rapidly rotating than galactic stars. It is cu-
tion distribution with a peak at Ω/Ωcrit = 0.6, which in the Roche rious to consider that even if LMC stars have faster equatorial
model corresponds to 3rot /3crit ≈ 0.4, and therefore is judged to be velocities, due to their compactness they have larger critical ve-
a similar initial rotation distribution to the one used in this work. locities, and hence perhaps the same initial 3rot /3crit as galactic
As demonstrated by Ekström et al. (2008) mass loss rates play a stars.
crucial role in the evolution of the surface rotation. In our models
above 10 M , a turn over in 3rot /3crit is caused by strong mass-
loss at late times (see Figure 2). When using the Kudritzki & 5.4. Comparing the Single and Binary Star Formation
Puls (2000) mass-loss scheme, Ekström et al. (2008) found this Channels of Be Sstars.
turn over not to occur, and so those models spend more time at
Whereas in this work we focus on the single star formation chan-
high 3rot /3crit values and so Be stars become more common. It is
nel for Be stars, it is evident that Be stars can also be formed
therefore concluded that mass-loss is just as important as angular
through close binary evolution. The main mechanism is spin-
momentum transport in producing stars which rotate close to the
up by accretion, which is expected to occur as a consequence
critical velocity at galactic metallicity. At Z = 0.002 (a metallic-
of mass transfer (Langer 2012). The class of Be/X-ray binaries
ity similar to our SMC models) Ekström et al. (2008) calculate a
(Reig 2011) provides strong support for this picture. Pols et al.
maximum in Be fraction at an age of 10Myr of around 10%, in
(1991) showed through detailed models that Be stars may be pro-
good agreement with the results presented here.
duced by mass accretion from a companion star in the course of
Both Ekström et al. (2008) and Granada et al. (2013) find that close binary evolution. Using simplified binary evolution calcu-
Be stars should become rarer at lower metallicities, which is in lations, Shao & Li (2014) demonstrated that potentially a large
contradiction to the results presented here and the general trend enough number of them could emerge from binary evolution to
seen by observers (Iqbal & Keller 2013; Maeder et al. 1999; Mar- explain the currently observed Be star populations. In the follow-
tayan et al. 2010). ing, we discuss several basic differences which can be expected
between Be stars formed through binary interaction compared to
5.3. Comparison with Further Observations those formed through the single star channel.
As we have seen above, when using the rotational veloc-
Iqbal & Keller (2013) observed Be star fractions within 1 vi- ity distribution of Dufton et al. (2013), it is difficult for single
sual magnitude of the turnoff in LMC and SMC clusters. Such star models to achieve very close to critical rotation (lower right
observations are directly comparable with Fig. 6, with the data panel of Fig. 6). To remedy this would require that a significant
over-plotted on our predictions. In clusters of ages from 7 to 8 fraction of stars is already born with near critical rotation. For
Myr the Be fraction decreases from 15 % to 0 %. When defining low enough mass or metallicity, this group of stars would the re-
the Be criterion to be 0.9 3rot /3crit , a similar behaviour is found main rapidly rotating throughout their main sequence life. Such
albeit at later ages. Iqbal & Keller (2013) also find that Be frac- a picture appears not to be supported by observations (McSwain
Article number, page 10 of 14
 Hastings, Wang, and Langer: The Single Star Path to Be Stars

& Gies 2005). Accordingly, single star evolution appears to be was found to be almost devoid of nitrogen lines and possessing a
able to explain significant Be star populations only if decretion spectrum inconsistent with single star evolution models (Lennon
discs could form also in stars rotating significantly below criti- et al. 2005), giving further evidence that Be stars can be formed
cal. by binary interactions.
In mass transferring binary evolution models, this is dif- Finally we note that a key difference between Be stars pro-
ferent. The angular-momentum of the mass-gainer increases duced by the two channels concerns their potential binary com-
quickly, such that critical rotation can be achieved after a rel- panions. Since the initial rotational velocity distribution for sin-
ative mass increase of the mass gainer of 10% or less (Packet gle stars and stars in binaries appears to be similar (Ramírez-
1981). During mass-transfer there is no fine-tuning mechanism Agudelo et al. 2015), we would expect a significant fraction
that switches off accretion when a given rotation rate is reached, of Be stars formed through the single star channel (i.e., with-
instead the only limit is critical rotation. Therefore, all mass out accretion-induced spin-up) to have unevolved main sequence
transferring binaries where tides do not limit the spin-up of the companions. However, essentially no such stars are known. Vice
mass gainer — which is the vast majority — are expected to versa, whereas massive binary-produced Be stars may be sin-
produce a critically rotating main sequence star (Langer 2012). gle since they lost their companion when it produced a super-
After the accretion phase, the two mechanisms which affect the nova explosion, the lower mass binary-produced Be stars should
single stars, i.e., spin-down by mass loss and spin-up due to core all have evolved companions: subdwarfs or white dwarfs. While
contraction, will also work in the spun-up mass gainer. Whereas those are very hard to detect (Schootemeijer et al. 2018), recent
the wind induced drain of angular momentum may spin down studies of disc truncation of apparently single Be stars suggest
some of the most massive mass gainers, the core contraction ac- that indeed unseen companions are present in the majority of
companying central hydrogen burning ensures that most of them cases (Klement et al. 2019).
remain at critical rotation for the rest of their main sequence evo-
lution.
Consequently, whereas single star evolution leads to an in- 6. Conclusions
crease of the rotation velocity compared to its critical value in
many cases, the binary channel can produce a much larger num- We have identified and discussed three factors which affect a
ber of stars living at critical rotation for a long time, compared to star’s evolution towards the critical velocity throughout main-
the single star channel. Furthermore, Chen et al. (2020) find that sequence evolution. Mass-loss through stellar winds has the ef-
the initial mass ratio limit for stable mass transfer leads naturally fect of removing angular momentum from a star, and so hinders
to a restriction of binary-produced Be stars to within about two the approach to the critical velocity. The fraction of convective
magnitudes of the cluster turn-off, which compares well with ob- core mass to the total stellar mass strongly affects the internal
servations when interpreting the Hα emitting stars in NGC 330 angular momentum transport, which is crucial for an expanding
as Be stars (Milone et al. 2018). envelope to maintain a fast rotational velocity. Lastly the occur-
A further important difference between the Be stars pro- rence of quasi-chemically homogeneous evolution prevents the
duced via single and binary evolution concerns their expected stellar envelope from expanding and thus critical velocity de-
surface abundances. As discussed above, the mass-gainer of a creasing, and also increases the angular momentum lost through
binary system may only accrete a small amount of mass to spin stellar winds.
up. In this case, only material from the outer envelope of the When using an observed distribution of B star rotational ve-
donor star is incorporated into the mass gainer. As this mate- locities for constructing synthetic stellar populations, we find
rial is generally not enriched in hydrogen burning products, one that our single star models predict few stars rotating at near crit-
would expect that Be stars formed through the binary channel are ical velocities, although we do predict as much as 35% of OB
not significantly polluted by accretion. Detailed binary evolution stars to rotate with more than 70% of their critical velocity. We
models with LMC metallicity (Langer & et al 2019) suggest that therefore conclude that if Be stars are near-critical rotators, then
the surface nitrogen mass fraction of spun-up mass gainers is at single star models cannot explain the observed numbers of Be
most tripled compared to the baseline nitrogen abundance. Ad- stars. In this case, most Be stars must be the product of mass-
ditionally, the spun-up mass gainers had ordinary rotation rates transfer in binary systems.
before the mass transfer episode. As such, they established a If Be stars instead only rotate at say 70-80% of their criti-
strong mean molecular weight barrier between core and enve- cal velocity, then the observed Be star fractions can be reason-
lope, which prevents any significant rotationally induced mixing ably described by single star evolution (see Fig 6). However, in
after their spin-up. the ∼40 Myr old SMC cluster NGC 330, Be stars are observed
In Sec. 3.5 we argue that single Be stars ought to have much in significant numbers down to almost two magnitudes below
larger surface nitrogen enhancements (by as much as a factor the main-sequence turn-off. Independent of the rotation thresh-
30; cf., Fig. 5). A diagnostic to discriminate single and binary Be old for the Be phenomenon, our single star models predict that
stars would then be their surface nitrogen abundances. Dunstall Be stars should be located only in a narrow luminosity range
et al. (2011) find that in the LMC cluster NGC 2004, only two near the turn-off (see Fig 8), which disagrees with observations
Be stars from a sample of 11 were measured to have a nitrogen of NGC 330.
abundance  = 12 + log(N/H) greater than 7.8, while the other 9 Whereas significant uncertainties remain, specifically in rec-
Be stars had  less than 7.4. This observed bimodal distribution onciling how stars that appear to rotate at 70% of their critical
supports the idea of the binary and single Be star formation chan- velocity can still form decretion discs and why so few Be stars
nels producing populations with different nitrogen enrichments, are observed to rotate near the critical velocity, it appears evident
and would suggest that in NGC 2004 the binary formation chan- that the observed Be star populations can not be explained by
nel dominates. Also Dunstall et al. (2011) found that the nitro- single star evolution alone, and that it may not be the dominant
gen abundances of the majority of the LMC Be stars observed channel for Be star formation. Nevertheless, single star evolu-
in the VLT-FLAMES Survey of Massive Stars are not consis- tion will contribute, most strongly so in the age range from 8 to
tent with single star evolution. Also, the Be star NGC 330-B 12 20 Myr, at least at sub-solar metallicity.
Article number, page 11 of 14
 A&A proofs: manuscript no. Be_paper

Furthermore, our single star models predict that the sur- Raguzova, N. V. & Popov, S. B. 2005, Astronomical and Astrophysical Transac-
faces of rapidly rotating single stars should be contaminated with tions, 24, 151
freshly synthesised nitrogen, the more the faster the rotation. In Ramírez-Agudelo, O. H., Sana, H., de Mink, S. E., et al. 2015, A&A, 580, A92
Ramírez-Agudelo, O. H., Simón-Díaz, S., Sana, H., et al. 2013, A&A, 560, A29
spun-up mass gainers of binary systems, this is not necessarily Reig, P. 2011, Ap&SS, 332, 1
so. The observations of non- or weakly nitrogen enriched sur- Rímulo, L. R., Carciofi, A. C., Vieira, R. G., et al. 2018, MNRAS, 476, 3555
faces in several groups of Be stars therefore strengthens the con- Rivinius, T., Baade, D., Štefl, S., et al. 2001, A&A, 369, 1058
clusion that the majority of these objects can not originate from Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&A Rev., 21, 69
Schootemeijer, A., Götberg, Y., de Mink, S. E., Gies, D., & Zapartas, E. 2018,
single star evolution. A&A, 615, A30
Secchi, A. 1866, Astronomische Nachrichten, 68, 63
Acknowledgements. We are grateful to Dietrich Baade, Shao, Y. & Li, X.-D. 2014, ApJ, 796, 37
Danny Lennon and Christophe Martayan for discussions and for Spruit, H. C. 2002, A&A, 381, 923
Struve, O. 1931, ApJ, 73, 94
pointing out the recent TESS result to us. We also thank our Suijs, M. P. L., Langer, N., Poelarends, A.-J., et al. 2008, A&A, 481, L87
anonymous referee for useful comments on an earlier version of Tarasov, A. E. 2017, in Astronomical Society of the Pacific Conference Series,
this manuscript. Vol. 510, Stars: From Collapse to Collapse, ed. Y. Y. Balega, D. O. Kudryavt-
sev, I. I. Romanyuk, & I. A. Yakunin, 132
Telting, J. H., Waters, L. B. F. M., Roche, P., et al. 1998, MNRAS, 296, 785
References Townsend, R. H. D., Owocki, S. P., & Howarth, I. D. 2004, MNRAS, 350, 189
van Bever, J. & Vanbeveren, D. 1997, A&A, 322, 116
Baade, D., Rivinius, T., Štefl, S., & Kaufer, A. 2002, A&A, 383, L31 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2000, A&A, 362, 295
Balona, L. A. & Ozuyar, D. 2019, arXiv e-prints, arXiv:1911.03068 von Zeipel, H. 1924, MNRAS, 84, 665
Bodensteiner, J., Sana, H., Mahy, L., et al. 2019, arXiv e-prints, Yoon, S.-C. & Langer, N. 2005, A&A, 443, 643
arXiv:1911.03477 Yoon, S.-C., Langer, N., & Norman, C. 2006, A&A, 460, 199
Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 Zorec, J. & Briot, D. 1997, A&A, 318, 443
Chen, W., Langer, N., Schootemaijer, A., et al. 2020, ApJL, submitted Zorec, J., Frémat, Y., Domiciano de Souza, A., et al. 2016, A&A, 595, A132
de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ,
764, 166
Dufton, P. L., Dunstall, P. R., Evans, C. J., et al. 2011, ApJ, 743, L22
Dufton, P. L., Langer, N., Dunstall, P. R., et al. 2013, A&A, 550, A109
Dunstall, P. R., Brott, I., Dufton, P. L., et al. 2011, A&A, 536, A65
Ekström, S., Meynet, G., Maeder, A., & Barblan, F. 2008, A&A, 478, 467
Espinosa Lara, F. & Rieutord, M. 2011, A&A, 533, A43
Girardi, L., Bertelli, G., Bressan, A., et al. 2002, A&A, 391, 195
Golden-Marx, J. B., Oey, M. S., Lamb, J. B., Graus, A. S., & White, A. S. 2016,
ApJ, 819, 55
Granada, A., Ekström, S., Georgy, C., et al. 2013, A&A, 553, A25
Granada, A., Georgy, C., Haemmerlé, L., Ekström, S., & Meynet, G. 2016, in
Astronomical Society of the Pacific Conference Series, Vol. 506, Bright Emis-
saries: Be Stars as Messengers of Star-Disk Physics, ed. T. A. A. Sigut & C. E.
Jones, 33
Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368
Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350
Huang, W., Gies, D. R., & McSwain, M. V. 2010, ApJ, 722, 605
Huat, A.-L., Hubert, A.-M., Baudin, F., et al. 2009, A&A, 506, 95
Hunter, I., Lennon, D. J., Dufton, P. L., et al. 2008, A&A, 479, 541
Iqbal, S. & Keller, S. C. 2013, MNRAS, 435, 3103
Keller, S. C. 2004, PASA, 21, 310
Kippenhahn & Weigert. 1990, Stellar Structure and Evolution (Springer-Verlag)
Klement, R., Carciofi, A. C., Rivinius, T., et al. 2019, arXiv e-prints,
arXiv:1909.12413
Kriz, S. & Harmanec, P. 1975, Bulletin of the Astronomical Institutes of
Czechoslovakia, 26, 65
Kudritzki, R.-P. & Puls, J. 2000, ARA&A, 38, 613
Lamers, H. J. G. L. M., Snow, T. P., & Lindholm, D. M. 1995, ApJ, 455, 269
Langer, N. 1998, A&A, 329, 551
Langer, N. 2012, ARA&A, 50, 107
Langer, N. & et al. 2019, in prep.
Lennon, D. J., Lee, J. K., Dufton, P. L., & Ryans, R. S. I. 2005, A&A, 438, 265
Li, G.-W., Shi, J.-R., Yanny, B., et al. 2018, ApJ, 863, 70
Liu, Q. Z., van Paradijs, J., & van den Heuvel, E. P. J. 2006, A&A, 455, 1165
Maeder, A. 1987, A&A, 178, 159
Maeder, A., Grebel, E. K., & Mermilliod, J.-C. 1999, A&A, 346, 459
Marchant, P. 2017, The impact of tides and mass transfer on the evolution of
metal-poor massive binary stars, PhD thesis, Rheinische Friedrich-Wilhelms-
Universität Bonn
Martayan, C., Baade, D., & Fabregat, J. 2010, A&A, 509, A11
Martayan, C., Frémat, Y., Hubert, A. M., et al. 2006, A&A, 452, 273
Martins, F. & Palacios, A. 2013, A&A, 560, A16
McSwain, M. V. & Gies, D. R. 2005, ApJS, 161, 118
Milone, A. P., Bedin, L. R., Piotto, G., & Anderson, J. 2009, A&A, 497, 755
Milone, A. P., Marino, A. F., Di Criscienzo, M., et al. 2018, MNRAS, 477, 2640
Neiner, C. & Mathis, S. 2014, in IAU Symposium, Vol. 301, Precision Astero-
seismology, ed. J. A. Guzik, W. J. Chaplin, G. Handler, & A. Pigulski, 465–
466
Packet, W. 1981, A&A, 102, 17
Penny, L. R., Sprague, A. J., Seago, G., & Gies, D. R. 2004, ApJ, 617, 1316
Pols, O. R., Cote, J., Waters, L. B. F. M., & Heise, J. 1991, A&A, 241, 419

Article number, page 12 of 14

 Hastings, Wang, and Langer: The Single Star Path to Be Stars

0.8

0.7
initial vrot /vcrit

0.6

0.5

0.4
5.0M¯ 15.0M¯ 25.0M¯
0.3 SMC SMC SMC

0.8

0.7
initial vrot /vcrit

0.6

0.5 0.6<vrot /vcrit <0.7

0.7<vrot /vcrit <0.8
0.4
0.8<vrot /vcrit <0.9
5.0M¯ 15.0M¯ 0.9<vrot /vcrit <0.98 25.0M
MW MW 0.98<vrot /vcrit <1.0 MW
¯

0.3
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
t/tMS t/tMS t/tMS
Fig. A.1. The critical velocity fractions, 3rot /3crit as a function of fractional hydrogen-burning lifetime, t/t MS and initial critical velocity fractions
for models with masses of 5, 15 and 25 M and SMC (top panels) and MW (bottom panels) metallicities, as indicated in the figure. The colours
indicate the critical velocity fraction, 3rot /3crit as given in the legend.

Appendix A: Initial Conditions Required to Reach Near Critical Velocities

In Section 3.1 we have discussed the evolution of rotational velocities of various models, when the models all have the same relative
initial rotation rates. Here we explore how the approach to the critical velocity depends on the initial rotation rate by performing
interpolations between the models. Fig. A.1 shows, for differing stellar masses and metallicities, the critical velocity fraction as a
function of time and initial rotation rate. The colour of each point on the plot shows the critical velocity fraction at a particular
fractional hydrogen-burning lifetime, t/t MS , and initial critical velocity fraction value. By following horizontal lines in the plot, one
traces the evolution of a single model through its evolution.
For 5 M models at both SMC and MW metallicities, the star evolves generally towards higher critical velocity fractions. This
can be seen as one traces a horizontal line, one moves always into regimes of higher critical velocity fractions. The exception is the
5 M MW models with 3rot /3crit > 0.7, which at the end of hydrogen burning spin down through increased rotationally enhanced
mass-loss. Looking at 15 M models at SMC metallicity, one sees also that there is a constant evolution towards higher critical
velocity fractions. On the other hand 15 M models with MW metallicity evolve to higher critical velocity fractions until around
80% of the hydrogen-burning lifetime, then they spin down due to angular momentum loss through winds (as discussed in Sec. 3.2).
The 25 M models behave in a more complicated way because the initially very fast rotating models can undergo quasi-chemically
homogeneous evolution (as discussed in Sec. 3.4). In the right panels of Fig. A.1 one can see the two regimes according to whether
or not the critical velocity fraction is increasing or decreasing. We can see that for 25 M MW models with initial 3rot /3crit ≈ 0.7,
they evolve at first to lower critical velocity fractions then after t/t MS ≈ 0.6 they begin to evolve towards higher critical velocity
fractions. This occurs because initially the star is evolving quasi-chemically homogeneously, during which time rotation rates and
hence rotationally induced mixing efficiency drops until quasi-chemically homogeneous evolution is shut down, at which point
the star begins to evolve with an expanding envelope and approaches the critical velocity. For similarly initially fast rotating SMC
models the same behaviour does not occur due to the weaker stellar winds at lower metallicities.
Fig. A.1 also shows us that the models only reach critical rotation (the black areas in the figure) very near core hydrogen
exhaustion. Furthermore, when one looks at the minimum initial rotation rate required to reach near critical rotation, it decreases
with increasing mass due to angular momentum transport efficiency ( as discussed in Sec. 3.3).

Appendix B: Synthetic Colour-magnitude Diagram Without Gravity Darkening

In Fig B.1 we present the results of our population synthesis to simulate the colour-magnitude diagram of NGC 330 while ignoring
the effects of gravity darkening. It is seen that along the main sequence there is a one-to-one relation between current critical velocity
fraction and the mF336W − mF814W colour, with faster rotators being redder. Such a relation is destroyed by gravity darkening (see
Fig. 7).

Article number, page 13 of 14

 A&A proofs: manuscript no. Be_paper

14.0 H α emitters
Stars 15.6
15.8 1.0
15.0
16.0 0.8
16.0 16.2
mF814W

mF814W
0.6

vrot /vcrit
16.4
17.0 0.4
16.6
16.8 0.2
18.0
17.0 0.0
-2.0 -1.0 0.0 -1.4 -1.2 -1.0
mF336W − mF814W mF336W − mF814W
Fig. B.1. The synthetic colour-magnitude diagram of a 35Myr star cluster at SMC metallicity, where each dot represents one single star and the
colour gives the critical velocity fraction, 3rot /3crit as indicated by the colour bar. Gravity darkening is not included. Over plotted are observations
of SMC cluster NGC 330 (Milone et al. 2018), with H α emitters marked by orange squares and normal stars as orange star symbols. The right
panel shows the region indicated by the red box in the left panel. To convert the models to apparent magnitudes a distance modulus of 18.92 mag
and a reddening of 0.06 mag were used.

Article number, page 14 of 14