11institutetext: Argelander-Institut für Astronomie, Universität Bonn, Auf dem
Hügel 71, 53121 Bonn, Germany
22institutetext: Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69,
53121 Bonn, Germany
33institutetext: Dipartimento di Fisica e Astronomia “Galileo Galilei” —Univ.
di Padova, Vicolo dell’Osservatorio 3, Padova, I-35122, Italy

# A stringent upper limit on Be star fractions produced by binary interaction

B. Hastings 1122 N. Langer 1122 C. Wang 11 A. Schootemeijer 11 A. P. Milone 33

###### Abstract

Context. Binary evolution can result in fast-rotating stars, predicted to be
observable as Be stars, through accretion of angular momentum during mass-
transfer phases. Despite numerous observational evidence pointing to this
possibly being the dominant Be formation channel, current models struggle to
produce a satisfactory description of Be star populations.

Aims. Given distinct uncertainties in detailed binary evolution calculations,
we investigate a rigorous and model independent upper limit for the production
of Be stars through binary interaction and aim to confront this limit with
observations of Be stars in young star clusters.

Methods. Using extreme assumptions, we calculate the number ratio of post-
interaction to pre-interaction binary systems in a coeval population, which
describes an upper limit to Be star formation through mass-transfer. A
detailed comparison is made between our derived upper limit and relevant
observations of Be stars, which allows us to probe several aspects of binary
star physics.

Results. We find that in coeval populations, binary interaction can at most
account for one third of all main-sequence stars being Be stars. Near the
cluster turn-off region, this limit appears to be realised in the clusters
studied. Away from the turn-off, applying simple assumptions about which
systems undergo unstable mass-transfer produces a good fit to the observed Be
fraction as a function of mass.

Conclusions.  We find that assuming distinct physics, binary evolution alone
can in principle match the high numbers of Be stars observed in open clusters.
Whether the required binary physics is realised in nature remains to be
investigated.

###### Key Words.:

stars: binaries: general / stars: emission line, Be / stars: evolution /
stars: massive / stars: rotation

##  1 Introduction

Be stars are massive main-sequence stars which display emission features in
their spectra. While, since their discovery over 150 years ago (Secchi, 1866),
we have advanced our understanding to explain the emission as a result of a
decretion-disc which is being ionised by the central star (Struve, 1931), it
is still not clear how a Be star gains its disc. Observations conclusively
show that Be stars rotate significantly faster than their B counterparts
(Struve, 1931; Porter, 1996; Huang et al., 2010; Zorec et al., 2016), such
that potentially, the centripetal force matches the gravitational force at the
equator (Collins & Truax, 1995; Townsend et al., 2004). However the
fundamental origin of this fast rotation is still unknown, although single and
binary star channels have been proposed.

One way to achieve such rotation is for a star to be spun up by mass-transfer
in a binary system (Kriz & Harmanec, 1975; Pols et al., 1991; Liu et al.,
2006; Langer, 2012). When a star accretes material, it also accretes angular
momentum, which in the absence of tidal forces can lead to critical rotation
of the accreting star, allowing material to become unbound and form a disc.
The accretion of angular momentum is an efficient process, with a star needing
to accrete typically a few percent of its own mass to rotate critically
(Packet, 1981).

In wide systems which initiate mass-transfer after the primary has exhausted
hydrogen in the core (so-called Case B mass-transfer), tidal forces are
generally weak for the accretor star and it can be spun up to near-critical
velocities. Furthermore, rapid rotators can also originate from close systems
which undergo mass-transfer while the donor is still core hydrogen burning
(so-called Case A mass-transfer). Although tides inhibit the spin-up of the
accretor during the initial mass-transfer phases, these phases cause a
widening of the binary (Petrovic et al., 2005) so that when the donor expands
to become a giant star (initiating Case AB mass-transfer), many systems are
wide enough to render tides ineffective, allowing the mass-gainer to rotate
super-synchronously (Sen et al., in prep.). Therefore rapidly rotating mass-
gainers can originate from both short and long period systems. What is common
between these cases is that the spun-up star is usually produced after the
initially more massive star in the system has exhausted its supply of core
hydrogen.

In Be star producing binary systems where the primary star is not massive
enough to undergo a supernova explosion, a short-lived helium star or long-
lived white-dwarf would be the companion to the Be star. Despite the
difficulty of detection, both of these types of systems have been observed (Li
et al., 2012; Schootemeijer et al., 2018; Shenar et al., 2020; Coe et al.,
2020). Furthermore, studies of Be star discs have found that many are
truncated, suggesting that they are acted upon by unseen companions (Klement
et al., 2017, 2019).

When the mass-donor does explode as a supernova, the majority of systems are
expected to become unbound (Brandt & Podsiadlowski, 1995) and the Be star will
probably have no companion. The fact that this does not occur in every case is
evidenced by large numbers of Be-Xray binaries (Raguzova & Popov, 2005), which
consist of a neutron star in an eccentric orbit around a Be star such that
Xrays are produced when the Be disc and neutron star interact. When the binary
is disrupted, the Be star would likely be a runaway star. Boubert & Evans
(2018) and Dorigo Jones et al. (2020) both find the peculiar space-velocities
of Be stars in the Gaia catalogue to be consistent with a binary origin of Be
stars.

Observations show the Be phenomenon to be more common at lower metallicities
(Maeder et al., 1999; Martayan et al., 2010; Iqbal & Keller, 2013), in good
agreement with predictions of single star models whereby metal-rich stars
suffer stronger angular momentum losses through winds, thus making fast-
rotators rarer at higher metallicities (Hastings et al., 2020). Naively, this
trend is difficult to explain in the binary framework. However, further
observational characteristics of Be stars have been uncovered that are
difficult to explain with a single star formation channel. Initially fast
rotating single stars are expected to exhibit enhanced surface nitrogen
abundances, as rotational mixing dredges up CNO processed material to the
photosphere. However, there appears to be an incompatibility between models of
rotating single stars and measurements of nitrogen abundances in Be stars,
with many Be stars showing much lower nitrogen abundances than expected
(Lennon et al., 2005; Dunstall et al., 2011; Ahmed & Sigut, 2017; Hastings et
al., 2020). On the other hand, spun-up mass-gainers might not be rich in
surface nitrogen. Although the physics governing the details of mass-transfer
remains uncertain, accretion may be limited by the angular momentum content of
the gainer, such that accretion becomes non-conservative once critical
rotation is achieved (Wang et al., 2020; Langer et al., 2020). Another factor
is the strong mean molecular weight barrier established from hydrogen burning
which prevents efficient rotational mixing in the critically rotating mass-
gainer (Kippenhahn, 1974; Pinsonneault et al., 1989).

As demonstrated by Ekström et al. (2008) and Hastings et al. (2020), single
stars may achieve near-critical rotation during the late stages of hydrogen
burning, in contrast to observations showing that Be stars have a range of
fractional main-sequence ages (Zorec et al., 2005; McSwain & Gies, 2005;
Milone et al., 2018). If Be stars are mostly single, one should expect pre-
interaction binaries to host Be primaries, as whatever proposed single star
mechanism causes the Be phenomenon should work for stars in a pre-interaction
binary just as well as for single stars. It is thus telling that almost no Be
stars with a main-sequence companion have been detected (Bodensteiner et al.,
2020b).

Despite the numerous pieces of evidence to support the dominance of a binary
formation channel, several uncertainties in binary evolution prevent a solid
and accurate theoretical description of Be star populations. Proof of the
difficulty in modelling the production of Be stars is given by the contrasting
results of previous authors. It has been concluded that binaries are
responsible for either all (Shao & Li, 2014), half (Pols et al., 1991) or only
a small minority (van Bever & Vanbeveren, 1997) of galactic Be stars. This
difference is mostly due to different assumptions on mass-transfer efficiency
and the stability of mass-transfer.

In light of these uncertainties, we find it useful to determine a model free
upper-limit to Be star production from mass-transfer in binary systems.
Assessment of this limit can provide insight into whether it is at all
possible for Be stars to be formed exclusively in binaries, and to what extent
other formation mechanisms must be invoked. Under the assumption that binary
evolution dominates the production of Be stars, we can also probe uncertain
binary physics. We shall use recent high quality observations of Be stars in
open clusters (Milone et al., 2018) to give a stringent test to our simple
picture.

In Section 2 we explain our procedure for calculating an upper-limit to Be
star production from mass-transfer in binary systems, with the results of this
endeavour presented in Section 3. In Section 4 we compare our results to the
numbers of Be stars observed in young open clusters. We infer the conditions
for stable mass-transfer that are required for our prescription to reproduce
the Be fractions along the main-sequences of young open clusters in Section 5.
Uncertainties and the implications of the upper-limit are discussed in Section
6. Concluding remarks are given in Section 7.

##  2 Method

###  2.1 A hypothetical population of interacting binary stars

In order to calculate an upper-limit to the numbers of Be stars that may be
produced, we take extreme assumptions. The first of which is that the initial
binary fraction in the population is 1; that is every star is born as a member
of a binary. Nextly, as the hydrogen-burning episode of a massive star makes
up around 90% of the star’s total lifetime, we shall assume that as soon as a
primary star leaves the main-sequence, stable mass-transfer will occur on a
very short timescale, instantly producing a Be star. In our model, a Be star
shall be produced regardless of the initial period, primary mass, or mass-
ratio of the system, so that every secondary star will at some point during
its lifetime become a Be star. In this framework, the orbital period
distribution becomes irrelevant. Furthermore we shall assume that once a Be
star is formed, it remains so for the rest of its lifetime.

