An analytic model for mass transfer in binaries with arbitrary eccentricity, with applications to triple-star systems

TL;DR

This paper introduces a new analytic model for mass transfer in binary systems with arbitrary eccentricity, enabling more accurate long-term evolution studies of such systems, including triples with Lidov-Kozai oscillations.

Contribution

The authors develop a smooth, phase-dependent mass transfer model valid for any eccentricity, overcoming previous high-eccentricity limitations and facilitating population synthesis applications.

Findings

Model valid for all eccentricities, including circular orbits.

Implementation in a publicly available code for population synthesis.

Application to systems with LK oscillations and eccentric binaries.

Abstract

Most studies of mass transfer in binary systems assume circular orbits at the onset of Roche lobe overflow. However, there are theoretical and observational indications that mass transfer could occur in eccentric orbits. In particular, eccentricity could be produced via sudden mass loss and velocity kicks during supernova explosions, or Lidov-Kozai (LK) oscillations in hierarchical triple systems, or, more generally, secular evolution in multiple-star systems. However, current analytic models of eccentric mass transfer are faced with the problem that they are only well defined in the limit of very high eccentricities, and break down for less eccentric and circular orbits. This provides a major obstacle to implementing such models in binary and higher-order population synthesis codes, which are useful tools for studying the long-term evolution of a large number of systems. Here, we present a new analytic model to describe the secular orbital evolution of binaries undergoing conservative mass transfer. The main improvement of our model is that the mass transfer rate is a smoothly varying function of orbital phase, rather than a delta function centered at periapsis. Consequently, our model is in principle valid for any eccentricity, thereby overcoming the main limitation of previous works. We implement our model in an easy-to-use and publicly available code that can be used as a basis for implementations of our model into population synthesis codes. We investigate the implications of our model in a number of applications with circular and eccentric binaries, and triples undergoing LK oscillations.

Figures

Equations

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}=-\frac{GM}{r^3}\mathbf{r}$$

$$+\frac{{\mathbf{f}}_{\mathrm{a}}}{M_{\mathrm{a}}}-\frac{{\mathbf{f}}_{\mathrm{d}}}{M_{\mathrm{d}}}$$

