© ESO, 2014

1. Introduction

Roche lobe overflow (RLOF) gives rise to a wide variety of phenomena, with implications for many areas of astrophysics. For instance, accretion onto a massive white dwarf may lead to Type Ia supernovae (Hoyle & Fowler 1960; Wang & Han 2012, for a review), or X-ray emission from accreting neutron stars and black holes (e.g. Zeldovich & Guseynov 1966). Additionally, RLOF is a viable formation channel for blue stragglers (McCrea 1964; Geller & Mathieu 2011; Leigh et al. 2013), sub-dwarf B stars (Mengel et al. 1976; Han et al. 2002; Chen et al. 2013), and for interacting compact binaries via the common envelope phase, which is triggered by dynamically unstable RLOF from a giant or asymptotic giant branch star (Paczynski 1976; Webbink 2008).

The exchange of mass and angular momentum during RLOF causes the orbital separation to change, dictating the subsequent fate of the binary system. For a primary (i.e. the initially more massive) star of mass M₁ and a secondary star of mass M₂, the rate of change of the semi-major axis, ȧ, is determined from the rates of change of the orbital angular momentum, , the primary mass, Ṁ₁, the secondary mass, Ṁ₂, and of the eccentricity, ė, according to $\frac{\mathit{ȧ}}{\mathit{a}}\mathrm{=}\mathrm{2}\frac{\mathit{J̇}\mathrm{orb}}{{\mathit{J}}_{\mathrm{orb}}}\mathrm{-}\mathrm{2}\frac{\mathit{Ṁ}\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{-}\mathrm{2}\frac{\mathit{Ṁ}\mathrm{2}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{+}\frac{\mathit{Ṁ}\mathrm{1}\mathrm{+}\mathit{Ṁ}\mathrm{2}}{{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}\mathit{e}\mathit{ė}}{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\mathrm{·}$(1)For conservative mass transfer within circular orbits (e = 0), it is assumed that mass and orbital angular momentum are conserved, i.e. Ṁ₂ = −Ṁ₁, and (e.g. Paczyński 1971). Here, mass transfer from the more massive primary to the less massive secondary (q = M₁/M₂> 1) leads to orbital shrinkage (ȧ< 0), while the reverse occurs for q< 1 (see, e.g. Pringle & Wade 1985). However, even if the mass remains in the system, orbital angular momentum may not be conserved () because of the exchange of angular momentum between the orbit and the stellar spins via tidal torques and mass transfer (Gokhale et al. 2007; Deschamps et al. 2013).

In the 60s, several works from Piotrowski (1964), Kruszewski (1964b) and Hadjidemetriou (1969a,b) evaluated the gravitational force between the material leaving the inner-Lagrangian, ℒ₁, point (the matter stream) and each star. As a particle travels, it generates a time-varying torque on the stars, allowing for angular momentum to be exchanged between the transferred material and the orbit. Subsequently, Luk’yanov (2008) demonstrated that orbital angular momentum is conserved only if the stars are point masses, if the gravitational force between the matter stream and the stars is neglected and if the velocity of the ejected (accreted) material is equal in magnitude but in the opposite direction to the orbital velocity of the mass loser (gainer; see Appendix B). In reality, the velocity of the accreted particle is determined by the initial ejection velocity, which in turn depends on the thermal sound speed in the primary’s atmosphere and its rotation rate (Kruszewski 1964a; Flannery 1975).

It is widely assumed that tides enforce the synchronous rotation of the primary with the orbit during RLOF. However, there is observational evidence for super- and sub-synchronously rotating stars in circular binaries (e.g. Habets & Zwaan 1989; Andersen et al. 1990; Meibom et al. 2006; Yakut et al. 2007). Furthermore, Pratt & Strittmatter (1976) and Savonije (1978) argued that if the ejected material removes angular momentum faster than tides can act to synchronize the rotation, the primary star will rotate sub-synchronously and its Roche lobe radius will be affected (e.g. Limber 1963; Sepinsky et al. 2007a).

A sub-synchronously rotating primary may cause the ejected material to be accreted back onto the primary (henceforth termed “self-accretion”), which causes the orbit to shrink even when q< 1 (Sepinsky et al. 2010). Super-synchronous rotation, on the other hand, may have the opposite effect (Kruszewski 1964b; Piotrowski 1964, 1967).

In this investigation, we apply the scheme of Hadjidemetriou (1969a,b), who derived the equations of motion of a binary system using the theory of osculating orbital elements. His scheme accounts for the transfer of linear momentum between the stars, and perturbations to the orbit due to the gravitational attraction between the stars and the mass transfer stream. Using our binary stellar evolution code BINSTAR, we calculate the resulting evolution of Algol binaries for a range of initial periods and masses. Currently, we assume that none of the exchanged mass leaves the system. Torques applied onto each star by tides and mass transfer are also included.

The paper is organized as follows. In Sect. 2 we introduce the formalism of Hadjidemetriou (1969a,b). Our results are presented in Sect. 3, and discussed in Sect. 4. We summarise and conclude our investigations in Sect. 5.

2. Computational method

BINSTAR is an extension of the single-star evolution code STAREVOL. Details on the stellar input physics can be found in Siess (2010), and references therein, while the binary input physics is described in Siess et al. (2013) and Deschamps et al. (2013). In this section we present our new implementation of the osculating scheme.

2.1. Variation of the orbital parameters

Consider a binary system with an eccentricity e. The stars orbit about their common centre of mass (see Fig. 1) with an orbital angular speed, ω, and orbital period P_orb.

Material is ejected from the primary star at a rate Ṁ₁ with a velocity V₁ relative to the primary’s mass centre, given by ${{V}}_{\mathrm{1}}\mathrm{=}{{W}}_{\mathrm{1}}\mathrm{-}{{v}}_{\mathrm{1}}\mathit{,}$(2)where W₁ is the absolute velocity of the ejected material, and v₁ is the orbital velocity of the primary’s mass centre. Similarly, the velocity of the accreted material V₂ with respect to the secondary’s mass centre is ${{V}}_{\mathrm{2}}\mathrm{=}{{W}}_{\mathrm{2}}\mathrm{-}{{v}}_{\mathrm{2}}\mathit{,}$(3)where W₂ is the absolute velocity of the accreted material, and v₂ is the secondary’s orbital velocity. Ejection occurs from the ℒ₁ point, located at a distance r_ℒ1 from the primary’s centre. The Roche radius R_ℒ1 and r_ℒ1, are calculated as in Davis et al. (2013), using the formalism described by Sepinsky et al. (2007a), that accounts for the donor’s rotation.

The impact site is located at r_acc with respect to the secondary’s centre of mass. The distance r_acc is either the secondary’s radius, R₂, for direct impact accretion or the accretion disc radius. The latter is estimated from the distance of closest approach of the accretion stream to the secondary, r_min, and is given by ${\mathit{r}}_{\mathrm{acc}}\mathrm{=}\mathrm{1.7}{\mathit{r}}_{\min }$(4)(Lubow & Shu 1975; Ulrich & Burger 1976), where r_min is determined from our ballistic calculations (see Sect. 2.3). The vector r_acc forms an angle $\mathit{\psi }\mathrm{=}\mathit{\pi }\mathrm{-}\mathit{\phi ,}$(5)with the line joining the two stars, ê_r (Fig. 1).

Fig. 1 Schematic of a binary system, consisting of a primary star of mass M1, and a secondary of mass M2 orbiting with an angular speed ω, where the centre of mass of the binary system is located at . The primary and secondary are respectively located at r1 and r2 with respect to . Mass is ejected from the inner-Lagrangian, ℒ1, point (yellow star), located at rℒ1 with respect to the primary’s mass centre, . Material falls towards the secondary (dashed line) and is accreted onto its surface (or at the edge of an accretion disc) at A, situated at racc with respect to its mass centre, . Alternatively, material falls onto the primary’s surface (dotted line), landing at B, located at ${{r}}_{\mathrm{acc}}^{\mathrm{\ast }}$ with respect to . The unit vectors êr and êt point along the line joining the two stars (towards the secondary), and perpendicular to this line, respectively, and ν is the true anomaly.

In the theory of osculating elements, the rate of change of the semi-major axis, ȧ, and of the eccentricity, ė, can be expressed as (see, e.g. Sterne 1960) $\mathit{ȧ}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{orb}}}{\mathit{\pi }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}}\left[\mathrm{𝒮}\mathit{e}\sin \mathit{\nu }\mathrm{+}\mathrm{𝒯}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{\nu }\mathrm{)}\right]$(6)and $\mathit{ė}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{orb}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\frac{\mathrm{1}}{\mathrm{2}}}}{\mathrm{2}\mathit{\pi a}}\left\{\mathrm{𝒮}\sin \mathit{\nu }\mathrm{+}\mathrm{𝒯}\left[\frac{\mathrm{2}\cos \mathit{\nu }\mathrm{+}\mathit{e}\mathrm{(}\mathrm{1}\mathrm{+}{\cos }^{\mathrm{2}}\mathit{\nu }\mathrm{)}}{\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{\nu }}\right]\right\}\mathit{,}$(7)where and are the perturbing forces per unit mass, acting along ê_r, and perpendicular to that line (along ê_t), respectively, and ν is the true anomaly.