For simplicity we ignore the effects of mass-loss through stellar winds, such
that every system remains at its initial mass-ratio, $q$, until mass-transfer
occurs (which may be either conservative or non-conservative). Also, given the
fact that the stellar mass-luminosity relation is very steep, we define each
binary system by its most luminous component, so that each binary can be
assigned an equivalent single-star mass. To facilitate comparison with open
cluster observations, our synthetic population is assumed to be coeval.

Other properties of our population are not designed to maximise the efficiency
of Be star formation, and are more or less standard in binary evolution
calculations. We denote the initial masses of the initially more massive star
as $M_{1}$, the initially less massive as $M_{2,i}$ and define the initial
mass-ratio, $q$, as

$\displaystyle q=\frac{M_{2,i}}{M_{1}},$ (1)

such that

$\displaystyle 0<q\leq 1.$ (2)

We consider a population of binary stars where the distribution of initial
primary mass follows a power law like

$\displaystyle\xi(M_{1})=\xi_{0}M_{1}^{\alpha},$ (3)

and the distribution of initial mass-ratios is described similarly as

$\displaystyle f(q)=f_{0}q^{\kappa},$ (4)

where $\xi_{0}$ and $f_{0}$ are normalising constants to ensure that the
integral over the whole parameter space is unity (as befitting a probability-
density function). For example, the value of $f_{0}$ is easily computed as

$\displaystyle f_{0}=\frac{\kappa+1}{1-q_{min}^{\kappa+1}},$ (5)

where $q_{min}$ is the minimum mass-ratio in our population and will be
nominally set to $q_{min}=0.1$ to match the observing campaigns of Sana et al.
(2012, 2013). It is assumed that systems born with mass-ratios smaller than
this value are likely to be unstable and merge either during their formation
or early in their evolution and hence are not considered.

Mass gain of the accretor shall be parameterised by assuming that a total mass
of $\Delta M$ is accreted, giving the relation between final and initial
masses of the accretor as

$\displaystyle M_{2,f}=M_{2,i}(1+\Delta M/M_{2,i}),$ (6)

with $\Delta M/M_{2,i}$ being a free parameter.

Our assumptions on the population are summarised in the list below.

  1. 1.

initial binary fraction is 1

  2. 2.

every system will undergo stable Case B mass-transfer and form a Be star,
irrespective of period or mass ratio

  3. 3.

once a primary star leaves the main-sequence a Be star is immediately formed

  4. 4.

once a Be star is formed, it remains so for the rest of its lifetime

  5. 5.

the accretor star gains mass $\Delta M$, with the relative mass-gain $\Delta
M/M_{2,i}$ being a free parameter.

  6. 6.

the effects of wind mass-loss are ignored so that a system remains at its
initial mass-ratio until mass-transfer occurs.

  7. 7.

the distribution of initial primary masses follows a power law,
$\xi(M_{1})\propto M_{1}^{\alpha}$

  8. 8.

the distribution of initial mass-ratios follows a power law, $f(q)\propto
q^{\kappa}$

  9. 9.

only considered are binaries with a mass-ratio greater than $q_{min}$, which
is set to 0.1.

  10. 10.

when both stars are hydrogen-burning, the luminosity of a binary system is
given by that of the primary. When the primary has evolved off the main-
sequence, the luminosity of the system is naturally that of the secondary.

According to our assumptions, every secondary star with a post-main-sequence
companion is a Be star, meaning the number of Be stars with given mass $M$ is

$\displaystyle n(\textrm{Be})=n(M_{2,f}=M\,\&\,M_{1}>M_{TO}),$ (7)

where $M_{TO}$ is the turn-off mass of our coeval population. In our model the
number of non-Be stars is given by the number of primaries at a given mass.
The Be fraction, $\phi_{\textrm{Be}}(M)$ shall be defined as the number
fraction of Be stars to all stars at a given mass. In a coeval population this
becomes

$\displaystyle\phi_{\textrm{Be}}(M)$
$\displaystyle=\frac{n(M_{2,f}=M\,\&\,M_{1}>M_{TO})}{n(M_{2,f}=M\,\&\,M_{1}>M_{TO})+n(M_{1}=M)}$
(8)
$\displaystyle=\left[1+\frac{n(M_{1}=M)}{n(M_{2,f}=M\,\&\,M_{1}>M_{TO})}\right]^{-1}.$
(9)

In our model, a Be star’s mass is related to its initial mass, $M_{2,i}$, and
the relative mass-gain, $\Delta M/M_{2,i}$ via Eq. 6, such that the expression
above becomes

$\displaystyle\phi_{\textrm{Be}}(M)$
$\displaystyle=\left[1+\frac{n(M_{1}=M)}{n(M_{2,i}=\frac{M}{1+\Delta
M/M_{2,i}}\,\&\,M_{1}>M_{TO})}\right]^{-1}.$ (10)

With the aid of Eqs. 1 and 6, the condition

$\displaystyle M_{1}>M_{TO}$ (11)

can be rewritten as

$\displaystyle q$ $\displaystyle<\frac{M_{2,i}}{M_{TO}},$ (12)

leading to

$\displaystyle q$ $\displaystyle<\frac{M_{2,f}}{M_{TO}(1+\Delta M/M_{2,i})}.$
(13)

This results in

$\displaystyle\phi_{\textrm{Be}}(M)$
$\displaystyle=\left[1+\frac{n(M_{1}=M)}{n(M_{2,i}=\frac{M}{1+\Delta
M/M_{2,i}}\,\&\,q<\frac{M}{M_{TO}(1+\Delta M/M_{2,i})})}\right]^{-1}.$ (14)

To study coeval populations, a more convenient approach is to find the Be
fraction as a function of the fractional main-sequence turn-off mass,
$M/M_{TO}$. This produces the expression

$\displaystyle\phi_{\textrm{Be}}(M/M_{TO})$
$\displaystyle=\left[1+\frac{n(M_{1}=M/M_{TO})}{n(M_{2,i}=\frac{M}{M_{TO}(1+\Delta
M/M_{2,i})}\,\&\,q<\frac{M}{M_{TO}(1+\Delta M/M_{2,i})})}\right]^{-1},$ (15)

which shall be our basis for exploring the Be fraction in coeval populations.

To evaluate the Be fraction, it is necessary to find the relative numbers of
primary stars to that of secondary stars at a given mass. We may write that
the number of primary stars with a given mass, $n(M_{1})$, is the integral of
the primary mass distribution across an infinitesimally small mass range,
$dM_{1}$, multiplied by the total number of stars in the population,
$n_{tot}$, as

$\displaystyle
n(M_{1})=n_{tot}\xi({M_{1}})dM_{1}=n_{tot}\xi_{0}M_{1}^{\alpha}dM_{1}.$ (16)

To tackle the number of secondary stars at a given mass is sightly more
involved as we do not have directly the distribution of secondary masses,
instead it is inferred from the primary mass and mass-ratio distributions.
First consider a population in which there exists only a single mass-ratio,
$q_{0}$, ie. the mass-ratio distribution is a delta-Dirac function. If one is
interested in the number of secondary stars with initial mass $M_{2,i}$, one
must count the number of primaries with mass $M_{2,i}/q_{0}$, so we have

$\displaystyle
n(M_{2,i}\,\&\,q_{0})=n_{tot}\xi\left(\frac{M_{2,i}}{q_{0}}\right)d\left(\frac{M_{2,i}}{q_{0}}\right),$
(17)

with $d\left(\frac{M_{2,i}}{q_{0}}\right)$ representing an infinitesimally
small change in $\frac{M_{2,i}}{q_{0}}$.

Any distribution may be expressed as an infinite sum of appropriately weighted
delta-Dirac distributions, with the weighting coming from the probability-
density function. Therefore for the general case we have

$\displaystyle
n(M_{2,i})=n_{tot}\int_{q_{min}}^{1}f(q)\xi\left(\frac{M_{2,i}}{q}\right)d\left(\frac{M_{2,i}}{q}\right)dq.$
(18)

It is then clear that the limits of the integral above place constraints on
the initial mass-ratios counted. The number of systems with a given initial
secondary mass $M_{2,i}$ and initial mass-ratios between $q_{min}$ and
$q_{max}$ can thus be written as

$\displaystyle
n(M_{2,i}\,\&\,q_{min}<q<q_{max})=n_{tot}\int_{q_{min}}^{q_{max}}f(q)\xi\left(\frac{M_{2,i}}{q}\right)d\left(\frac{M_{2,i}}{q}\right)dq,$
(19)

The differential $d\left(\frac{M_{2,i}}{q}\right)$ in Eq.19 is quite
cumbersome so we chose to let

$\displaystyle r=\frac{M_{2,i}}{q}.$ (20)

We may now write

$\displaystyle
n(M_{2,i}\,\&\,q_{min}<q<q_{max})=n_{tot}\int_{q_{min}}^{q_{max}}f(q)\xi(r)drdq.$
(21)

We have

$\displaystyle M_{2,i}=qM_{1},$ (22)

thus

$\displaystyle dM_{2,i}=qdM_{1}+M_{1}dq.$ (23)

As we are interested in the number of secondary stars at a fixed mass,
$d{M_{2,i}}=0$, so

$\displaystyle dq=-\frac{q}{M_{1}}dM_{1}.$ (24)

Differentiating $r$ gives

$\displaystyle dr=\frac{-M_{2,i}}{q^{2}}dq.$ (25)

Combining Eqs. 24 and 25 results in