$$+\frac{{\dot{M}}_{\mathrm{a}}}{M_{\mathrm{a}}}\left({\mathbf{w}}_{\mathrm{a}}+{\mathbf{\omega }}_{\mathrm{orb}}\times {\mathbf{r}}_{{\mathrm{A}}_{\mathrm{a}}}\right)-\frac{{\dot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}\left({\mathbf{w}}_{\mathrm{d}}+{\mathbf{\omega }}_{\mathrm{orb}}\times {\mathbf{r}}_{{\mathrm{A}}_{\mathrm{d}}}\right)$$

$$+\frac{{\ddot{M}}_{\mathrm{a}}}{M_{\mathrm{a}}}{\mathbf{r}}_{{\mathrm{A}}_{\mathrm{a}}}-\frac{{\ddot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}{\mathbf{r}}_{{\mathrm{A}}_{\mathrm{d}}}.$$

$$\frac{\dot{a}}{a}=-2\frac{\dot{M_{\mathrm{d}}}}{M_{\mathrm{d}}}\left(1-\frac{M_{\mathrm{d}}}{M_{\mathrm{a}}}\right).$$

$$\frac{{\dot{L}}_{\mathrm{orb}}}{L_{\mathrm{orb}}}=\frac{{\dot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}+\frac{{\dot{M}}_{\mathrm{a}}}{M_{\mathrm{a}}}-\frac{1}{2}\frac{\dot{M}}{M}+\frac{1}{2}\frac{\dot{a}}{a}-\frac{e\dot{e}}{1-e^2}$$

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}$$

$$⟨\frac{\dot{a}}{a}⟩$$

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}$$

$$-\frac{{\ddot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}\left({\mathbf{r}}_{{\mathrm{A}}_{\mathrm{d}}}+q{\mathbf{r}}_{{\mathrm{A}}_{\mathrm{a}}}\right)$$

$$\equiv -\frac{GM}{r^3}\mathbf{r}-\frac{{\dot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}\left(1-q\right)\dot{\mathbf{r}}-\frac{{\dot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}{\mathbf{\omega }}_{\mathrm{orb}}\times {\mathbf{r}}_{\mathrm{A}}-\frac{{\ddot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}{\mathbf{r}}_{\mathrm{A}},$$

$${\mathbf{r}}_{\mathrm{A}}\equiv {\mathbf{r}}_{{\mathrm{A}}_{\mathrm{d}}}+q{\mathbf{r}}_{{\mathrm{A}}_{\mathrm{a}}}.$$

$$\frac{q}{X_{\mathrm{L}}^2}-\frac{1}{(1-X_{\mathrm{L}}{)}^2}-X_{\mathrm{L}}(1+q)\mathcal{A}(\hat{\Omega },e,\theta )+1=0.$$

$$\mathcal{A}(\hat{\Omega },e,\theta )=\frac{{\hat{\Omega }}^2(1+e{)}^4}{(1+e\cos \theta {)}^3}={\hat{\Omega }}^2\frac{1+e}{(1-e{)}^3}{\left(\frac{r}{a}\right)}^3.$$

$$X_{\mathrm{L}}={\mathcal{A}}^{-1/3}={\hat{\Omega }}^{-2/3}\frac{1-e}{(1+e{)}^{1/3}}\frac{a}{r}\equiv X_{\mathrm{L},0}(e,\hat{\Omega })\frac{a}{r},$$

$$X_{\mathrm{L},0}(e,\hat{\Omega })\equiv {\hat{\Omega }}^{-2/3}\frac{1-e}{(1+e{)}^{1/3}}.$$

$$\displaystyle\mathbf{r}_{\mathrm{A}}=\left\{\begin{array}[]{cc}\displaystyle\left[r(t)X_{\mathrm{L,0}}(q)-qr_{\mathrm{A}_{\mathrm{a}}}\right]\hat{\mathbf{r}}(t),&\text{(Case 1)}\\ \\ \displaystyle\left[aX_{\mathrm{L,0}}(e,\hat{\Omega})-qr_{\mathrm{A}_{\mathrm{a}}}\right]\hat{\mathbf{r}}(t).&\text{(Case 2)}\end{array}\right.$$

$$R_{\mathrm{L}}(t)=\frac{R_{\mathrm{L}}^{\mathrm{c}}}{a}r(t),$$

$$\frac{R_{\mathrm{L}}^{\mathrm{c}}}{a}=\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln \left(1+q^{1/3}\right)}$$

$$x\equiv \frac{R_{\mathrm{L}}^{\mathrm{c}}}{R}=\frac{a}{R}\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln \left(1+q^{1/3}\right)},$$

$$r(\mathcal{E})=a(1-e\cos \mathcal{E}),$$

$$\cos {\mathcal{E}}_0=\frac{1}{e}\left(1-\frac{1}{x}\right).$$

$$\dot{M_{\mathrm{d}}}\propto \exp \left(\frac{R-R_{\mathrm{L}}}{H_{\mathrm{P}}}\right),$$

$$\dot{M_{\mathrm{d}}}\propto {\left(\frac{R-R_{\mathrm{L}}}{R}\right)}^{n_{\mathrm{p}}+3/2},$$

$${\dot{M}}_{\mathrm{d}}={\dot{M}}_{\mathrm{d},0}{\left(\frac{R-R_{\mathrm{L}}}{R}\right)}^3={\dot{M_{\mathrm{d}}}}_{,0}{\left[1-x(1-e\cos \mathcal{E})\right]}^3,$$

$${\ddot{M}}_{\mathrm{d}}=-3nxe{\dot{M}}_{\mathrm{d},0}{\left[1-x(1-e\cos \mathcal{E})\right]}^2\frac{\sin \mathcal{E}}{1-e\cos \mathcal{E}},$$

$${\dot{M}}_{\mathrm{d}}$$

$$={\dot{M_{\mathrm{d}}}}_{,0}{\left[1-x\left\{1-e\cos (\mathcal{E}-{\mathcal{E}}_{\tau })\right\}\right]}^3,$$

$$n\tau ={\mathcal{E}}_{\tau }-e\sin {\mathcal{E}}_{\tau }.$$

$$\frac{{\mathrm{d}}^2\mathbf{r}}{\mathrm{d}{t}^2}$$

$$\mathbf{f}=-\frac{{\dot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}\left(1-q\right)\dot{\mathbf{r}}-\frac{{\dot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}{\mathbf{\omega }}_{\mathrm{orb}}\times {\mathbf{r}}_{\mathrm{A}}-\frac{{\ddot{M}}_{\mathrm{d}}}{M_{\mathrm{d}}}{\mathbf{r}}_{\mathrm{A}},$$

$$\dot{E}=-E(\dot{a}/a)=\dot{\mathbf{r}}\cdot \mathbf{f}.$$

$$\mathbf{e}=\frac{1}{GM}\dot{\mathbf{r}}\times \mathbf{h}-\hat{\mathbf{r}}=\frac{1}{GM}\left[\mathbf{r}(\dot{\mathbf{r}}\cdot \dot{\mathbf{r}})-\dot{\mathbf{r}}(\mathbf{r}\cdot \dot{\mathbf{r}})\right]-\hat{\mathbf{r}},$$

Full text

1. Introduction

A common evolutionary process in binary and multiple-star systems is the transfer of mass between stars. Mass transfer is thought to be responsible for a wide range of phenomena such as X-ray emission in low- and high-mass X-ray binaries (see, e.g., Verbunt 1993; Remillard & McClintock 2006 for reviews), spin-up of neutron stars (e.g., Lorimer 2008), and cataclysmic variables (e.g., Knigge 2011).

Since tides are generally thought to be efficient in close binary systems (Zahn, 1977; Shu & Lubow, 1981), most theoretical studies of mass transfer assume that the orbit has circularized by the time of onset of mass transfer (Hurley et al., 2002; Pols et al., 2003). However, various processes are known to be able to excite significant eccentricity, even when tides are taken into account during the evolution of the binary. These processes include sudden mass loss and an imparted velocity kick during supernova explosions (e.g., Hills 1983; Brandt & Podsiadlowski 1995; Kalogera 1996), enhanced mass loss at periapsis (Soker, 2000; Bonačić Marinović et al., 2008), interactions with a massive circumbinary disk (Dermine et al., 2013; Antoniadis, 2014; Rafikov, 2016; Muñoz et al., 2018), or Lidov-Kozai (LK) oscillations (Lidov 1962; Kozai 1962; see Naoz 2016 for a review) if the binary is orbited by a third star (or, more generally, secular evolution in higher-order systems including, but not limited to, quadruple systems, e.g., Pejcha et al. 2013; Hamers et al. 2015; Hamers & Portegies Zwart 2016; Hamers & Lai 2017; Grishin et al. 2018). In particular, population synthesis studies of triple stars find that of the order of 10% of systems undergo mass transfer in eccentric orbits at some point in their evolution (Toonen et al., 2016, 2018). Similarly, mass transfer in eccentric orbits can be triggered by secular evolution in quadruple-star systems (Hamers, 2018, 2019). From an observational side, semidetached binaries are known to have nonzero eccentricities (Petrova & Orlov, 1999; Bonačić Marinović et al., 2008; Vos et al., 2013; Boffin et al., 2014), as well as high-mass X-ray binaries (Raguzova & Popov, 2005), and post-AGB binaries (Van Winckel et al., 1995).

Despite the relevance of nonzero eccentricity in mass transfer processes, binary population synthesis codes such as StarTrack (Belczynski et al., 2008), BSE (Hurley et al., 2002) and the updated binary_c (Izzard et al., 2004, 2006, 2009; Claeys et al., 2014), and SeBa (Portegies Zwart & Verbunt, 1996; Toonen et al., 2012), enforce circular orbits at the onset of mass transfer. Nevertheless, the problem of mass loss/mass transfer has been studied for over half a century (e.g., Huang 1956; Hadjidemetriou 1963; Kruszewski 1964; Piotrowski 1964; Matese & Whitmire 1983, 1984), and has received more recent attention in numerical studies (e.g., Regös et al. 2005; Church et al. 2009; Sepinsky et al. 2010; Lajoie & Sills 2011a, b; van der Helm et al. 2016; Bobrick et al. 2017), as well as in (semi)analytical work (Sepinsky et al., 2007b, 2009, 2010; Veras et al., 2011; Veras & Tout, 2012; Veras et al., 2013, 2014; Dosopoulou & Kalogera, 2016a, b).

Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b) in particular derived equations for the secular (i.e., orbit-averaged) changes of the orbital elements due to mass transfer in eccentric binaries. For the case of Roche Lobe overflow (RLOF), they assumed that the mass transfer rate is a delta function centered at periapsis, i.e., the donor star transfers its mass in a burst at its closest approach to its companion. This assumption is physically reasonable in the limit of very high eccentricity. However, it is clearly no longer reasonable for less eccentric or even circular orbits. In the latter case, the mass transfer rate is expected to be constant during the orbit. Moreover, periapsis is no longer defined for a circular orbit.

In the delta function model of eccentric RLOF of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b), the equations of motion state that the eccentricity time derivative is negative at zero eccentricity, provided that the mass ratio $q$ of the donor to the accretor mass is $q>1$ (and assuming point masses). This has the practical implication that, when solving the equations of motion (a set of first-order ordinary differential equations [ODEs]) numerically, the eccentricity becomes negative as the system evolves toward circularization (see Fig. 3 in Section 3.1 below for an example). Clearly, this is an undesirable property, especially when the equations are to be implemented in population synthesis codes. The case $q<1$ is problematic as well, since for $q<1$ the eccentricity time derivative is positive at zero eccentricity, and the equations of motion predict a growing eccentricity, which ultimately leads to significant deviation from the expected evolution of the semimajor axis for strictly circular orbits (see equation 2 below and Fig. 4 for an example).

In this paper, we present an analytic model for mass transfer in eccentric orbits that shares some of the same basic assumptions as those of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b), but assumes a physically motivated and more realistic model for the mass transfer rate. Consequently, our model is in principle valid for any eccentricity, including zero, eliminating the issues described above. We give explicit expressions for the orbit-averaged equations of motion and apply them to mass-transferring isolated binaries and triple systems. We also make an easy-to-use Python code publicly available (see Section 2.6 for the url) to quickly solve the equations of motion numerically, and which can be used as a basis for implementations of our model into binary (and higher-order multiplicity) population synthesis codes.

The plan of the paper is as follows. The analytic model is presented in Section 2. We give applications of the model in isolated binaries in Section 3, and in triple systems undergoing LK oscillations in Section 4. We discuss limitations of our model in Section 5, and conclude in Section 6.

2. The analytic model

Consider a binary with a donor star with mass $M_{\mathrm{d}}$ and radius $R$, and an accreting star with mass $M_{\mathrm{a}}$. Let the total mass be denoted with $M\equiv M_{\mathrm{d}}+M_{\mathrm{a}}$, the mass ratio with $q\equiv M_{\mathrm{d}}/M_{\mathrm{a}}$, the reduced mass with $\mu\equiv M_{\mathrm{d}}M_{\mathrm{a}}/M$, and the relative separation vector between the centers of mass of the stars with $\mathbf{r}\equiv\mathbf{R}_{\mathrm{a}}-\mathbf{R}_{\mathrm{d}}$, where $\mathbf{R}_{\mathrm{d}}$ and $\mathbf{R}_{\mathrm{a}}$ are the absolute position vectors of the centers of mass of the donor and accretor star, respectively (here, ‘absolute’ means relative to an inertial reference frame). The orbital angular frequency vector is $\bm{\omega}_{\mathrm{orb}}$, which is directed along the orbital angular momentum vector, and which has a magnitude $\omega_{\mathrm{orb}}=\dot{\theta}=n\sqrt{1-e^{2}}\,(a/r)^{2}$, where $\theta$ is the true anomaly, the dot denotes the derivative with respect to time $t$, $n=\sqrt{GM/a^{3}}$ is the mean motion, and $a$ and $e$ are the orbital semimajor axis and eccentricity, respectively.

Assume that the donor loses mass with a mass-loss rate $\dot{M}_{\mathrm{d}}$, and an absolute velocity $\mathbf{W}_{\mathrm{d}}$ at the position $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ relative to the center of mass of the donor (i.e., at the Lagrangian point $L_{1}$ in the case of RLOF). The accretor accretes mass with a rate $\dot{M}_{\mathrm{a}}$, and absolute velocity $\mathbf{W}_{\mathrm{a}}$ at the position $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ relative to the center of mass of the accretor ($||\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}||$ could be the accretor’s radius in the case of direct impact accretion, or, e.g., the size of the accretion disk if it is present). The ejection/accretion velocities relative to the donor and accretor are $\mathbf{w}_{\mathrm{d}}\equiv\mathbf{W}_{\mathrm{d}}-\mathbf{V}_{\mathrm{d}}$ and $\mathbf{w}_{\mathrm{a}}\equiv\mathbf{W}_{\mathrm{a}}-\mathbf{V}_{\mathrm{a}}$, respectively, where $\mathbf{V}_{\mathrm{d}}\equiv\mathrm{d}\mathbf{R}_{\mathrm{d}}/\mathrm{d}t$ and $\mathbf{V}_{\mathrm{a}}\equiv\mathrm{d}\mathbf{R}_{\mathrm{a}}/\mathrm{d}t$ are the absolute velocities of the centers of mass of the donor and accretor, respectively.

2.1. Perturbing acceleration due to mass transfer

2.1.1 General equations of motion

We adopt the approach in which the effects of mass transfer are treated as perturbations to the instantaneous (osculating) Kepler orbit of the binary. Note that another approach is to consider changes in the total binary orbital energy and angular momentum (e.g., Huang 1956; Bonačić Marinović et al. 2008). As derived by Hadjidemetriou (1969b) and re-derived by Sepinsky et al. (2007b), the acceleration of the relative position vector $\mathbf{r}$ can be written as

[TABLE]

The first term after the equality in the first line (1a) is the Keplerian acceleration, and all other terms represent perturbations associated with mass transfer. The terms in the second line (1b) represent perturbations from the ejected mass on the orbit. These terms are generally hard to calculate analytically, and typically require numerical integration of the trajectories of the particles in the mass transfer stream (e.g., Hadjidemetriou 1969a; Sepinsky et al. 2010; Davis et al. 2014). The terms in the third line (1c) are associated with the change of linear momentum of the accretor and donor due to mass transfer. The terms in the fourth line (1d) are due to the acceleration of the centers of mass of the accretor and donor.

We note that equation (1) was also derived by Matese & Whitmire (1983, 1984), although, due to an error made in the equations for the absolute position vectors of the components in terms of the relative separation, the acceleration term (line 1d) was missing, as was pointed out in section 3.3 of Sepinsky et al. (2007b).

2.1.2 Simplified equations – overview

The general equations of motion, equation (1), are hard to model analytically without detailed (hydrodynamical) simulations (see, e.g., Edwards & Pringle 1987; Regös et al. 2005; Church et al. 2009; Sepinsky et al. 2010; Lajoie & Sills 2011a, b; van der Helm et al. 2016; Bobrick et al. 2017 for numerical/hydrodynamical studies). We make a number of simplifying assumptions in order to arrive at equations that are analytically tractable. These assumptions are summarized below.

We assume that the effects of the mass stream on the orbit are negligible, i.e., we set $\mathbf{f}_{\mathrm{a}}=\mathbf{f}_{\mathrm{d}}=\mathbf{0}$. 2. 2.

We assume conservative mass transfer, i.e., $\dot{M}_{\mathrm{d}}=-\dot{M}_{\mathrm{a}}$, such that $\dot{M}=0$. 3. 3.

We assume that the donor ejects mass at a relative velocity of $\mathbf{w}_{\mathrm{d}}=\dot{\mathbf{r}}$, and that the accretor accretes mass at a relative velocity of $\mathbf{w}_{\mathrm{a}}=-\dot{\mathbf{r}}$. 4. 4.

We assume that $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ and $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ corotate with the orbit, i.e., they are proportional to $\hat{\mathbf{r}}$. In addition, we take $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ to have a constant magnitude, whereas we make two limiting assumptions of the magnitude of $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$: proportional to $r(t)$, and constant as a function of orbital phase.

The third assumption in particular may seem peculiar. However, we show below in Section 2.1.3 that these assumptions, in the case of point masses with zero-size ejection/accretion radii ($\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{0}$), lead to the canonical relation for the change of the semimajor axis due to conservative mass transfer in circular orbits.

2.1.3 Simplified equations – point masses and circular orbits

Here, we show that the assumptions of Section 2.1.2, combined with $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{0}$, yield the canonical relation for mass transfer in circular orbits, i.e.,

[TABLE]

This relation is usually derived using conservation of the orbital angular momentum, $L_{\mathrm{orb}}=\mu\sqrt{GMa(1-e^{2})}$, i.e., by setting

[TABLE]

to zero, combined with the assumptions of conservative mass transfer ($\dot{M}=0$), and circular orbits ($e=0$). Equation (2) states that the orbit shrinks when the donor is more massive than the accretor; when the mass ratio has reversed, the orbit expands (note that $\dot{M}_{\mathrm{d}}<0$).

With the assumptions outlined above in Section 2.1.2 and setting $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{0}$, equation (1) reduces to

[TABLE]

Consider the secular orbital evolution implied by equation (4). If the orbits are circular, then it is reasonable to assume that the mass transfer rate $\dot{M}_{\mathrm{d}}$ is constant as well. In addition, if the orbital timescale is short compared to the mass transfer timescale, $P_{\mathrm{orb}}\ll M_{\mathrm{d}}/\dot{M}_{\mathrm{d}}$, then we can assume that $M_{\mathrm{d}}$ and $M_{\mathrm{a}}$ are approximately constant during the orbit (i.e., the adiabatic approximation). The resulting secular semimajor axis change (see Section 2.4.1 below for details) is

[TABLE]

We thus arrive at the canonical relation, equation (2).

We note that equation (4) is also commonly used in other studies of mass transfer (e.g., Eggleton 2006, eq. C98; Kashi & Soker 2018).

2.1.4 Simplified equations – extended bodies with nonzero ejection/accretion radii

With the assumptions described in Section 2.1.2 and not setting the ejection/accretion radii to zero, the equations of motion (equation 1) can be written as

[TABLE]

where we introduced, for convenience, the short-hand notation

[TABLE]

When averaging over the orbit below in Section 2.4, we follow the adiabatic approximation and set $M_{\mathrm{d}}$ and the mass ratio $q$ to be constant. We do not assume a constant $\mathbf{r}_{\mathrm{A}}$, as is described in more detail below in Section 2.2.

Equation (6) is consistent with the equations of motion assumed by Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b). In further steps below, we deviate from these works. In particular, we make different assumptions of the ejection/accretion radii (Section 2.2), and the instantaneous mass transfer rate (Section 2.3).

2.2. Assumptions of the ejection/accretion radii

The ejection/accretion positions $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ and $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ describe the locations of the ejected/accreted mass relative to the donor/accretor star. Seen from an inertial frame, the vectors $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ and $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ rotate. Here, we assume that $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ is aligned with the separation vector $\mathbf{r}$ between the two stars, and points towards the accretor, i.e., $\hat{\mathbf{r}}_{\mathrm{A}_{\mathrm{d}}}=\hat{\mathbf{r}}$. We assume that the mass is accreted by the accretor on the same axis, but in the opposite direction, i.e., $\hat{\mathbf{r}}_{\mathrm{A}_{\mathrm{a}}}=-\hat{\mathbf{r}}$.

In general, the magnitudes of both $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ and $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ could be functions of the orbital phase. In the case of RLOF, $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ is the location of the first Lagrangian point $L_{1}$, interpreted to be a function of orbital phase, and insofar as this point can be defined for eccentric orbits. In the limit that the dynamical timescale of the donor is much shorter than the timescale associated with the orbital angular velocity and donor rotation (also known as the first approximation, Limber 1963), Sepinsky et al. (2007a) showed that the stationary point $L_{1}$ between the two stars can be determined according to the equation

[TABLE]

Here, $X_{\mathrm{L}}\equiv X/r$ is the location of $L_{1}$ relative to the donor’s center of mass normalized to the orbital separation, $\hat{\Omega}\equiv\Omega/\omega_{\mathrm{orb},\,\mathrm{P}}$ is the donor’s spin frequency $\Omega$ normalized to $\omega_{\mathrm{orb},\,\mathrm{P}}=n(1+e)^{1/2}/(1-e)^{3/2}$, the orbital angular frequency at periapsis, and

[TABLE]

Unfortunately, no (simple) analytic solutions exist for $X_{\mathrm{L}}$ in equation (8) as a function of $q$ and $\mathcal{A}$. Instead, we make two limiting assumptions.

Negligible donor spin: $\hat{\Omega}\approx 0$. 2. 2.

Large mass ratio: $q\gg 1$.

In the first case, $\mathcal{A}\approx 0$, and we can neglect the associated term in equation (8). Consequently, $X_{\mathrm{L}}$ is a function only of $q$. An analytic (although not very illuminating) solution for $X_{\mathrm{L}}=X_{\mathrm{L,0}}(q)$ exists in this case, and is given explicitly in Appendix  A. Note that, although $X_{\mathrm{L}}$ is not a function of orbital phase, the location of $L_{1}$ itself, $X=X_{\mathrm{L,0}}(q)r(t)$, varies along the orbit.

In the second case, we can ignore any terms not involving $q$ in equation (8), giving the simple solution

[TABLE]

where we defined the function

[TABLE]

The magnitude of $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$, $r_{\mathrm{A}_{\mathrm{a}}}$, could be the radius of the accretor star, or, if present, the size of the accretion disk. Generally, $r_{\mathrm{A}_{\mathrm{a}}}$ could vary along the orbit. However, for simplicity, we assume that $r_{\mathrm{A}_{\mathrm{a}}}$ is constant.

In summary, for the ejection/accretion locations we assume

[TABLE]

Here, under the adiabatic approximation, $e$, $q$, $\hat{\Omega}$, and $r_{\mathrm{A}_{\mathrm{a}}}$ are taken to be constant along one orbit when orbit-averaging (see Section 2.4).

2.3. Model for the mass transfer rate

In order to proceed to find the orbit-averaged rates of change of the orbital elements implied by the simplified equations of motion (equation 6), we need to specify the mass transfer rate as a function of orbital phase, i.e., $\dot{M}_{\mathrm{d}}=\dot{M}_{\mathrm{d}}(t)$. We recall from the Introduction that Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b) assumed instantaneous transfer at periapsis, i.e., a delta function $\dot{M}_{\mathrm{d}}\propto\delta(\theta)$. Here, we assume a more general model that resembles a delta function at very high eccentricity, but is also well defined for lower eccentricities, in particular, for circular orbits in which mass transfer is expected to occur continuously during the orbit.

