© ESO, 2013

1. Introduction

Algols are short orbital period (from several hours to tens of days) semi-detached binaries composed of a hot B-A main sequence star and a cool F-K giant star (Kopal 1955; Giuricin et al. 1983). Because of the small separation, the initially more massive star (henceforth, the donor) overfills its Roche lobe and matter then channels to the companion (the gainer) through the inner Lagrangian point ℒ₁, by Roche lobe over-flow (RLOF). The duration of the mass-transfer phase (with rates reaching ≈10^-4  M_⊙ yr^-1) is long enough (several 10⁵ yrs) that a significant amount of mass can be accreted by the gainer (up to several solar masses). Observations clearly attest to the reversal of the mass ratio during RLOF with a donor star (now a giant) that is less massive than its main-sequence companion. The global evolution is well established, but many open questions, such as discrepancies between canonical simulations and observations about the final orbital periods and mass ratios (van Rensbergen et al. 2010b, 2011), still remain.

Previous numerical simulations of Algols (Nelson & Eggleton 2001; van Rensbergen et al. 2011) have provided a broad understanding of short-period binaries, focusing on some specific aspects, such as non-conservative mass transfer (Giannuzzi 1981) or the properties of the mass-transfer stream (Flannery 1975). Except for the study of Dervişoğlu et al. (2010), there are, however, basically no simulations that follow the angular-momentum evolution for both the stars and the orbit in detail, and that account for torques arising from tides, mass transfer, magnetic-field effects, disc accretion or direct impact, and spin-down mechanisms.

Moreover, the evolution of critically-rotating accretors is poorly understood: it is well known that the transferred mass carries enough angular momentum (Kruszewski 1967) to rapidly spin the gainer up to its critical spin-angular velocity (Packet 1981): ${\mathrm{\Omega }}_{\mathrm{crit}\mathit{,}\mathrm{\ast }}\mathrm{=}\mathrm{(}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{/}{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{3}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(1)For a typical 6 + 3.6  M_⊙ system with initial period P_init = 2.5 days in the absence of spin-down mechanisms, only 3 percent (0.12 M_⊙) of the total amount of matter transferred by RLOF (more than 5 M_⊙) is enough to spin the gainer up to the critical rotation. However, it is not clear whether the gainer reaches the critical velocity, whether it is spun down due to any braking mechanisms, or if accretion stops (de Mink et al. 2007). Many physical processes have been proposed to spin down the gainer, including tides (Zahn 1977), magnetic braking (Armitage & Clarke 1996; Stȩpień 2000; Dervişoğlu et al. 2010), or the limitation of the accretion of angular momentum through the interaction with an accretion disc (Popham & Narayan 1991; Paczynski 1991; Colpi et al. 1991; Bisnovatyi-Kogan 1993). However, none of these processes have been applied to a binary-star evolution code that consistently follows the stellar spin velocity and all torques owing to mass accretion and mass loss.

It has also been shown that some Algol systems must be non-conservative to comply with observations of mass ratios (for example Refsdal et al. 1974; Massevitch & Yungelson 1975; Mezzetti et al. 1980; Sarna 1993; Nelson & Eggleton 2001). Yet, only mass loss through the outer Lagrangian point ℒ₃ (Sytov et al. 2007), bipolar jets (Ak et al. 2007) or hotspot formation (van Rensbergen et al. 2011) have been proposed to account for this kind of evolution. Non-conservative mass transfer in stellar evolution codes is generally treated via a free parameter, defining the amount of matter accreted by the gainer star compared to the amount of matter passing through the inner Lagrangian point (de Greve 1993; Sarna et al. 1998; Eggleton & Kiseleva-Eggleton 2002; Dervişoğlu et al. 2010).

In this article, we focus on the evolution of the gainer’s spin-angular velocity and non-conservative evolution. To do so, we study a specific case B mass-transfer system (M₁ = 3.6  M_⊙, M₂ = 6  M_⊙, and P_orb = 2.5 days) in which the donor is a hydrogen-shell burning star. Mass transfer starts just after the star leaves the main sequence, and the rate is characterised by a short phase of rapid mass transfer (hereafter rapid phase) where most of the mass is transferred, followed by a long-lasting slow phase (hereafter quiescent phase) where most of the orbital-separation changes take place.

To address the aforementioned problems (i.e., critical rotation of the gainer and non-conservative evolution), we used the binary-star evolution code Binstar (Siess et al. 2013). Section 2 presents the implemented binary input physics; our treatment of rotation and the various spin-down mechanisms and the torques due to magnetic fields, star-disc interactions, and tides. We also describe our implementation of the hotspot physics and its consequences on the mass accretion rate. In Sect. 3, we apply these prescriptions to our specific Algol system and analyse the effect of each torque on the gainer’s spin. A calculation with a hotspot configuration is also presented. In Sect. 4, observational constraints for Algol systems and related objects are discussed, and their impact on the spin-down mechanism investigated. Finally, we discuss the effect of non-conservative evolution on the close environment of the system (i.e., the region impacted by the system mass-loss) and examine the limitations of our assumptions regarding solid-body rotation and disc formation. The conclusions are given in Sect. 5.

2. The Binstar code

Binstar is an extension of the 1-dimensional stellar evolution code Starevol (Siess 2006) that handles the simultaneous calculation of the binary orbital parameters (separation and eccentricity) and the two stellar components. A complete description of the code can be found in Siess et al. (2013).

In this section, we outline our treatment of stellar rotation and the effects of mass transfer on the stellar spins. We specially focus on the torque due to magnetic fields and the formation of a star-disc boundary layer. Hereafter, the subscripts “g” and “d” refer to the gainer (primary star in Algol systems) and donor (secondary star), respectively. The models presented hereafter assume solid-body rotation. The limitations associated with this assumption are discussed in Sect. 4.3.

The RLOF mass-transfer rate follows the prescriptions described by Ritter (1988) and Kolb & Ritter (1990). In their formalism, two regimes have been considered, depending on whether the Roche radius lies in the optically thick or thin layers of the donor star.

Let us finally emphasize that the (solid) rotation profile and all the torques are computed at each iteration during the convergence process, along with the evolution of the separation and mass transfer rates. This procedure brings numerical stability and ensures that the calculation of the orbital elements are in step with the stellar structures.

2.1. Treatment of angular momentum in BINSTAR

Binstar follows the angular-momentum evolution of each star and the orbit. The total angular momentum of the system (J_Σ) is given by ${\mathit{J}}_{\mathrm{\Sigma }}\mathrm{=}{\mathit{J}}_{\mathrm{d}}\mathrm{+}{\mathit{J}}_{\mathrm{g}}\mathrm{+}{\mathit{J}}_{\mathrm{orb}}\mathit{,}$(2)where J_d,g are the stellar spin-angular momenta and J_orb is the orbital component, such that ${\mathit{J}}_{\mathrm{orb}}\mathrm{=}{\mathit{M}}_{\mathrm{d}}{\mathit{M}}_{\mathrm{g}}{\left[\frac{\mathit{Ga}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}{{\mathit{M}}_{\mathrm{\Sigma }}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(3)where a is the binary separation, e is the eccentricity, and M_Σ = M_d + M_g is the total mass of the system. Evolution is conservative if M_Σ and J_Σ are constant during the evolution. The conservation of angular momentum implies that $\mathit{J̇}\mathrm{orb}\mathrm{=}\mathit{J̇}\mathrm{\Sigma }\mathrm{-}\hspace{0.17em}\mathit{J̇}\mathrm{d}\mathrm{-}\hspace{0.17em}\mathit{J̇}\mathrm{g}\mathit{,}$(4)where are the torques applied on each star, is the rate of change of the orbital angular momentum, and represents the angular-momentum loss rate by the system in non-conservative evolution. Torques from mass transfer and tides are explained in detail in Siess et al. (2013) and summarised in Table 1.

2.2. Magnetic torques

Depending on its configuration, a magnetic field can either spin the star up or down. We consider two effects: magnetic-wind braking and disc-locking.

2.2.1. Magnetic-wind braking

Stars expel matter through winds. Assuming the specific angular momentum of the surface layers of the star to be ${\mathit{j}}_{\mathrm{surf}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{\ast }}{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}\mathit{,}$(5)this matter removes angular momentum at a rate $\mathit{J̇}\begin{array}{c}\mathrm{wind}\mathit{,}{B}\\ \mathrm{loss}\mathit{,}\mathrm{\ast }\end{array}\mathrm{=}\mathrm{-}\mathrm{\left|}\mathit{Ṁ}\mathrm{W}\mathrm{\right|}{\mathrm{\Omega }}_{\mathrm{\ast }}{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}\mathit{<}\mathrm{0}\mathit{,}$(6)where Ṁ_W is the wind mass-loss rate, R_∗ is the stellar radius, and Ω_∗ is its spin-angular speed. If the magnetic field co-rotates with the star, particles emitted by the wind can either be trapped within the dead zone (where the magnetic-field lines are closed) or leave the system by flowing along open magnetic-field lines up to the Alfvèn radius R_A (Weber & Davis 1967; Mestel 1968). Beyond R_A where the wind speed equals the Alfvèn speed, matter decouples from the magnetic field and particles can freely escape the system. The angular-momentum loss rate at R_A is given by $\mathit{J̇}\begin{array}{c}\mathrm{wind}\mathit{,}{B}\\ \mathrm{loss}\mathit{,}\mathrm{\ast }\end{array}\mathrm{=}\mathrm{-}\mathrm{\left|}\mathit{Ṁ}\mathrm{W}\mathrm{\right|}{\mathrm{\Omega }}_{\mathrm{\ast }}{\mathit{R}}_{\mathrm{A}}^{\mathrm{2}}\mathit{.}$(7)The torque is efficient because R_A ≫ R_∗ and can be written as (e.g. Dervişoğlu et al. 2010) $\mathit{J̇}\begin{array}{c}\mathrm{wind}\mathit{,}{B}\\ \mathrm{loss}\mathit{,}\mathrm{\ast }\end{array}\mathrm{=}\mathrm{-}{\mathrm{\left[}\mathrm{\right|}\mathit{Ṁ}\mathrm{W}{\mathrm{|}}^{\mathrm{(}\mathrm{4}\mathit{m}\mathrm{-}\mathrm{9}\mathrm{)}}{\mathit{B}}_{\mathrm{\ast }}^{\mathrm{8}}\mathrm{(}\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}{\mathrm{)}}^{-2}{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{8}\mathit{m}}}^{\mathrm{]}}}^{\mathrm{1}\mathit{/}\mathrm{(}\mathrm{4}\mathit{m}\mathrm{-}\mathrm{5}\mathrm{)}}{\mathrm{\Omega }}_{\mathrm{\ast }}\mathit{,}$(8)where m is a parameter characterising the topology of the magnetic field of magnitude B_∗ at the stellar surface. For our purpose, we assume a bipolar configuration (m = 3), which is a reasonable approximation around the Alfvèn radius where the matter is released (Livio & Pringle 1992). In our models, only the magnetic-field strength B_∗ is left as a free parameter. The wind mass-loss rate is determined using the prescription of Reimers (1975). In a case study, we also considered a rotationally enhanced wind mass loss where the modified rate is given by (Maeder & Meynet 2001) $\frac{\mathit{Ṁ}\mathrm{wind}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}}{\mathit{Ṁ}\mathrm{wind}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{=}{\left(\frac{\mathrm{1}\mathrm{-}\mathrm{\Gamma }}{\mathrm{1}\mathrm{-}\frac{{\mathrm{\Omega }}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi G}{\mathit{\rho }}_{\mathrm{m}}}\mathrm{-}\mathrm{\Gamma }}\right)}^{\frac{\mathrm{1}}{\mathit{\alpha }}\mathrm{-}\mathrm{1}}\mathit{,}$(9)where Γ is the ratio of the stellar luminosity over the Eddington luminosity (Eq. (35)), ρ_m the stellar mean density ${\mathit{\rho }}_{\mathrm{m}}\mathrm{=}\frac{\mathrm{3}{\mathit{M}}_{\mathrm{\ast }}}{\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{3}}}\mathit{,}$(10)and α a force multiplier such that α = 0.52, 0.24, 0.17 and 0.15 for log (T_eff ≥ 4.35, 4.30, 4.00 and 3.90) respectively (Lamers et al. 1995).

2.2.2. Disc-locking

If an accretion disc is present, the gainer’s magnetic field anchors into it. Since the disc does not co-rotate with the star, magnetic-field lines are twisted, creating a toroidal magnetic-field component that generates a torque on the star. If the magnetic field is anchored beyond the co-rotation radius R_co, defined as ${\mathit{R}}_{\mathrm{co}}\mathrm{=}{\left(\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathrm{\Omega }}_{\mathrm{\ast }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(11)then the torque is negative and the star is spun down. The inner disc is truncated at the Alfvén radius R_A, because the magnetic pressure is higher than the gas pressure below this region.

We assume a simple power law for the magnetic-field strength of the form B(r) = B_∗(R_∗/r)^m for r > R_A, considering the same magnetic-field geometry parameter as before (m = 3 for a bipolar configuration). The resulting torque is given by (Armitage & Clarke 1996) $\mathit{J̇}\begin{array}{c}\mathrm{disc}\mathit{,}{B}\\ \mathrm{loss}\mathit{,}\mathrm{\ast }\end{array}\mathrm{=}\frac{{\mathit{\mu }}^{\mathrm{2}}}{\mathrm{3}}\mathrm{(}{\mathit{R}}_{\mathrm{A}}^{-3}\mathrm{-}\mathrm{2}{\mathit{R}}_{\mathrm{co}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{R}}_{\mathrm{A}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{)}$(12)with $\mathit{\mu }\mathrm{=}{\mathit{B}}_{\mathrm{\ast }}{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{3}}$. Depending on the location of R_co and R_A, the global torque can be either positive or negative and is most efficient at spinning-down the gainer when R_A = R_co because there are no regions of the disc that rotate faster than the star (remember that below R_A the disc is truncated) in this configuration. The expression of this maximum torque (Stȩpień 2000) is given by $\mathit{J̇}\begin{array}{c}\mathrm{disc}\mathit{,}{B}\\ \mathrm{loss}\mathit{,}\mathrm{\ast }\end{array}\mathrm{=}\mathrm{-}\frac{{\mathit{\mu }}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{\ast }}^{\mathrm{2}}}{\mathrm{3}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}$(13)and does not depend on the disc properties, such as its mass.

2.3. Tidal effects

Tides have two effects on a binary system: they circularise the orbit of eccentric systems and synchronise the stellar spin with the orbital period. The implementation of tides in Binstar follows the prescription described by Zahn (1977) and includes the Zahn (1989) refinement for convective stars in short-period systems (see Siess et al. 2013, for details).

2.4. Accretion disc and star-disc boundary layer

The matter-stream escaping the donor through the ℒ₁ point will free-fall into the gainer’s potential well. If the stream’s minimum distance of approach to the gainer is such that R_min > R_g, then the stream will orbit the gainer until it collides with itself. Subsequently, collisions between particles lead to the formation of a ring at the circularisation radius. The matter will then spread out, and an accretion disc forms (Kruszewski 1967; Lubow & Shu 1975).

The properties of accretion discs around fast-rotating stars have been widely studied (e.g., Paczynski 1991; Popham & Narayan 1991; Colpi et al. 1991; Bisnovatyi-Kogan 1993). Their 1D calculations show that accretion does not stop when the gainer reaches its critical velocity because of the presence of a boundary layer between the star and the disc. The physics of this boundary and its implementation in Binstar are presented in the next section.

2.4.1. The Paczynski model

The main idea of this model (Paczynski 1991) is to assume that the star and the disc behave as one fluid. A polytropic equation of state is used for both the disc and the star with the pressure given by $\mathit{P}\mathrm{=}\mathit{K}{\mathit{\rho }}^{\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}\mathit{n}}\mathit{,}$(14)where K and n are two constants. The hydrostatic equilibrium of the disc is ensured by balancing the gravitational acceleration, centrifugal acceleration, and pressure gradient (in a cylindrical coordinate system in which r is the distance from the rotation axis and z the distance from the equatorial plane); that is, $\mathit{r}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{=}{\frac{\mathit{\partial }\mathrm{\Psi }}{\mathit{\partial r}}|}_{\mathit{z}}\mathrm{+}{\frac{\mathit{\partial }\mathrm{\Psi }}{\mathit{\partial z}}|}_{\mathit{r}}\frac{\mathrm{d}\mathit{h}}{\mathrm{d}\mathit{r}}\mathit{,}$(15)where z = h(r) is the disc thickness as a function of r, Ω is the disc angular velocity and $\mathrm{\Psi }\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}$(16)is the gravitational potential of a point mass with M_∗ being the stellar mass. For a thin disc [h(r) ≪ r], the equation of hydrostatic equilibrium in the z-direction is $\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathit{\partial P}}{\mathit{\partial r}}\mathrm{=}{\frac{\mathit{\partial }\mathrm{\Psi }}{\mathit{\partial z}}|}_{\mathit{r}}\mathrm{\approx }\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{r}}^{\mathrm{3}}}\mathit{h}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}$(17)After some algebra and using Eq. (15), the condition for hydrostatic equilibrium may be re-written as $\frac{\mathit{h}}{\mathit{r}}\frac{\mathrm{d}\mathit{h}}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{h}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}}{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}\mathrm{-}\mathrm{1.}$(18)The second equation expressing the angular-momentum conservation is written as $\mathit{Ṁ}\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{\mathrm{d\Gamma }}{\mathrm{d}\mathit{r}}\mathit{,}$(19)where Ṁ is the mass flux in the disc, j is the specific angular momentum carried with the disc material, and Γ is the torque due to viscous stresses inside the disc. Using the polytropic equation of state, the α-prescription for the disc viscosity (Shakura & Sunyaev 1973), and the equation of the torque acting between two cylinders in a differentially rotating disc (Pringle 1981), Eq. (19) is rewritten as (Paczynski 1991) $\left[\mathrm{4}\mathit{\pi \alpha }\frac{\mathrm{(}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}{\mathrm{)}}^{\mathit{n}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}}{{\mathit{K}}^{\mathit{n}}}\right]\frac{\mathrm{d\Omega }}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{{\mathit{r}}^{\mathrm{3}\mathit{n}\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}}{{\mathit{z}}^{\mathrm{2}\mathit{n}\mathrm{+}\mathrm{3}}}\mathrm{(}\mathit{J̇}\mathrm{-}\mathit{Ṁ}{\mathit{r}}^{\mathrm{2}}\mathrm{\Omega }\mathrm{)}\mathit{.}$(20)