To account for the fact that the primary’s sub-synchronous rotation may cause a fraction α_self of ejected matter to be self-accreted, we decompose and into two contributions; that resulting from accretion onto the secondary, and that arising from self-accretion. If $\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}$ and $\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}$ are the mass transfer rates onto the primary and towards the secondary, respectively, then $\mathit{Ṁ}\mathrm{1}\mathit{,}\mathrm{ej}\mathrm{=}\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}\mathrm{+}\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}\mathrm{=}{\mathit{\alpha }}_{\mathrm{self}}\mathit{Ṁ}\mathrm{1}\mathit{,}\mathrm{ej}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\alpha }}_{\mathrm{self}}\mathrm{)}\mathit{Ṁ}\mathrm{1}\mathit{,}\mathrm{ej}\mathit{.}$(8)The mass accretion rate onto the secondary is $\mathit{Ṁ}{}_{\mathrm{2}\mathit{,}\mathrm{acc}}\mathrm{=}\mathrm{-}\mathit{\beta }\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}$, where β is the accretion efficiency (β = 1 for conservative mass transfer). The corresponding perturbing forces, and , are respectively $\begin{array}{ccc}{\mathrm{𝒮}}_{\mathrm{comp}}& \mathrm{=}& {\left(\frac{{\mathit{f}}_{\mathrm{2}\mathit{,r}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}\mathit{,r}}}{{\mathit{M}}_{\mathrm{1}}}\right)}_{\mathrm{comp}}\mathrm{+}\frac{\begin{array}{c}\mathrm{¨}\\ {\mathit{M}}_{\mathrm{2}\mathit{,}\mathrm{acc}}\end{array}}{{\mathit{M}}_{\mathrm{2}}}{\mathit{r}}_{\mathrm{acc}}\cos \mathit{\psi }\mathrm{-}\frac{\begin{array}{c}\mathrm{¨}\\ {\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{comp}}\end{array}}{{\mathit{M}}_{\mathrm{1}}}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\\ & & \end{array}$(9)and $\begin{array}{ccc}{\mathrm{𝒯}}_{\mathrm{comp}}& \mathrm{=}& \begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\left(\frac{{\mathit{f}}_{\mathrm{2}\mathit{,t}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}\mathit{,t}}}{{\mathit{M}}_{\mathrm{1}}}\right)}_{\mathrm{comp}}}\\ {\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}\end{array}\mathrm{+}\begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\frac{\begin{array}{c}\mathrm{¨}\\ {\mathit{M}}_{\mathrm{2}\mathit{,}\mathrm{acc}}\end{array}}{{\mathit{M}}_{\mathrm{2}}}{\mathit{r}}_{\mathrm{acc}}\sin \mathit{\psi }}\\ {\mathrm{\mathscr{H}}}_{\mathrm{comp}}\end{array}\\ & & \end{array}$(10)(Hadjidemetriou 1969b; Sepinsky et al. 2007b). The subscripts “r” and “t” refer to the radial and transverse components (i.e. along ê_r and ê_t) of the vector quantities, respectively. In Eqs. (9) and (10), terms proportional to f_i are the gravitational forces per unit mass acting on the ith star due to the mass transfer stream (see Sect. 2.4), and terms proportional to Ṁ_i are associated with the change of linear momentum for the ith star, while terms proportional to $\begin{array}{c}{\mathrm{¨}}_{}\\ {\mathit{M}}_{\mathit{i}}\end{array}$ represent the acceleration of the mass centre of the ith star arising from asymmetric mass loss or gain. In Eq. (10), ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$, and ℋ_comp correspond respectively to the gravitational force acting on the secondary by the mass transfer stream, the linear momentum transferred to the secondary and the acceleration of its mass centre, all with respect to the primary (see Sect. 3).

The self-accretion rate back onto the primary is $\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{acc}\end{array}\mathrm{=}\mathrm{-}\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}$. The associated perturbing forces, and , are found by considering the total acceleration experienced by the primary resulting from the ejection and re-capture of material (see Appendix A), giving $\begin{array}{ccc}{\mathrm{𝒮}}_{\mathrm{self}}& \mathrm{=}& {\left(\frac{{\mathit{f}}_{\mathrm{2}\mathit{,r}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}\mathit{,r}}}{{\mathit{M}}_{\mathrm{1}}}\right)}_{\mathrm{self}}\mathrm{-}\frac{\stackrel{}{{{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}_{\mathrm{¨}}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{-}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\cos {\mathit{\psi }}^{\mathrm{\ast }}\mathrm{)}\\ & & \end{array}$(11)and $\begin{array}{ccc}{\mathrm{𝒯}}_{\mathrm{self}}& \mathrm{=}& \begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\left(\frac{{\mathit{f}}_{\mathrm{2}\mathit{,t}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}\mathit{,t}}}{{\mathit{M}}_{\mathrm{1}}}\right)}_{\mathrm{self}}}\\ {\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{self}}\end{array}\begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\mathrm{-}\frac{\stackrel{}{{{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}_{\mathrm{¨}}}}{{\mathit{M}}_{\mathrm{1}}}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\sin {\mathit{\psi }}^{\mathrm{\ast }}}\\ {\mathrm{\mathscr{H}}}_{\mathrm{self}}\end{array}\\ & & \end{array}$(12)where the asterisks indicate quantities calculated at self-accretion. Here, ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{self}}$ is similar to ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$ but now for the self-accreted material, while and ℋ_self correspond respectively to the net momentum transferred to the primary by the ejected and self-accreted material, and the acceleration of its centre of mass. The radius ${\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}$ is determined from the ballistic calculations, and corresponds to the location of the particle where the Roche potential is equal to the potential at the ℒ₁-point. Using Eqs. (9) to (12), the total perturbing forces are and .

For circular orbits, ȧ is only a function of . For infinitesimal changes of a and e over one orbital period, the term in braces in Eq. (7) averages to zero over one orbit, so ė = 0 (Hadjidemetriou 1969b). Therefore, in the remainder of Sect. 2, we will just describe the quantities pertinent to the calculation of .

The total angular momentum of the binary system, J, is the sum of the spin angular momenta of each star, J_1,2, the orbital angular momentum, J_orb, and the angular momentum carried by the mass that is not attached to the stars (i.e. the mass in the wind and in the mass transfer stream), J_MT, i.e. $\mathit{J}\mathrm{=}{\mathit{J}}_{\mathrm{1}}\mathrm{+}{\mathit{J}}_{\mathrm{2}}\mathrm{+}{\mathit{J}}_{\mathrm{orb}}\mathrm{+}{\mathit{J}}_{\mathrm{MT}}\mathrm{=}\mathrm{const}\mathit{.}$(13)The rate of change of the orbital angular momentum, , is then determined by taking the time derivative of Eq. (13), and solving for , giving $\mathit{J̇}\mathrm{orb}\mathrm{=}\mathrm{-}\mathit{J̇}\mathrm{1}\mathrm{-}\mathit{J̇}\mathrm{2}\mathrm{-}\mathit{J̇}\mathrm{MT}\mathit{,}$(14)where $\begin{array}{ccc}\mathit{J̇}\mathrm{MT}& \mathrm{=}& \begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\mathrm{-}{\mathit{J}}_{\mathrm{orb}}\left[\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{\mathit{M}}\left(\frac{\mathrm{1}}{\mathit{q}}\mathrm{-}\mathit{\beta q}\right)\right]\mathrm{-}\mathit{ma}\mathrm{𝒯}}\\ \mathit{J̇}\mathrm{RLOF}\end{array}\\ & & \end{array}$(15)for circular orbits (see Appendix B). Here, m = M₁M₂/ (M₁ + M₂) is the reduced mass, Ṁ_1,2,loss< 0 is the systemic mass loss rate from each star (either via winds and/or non-conservative evolution), is the torque resulting from the material transferred between the stars, and is the torque generated by the material leaving the system (and thus associated with Ṁ_1,2,loss). The expression for is identical to that of Bonačić Marinović et al. (2008), which considers that the escaping material carries the specific orbital angular momentum of the star.

For brevity, we term our formalism the osculating scheme. If the stars are point masses, and the gravitational attraction exerted by the accretion stream is neglected, as is usually assumed, then for conservative mass transfer in Eq. (15) and we recover the classical formulation¹.

2.2. Stellar torques