We assume that the instantaneous Roche lobe radius of the donor is given by

[TABLE]

where $r(t)$ is the instantaneous orbital separation, and

[TABLE]

is a fit to the Roche lobe radius in a circular orbit, provided by Eggleton (1983). A more accurate but more complicated expression for $R_{\mathrm{L}}(t)$ is given by Sepinsky et al. (2007a). For the moment, let the semimajor axis $a$ and the eccentricity $e$ be fixed. We consider three cases with regard to RLOF.

$R<R_{\mathrm{L}}(t)$: the donor does not fill its Roche lobe during any orbital phase (‘no RLOF’). With $\dot{M_{\mathrm{d}}}=0$, equation (6) trivially reduces to the unperturbed two-body problem, and the orbit remains unchanged ($\dot{a}=\dot{e}=0$). 2. 2.

$R\geq R_{\mathrm{L}}(t)$ for any orbital phase: the donor fills its Roche lobe during the entire orbit (‘full RLOF’). This is the case if the orbit is (close to) circular. 3. 3.

$R\geq R_{\mathrm{L}}(t)$ for a range of orbital phases: the donor fills its Roche lobe during part of the orbit (‘partial RLOF’). This occurs if the orbit is (sufficiently) eccentric.

Defining

[TABLE]

these conditions can be written as