2.4.2. Dimensionless equations

For convenience, the dimensionless variables are introduced:

$\begin{array}{ccc}& & \mathit{x}\mathrm{=}\frac{\mathit{r}}{{\mathit{R}}_{\mathrm{\ast }}}\mathit{,}\mathit{y}\mathrm{=}\frac{\mathit{h}}{{\mathit{R}}_{\mathrm{\ast }}}\mathit{,}\\ & & {\mathit{\omega }}^{\mathrm{2}}\mathrm{=}{\mathrm{\Omega }}^{\mathrm{2}}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{3}}}{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}\mathit{,}{\mathit{\omega }}_{\mathrm{\ast }}^{\mathrm{2}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{\ast }}^{\mathrm{2}}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{3}}}{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}\mathit{,}{\mathit{\omega }}_{\mathrm{crit}}^{\mathrm{2}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{crit}}^{\mathrm{2}}\frac{{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{3}}}{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}\mathit{,}\end{array}$where Ω_∗ and Ω_crit are the stellar and critical angular velocities, respectively. We also define the two constants: $\mathit{\varsigma }\mathrm{\equiv }\mathrm{4}\mathit{\pi \alpha }\frac{{\mathit{G}}^{\mathrm{2}}{\mathit{M}}^{\mathrm{2}}}{\dot{{\mathit{M}}_{\mathrm{\ast }}}{\mathit{K}}^{\mathrm{1.5}}}\mathit{,}$(23)and ${\mathit{j}}_{\mathrm{\ast }}\mathrm{\equiv }\frac{\mathit{J̇}}{\mathit{Ṁ}\mathrm{(}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}{\mathit{R}}_{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathit{,}$(24)where j_∗ represents the specific angular momentum of the disc material at the effective stellar radius R_∗, which is defined as (Chandrasekhar 1939) ${\mathit{R}}_{\mathrm{\ast }}\mathrm{=}\frac{\mathit{K}}{\mathrm{0.4242}\mathit{G}{\mathit{M}}_{\mathrm{\ast }}^{\mathrm{1}\mathit{/}\mathrm{3}}}\mathrm{·}$(25)Following Paczynski (1991), we use ς = 10⁶ throughout the computations. In dimensionless form, Eqs. (18) and (20) can be re-expressed as $\begin{array}{ccc}& & \frac{\mathit{y}}{\mathit{x}}\frac{\mathrm{d}\mathit{y}}{\mathrm{d}\mathit{x}}\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{y}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{-}\mathrm{1}\mathit{,}\\ & & \frac{\mathrm{d}\mathit{\omega }}{\mathrm{d}\mathit{x}}\mathrm{=}\frac{{\mathit{x}}^{\mathrm{3}}}{\mathit{\varsigma }{\mathit{y}}^{\mathrm{6}}}\mathrm{(}{\mathit{j}}_{\mathrm{\ast }}\mathrm{-}\mathit{\omega }{\mathit{x}}^{\mathrm{2}}\mathrm{)}\mathrm{·}\end{array}$This system of equations does not depend on the disc properties, such as mass or radial extent. This system is solved using a two-point boundary value solver (Capper et al. 2007). For numerical reasons, two sets of boundaries are used depending on the stellar spin-angular velocity.

For rapidly rotating stars with ω_∗/ω_crit ≈ 1, the structure of the disc must satisfy the boundary condition $\mathit{\omega }\mathrm{=}{\mathit{\omega }}_{\mathrm{\ast }}\mathrm{at}\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{tr}}\mathit{,}$(28)where r_tr = 0.8r_polar is an arbitrary transition radius between the star and the disc, and r_polar is the star’s polar radius.

For slowly rotating stars, the spin-angular velocity profile presents a local maximum in the boundary layer. This maximum, which is located at a radius r_nt, corresponds to a zero-torque condition between the star and the disc. The inner boundary condition is written as ${\left(\frac{\mathrm{d}\mathit{\omega }}{\mathrm{d}\mathit{r}}\right)}_{{\mathit{r}}_{\mathrm{nt}}}\mathrm{=}\mathrm{0.}$(29)A Runge-Kutta integrator (Press et al. 2007) is used to solve the inner part of the system between r_nt and r_tr.

For both rapidly and slowly rotating systems, the outer boundary condition is $\frac{\mathrm{d}\mathit{\omega }}{\mathrm{d}\mathit{r}}\mathrm{=}\frac{\mathrm{d}{\mathit{\omega }}_{\mathrm{Kepler}}}{\mathrm{d}\mathit{r}}\mathit{,}$(30)where ω_Kepler is the Keplerian angular speed at the outer edge of the disc. Details on the method of calculation can be found in Paczynski (1991).

The solution of this system of equations provides the value of the specific angular momentum j_disc = j_∗ accreted by the star as a function of its stellar spin. The torque exerted on the star is then $\mathit{J̇}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}\mathrm{=}\mathrm{\left|}\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}\mathrm{\right|}{\mathit{j}}_{\mathrm{disc}}\mathit{,}$(31)where $\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}$ is the mass accretion rate for the gainer star. Note that j_disc can be negative.