The torques acting on the ith star can be decomposed into the tidal torque , and the torque arising from mass ejection or accretion, , to give $\mathit{J̇}\mathit{i}\mathrm{=}\mathit{J̇}\mathrm{tide}\mathit{,i}\mathrm{+}\mathit{J̇}\mathit{Ṁ}\mathit{i}\mathit{.}$(16)For stars with radiative envelopes, we apply the prescription for dynamical tides described by Zahn (1989). For convective stars, we use the formalism of Zahn (1977) describing equilibrium tides (see Siess et al. 2013, for further details). For we have $\mathit{J̇}\mathit{Ṁ}\mathit{i}\mathrm{=}\mathit{Ṁ}\mathit{i}\mathrm{(}\mathit{\omega }{\mathit{r}}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}_{\mathit{i,t}}\hspace{0.17em}{\mathit{r}}_{\mathit{i}}\cos \mathit{\phi }\mathrm{-}{\mathit{U}}_{\mathit{i,r}}\hspace{0.17em}{\mathit{r}}_{\mathit{i}}\sin \mathit{\phi }\mathrm{)}\mathit{,}$(17)(Piotrowski 1964; Flannery 1975) where U_i,t and U_i,r are the tangential and radial components of the ejection or accretion velocity with respect to a frame of reference co-rotating with the binary, r_i is the distance from the ith star’s mass centre to the mass ejection/accretion point (i.e. r_ℒ1 or r_acc), φ is the angle between the ejection/accretion point and the line joining the two stars, and Ṁ_i is the mass loss/accretion rate. From Eq. (17), the torque applied onto the primary because of mass ejection is $\mathit{J̇}\mathit{Ṁ}\mathrm{1}\mathrm{=}\mathrm{-}\mathrm{\left|}\mathit{Ṁ}\mathrm{1}\mathit{,}\mathrm{ej}\mathrm{\right|}\mathrm{(}\mathit{\omega }{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}_{\mathrm{1}\mathit{,t}}\hspace{0.17em}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{)}\mathit{,}$(18)where we have used the fact that φ = 0. The torque arising from self-accretion is $\mathit{J̇}\begin{array}{c}\mathrm{self}\\ \mathit{Ṁ}\mathrm{1}\end{array}\mathrm{=}\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{acc}\end{array}\mathrm{\left[}\mathit{\omega }\mathrm{\right(}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}_{\mathrm{1}\mathit{,t}}^{\mathrm{\ast }}\hspace{0.17em}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\cos {\mathit{\phi }}^{\mathrm{\ast }}\mathrm{-}{\mathit{U}}_{\mathrm{1}\mathit{,r}}^{\mathrm{\ast }}\hspace{0.17em}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\sin {\mathit{\phi }}^{\mathrm{\ast }}\mathrm{]}\mathit{,}$(19)while the torque applied onto the secondary is $\mathit{J̇}\mathit{Ṁ}\mathrm{2}\mathrm{=}\mathit{Ṁ}\mathrm{2}\mathit{,}\mathrm{acc}\mathrm{(}\mathit{\omega }{\mathit{r}}_{\mathrm{acc}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}_{\mathrm{2}\mathit{,t}}\hspace{0.17em}{\mathit{r}}_{\mathrm{acc}}\cos \mathit{\phi }\mathrm{-}{\mathit{U}}_{\mathrm{2}\mathit{,r}}\hspace{0.17em}{\mathit{r}}_{\mathrm{acc}}\sin {\mathit{\phi }}^{\mathrm{)}}\mathit{.}$(20)As shown by Packet (1981), accretion may rapidly spin up the secondary to its critical angular velocity ${\mathrm{\Omega }}_{\mathrm{2}}^{\mathrm{crit}}\mathrm{=}\mathrm{(}\mathit{G}{\mathit{M}}_{\mathrm{2}}\mathit{/}{\mathit{R}}_{\mathrm{2}}^{\mathrm{3}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. Deschamps et al. (2013) argued that super-critical rotation can be avoided either by the interaction between the secondary and an accretion disc, or magnetic braking (wind braking and disc-locking). For simplicity, we mimic these mechanisms by forcing the secondary’s spin to remain below ${\mathrm{\Omega }}_{\mathrm{2}}^{\mathrm{crit}}$. This translates into an effective torque on the gainer given by $\mathit{J̇}\begin{array}{c}\mathrm{eff}\\ \mathit{Ṁ}\mathrm{2}\end{array}\mathrm{=}\begin{array}{c}\min \end{array}\hspace{0.17em}(\frac{{\mathit{I}}_{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{2}}^{\mathrm{crit}}\mathrm{-}{\mathit{J}}_{\mathrm{2}\mathit{,}\mathrm{0}}}{\mathrm{\Delta }\mathit{t}}\mathit{,}\mathit{J̇}\mathit{Ṁ}\mathrm{2})$(21)where Δt is the evolutionary time step, which is constrained using the nuclear burning timescales, changes in the stars’ structures, the rates of change of the orbital parameters and the mass transfer rate (see Siess et al. 2013, for further details). Also, I₂ is the secondary’s moment of inertia and J_2,0 is the secondary’s angular momentum at the previous time step. For simplicity, we assume that each star rotates as a solid body, since the treatment of differential rotation is beyond the scope of this investigation (but see Sect. 4).

2.3. Ejection and accretion velocities

The components of V₁ and V₂ along ê_t, i.e. V_1,t and V_2,t, respectively are given by (Hadjidemetriou 1969b) ${\mathit{V}}_{\mathrm{1}\mathit{,t}}\mathrm{=}{\mathit{U}}_{\mathrm{1}\mathit{,t}}\mathrm{+}\mathit{\omega }{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathit{,}$(22)and ${\mathit{V}}_{\mathrm{2}\mathit{,t}}\mathrm{=}{\mathit{U}}_{\mathrm{2}\mathit{,t}}\mathrm{+}\mathit{\omega }{\mathit{r}}_{\mathrm{acc}}\cos \mathit{\psi }\mathit{.}$(23)For particles ejected from the ℒ₁ point, we set U_1,r to the sound speed, c_s, at the primary’s photosphere, and U_1,t is calculated from ${\mathit{U}}_{\mathrm{1}\mathit{,t}}\mathrm{=}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{1}}\mathrm{-}\mathit{\omega }\mathrm{)}\mathit{,}$(24)where Ω₁ is the spin angular velocity of the primary. We evaluate ψ,U_2,r and U_2,t by solving the restricted three-body equations of motion (e.g. Hadjidemetriou 1969a; Flannery 1975).

2.4. Perturbing forces from the accretion stream

We discretize the mass transfer stream as a succession of individual particles of mass δm_i. The force exerted onto the jth star by the ith particle is ${{f}}_{\mathit{ji}}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathit{j}}\hspace{0.17em}\mathit{\delta }{\mathit{m}}_{\mathit{i}}}{{\mathit{r}}_{\mathit{ji}}^{\mathrm{3}}}{{r}}_{\mathit{ji}}\mathit{,}$(25)where r_ji is the distance between the ith particle and star j. In the transverse direction, we have for the primary ${\mathit{r}}_{\mathrm{1}\mathit{i}}\mathrm{=}{\mathrm{\left\{}\mathrm{\right[}{\mathit{x}}_{\mathit{i}}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\mu }\mathrm{)}\mathit{a}{\mathrm{]}}^{\mathrm{2}}\mathrm{+}{{\mathit{y}}_{\mathit{i}}^{\mathrm{2}}}^{\mathrm{\}}}}^{\frac{\mathrm{1}}{\mathrm{2}}}$(26)and for the secondary ${\mathit{r}}_{\mathrm{2}\mathit{i}}\mathrm{=}{\mathrm{\left[}\mathrm{\right(}{\mathit{x}}_{\mathit{i}}\mathrm{-}\mathit{\mu a}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{{\mathit{y}}_{\mathit{i}}^{\mathrm{2}}}^{\mathrm{]}}}^{\frac{\mathrm{1}}{\mathrm{2}}}$(27)where μ ≡ M₁/ (M₁ + M₂) and (x_i,y_i) are the coordinates of particle i in the co-rotating frame. Summing Eq. (25) over all stream particles, and taking the difference between the gravitational force (per unit mass) acting on the secondary and on the primary star gives ${\mathrm{ℱ}}_{\mathrm{stream}}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{2}\mathit{,t}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}\mathit{,t}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{=}\mathit{G}\sum _{\mathit{i}}\mathit{\delta }{\mathit{m}}_{\mathit{i}}{\mathit{y}}_{\mathit{i}}\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{2}\mathit{i}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{1}\mathit{i}}^{\mathrm{3}}}\right)\mathrm{·}$(28)During a time-step δt of the ballistic calculation, the mass contained in the stream that falls towards the secondary is given by $\mathit{\delta }{\mathit{m}}_{\mathit{i}}^{\mathrm{comp}}\mathrm{=}\mathrm{-}\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}\mathit{\delta t}\mathit{.}$(29)Inserting Eq. (29) into Eq. (28), then in the limit δt → 0, ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$ writes as ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}\mathrm{=}\mathrm{-}\mathit{G}\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}{\mathrm{\int }}_{\mathrm{0}}^{\begin{array}{c}\mathrm{˜}\\ {\mathit{t}}_{\mathrm{comp}}\end{array}}{\mathit{y}}_{\mathrm{comp}}{\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{2}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{1}}^{\mathrm{3}}}\right)}_{\mathrm{comp}}\hspace{0.17em}\mathrm{d}\mathit{t}\mathit{.}$(30)Similarly for the stream falling back onto the primary, we have $\mathit{\delta }{\mathit{m}}_{\mathit{i}}^{\mathrm{self}}\mathrm{=}\mathrm{-}\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}\mathit{\delta t}$(31)giving $\begin{array}{ccc}{\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{self}}\mathrm{=}\mathrm{-}\mathit{G}\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}{\mathrm{\int }}_{\mathrm{0}}^{\begin{array}{c}\mathrm{˜}\\ {\mathit{t}}_{\mathrm{self}}\end{array}}{\mathit{y}}_{\mathrm{self}}{\left(\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{2}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{1}}^{\mathrm{3}}}\right)}_{\mathrm{self}}\hspace{0.17em}\mathrm{d}\mathit{t}\mathit{.}& & \end{array}$(32)In Eqs. (30) and (32), $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{t}\end{array}$ is the particle’s travel time between the ℒ₁ point and the point of impact, and r₁, r₂, and y describe the position of the particle at time t, given by our integration of the ballistic trajectories. The subscripts “comp” and “self” indicate quantities pertaining to material falling onto the secondary and onto the primary, respectively.

2.5. Treatment of the mass transfer stream

The finite width of the matter stream that is ejected from the ℒ₁-point is accounted for when calculating the quantities found in Eqs. (10) and (12) as described below. For a primary with an effective temperature T_eff,1, mean molecular weight at the photosphere μ_ph,1, the surface area, , of the stream at the ℒ₁ point is $\mathrm{𝒜}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }\mathrm{ℛ}{\mathit{T}}_{\mathrm{eff}\mathit{,}\mathrm{1}}}{{\mathit{\mu }}_{\mathrm{ph}\mathit{,}\mathrm{1}}}\frac{{\mathit{a}}^{\mathrm{3}}\mathit{q}}{\mathit{G}{\mathit{M}}_{\mathrm{1}}}{\left\{\mathit{g}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{\left[}\mathit{g}\mathrm{\right(}\mathit{q}\mathrm{)}\mathrm{-}\mathit{f}\mathrm{-}\mathit{qf}\mathrm{]}\right\}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{,}$(33)(Davis et al. 2013) where ℛ is the ideal gas constant. In Eq. (33), $\mathit{g}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\frac{\mathit{q}}{\mathrm{(}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathit{/}\mathit{a}{\mathrm{)}}^{\mathrm{3}}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathit{/}\mathit{a}\mathrm{)}{\mathrm{]}}^{\mathrm{3}}}\mathit{,}$(34)(Kolb & Ritter 1990) and $\mathit{f}\mathrm{=}\frac{{\mathrm{\Omega }}_{\mathrm{1}}}{\mathit{\omega }}$(35)corrects for any effects arising from the primary’s rotation.

Fig. 2 Evolution in time (since the start of mass transfer at tRLOF = 9.62034 × 107 yr) of the primary’s spin angular velocity in units of the orbital velocity, Ω1/ω (dotted green curve, left axis), for the 5 + 3 M⊙ system, initial period Pi = 7 d, using the osculating formalism. Panel a): as the primary evolves off the main sequence up to the onset of mass transfer; panel b): during the self-accretion phase; panel c): during the slow mass transfer phase. Panel b) also shows the fraction of ejected material that falls back onto the primary, αself (solid black curve, left axis), and the ratio of the tidal synchronization timescale to the mass transfer timescale, τsync/τṀ (dashed blue curve, right axis). Note the change in scales along the x-axis in each panel. The top panels indicate the evolution of q = M1/M2, and the open red square indicates where q = 1.

We then calculate the ballistic trajectory of N particles, which are ejected at evenly spaced intervals, Δℓ, along the width of the ℒ₁-point, $\mathit{\ell }\mathrm{\approx }\sqrt{\mathrm{𝒜}}$. To account for the Gaussian density distribution of the particles at the ℒ₁ nozzle (see, e.g. Lubow & Shu 1975; Edwards & Pringle 1987; Raymer 2012) we weight each trajectory using $\mathit{\xi }\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ \mathit{\ell }\end{array}\mathrm{\right)}\mathrm{=}\mathit{\eta }{\mathrm{e}}^{\mathrm{-}\frac{\stackrel{}{{{\mathit{\ell }}^{\mathrm{2}}}_{\mathrm{˜}}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\mathit{,}$(36)where $\mathrm{-}\mathrm{1}\mathrm{\le }\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\ell }\end{array}\mathrm{\le }\mathrm{1}$ is the normalised distance along the finite width of the stream (the ℒ₁-point is located at $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\ell }\end{array}\mathrm{=}\mathrm{0}$), σ is the standard deviation and η is a normalisation constant set by the requirement that ${{}^{\mathrm{\int }}}_{-1}^{\mathrm{1}}\mathit{\xi }\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\ell }\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\ell }\end{array}\mathrm{=}\mathrm{1}$. We use σ = 0.4 so that the density at $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\ell }\end{array}\mathrm{=}\mathrm{-}\mathrm{1}$ and $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{\ell }\end{array}\mathrm{=}\mathrm{1}$ is equal to the donor’s photospheric density, and use N = 128 particles. We checked that increasing N any further or slightly varying σ has negligible impact on the calculations.