$\displaystyle dr=\frac{1}{q}\frac{M_{2,i}}{M_{1}}dM_{1}.$ (26)

Inserting this into our expression for $n(M_{2,i}\,\&\,q_{min}<q<q_{max})$
(Eq. 21) gives

$\displaystyle n(M_{2,i}\,\&\,q_{min}<q<q_{max})$
$\displaystyle=n_{tot}\int_{q_{min}}^{q_{max}}f(q)\xi\left(\frac{M_{2,i}}{q}\right)\frac{1}{q}\frac{M_{2,i}}{M_{1}}dM_{1}dq.$
(27)

We now divide Eq. 16 by Eq. 27 leaving

$\displaystyle\frac{n(M_{1})}{n(M_{2,i}\,\&\,q_{min}<q<q_{max})}=\frac{M_{1}\xi(M_{1})}{M_{2,i}\int_{q_{min}}^{q_{max}}f(q)\xi\left(\frac{M_{2,i}}{q}\right)\frac{1}{q}dq}.$
(28)

When the distributions for initial primary mass and mass-ratio, Eqs. 3 and 4,
are inserted, Eq. 28 simplifies further to

$\displaystyle\frac{n(M_{1})}{n(M_{2,i}\,\&\,q_{min}<q<q_{max})}=\left(\frac{M_{1}}{M_{2,i}}\right)^{\alpha+1}\frac{1}{\int_{q_{min}}^{q_{max}}f_{0}q^{\kappa-\alpha-1}dq}.$
(29)

This result may be readily checked against Monte-Carlo sampling of the primary
mass and mass-ratio distributions.

Equation 29 can be used to directly determine the Be fraction (Eq. 15) by
setting $M_{1}=M/M_{TO}$ , $M_{2,i}=\frac{M}{M_{TO}(1+\Delta M/M_{2,i})}$ and
$q_{max}=\frac{M}{M_{TO}(1+\Delta M/M_{2,i})}$. This leaves

$\displaystyle\phi_{Be}(M/M_{TO})=\Bigg{[}1+\frac{(1+\Delta
M/M_{2,i})^{(\alpha+1)}}{\int_{q_{min}}^{\frac{M}{(1+\Delta
M/M_{2,i})M_{TO}}}f_{0}q^{\kappa-\alpha-1}dq}\bigg{]}^{-1},$ (30)

where the integral has a simple analytic solution. Eq. 30 describes the Be
star fraction as a function of the fractional turn-off mass, $M/M_{TO}$ for
our model open cluster. As the mass-dependence in Eq. 30 is expressed by the
fractional turn-off mass, is not necessary to specify the turn-off mass.

###  2.2 Limits for remaining parameters

All that remains is to explore Eq. 30 in a suitable parameter space. The
parameters we have are $\alpha$, $\kappa$ and $\Delta M/M_{2,i}$, the primary
mass distribution exponent, initial mass-ratio distribution exponent and
relative accretor mass-gain respectively.

The canonical value for the initial-mass-function (IMF) exponent, $\alpha$, is
given by the Salpeter IMF, $\alpha=-2.35$ (Salpeter, 1955). However recent
observations of young stars in the 30 Doradus starburst region suggest instead
$\alpha=-1.90^{+0.37}_{-0.26}$ (Schneider et al., 2018). Similarly in the R136
star forming region an exponent of $\alpha=-2.0{\pm 0.3}$ was found
(Bestenlehner et al., 2020). On the other hand it has also been proposed that
the IMF follows an even steeper law with $\alpha=-2.7$ (Scalo, 1986).
Therefore we consider the range $-1.9<\alpha<-2.7$ .

Observations of Galactic O-type stars show that the mass-ratio distribution
follows a power law with exponent $\kappa=-0.1\pm 0.6$ (Sana et al., 2012) for
$0.1<q<1$. In the Large Magellanic Cloud, the mass-ratios of massive binaries
appear to be distributed differently with $\kappa=-1.0\pm 0.4$ again in the
range $0.1<q<1$ (Sana et al., 2013). There are many claims that mass-ratios of
binaries favour either low values (Trimble, 1990; Tout, 1991; Hogeveen, 1991)
or follow a uniform distribution(Kobulnicky & Fryer, 2007; Kouwenhoven et al.,
2007). In light of these findings we shall consider $\kappa$ values in the
range $-1<\kappa<0$.

Estimates of the accretor mass-gain, $\Delta M/M_{2,i}$, obtained by demanding
that mass-transfer stops once the mass-gainer reaches critical rotation tell
us that $\Delta M/M_{2,i}$ is at the very most 0.1 and in most cases around
0.02, depending on the angular momentum content and physical structure of the
mass-gainer before accretion (Packet, 1981; Petrovic et al., 2005; Wang et
al., 2020). It has been found that around 70% of the mass leaving the donor
must be ejected from the system to explain observed distributions of Be star
masses in Be X-ray binaries (Vinciguerra et al., 2020). However it must be
noted that because it is expected that up to 90% of massive binary systems are
broken apart by a supernova kick (Brandt & Podsiadlowski, 1995), Be X-ray
binaries represent a small fraction of the population, and hence may well
contain strong biases. Furthermore, it is believed that mass-transfer must be
highly non-conservative to explain observed populations of Wolf-Rayet O-star
binaries (Petrovic et al., 2005; Shao & Li, 2016). On the other hand, several
systems exist which show evidence of near-conservative mass-transfer having
taken place (de Mink et al., 2007; Schootemeijer et al., 2018; Brož et al.,
2021). To fully explore the effects of mass-transfer efficiency on Be star
populations, we take the full range $0<\Delta M/M_{2,i}<1$.

##  3 Results

The results of Eq. 30 are plotted in Fig. 1 for the extremal parameters
outlined in Sect. 2.2. The primary mass distribution affects the absolute
numbers of Be stars because for a shallower distribution ($\alpha$ closer to
0), there is an abundance of massive binaries, such that many systems contain
post main-sequence primaries and therefore the number of Be stars increases.
Conversely when $\alpha<<0$, the population contains fewer primaries of mass
greater than the turn-off mass and the Be count decreases.

The effect of the mass-ratio distribution can be understood by considering a
population with a high value of $\kappa$ such that secondary stars have a
similar mass to their companion. In this case, when the primary leaves the
main-sequence, the secondary will be rather evolved, and hence most Be stars
will be found near the turn-off. On the other hand in a population with a low
$\kappa$, the opposite is true; the secondary stars will have low masses
compared to the turn-off mass and Be stars will be more evenly distributed
along the main-sequence, as seen in Fig. 1.

Figure 1 shows how a varying mass-gain changes the Be count, with accretors
that gain more mass producing fewer Be stars. This can be understood by
considering a Be star of mass 0.9$M_{TO}$ that is produced by inefficient
mass-transfer. For the primary to exceed the turn-off mass, the initial mass-
ratio of the system must be less than 0.9. Now if this star had gained
0.1$M_{TO}$, the initial mass would be 0.8$M_{TO}$ and the initial mass-ratio
must be less than 0.8. Therefore, mass gain restricts the number of systems
that are able to produce Be stars of a given mass, and low mass-transfer
efficiency leads to higher numbers of Be stars being produced at a given mass.

It is seen in Fig. 1 that the largest Be fractions are produced when mass-
transfer is inefficient ($\Delta M/M_{2,i}=0$) and the initial mass function
is shallow. The mass-ratio distribution then tunes the distribution of Be
stars along the main-sequence. Therefore it is judged that the most Be stars
are produced with the parameters $\Delta M/M_{2,i}=0$, $\alpha=-1.9$ and
$-1<\kappa<0$. Depending on the chosen parameters, the maximum Be fraction is
in the range 0.2-0.35 near the main-sequence turn-off.

Figure 1: The maximum Be fraction, $\Phi_{Be}$, in a coeval population as
defined by Eq. 30 plotted as a function of fractional main-sequence turn-off
mass, $M/M_{TO}$ for varying parameters. The left and right Panels show
$\Delta M/M_{2,i}=$ 0, and 1 respectively. The colour of the lines represents
differing $\kappa$ values with red being $\kappa=0$ and blue $\kappa=-1$.
Dashed lines show $\alpha=-1.9$ and solid lines $\alpha=-2.7$ as indicated by
the annotations.

##  4 Comparison to observations

To contextualise to our results, we attempt here a comparison with
observations using high quality Hubble Space Telescope photometry of young
Small and Large Magellanic Cloud open clusters in which Be stars are revealed
as bright objects in a narrow-band filter centred on H$\alpha$ (Milone et al.,
2018). Photometry was performed with Hubble wide-band filters $F814W$ and
$F336W$ and the narrow band $F656N$ filter, allowing one to produce colour-
magnitude diagrams in which Be stars are identified from H$\alpha$ photometry.

As many spectroscopically confirmed Be stars in NGC 330 are bright in
H$\alpha$ (Bodensteiner et al., 2020a), we judge H$\alpha$ emission to be a
good proxy for Be stars. It is possible for the accretion discs of Algol-type
binaries to exhibit H$\alpha$ emission (Peters, 1989) however, such systems
are expected to contribute around 3% to the total population (de Mink et al.
2014, Sen et al. in prep.). Furthermore, some field stars may be H$\alpha$
emitters, with Milone et al. (2018) noting that no more than one-tenth of
stars in the cluster field are suspected field stars. Field stars will also
contaminate the population of stars not emitting in H$\alpha$, therefore their
presence is not expected to significantly alter the relative fractions of
H$\alpha$ emitters and non-emitters.