To ease the computation when the star is not critically rotating (less than 80% of the critical spin), the accreted specific angular momentum can be set as a fraction f_jdisc of the Keplerian specific angular momentum at the surface of the star; that is, ${\mathit{j}}_{\mathrm{disc}}\mathrm{=}{\mathit{f}}_{{\mathit{j}}_{\mathrm{disc}}}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{g}}{\mathit{R}}_{\mathrm{g}}}\mathit{.}$(32)We use f_jdisc = 1, which is the mean value computed with the “boundary-layer” mechanism for non-critically rotating accretors.

In the disc, advection of matter occurs due to the transport of angular momentum by viscosity. By treating the star-disc boundary layer as one fluid, the same mechanism also applies at the stellar surface: when the star reaches critical rotation, the accretion of angular momentum is no longer possible, and viscous processes remove spin-angular momentum from the star, keeping it at a critical value but not exceeding it. This mechanism allows the star to accrete large amounts of mass, while giving angular momentum back to the disc. This extra angular momentum leads to the spreading of the disc. The outer disc radius is truncated because of tidal interaction between the disc and the donor star (Lin & Papaloizou 1979; Whitehurst 1988; Ichikawa & Osaki 1994) and, consequently, some angular momentum is returned to the orbit.

2.5. Direct impact and hotspot formation

A second mode of mass transfer occurs when R_min < R_g, such that the stream impacts the gainer’s surface (hereafter “direct impact”). In this case, we follow the evolution of the stream through a ballistic approach to determine the precise location, velocity, and angular momentum of the impacting stream (Huang 1963; Kruszewski 1964; Flannery 1975).

When the particle collides with the star, its specific angular momentum is ${\mathit{j}}_{\mathrm{impact}}\mathrm{=}\mathrm{\parallel }{r}\mathrm{\wedge }{v}\mathrm{\parallel }\mathit{,}$(33)where r and v are respectively the radius vector of the particle from the gainer’s stellar centre, and the particle velocity at the impact location. These are determined by integrating the equations of motion of the particle using a Runge-Kutta scheme. Using Eq. (33), the torque applied on the accretor is given by $\mathit{J̇}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}\mathrm{=}\mathrm{\left|}\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}\mathrm{\right|}{\mathit{j}}_{\mathrm{impact}}\mathit{.}$(34)When the stream impacts the star (or the disc’s outer radius), its kinetic energy is converted into thermal energy and a hotspot forms (Peters & Polidan 2004; Lomax et al. 2012). If the hotspot luminosity exceeds the Eddington value such that ${\mathit{L}}_{\mathrm{edd}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi cG}{\mathit{M}}_{\mathrm{\ast }}}{\mathit{\kappa }}\mathit{,}$(35)where κ is the plasma opacity, the matter escapes the system. We show in Appendix A (all quantities entering Eq. (36) being defined there) that if the mass-transfer rate exceeds $\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}\mathrm{=}\frac{\mathrm{2}}{{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}}\left(\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}\mathrm{-}{\mathit{L}}_{\mathrm{g}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}}\mathrm{+}\frac{{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}}{\mathrm{\Delta }\mathit{t}}\right)\mathit{,}$(36)the evolution becomes non-conservative (Eggleton 2000; van Rensbergen et al. 2008). We use a prescription for $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$, depending only on the total mass of the system and based on observations (see Appendix A, Eq. (A.20)).

The parameter β defining the fraction of the transferred mass that escapes the system is given by $\mathit{\beta }\mathrm{=}\left|\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{g}\mathit{,}\mathrm{acc}\end{array}}{\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{d}\mathit{,}\mathrm{loss}\end{array}}\right|\mathit{,}$(37)where $\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{g}\mathit{,}\mathrm{acc}\end{array}\mathrm{=}\min \mathrm{\left(}\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{loss}\mathit{,}\mathrm{d}\end{array}\mathit{,}\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}\mathrm{\right)}$. Note that the wind mass-loss contribution is not included.

In short-period binaries, matter may also escape the system through the third Lagrangian point (Sytov et al. 2007) at a rate of approximately 10^-10  M_⊙ yr^-1, which is higher than stellar winds but still negligible compared to the RLOF mass-transfer rate. We neglect this aspect in our present investigation.

3. Results

In this section, we present new calculations of Algol evolution, considering the previously described spin-down mechanisms and non-conservative evolution. Our reference system is a 6 M_⊙ + 3.6  M_⊙ with an initial period of 2.5 days. Both stars are initially relaxed on the zero-age main sequence with composition X_g = 0.735 and X_d = 0.695 and are synchronised with the orbit. Initial effective temperatures (luminosities) are ≈10 500 K (≈130 L_⊙) for the gainer and ≈19 500 K (≈1030 L_⊙) for the donor. In this study, we do not consider Bondi-Hoyle accretion ($\mathit{J̇}\begin{array}{c}\mathrm{wind}\\ \mathrm{acc}\mathit{,}\mathrm{\left[}\mathrm{g}\mathit{,}\mathrm{d}\mathrm{\right]}\end{array}$, see Table 1) because of its negligible impact compared to RLOF mass transfer. By default, no spin-down mechanisms are used unless mentioned.

3.1. On the spin up of the gainer

Fig. 1 Evolution of the mass-transfer rate (top), stellar masses (middle), surface velocity (bottom) for a 6 M⊙ donor (solid black line), and 3.6 M⊙ gainer (dotted red line) system with initial period Pinit = 2.5 days.

Figure 1 shows the evolution of the mass-transfer rate (top panel), stellar masses (middle panel), and the surface-spin velocity of the gainer (in units of the critical Keplerian velocity; bottom panel) for our reference system, undergoing case B mass transfer. Here, we have not considered any spin-down mechanisms and have assumed fully conservative mass transfer. Mass and angular-momentum losses from the system only occur via winds (quasi-conservative evolution).

During the rapid phase, the mass-accretion rate reaches values up to 10^-4  M_⊙ yr^-1, and by about 3 × 10⁴ years, the gainer reaches critical velocity (blue vertical line in Fig. 1). At this point, we stop the evolution. During this period, only a low amount of matter is transferred, and the gainer’s mass increases by 0.12 M_⊙, which represents 3% of its initial mass (Fig. 1, middle panel). However, more than 5 M_⊙ would have been transferred if the evolution was continued.

3.2. Effect of magnetic fields

The efficiency of magnetic-wind braking strongly depends on the mass-loss rate of the gainer (and on the parameter β for non-conservative evolution to a certain extent), and on | B_∗ |. In the following simulations, we assume that the disc-locking activates as soon as the gainer is spun up and is at work even if direct impact is expected. Although contradictory, this assumption does not conflict with an observation of binaries in which the gainer is surrounded by a disc, even though the Lubow & Shu (1975) criterion¹ implies direct impact (see Sect. 4.2.2).

The top panel of Fig. 2 shows the evolution of the gainer’s surface spin-angular velocity (normalised to the critical value) for different values of | B_∗ |, where we consider both magnetic-wind braking and disc-locking. Only a magnetic field of intensity | B_∗ |  ≳ 3 kG can spin down the star sufficiently to prevent critical rotation.