Since we are dealing with more than one particle landing on each star, we calculate the average stream properties as follows. Let ${\mathit{Q}}_{\mathit{k}}^{\mathrm{comp}}\mathrm{\in }\mathrm{\left\{}{\mathrm{ℱ}}_{\mathrm{stream}\mathit{,k}}^{\mathrm{comp}}\mathit{,}{\mathit{U}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,t,k}}\mathit{,}{\mathit{r}}_{\mathrm{acc}\mathit{,k}}\mathit{,}{\mathit{\psi }}_{\mathit{k}}\mathrm{\right\}}$ be some quantity related to the kth particle that is ejected from $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\ell }}_{\mathit{k}}\end{array}$, which is subsequently accreted by the companion. Then the mean value of this quantity for a given model is calculated using $\mathrm{⟨}{\mathit{Q}}^{\mathrm{comp}}\mathrm{⟩}\mathrm{=}\frac{{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{comp}}}{\mathit{Q}}_{\mathit{k}}^{\mathrm{comp}}\mathit{\xi }\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ {\mathit{\ell }}_{\mathit{k}}\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{\Delta }\begin{array}{c}\mathrm{˜}\\ \mathit{\ell }\end{array}}{{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{comp}}}\mathit{\xi }\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ {\mathit{\ell }}_{\mathit{k}}\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{\Delta }\begin{array}{c}\mathrm{˜}\\ \mathit{\ell }\end{array}}\mathit{,}$(37)where N_comp is the number of particles landing on the companion. A similar expression for the particles landing on the donor is used, where we replace “comp” with “self”. Hence, we replace ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{self}\mathit{,}\mathrm{comp}}$, U_1,2,t, ${\mathit{U}}_{\mathrm{1}\mathit{,t}}^{\mathrm{\ast }}$, ${\mathit{r}}_{\mathrm{acc}}^{\mathrm{(}\mathrm{\ast }\mathrm{)}}$ and ψ^(∗) in Eqs. (10) and (12) with the corresponding mean values, as determined by Eq. (37). Finally, the fraction of particles landing back on the donor star is ${\mathit{\alpha }}_{\mathrm{self}}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{self}}}\mathit{\xi }\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ {\mathit{\ell }}_{\mathit{k}}\end{array}\mathrm{\right)}\hspace{0.17em}\mathrm{\Delta }\begin{array}{c}\mathrm{˜}\\ \mathit{\ell }\end{array}\mathit{,}$(38)where N_self is the number of particles undergoing self-accretion (N = N_self + N_comp).

Fig. 3 Left: evolution for the 5 + 3 M⊙ binary system with Pi = 7 d of (solid black curve), (long-dashed green) and ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}\mathit{,}\mathrm{self}}$ (short-dashed red) arising from mass transfer onto the companion (panel a)) and self-accretion (panel b)). Panel c) compares (long-dashed magenta), (solid cyan) and the total (short-dashed blue). Right: the same but for the 6 + 4 M⊙ system, Pi = 2.5 d.

2.6. The binary models

To analyse the impact of this new formalism on the evolution of Algols, we consider two systems. The first configuration is a 5 + 3 M_⊙ binary, with an initial period P_i = 7 days, undergoing case B mass transfer (during shell H-burning). Our second system is a 6 + 4 M_⊙ binary, P_i = 2.5 days, which commences mass exchange during core H-burning (case A). We assume a solar composition (Z = 0.02), and apply moderate convective core overshooting, with α_ov = Λ /H_p = 0.2, where Λ is the mixing length and H_p is the pressure scale height. The initial spin periods of each star are set to the initial orbital periods at the start of the simulation (i.e. f = 1 in Eq. (35)) and the orbits are circular.

3. Results

Section 3.1 presents our osculating calculations for the case B system, which are compared to calculations determined from the classical scheme (Sect. 3.1.1), followed by case BB (shell He-burning) mass transfer in Sect. 3.1.2. Calculations for the case A model are given in Sect. 3.2.

3.1. 5 + 3 M_⊙, P_i = 7 days

3.1.1. Case B mass transfer

As the primary evolves off the main sequence, its radius and moment of inertia increase on a much shorter timescale (τ_R = R/Ṙ ≈ 10⁸ yr) than tides can maintain synchronization ( yr) and Ω₁ declines (Fig. 2a). Therefore, mass transfer starts while the primary star is rotating significantly sub-synchronously, and this has two important consequences. First, the Roche radius calculated using the Sepinsky et al. (2007a) formalism is 7 per cent larger than the Eggleton (1983) prescription. Consequently, the primary must evolve further along the sub-giant branch before it can start transferring mass. Secondly, at the onset of RLOF, all the material ejected by the primary star is initially self-accreted (α_self = 1, panel b, solid black curve). This process induces a positive torque onto the primary (Eq. (19)), increasing Ω₁ (dotted green curve) and thus U_1,t, causing material to progressively flow onto the companion (α_self< 1). In parallel, mass ejected through the ℒ₁-point applies a negative torque onto the primary (Eq. (18)). Since, initially, more material is falling onto the primary than onto the companion, the net effect is an acceleration of the primary’s rotation rate, i.e. $\mathrm{\left|}\mathit{J̇}{}_{\mathit{Ṁ}\mathrm{1}}\mathrm{\right|}\mathit{<}\mathit{J̇}\begin{array}{c}\mathrm{self}\\ \mathit{Ṁ}{}_{\mathrm{1}}\end{array}$, and a reduction in the amount of self-accreted material. As Ω₁ increases, $\mathit{J̇}\begin{array}{c}\mathrm{self}\\ \mathit{Ṁ}{}_{\mathrm{1}}\end{array}$ drops until eventually the net torque applied onto the primary via the ejection and self-accretion process is zero, i.e. $\mathrm{\left|}\mathit{J̇}{}_{\mathit{Ṁ}\mathrm{1}}\mathrm{\right|}\mathrm{\approx }\mathit{J̇}\begin{array}{c}\mathrm{self}\\ \mathit{Ṁ}{}_{\mathrm{1}}\end{array}$. This situation is characterized by a plateau in α_self ≈ 0.2. Throughout the self-accretion process, the mass loss timescale τ_Ṁ = M₁/ | Ṁ₁ | <τ_sync, so tides are unable to enforce synchronous rotation (dashed, blue curve).

For circular orbits, Eq. (6) reduces to $\mathit{ȧ}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{orb}}}{\mathit{\pi }}\mathrm{𝒯}\mathrm{=}\frac{{\mathit{P}}_{\mathrm{orb}}}{\mathit{\pi }}\mathrm{(}{\mathrm{𝒯}}_{\mathrm{comp}}\mathrm{+}{\mathrm{𝒯}}_{\mathrm{self}}\mathrm{)}\mathit{,}$(39)and so the sign of ȧ depends on the signs of and . To understand the evolution of the separation, we show in Fig. 3 the various contributions entering the expression for , i.e. the gravitational force applied onto the secondary (primary) by the ejected particle ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}\mathit{,}\mathrm{\left(}\mathrm{self}\mathrm{\right)}}$, and the linear momentum transferred to the companion (donor), , all with respect to the donor. Throughout mass transfer, the terms ℋ_self and ℋ_comp, corresponding to the acceleration of the donor’s or companion’s mass centre, are negligible because $\begin{array}{c}{}_{\mathrm{¨}}\\ {\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}\end{array}\mathrm{\approx }\mathrm{0}$, and for this reason it is not displayed in Fig. 3. At the start of mass transfer, the contribution to partially comes from ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$ (Fig. 3a, left panels, red short-dashed curves), which is negative for two reasons. Firstly, the particles are located at y_comp< 0 and secondly they are typically situated in the vicinity of the secondary i.e. 1 /r₂> 1 /r₁ (see Eq. (30)).

The term (green, long-dashed curves) is small for q> 1, but increases as q declines. Indeed, as mass transfer proceeds, the primary’s Roche lobe radius and therefore r_ℒ1 shrink, and so the ℒ₁ point moves further away from the companion’s surface. A particle’s travel time thus rises and it can accelerate to a larger impact velocity (V_2,t). Even though also depends on V_1,t + ωr_ℒ1, its magnitude is about a factor of 10 smaller than the | V_2,t + ωr_acccosψ | term.

During self-accretion, is dictated by (Fig. 3b). We find that V_1,t ≈ − 1 × 10⁷ cm s^-1 and ${\mathit{V}}_{\mathrm{1}\mathit{,t}}^{\mathrm{\ast }}\mathrm{\approx }\mathrm{8}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}$ cm s^-1, which, for small ψ^∗, gives (Eq. (12)). Similarly to ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$, ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{self}}\mathit{<}\mathrm{0}$ because y_self< 0 and, owing to the large value of r_ℒ1 for q> 1, the particle is closer to the secondary than to the primary.

Hence, self-accretion enhances the rate of shrinkage of both the orbital separation (Fig. 3c) and the primary’s Roche radius. This faster contraction of R_ℒ1 increases the overfilling factor resulting in a mass transfer rate that is a factor of ~1.6 higher than in the classical scheme. In response to this higher mass loss rate, the primary’s radiative layers further contract and it attains a lower surface luminosity on the Hertzsprung-Russel (HR) diagram compared to the classical calculation (Fig. 4a).

To quantify the impact of self-accretion, we re-ran a simulation by setting , as indicated by the short-dashed green curve in Fig. 4d. Differences in the mass transfer rate and the evolution along the HR diagram are negligible. However, the post-mass transfer orbital period (about 21 d) is slightly longer than when self-accretion is included (17 d); a relative difference of about 20 per cent.

Fig. 4 Evolution of the 5 + 3 M⊙ system, Pi = 7 d, calculated using the osculating (solid black curves) and the classical (dotted red curves) prescriptions. The short-dashed green and the long-dashed blue curves also use the osculating scheme, but and ℱstream have been set to zero, respectively. For clarity, the dashed green curve has been omitted from panels a) to c) as it is indistinguishable from the black curve. The long-dashed blue curve has been truncated in panel a), since it follows the same track as the black and red curves during core-He burning. Panels a) and b) show the evolutionary path in the Hertzsprung-Russel diagram of the primary and secondary. Panel c) shows the mass loss rate, Ṁ1, as a function of time since mass transfer started and panel d) the evolution of the orbital period, Porb, as a function of the mass ratio, q = M1/M2. The different evolutionary phases are indicated by the open green squares, and are labelled as follows: A: start of case B mass transfer; B: q = 1; C: end of case B mass transfer (coincides with the end of H-shell burning); D: start of case BB mass transfer (onset of He-core burning); E: end of case BB mass transfer (onset of He-shell burning).

As mass transfer proceeds, τ_sync decreases as a result of the deepening of the primary’s surface convection zone in response to mass loss (Webbink 1977a). Eventually, τ_sync<τ_Ṁ (Fig. 2b), and tides can re-synchronize the primary. Moreover, as the primary’s mass declines, it exerts a smaller gravitational attraction onto the ejected particle, and evermore material falls onto the secondary star, as indicated by the decrease in α_self.

Once q ≲ 1, self-accretion shuts off and soon after (when q ≃ 0.9) becomes positive, allowing the orbit to expand. The rise in is because of the aforementioned growth of the term, resulting from the higher impact velocity (V_2,t). Eventually, the primary restores thermal equilibrium, Ṁ₁ declines (Fig. 4c) and the evolution enters the slow phase (around q ≈ 0.16), where mass transfer occurs on the nuclear timescale of the hydrogen-burning shell (Paczyński 1971). The calculated mass transfer rates for the classical and osculating models are virtually identical, since the shell-burning properties in both cases are the same². Also note that throughout the slow phase, τ_sync ≪ τ_Ṁ and so tides can enforce synchronous rotation of the primary (Fig. 2c).