Be star fractions have previously been measured as a function of magnitude
(Keller et al., 1999; Milone et al., 2018; Bodensteiner et al., 2020a).
However we find it worthwhile to repeat this exercise, including several
factors which were previously overlooked.

###  4.1 Counting Be stars

Our goal is to measure the observed Be fraction as a function of mass along
the main-sequence of an open cluster. To do this we must note the two major
differences between a Be star and a ”normal” B star; fast rotation and the
presence of a decretion disc. Due to the effect of the centrifugal force, a
fast-rotating star suffers from reduced effective gravity at the equator, and
according to the von Zeipel theorem (von Zeipel, 1924), this results in a
lower effective temperature. Therefore fast-rotating stars are cooler and
redder than their non-rotating counterparts. Furthermore, light from a Be star
consists of radiation from the star itself and also light from the decretion
disc. Typically the average temperature of the disc is around 70% that of the
star’s effective temperature (Sigut et al., 2009), and so the disc is expected
to emit mostly in visible and infra-red wavelengths.

In a colour-magnitude diagram, the magnitude in a red filter is plotted on the
y-axis and a colour defined by the blue and red filter (B - R) on the x-axis.
When a star becomes brighter in the red filter, it will therefore move to the
right and upwards in the colour-magnitude diagram. This effect means that to
count the Be stars as a function of mass, we must do so in bins that are
sloped with respect to the x-axis. The gradient of this slope depends on how
much redder a near-critically rotating star is than a slow rotator at the same
mass, and on how much light the decretion disc radiates.

As no reliable numerical models exist of stars rotating at the critical
velocity, we shall adopt a simple model to relate the luminosity and
temperature of a critical rotator to an equivalent non-rotating star. After
having been spun up, a star will change its shape, becoming oblate. At the
same time we do not expect a great difference in luminosity between a star
before and after the spin-up. This is because stars are generally very
centrally condensed, such that the centrifugal force is small compared to
gravity in the regions where nuclear burning occurs, meaning that (excluding
the effects of rotational mixing) central temperatures and thus luminosities
are not very sensitive to rotation, in agreement with models (Brott et al.,
2011; Paxton et al., 2019). Following the Steffan-Boltzmann law, because of
the increased surface area of a critical rotator, the effective temperature
decreases. Using the Roche model (see Appendix A) , one can show that a
critically rotating star has a surface area of approximately $1.58\times 4\pi
R_{p}^{2}$, with $R_{p}$ being the polar radius, which corresponds to a
decrease in effective temperature by a factor of ${1.58^{-\frac{1}{4}}}\approx
0.89$. Knowing this, we can construct isochrones describing the intrinsic
properties of critically rotating stars from non-rotating isochrones.

A further complication that is brought about by gravity darkening is that a
fast rotating star appears cooler and dimmer when viewed equator-on as
compared to pole-on. Assuming a random orientation of the inclination axis,
the mean value of the sine of the inclination angle is $\pi/4$, corresponding
to a mean inclination angle of 51.8 $\mathrm{\SIUnitSymbolDegree}$. To take
into account the mean effect of gravity darkening, we employ the model of
Espinosa Lara & Rieutord (2011) as implemented in MESA (Paxton et al., 2019).
Here, the projected luminosity and effective temperature, $L_{\textrm{proj}}$,
$T_{\textrm{eff,proj}}$ are related to the intrinsic luminosity and effective
temperature, $L$, $T_{\textrm{eff}}$ by

$\displaystyle L_{\textrm{proj}}$ $\displaystyle=C_{T}(\omega,i)L$ (31)
$\displaystyle T_{\textrm{eff,proj}}$
$\displaystyle=C_{L}(\omega,i)T_{\textrm{eff}},$ (32)

with $C_{T}$ and $C_{L}$ depending on the fraction of critical velocity,
$\omega$ and inclination angle $i$.

The temperatures and luminosities of critically rotating stars are found by
using the coefficients
$C_{T}(\omega=1,i=51.8$\mathrm{\SIUnitSymbolDegree}$)=1.02$ and
$C_{L}(\omega=1,i=51.8$\mathrm{\SIUnitSymbolDegree}$)=1.22$. It is a rather
curious feature of the gravity darkening model that at the mean inclination,
the coefficients exceed unity, meaning the average effect of gravity darkening
is not darkening at all, but brightening. Finally, by interpolating tables of
synthetic stellar spectra (Choi et al., 2016) to produce magnitudes in Hubble
filters, we are able to produce an isochrone of critical rotators, as shown in
Fig. 2.

The contribution of a Be star’s disc to its total flux is more difficult to
assess. It has been noticed that a loss of spectral emission features in
certain Be stars coincides with a dimming of around 0.3-0.5 magnitudes in the
R and V filters (Carciofi et al., 2012; Labadie-Bartz et al., 2017; Rímulo et
al., 2018). If the loss of emission features is interpreted as the
disappearance of the disc, one can take this change in brightness to equal the
flux contribution of the disc. By comparing the colour of our isochrones with
the colours of Be stars in NGC 330, we can assess how much the Be disc shines,
as in Fig. 2. After assuming that the disc shines in the F814W filter but not
in the F336W filter, we find a reasonable fit to the H$\alpha$ emitters when a
disc brightness of 0.25${m}_{\textrm{F814W}}\thinspace$is adopted, as shown by
the solid and dashed purple lines in Fig. 2.

Figure 2: Colour-magnitude diagram of NGC 330 focused on the turn-off region
with H$\alpha$ emitters marked in red. An isochrone of non-rotating stars is
plotted in blue (see App. B for model details). Green and red isochrones
depict critically rotating stars viewed pole-on and equator-on respectively,
as derived from a simple model of critical rotators (see App. A). The solid
purple isochrone represents critically rotating stars viewed at the mean
inclination angle when the rotation axis is randomly oriented
(51.8$\mathrm{\SIUnitSymbolDegree}$) and for the dotted purple
0.25${m}_{\textrm{F814W}}\thinspace$has been added to simulate the decretion
disc. The isochrone age is 30 Myr, distance modulus $\mu$ = 18.8 mag and
reddening of E(B-V)=0.1 mag. Data from Milone et al. (2018)

For NGC 330, in the regions of the colour-magnitude diagram containing Be
stars, we find stars of equal mass on the non-rotating and Be star isochrones
to be connected by lines of gradient
$\frac{-d{m}_{\textrm{F814W}}}{d(m_{\textrm{F336W}}-{m}_{\textrm{F814W}})}=$
2\. For NGC 2164, again assuming a constant disc magnitude of
0.25${m}_{\textrm{F814W}}\thinspace$, the gradient is found to be 1.8. These
differing values are caused by the ways in which stellar spectra, and hence
magnitude in a given filter, vary with luminosity and effective temperature.

Figure 3 shows the colour-magnitude diagrams of NGC 330 and NGC 2164 with the
Be fraction as counted in slanted bins with gradients of 2.0 and 1.8
respectively. It is noted that as compared to counting the Be fraction in bins
of constant ${m}_{\textrm{F814W}}\thinspace$magnitude (ie. horizontal bins)
the values measured here are lower because in the horizontal bin case, one is
counting B stars with a higher mass than the Be stars in the same bin.
According the initial-mass-function, the higher mass stars are less populous
and hence the Be fraction increases solely because there are fewer B stars
being counted.

We use isochrones of rotating single stars based on an extended model grid of
Schootemeijer et al. (2019) (see Appendix B for a thorough description) with
an initial rotation rate of
${v_{\textrm{rot}}}/{v_{\textrm{crit}}}\thinspace$= 0.6 to assign mass-ranges
to each bin, so that the Be fraction can be evaluated as a function of mass.
The bins are placed so that the outer edge of the last bin is at the point
where hydrogen has been exhausted in the stellar core. The value of
${v_{\textrm{rot}}}/{v_{\textrm{crit}}}\thinspace$= 0.6 is chosen as suggested
by Gossage et al. (2019) and Wang et al. (in prep.), and produces equatorial
rotation velocities that are in broad agreement with spectroscopic
observations (Dufton et al., 2013; Marino et al., 2018; Sun et al., 2019;
Kamann et al., 2020). The reddening and distance modulus values are tailored
to give the best fit to the cluster and are in good agreement with previous
isochrone fittings for these clusters (Milone et al., 2018). The isochrone
fits are shown in Fig 3.

Figure 3: Colour-magnitude diagrams with isochrone fits and Be star counts for
Small Magellanic Cloud cluster NGC 330 (left) and Large Magellanic Cloud
cluster NGC 2164 (right). H$\alpha$ emitters are marked in red. Bottom panels
show the Be fraction counted in bins as defined in the top panels with the
errors given by the binomial counting error. The bins have a gradient of 2.0
and 1.8 for NGC 330 and NGC 2164 respectively. Mass values associated with the
bins are provided by the isochrone fit. For both clusters the isochrone
depicts stars with initial rotation equal to 0.6
${v_{\textrm{rot}}}/{v_{\textrm{crit}}}\thinspace$. For NGC 330 the isochrone
age is 30 Myr, distance modulus $\mu$ = 18.8 mag and reddening of E(B-V)=0.1
mag. For NGC 2164 the age is 80 Myr, $\mu$= 18.3 mag and E(B-V)=0.12 mag. Data
from Milone et al. (2018)

The isochrones allow us to measure the turn-off mass and the masses associated
with each bin, thereby a direct comparison between the theory presented in
Sec. 3 and observations is possible. Figure 4 shows this comparison, with
counting uncertainties on the Be fraction given by the standard error,
$\sigma$, assuming a binomial distribution as