Fig. 2 Evolution of the gainer’s surface velocity (normalised to the critical velocity) for different spin-down mechanisms for our 6 + 3.6 M⊙ system with an initial period Pinit = 2.5 days. The horizontal line represents the critical spin-angular rotation of the star. All simulations are stopped when Ωsurf ≈ Ωcrit. Top panel: solid black line (labelled “Ref.”): | B∗ | = 0 kG, no braking mechanism; dotted red line (labelled “B.-L.”): boundary-layer mechanism, no magnetic field; dashed green line: Zahn (1989) tidal effect treatment only; long-dashed blue line: magnetic braking with | B∗ |  = 3 kG. Bottom panel: solid black line: magnetic braking with | B∗ | = 1 kG; dotted red line: magnetic braking with | B∗ | = 2 kG; dashed green line (labelled M.-W.): magnetic-wind braking only with | B∗ |  = 3 kG; long-dashed blue line (labelled D.-L.): disc locking only with | B∗ |  = 3 kG.

In the lower panel of Fig. 2, we present the evolution of the normalised spin-angular velocity for each magnetic-braking mechanism individually. As shown, wind braking is very inefficient, but various effects can contribute to strengthen its action. First, rotation can substantially enhance the mass loss rate because of the reduced gravity. To quantify this effect, we ran additional simulations using Eq. (9) for the wind loss rate. Although the mass-loss rate is increased by a factor of almost 3, the 2 kG configuration still produces a critically rotating gainer, and in the 3 kG case, the difference in stellar spin is insignificant, of the order of 2 percents. Therefore, change in the wind rate due to rotation has a negligible impact because disc-locking is so efficient (and independent of mass loss rate). Alternatively, Tout & Pringle (1992) argue that the dynamo process at the origin of the magnetic-field generation can enhance the mass loss rate and hence, the magnetic torque. Moreover, the evolution may not be conservative and some mass may leave the system from the hotspot (see Sects. 3.6 and 3.7) or from the outer Lagrangian points ℒ₂ and ℒ₃ if, for example, the Roche lobe geometry is modified by radiation pressure (Dermine et al. 2009). On the other hand, the disc-locking torque has been overestimated, because we set R_co = R_A in Eq. (12) and considered that the whole outer disc is contributing to spin down the star. This is not necessarily the case because the disc magnetosphere can screen the stellar magnetic field, preventing it to anchor into the disc and therefore to spin the star down (Ghosh & Lamb 1979). Note also that not all Algol systems possess a disc during the rapid mass-transfer phase, and therefore, the disc-locking mechanism may not operate universally.

Observations suggest that magnetic fields of 3 kG for main sequence stars are rare even for Ap and Bp stars, which are known to be more magnetically active than standard A-B stars typically found in Algols (Hubrig et al. 2006, 2007; Bychkov et al. 2003, 2009). However, these observations relate to slow rotators, and we may expect the surface magnetic-field strength to have been much larger when the star was spinning faster due to dynamo generation during the rapid phase (Spruit 1999). On the other hand, the (initial) 3.6 M_⊙ gainer star has an extended radiative envelope, which limits the dynamo processes (Piddington 1983), so dynamo effects might not be efficient. In addition, it is difficult to observationally determine the magnetic-field strength of the main-sequence gainer in Algols, because it is likely weaker than that of the red-giant gainer, which possesses an expended convective envelope (Retter et al. 2005). Given the large uncertainties surrounding the modelling of these braking mechanisms, it is difficult to state whether or not a 3 kG magnetic field is a plausible value. Additional magnetic field determinations are clearly needed to better understand the process(es) at work in Algols.

3.3. Effects of tides

Tidal forces steeply increase with decreasing orbital separation and therefore may be efficient in short-period Algols. The dashed green curve of Fig. 2 shows the evolution of the surface velocity (in units of the critical velocity) for our 6 + 3.6 M_⊙ system, which includes tidal braking. We confirm the results of Dervişoğlu et al. (2010) that tides are not efficient enough to spin-down the gainer and compensate for the spin-up produced by angular-momentum accretion because of the much longer tidal timescale (τ_sync ≈ 10⁸ yr) compared to the mass-transfer timescale (${\mathit{\tau }}_{\mathrm{acc}}\mathrm{=}{\mathit{M}}_{\mathrm{g}}\mathit{/}\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}\mathrm{\approx }{\mathrm{10}}^{\mathrm{4}}\mathrm{-}{\mathrm{10}}^{\mathrm{5}}$ yr).

Observations of Algols during the slow accretion phase indicate that they are spinning slower than their critical velocity. Most Algols are synchronised with the orbital period (Wilson 1989; Mukherjee et al. 1996). Thus, tides may help to spin-down the gainer during the long-lasting quiescent phase, but they are too weak to maintain the gainer below the critical spin velocity during the entire mass-transfer episode.

Fig. 3 Evolution of the torque on the gainer star $\mathit{J̇}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}$. Solid black line: no spin down mechanism (the simulation is stopped, when the gainer reaches the critical rotation). Dashed red line: boundary-layer model. The horizontal line represents the limit where the net torque spins the star up (above) or spins it down (below).

3.4. Star-disc boundary-layer treatment

Although disc formation in semi-detached binaries is determined by the stellar and orbital properties of the system (Lubow & Shu 1975), we assume the formation of an accretion disc whenever the star reaches its critical spin rate. Our paradigm is that the stellar matter which, cannot be accreted by the fast-rotating and distorted star, will surround the gainer, eventually forming a disc (see Sect. 4.2.2). In practice, the treatment of the boundary layer activates once the star has reached 80% of its critical spin-angular velocity. As a result, the gainer is initially spun up in the same way as if no braking mechanism was present and is then kept at the critical spin-angular velocity when the boundary layer is at work.

Figure 3 shows the evolution of the torque applied on the star due to mass accretion. When the boundary-layer mechanism is activated, the torque quickly drops to zero when the star reaches the critical spin velocity. Since the momentum of inertia of the donor keeps rising due to the increase in mass and radius, angular momentum needs to be evacuated to maintain the rotation at the critical value, hence the temporarily negative value of $\mathit{J̇}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}$ (Fig. 2). All the transferred matter is accreted, and the evolution of the system remains conservative.

When a disc forms because the separation is large enough or when the gainer has a small filling-factor, the boundary model is activated, and as long as the gainer’s surface velocity remains below 0.9Ω_crit, the accreted specific angular momentum is basically equal to the Keplerian value of ${\mathit{j}}_{\mathrm{crit}}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}{\mathit{R}}_{\mathrm{\ast }\mathit{,}\mathrm{eq}}}\mathit{,}$(38)where R_∗ ,eq is the stellar equatorial radius. Note that our computed values for j_crit are in very good agreement with those of Popham & Narayan (1991), who use a different set of equations.

3.5. Evolution of the system

Fig. 4 Evolution of our 6 + 3.6 M⊙ system with different braking mechanisms: (solid black line) rotation-free model; (dotted red line) magnetic field 3 kG; (dashed green line) boundary-layer treatment. Top left: mass accretion rate on the gainer. The vertical line represents the epoch when the mass ratio q = 1. Top right: orbital separation. Bottom left: HR diagram for the three configurations: (solid black) rotation-free model; (dotted red) 3 kG magnetic field; (dashed blue) boundary-layer treatment; (dot-dashed green) hotspot (see Appendix A). The open circles correspond to the beginning of mass transfer, triangles to the transition between the rapid and quiescent phase of mass transfer, and squares to the end of mass transfer (only on the rotation-free model). Bottom right: extension of the donor (top cyan hatched area) and gainer (bottom magenta hatched area) within their Roche lobes RL (donor (solid blue lines) gainer (solid red lines)) for the three cases. The y-axis represents the position relative to the system’s centre of mass. The gap between the Roche lobe radii RL of the two stars is because RL is defined as the radius of a sphere with a volume equivalent to the Roche lobe and therefore differs from the true location of the ℒ1 point.

The evolution of the binary parameters strongly depends on the spin-down mechanism, because it controls how much mass and angular momentum are transferred.

The top right panel of Fig. 4 shows the evolution of the separation for the considered spin-down mechanisms. The “boundary-layer” simulation ends up with a shorter orbital separation (about 100 R_⊙, orbital period of 36 days) than the magnetic-field case (about 135 R_⊙, orbital period of 59 days), which is slightly less than the canonical “rotation-free” model.

When a boundary layer forms, the star accretes the maximum allowed amount of angular momentum necessary to maintain its rotational velocity at the critical value. This process minimises the angular momentum returned to the orbit, hence leading to a smaller orbital separation compared to the other mechanisms. For the disc-locking mechanism, the gainer’s spin is kept below the critical rate, so more angular momentum is evacuated leading to longer period systems. For the stronger wind-braking mechanism (as compared to disc-locking), more angular momentum will be lost from the system, therefore reducing the orbital separation.

The differences in the mass-transfer rates between the three aforementioned models are small, (see Fig. 4, top left panel) since the separation does not vary between the considered models during the rapid phase of mass transfer. Therefore, the Roche radius – and in turn the mass-transfer rate – of the donor remain unchanged. Then, the differences in the separation are negligible on the (weak) mass-transfer rate during the quiescent phase. As a result, the final masses are exactly the same. The donor star ends up with a mass of 0.8 M_⊙ and the gainer at about 8.75 M_⊙. The final mass ratios are q = 0.096 for the star-disc spin-down mechanism and q = 0.098 for the wind-braking mechanism.

During the mass-transfer phase, the donor radius equals the Roche radius, which depends on the separation and thus on the braking mechanism. As the Roche radius is smaller for the boundary-layer model than for the magnetic-braking mechanism, so is the donor’s radius (by about 30%; bottom-right panel of Fig. 4).

The bottom-left panel of Fig. 4 presents the Hertzsprung-Russell diagram (HRD) for the same three configurations. There are no differences between these models for the gainer star since roughly the same amount of matter is accreted: the radiative accretor ascends the main sequence at a slightly higher luminosity. On the other hand, the donor’s evolution exhibits small differences during the quiescent mass-transfer phase due to the different evolutionary histories of the separation and stellar radius. By the end of the simulation, all tracks converge to the same point.