By the end of the simulations, the difference in the periods is significant with 17 days for the osculating scheme compared to 71 days for the classical model. The primary’s radius is correspondingly smaller, since it keeps track of its Roche lobe, which is a function of a. This explains why, for a given luminosity, the osculating model gives a hotter primary than in the classical case (Fig. 4a).

The shorter orbital period predicted by the osculating scheme is caused by the negative ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$ contribution. Indeed, neglecting this term gives , yielding the longest post-mass transfer orbital periods out of all our calculations (Fig. 4d, long-dashed blue curve). In this case, the Roche filling primary has a larger radius and a lower effective temperature, causing a substantial surface convection zone to develop (radial extent of about 60 R_⊙). This enhances the mass transfer rate because convective layers expand upon mass loss, causing the star to over-fill its Roche lobe further.

When He ignites in the primary, mass transfer terminates and the final masses are virtually identical between the classical and osculating schemes (M₁ ≈ 0.87 M_⊙, M₂ ≈ 7.3 M_⊙).

Both the osculating and classical formalisms predict similar evolutionary tracks of the secondary on the HR diagram (Fig. 4b). Once mass transfer enters the slow phase, the secondary relaxes towards thermal equilibrium, establishing a new effective temperature and luminosity along the main sequence which is appropriate for its new mass.

3.1.2. Case BB mass transfer

Once mass transfer has stopped, the structure of the primary consists of a convective He-burning core of 0.27 M_⊙, surrounded by a radiative envelope of ~ 0.60 M_⊙. The H-burning shell is located at mass coordinate M_r ≈ 0.73 M_⊙, and has a mass of about 0.06 M_⊙. When He ignites in the core, the primary contracts within its Roche lobe on a timescale much shorter than τ_sync leading to super-synchronous rotation (Fig. 5, left panel, solid green curve). The star accelerates up to approximately 3 per cent of the critical velocity when Ω₁/ω ≈ 12. These calculations therefore predict the presence of rapidly rotating core-He burning stars in detached binaries. To the best of our knowledge, there are no available observations of such systems in this evolutionary phase.

At about 2.72 × 10⁷ yr since the start of case B mass transfer, the activation of shell-He burning produces a rapid expansion of the primary. As the star fills more of its Roche lobe (R₁/R_ℒ,1, dashed cyan curve) the tidal forces strengthen and within ~ 1.5 × 10⁶ yr, the primary is re-synchronized. The ensuing case BB (Fig. 5, shaded region) is characterized by a constant mass exchange rate of about 10^-7M_⊙ yr^-1, and an expansion of the orbit because q< 1. Synchronous rotation is maintained during the entire phase of mass transfer, so no self-accretion occurs. In both schemes, mass transfer lasts for about 3 × 10⁵ yr.

Mass transfer ceases as a result of re-ignition of the H-burning shell, which causes the primary’s radius to shrink. At this point, the binary consists of a 0.8 M_⊙ CO star, undergoing H-shell burning in the surface layers, and a 7.2 M_⊙ main sequence companion. The final orbital periods are approximately 80 days and 19 days for the classical and osculating calculations, respectively.

Fig. 5 Similar to Fig. 2, but now during core He burning (left panel), and shell He-burning (right panel), the onset of which is indicated by the arrow. The dashed cyan curve shows R1/Rℒ,1 (right axis), and the shaded region marks case BB mass transfer.

To calculate the subsequent evolution of this system is beyond the scope of the investigation. Nonetheless, we can infer its fate based on the binary parameters. Eventually, the secondary star will evolve off the main sequence, fill its own Roche lobe and transfer mass back to the CO primary star. For the osculating scheme, we estimate that the secondary will fill its Roche lobe with a radius of about 34 R_⊙ as it is crossing the sub-giant branch, and for the classical scheme, when the star approaches the base of the giant branch with a radius of ~ 60 R_⊙. Because of the high mass ratio M₂/M₁ ≈ 8, which lies well above the critical value of between 1.2 and 1.3 for dynamically unstable mass transfer (Webbink 2008), we therefore expect this system to enter common envelope evolution.

Fig. 6 a) Evolution of αself as a function of time, tself, since the start of self-accretion for a 5 + 3 M⊙ system, with an initial orbital period, Pi of 3 d (dotted red curve), 7 d (solid black), 10 d (green short dashed) and 13 d (blue long-dashed). b) Evolution of the orbital period, Porb, as a function of q = M1/M2, for different initial orbital periods: 3 d (dot-dashed); 7 d (solid); 10 d (short-dashed) and 13 d (dotted). Black curves refer to classical calculations, and the red lines to the osculating scheme.

3.1.3. Effect of changing the initial orbital period

Figure 6a shows that α_self levels off to between 0.15 and 0.2 for the considered values of P_i. In addition, the duration of the self-accretion phase progressively decreases, from 5.2 × 10⁴ to 2.6 × 10⁴ yr when the initial period, P_i, rises from 7 to 13 d. By increasing P_i, the primary evolves further along the sub-giant branch before mass transfer starts, and can develop a deeper convection zone. When RLOF occurs, tides are more efficient at re-synchronizing the primary, and self-accretion is stopped earlier.

For P_i between 7 and 13 days, α_self ≈ 1 at the start of mass transfer. By contrast, for P_i = 3 d, α_self rises from zero to a constant value of approximately 0.2. In this case, the primary is still very close to its main sequence location in the HR diagram, and it is initially not rotating sufficiently slowly to trigger self-accretion. Only once enough angular momentum has been removed from the primary does self-accretion occur.

The final primary masses are between 0.85 and 0.88 M_⊙ and the secondary masses between 7.15 and 7.12 M_⊙, irrespective of whether the classical or the osculating scheme is used. However, the osculating scheme systematically yields shorter post-mass transfer orbital periods, by a factor of about 4 (Fig. 6b), and for the model with P_i = 3 d, it gives rise to a contact system, in contrast to the classical scheme.

3.2. 6 + 4 M_⊙, P_i = 2.5 d

Since the initial period is shorter for this model, tides are able to establish synchronous rotation by the time the primary fills its Roche lobe (Fig. 7a). Therefore, all ejected mass is initially accreted onto the companion (α_self = 0) and only once mass ejection has removed a sufficient amount of angular momentum from the primary does self-accretion occur, with α_self ≈ 0.15.

As we remarked, self-accretion enhances the orbital contraction because of the negative contribution from (right panel, Fig. 3b). In contrast to the case B model, this process continues after the mass ratio has reversed and only ceases when q ≈ 0.6. The reason for this difference stems from the time delay associated with the appearance of a surface convection zone, which reinforces the tidal interaction and accelerates the primary’s rotation velocity back to synchronous rotation. This persisting self-accretion episode maintains for a longer period of time and the orbit contracts until q ≈ 0.8 (Fig. 3c). Also note that ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{self}}\mathit{>}\mathrm{0}$ when q< 1 because the ejected material is typically located in the vicinity of the primary (1 /r₁> 1 /r₂ in Eq. (32)), owing to the close proximity of the ℒ₁-point to the primary. The positive value for ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$, on the other hand, is because the primary’s sub-synchronous rotation deflects the particles such that they have y_comp> 0.

Fig. 7 Similar to Fig. 2, but now for the 6 + 4 M⊙ system, Pi = 2.5 days. Here, tRLOF = 6.12305 × 107 yr.

Fig. 8 Similar to Fig. 4 but for the 6 + 4 M⊙, Pi = 2.5 d system. A: start of case A mass transfer; B: q = 1; C: end of core H-burning; D: start of shell H-burning; E: end of case A mass transfer (H-shell burning ceases).

At t ≈ t_RLOF + 1 × 10⁶ yr (q ≈ 0.3), the convection zone recedes and τ_sync rises again (Fig. 7b). Subsequent mass ejection brings the primary back into sub-synchronous rotation, triggering a second self-accretion episode, with about 10 per cent of the ejected material falling back onto the star. This second self-accretion event is absent in the case B model because of the higher sound speed (U_1,r) for the primary. The impact of this second occurrence on the orbital evolution, however, is negligible. This is because $\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}$ and ${\mathit{V}}_{\mathrm{1}\mathit{,t}}\mathrm{-}{\mathit{V}}_{\mathrm{1}\mathit{,t}}^{\mathrm{\ast }}$ are respectively a factor of about 10 and 6 smaller than the values for the first self-accretion episode. Neglecting (short-dashed green curve, Fig. 8d) shows that the relative difference between the post-mass transfer orbital periods with and without self-accretion is about 15 per cent.

As the mass transfer rate decelerates, τ_sync/τ_Ṁ correspondingly declines, Ω₁ increases back towards synchronous rotation and self-accretion shuts off. Eventually, hydrogen is exhausted in the core, the primary shrinks within its Roche lobe, and mass transfer is terminated (giving rise to the hook feature around log ₁₀(L/L_⊙) ≃ 2.5 in Fig. 8a). Mass exchange resumes once the primary re-expands because of H-shell ignition. Since now τ_sync/τ_Ṁ ≪ 1, synchronous rotation is maintained throughout the subsequent mass transfer episode which starts at t = t_RLOF + 5 × 10⁶ yr (Fig. 7c).

Upon He-core ignition, the final mass of the osculating primary is marginally less massive (≈ 0.9 M_⊙) than in the classical models (1.0 M_⊙), the total system mass being held constant. As with the case B model, however, the osculating scheme gives a final orbital period which is a factor of 5.5 smaller (≈6 d) than the classical scheme (≈33 d, Fig. 8d).

As for the case B model, ${\mathrm{ℱ}}_{\mathrm{stream}}^{\mathrm{comp}}$ is responsible for the shorter post-mass transfer orbital period. The long-dashed blue curve in Fig. 8 presents the evolution if the ℱ_stream terms are neglected. When q> 1, the orbital separation slightly increases (Fig. 8d) because of the dominant, positive contribution from the ωr_acccosψ term in Eq. (10). The expansion of the orbit reduces the amount that the primary star over-fills its Roche lobe, giving a correspondingly smaller mass loss rate peaking at about 2 × 10^-7M_⊙ yr^-1 (Fig. 8c). For a given t − t_RLOF, the primary is therefore more massive than found by the other models, and its core-hydrogen burning timescale is, in turn, shorter. By consequence, both shell-hydrogen burning and core-helium burning commence sooner in the binary’s evolution and so the duration of mass transfer is shorter.

Eventually, the secondary will fill its Roche lobe as it evolves off the main sequence, before He-core burning has ceased in the primary. This second mass transfer episode starts when the secondary’s radius is 17 R_⊙ for the osculating scheme and 53 R_⊙ for the classical model. As for the case B system, we expect common envelope evolution to follow because of the extreme mass ratio.

4. Discussion

4.1. Consequences of the osculating scheme on the orbital evolution