$\displaystyle\sigma=\sqrt{\Phi(1-\Phi)/N},$ (33)

with $\Phi$ being the measured Be fraction and $N$ the total number of stars
in a given bin.

We find that despite the two clusters being of different metallicities and
ages, they seem to have similar Be fractions as a function of relative turn-
off mass. This may be an indication that whatever the dominant Be production
channel is, it is universal.

In both clusters the Be fraction steadily increases from zero to around 0.4 in
the range 60-80% of the turn-off mass. Near the turn-off, the Be fraction is
found to be approximately 0.4 with significant counting uncertainty due to the
relatively small numbers of stars occupying this region. Taking into account
these uncertainties, it is seen that our upper limit can describe the numbers
of Be stars in the upper part of the main-sequence. It is important to note
that because of the difficulty in performing an isochrone fit, the Be
fractions near the turn-off are particularly uncertain, with a small change in
the isochrone fit resulting in a large change in the measured Be fraction (see
Sec. 6.1 for a quantitative discussion). Therefore, despite the measured Be
fraction in NGC 330 at times exceeding our upper limit, it is reasonable to
conclude that the upper limit does provide a reasonable fit to the Be star
numbers near the turn-off. However it does fail to explain the lack of Be
stars below $M/M_{TO}\approx 0.7$. This may be the result of certain systems
not forming Be stars but instead merging, as shall be discussed in the next
section.

Figure 4: Comparison between theory and observations. Be fraction as a
function of fractional main-sequence turn-off mass in NGC 330 and NGC 2164 as
shown in Fig. 3. Dashed lines show theoretical upper limit given by Eq. 30
with $\alpha=-1.9$, $\Delta M/M_{2,i}=0$ and $\kappa=-1.0,0$ (see Fig. 1) as
given by legend.

##  5 Inferring the initial conditions for stable mass-transfer

The observations presented in Sec. 4 show that our upper limit can
approximately describe the numbers of Be stars near the turn off, but fails to
reproduce the Be sequence’s sharp cut-off. Here, we shall investigate how our
prescription will need to be changed in order to reproduce this feature.

In reality, not every binary system will undergo stable mass-transfer to form
a Be star. For the specific case of the donor being in the Hertzsprung gap, as
it is in Case AB or Case B mass-transfer, the mass-transfer proceeds at the
Kelvin-Helmholtz (or thermal) timescale (Tout et al., 1997; Wellstein et al.,
2001), meaning that if there is a large discrepancy in the Kelvin-Helmholtz
timescales of the donor and accretor, the mass-transfer will become unstable
and a common envelope situation will ensue, most likely leading to a stellar
merger.

To model the occurrence of mergers, it is often assumed in simplified binary
evolution calculations (Pols et al., 1991; Hurley et al., 2002; Schneider et
al., 2015) that systems below a certain mass-ratio will merge, however such a
simple criterion is unsuitable to reproduce the observations shown in Fig. 3.
Equation 30 gives the Be star fraction as an integral quantity, such that the
Be star fraction at the main-sequence turn-off is the accumulation of systems
with mass-ratios from $q_{min}$ to 1. This may be understood intuitively by
noting that a Be star of mass near the main-sequence turn-off mass can
originate from either an extreme mass-ratio system with a very massive
primary, or from a system with mass-ratio close to unity. Therefore, when we
demand that all systems below a given mass-ratio merge, we will naturally
decrease the Be fraction at the turn-off, which we must avoid to retrieve high
numbers of Be stars at the turn-off.

To keep the Be fraction near the turn-off high and produce a sharp break in
the Be fraction at $M/M_{TO}\approx 0.7$, more sophisticated criteria are
needed, namely with dependence on primary mass and mass-ratio. We propose that
the systems most likely to suffer unstable mass-transfer are those with an
extreme mass-ratio and low primary mass, as the components of such systems
have the largest difference in Kelvin-Helmholtz timescales. This can be
visualised in a grid of primary mass against mass-ratio, with the bottom
corner consisting of systems that merge. In such a grid, systems with a fixed
secondary mass are represented by parabolae, as depicted in Fig. 5a. If the
parabola representing a secondary with the turn-off mass can avoid the region
containing merger progenitors, the Be fraction at the turn-off will remain
close to the maximum theoretical prediction. Then as the secondary mass
decreases, the parabolae will move into the corner with low mass-ratio and low
primary mass and consequently the Be fraction will decrease.

To make a test of our hypothesis, we perform a Monte-Carlo simulation, whereby
systems are picked randomly from given distributions of initial primary mass
and initial mass-ratio. As before, we shall assume that mass-transfer is
completely non-conservative ($\Delta M/M_{2,i}=0$). By choosing a turn-off
mass, we can calculate the masses of Be stars in the simulation and therefore
assess the Be fraction. The occurrence of mergers is decided using the stable
mass-transfer region depicted in Fig. 5 a. The motivation for selecting this
region will now be explained. Analysis of mass-transfer from giant donors
(Pavlovskii & Ivanova, 2015) has indicated that mass-transfer from
Hertzsprung-gap stars is stable at mass-ratios greater than around 0.6. We
will therefore assume that all systems with initial mass-ratios greater than
0.6 will undergo stable mass-transfer. The stability of mass-transfer is
determined by the donor’s reaction to mass-loss, where stars with radiative
envelopes generally tend to contract as the envelope is being stripped
(Hjellming & Webbink, 1987). This is reversed for convective-envelope stars,
which typically expand in response to mass loss (Hjellming & Webbink, 1987).
Stellar structure calculations suggest that stars with a mass greater than
around 60$\,\mathrm{M_{\odot}}$ spend very little time as red-giants, meaning
that they mostly have radiative envelopes (Schootemeijer et al., 2019; Klencki
et al., 2021) such that mass-transfer is much more likely to occur when the
donor has a radiative envelope. We therefore propose that mass-transfer will
be stable for all systems with a primary mass exceeding
60$\,\mathrm{M_{\odot}}$. The region of instability is then defined by a
linear interpolation between systems with $M_{1}=60\,\mathrm{M_{\odot}},q=0.1$
and $M_{1}=5\,\mathrm{M_{\odot}},q=0.6$, as depicted in Fig. 5a. We shall
again assume that the orbital period plays no role in determining the
stability of mass-transfer, hence allowing us to not specify an orbital period
distribution.

Merger products change the distribution of masses in a population, hence
affect the Be fraction as a function of mass. We will assume that Be stars are
not merger products for two reasons. Firstly, it is believed that, although
merger products are fast-rotators initially, while thermal equilibrium is
returned, internal angular momentum redistribution causes a rapid spin-down
(Schneider et al., 2019). What is more, stellar mergers may produce strongly
magnetised stars (Ferrario et al., 2009; Wickramasinghe et al., 2014;
Schneider et al., 2019) which would further spin-down due to magnetic braking.
Secondly, a merger between a star with a helium core and a main-sequence
object will not produce a hydrogen burning star owing to the higher mean
molecular weight and lower entropy of the helium-core star (Langer, 2012;
Justham et al., 2014). As observations (Milone et al., 2018) show Be stars to
be concentrated on the main-sequence, we assume that Be stars are unlikely to
be produced from the merging of two stars. For mergers, we assume the fraction
of mass lost during the merging process to the total binary mass to be equal
to

$\displaystyle\mu_{loss}=\frac{0.3q}{1+q^{2}},$ (34)

(Glebbeek & Pols, 2008) which equates to between 2 and 15% over the range
$0.1<q<1$. As the mass lost during the merging process is assumed to be low,
mergers will always have a mass exceeding the turn-off mass and will not
affect the Be fractions on the main-sequence.

Figure 5 b shows the results of the simulation for clusters with turn-off
masses of 9 and 6$\,\mathrm{M_{\odot}}$, which roughly correspond to NGC 330
and NGC 2164, respectively. The chosen criteria have maintained a high Be
fraction near the turn-off and also produced a sudden end to the Be-sequence,
and provide a reasonable fit to the measured Be fractions in NGC 330 and NGC
2164. It is remarkable that such simple, although physically motivated, stable
mass-transfer criteria can successfully reproduce the numbers of Be stars in
the open clusters studied. Our empirical mass-transfer stability criteria
could be tested in the next generation of detailed binary evolution models.

Figure 5: a) Adopted region of stable mass-transfer in the primary mass-mass-
ratio plane. Regions coloured red experience unstable mass-transfer and merge,
while for green regions, mass-transfer is stable and a Be star is formed. The
black and orange lines show systems with secondary masses of 9 and 6
$\,\mathrm{M_{\odot}}$ respectively. b) Results of a Monte-Carlo simulation
showing the Be fraction, $\Phi_{Be}$, when the stable mass-transfer region in
a) is applied. Binary systems have a flat mass-ratio distribution
($\kappa=0$), a primary mass distribution $\xi(M_{1})\propto M^{-1.9}$ and we
have assumed inefficient accretion ($\Delta M/M_{2,i}=0$). The black line
shows a simulation with a turn-off mass of 9$\,\mathrm{M_{\odot}}$, and the
orange line of 6$\,\mathrm{M_{\odot}}$. The dashed grey line shows the
theoretical upper limit, as given by Eq. 30. Measured Be fractions of NGC 330
and NGC 2164 according to Fig. 3 are plotted as black and orange crosses
respectively.

##  6 Discussion

###  6.1 Uncertainties