Fig. 5 Evolution of our 6 + 3.6 M⊙ system with hotspot formation. Top left: normalised hotspot Lhotspot (solid black line), accretion ${\mathit{L}}_{\mathrm{acc}}^{\mathrm{std}}$ (dashed red) and Eddington (dotted green) luminosities. Top right: evolution of the masses with (solid black) and without (dashed blue) our hotspot formalism compared to van Rensbergen et al. (2008) hotspot formalism (dotted red). Bottom left: evolution of β for the gainer using our (solid black) and van Rensbergen et al. (2008) hotspot formalism (dotted red). Bottom right: evolution of the orbital separation (same colours as top right panel).

3.6. Hotspot

The hotspot formalism derived in Appendix A provides a self-consistent determination of the β parameter. In this section, it is applied in conjunction with the boundary-layer formalism to avoid super-critical velocities. During the rapid phase, the hotspot luminosity L_hotspot (Fig. 5, top left) is higher than the critical Eddington value $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}$ (Eq. (A.15)) and largely exceeds the accretion luminosity${\mathit{L}}_{\mathrm{acc}}^{\mathrm{std}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{g}}\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}}{{\mathit{R}}_{\mathrm{g}}}$(39)because of the small impact area (Eq. (A.17)). This hotspot formalism leads to a non-conservative evolution and small values of β during both the rapid and quiescent phases (see bottom-left panel of Fig. 5). While almost the same amount of mass is lost by the donor² (top-right panel of Fig. 5), the gainer’s final mass is different: 7.4 M_⊙ versus 8.8 M_⊙ for the fully conservative case. In total, 1.4 M_⊙ are expelled from the system, but this barely affects the final mass ratio; 0.1 for the conservative model compared to 0.12 for the hotspot model.

There are two regimes for the evolution of β (see Fig. 5, bottom left). The first one occurs during the onset of mass transfer, when the gainer is still slowly rotating (Ω_∗/Ω_crit ≲ 0.2). This phase is characterised by a drop in β as the first term of Eq. (36), $\mathrm{(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}\mathrm{-}{\mathit{L}}_{\mathrm{g}}\mathrm{)}\mathit{/}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{K}\end{array}$, dominates the second (rotational) term, $\mathrm{\Delta }{\mathit{E}}_{\mathrm{rot}}\mathit{/}\mathrm{\Delta }\mathit{t}\mathrm{=}\mathrm{(}{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}\mathrm{)}\mathit{/}\mathrm{\Delta }\mathit{t}$. With increasing velocity, the star enters the second regime where the rotational term, which is proportional to Ω_∗, becomes dominant. During this phase, β starts to increase again, because $\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}$ gets higher. Eventually, the star reaches its critical spin velocity, and β saturates around 0.6−0.7. The noise in the β profile is due to the explicit dependence of the rotational term on the time-step Δt.

Once the mass ratio has been reversed, the orbital separation starts to increase, leaving space for the formation of an accretion disc. At time t ≈ 5.486 × 10⁷ yr, the stream no longer impacts the star, and our hotspot formalism stops applying, terminating the non-conservative evolution.

In the Hertzsprung-Russell diagram (Fig. 4, bottom-left panel), the gainer is less luminous, because it accretes less mass. For the donor, on the other hand, the presence of a hotspot does not alter its evolution in the HRD. The main effects of this non-conservative evolution are the slight decrease in the final orbital separation (94 R_⊙ vs. 100 R_⊙ in the conservative case) and a slower acceleration of the gainer’s rotational velocity.

The parameter $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$, defining the hotspot luminosity (Eq. (A.18)), has a strong impact on the critical mass-accretion rate and thus on β. Given the range of values for $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$ derived from observations (from $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{K}\end{array}\mathrm{\approx }\mathrm{1}$ to $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{K}\end{array}\mathrm{\approx }\mathrm{700}$; van Rensbergen et al. 2011), a fully conservative or a significantly non-conservative evolution can result from the simulations. In our calculations, $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$ varies between 125 and 161, and the total amount of mass ejected is 1.4 M_⊙. If we were to use $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{K}\end{array}\mathrm{=}\mathrm{1}$, we would have a fully conservative evolution for the same system.

The opacity of the material at the impact location (where the Eddington luminosity is computed) is set to the photospheric value that is higher than the Thomson electron-scattering opacity κ_es generally used. Given the sensitivity of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}$ to this quantity, we ran a simulation using κ_es for the computation of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}$ and found a fully conservative evolution. The opacity has thus a strong impact on β. It is not straightforward (and beyond the scope of this article) to precisely compute the physics (opacity, shock formation) at the impact location, since this depends on the temperature of the disturbed material and the penetration depth of the stream. This depth can be estimated by equating the stream ram pressure ρ | v² | to the stellar gas pressure (Ulrich & Burger 1976). However, complications arise because the rotation of the star makes the stream impact upon a new undisturbed part of the star at each instant. Moreover, the matter can be ejected directly from the impact location but also from its trail (i.e., from the perturbed material no longer at the impact location) where the opacity may be different.

Due to the large amount of mass transferred in massive short-period binaries, the gainer’s radius may increase until it fills its own Roche radius and the system may then evolve into contact. However, the presence of a hotspot can prevent this evolution, because less mass is accreted onto the gainer star.

In contrast to our simulation, the model of van Rensbergen et al. (2008) for the same system (with initial masses 6 + 3.6  M_⊙ and initial period P_init = 2.5 d) is fully conservative (neglecting winds). It is not straightforward to disentangle the differences between the two runs because of the different approaches used to compute various physical quantities. For instance, the stream velocity is computed in Binstar via ballistic motion, while van Rensbergen et al. (2008) use geometrical arguments. They also use a parametric tidal spin-down prescription, which is more efficient than the one computed in Binstar, derived from the stellar structure. Therefore, the gainer star is more efficiently spun down and, in turn, has a lower mass-loss rate due to the higher $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}$ (Eq. (A.16)) in their models. The computation of the opacity also differs; we use the photospheric opacity (at τ = 2/3), while they use a parametric expression based on the electron-scattering opacity. For comparison, we implemented their formulation in Binstar and found that less mass is lost from the system (≈0.9 M_⊙) mainly, because they neglect the term ${\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}$ in Eq. (36), which is not negligible (${\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\mathit{/}{\mathit{v}}^{\mathrm{2}}\mathrm{\approx }\mathrm{0.2}$) and contributes to a greater decrease $\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}$. As a result, the final separation (bottom-right panel of Fig. 5) also changes between the two hotspot formalisms since it is linked with the masses of the two stars. Our less conservative model leads to a shorter-period system. Nevertheless, both formalisms show the same two phases for the evolution of β.

Fig. 6 Evolution of our 6 + 3.6  M⊙ system with Pinit = 2.5 days with hotspot and different magnetic field strengths. Black: hotspot with boundary layer (BL) and no magnetic field. Red: hotspot (no BL) + magnetic field | B∗ |  = 1 kG. Green: hotspot (no BL) + magnetic field | B∗ |  = 2 kG. The red simulation stops when the gainer reaches its critical spin-angular momentum (indicated by a red circle). Top left: surface spin-angular velocity. Top right: orbital separation. Bottom left: accretion efficiency β. Bottom right: masses.

3.7. Hotspot and magnetic braking

We mentioned that β is a key parameter for the magnetic-wind braking. In Sect. 3.2, only stellar winds have been considered as a mass-loss mechanism from the system, but the matter emanating from the hotspot might contribute to the magnetic-wind braking by increasing the mass-loss rate Ṁ_W in Eq. (8). As a result, a lower magnetic-field strength is needed to keep the gainer’s spin velocity below the critical value.

In the configuration with a hotspot plus a 1 kG magnetic field and no boundary-layer mechanism, the gainer reaches the critical rate, although later than in the conservative case (Fig. 6). Above 2 kG, the critical rotation can be avoided but the required field strength still remains higher than measured in Algols.

As expected, coupling the expelled material to the magnetic field will increase the system angular momentum loss and thus the separation. As a consequence, the impact velocity will be greater and according to Eq. (36), the threshold for mass ejection $\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}$ will be lower. As explained previously, the value of β drops until the rotational term $\frac{{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}}{\mathrm{\Delta }\mathit{t}}$ starts to dominate. The β profiles diverge more as time goes on because of the increasing difference in the gainer’s luminosities due to their different masses. In the case of a strong magnetic field (2 kG), β drops to lower values (β ≈ 0.1) and for a much longer duration. The donor reaches the same final mass (M_d ≈ 0.8  M_⊙), but the final mass ratio increases significantly from 0.1 to 0.25 because 3.6 M_⊙ are removed from the system.

In their study, Dervişoğlu et al. (2010) use a constant value of β ranging from 0.1 to 0.9 to account for the system mass-loss. This is in contrast to our approach where β is computed consistently. Our simulations cannot maintain values as low as those used by these authors during the entire phase of mass transfer and this explains why in our case the magnetic-wind braking is not as efficient. Because of our poor knowledge of the hotspot outflow geometry and its link with the magnetic field, that coupling should be regarded as exploratory.

4. Can observations help constraining Algol evolution?

4.1. Predicting the mass ejected from the hotspot

Fig. 7 Density profile of the matter expelled by the hotspot and surrounding the system. The x-axis represents the angular extent for a system located at 100 pc. The profile starts at the stellar radius of the gainer (R∗ = 2.2 × 1011 cm) and is cut at 2 arcmin. The density drops below the mean density of the ISM value at 12 arcmin (1.05 × 1018 cm) from the star.