Our simulations show that the osculating scheme yields a significantly shorter post-mass transfer orbital period than the classical formalism. Alternatively expressed, to obtain an Algol with a given orbital period, the osculating prescription requires a longer initial period. Consequently, the progenitor primary may fill its Roche lobe once it has already developed a deep surface convection zone near the base of the giant branch. Tout & Eggleton (1988a) suggested this was the case for a number of observed Algols, for example TW Dra and AR Mon. Mass transfer would therefore proceed on the dynamical – rather than the thermal – timescale, possibly causing common envelope evolution. They proposed, however, that such a fate can be avoided if the primary loses sufficient mass via an enhanced stellar wind (companion-reinforced attrition process; Tout & Eggleton 1988b) to reduce the mass ratio close to unity, before RLOF starts.

We currently consider conservative evolution, which may be a reasonable assumption during the slow mass transfer phase. van Rensbergen et al. (2008) suggested that non-conservative mass transfer in Algols is triggered by a hotspot located at the secondary’s surface, or at the edge of an accretion disc. They further found that this mechanism typically operates during the rapid mass transfer phase, and becomes quiescent within the slow regime. Other suggested mechanisms – albeit poorly studied – include mass escaping through the outer Lagrangian ℒ₃ point (Sytov et al. 2007) or by bipolar jets (Ak et al. 2007). The removal of orbital angular momentum via systemic mass losses was also invoked by De Greve et al. (1985) and De Greve & Linnell (1994) to explain the observed orbital periods of TV Cas, β Lyrae and SV Cen. For TV Cas, for example, the authors estimated that 80 per cent of the transferred mass is ejected from the system, removing about 40 per cent of the progenitor’s orbital angular momentum. However, given the effect of ℱ_stream on the orbital evolution, we argue that a fraction of the estimated orbital angular momentum loss can be attributed to the gravitational interaction between the stars and the mass transfer stream. Therefore, the amount of orbital angular momentum carried by the expelled mass may be lower than quoted by De Greve et al. (1985). We will consider the impact of non-conservative evolution in a future investigation.

Our simulations indicate that, contrarily to what is usually assumed, the primary does not always rotate synchronously throughout mass exchange. As shown in Figs. 2 and 7, the primary is rotating significantly sub-synchronously during the rapid phase. Only once the mass transfer rate decelerates and convection develops in the surface layers, are tides effective enough to re-synchronize the primary.

Unfortunately, published spin rates of the primary during the rapid phase are, to the best of our knowledge, not available, but this is likely a result of the short duration of this phase (between about 10⁵ to 10⁶ years in our models). Existing studies of rapid mass transfer systems, such as UX Mon (Sudar et al. 2011), assume that the primary is synchronous, based upon the expectations that tides are always efficient enough to maintain synchronism, although this has never been proven observationally. On the other hand, evidence for synchronous primaries where the mass ratio has reversed are relatively abundant, such as the eclipsing δ Scuti star KIC 10661783 (Lehmann et al. 2013), TX Uma (Glazunova et al. 2011), KZ Pav (Sürgit et al. 2010) and RX Cas (Andersen et al. 1989).

4.2. Contact evolution

The enhanced orbital contraction that self-accretion generates, combined with expanding radiative secondaries in response to mass accretion (e.g. Neo et al. 1977), could lead to more contact systems. The situation may be particularly severe for case A systems when self-accretion operates even after the mass ratio has reversed. Indeed, even though the classical scheme predicts orbital expansion when q< 1, self-accretion still causes the orbit to contract down to q ≃ 0.8 in our case A simulation. However, contact may be avoided if, for a given initial period, the mass ratio is initially close to unity. In this configuration, the mass ratio is reversed sooner, limiting the orbital shrinkage.

The spin-up of the secondary, as it accretes angular momentum from the transferred material (Packet 1981), may also lead to a contact binary. Sepinsky et al. (2007a) demonstrated that the Roche lobe radius of a super-synchronously rotating star will be smaller than the value determined using the standard Eggleton (1983) formula. So, if we also consider the impact of the secondary’s rotation on its Roche radius (as should be done), we indeed find that all our models enter a contact configuration during the rapid phase. However, magnetic braking may prevent significant spin-up of the secondary (Dervişoǧlu et al. 2010; Deschamps et al. 2013) although, in light of the expanding radiative secondary, these mechanisms may only delay the onset of contact.

Contact evolution requires detailed modelling of energy transport between the stars within their common envelope (e.g. Webbink 1977b; Stȩpień 2009), which is beyond the scope of the present paper. Nonetheless, observations of the contact systems with early spectral-type stars, such as LY Aur (Zhao et al. 2014), V382 Cyg and TU Mus (Qian et al. 2007) indicate that their evolution is similar to that of a semi-detached Algol, but much shorter-lived. The authors suggest that the contact configuration was most likely triggered during a rapid case A mass transfer phase, and that the observed period increase is caused by mass transfer from the less massive to the more massive star. They also expect that the increasing orbital separation will break the contact configuration, giving a semi-detached Algol, suggesting that not all contact systems necessarily merge.

Further complications arise from the fact that the rotation rate of each star is different. The concept of the Roche model, which assumes that both stars rotate synchronously with the orbit, is therefore incorrect. As outlined by Vanbeveren (1977), asynchronous rotation of both components gives Roche lobes that do not necessarily coincide at a common inner-Lagrangian point. This situation may greatly complicate the mass flow structure between the stars, possibly leading to non-conservative evolution. We appeal to smooth particle hydrodynamical (SPH) simulations to investigate this in further detail.

4.3. Physical considerations

4.3.1. Stellar rotation

By adopting the solid-body approximation, we assume that torques will spin up or spin down each star as a whole. However, it is more likely that only the outer-most layers will be affected, thereby triggering differential rotation. Subsequently, angular momentum is re-distributed within the stellar interior via meridional circulation and shear instabilities (see Maeder 2009, for a review).

Song et al. (2013) studied the affect of angular momentum transport within a differentially rotating 15 M_⊙ primary, with a 10 M_⊙ companion as the primary evolved from the ZAMS until the onset of RLOF. They found that meridional circulation always counteracts the impact of tides; spinning up the surface layers when tides spin them down and vice versa, increasing the time for the star to rotate synchronously with the orbit. Consequently, some of their models commence RLOF before the primary is synchronised. In addition, differential rotation significantly enhances nitrogen abundances at the stellar surface.

We can speculate that significant differential rotation will occur while our donors become sub-synchronous, with meridional circulation initially opposing the spin-down triggered by rapid mass loss, and then the subsequent tidal forces which act to re-synchronise the donor. In turn, this may affect both the duration of the self-accretion phase and the value of α_self.

4.3.2. Tides

Tassoul (1987, 1988) proposed an alternative theory that invokes large-scale hydrodynamical flows as a means to dissipate kinetic energy, which are more efficient than the dynamical tide model described by Zahn (1977). Indeed, the synchronization and circularisation timescales between the two approaches vary by up to three orders of magnitude (Khaliullin & Khaliullina 2010).

We can speculate on the impact of Tassoul’s formalism on our calculations as follows. If τ_sync is 1000 times shorter, then for our case B system τ_sync<τ_R, i.e. tides keep the donor in synchronous rotation by the time RLOF starts. Additionally, for both the case A and case B models, an inspection of Figs. 2 and 7 shows that we would have τ_sync<τ_Ṁ. Hence, even during the fast phase, mass loss would not be sufficiently rapid to spin down the donor star to sub-synchronous rates and, in turn, self-accretion would not be triggered. In this case, the orbital evolution will correspond to the short-dashed green curves given in Figs. 4 and 8, where we neglected the term. Therefore, using Tassoul’s formalism will not significantly affect our results, namely that our osculating scheme still yields much shorter orbital periods.

However, the hydrodynamical model has since been criticised on theoretical grounds by Rieutord (1992) and Rieutord & Zahn (1997; but see Tassoul & Tassoul 1997, for a counter-argument), while observations attempting to constrain the mechanism underpinning tidal interactions are inconclusive. On the one hand, Claret et al. (1995) and Claret & Cunha (1997) found that both the Zahn and Tassoul formalisms can adequately account for the observed eccentricities and rotational velocities of early-type eclipsing binaries. On the other hand, using an updated sample of such binaries, Khaliullin & Khaliullina (2007, 2010) found that Tassoul’s theory is in contradiction with observations, which are better reproduced by Zahn’s formalism. However, Meibom & Mathieu (2005) derived, for stellar populations with a range of ages, the tidal circularization period (i.e. the orbital period at which the binary orbit circularises at the age of the population). Their results indicated that, at a given age, the observed circularisation period is larger than the predicted value from the dynamical tide model, suggesting that it is too inefficient. Clearly, more observational and theoretical work is required in this area.

4.3.3. Self-accretion

The phenomenon of self-accretion has also been found in the studies of Kruszewski (1964a) and Sepinsky et al. (2010), who use the same ballistic approach, which assumes that no collisions between particles occur. Clearly, this is not realistic, as there are indeed multiple collisions along the trajectory.

Nonetheless, our results are in qualitative agreement with the SPH calculations of Belvedere et al. (1993), who find that below some critical rotational velocity, the entire stream is deflected towards the primary. For larger spin rates, material falls onto both components. However, their study focused on the impact of asynchronous rotation on disc formation around the secondary, and they did not quantify self-accretion. We hope our work will motivate future SPH studies in this area.

5. Summary and conclusions

We use our stellar binary evolution code BINSTAR to calculate the evolution of Algol systems using the theory of osculating orbital elements. By calculating the ballistic trajectories of ejected particles from the mass losing star (the donor), we determine the change of linear momentum of each star, and the gravitational perturbation applied to the stars by the mass transfer stream. As a consequence of the latter, the osculating formalism predicts signicantly shorter post-mass transfer orbital periods, typically by a factor of 4, than the widely applied classical scheme.

Also contrary to widely held belief, the donor star does not remain in synchronous rotation with the orbital motion throughout mass exchange. The initially rapid mass ejection spins down the donor on a shorter timescale than the tidal synchronization timescale, enforcing sub-synchronous rotation and causing about 15 to 20 per cent of the ejected material to fall back onto the donor during these episodes of self-accretion. Self-accretion, combined with the sink of orbital angular momentum that mass transfer provides, may lead to the formation of more contact binary systems.

While we have mainly focused on conservative Algol evolution, the osculating prescription clearly applies to all varieties of interacting binaries. In the future, we will apply our osculating formalism to investigate the evolution of eccentric systems.

Online material

Appendix A: The perturbing forces and