The largest uncertainty in our procedure comes from the isochrone fits. Most,
if not all open clusters display an extended main-sequence turn-off, making
the choice of a suitable isochrone age difficult. This is illustrated in Fig.
6, where isochrones of two different ages are fitted to NGC 330 and the Be
fractions are evaluated. It is seen that a small variation in the adopted age
can cause the Be count in some bins to vary by up to 0.2, with the end of the
Be sequence being particularly affected. A similar sensitivity is also found
for small differences in the distance modulus, reddening and isochrone
rotation rates. From Fig. 6, it is judged that the uncertainty on the measured
Be fractions is approximately 0.1 without including the counting error.

Figure 6: Left panel: colour-magnitude diagram of NGC 330 with H$\alpha$
emitters marked in red. Isochrone fits with ages 30 and 36Myr are plotted in
purple and green respectively. Both isochrones have initial rotation of
${v_{\textrm{rot}}}/{v_{\textrm{crit}}}\thinspace$=0.6, a distance modulus of
18.8 mag and reddening of E(B-V)=0.1 mag. Right panel: Be fraction as a
function of fractional turn-off mass as measured by the 30Myr isochrone (in
purple) and the 36Myr isochrone (in green). Dashed lines show theoretical
upper limit given by Eq. 30 with $\alpha=-1.9$, $\Delta M/M_{2,i}=0$ and
$\kappa=-1.0,0$ as given by legend.

To measure the observed Be fraction as a function of mass, we must use slanted
bins. In calculating the gradient of these bins, we have assumed that the Be
star disc always adds 0.25${m}_{\textrm{F814W}}\thinspace$to the magnitude of
the star. This may be an over-simplification, with Be stars of differing mass
or evolutionary status hosting relatively brighter or dimmer discs.
Unfortunately this effect is difficult to observe and characterise and is also
compounded by the fact that Be stars can display spectral and photometric
variability (Porter & Rivinius, 2003). Sigut et al. (2009) report that the
ratio of stellar effective temperature to mean disc temperature and infrared
excess are indeed functions of spectral type.

Far older clusters, such as the 300Myr old NGC 1856 have much lower Be
fractions than their younger counterparts (see Milone et al., 2018, Fig. 17).
Our simple model and mass-transfer stability criteria predict that the Be
fraction does not vary strongly with turn-off mass and therefore is unable to
explain the turn-off Be fraction in NGC 1856 of around 0.2. However this
discrepancy may be partly explained by a changing binary fraction with mass,
as it is known that more massive stars display a stronger preference for
binary companionship (Köhler et al., 2006; Kouwenhoven et al., 2009), thus
older clusters which contain fewer binaries will naturally have fewer Be
stars. Another aspect behind the emission line phenomenon is the ionising
power of the star, because to produce an emission line, the central Be star
must ionise its decretion disc. Without sufficient ionising power, even if a
decretion disc is present, no emission line will be observable and hence the
star will seem ordinary. The ionising photon emission rate is known to be
strongly dependant on effective temperature (Sternberg et al., 2003), so at
some limiting mass one would expect the central star to be unable to ionise a
disc. This effect may play a role in lowering the Be fraction in older
clusters and causing a dearth of Be stars at low magnitudes in the colour-
magnitude diagram. Both of the clusters studied here have several stars that
are very red, despite not being marked as H$\alpha$ emitters, which could in
fact be such ”dormant Be stars”.

Lastly we note that in our work we have assumed that the properties of binary
systems are distributed according to very simple laws. In reality however, the
distributions may well be complex functions of one another, for example the
mass-ratio distribution might be a function of the primary mass. The nature of
these distributions is set by poorly understood binary star formation
mechanisms as outlined by Tokovinin & Moe (2020).

###  6.2 Mass-transfer efficiency

To construct our prediction of Be star fractions, we have assumed that no mass
is accreted during mass-transfer and seen that this scenario fits observations
reasonably well. It has been demonstrated that for efficient mass-transfer,
the Be star fraction decreases, meaning that the theoretical framework
presented here would not fit observations if mass-transfer were efficient.
This leads us to propose that, if the binary Be formation channel is the
dominant one, mass-transfer is, on average, far from conservative.

Binary models with conservative mass-transfer predict Be stars to be blue
stragglers, after having gained a lot of mass (van Bever & Vanbeveren, 1997).
The observations presented in Fig. 3 contradict this prediction, with the vast
majority of Be stars lying either on the main-sequence or being slightly
redder than it, strengthening our conclusion about mass-transfer being highly
non-conservative.

###  6.3 On the initial binary fraction

In obtaining our results we have assumed an initial binary fraction of 1,
which could be criticised as too extreme. We have demonstrated that in a
coeval population of binary systems, at most 30% of systems are post-
interaction binaries (see Sec. 3). Pre-interaction systems would therefore
make up no less than 70% of this population. Dedicated models show that under
the assumption of a constant star formation rate, $30^{+10}_{-15}$ % of
massive stars are the products of binary interaction (de Mink et al., 2014),
in broad agreement with this work.

Post-interaction binaries are either merger products, contain a relatively low
mass post-main-sequence object (helium star, black-hole, neutron star or
white-dwarf) with a main-sequence (possibly emission line) star or form a
runaway star ejected from the binary orbit after a supernova. Such objects
manifest themselves either as single stars, or would be difficult to detect as
binaries (de Mink et al., 2011b, a). Thus, even in a population whose initial
binary fraction is 1, apparently single stars are present in the proportions
described above.

By examining radial velocity variations of very massive stars, one may only
measure the pre-interaction binary fraction (as supernova kicks are believed
to disrupt almost all binary systems (Brandt & Podsiadlowski, 1995)), which
has been observed to be around 0.7 for O-type stars (Sana et al., 2012). We
therefore argue that the initial binary fraction is certainly greater than the
observed pre-interaction binary fraction, such that at this stage, we must
remain open to the possibility that an initial binary fraction very close to
one is indeed realised in nature.

##  7 Conclusions

In light of various uncertainties plaguing binary evolution calculations, we
have investigated whether binary evolution can possibly reflect the large
numbers of Be stars observed in open clusters. Starting from the premise that
any binary system, regardless of primary mass, orbital period or mass-ratio,
will undergo stable mass-transfer to form a Be star, we have calculated a
rigorous upper limit to the formation of Be stars through this channel. It has
been demonstrated that such binary evolution does not allow more than around
30% of stars to have been spun up though binary interaction and become
emission line objects.

After using isochrone fits to assign stars in the colour-magnitude diagram
masses, a count of the H$\alpha$ emitters in two open clusters reveals that
for objects near the turn-off, our upper-limit provides a reasonable
description of the numbers of Be stars, especially when uncertainties arising
from the counting method are taken into account. The upper limit does however
fail to describe the sudden decrease in Be fraction that both clusters exhibit
at a mass approximately 70-80% of the turn-off mass.

This problem can be rectified by assuming that systems of low mass-ratio and
low primary mass merge. By adopting simple, although physically justified
stable mass-transfer criteria, we have shown that a good fit to the
observational data is produced by this postulate.

It has been demonstrated in a qualitative way that in coeval populations, a
larger mass-gain of the donor results in a smaller Be fraction at a given
mass. Given that the observed Be fractions are very close to our upper limit
when assuming totally inefficient mass-transfer, it follows that to be able to
explain such high Be fractions, mass-transfer must be non-conservative.

We have highlighted the distinction between the initial binary fraction and
the binary fraction that one is able to observe, and argued that these two
quantities are not equal. This is so because a population of binary stars will
always contain post-interaction systems that will appear to be single stars.
The calculations outlined in this work provide rough constraints on this
discrepancy, suggesting that the initial binary fraction is much higher than
previously thought.

In conclusion, our theoretical argument serves to reinforce numerous
observational arguments that suggest binary interactions to be responsible for
Be stars. We conclude that observations of Be stars in young open clusters
(Milone et al., 2018; Bodensteiner et al., 2020a) do not contradict the
hypothesis that Be stars originate exclusively from mass-transfer in binary
systems. We have shown that if all Be stars are binary interaction products,
somewhat extreme assumptions must be realised such as an initial binary
fraction very close to unity, a shallow initial mass function and very non-
conservative mass transfer. Whether or not these conditions can be met by the
stars in the sky remains to be seen.

Acknowledgements: The authors extend gratitude to an anonymous referee for
useful comments on an earlier version of this manuscript.