It is not clear what happens to the matter ejected via the hotspot, whether it freely escapes the system as assumed in our simulation or if it remains in the gainer’s Roche lobe to be eventually re-accreted. Since the hotspot luminosity exceeds the Eddington luminosity in our model, we expect the matter to reach the escape velocity ${\mathit{v}}_{\mathrm{esc}}\mathrm{\left(}{\mathit{R}}_{\mathrm{\ast }}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{eff}}}{{\mathit{R}}_{\mathrm{\ast }}}}\mathit{,}$(40)where M_eff is the effective gravitational mass felt by the wind (Castor et al. 1975): ${\mathit{M}}_{\mathrm{eff}}\mathrm{=}{\mathit{M}}_{\mathrm{\ast }}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{L}}\mathrm{)}\mathrm{=}{\mathit{M}}_{\mathrm{\ast }}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{\kappa }}_{\mathrm{es}}{\mathit{L}}_{\mathrm{\Sigma }}}{\mathrm{4}\mathit{\pi cG}{\mathit{M}}_{\mathrm{\ast }}}\right)\mathit{,}$(41)where κ_es is the Thomson electron-scattering opacity. In this equation, the luminosity L_Σ is due to both the star and the hotspot. We further assume that the flow, once ejected, is accelerated by radiation pressure. Following Castor et al. (1975)’s prescription for radiation-driven winds, the material velocity at a distance r from the gainer is then given by $\mathit{v}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathit{\alpha }}{\mathrm{1}\mathrm{-}\mathit{\alpha }}}\hspace{0.17em}{\mathit{v}}_{\mathrm{esc}}\mathrm{\left(}{\mathit{R}}_{\mathrm{\ast }}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\sqrt{\mathrm{1}\mathrm{-}\frac{{\mathit{R}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathit{r}}}\mathit{,}$(42)where α is a force multiplier parameter (Lamers & Cassinelli 1999) that has been set to three different values: α = 0.465 (for T_Σ = 6 × 10³ K³); α = 0.5 (for T_Σ = 3 × 10⁴ K typical of hotspot temperatures; van Rensbergen et al. 2010b); and α = 0.640 (for T_Σ = 5 × 10⁴ K). For the sake of simplicity, the influence of the secondary star and the internal energy of the flow have been neglected.

The density of the material surrounding the system at a given time t and distance r is calculated using the continuity equation $\mathit{\rho }\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{Ṁ}\mathrm{flow}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\begin{array}{c}\mathrm{˜}\\ \mathit{\sigma }\end{array}{\mathit{r}}^{\mathrm{2}}\mathit{v}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}\mathit{,}$(43)where $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\sigma }\end{array}$ is a solid angle over which the material is ejected. The geometry of the mass ejection is badly constrained but will probably be collimated due to the small area of the hotspot. For simplicity, we use $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\sigma }\end{array}\mathrm{=}\mathit{\pi }$. The mass flow Ṁ_flow(r,t) at time t and position r is equal to the system mass-loss rate $\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}\mathit{Ṁ}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{loss}\mathit{,}\mathrm{d}\end{array}\mathrm{\left(}\mathit{t}\mathrm{\right)}$ estimated at the time such that $\mathit{r}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{t}}\mathit{v}\mathrm{\left(}\mathit{r,\tau }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\tau }$. In other words, r is the distance travelled by the ejected layer, since the beginning of the non-conservative phase.

In this simplistic model, it is also possible to follow the chemical stratification of the expelled material because we know from the stellar model the composition of the ejected layers at each time-step. The chemical changes are small until the gainer’s mass has dropped below ≈3  M_⊙, and the CNO-processed layers start to be ejected. By the end of the non-conservative evolution, the abundance of ¹⁴N has increased by a factor of 5, while ¹²C has been almost completely depleted. Since C+N+O remains constant, the global metallicity of the ejected material is almost constant.

The final density profile (when the non-conservative phase ends) is shown in Fig. 7. The x-axis represents the angle subtended by the flow for a system located at 100 pc (to fix the ideas, we remind that β Persei is situated at ≈30 pc and β Lyræ at ≈300 pc). The density drops below the mean interstellar medium (ISM) density (≈10^-24 g cm^-3, Ferrière 2001)⁴ at 12 arcmin from the star, or at 1.05 × 10¹⁸ cm. Closer to the star, the density rises up to 10^-9 g cm^-3 in the thin region of the wind acceleration. Finally, we report that the parameter α has almost no impact on the density profile in the inner region (inside the 2 arcmin region).

The calculation of the radiative transfer through this surrounding material would allow a meaningful comparison with infrared flux observations, but it is beyond the scope of this article and deferred to a future work.

4.2. Observational constraints on the spin-down mechanism

In this section, we attempt to identify observational constraints that can corroborate or invalidate spin-down mechanisms and non-conservative evolution.

4.2.1. Classification of Algol systems

Fig. 8 Evolution of the mass-transfer rate (top-panel), mass ratio (q = Md/Mg), and separation for our 6 + 3.6 M⊙ system with Pinit = 2.5 days. The different evolutionary states are coloured in red (Class I), green (Class II), and blue (Class III).

Depending on their observable properties (presence of a disc, mass ratio, and mass-transfer rate, etc.), binary systems in a given and the same evolutionary state from a theoretical standpoint may end up being associated with different observational classes. Therefore, we define theoretical classes based on our simulation results and associate a prototype with each class to ease the comparison between observations and models.

The first class (Class I, see Fig. 8) refers to the state before the mass-ratio reversal when the mass ratio and the orbital period decrease. Observational prototypes are UX Mon (Sudar et al. 2011) or SV Cen (see Table 2). During this phase, mass transfer is increasing until the system enters the rapid mass-transfer phase. These objects can host accretion discs and exhibit bipolar jet and/or show strong UV emission lines likely due to the formation of a hotspot (Guinan 1989). Taranova & Shenavrin (1997) also report excess IR emission attributed to dusty material surrounding the system.

For Class II (Fig. 8), the mass ratio is reversed and the mass-transfer rate is at a maximum. This stage is the realm of β Lyræ or W Ser objects (see Table 2). The same features as for Class I objects (discs, jets, hotspots, and circumbinary gas) may possibly be present, but Class II objects differ from Class I by their longer and increasing orbital periods. They are often referred as “active” or “hyper-active” Algols.

The last class (Class III, Fig. 8) includes quiescent objects with their mass ratio reversed. It is during this phase that the separation evolves the most, but the mass ratio remains roughly constant. This class mainly corresponds to genuine Algol systems (including the prototype β Per; see Table 2).

4.2.2. Observational constraints from Class I-II objects

Spin-down mechanisms and system-mass loss may act more efficiently during these phases (Classes I-II) and lead to clearer observational signatures, making these systems more interesting despite that they do not last long compared to Class III objects.

The observationally defined W Ser systems that we associate in our theoretical system to Class I-II objects were first detected and classified owing to their UV emission, which is attributed to the impact of the accretion stream on the edge of the accretion disc. The presence of gas around these systems (see Sect. 4.2.1) is consistent with mass ejection from a hotspot.

Accretion discs are frequently detected in Class I-II systems (Guinan 1989; van Rensbergen et al. 2011) and contribute to the regulation of the gainer’s spin velocity. The Doppler tomography method (mapping the circumstellar matter in a system) has been extensively applied to β Lyræ systems (Richards et al. 1995; Albright & Richards 1996; Richards 2004; Miller et al. 2007; Richards et al. 2010). These studies reveal the presence of disc-like structures, even in the case where the gainer has a large filling factor and where direct impact is expected instead (Barai et al. 2004). This supports our paradigm of disc formation around critically-rotating gainers and the activation of the boundary-layer mechanism (Sect. 2.4) or the disc-locking mechanism (Sect. 2.2.2), especially for Class I-II objects where the separation is short and the mass-transfer rate high.

4.3. On the rotation of Algols

Most observed Algol systems have an accretor with a rotational velocity in the range 0.1 < Ω_∗/Ω_Kep < 0.4 (van Hamme & Wilson 1990), but some of them are very rapidly rotating with a Keplerian (or close to Keplerian) surface spin velocity (McNamara 1957; Hansen & McNamara 1959; Wilson 1989; Mukherjee et al. 1996). Unfortunately, there are no observations of spin rates for the gainer in the interesting W Ser objects, because they are completely embedded in their accretion discs. The very high spin velocities in genuine Algols can be regarded as the signature of recent angular-momentum transfer or as evidence of the braking mechanism that did not have the time to yet operate. This is consistent with the boundary-layer paradigm that necessarily leads the star to critical rotation. On the other hand, another mechanism must be advocated to spin the gainer down after the rapid phase of mass transfer to explain the slow rotators. This can be due to tides or magnetic braking (with a low magnetic field strength <3 kG), which are processes that were inefficient during the earlier evolution of rapid mass transfer.

If the star is not rotating as a solid body, the outermost layers will first be spun up and, because of the differential rotation, instabilities will develop and allow angular momentum to be transported in the interior. In this case, the outermost layers will be accelerated to the critical velocity after a low amount of mass has been accreted, and we may expect the boundary-layer mechanism to be activated earlier during the evolution. If the angular-momentum redistribution inside the star is not efficient enough for the gainer to reach solid-body rotation by the end of the mass-transfer phase, less angular momentum will then be stored into the star and a larger amount given back to the orbit, thus increasing the orbital separation.

However, rotation can also affect the mechanical and thermal structure of the star, because it brakes the spherical symmetry. Considering these effects in the stellar structure equations (Endal & Sofia 1976, hereafter referred to as centrifugally supported models) leads to larger effective radii and cooler surface temperatures. Because of the short separation, the result is that our system runs into contact just before the reversal of the mass ratio. However, this contact phase can be avoided if the initial period is increased to ≈5 days (instead of 2.5 days). To reproduce the evolution of canonical models, for which the effects of departure from spherical symmetry are not accounted, our centrifugally supported models must start with a longer initial period. It also implies that neglecting these structural effects of rotation will produce less contact systems. Finally, we also point out that less room is left for an accretion disc to form due to the growth of the stellar equatorial radius, and mass loss through the third Lagrangian point ℒ₃ may be enhanced due to the growth of the stellar equatorial radius.

Finally, rotational instabilities also contribute to the transport of chemicals, and the mixing is likely to be very efficient in the gainer star because of the strong torques applied at the surface (Decressin et al. 2009). Similar to fast rotating massive stars, we may expect the accretor to evolve almost chemically homogeneously (Maeder 1987; de Mink et al. 2010), showing surface enhancements in helium and nitrogen. These alterations of the surface composition can decrease the envelope opacity, making the gainer star more compact and more luminous, thus counterbalancing the radius growth associated with the rotational deformation of the structure described above. Clearly, stellar calculations, including angular momentum transport and rotationally induced mixing, are required to clarify the situation.

5. Conclusions