Consider a primary star of mass M₁. At time t = t₀, it has an orbital velocity v₁ = v₀, and an orbital angular velocity ω. At a later time t = t′, a particle is ejected from the primary at the inner-Lagrangian point, located at a distance r_ℒ1 from the primary’s mass centre (see Fig. A.1), and so the primary’s mass becomes ${\mathit{M}}_{\mathrm{1}}\mathrm{+}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}$, where $\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}\mathit{<}\mathrm{0}$. The particle’s absolute velocity is ${{W}}_{\mathrm{ej}}^{\mathrm{\prime }}$. As a result, the centre of mass is shifted by $\mathit{\delta }{{r}}_{\mathrm{1}}\mathrm{=}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}{{\mathit{M}}_{\mathrm{1}}}{{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}$(A.1)with respect to its unperturbed location, and its new orbital velocity is ${{v}}_{\mathrm{1}}^{\mathrm{\prime }}$. During a time interval δt, the change in the primary’s orbital velocity is (see Hadjidemetriou 1969b; Sepinsky et al. 2007b, for further details) $\begin{array}{ccc}\frac{{{v}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{-}{{v}}_{\mathrm{1}}}{\mathit{\delta t}}& \mathrm{=}& \frac{{{Q}}_{\mathrm{1}}^{\mathrm{\prime }}\mathrm{-}{{Q}}_{\mathrm{1}}}{\mathit{\delta t}}\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{+}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}{\mathit{\delta t}}\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{{V}}_{\mathrm{ej}}\mathrm{+}{{\omega }}^{\mathrm{\prime }}\mathrm{\times }{{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{)}\\ & & \end{array}$(A.2)where Q₁ is the primary’s momentum, primed quantities indicate values at time t′, and V_ej = W_ej − v₁ is the relative velocity of the ejected material with respect to the primary’s mass centre.

At self-accretion, the primary accretes the particle of mass $\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}\mathrm{=}\mathrm{-}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}\mathit{>}\mathrm{0}$. Just before self-accretion occurs at time t = t′′, the momentum of the primary and ejected particle, ${{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime }}$, is ${{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{=}\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{-}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}\mathrm{)}{{v}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{+}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}{{W}}_{\mathrm{acc}}$(A.3)where W_acc is the absolute velocity of the self-accreted particle. The orbital velocity, ${{v}}_{\mathrm{1}}^{\mathrm{\prime \prime }}$, is the sum of the non-perturbed orbital velocity at time t′′, ${{v}}_{\mathrm{0}}^{\mathrm{\prime \prime }}$, and the perturbation to the velocity because of the original ejection episode, so ${{v}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{=}{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime }}\mathrm{+}\mathrm{(}{{\omega }}^{\mathrm{\prime }}\mathrm{+}\mathit{\delta }{\omega }\mathrm{)}\mathrm{\times }\mathit{\delta }{{r}}_{\mathrm{1}}\mathrm{=}{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime }}\mathrm{+}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{{\omega }}^{\mathrm{\prime }}\mathrm{\times }{{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{)}\mathit{,}$(A.4)

Fig. A.1 Illustration of the self-accretion process. At time t = t0, the donor star moves along its orbital path (solid black curve), with an orbital velocity v0 (cyan arrow). At t′, a particle of mass $\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}$ is ejected from the inner-Lagrangian point located at rℒ1 with respect to the donor’s mass centre. The ejection shifts the centre of mass by δr1, and the donor follows a new orbit (long dashed curve) with velocity v′. At t′′ just before the particle is re-accreted (at ${{r}}_{\mathrm{acc}}^{\mathrm{\ast }}$), the orbital velocity is v′′. Subsequently, self-accretion shifts the donor’s mass centre by δrself, and it follows the orbit indicated by the dot-dashed curve, with a velocity v′′′. The dashed circles represent the locations of the donor if no mass ejection had taken place, while the dotted circle indicates the donor’s location had no self-accretion occurred.