## References

  * Ahmed & Sigut (2017) Ahmed, A. & Sigut, T. A. A. 2017, MNRAS, 471, 3398
  * Bestenlehner et al. (2020) Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS [arXiv:2009.05136]
  * Bodensteiner et al. (2020a) Bodensteiner, J., Sana, H., Mahy, L., et al. 2020a, A&A, 634, A51
  * Bodensteiner et al. (2020b) Bodensteiner, J., Shenar, T., & Sana, H. 2020b, A&A, 641, A42
  * Boubert & Evans (2018) Boubert, D. & Evans, N. W. 2018, MNRAS, 477, 5261
  * Brandt & Podsiadlowski (1995) Brandt, N. & Podsiadlowski, P. 1995, MNRAS, 274, 461
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115
  * Brož et al. (2021) Brož, M., Mourard, D., Budaj, J., et al. 2021, A&A, 645, A51
  * Carciofi et al. (2012) Carciofi, A. C., Bjorkman, J. E., Otero, S. A., et al. 2012, ApJ, 744, L15
  * Castro et al. (2014) Castro, N., Fossati, L., Langer, N., et al. 2014, A&A, 570, L13
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Claret & Torres (2016) Claret, A. & Torres, G. 2016, A&A, 592, A15
  * Coe et al. (2020) Coe, M. J., Kennea, J. A., Evans, P. A., & Udalski, A. 2020, MNRAS, 497, L50
  * Collins & Truax (1995) Collins, George W., I. & Truax, R. J. 1995, ApJ, 439, 860
  * de Mink et al. (2011a) de Mink, S. E., Langer, N., & Izzard, R. G. 2011a, Bulletin de la Societe Royale des Sciences de Liege, 80, 543
  * de Mink et al. (2011b) de Mink, S. E., Langer, N., & Izzard, R. G. 2011b, in Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits, ed. C. Neiner, G. Wade, G. Meynet, & G. Peters, Vol. 272, 531–532
  * de Mink et al. (2007) de Mink, S. E., Pols, O. R., & Hilditch, R. W. 2007, A&A, 467, 1181
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7
  * Dorigo Jones et al. (2020) Dorigo Jones, J., Oey, M. S., Paggeot, K., Castro, N., & Moe, M. 2020, ApJ, 903, 43
  * Dufton et al. (2013) Dufton, P. L., Langer, N., Dunstall, P. R., et al. 2013, A&A, 550, A109
  * Dunstall et al. (2011) Dunstall, P. R., Brott, I., Dufton, P. L., et al. 2011, A&A, 536, A65
  * Ekström et al. (2008) Ekström, S., Meynet, G., Maeder, A., & Barblan, F. 2008, A&A, 478, 467
  * Espinosa Lara & Rieutord (2011) Espinosa Lara, F. & Rieutord, M. 2011, A&A, 533, A43
  * Ferrario et al. (2009) Ferrario, L., Pringle, J. E., Tout, C. A., & Wickramasinghe, D. T. 2009, MNRAS, 400, L71
  * Glebbeek & Pols (2008) Glebbeek, E. & Pols, O. R. 2008, A&A, 488, 1017
  * Gossage et al. (2019) Gossage, S., Conroy, C., Dotter, A., et al. 2019, ApJ, 887, 199
  * Hastings et al. (2020) Hastings, B., Wang, C., & Langer, N. 2020, A&A, 633, A165
  * Hjellming & Webbink (1987) Hjellming, M. S. & Webbink, R. F. 1987, ApJ, 318, 794
  * Hogeveen (1991) Hogeveen, S. J. 1991, The mass-ratio distribution of binary stars, PhD thesis, -
  * Huang et al. (2010) Huang, W., Gies, D. R., & McSwain, M. V. 2010, ApJ, 722, 605
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
  * Iqbal & Keller (2013) Iqbal, S. & Keller, S. C. 2013, MNRAS, 435, 3103
  * Justham et al. (2014) Justham, S., Podsiadlowski, P., & Vink, J. S. 2014, ApJ, 796, 121
  * Kamann et al. (2020) Kamann, S., Bastian, N., Gossage, S., et al. 2020, MNRAS, 492, 2177
  * Keller et al. (1999) Keller, S. C., Wood, P. R., & Bessell, M. S. 1999, A&AS, 134, 489
  * Kippenhahn (1974) Kippenhahn, R. 1974, in Late Stages of Stellar Evolution, ed. R. J. Tayler & J. E. Hesser, Vol. 66, 20
  * Klement et al. (2019) Klement, R., Carciofi, A. C., Rivinius, T., et al. 2019, ApJ, 885, 147
  * Klement et al. (2017) Klement, R., Carciofi, A. C., Rivinius, T., et al. 2017, A&A, 601, A74
  * Klencki et al. (2021) Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2021, A&A, 645, A54
  * Kobulnicky & Fryer (2007) Kobulnicky, H. A. & Fryer, C. L. 2007, ApJ, 670, 747
  * Köhler et al. (2006) Köhler, R., Petr-Gotzens, M. G., McCaughrean, M. J., et al. 2006, A&A, 458, 461
  * Kouwenhoven et al. (2009) Kouwenhoven, M. B. N., Brown, A. G. A., Goodwin, S. P., Portegies Zwart, S. F., & Kaper, L. 2009, A&A, 493, 979
  * Kouwenhoven et al. (2007) Kouwenhoven, M. B. N., Brown, A. G. A., Portegies Zwart, S. F., & Kaper, L. 2007, A&A, 474, 77
  * Kriz & Harmanec (1975) Kriz, S. & Harmanec, P. 1975, Bulletin of the Astronomical Institutes of Czechoslovakia, 26, 65
  * Labadie-Bartz et al. (2017) Labadie-Bartz, J., Pepper, J., McSwain, M. V., et al. 2017, AJ, 153, 252
  * Langer (2012) Langer, N. 2012, ARA&A, 50, 107
  * Langer et al. (2020) Langer, N., Schürmann, C., Stoll, K., et al. 2020, A&A, 638, A39
  * Lennon et al. (2005) Lennon, D. J., Lee, J. K., Dufton, P. L., & Ryans, R. S. I. 2005, A&A, 438, 265
  * Li et al. (2012) Li, K. L., Kong, A. K. H., Charles, P. A., et al. 2012, ApJ, 761, 99
  * Liu et al. (2006) Liu, Q. Z., van Paradijs, J., & van den Heuvel, E. P. J. 2006, A&A, 455, 1165
  * Maeder et al. (1999) Maeder, A., Grebel, E. K., & Mermilliod, J.-C. 1999, A&A, 346, 459
  * Marino et al. (2018) Marino, A. F., Przybilla, N., Milone, A. P., et al. 2018, AJ, 156, 116
  * Martayan et al. (2010) Martayan, C., Baade, D., & Fabregat, J. 2010, A&A, 509, A11
  * McSwain & Gies (2005) McSwain, M. V. & Gies, D. R. 2005, ApJS, 161, 118
  * Milone et al. (2018) Milone, A. P., Marino, A. F., Di Criscienzo, M., et al. 2018, MNRAS, 477, 2640
  * Packet (1981) Packet, W. 1981, A&A, 102, 17
  * Pavlovskii & Ivanova (2015) Pavlovskii, K. & Ivanova, N. 2015, MNRAS, 449, 4415
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Peters (1989) Peters, G. J. 1989, Space Sci. Rev., 50, 9
  * Petrovic et al. (2005) Petrovic, J., Langer, N., & van der Hucht, K. A. 2005, A&A, 435, 1013
  * Pinsonneault et al. (1989) Pinsonneault, M. H., Kawaler, S. D., Sofia, S., & Demarque, P. 1989, ApJ, 338, 424
  * Pols et al. (1991) Pols, O. R., Cote, J., Waters, L. B. F. M., & Heise, J. 1991, A&A, 241, 419
  * Porter (1996) Porter, J. M. 1996, MNRAS, 280, L31
  * Porter & Rivinius (2003) Porter, J. M. & Rivinius, T. 2003, PASP, 115, 1153
  * Raguzova & Popov (2005) Raguzova, N. V. & Popov, S. B. 2005, Astronomical and Astrophysical Transactions, 24, 151
  * Rímulo et al. (2018) Rímulo, L. R., Carciofi, A. C., Vieira, R. G., et al. 2018, MNRAS, 476, 3555
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Sana et al. (2013) Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Scalo (1986) Scalo, J. M. 1986, Fund. Cosmic Phys., 11, 1
  * Schneider et al. (2015) Schneider, F. R. N., Izzard, R. G., Langer, N., & de Mink, S. E. 2015, ApJ, 805, 20
  * Schneider et al. (2019) Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211
  * Schneider et al. (2018) Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science, 359, 69
  * Schootemeijer et al. (2018) Schootemeijer, A., Götberg, Y., de Mink, S. E., Gies, D., & Zapartas, E. 2018, A&A, 615, A30
  * Schootemeijer et al. (2019) Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132
  * Secchi (1866) Secchi, A. 1866, Astronomische Nachrichten, 68, 63
  * Shao & Li (2014) Shao, Y. & Li, X.-D. 2014, ApJ, 796, 37
  * Shao & Li (2016) Shao, Y. & Li, X.-D. 2016, ApJ, 833, 108
  * Shenar et al. (2020) Shenar, T., Bodensteiner, J., Abdul-Masih, M., et al. 2020, A&A, 639, L6
  * Sigut et al. (2009) Sigut, T. A. A., McGill, M. A., & Jones, C. E. 2009, ApJ, 699, 1973
  * Sternberg et al. (2003) Sternberg, A., Hoffmann, T. L., & Pauldrach, A. W. A. 2003, ApJ, 599, 1333
  * Struve (1931) Struve, O. 1931, ApJ, 73, 94
  * Sun et al. (2019) Sun, W., de Grijs, R., Deng, L., & Albrow, M. D. 2019, ApJ, 876, 113
  * Tokovinin & Moe (2020) Tokovinin, A. & Moe, M. 2020, MNRAS, 491, 5158
  * Tout (1991) Tout, C. A. 1991, MNRAS, 250, 701
  * Tout et al. (1997) Tout, C. A., Aarseth, S. J., Pols, O. R., & Eggleton, P. P. 1997, MNRAS, 291, 732
  * Townsend et al. (2004) Townsend, R. H. D., Owocki, S. P., & Howarth, I. D. 2004, MNRAS, 350, 189
  * Trimble (1990) Trimble, V. 1990, MNRAS, 242, 79
  * van Bever & Vanbeveren (1997) van Bever, J. & Vanbeveren, D. 1997, A&A, 322, 116
  * Vinciguerra et al. (2020) Vinciguerra, S., Neijssel, C. J., Vigna-Gómez, A., et al. 2020, MNRAS, 498, 4705
  * von Zeipel (1924) von Zeipel, H. 1924, MNRAS, 84, 665
  * Wang et al. (2020) Wang, C., Langer, N., Schootemeijer, A., et al. 2020, ApJ, 888, L12
  * Wellstein et al. (2001) Wellstein, S., Langer, N., & Braun, H. 2001, A&A, 369, 939
  * Wickramasinghe et al. (2014) Wickramasinghe, D. T., Tout, C. A., & Ferrario, L. 2014, MNRAS, 437, 675
  * Zahn et al. (2010) Zahn, J. P., Ranc, C., & Morel, P. 2010, A&A, 517, A7
  * Zorec et al. (2005) Zorec, J., Frémat, Y., & Cidale, L. 2005, A&A, 441, 235
  * Zorec et al. (2016) Zorec, J., Frémat, Y., Domiciano de Souza, A., et al. 2016, A&A, 595, A132