In this study, we performed the first binary-evolution calculations that consistently consider the torques arising from magnetic fields and star-disc interaction. This is the first attempt to confront several braking mechanisms with the use of the state-of-the-art binary-star evolutionary code Binstar. We show that the torques due to magnetic field and star-disc interaction can prevent the accreting star from reaching a super-critical velocity, although the magnetic-field strength required are stronger than that of typical Algols. Tides are shown to be inefficient, but we do not exclude them to slow down the gainer after the mass-transfer phase. We show that some orbital key parameters of the system (e.g. the orbital separation) strongly depend on the spin-down mechanism.

Our new hotspot prescription leads to a less conservative evolution than the van Rensbergen et al. (2008) formalism. We showed that the opacity of the impacted material and the geometry of the accretion (through the $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$ parameter) are two highly-sensitive, badly-constrained key parameters for the computation of non-conservative models, stressing the need for observations and/or hydrodynamic simulations of hotspots.

A better understanding of disc formation around critically-rotating accretors and observational constraints on the magnetic-field strength of these stars are needed. In the future, we wish to establish a grid of non-conservative evolutions using Binstar to perform a statistical study, comparing these data to observational data and evaluate the impact of the various spin-down mechanisms and the hotspot formalism on the properties of the binary remnant.

Acknowledgments

We are most grateful to Walter van Rensbergen and Jean-Pierre de Grève for sharing their great expertise and for many fruitful discussions. R.D. and P.J.D. acknowledge support from the Communauté française de Belgique – Actions de Recherche Concertées. LS is an FNRS Researcher. P.J.D. is Chargé de Recherche (FNRS-F.R.S.).

References

Appendix A: Hotspots and critical mass-transfer rate

In this Appendix, we derive a limit for the gainer’s mass-accretion rate by imposing that the luminosity at the hotspot location always remains below the critical Eddington value (Eq. (35)). In the case of direct impact, the kinetic energy of the gas stream emerging from the inner Lagrangian point ℒ₁ is partly converted into thermal energy causing a hotspot and is partly radiated away.

To calculate the accretion luminosity (L_acc), we estimate the Jacobi constant (only conserved quantity in our problem) of a test particle at ℒ₁. In the rotating frame, $\begin{array}{ccc}{\mathit{C}}_{{\mathrm{ℒ}}_{\mathrm{1}}}& \mathrm{=}& \mathrm{-}\mathit{G}\left(\frac{{\mathit{M}}_{\mathrm{d}}}{{\mathit{R}}_{\mathrm{d}}}\mathrm{+}\frac{{\mathit{M}}_{\mathrm{g}}}{{\mathit{l}}_{\mathrm{g}}}\right)\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\hspace{0.17em}\left(\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{{\mathrm{ℒ}}_{\mathrm{1}}}^{\mathrm{2}}\right)\mathit{,}\end{array}$(A.1)where l_g is the distance between the centre of the gainer and ℒ₁; r_ℒ1 is the distance between ℒ₁ and the centre of mass of the system; v_ℒ1 is the initial speed of the stream at ℒ₁ in the rotating frame; and R_d is the radius of the donor filling its Roche lobe. The last term is usually negligible because, the velocity of the particle at ℒ₁ (assumed to be equal to the sound speed) is much lower than the orbital velocity (v_ℒ1 ≪ r_ℒ1Ω_orb). Similarly, we define C_imp at the impact location: $\begin{array}{ccc}{\mathit{C}}_{\mathrm{imp}}& \mathrm{=}& \mathrm{-}\mathit{G}\left(\frac{{\mathit{M}}_{\mathrm{d}}}{{\mathit{l}}_{\mathrm{d}}}\mathrm{+}\frac{{\mathit{M}}_{\mathrm{g}}}{{\mathit{R}}_{\mathrm{g}}}\right)\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{r}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\\ & \mathrm{=}& {\mathit{e}}_{\mathrm{imp}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathit{,}\end{array}$where l_d is the distance between the centre of the donor and the hotspot; v_imp is the stream velocity in the rotating frame at the impact location; and r_g is the distance between the centre of mass of the system and the hotspot. Since C_ℒ1 = C_imp, the energy gained by the stream is written as $\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{=}{\mathit{C}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{-}{\mathit{e}}_{\mathrm{imp}}$(A.4)and is partially radiated away by the hotspot and imparted to the star in the form of internal energy. Note that the stream velocity at impact v_imp is consistently computed by following the ballistic motion of the stream (see Sect. 2.5) in the rotating frame in Binstar. Before the accretion starts, the gainer’s total energy (in the rotating frame, e.g. Padmanabhan 2000) is given by ${\mathit{U}}_{\mathrm{g}}\mathrm{=}{\mathit{E}}_{\mathrm{kin}}\mathrm{+}{\mathit{E}}_{\mathrm{grav}}\mathrm{+}{\mathit{E}}_{\mathrm{i}}\mathrm{+}{\mathit{E}}_{\mathrm{nuc}}\mathrm{+}{\mathit{E}}_{\mathrm{rot}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{g}}{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\mathit{,}$(A.5)where a_g is the distance between the centre of mass of the system and the centre of mass of the gainer. After the mass deposition, ${\mathit{U}}_{\mathrm{g}}^{\mathrm{\prime }}\mathrm{=}{\mathit{E}}_{\mathrm{kin}}^{\mathrm{\prime }}\mathrm{+}{\mathit{E}}_{\mathrm{grav}}^{\mathrm{\prime }}\mathrm{+}{\mathit{E}}_{\mathrm{i}}^{\mathrm{\prime }}\mathrm{+}{\mathit{E}}_{\mathrm{nuc}}^{\mathrm{\prime }}\mathrm{+}{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{g}}^{\mathrm{\prime }}{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\mathit{,}$(A.6)where ${\mathit{E}}_{\mathrm{rot}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{I}}_{\mathrm{g}}{\mathrm{\Omega }}_{\mathrm{g}}^{\mathrm{2}}$ and E_kin are the rotational and radial motion kinetic energy of the star, and E_grav, E_i, and E_nuc are the gravitational, internal, and nuclear energy content of the star, respectively (Kippenhahn & Weigert 1990). Conservation of the total energy of the system composed of the gainer and accreted material (we neglect changes in the donor’s energy) is written as $\mathrm{(}{\mathit{U}}_{\mathrm{g}}^{\mathrm{\prime }}\mathrm{-}{\mathit{U}}_{\mathrm{g}}\mathrm{)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{acc}}{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{+}\mathrm{d}{\mathit{U}}_{\mathrm{rad}}\mathrm{+}\mathrm{d}{\mathit{U}}_{\mathrm{acc}}\mathrm{=}\mathrm{0}\mathit{,}$(A.7)where M_acc = Ṁ_acc × Δt is the accreted mass; dU_rad is the energy radiated away by the star; and dU_acc the energy lost (by radiation) during the accretion process. For a star in hydrostatic equilibrium, we have $\mathit{L}\mathrm{=}\frac{\mathrm{d}{\mathit{U}}_{\mathrm{rad}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\mathrm{(}{\mathit{E}}_{\mathrm{kin}}\mathrm{+}{\mathit{E}}_{\mathrm{grav}}\mathrm{+}{\mathit{E}}_{\mathrm{i}}\mathrm{+}{\mathit{E}}_{\mathrm{nuc}}\mathrm{)}\mathit{.}$(A.8)In the latter expression, we did not consider the rotational kinetic energy, because it is not accounted for in the stellar structure equations (there is no feed back of rotation on the star’s structure). Inserting Eqs. (A.4)−(A.6) into (A.7) with the help of (A.8) we obtain $\begin{array}{ccc}\mathrm{d}{\mathit{U}}_{\mathrm{acc}}& \mathrm{=}& {\mathit{L}}_{\mathrm{acc}}\mathrm{\times }\mathrm{\Delta }\mathit{t}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{acc}}{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{+}{\mathit{E}}_{\mathrm{rot}}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\mathrm{(}{\mathit{M}}_{\mathrm{g}}\mathrm{-}{{\mathit{M}}_{\mathrm{g}}^{\mathrm{\prime }}}^{\mathrm{)}}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{acc}}{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{+}{\mathit{E}}_{\mathrm{rot}}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{acc}}{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\mathit{.}\end{array}$In a binary system with rotating stars, the accretion luminosity is finally written as ${\mathit{L}}_{\mathrm{acc}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{\Delta }\mathit{t}}\left(\frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{acc}}{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{+}{\mathit{E}}_{\mathrm{rot}}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{M}}_{\mathrm{acc}}{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\right)\mathit{.}$(A.12)Considering that not all the accretion luminosity is radiated away⁵ (since some fraction is used to increase the internal energy), we introduce the factor α_acc < 1 to evaluate the actual luminosity associated with accretion: ${\mathit{L}}_{\mathrm{acc}\mathit{,}\mathrm{rad}}\mathrm{=}{\mathit{\alpha }}_{\mathrm{acc}}{\mathit{L}}_{\mathrm{acc}}\mathit{.}$(A.13)

The matter at the hotspot location will be ejected if the hotspot luminosity is higher than the gainer’s Eddington luminosity L_edd,g: ${\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi cG}{\mathit{M}}_{\mathit{g}}}{\mathit{\kappa }}\mathit{,}$(A.14)where κ is the plasma opacity. However, there are a few subtleties to handle properly. First, we should use the Eddington luminosity corrected for the effects of stellar rotation (Maeder 2009) for rotating stars $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}& \mathrm{=}& \mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Gamma }}^{\mathrm{)}}{\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\\ & \mathrm{=}& \end{array}$Second, the stellar luminosity contributes to the accretion luminosity in powering the mass ejection from the hotspot. However, the surfaces involved are different, and this must be considered, especially in the comparison to the Eddington luminosity, whose 4π factor in expression Eq. (A.14) implicitly assumes that the whole sphere contributes to the luminosity. For the purpose of this comparison, it is therefore necessary to multiply L_acc,rad by the factor S_star/S_hotspot to make all luminosities correspond to the same surface: ${\mathit{L}}_{\mathrm{hotspot}}\mathrm{=}{\mathit{L}}_{\mathrm{acc}\mathit{,}\mathrm{rad}}\mathrm{\times }\frac{{\mathit{S}}_{\mathrm{star}}}{{\mathit{S}}_{\mathrm{hotspot}}}\mathrm{=}{\mathit{\alpha }}_{\mathrm{acc}}\mathrm{\times }{\mathit{L}}_{\mathrm{acc}}\mathrm{\times }\frac{{\mathit{S}}_{\mathrm{star}}}{{\mathit{S}}_{\mathrm{hotspot}}}\mathit{,}$(A.17)or in other words

$\begin{array}{ccc}& & {\mathit{L}}_{\mathrm{hotspot}}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}{\mathit{L}}_{\mathrm{acc}}\\ & & \begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}\mathrm{=}{\mathit{\alpha }}_{\mathrm{acc}}\frac{{\mathit{S}}_{\mathrm{star}}}{{\mathit{S}}_{\mathrm{hotspot}}}\mathit{,}\end{array}$where S_hotspot and S_star are the hotspot and gainer’s surface areas, respectively. Values of $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{K}\end{array}$ have been calibrated by van Rensbergen et al. (2011), and we use their model $\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}\mathrm{=}\mathrm{3.9188}{\left(\frac{{\mathit{M}}_{\mathrm{d}}\mathrm{+}\mathit{M}\mathrm{g}}{{\mathit{M}}_{\mathrm{\odot }}}\right)}^{\mathrm{1.645}}\mathit{,}$(A.20)which only depends on the total mass of the system. Thus, matter at the hotspot location will be ejected if ${\mathit{L}}_{\mathrm{g}}\mathrm{+}{\mathit{L}}_{\mathrm{hotspot}}\mathit{>}\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}\mathit{,}$(A.21)where $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}$ is the gainer’s Eddington luminosity and L_g is its luminosity.