$\displaystyle x>\frac{1}{1-e}$;  (no RLOF) 2. 2.

$\displaystyle x\leq\frac{1}{1+e}$;  (full RLOF) 3. 3.

$\displaystyle\frac{1}{1+e}<x\leq\frac{1}{1-e}$.  (partial RLOF)

The three cases are illustrated graphically in Fig. 1 with the boundaries plotted in the $(e,x)$ plane, and with the hatched regions corresponding to RLOF (partial or full). For $e=0$, the three cases reduce to two cases consistent with the standard picture of RLOF in circular orbits: no RLOF if $x>1$, and (full) RLOF if $x\leq 1$. In the circular case, no RLOF is possible if $x>1$ (i.e., $R<R_{\mathrm{L}}^{\mathrm{c}}$). In the eccentric case, however, partial RLOF is possible provided that $e$ is large enough, i.e., $e>1-1/x$.

The range of orbital phases for which RLOF occurs in case (3) is given by $-\mathcal{E}_{0}<\mathcal{E}<\mathcal{E}_{0}$, where $\mathcal{E}$ is the eccentric anomaly (of course, the mean and true anomalies are equally valid variables to describe this range). Using the canonical relation

[TABLE]

$\mathcal{E}_{0}$ is given by

[TABLE]

Note that, strictly speaking, case (2) can be considered as a special case of case (3), i.e., case (3) reduces to case (2) if $\mathcal{E}_{0}=\pi$.

One approach to model the phase-dependent mass transfer rate might be to assume a step function that is zero in case (1), and constant in cases (2) and (3) in the range $-\mathcal{E}_{0}<\mathcal{E}<\mathcal{E}_{0}$. However, a physically more realistic model should take into account the known property that the mass transfer rate is highly sensitive to the ‘radius excess’ $\Delta R\equiv R-R_{\mathrm{L}}$, i.e., the degree at which the donor overflows its Roche lobe (e.g., Lajoie & Sills 2011a). Therefore, the mass transfer rate is expected to be higher closer to periapsis, where the (instantaneous) Roche lobe radius is smaller.

For a donor with an atmospheric scale height $H_{\mathrm{P}}$, Ritter (1988) derived that the mass transfer rate should be

[TABLE]