##  Appendix A Stars at the critical velocity

Here we investigate how critically rotating and slowly rotating stars vary in
terms of their effective temperatures and luminosities. Owing to the fact that
stable and reliable numerical models of very fast rotating stars are difficult
to produce, we shall employ a simple analytic model.

The luminosity of a main-sequence star is generated from nuclear reactions in
the central region. As stars are centrally condensed, when a given star is
spun up to the critical velocity the centripetal forces acting in the central
regions are much weaker than the force of gravity, meaning the structure of
the core is largely unchanged. Hence, the intrinsic luminosity is constant to
first order.

What does however change when a star is spun up is the outer structure. The
equatorial radius increases, and therefore so does the surface area of the
star, $S$. This then causes the effective temperature to decrease, as
evidenced by the Stefan–Boltzmann law

$\displaystyle L\propto ST_{\rm{eff}}^{4}.$ (35)

To characterise this change in effective temperature we use the Roche model,
which describes a star with all mass concentrated at the centre, and rotating
with constant angular velocity $\Omega$. In this framework, the effective
potential, with respect to the radial coordinate $r$ and latitude $\theta$, is

$\displaystyle\Psi(r,\theta)=-\frac{GM}{r}-\frac{1}{2}\Omega^{2}r^{2}sin^{2}(\theta).$
(36)

At the critical angular velocity, the polar radius, $R_{p}$ is equal to the
radius of an equivalent non-rotating star of the same mass, whereas the
equatorial radius, $R_{e}$ is given by

$\displaystyle R_{e}=\frac{3}{2}R_{p}.$ (37)

Taking the $x-y$ plane to be parallel to the axis of rotation, where $y$
represents the distance along the rotation axis and $x$ the perpendicular
distance from the rotation axis, the surface of a critically rotating star is
described by

$\displaystyle\left(\frac{y}{R_{e}}\right)^{2}=\left(\frac{2}{3-x^{2}/R_{e}^{2}}\right)^{2}-\left(\frac{x}{R_{e}}\right)^{2},$
(38)

(Zahn et al., 2010) The surface area of a star is (see Paxton et al., 2019,
Appendix B) , in general

$\displaystyle
S=4\pi\int_{0}^{R_{e}}x\sqrt{\left(\frac{dy}{dx}\right)^{2}+1}\,dx.$ (39)

After choosing units such that $R_{e}=1$, Eqs. 38 and 39 can be solved
numerically to give the surface area of a critical rotator, $S_{c}$, as

$\displaystyle S_{c}\approx 4\pi\times 0.7028.$ (40)

Compare this to the surface areas of non-rotating stars with radii $R_{e}$ and
$R_{p}$,

$\displaystyle S_{0}(r=R_{e})=4\pi$ (41)

and

$\displaystyle S_{0}(r=R_{p})=4\pi\left(\frac{2}{3}\right)^{2}\approx
4\pi\times 0.4444.$ (42)

As expected we have

$\displaystyle S_{o}(r=R_{p})<S_{c}<S_{o}(r=R_{e}).$ (43)

We see that the surface area of a critically rotating star is around 1.58
times larger than that of its non-rotating counterpart.

Therefore from Eq. 35, the temperature of a star after having been spun up to
critical decreases by a factor of 0.89.

##  Appendix B Stellar isochrones

As our model predictions give the Be fraction as a function of mass, to make
an effective comparison with observations, we must extract masses from stars
in the colour-magnitude diagram. To this end, we employ isochrones of single
rotating stars to assign mass ranges to different areas of the colour
magnitude diagram.

We use the grid of Schootemeijer et al. (2019) that has been extended to
masses between 2 and 20$\,\mathrm{M_{\odot}}$ with slight changes to internal
mixing - see below. The code used was MESA (Paxton et al., 2011, 2013, 2015,
2018, 2019). Models were computed at initial rotation rates between 0 and 80%
of critical velocity in steps of 10%. As is standard in MESA the critical
velocity is defined as

$\displaystyle v_{\textrm{crit}}=\sqrt{\frac{GM}{R}(1-\Gamma)},$ (44)

where $\Gamma$ is the ratio of luminosity to Eddington luminosity, and is
negligible for the models presented here. During early times, the models
undergo a relaxation period, during which the critical velocity fraction can
oscillate wildly. To circumvent this, we define the initial critical velocity
fraction at the point when the model has burnt 3% of its initial hydrogen
content by mass.

The physics employed in the models is mostly identical to that of Brott et al.
(2011), except for the treatment of two mixing processes. Stepped convective
overshooting is adopted that extends the convective zone by
$\alpha_{\textrm{OV}}$ times the local pressure scale-height. A dependence of
$\alpha_{\textrm{OV}}$ with mass accounts for observational trends (Castro et
al., 2014; Claret & Torres, 2016; Schootemeijer et al., 2019), whereby
$\alpha_{\textrm{OV}}$ increases linearly from 0.1 at
1.66$\,\mathrm{M_{\odot}}$ Claret & Torres (2016) to 0.3 at
20$\,\mathrm{M_{\odot}}$ Brott et al. (2011). Furthermore, time smoothing in
rotational mixing is turned off to avoid unrealistically strong mixing.

Isochrones are generated through a series of linear interpolations and are
split up into two equivalent-evolutionary-phases (EEPs). The first phase lasts
until core hydrogen depletion, and the second phase from core hydrogen
depletion until core helium depletion. To compute the parameters of a star in
the first EEP with initial mass $M_{i}$, initial critical velocity fraction
$\varv_{i}$ at time $t$, we first find the time at which this star would
experience core hydrogen exhaustion, $T_{MS}$. In total, four models are used
for the interpolation, two models with initial masses $M_{1}$ and initial
critical velocity fractions $\varv_{1,a}$ and $\varv_{1,b}$ and similarly two
models with initial masses $M_{2}$ and initial critical velocity fractions
$\varv_{2,a}$ and $\varv_{2,b}$. The models are selected such that
$M_{1}<M_{i}<M_{2}$ and $\varv_{1,a}<\varv_{i}<\varv_{1,b}$ and similarly for
$\varv_{2,a}$, $\varv_{2,b}$. For $M_{1}$ and $M_{2}$ we interpolate the
lifetime when initial ${v_{\textrm{rot}}}/{v_{\textrm{crit}}}\thinspace$=
$\varv_{i}$ from these models as is shown in Fig. 7 a. The hydrogen burning
lifetime is then computed as an interpolation in mass between the values of
$M_{1}$ and $M_{2}$, as depicted in Fig. 7 b. For this step, the most accurate
results are obtained when the logarithm of hydrogen burning lifetime is
interpolated against the logarithm of initial mass. Using the interpolated
lifetime, $T_{i}$, of this star with initial mass $M_{i}$ and initial critical
velocity fraction $\varv_{i}$, we define its fractional lifetime as $t/T_{i}$.
This fractional lifetime is the value at which all further interpolations will
be carried out. Next, a given quantity (for the purposes of making isochrones,
the quantities of interest are effective temperature and luminosity), $Q$, is
interpolated at a fractional lifetime of $t/T_{i}$ for the four selected
models, as in 7 c. The penultimate step is to calculate the quantities
$Q_{M_{1}}$, $Q_{M_{2}}$ which represent the values of $Q$ of a star with mass
$M_{1}$, $M_{2}$, initial critical velocity fraction $\varv_{i}$ and
fractional lifetime $t/T_{i}$ by interpolating across initial critical
velocity fraction like in Fig. 7 d. Finally, an interpolation in initial mass
between the quantities $Q_{M_{1}}$ and $Q_{M_{2}}$ is done to produce the
value of the chosen parameter for a star of given mass, initial rotation rate
and age.

To generate the second EEP, the same procedure is used but only with a
different fractional lifetime, namely the fractional helium-burning lifetime,
$t/T_{He}$ such that at core hydrogen exhaustion $t/T_{He}=0$ and at core
helium exhaustion $t/T_{He}=1$.

Absolute magnitudes in Hubble Space Telescope filters are computed by
interpolating tables of synthetic stellar spectra provided by the MIST project
Choi et al. (2016). Apparent magnitudes are then calculated as

$\displaystyle m_{\textrm{F814W}}=M_{\textrm{F814W}}+A_{F814W}+\mu,$ (45)
$\displaystyle m_{\textrm{F336W}}=M_{\textrm{F336W}}+A_{F336W}+\mu,$ (46)

with $\mu$ being the distance modulus and absorption coefficients
$A_{F814W}=2.04E(B-V)$ and $A_{F336W}=5.16E(B-V)$ (Milone et al., 2018).

Figure 7: Schematic representation of the interpolation procedure employed to
produce isochrones. See text for a thorough explanation.