The critical rate above which the hotspot luminosity prevents matter from being accreted is given by solving for Ṁ_acc in Eq. (A.12) combined with Eqs. (A.18), (A.21), and (A.13), which then yields (Eq. (36)) $\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}\mathrm{=}\frac{\mathrm{2}}{{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}}\left(\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}\mathrm{-}{\mathit{L}}_{\mathrm{g}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}}\mathrm{+}\frac{{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}}{\mathrm{\Delta }\mathit{t}}\right)\mathrm{·}$(A.22)When the mass-transfer rate exceeds the critical value, no mass can be accreted.

van Rensbergen et al. (2008) derived a different expression for $\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}$ (their Eq. (18)), which is written using our notations as $\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}\mathrm{=}\frac{\mathrm{2}}{{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}}\left[\left(\mathrm{\Delta }\mathit{E}\frac{{\mathit{R}}_{\mathrm{g}}{\mathit{L}}_{\mathrm{edd}\mathit{,g}}}{\mathit{G}{\mathit{M}}_{\mathrm{g}}}\mathrm{-}{\mathit{L}}_{\mathrm{g}}\right)\frac{\mathrm{1}}{\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}}\mathrm{+}\frac{{\mathit{E}}_{\mathrm{rot}}^{\mathrm{\prime }}\mathrm{-}{\mathit{E}}_{\mathrm{rot}}}{\mathrm{\Delta }\mathit{t}}\right]\mathit{,}$(A.23)where $\begin{array}{ccc}\mathrm{\Delta }\mathit{E}& \mathrm{=}& \mathrm{(}\mathrm{-}{\mathit{E}}_{\mathrm{grav}}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{-}{\mathit{E}}_{\mathrm{orb}}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{-}{\mathit{E}}_{\mathrm{eq}}\mathrm{)}\\ & \mathrm{=}& \mathit{G}\left(\frac{{\mathit{M}}_{\mathrm{d}}}{{\mathit{S}}_{\mathrm{d}}}\mathrm{+}\frac{{\mathit{M}}_{\mathrm{g}}}{{\mathit{R}}_{\mathrm{g}}}\right)\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{R}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{spin}\mathit{,}\mathrm{g}}^{\mathrm{2}}\mathit{.}\end{array}$Based on our simulations of a 6 + 3.6 M_⊙ system with initial period P_init = 2.5 days, we find that E_grav ≈ 2 − 3 ×  | E_eq |  ≫ E_orb which allows us to simplify Eq. (A.23) (neglecting the E_rot terms) $\begin{array}{ccc}\mathit{Ṁ}\begin{array}{c}\mathrm{crit}\\ \mathrm{acc}\end{array}& \mathrm{\approx }& \frac{\mathrm{2}}{{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}}\left[\left(\mathrm{1}\mathrm{+}\frac{{\mathit{E}}_{\mathrm{eq}}}{{\mathit{E}}_{\mathrm{pot}}}\right){\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\mathrm{-}{\mathit{L}}_{\mathrm{g}}\right]\\ & \mathrm{\approx }& \frac{\mathrm{2}}{{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}}\left[\left(\mathrm{1}\mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{spin}\mathit{,}\mathrm{g}}^{\mathrm{2}}{\mathit{R}}_{\mathrm{g}}^{\mathrm{3}}}{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathrm{g}}}\right){\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\mathrm{-}{\mathit{L}}_{\mathrm{g}}\right]\\ & \mathrm{\approx }& \frac{\mathrm{2}}{{\mathit{v}}_{\mathrm{imp}}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{K}\end{array}}\mathrm{(}\begin{array}{c}\mathrm{˘}\\ {\mathit{L}}_{\mathrm{edd}\mathit{,}\mathrm{g}}\end{array}\mathrm{-}{\mathit{L}}_{\mathrm{g}}\mathrm{)}\mathit{.}\end{array}$This formulation is very comparable to our Eq. (36) in the case that E_rot is small, apart from some numerical factors of order of unity and the presence of ${\mathit{a}}_{\mathrm{g}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{orb}}^{\mathrm{2}}$ in our formalism.

All Tables

All Figures

	Fig. 1 Evolution of the mass-transfer rate (top), stellar masses (middle), surface velocity (bottom) for a 6 M⊙ donor (solid black line), and 3.6 M⊙ gainer (dotted red line) system with initial period Pinit = 2.5 days.
In the text

Fig. 2 Evolution of the gainer’s surface velocity (normalised to the critical velocity) for different spin-down mechanisms for our 6 + 3.6 M⊙ system with an initial period Pinit = 2.5 days. The horizontal line represents the critical spin-angular rotation of the star. All simulations are stopped when Ωsurf ≈ Ωcrit. Top panel: solid black line (labelled “Ref.”): | B∗ | = 0 kG, no braking mechanism; dotted red line (labelled “B.-L.”): boundary-layer mechanism, no magnetic field; dashed green line: Zahn (1989) tidal effect treatment only; long-dashed blue line: magnetic braking with | B∗ |  = 3 kG. Bottom panel: solid black line: magnetic braking with | B∗ | = 1 kG; dotted red line: magnetic braking with | B∗ | = 2 kG; dashed green line (labelled M.-W.): magnetic-wind braking only with | B∗ |  = 3 kG; long-dashed blue line (labelled D.-L.): disc locking only with | B∗ |  = 3 kG.

In the text

Fig. 3 Evolution of the torque on the gainer star $\mathit{J̇}\begin{array}{c}\mathrm{RLOF}\\ \mathrm{acc}\mathit{,}\mathrm{g}\end{array}$. Solid black line: no spin down mechanism (the simulation is stopped, when the gainer reaches the critical rotation). Dashed red line: boundary-layer model. The horizontal line represents the limit where the net torque spins the star up (above) or spins it down (below).

In the text

Fig. 4 Evolution of our 6 + 3.6 M⊙ system with different braking mechanisms: (solid black line) rotation-free model; (dotted red line) magnetic field 3 kG; (dashed green line) boundary-layer treatment. Top left: mass accretion rate on the gainer. The vertical line represents the epoch when the mass ratio q = 1. Top right: orbital separation. Bottom left: HR diagram for the three configurations: (solid black) rotation-free model; (dotted red) 3 kG magnetic field; (dashed blue) boundary-layer treatment; (dot-dashed green) hotspot (see Appendix A). The open circles correspond to the beginning of mass transfer, triangles to the transition between the rapid and quiescent phase of mass transfer, and squares to the end of mass transfer (only on the rotation-free model). Bottom right: extension of the donor (top cyan hatched area) and gainer (bottom magenta hatched area) within their Roche lobes RL (donor (solid blue lines) gainer (solid red lines)) for the three cases. The y-axis represents the position relative to the system’s centre of mass. The gap between the Roche lobe radii RL of the two stars is because RL is defined as the radius of a sphere with a volume equivalent to the Roche lobe and therefore differs from the true location of the ℒ1 point.

In the text

Fig. 5 Evolution of our 6 + 3.6 M⊙ system with hotspot formation. Top left: normalised hotspot Lhotspot (solid black line), accretion ${\mathit{L}}_{\mathrm{acc}}^{\mathrm{std}}$ (dashed red) and Eddington (dotted green) luminosities. Top right: evolution of the masses with (solid black) and without (dashed blue) our hotspot formalism compared to van Rensbergen et al. (2008) hotspot formalism (dotted red). Bottom left: evolution of β for the gainer using our (solid black) and van Rensbergen et al. (2008) hotspot formalism (dotted red). Bottom right: evolution of the orbital separation (same colours as top right panel).

In the text

Fig. 6 Evolution of our 6 + 3.6  M⊙ system with Pinit = 2.5 days with hotspot and different magnetic field strengths. Black: hotspot with boundary layer (BL) and no magnetic field. Red: hotspot (no BL) + magnetic field | B∗ |  = 1 kG. Green: hotspot (no BL) + magnetic field | B∗ |  = 2 kG. The red simulation stops when the gainer reaches its critical spin-angular momentum (indicated by a red circle). Top left: surface spin-angular velocity. Top right: orbital separation. Bottom left: accretion efficiency β. Bottom right: masses.

In the text

	Fig. 7 Density profile of the matter expelled by the hotspot and surrounding the system. The x-axis represents the angular extent for a system located at 100 pc. The profile starts at the stellar radius of the gainer (R∗ = 2.2 × 1011 cm) and is cut at 2 arcmin. The density drops below the mean density of the ISM value at 12 arcmin (1.05 × 1018 cm) from the star.
In the text

	Fig. 8 Evolution of the mass-transfer rate (top-panel), mass ratio (q = Md/Mg), and separation for our 6 + 3.6 M⊙ system with Pinit = 2.5 days. The different evolutionary states are coloured in red (Class I), green (Class II), and blue (Class III).
In the text