where $H_{\mathrm{P}}$ is the pressure scale height of the donor’s atmosphere. Equation (21) states that the transfer rate increases exponentially with $\Delta R$. Another result can be obtained by assuming a polytropic equation of state for the donor coupled with Bernoulli’s equation, giving (Paczyński & Sienkiewicz 1972; Edwards & Pringle 1987 111See also section 7.1, pages 8-9, of lecture notes on binary star evolution by Onno Pols, which are downloadable at the url https://www.astro.ru.nl/$\sim$onnop/education/binaries_utrecht_notes/Binaries_ch6-8.pdf.)

[TABLE]

where $n_{\mathrm{p}}$ is the polytropic index.

Ideally, one should consider both cases: equation (21), and equation (22) with a general $n_{\mathrm{p}}$. However, with these functional relations for $\dot{M}_{\mathrm{d}}$, we were unable to analytically compute the resulting integrals required to derive the orbit-averaged equations of motion (see Section 2.4 below). Fortunately, however, we were able to do so in the specific case of equation (22) with $n_{\mathrm{p}}=3/2$, which is a reasonable approximation for a wide range of objects such as convective stars and low-mass white dwarfs (e.g., Chandrasekhar 1939), and gas giant planets (e.g., Weppner et al. 2015).

Therefore, we choose to adopt the power-law dependence of $\dot{M}_{\mathrm{d}}$ with $n_{\mathrm{p}}=3/2$, i.e., in our model, we set

[TABLE]

where $\dot{M}_{\mathrm{d},\,0}$ is a phase-independent mass transfer rate. We note that equation (23) states that $\dot{M}_{\mathrm{d}}$ is a strong function of $\Delta R$, similarly to equation (21). Also, we emphasize that equation (23) only gives a qualitative description of the mass transfer rate; other choices might be equally valid.

Below, we relate $\dot{M}_{\mathrm{d},\,0}$ to the orbit-averaged mass transfer rate $\langle\dot{M}_{\mathrm{d}}\rangle$; we assume that $\langle\dot{M}_{\mathrm{d}}\rangle$ is known. For illustration, we plot in Fig. 2 the mass transfer rate $\dot{M}_{\mathrm{d}}$ normalized to $\langle\dot{M}_{\mathrm{d}}\rangle$ as a function of the eccentric anomaly, $\mathcal{E}$, for different values of $x$ and $e$. For low eccentricity and $x$ close to zero, the mass transfer rate is nearly constant. As $x$ and $e$ increase, the mass transfer rate becomes increasingly peaked around $\mathcal{E}=0$ (with a decreasing width, and increasing height).

From the equations of motion, equation (1), it is clear that the double time derivative $\ddot{M}_{\mathrm{d}}$ should also be specified. From equation (23), we directly obtain $\ddot{M}_{\mathrm{d}}$ assuming that $\dot{M_{\mathrm{d}}}_{,0}$, $x$, and $e$ are constant during the orbit, i.e.,

[TABLE]

where we used the canonical relation $\dot{\mathcal{E}}=na/r$.

As a further sophistication, we introduce a time lag $\tau$ (with dimensions of time) to take into account the known phenomenon that the mass transfer rate is not symmetric around the apsidal line, but peaks just after periapsis (e.g., van der Helm et al. 2016). In this case, we evaluate the Roche lobe radius in the expression for the mass transfer rate at $t-\tau$, i.e.,

[TABLE]

where $\mathcal{E}_{\tau}$ is the eccentric anomaly corresponding to the time interval $\Delta t=\tau$ as measured from periapsis. The quantity $\mathcal{E}_{\tau}$ is a function of $\tau$, the mean motion $n$, and the eccentricity, and is found using the Kepler equation, i.e.,

[TABLE]

For small $\tau$ ($n\tau\ll 1$), equation (26) yields the approximate solution $\mathcal{E}_{\tau}\approx n\tau/(1-e)$. Below, we distinguish between the cases $\tau=0$, and $\tau\neq 0$. It turns out that the expressions associated with the ejection/accretion radii in the orbit-averaged equations become excessively complicated if $\tau\neq 0$ (i.e., with hundreds of terms appearing). Therefore, when presenting the equations below in Section 2.4, in the case $\tau\neq 0$ we include only the terms that appear if $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{0}$.

2.4. Orbit-averaged equations of motion

2.4.1 Preliminaries

For completeness, we first give the standard equations that describe the evolution of the orbital elements of the binary (see, e.g., Appendix C of Eggleton 2006). Write the relative acceleration as

[TABLE]

where $\mathbf{f}$ is the perturbing acceleration (force per unit mass). In our case, the perturbing acceleration is

[TABLE]

where $\dot{\mathbf{r}}$, $\bm{\omega}_{\mathrm{orb}}=n\sqrt{1-e^{2}}(a/r)^{2}\,\hat{\mathbf{h}}$, $\dot{M}_{\mathrm{d}}$, $\ddot{M}_{\mathrm{d}}$, and $\mathbf{r}_{\mathrm{A}}$ are all functions of the orbital phase (see equations 23, 24 and 15 for the latter three terms, respectively), whereas $M_{\mathrm{d}}$ and $q$ are assumed to be phase independent (adiabatic regime).

The perturbation $\mathbf{f}$ gives rise to a change in the specific orbital energy $E=\frac{1}{2}(\dot{\mathbf{r}}\cdot\dot{\mathbf{r}})-GM/r$ given by

[TABLE]

Furthermore, the change in the eccentricity vector,

[TABLE]

where $\mathbf{h}\equiv\mathbf{r}\times\dot{\mathbf{r}}$ is the specific angular momentum vector, is given by

[TABLE]

The specific angular momentum vector changes as a result of the torque, i.e.,

[TABLE]

It is important to note that equations (29) and (31) assume that the total mass $M$ is constant. This is the case for conservative mass transfer ($\dot{M}=0$), but not necessarily for nonconservative mass transfer (the latter is beyond the scope of this paper).

The first term of equation (28) is proportional to the relative velocity, $\dot{\mathbf{r}}$. For a perturbation of the form $\tilde{\mathbf{f}}=C\,\dot{\mathbf{r}}$, where $C$ is a scalar quantity that can depend on time/orbital phase, the above expressions imply

[TABLE]

where $\hat{\mathbf{q}}\equiv\hat{\mathbf{h}}\times\hat{\mathbf{e}}$. Equation (33b) implies a scalar eccentricity change of

[TABLE]

and the argument of periapsis $\omega$ changes according to

[TABLE]

Under the influence of $\tilde{\mathbf{f}}=C\,\dot{\mathbf{r}}$, the inclination $i$ and longitude of the ascending node $\Omega$ remain constant, since $\dot{\mathbf{h}}=C\,\mathbf{r}\times\dot{\mathbf{r}}=C\,\mathbf{h}$, implying that $\mathbf{h}$ does not change its direction.

More generally, for any phase-dependent perturbation $\tilde{\mathbf{f}}$, the orbital elements change according to

[TABLE]

We assume that $\mathbf{r}_{\mathrm{A}}$ in the second and third terms of equation (28) is directed along $\hat{\mathbf{r}}$. Therefore, the second term in equation (28) is $\propto\hat{\mathbf{h}}\times\hat{\mathbf{r}}$, implying that the associated $\dot{\mathbf{h}}\propto\mathbf{r}\times\left(\hat{\mathbf{h}}\times\hat{\mathbf{r}}\right)=r\,\hat{\mathbf{h}}$, i.e., $\mathbf{h}$ does not change direction as a result of the second term. Lastly, due to the third term in equation (28) $\dot{\mathbf{h}}\propto\mathbf{r}\times\hat{\mathbf{r}}=\mathbf{0}$, i.e., $\mathbf{h}$ does not change at all as a result of the third term. We conclude that, for our assumed perturbation in equation (28), $i$ and $\Omega$ remain constant.

Next, we orbit-average the remaining nontrivial equations for $\dot{a}$, $\dot{e}$ and $\dot{\omega}$. We define orbit-averaged quantities in the usual way and formulated in terms of the eccentric anomaly, i.e.,

[TABLE]

where $(...)$ denotes the quantity to be averaged. In equation (37), we assume that the orbital elements (most importantly, $a$ and $e$) are constant during the orbit (adiabatic approximation). In our model, the mass transfer rate is zero for $\mathcal{E}<-\mathcal{E}_{0}$, and $\mathcal{E}>\mathcal{E}_{0}$. Therefore, in practice, the range of the integral in equation (37) is taken to be $-\mathcal{E}_{0}<\mathcal{E}<\mathcal{E}_{0}$. Note that $\mathcal{E}_{0}$ is a function of $e$ and $x$ (see equation 20).

Normalization

As mentioned in Section 2.3, we relate the quantity $\dot{M_{\mathrm{d}}}_{,0}$ in equation (23) to the orbit-averaged mass transfer rate $\langle\dot{M}_{\mathrm{d}}\rangle$, and assume that $\langle\dot{M}_{\mathrm{d}}\rangle$ is known. For $\tau=0$, equations (23) and (37) imply

[TABLE]

where the dimensionless function $f_{\dot{M}}(e,x)$ is given explicitly by equation (B.1) in Appendix B. When $\tau\neq 0$, the dimensionless function should be replaced with $f_{\dot{M}}(e,x,\mathcal{E}_{\tau})$, which is given explicitly by equation (B.1) in Appendix B.

We remark that $\langle\dot{M}_{\mathrm{d}}\rangle$ is assumed to be known and constant during the orbit. However, $\langle\dot{M}_{\mathrm{d}}\rangle$ can change owing to the changing structure of the donor as a result of mass transfer, and/or owing to stellar evolution. For simplicity and to separate the orbital evolution from the donor’s structure and stellar evolution, we assume that $\langle\dot{M}_{\mathrm{d}}\rangle$ is constant in our applications (Sections 3 and 4). However, in other situations, e.g., when modeling the long-term evolution in population synthesis studies, $\langle\dot{M}_{\mathrm{d}}\rangle$ should be calculated self-consistently and allowed to vary over timescales much longer than the orbital timescale.

2.4.2 Secular change of the orbital elements with $\tau=0$ — case (1): negligible donor spin

Setting $\tau=0$ and in case (1) of Section 2.2 for $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ (negligible donor spin), the orbit-averaged equations of motion are given by

[TABLE]

The dimensionless quantities $g_{a}(e,x)$, $h_{a}(e,x)$, $g_{e}(e,x)$, and $h_{e}(e,x)$ are closed-form analytic functions given explicitly in Appendix B. In presenting these functions, we keep the explicit dependence on $\mathcal{E}_{0}=\mathcal{E}_{0}(e,x)$ (see equation 20) which appears through the integration limits. In our practical implementation, $\mathcal{E}_{0}$ is replaced using equation (20) instead of carrying out this replacement analytically. This approach turns out to be numerically favorable, since situations can otherwise occur in which the ODE integrator evaluates the functions at illegal combinations of $e$ and $x$ (i.e., such that $\cos\mathcal{E}_{0}$ does not lie between -1 and 1).

Note that there is no secular change of the argument of periapsis, $\omega$. This is no longer the case when $\tau\neq 0$ (see Section 2.4.4).

2.4.3 Secular change of the orbital elements with $\tau=0$ — case (2): large mass ratio

Setting $\tau=0$ and in case (2) of Section 2.2 for $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ (large mass ratio $q$), the orbit-averaged equations of motion are given by

[TABLE]

2.4.4 Secular change of the orbital elements with $\tau\neq 0$ and zero ejection/accretion radii

As mentioned in Section 2.3, the terms associated with $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ and $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ are excessively complicated if $\tau\neq 0$. Therefore, we here restrict to the secular orbital evolution associated with the term $\propto\dot{\mathbf{r}}$ in equation (28) only, i.e., setting $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{0}$. In this case, the orbit-averaged equations of motion are given by

[TABLE]

The functions $f_{a}(e,x,\mathcal{E}_{\tau})$, $f_{e}(e,x,\mathcal{E}_{\tau})$, and $f_{\omega}(e,x,\mathcal{E}_{\tau})$ are given explicitly in Appendix B. We recall that $\mathcal{E}_{\tau}$ is a function of $\tau$, $n$, and $e$ (see equation 26).

2.5. Properties of the orbit-averaged equations of motion

Although the general expressions for the dimensionless functions in the orbit-averaged equations of motion are cumbersome (see Appendix B), some insight can be gained in limiting cases. Here, we consider the functions $f_{a}(e,x,\mathcal{E}_{\tau})$, $f_{e}(e,x,\mathcal{E}_{\tau})$, and $f_{\omega}(e,x,\mathcal{E}_{\tau})$, all normalized to $f_{\dot{M}}(e,x,\mathcal{E}_{\tau})$.

In the limit of $\mathcal{E}_{0}\rightarrow\pi$, RLOF occurs during the entire orbit. In this case, setting $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ and $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$ to zero and to second order in the eccentricity, the equations of motion are given by (assuming $x\neq 1$)

[TABLE]

As the orbit circularizes ($e\rightarrow 0$), equation (42a) reduces to the canonical relation for conservative mass transfer in circular orbits, equation (2). Furthermore, equation (42b) states that, as $e\rightarrow 0$, $\langle\dot{e}\rangle\propto e$, i.e., $e$ decays exponentially to zero. In contrast, according to the delta function model of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b), for low eccentricities $\langle\dot{e}\rangle\propto\sqrt{1-e^{2}}(1-e)\approx 1-e$, i.e., $\langle\dot{e}\rangle$ is nonzero at $e=0$, and this can lead to undesirable properties (see Section 3 below for explicit examples). This is a consequence of the assumption of a delta function mass transfer rate at periapsis, whereas periapsis is not defined for circular orbits.