where we have used Eq. (A.1) and ignored terms larger than first-order. Inserting Eq. (A.4) into Eq. (A.3) gives ${{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime }}\mathrm{=}{\mathit{M}}_{\mathrm{1}}{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime }}\mathrm{+}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}\mathrm{(}{{\omega }}^{\mathrm{\prime }}\mathrm{\times }{{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{)}\mathrm{-}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime }}\mathrm{+}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}{{W}}_{\mathrm{acc}}\mathit{.}$(A.5)The particle is self-accreted at time t^′′′, at a location ${{r}}_{\mathrm{acc}}^{\mathrm{\ast }}$ with respect to the primary’s mass centre. This shifts the mass centre by an amount $\mathit{\delta }{{r}}_{\mathrm{self}}\mathrm{=}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}}{{\mathit{M}}_{\mathrm{1}}}{{r}}_{\mathrm{acc}}^{\mathrm{\ast }}$(A.6)and the primary’s momentum is now ${{Q}}^{\mathrm{\prime \prime \prime }}\mathrm{=}{\mathit{M}}_{\mathrm{1}}{{v}}_{\mathrm{1}}^{\mathrm{\prime \prime \prime }}\mathit{.}$(A.7)The new orbital velocity, ${{v}}_{\mathrm{1}}^{\mathrm{\prime \prime \prime }}$, is the sum of the unperturbed orbital velocity, ${\mathit{v}}_{\mathrm{0}}^{\mathrm{\prime \prime \prime }}$, and the perturbations to the velocity arising from the ejection and self-accretion processes, i.e. $\begin{array}{ccc}{{v}}_{\mathrm{1}}^{\mathrm{\prime \prime \prime }}& \mathrm{=}& {{v}}_{\mathrm{0}}^{\mathrm{\prime \prime \prime }}\mathrm{+}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{{\omega }}^{\mathrm{\prime }}\mathrm{\times }{{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{)}\mathrm{+}\mathrm{(}{{\omega }}^{\mathrm{\prime \prime }}\mathrm{+}\mathit{\delta }{\omega }\mathrm{)}\mathrm{\times }\mathit{\delta }{{r}}_{\mathrm{self}}\\ & \mathrm{=}& \end{array}$(A.8)where we have used Eq. (A.6), and once again ignored terms higher than first-order. Inserting Eq. (A.8) into Eq. (A.7) yields ${{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime \prime }}\mathrm{=}{\mathit{M}}_{\mathrm{1}}{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime \prime }}\mathrm{+}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}\mathrm{(}{{\omega }}^{\mathrm{\prime }}\mathrm{\times }{{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{)}\mathrm{+}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}\mathrm{(}{{\omega }}^{\mathrm{\prime \prime }}\mathrm{\times }{{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\mathrm{)}\mathit{.}$(A.9)Taking the difference between Eqs. (A.9) and (A.5), dividing the result by δt and using the fact that $\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}\mathrm{=}\mathrm{-}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}$ yields $\begin{array}{ccc}\frac{{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime \prime }}\mathrm{-}{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime }}}{\mathit{\delta t}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\frac{{{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime \prime }}\mathrm{-}{{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime }}}{\mathit{\delta t}}\mathrm{-}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}}{\mathit{\delta t}}\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{{\omega }}^{\mathrm{\prime \prime }}\mathrm{\times }{{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\mathrm{)}\\ & & \end{array}$(A.10)Following Sepinsky et al. (2007b), the absolute acceleration of the primary’s mass centre is the sum of Eq. (A.10), the relative acceleration of the primary’s mass centre, $\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{acc}\end{array}}{\mathrm{(}\mathit{\delta t}{\mathrm{)}}^{\mathrm{2}}}\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}$, and the Coriolis acceleration, $\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}}{\mathit{\delta t}}\mathrm{2}\mathrm{(}{{\omega }}^{\mathrm{\prime \prime }}\mathrm{\times }{{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\mathrm{)}$ to give $\begin{array}{ccc}\frac{{{v}}_{\mathrm{1}}^{\mathrm{\prime \prime \prime }}\mathrm{-}{{v}}_{\mathrm{1}}^{\mathrm{\prime \prime }}}{\mathit{\delta t}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\frac{{{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime \prime }}\mathrm{-}{{Q}}_{\mathrm{1}}^{\mathrm{\prime \prime }}}{\mathit{\delta t}}\mathrm{+}\frac{\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}}{\mathit{\delta t}}\frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{\left[}\mathrm{\right(}{{\omega }}^{\mathrm{\prime \prime }}\mathrm{\times }{{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\mathrm{)}\mathrm{+}{{V}}_{\mathrm{acc}}\mathrm{]}\\ & & \end{array}$(A.11)where ${{V}}_{\mathrm{acc}}\mathrm{=}{{W}}_{\mathrm{acc}}\mathrm{-}{{v}}_{\mathrm{0}}^{\mathrm{\prime \prime }}$ is the relative velocity of the self-accreted particle with respect to the primary’s mass centre.

Summing Eqs. (A.11) and (A.2) gives the acceleration of the primary’s mass centre from both the ejection and self-accretion process. In the limit δt → 0, and remembering that $\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{acc}}^{\mathrm{self}}\mathrm{=}\mathrm{-}\mathit{\delta }{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}$, then the acceleration of the primary’s mass centre is $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}{{R}}_{\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}& \mathrm{=}& \frac{{{F}}_{\mathrm{1}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{+}\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{[}\mathrm{(}{{V}}_{\mathrm{ej}}\mathrm{-}{{V}}_{\mathrm{acc}}\mathrm{)}\mathrm{+}\mathrm{(}{{\omega }}^{\mathrm{\prime }}\mathrm{\times }{{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{-}{{\omega }}^{\mathrm{\prime \prime }}\mathrm{\times }{{r}}_{\mathrm{acc}}^{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{]}}\\ & & \end{array}$(A.12)where R₁ is the position vector of the primary with respect to an inertial reference frame, and F₁ is the sum of all external forces acting on the primary, which writes as ${{F}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}\frac{{{R}}_{\mathrm{1}}}{{\mathit{R}}_{\mathrm{1}}}\mathrm{+}{{f}}_{\mathrm{1}}\mathit{,}$(A.13)where r = R₂ − R₁, R₂ is the position vector of the secondary, and f₁ is the force acting on the primary via the matter stream. Since the secondary is not accreting, its equation of motion is $\frac{{\mathrm{d}}^{\mathrm{2}}{{R}}_{\mathrm{2}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{1}}}{{\mathit{r}}^{\mathrm{2}}}\frac{{{R}}_{\mathrm{2}}}{{\mathit{R}}_{\mathrm{2}}}\mathrm{+}\frac{{{f}}_{\mathrm{2}}}{{\mathit{M}}_{\mathrm{2}}}\mathit{,}$(A.14)where f₂ is the force acting on the secondary by the accretion stream. Subtracting Eq. (A.12) from Eq. (A.14) gives the equation of motion of the secondary with respect to the primary, which is $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}{r}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}& \mathrm{=}& \mathrm{-}\frac{\mathit{G}\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{)}}{{\mathit{r}}^{\mathrm{3}}}{r}\mathrm{+}\frac{{{f}}_{\mathrm{2}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{{f}}_{\mathrm{1}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{-}\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{{\mathit{M}}_{\mathrm{1}}}[\mathrm{(}{{V}}_{\mathrm{ej}}\mathrm{-}{{V}}_{\mathrm{acc}}\mathrm{)}\\ & & \end{array}$(A.15)which takes the form $\frac{{\mathrm{d}}^{\mathrm{2}}{r}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\mathrm{-}\frac{\mathit{G}\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{)}}{{\mathit{r}}^{\mathrm{3}}}{r}\mathrm{+}{\mathrm{𝒮}}_{\mathrm{self}}\mathit{ê}\mathit{r}\mathrm{+}{\mathrm{𝒯}}_{\mathrm{self}}\mathit{ê}\mathit{t}\mathit{.}$(A.16)Here, ê_r is a unit vector pointing along r, and ê_t is a unit vector perpendicular to ê_r in the direction of the orbital motion. Taking the dot product of Eq. (A.16) with ê_r and ê_t respectively, yields $\begin{array}{ccc}{\mathrm{𝒮}}_{\mathrm{self}}& \mathrm{=}& \frac{{\mathit{f}}_{\mathrm{2}\mathit{,r}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}\mathit{,r}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{-}\frac{\stackrel{}{{{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}_{\mathrm{¨}}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{-}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\cos {\mathit{\psi }}^{\mathrm{\ast }}\mathrm{)}\\ & & \mathrm{-}\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{\mathit{V}}_{\mathrm{ej}\mathit{,r}}\mathrm{-}{\mathit{V}}_{\mathrm{acc}\mathit{,r}}\mathrm{+}{\mathit{\omega }}^{\mathrm{\prime \prime }}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\sin {\mathit{\psi }}^{\mathrm{\ast }}\mathrm{)}\mathit{,}\end{array}$(A.17)and $\begin{array}{ccc}{\mathrm{𝒯}}_{\mathrm{self}}& \mathrm{=}& \frac{{\mathit{f}}_{\mathrm{2}\mathit{,t}}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{{\mathit{f}}_{\mathrm{1}\mathit{,t}}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{-}\frac{\stackrel{}{{{\mathit{M}}_{\mathrm{1}\mathit{,}\mathrm{ej}}^{\mathrm{self}}}_{\mathrm{¨}}}}{{\mathit{M}}_{\mathrm{1}}}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\sin {\mathit{\psi }}^{\mathrm{\ast }}\\ & & \mathrm{-}\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{self}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}{\mathit{V}}_{\mathrm{ej}\mathit{,t}}\mathrm{-}{\mathit{V}}_{\mathrm{acc}\mathit{,t}}\mathrm{+}{\mathit{\omega }}^{\mathrm{\prime }}{\mathit{r}}_{{\mathrm{ℒ}}_{\mathrm{1}}}\mathrm{-}{\mathit{\omega }}^{\mathrm{\prime \prime }}{\mathit{r}}_{\mathrm{acc}}^{\mathrm{\ast }}\cos {\mathit{\psi }}^{\mathrm{\ast }}\mathrm{)}\mathit{,}\end{array}$(A.18)which are the same as Eqs. (11) and (12), noting that ω′ = ω′′ = ω for circular orbits. The quantity ψ^∗ is the angle between ê_r and the impact site, and the subscripts “r” and “t” indicate components along ê_r and ê_t respectively.

Appendix B: Torque arising from mass transfer,

Consider a primary star of mass M₁, and a secondary of mass M₂, separated by a distance r. They respectively orbit the common centre of mass with an orbital velocity v₁ and v₂. The velocity of the secondary with respect to the primary is v = v₂ − v₁, and the orbital angular momentum is given by ${\mathit{J}}_{\mathrm{orb}}\mathrm{=}\mathit{m}{\mathrm{[}\mathit{GMa}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{]}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathit{mr}{\mathit{v}}_{\mathit{t}}\mathit{,}$(B.1)where e is the eccentricity, m = M₁M₂/M is the reduced mass, M = M₁ + M₂, and v_t is the orbital velocity along ê_t, given by ${\mathit{v}}_{\mathit{t}}\mathrm{=}{\left[\frac{\mathit{GM}}{\mathit{a}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}\right]}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{\nu }\mathrm{)}$(B.2)and ν is the true anomaly. Taking the time derivative of the last equality in Eq. (B.1) and noting that ṙ = 0 for an osculating orbit (see, e.g. Bonačić Marinović et al. 2008), yields $\frac{\mathit{J̇}\mathrm{orb}}{{\mathit{J}}_{\mathrm{orb}}}\mathrm{=}\frac{\mathit{Ṁ}\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{+}\frac{\mathit{Ṁ}\mathrm{2}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{\mathit{Ṁ}}{\mathit{M}}\mathrm{+}\frac{\mathit{v̇}\mathit{t}}{{\mathit{v}}_{\mathit{t}}}\mathrm{·}$(B.3)Similarly to Eq. (A.16), the equation of motion of a binary acted on by perturbing forces and reads $\frac{{\mathrm{d}}^{\mathrm{2}}{r}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}\frac{\mathrm{d}{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{G}\mathrm{(}{\mathit{M}}_{\mathrm{1}}\mathrm{+}{\mathit{M}}_{\mathrm{2}}\mathrm{)}}{\mathrm{|}\mathit{r}{\mathrm{|}}^{\mathrm{3}}}{r}\mathrm{+}\mathrm{𝒮}\mathit{ê}\mathit{r}\mathrm{+}\mathrm{𝒯}\mathit{ê}\mathit{t}\mathit{.}$(B.4)Taking the dot product of Eq. (B.4) with ê_t gives $\frac{\mathrm{d}{\mathit{v}}_{\mathit{t}}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{𝒯}\mathit{.}$(B.5)Inserting Eqs. (B.5) and (B.2) into Eq. (B.3), and using the first equality in Eq. (B.1), gives the torque applied onto the orbit from mass transfer, $\mathit{J̇}\mathrm{orb}\mathit{,}\mathrm{MT}\mathrm{=}{\mathit{J}}_{\mathrm{orb}}\left(\frac{\mathit{Ṁ}\mathrm{1}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{+}\frac{\mathit{Ṁ}\mathrm{2}}{{\mathit{M}}_{\mathrm{2}}}\mathrm{-}\frac{\mathit{Ṁ}}{\mathit{M}}\right)\mathrm{+}\mathit{m}\frac{\mathit{a}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}{\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{\nu }}\mathrm{𝒯}\mathit{.}$(B.6)The net change of the primary’s mass is the sum of mass transferred to the companion via RLOF, $\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}\mathit{<}\mathrm{0}$, and mass ejected by the wind, Ṁ_1,loss< 0, i.e. $\mathit{Ṁ}\mathrm{1}\mathrm{=}\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}\mathrm{+}\mathit{Ṁ}\mathrm{1}\mathit{,}\mathrm{loss}\mathit{.}$(B.7)Similarly, for the secondary $\mathit{Ṁ}\mathrm{2}\mathrm{=}\mathrm{-}\mathit{\beta }\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}\mathrm{+}\mathit{Ṁ}\mathrm{2}\mathit{,}\mathrm{loss}\mathit{,}$(B.8)where the first term on the right hand side gives the accretion rate and Ṁ_2,loss includes the mass ejected from the system during non-conservative mass transfer. Substituting Eqs. (B.7) and (B.8) into Eq. (B.6) gives $\begin{array}{ccc}\mathit{J̇}\mathrm{orb}\mathit{,}\mathrm{MT}& \mathrm{=}& \begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{{\mathit{J}}_{\mathrm{orb}}\left[\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{\mathit{M}}\left(\frac{\mathrm{1}}{\mathit{q}}\mathrm{-}\mathit{\beta q}\right)\right]\mathrm{+}\mathit{m}\frac{\mathit{a}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}{\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{\nu }}\mathrm{𝒯}}\\ \mathit{J̇}\mathrm{orb}\mathit{,}\mathrm{RLOF}\end{array}\\ & & \end{array}$(B.9)where q = M₁/M₂, is the torque acting on the orbit as a consequence of RLOF, while is the torque applied by the material leaving the system. The corresponding torque applied onto the transferred mass is just $\mathit{J̇}\mathrm{MT}\mathrm{=}\mathrm{-}\mathit{J̇}\mathrm{orb}\mathit{,}\mathrm{MT}\mathit{.}$(B.10)Using Eqs. (B.9) and (B.10) with e = 0 gives Eq. (15).

Next, we demonstrate the consistency of Eq. (B.9) by showing that in the classical formalism for conservative mass transfer . If all material is transferred to the secondary (α_self = 0, β = 1), if the stars are treated as point masses (r_ℒ1 = 0, r_acc = 0), and if we neglect the gravitational attraction by the accretion stream (f₁ = 0, f₂ = 0), Eq. (10) reduces to ${\mathrm{𝒯}}_{\mathrm{comp}}\mathrm{=}\mathrm{-}\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{{\mathit{M}}_{\mathrm{1}}}\mathrm{(}\mathit{q}{\mathit{V}}_{\mathrm{2}\mathit{,t}}\mathrm{+}{\mathit{V}}_{\mathrm{1}\mathit{,t}}\mathrm{)}\mathit{,}$(B.11)where $\mathit{Ṁ}{}_{\mathrm{2}\mathit{,}\mathrm{acc}}\mathrm{=}\mathrm{-}\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}$ for conservative mass transfer. Sepinsky et al. (2007b) and Luk’yanov (2008) demonstrated that, if the

orbital angular momentum is conserved, in a circular orbit V_1,t and V_2,t are related by $\mathit{q}{\mathit{V}}_{\mathrm{2}\mathit{,t}}\mathrm{+}{\mathit{V}}_{\mathrm{1}\mathit{,t}}\mathrm{=}{\left(\frac{\mathit{GM}}{\mathit{a}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{q}\mathrm{)}\mathit{.}$(B.12)Substituting Eq. (B.12) into Eq. (B.11), and that result into Eq. (B.6) gives for a circular orbit $\begin{array}{ccc}\mathit{J̇}\mathrm{orb}\mathit{,}\mathrm{MT}& \mathrm{=}& \mathit{J̇}\mathrm{orb}\mathit{,}\mathrm{RLOF}\mathrm{=}{\mathit{J}}_{\mathrm{orb}}\frac{\mathit{Ṁ}\begin{array}{c}\mathrm{comp}\\ \mathrm{1}\mathit{,}\mathrm{ej}\end{array}}{{\mathit{M}}_{\mathrm{1}}}\left(\frac{\mathit{q}}{\mathrm{1}\mathrm{+}\mathit{q}}\right)\left(\frac{\mathrm{1}\mathrm{-}{\mathit{q}}^{\mathrm{2}}}{\mathit{q}}\right)\\ & & \end{array}$(B.13)where we have used M₁/M = q/ (1 + q). Using Eq. (B.1), Eq. (B.13) reduces to zero, as required.

Acknowledgments

P.J.D. acknowledges financial support from the FNRS Research Fellowship – Chargé de Recherche. L.S. is a FNRS Research Associate. We thank the anonymous referee whose constructive comments helped to improve the quality of the manuscript.

References