Lastly, note that $\langle\dot{\omega}\rangle\propto\sin\mathcal{E}_{\tau}$ vanishes for $\mathcal{E}_{\tau}=0$ ($\tau=0$), i.e., the orbit only secularly precesses if mass is transferred asymmetrically with respect to periapsis.

2.6. Numerical implementation

To numerically integrate the orbit-averaged equations of motion as presented in Sections 2.4.2, 2.4.3, and 2.4.4, we implemented them into a code called emt (Eccentric Mass Transfer) which is freely available at https://github.com/hamers/emt. The code is an easy-to-use and standalone Python script (using the standard NumPy and SciPy libraries). The functions, as given in Appendix B, are implemented in C and interfaced with Python using ctypes. These functions can therefore be easily ported to other codes such as binary or higher-order multiplicity population synthesis codes, which are typically written in C or similar languages. In addition, we implemented the secular equations of motion associated with a third body, using the standard expressions to quadrupole and octupole order (e.g., Harrington 1968; Ford et al. 2000; Naoz et al. 2013). We also included the first-order post-Newtonian (PN) terms in the inner orbit that give rise to orbital precession. However, we did not include higher-order PN terms or tidal effects (tidal dissipation and tidal bulges).

Below, in Sections 3 and 4, we use the emt code to investigate mass transfer in binary and triple systems, respectively.

3. Applications: mass transfer in circular and eccentric binaries

We apply our model to a number of systems using the numerical implementation described in Section 2.6. In Section 3, we restrict to the case of isolated binaries and concentrate on comparing our results to those of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b). In Section 4, we focus on hierarchical triple systems undergoing LK oscillations.

3.1. Circular orbit

We first consider the simplest case of mass transfer in a circular binary, also setting $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{0}$ (point masses). As mentioned above, the expected orbital evolution in this case is described by equation (2). From the fact that, for a circular orbit and conservative mass transfer, $L_{\mathrm{orb}}^{2}=\mu^{2}GMa$ and $M$ are constant, it immediately follows that $M_{\mathrm{d}}^{2}M_{\mathrm{a}}^{2}a$ is conserved during mass transfer.

In our numerical example, we set the initial parameters to $M_{\mathrm{d}}=1\,\mathrm{M}_{\odot}$, $M_{\mathrm{a}}=0.8\,\mathrm{M}_{\odot}$, $\langle\dot{M}_{\mathrm{d}}\rangle=-10^{-8}\,\mathrm{M}_{\odot}\,\mathrm{yr^{-1}}$, $a=1\,\,\textsc{au}$, and $R=1000\,\mathrm{R}_{\odot}$, i.e., RLOF triggered in a relatively tight binary when a solar-type primary evolves to a giant. We emphasize that our parameters are chosen for illustration purposes only, and may not be entirely realistic. Nonetheless, in this case the qualitative behavior of mass transfer is independent of the exact choice of parameters.

We show the evolution of the semimajor axis and eccentricity in the top and middle panels of Fig. 3, respectively. The bottom panel shows the mass ratio $q$ and $\mathcal{E}_{0}$ (if applicable) as a function of time. We include results according to our model (labeled ‘emt’ in the figure), and according to the equations of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b) (labeled ‘Sepinsky’). The latter equations, for zero ejection/accretion radii, are of the form of equations (41), with $f_{\dot{M}}=1$, $f_{a}=\sqrt{1-e^{2}}$, $f_{e}=\sqrt{1-e^{2}}(1-e)$, and $f_{\omega}=0$. The canonical analytic expectation, $M_{\mathrm{d}}^{2}M_{\mathrm{a}}^{2}a$ is constant, is represented in the top panel of Fig. 3 with the red dotted line.

In our example in Fig. 3, $q>1$ initially, and the orbit shrinks accordingly. After $\approx 10\,\mathrm{Myr}$, $q$ reaches unity (see the black horizontal dotted line in the bottom panel of Fig. 3). Subsequently, the orbit expands. The middle panel illustrates the issue with the equations of motion of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b), as mentioned in the Introduction: in the ‘Sepinsky’ model, the eccentricity becomes negative, and subsequently grows to $\gtrsim 0.6$ by 70 Myr. Since $\langle\dot{a}\rangle$ depends on $e$, this also slightly affects the evolution of $a$ after $\approx 50\,\mathrm{Myr}$. In contrast, the eccentricity in the ‘emt’ model remains zero the entire time, and the semimajor axis evolution is in accordance with conservation of $M_{\mathrm{d}}^{2}M_{\mathrm{a}}^{2}a$.

Note that, as shown in the bottom panel of Fig. 3, $\mathcal{E}_{0}=\pi$ at all times in the ‘emt’ model. This reflects the fact that the donor is transferring mass in a circular orbit during all orbital phases (initially, $x\simeq 0.086$).

We emphasize that $\langle\dot{e}\rangle<0$ for $e=0$ according to the Sepinsky model if $q>1$. If $q<1$, then the Sepinsky model gives $\langle\dot{e}\rangle>0$, but this still gives evolution inconsistent with the canonical equation (2) for circular orbits. We show an example in Fig. 4, in which the same system is taken as in Fig. 3, except that the initial masses are now $M_{\mathrm{d}}=1.5\,\mathrm{M}_{\odot}$ and $M_{\mathrm{a}}=2\,\mathrm{M}_{\odot}$, such that initially $q=0.75<1$. In this case, the Sepinsky model predicts a growth of eccentricity, with $e$ increasing to $\sim 0.7$ by 100 Myr. Similarly to the case $q>1$, the semimajor axis evolution according to the Sepinsky model deviates from equation (2). The ‘emt’ model yields zero eccentricity at all times, with the semimajor axis increasing in accordance with the canonical expectation.

3.2. Eccentric orbit

3.2.1 Zero ejection/accretion radii

Next, we consider an eccentric orbit, still setting $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{0}$. The initial parameters are $M_{\mathrm{d}}=8\,\mathrm{M}_{\odot}$, $M_{\mathrm{a}}=1.4\,\mathrm{M}_{\odot}$, $\langle\dot{M}_{\mathrm{d}}\rangle=-10^{-8}\,\mathrm{M}_{\odot}\,\mathrm{yr^{-1}}$, $a=1\,\,\textsc{au}$, $e=0.92$, and $R=10\,\mathrm{R}_{\odot}$. This can be representative of a binary in which the primary star underwent a supernova explosion, leaving a neutron star (the accretor). In our scenario, the sudden mass loss and velocity kick produced an eccentric (but still bound) orbit, triggering RLOF of the companion star (the donor). The chosen initial parameters are just sufficient to trigger RLOF near periapsis in our model, i.e., $\mathcal{E}_{0}\simeq 0.13>0$.

The orbital evolution is shown in Fig. 5. According to the ‘Sepinsky’ model, $a$ decreases to $\simeq 0.5\,\,\textsc{au}$ after $\approx 400\,\mathrm{Myr}$, whereas the eccentricity decreases slightly. In contrast, according to the ‘emt’ model, the semimajor axis decreases significantly, and the orbit circularizes. Gradually, $\mathcal{E}_{0}$ increases, until around 60 Myr, $\mathcal{E}_{0}=\pi$, and RLOF occurs during the entire orbit. At $t\approx 200\,\mathrm{Myr}$, the periapsis distance $r_{\mathrm{p}}=a(1-e)$ reaches $r_{\mathrm{p}}\leq R$, in which case we no longer integrate the equations of motion (subsequently, all quantities are taken to be constant). In reality, strong interactions would likely play an important role well before $r_{\mathrm{p}}\leq R$ is reached. This shows an important qualitative difference between the ‘Sepinsky’ and ‘emt’ models: the ‘emt’ model predicts a strong interaction to occur such as a merger, collision, or strong tidal effects (the latter are not the focus here and so are not included), whereas according to the ‘Sepinsky’ model, no such strong interactions occur.

In the middle panel of Fig. 5, we also show with blue and brown lines the eccentricity evolution for two modified cases with smaller initial eccentricities, and where $a$ was decreased as well to ensure RLOF still occurs during (part of) the orbit. These cases correspond to $a=0.1\,\,\textsc{au}$, $e=0.45$ ($a=0.08\,\,\textsc{au}$, $e=0.2$) for the blue and (thinner) brown lines, respectively. In the ‘emt’ model, the eccentricity and semimajor axis decrease smoothly until $r_{\mathrm{p}}\leq R$. In the ‘Sepinsky’ model, the decrease in $a$ is smaller, and no collision occurs initially. The eccentricity decreases until $e=0$, after which $e<0$. The integration is halted (constant $e$) when $r_{\mathrm{p}}\leq R$.

These examples show that our model deviates significantly from the ‘Sepinsky’ model not only for circular orbits, but for more eccentric orbits as well. The latter can be understood by noting that, even for eccentric orbits, mass transfer in the ‘emt’ model occurs during a finite range of orbital phases, whereas the ‘Sepinsky’ delta function mass transfer rate applies strictly in the limit $e\rightarrow 1$.

We show the possible effect of a nonzero $\tau$ in Fig. 6, where, for the same parameters as above (with $a=1\,\,\textsc{au}$ and $e=0.92$), we compare the cases $\tau=0$, and $\tau\neq 0$. The mass-loss delay time $\tau$ is expected to be of the order of the hydrostatic timescale $\tau_{\mathrm{hyd}}$ (e.g., van der Helm et al. 2016), and the latter can be estimated by (Kippenhahn et al., 2012) $\tau_{\mathrm{hyd}}\approx\sqrt{R^{3}/(GM_{\mathrm{d}})}\simeq 0.026\,\mathrm{d}$. Here, we set $\tau$ to $\tau=10\,\tau_{\mathrm{hyd}}$ (setting $\tau=\tau_{\mathrm{hyd}}$ did not yield noticeable differences between the cases $\tau=0$ and $\tau\neq 0$). As shown in Fig. 6, for a nonzero $\tau$, the orbit starts expanding again after significant shrinkage of the orbit, and the eccentricity increases. At $\approx 70\,\mathrm{Myr}$, $r_{\mathrm{p}}\leq R$, which is earlier compared to the case $\tau=0$. Also, note that the apsidal line has advanced by nearly $2\pi$ by the time at which $r_{\mathrm{p}}\leq R$ (see the blue dashed line in the bottom panel of Fig. 6).

3.2.2 Nonzero ejection/accretion radii

Next, we focus solely on the ‘emt’ model, and briefly investigate the effect of terms associated with nonzero ejection/accretion radii (i.e., extended bodies). For the eccentric system in Section 3.2.1 with $a=1\,\,\textsc{au}$ and $e=0.92$, we compare in Fig. 7 the case with $\mathbf{r}_{\mathrm{A}}=\mathbf{0}$ (zero ejection and accretion radii), to $\mathbf{r}_{\mathrm{A}}\neq\mathbf{0}$ (nonzero ejection and accretion radii). In the latter case, we include both a nonzero ejection and accretion radius, where we adopt case (1) from Section 2.2 for the ejection radius (i.e., slow donor spin), and set the accretion radius to $r_{\mathrm{A}_{\mathrm{a}}}=0.01\,\mathrm{R}_{\odot}$ (note that the effect of $r_{\mathrm{A}_{\mathrm{a}}}$ for our chosen system is very small, unless $r_{\mathrm{A}_{\mathrm{a}}}\gtrsim 1\,\mathrm{R}_{\odot}$).

As shown in Fig. 7, the differences between zero and nonzero $\mathbf{r}_{\mathrm{A}}$ are not very large. Nevertheless, with $\mathbf{r}_{\mathrm{A}}\neq\mathbf{0}$, the orbit shrinks less, and the collision criterion $r_{\mathrm{p}}\leq R$ is avoided. Instead, the system evolves to reach mass ratio reversal at $t\approx 330\,\mathrm{Myr}$, around which time the orbit starts expanding.

Lastly, we compare in Fig. 8, for the same system as above, the cases (1) and (2) for the ejection radius as described in Section 2.2, i.e., either low donor spin (case 1), or high mass ratio (case 2). In the high-$q$ case, the orbit expands more rapidly after mass ratio reversal. Consequently, mass transfer transitions to partial RLOF around 350 Myr, and the eccentricity starts increasing again. Similar behavior occurs in the low donor spin case, although at later times.

4. Applications: triple-star systems

As demonstrated in Section 3.2, our model already shows more complicated behavior compared to the canonical relation for mass transfer in circular orbits with $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}=\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}=\mathbf{0}$, i.e., equation (2). The latter relation was shown with the red dotted lines in the top panels of the figures. If a tertiary object is included, the realm of possibilities increases even more. Here, we show two illustrative cases, but we emphasize that the parameter space is large; a full investigation is beyond the scope of this paper. Specifically, we consider cases when mass transfer occurs quickly (Section 4.1), and slowly (Section 4.2). In both cases, we use the code as described in Section 2.6, in which we include the standard secular quadrupole- and octupole-order terms, as well as 1PN terms in the inner orbit. Apart from the 1PN terms, we here focus on mass transfer effects only, and do not include dissipative higher-order PN terms or tides.

4.1. Fast mass transfer

We choose a triple system with inner binary parameters $M_{\mathrm{d}}=1\,\mathrm{M}_{\odot}$, $M_{\mathrm{a}}=0.1\,\mathrm{M}_{\odot}$, $\langle\dot{M}_{\mathrm{d}}\rangle=-10^{-8}\,\mathrm{M}_{\odot}\,\mathrm{yr^{-1}}$, $a=1\,\,\textsc{au}$, $e=0.001$, $\omega=45^{\circ}$, and $R=1\,\mathrm{R}_{\odot}$. The tertiary star has a mass $M_{\mathrm{t}}=1\,\mathrm{M}_{\odot}$; the outer orbit has parameters $a_{\mathrm{out}}=150\,\,\textsc{au}$, $e_{\mathrm{out}}=0.6$, $\omega_{\mathrm{out}}\simeq 5.7^{\circ}$, and the mutual inclination with respect to the inner orbit is $i_{\mathrm{rel}}=85^{\circ}$.

We carry out two integrations: one with mass transfer included ($\langle\dot{M}_{\mathrm{d}}\rangle=-10^{-8}\,\mathrm{M}_{\odot}\,\mathrm{yr^{-1}}$), and one (effectively) without mass transfer ($\langle\dot{M}_{\mathrm{d}}\rangle=-10^{-30}\,\mathrm{M}_{\odot}\,\mathrm{yr^{-1}}$). In Fig. 9, we compare the two cases. In the absence of mass transfer, the inner orbit eccentricity oscillates with quasi-regular LK oscillations with high amplitude. In the mass transfer case, mass transfer does not occur initially ($\mathcal{E}_{0}=0$), but after a few cycles, the eccentricity becomes high enough to trigger partial RLOF ($\mathcal{E}_{0}>0$). The semimajor axis quickly decreases (since $q>1$), and the orbit shrinks and circularizes. From the lack of further oscillations after the high-$e$ peak at $\approx 10\,\mathrm{Myr}$, it can be deduced that the inner binary quickly decouples from the tertiary after mass transfer has ensued.

This type of evolution is analogous to the ‘fast’ type of mergers occurring in black hole triple systems (e.g., Randall & Xianyu 2018), or stellar/planetary systems in which tidal friction is important (e.g., Petrovich 2015).

4.2. Slow mass transfer

In Fig. 10, we show the evolution of a very similar system to that in Section 4.1 – the only modification is that the average mass transfer rate has been decreased (in absolute value) by a factor of 10, i.e., $\langle\dot{M}_{\mathrm{d}}\rangle=-10^{-9}\,\mathrm{M}_{\odot}\,\mathrm{yr^{-1}}$. The resulting evolution is significantly different, however: the inner binary still shrinks and circularizes, but this process takes $\sim 10$ times longer. In particular, the inner binary does not immediately decouple from the tertiary star, as is illustrated with the insets in Fig. 10. Initially, around 10 Myr, mass transfer ensues in short bursts at which $\mathcal{E}_{0}$ increases to $\approx 0.05\,\mathrm{rad}$ (see the bottom panel of Fig. 10). These bursts are associated with stair-wise jumps in the semimajor axis, and gradually locking the eccentricity to a high value. After $\approx 15\,\mathrm{Myr}$, the LK oscillations stop and the orbit circularizes, although partial RLOF continues. After $\approx 60\,\mathrm{Myr}$, full RLOF occurs ($\mathcal{E}_{0}=\pi$ in the bottom panel of Fig. 10).

5. Discussion

In Section 2, we made a number of simplifying assumptions in order to arrive at an analytically tractable model. These simplifications include the neglect of the force of the mass transfer steam on the stars, conservative mass transfer, imposed ejection and accretion velocities, and specific assumed directions and magnitudes of the ejection/accretion points relative to the donor/accretor. Although most of these assumptions have at least some physical basis, their validity should be evaluated in future work using detailed hydrodynamic simulations. Of course, this endeavor is hampered by the fact that such simulations are typically limited to no more than $\mathcal{O}(10)$ orbits, whereas the secular evolution takes place on much longer timescales. Alternatively, $N$-body simulations tailored for mass transfer to compute the mass transfer stream trajectories and the effect on the orbit could be used to test the long-term evolution, similarly to, e.g., Sepinsky et al. (2010) and Dosopoulou et al. (2017).

Of highest importance is to evaluate the correctness of our assumption of the ejection/accretion velocities, $\mathbf{w}_{\mathrm{d}}$ and $\mathbf{w}_{\mathrm{a}}$, which leads to the (commonly used) expression equation (4) (ignoring back-reaction forces of the matter stream on the stars, and finite sizes of the ejection/accretion radii). More fundamentally, even the physical motivation behind the ‘canonical’ relation for conservative mass transfer in circular binaries, equation (2), has been called into question (Luk’yanov, 2008).

6. Conclusions

We presented an analytic model to describe the long-term orbital evolution of mass-transferring binaries. Our model applies to conservative mass transfer, in which no net mass is lost from the system. In contrast to previous works, our model is, in principle, applicable for any eccentricity $e$, and it gives qualitatively and quantitively different evolution for both circular and eccentric orbits. We implemented the model in a code that is publicly available (see Section 2.6 for the url), and we investigated several cases of Roche Lobe overflow mass transfer in circular and eccentric binaries, and in triples undergoing LK oscillations. Our main conclusions are given below.

1. The model of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b) assumed that the mass transfer rate is a delta function of the orbital phase, centered at periapsis. Although physically reasonable in the limit $e\rightarrow 1$, this assumption breaks down for smaller eccentricities, when mass transfer is expected to occur for a finite range of orbital phases. In the limit of a circular orbit (in which RLOF occurs during the entire orbit), the mass transfer rate is expected to be constant. As a consequence of the delta function assumption, the equations of motion of the model of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b) state that the eccentricity time derivative is negative at zero eccentricity (assuming point masses and a donor-to-accretor mass ratio $q>1$), which, when solved (numerically) as a function of time, yields unphysical negative eccentricities (see, e.g., the middle panel of Fig. 3). We showed that the delta function model is problematic in the case $q<1$ as well (see Fig. 4). We remedied these issues by assuming that the mass transfer rate is a smooth function of orbital phase, reducing to a delta-like function at high eccentricities, and a flatter function at lower eccentricities. In our model, the eccentricity decays exponentially to zero for small eccentricities (see Section 2.5).

2. We derived explicit expressions for the orbit-averaged equations of motion in our model. We also included the effect of nonzero ejection/accretion radii, making two limiting assumptions related to the location of the ejected material relative to the donor (see Section 2.2). For zero ejection and accretion radii, we also derived expressions for the orbital evolution when mass transfer occurs with a delay after periapsis. We implemented our expressions into an easy-to-use and freely available Python code to quickly solve the equations of motion numerically, and which can be used as a basis for implementations of our model into binary (and higher-order multiplicity) population synthesis codes.

3. We applied our model to circular and eccentric binaries. We showed that, in the limit of circular orbits, our model is in agreement with the canonical relation describing conservative mass transfer in circular orbits. In contrast, in our example with a circular system, the model of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b) yielded an unphysical negative eccentricity, and different semimajor axis evolution compared to the canonical relation (the latter applied to both cases $q>1$, and $q<1$). In examples of eccentric binaries, we found that our model predicts faster orbital shrinkage and circularization compared to the model of Sepinsky et al. (2007b) and Dosopoulou & Kalogera (2016b). Also, assuming a delay time between periapsis passage and mass transfer resulted in orbital expansion and eccentricity driving. Furthermore, we found that effects associated with nonzero ejection and accretion radii are somewhat, but not extremely, important. However, we emphasize that we considered only a few examples. A comprehensive parameter space study is beyond the scope of this work.

4. The implementation of our model also allows for inclusion of the (secular) perturbation by a distant orbiting third body. We illustrated the process of mass transfer in triples undergoing LK oscillations in two examples, and showed that mass transfer can act ‘quickly’ (i.e., mass transfer effects are rapidly decoupled from LK oscillations), and ‘slowly’ (in which case many oscillations can occur before the inner binary becomes decoupled). This dichotomy is analogous to triple systems in which dissipative effects are important, such as tidal evolution in systems with planets and/or stars, and gravitational wave emission in systems containing compact objects.

Acknowledgements

We thank Scott Tremaine for insightful discussions and comments on the manuscript. We also thank Dimitri Veras for feedback, and the anonymous referee for helpful comments. A.S.H. gratefully acknowledges support from the Institute for Advanced Study, the Peter Svennilson Membership, and the Martin A. and Helen Chooljian Membership. F.D. acknowledges support from PCTS and Lyman Spitzer Jr fellowships.

Appendix A A. Explicit expression for the first Lagrange point in the limit of small spin

If $\mathcal{A}=0$ is set in equation (8), then the (scaled) first Lagrange point $X_{\mathrm{L}}$ has an analytic solution of $q$, given by

[TABLE]

Here,

[TABLE]

Appendix B B. Explicit expressions for the functions appearing in the orbit-averaged equations of motion

Here, we give explicit expressions for the various functions that describe the orbit-averaged equations of motion, and which are referred to in Section 2.4.

B.1. Normalization

The normalization function is obtained by orbit-averaging $\dot{M}_{\mathrm{d}}$, and reads, for $\tau=0$,

[TABLE]

For $\tau\neq 0$, the normalization function reads

[TABLE]

B.2. Terms associated with zero ejection/accretion radii

For $\tau=0$, the functions associated with the terms in $\bf{f}$ proportional to $\dot{\mathbf{r}}$, which appear regardless of the assumed ejection/accretion locations $\mathbf{r}_{\mathrm{A}_{\mathrm{d}}}$ and $\mathbf{r}_{\mathrm{A}_{\mathrm{a}}}$, are given by

[TABLE]

For $\tau\neq 0$, these functions are

[TABLE]

B.3. Terms associated with nonzero ejection/accretion radii

The following functions apply only to the case $\tau=0$. The semimajor axis functions associated with the terms involving nonzero ejection/accretion radii are given by

[TABLE]

Lastly, the eccentricity functions associated with the terms involving nonzero ejection/accretion radii are given by

[TABLE]