Orbital evolution of mass-transferring eccentric binary systems

# Orbital evolution of mass-transferring eccentric binary systems. I. Phase-dependent evolution

Fani Dosopoulou and Vicky Kalogera Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astrophysics, Northwestern University, Evanston, IL 60208
###### Abstract

Observations reveal that mass-transferring binary systems may have non-zero orbital eccentricities. The time-evolution of the orbital semi-major axis and eccentricity of mass-transferring eccentric binary systems is an important part of binary evolution theory and has been widely studied. However, various different approaches and assumptions on the subject have made the literature difficult to comprehend and comparisons between different orbital element time-evolution equations not easy to make. Consequently, no self-consistent treatment of this phase has been ever included in binary population synthesis codes. In this paper, we present a general formalism to derive the time-evolution equations of the binary orbital elements, treating mass-loss and mass-transfer as perturbations to the general two-body problem. We present the self-consistent form of the perturbing acceleration and the phase-dependent time-evolution equations for the orbital elements under different mass-loss/transfer processes. First, we study the cases of isotropic and anisotropic wind mass-loss. Then, we proceed with the non-isotropic ejection and accretion in conservative as well as non-conservative manner for both point masses and extended bodies. Comparison of the derived equations with similar work in the literature is made and explanation of the existing discrepancies is provided.

## 1. Introduction

In modeling of binary populations, most interacting mass-transferring binary systems are assumed to be circular (e.g. Belczynski et al., 2008). However, observations reveal that semi-detached binary systems might have non-zero orbital eccentricities (e.g., Petrova & Orlov, 1999; Bonačić Marinović et al., 2008; Vos et al., 2013; Boffin et al., 2014). In addition, recent observations indicate that High-Mass X-ray Binaries (HMXB) are eccentric, with relatively high eccentricities in the case of Be/X-ray binaries (e.g., Raguzova & Popov, 2005; Walter et al., 2015). The residual eccentricity that has been observed in mass-transferring binary systems is quite surprising, since tidal dissipation might be expected to circularize the orbit and drive the star towards co-rotation (e.g., Hut, 1981; Eggleton et al., 1998; Bonačić Marinović et al., 2008). This implies that the efficiency of tidal circularization in interacting eccentric binaries is not as high as typically assumed. Previous work has shown that especially for long-period systems and solar-type binaries in open clusters tides are not sufficient to circularize the orbit, leaving the systems with a considerable non-zero eccentricity (e.g., Vos et al., 2013, 2012, 2015; Meibom & Mathieu, 2005). Possible mechanisms that have been proposed to help the system retain a finite eccentricity are mass-loss and mass-transfer processes.

Mass-loss/transfer can indeed enhance eccentricity (e.g., Vos et al., 2015; Bonačić Marinović et al., 2008; Soker, 2000) and is a common evolutionary phase in a binary system, being responsible for many observed astrophysical phenomena such as Type Ia supernovae, X-ray binaries, neutron-stars spin up as well as orbital contraction or expansion. Mass can be ejected or accreted in a binary system in many different and complicated ways and is phase-dependent in an eccentric orbit. Commonly studied cases include isotropic wind mass-loss (e.g., Hadjidemetriou, 1963), Roche-Lobe-Overflow (RLOF) (e.g., Sepinsky et al., 2007, 2009) and Bondi-Hoyle accretion (e.g., Bondi & Hoyle, 1944; Boffin & Jorissen, 1988; Hurley et al., 2002). All the aforementioned processes change the total mass, energy and orbital angular momentum of the binary, which corresponds to a change in the binary orbital elements. The effect of mass-transfer on the orbital elements of the orbit goes back a while ago (e.g., Hadjidemetriou, 1963; Huang, 1956) and continues to recent years with work studying various and different cases of mass-loss and mass-transfer (e.g., Sepinsky et al., 2007, 2009, 2010; Veras et al., 2011, 2013).

In the existing literature, to our knowledge, the time-evolution of the orbital elements of a binary undergoing mass-loss and mass-transfer has been derived following two different approaches. One involves calculating time variations of the orbital elements based on changes in the total orbital energy and angular momentum of the system (e.g., Huang, 1956; Boffin & Jorissen, 1988; Bonačić Marinović et al., 2008). The other is treating mass-loss/transfer as a perturbation in the two-body problem and using the latter to calculate the time-evolution of the orbital elements (e.g., Hadjidemetriou, 1963; Sepinsky et al., 2009; Veras et al., 2011, 2013). However, consistency between these two different approaches is not always clear and many assumptions that have implicitly been made following either technique are not always clearly stated. Thus, having a general mathematical formalism to describe the orbital evolution of a binary system undergoing mass-transfer of various types is both useful and important in binary evolution theory. In the first part of this paper, we present how to derive the time-evolution equations of the orbital elements based on the Lagrange perturbation formalism (Lagrange ) and then apply this to specific cases of mass-loss and mass-transfer in eccentric binaries.

A binary system exposed to different kind of perturbations follows an orbit that is continuously evolving in time, i.e., an open orbit. This means that much attention must be paid to the physical interpretation of the orbital elements we assign to the evolving orbit (e.g., Veras et al., 2011; Efroimsky & Goldreich, 2004). In this work, we refer to the importance of this interpretation which is connected both to the freedom in the choice of orbital elements and the gauge freedom in orbital mechanics.

The paper is organized as follows. In Section 2 we present the general formalism we use to describe the perturbed two-body problem. In section 3 we comment on the importance of the choice of reference frame and present the time-evolution equations of the orbital elements in three different reference frames. In sections 4 and 5 we apply this formalism to the cases of isotropic and anisotropic wind mass-loss respectively. In Section 6 we study conservative as well as non-conservative non-isotropic mass ejection and accretion in the system for extended binary components. In Section 7 we present the form the time-evolution equations of the orbital elements reduce to under the point-mass approximation. Throughout all the aforementioned sections, comparison to previous work in the literature that has considered similar cases is made. Finally, an explanation of any existing discrepancies is provided. In Section 8 we summarize our results.

The perturbing acceleration is in principle phase-dependent. In this paper, we present the general form of the perturbing acceleration and the orbital element time-evolution equations, keeping this phase-dependence. In a follow-up paper Dosopoulou Kalogera 2016b (from now on Paper II), we orbit-average the phase-dependent time-evolution equations of the orbital elements relevant to different types of mass-loss/transfer in the binary under either the assumption of adiabaticity or for delta-function mass-loss/transfer at periastron.

## 2. The perturbed two-body problem

Two-body systems like binary stars can be exposed to gravitational perturbations from different agents. These include tidal dissipation, relativistic corrections, gravitational-wave radiation, magnetic fields, inertial forces, mass-loss/mass-transfer phenomena and other. All these physical processes act like perturbing forces to the general two-body problem and need a new mathematical way of treatment compared to the unperturbed case. Due to the various perturbations, each body is no longer moving in the actual physical trajectory (conic section) it would if no perturbations existed, but its physical orbit is slowly changing with time.

Either called Variation of constants in pure mathematical language or Varying-Conic method in purely geometric terms, the method that was proposed to study these systems was advanced and completed by Lagrange. This method involves an approximation of the true physical orbit by a family of evolving instantaneous conic sections.

In subsection , we will prove that the latter representation is not unique, thus one should be very meticulous when applying this mathematical formulation. Specifically, the actual physical orbit of the body is open, the orbital elements derived though this method are not the classic orbital elements of a closed orbit. Rather, they either can have a possible physical interpretation under certain assumptions or they should be treated just as mathematical tools otherwise.We shall discuss this issue in more detail in the next sections.

### 2.1. Variation of Constants method

The general reduced two-body problem is described by the equation

 ¨r=−GMr3r (1)

where is the gravitational constant, is the total mass of the system and is the relative position between the two bodies. The general solution to this problem is a Keplerian conic section defined by six constant orbital elements .

In a specified inertial coordinate system, this conic is described by

 r=r(C1,...,C6,t),˙r=(C1,...,C6,t) (2)

where by definition

 v=(∂r∂t)Ci=const. (3)

is the relative velocity between the two bodies in the system.

The perturbed two-body problem can then be written as

 ¨r=−GMr3r+F(r,˙r) (4)

where the perturbing force depends upon both the relative position and velocity.

Equation can be solved assuming that at each instant of time, the true orbit can be approximated by an instantaneous Keplerian Conic section which is changing over time through the now time-dependent orbital elements , i.e.,

 r=r(C1(t),...,C6(t),t) (5)

and the velocity of the body will now be given by

 ˙r=∂r∂t+6∑i=1∂r∂CidCidt. (6)

Equations and constitute a system of three independent equations, one independent variable (time ) and six dependent variables (orbital elements ). This means that the parametrization in terms of the orbital elements is not unique and additional constraints need to be applied.

#### 2.1.1 Lagrange constraint - Osculating elements

Imposing the Lagrange constraint

 6∑i=1∂r∂CidCidt=0 (7)

makes the orbital elements osculating. This means that the instantaneous conics are tangential to the perturbed actual physical orbit, as it can be seen by equation and is also depicted in Figure .

Another way to think about osculating orbits is that these are the orbits that the body would follow if perturbations were to cease instantaneously.

Equations , and constitute now a well-defined system of six equations and six variables.

Following the Lagrange constraint, the perturbed equation of motion can be written as

 ∂2r∂t2+GMr3r+6∑i=1∂v∂CidCidt=F (8)

where refers to the solution of the unperturbed two-body problem, i.e., satisfies equation .

This simplifies equation to

 6∑i=1∂v∂CidCidt=F (9)

which using the definition of Lagrange and Poisson brackets can be decoupled to read

 dCidt=−6∑j=1{CiCj}∂r∂CjF. (10)

Poisson brackets are defined by

 (11)

they are the negative inverse to Lagrange brackets

 [CkCi]=∂r∂Ck∂v∂Ci−∂r∂Ci∂v∂Ck, (12)

and do not depend on time explicitly.

Equations give the time-evolution equations of the six osculating orbital elements, which are invariant with respect to the change of coordinates and/or orbital parametrization.

We choose the six Keplerian elements , namely semi-major axis, eccentricity, inclination, argument of periapsis, longitude of the ascending node, and true anomaly . Also, we define as the mean motion.

Making use of the following Poisson brackets

 {σa} =−{aσ}=2na (13) {σe} =1−e2na2e,{ωe}=−√1−e2na2e (14) {σa} =−√1−e2na2e (15) {Ωi} =−1na2√1−e2sini (16) {ωi} =−cosina2√1−e2sini (17)

yields the following Lagrange time-evolution equations for the aforementioned six orbital elements

 dadt =2naF∂r∂σ (18) dedt =1−e2na2eF∂r∂σ−√1−e2na2eF∂r∂ω (19) didt (20) dωdt =√1−e2na2eF∂r∂e−cosina2√1−e2siniF∂r∂i. (21) dΩdt =1na2√1−e2siniF∂r∂i (22) dνdt =n(1+ecosν)2(1−e2)3/2−dωdt−cosidΩdt. (23)

We note here that the orbital elements appearing in equations - are, as we will explain in detail in the next sections, the so-called osculating elements. The physical meaning of these elements as well as their uniqueness with respect to any other choice of orbital elements will be explained in the next sections as well.

Equations - describe the time-evolution of all the six orbital elements. Because in this paper we are interested only in perturbations induced by mass-loss/transfer processes and for reasons we describe in detail in the last paragraph of the following section, in this paper we present only the time-evolution of the semi-major axis and the eccentricity [see equations and ]. However, as we can see from equations - for a general perturbation the orbit undergoes various precessions. At the end of the following section we discuss briefly these precessions as well. In Paper II we discuss in more detail the importance of these precessions within the context of secular evolution.

#### 2.1.2 Introducing Gauge freedom in orbital mechanics - Non-osculating elements

According to Efroimsky & Goldreich (2004), imposing the Lagrange constraint is not necessary and one is free to set

 6∑i=1∂r∂CidCidt=Φ (24)

where is called the gauge and represents the gauge freedom in orbital mechanics.

Generalizing more the formulation by writing the perturbing force in terms of a perturbing Lagrangian

 F=∂ΔL∂r−ddt(∂ΔL∂˙r) (25)

one can re-write equations and as follows (e.g., Efroimsky & Goldreich, 2004)

where is the perturbing Hamiltonian defined by

 ΔH=−[ΔL+12(∂ΔL∂˙r)2] (28)

and functions and have the same functional form as the unperturbed-case position and velocity vector but are now time dependent through the time-dependent orbital elements . We note here that equation seems to have a singularity in the case of circular orbits (). This singularity is not real since to take the limit of equation for one needs first to calculate explicitly the form of the partial derivatives included in this equation as functions of the orbital elements.

For a general gauge they are connected to the perturbed velocity and position by the relations

 r =f(C1,...,C6,t) (29) ˙r =∂f∂t+Φ (30) =g+Φ. (31)

Equations and constitute a physical interpretation of the gauge freedom.

Choosing the Lagrange constraint , one enforces the family of the approximating conics to be always tangent to the physical orbit of the body. In other words, this means that the orbital elements under this constraint are osculating elements and describe the conic the body would follow if the perturbation ceased instantaneously.

Any other non-zero constraint would lead to orbital elements that do not describe conics that are tangential to the real orbit, but differ from the tangent velocity vector by [see Figure 2]. These elements are the so-called non-osculating elements and can be used in some cases as useful mathematical tools to describe the problem (e.g., Efroimsky, 2002).

It is obvious that a real physical interpretation can be assigned only to the osculating elements (Lagrange constraint), since only those describe conics that are always and at any time tangent to the actual path of the body.

We should point out here that even the osculating orbital elements derived by this formulation are mathematical tools that can be used to approximately describe the real trajectory of the body. This means that one should be very meticulous when interpreting them. The physical time-evolution of the system is precisely expressed with the time-derivative of the position vector of the body .

When applying equations and or alternatively and to any kind of problem, one needs to know the various partial derivatives of and with respect to the orbital elements. For this reason, we calculated and we present in Appendix A all the needed mathematical quantities as functions of the true anomaly .

Before proceeding to the application of the zero-gauge and gauge-free aforementioned formulations to specific physical problems, we would like to make a comment on the usefulness and advantages of each of them.

Based on what we presented here there are two ways one can derive the time-evolution equations of the orbital elements. One is using the perturbing force directly applying the osculating equations and and the other is using the equivalent Lagrangian defined by equation and applying equations and . Although these two different approaches are completely equivalent eventually leading to the same equations for the osculating orbital elements, the second method which makes use of the Lagrangian is more general because it introduces the gauge freedom. This means that the latter is the only method one could use to choose a gauge different than zero. There are cases where the use of a gauge other than zero is needed, e.g., when the Delaunay variables need to remain canonical for a general perturbation depending both on positions and velocities of the bodies (e.g., Efroimsky & Goldreich, 2003, 2004).

In the next sections, we will apply the gauge-free method to derive and compare the orbital element time-evolution equations between a zero-gauge and a different gauge , called the generalized Lagrange constraint, which is commonly used in the literature (e.g., Efroimsky, 2002; Gurfil & Belyanin, 2008) and is defined by

 Φ=−∂ΔL∂˙r. (32)

Using this gauge in equations and gives the time-evolution equations of the so-called contact orbital elements. As mentioned above, although equations and are simplified by using this gauge, the derived contact orbital elements have no physical meaning compared to the zero-gauge and the relevant osculating orbital elements.

## 3. The importance of reference frame choice

It is important to note that equations and are in vector form and thus independent of the reference frame used in each specific problem. After decomposing the perturbing force in these equations into its components in a chosen reference frame, one should be careful in interpreting the resulting equations for the semi-major axis and the eccentricity. In this section we discuss the importance of the reference frame choice pointing out that a vanishing component of the perturbing force in one specific reference frame does not guarantee any other vanishing component of the perturbing force in a different reference frame [see the rotation matrices -].

The most commonly used reference frames to describe a two-body problem are depicted in Figure . Reference system is the inertial reference system. In the reference system , the unit vector is along the eccentricity vector and points always into the direction of the periastron (rotated by the two angles and relative to ). Finally, in the reference system , the unit vector is along the relative position vector between the two bodies in the binary (rotated relative to by the angle of true anomaly ).

Acceleration can be written in these reference frames as follows

 A =Ax^x+Ay^y+Az^z (33) A =Ap^p+Aq^q+Al^l (34) A =Ar^r+Aτ^τ+An^z. (35)

The form of equations and will be different based on the reference frame chosen each time in the problem. The use of the correct set of equations for the chosen reference frame is important since using the incorrect set of equations can lead to misinterpretation of the system’s orbital evolution.Here we present this set of equations for three different refernce frames.

For decomposition in the reference system equations and take the form

where the different are defined by

 C1 =ecosω+cos(ν+ω) (38) C2 =esinω+sin(ν+ω) (39) C5 =(2ecosν+1+cos2ν)sinω (40) +2(e+cosν)cosωsinν (41) C6 =(2ecosν+1+cos2ν)cosω (42) −2(e+cosν)sinωsinν. (43)

For decomposition in the reference system equations and become

For decomposition in the reference system equations and yield

Transformation between the different reference systems is done using the appropriate rotation matrices. Specifically, the cyclic transformation is done using the relations

 QTKI→KJ (48) RKJ→KR (49) QTKI→KR,ω→ω+ν (50)

where the right arrow in equation means substituting with in the rotation matrix . The rotation matrices and are defined by

 Q=⎛⎜⎝Q11Q12Q13Q21Q22Q23Q31Q32Q33⎞⎟⎠ (51)

and

 Q11 =cosΩcosω−sinΩsinωcosi (52) Q12 =−cosΩsinω−sinΩcosωcosi (53) Q13 =sinΩsini (54) Q21 =sinΩcosω−cosΩsinωcosi (55) Q22 =−sinΩsinω−cosΩcosωcosi (56) Q23 =−cosΩsini (57) Q31 =sinωsini (58) Q32 =cosωsini (59) Q33 =cosi (60)

and

 R=⎛⎜⎝cosνsinν0−sinνcosν0001⎞⎟⎠. (61)

In an attempt to clarify the meaning of equations and we mention here that the latter do not refer to inertial reference frames. The reference frames and are not inertial since they are rotating. Equations and constitute a decomposition of the perturbing force into its components in the relevant reference frames and are true at any instant of time. The matrices and are time-dependent since the orbital elements are time-dependent. Thus, the transformations between the different reference frame components of the force are true at any instant of time. These transformations are constructed in such a way that they include the time-derivative of the base unit vectors in and due to the various precessions (). Consequently, the various inertial forces that appear when we transform into a rotating reference frame, namely and , are already embedded in the form of the time-dependent rotation matrices and . Inertial forces like Coriolis acceleration, centripetal acceleration and the Euler acceleration will appear applying the relevant time-dependent rotation matrix to the position and force vectors and substitute the transformed position and force vectors into the form of the equation of motion in the inertial system. It is after this substitution that the inertial forces appear in their known form.

We note here that Veras et al. (2013) presented equations and as well as and and were the first to apply them to study the case of anisotropic wind mass-loss [see Section 5 here and Section 5 in Paper II]. However, choosing the reference system they incorrectly used in their equations the components of the perturbation along the inertial Cartesian system instead of the components of the perturbation in the reference system [compare to our equations and ]. As previously mentioned, this is important since the incorrect use of the right set of time-evolution equations of the orbital elements for the chosen reference system can lead to a misinterpretation of the evolution changes in the semi-major axis and eccentricity .

As we mentioned before in this paper we are only interested in the time-evolution of the semi-major axis and the eccentricity. The basic reason for this is the type of perturbations induced by mass-loss/transfer processes. To explain this in more detail we choose to decompose equations in the reference frame

 dadt =2n√1−e2[Fresinν+Fτ(1+ecosν)] (62) dedt =√1−e2na[Frsinν+Fτ(cosν+e+cosν1+ecosν)] (63) didt =Fncos(f+ω)√1−e2na(1+ecosf) (64) dΩdt =1na2√1−e2siniFna(1−e2)sin(ν+ω)1+ecosν (65) dωdt =√1−e2nae(−Frcosν+Fτsinν2+ecosν1+ecosν) −cosinasiniFn√1−e2sin(ν+ω)1+ecosν (66) dνdt =n(1+ecosf)2(1−e2)3/2−dωdt−cosidΩdt (67)

As we will see in the following sections the perturbations induced by the mass-loss/transfer processes studied in this paper have the general form . From equations and we can see that the inclination and the longitude of the ascending node evolve only for a non-zero vertical to the orbital plane component of the perturbing force, i.e., . However, perturbations of the form do not contain such a vertical perturbing force component (i.e., ) and thus and do not evolve in the cases we study here. On the contrary, from equation we see that the argument of periapsis is indeed precessing due to mass-loss/transfer but we assume this precession is at first order approximation small (i.e., ). Under these considerations, we choose to focus and present only the time-evolution of the semi-major axis and the eccentricity .

## 4. Isotropic wind mass-loss

The components of a binary system can lose mass in many different ways. Due to the energy and angular momentum conservation laws mass-loss procedures expose the system to additional perturbations including linear and angular momentum recoils. One characteristic way through which stars can lose mass is stellar winds. However, the actual form and structure of these winds is a complicated phenomenon and depends strongly on the structure and rotation of the star as well as on the potential existence of a dynamically important magnetic field in the star. In many cases, e.g., in a supernova explosion, the mass-loss can be assumed to be spherically symmetric. This means that the wind velocity is the same for all points along the surface of the star and thus we refer to this case of wind mass-loss as isotropic. We note here that isotropic mass-loss does not produce any momentum kick on the mass-losing body. However, the total mass of the system changes with time and this induces in the system a perturbation which acts as we will see below as a linear drag force.

As mentioned above, in the case of isotropic wind mass-loss the perturbation in equation comes implicitly from the fact that the mass of one or both stars is changing leading to a time-dependent total mass of the system. Figure 4 describes schematically the case of isotropic mass-loss considered in this section. Let’s assume that at some point in time the two-body system has a total mass and the body has exactly the needed velocity to follow the unperturbed orbit (black ellipse in Figure 4). If the total mass of the system momentarily changes by introducing no additional momentum kicks to the bodies (what we refer to as isotropic) then the period increases and the velocity is more than needed to follow the unperturbed orbit. The body will follow the perturbed orbit instead (red ellipse in FIgure 4). As we will see below we can describe this perturbation as a linear drag force acting to the body , forcing it to follow the osculating orbit. Before deriving the actual form of the perturbation we remind here that the perturbed orbit is an osculating orbit. This means that the body will never follow the osculating orbit since after time the total mass of the system will continue changing and the actual physical orbit of the body will never be a closed orbit in time. From equations - we know that because of the perturbation induced, the osculating orbit will be both evolved and precessed relative to the unperturbed orbit. In Paper II, we will prove that the osculating orbit as already indicated in Figure 4 has on average the same eccentricity, a greater semi-major axis and it precesses compared to the perturbed orbit.

### 4.1. Perturbing force treatment

Hadjidemetriou (1963) and Omarov (1964) independently derived the appropriate perturbation for isotropic wind mass-loss. For a donor of mass and a companion of mass , if the total mass of the two-body system is changing in time with a rate , making use of the energy equation leads to the following perturbation for isotropic wind mass-loss (Hadjidemetriou, 1963; Omarov, 1964)

 F=−12˙M(t)M(t)˙r. (68)

Perturbation has the form of a linear drag acceleration acting always in the opposite to the motion of the body direction. This linear drag is proportional to the fractional mass-loss rate and is responsible for the reduction of the instantaneous velocity magnitude described above and depicted in Figure 4. This reduction leads to the new osculating orbit which will continue to evolve in time as long as the total mass is continuously changing in time or will be the actual physical orbit of the body if the perturbation ceases instantaneously.

Using the perturbing force , the equation of motion becomes

 ¨r=−GMr3r−12˙M(t)M(t)˙r. (69)

For the perturbation defined by equation , making use of the relations in equations and , one derives the following phase-dependent time-evolution equations for the osculating semi-major axis and eccentricity in the case of isotropic wind mass-loss (Hadjidemetriou, 1963)

 (dadt)iso =−(1+e2+2ecosν)1−e2˙MMa (70) (dedt)iso =−(e+cosν)˙MM. (71)

Equations and show that the osculating semi-major axis and eccentricity are in general phase-dependent. Using the identity we can re-write the nominator in equation as , which is always positive. This proves that for any isotropic mass-loss the semi-major axis can never decrease throughout the orbit and is not oscillating. Following the same procedure and given that the periastron position is given by one can prove that . Since and always for any isotropic mass-loss the periastron position can never decrease throughout the orbit and is not oscillating as well. However, the eccentricity does not follow a similar behaviour and can both increase and decrease throughout the orbit. The eccentricity undergoes oscillations on the orbital timescale (see e.g., Veras et al., 2011). The amplitude of these oscillations is proportional to , which is a scaled ratio of the orbital period to the mass-loss timescale. These oscillations can be smoothed out through orbit-averaging procedures when deriving the secular time-evolution equations of the osculating orbital elements. These secular time-evolution equations can be even more simplified in the adiabatic regime [see subsection 2.1.2 and Paper II for more details]. Over time the amplitude of these oscillations is increasing, since eventually the orbital period increases and becomes comparable to the mass-loss time scale. This is the point where adiabaticity is broken and the “adiabatic” equations no longer apply.

### 4.2. Perturbing Lagrangian treatment

Equations and can alternatively be derived from equations and for a zero-gauge by applying the perturbing Lagrangian

 ΔL=f(M1,M2)(12˙r2+GMr) (72)

where , with being a constant equal to the the total mass of the system at .

We note here that the Lagrangian was derived in such a way that it satisfies equation with the force given by equation . Since equation cannot uniquely determine the relative Lagrangian for a given force, the Lagrangian is not unique, has no physical meaning and it acts mainly as a mathematical tool to derive the time-evolution of the osculating or non-osculating orbital elements through equations and , also allowing for the gauge freedom. This freedom in the choice of the Lagrangian is depicted in the freedom in the functional dependence of the function on and . This means that the osculating orbital elements which have a specific physical interpretation cannot depend on the function and one can prove that equations and are constructed in such a way that for a zero-gauge that leads to the osculating orbital elements the functional dependence on diminishes. However, for a non-zero gauge and thus non-osculating orbital elements with no physical meaning this is no longer true. For the generalized Lagrange constraint given in equation and isotropic wind mass-loss the phase-dependent time-evolution equations for the contact orbital elements [see subsection 2.1.2] are given by

 (dacondt)iso =f(f+2)f+1[−2naesinν(1+ecosν)2(1−e2)5/2] (73) (decondt)iso =f(f+2)f+1[−nsinν(1+ecosν)2(1−e2)3/2]. (74)

We note here that equations and depend on the arbitrary function since they describe the time-evolution of the so-called contact elements which as we mentioned before do not have physical meaning. In the case of the non-zero gauge both the contact semi-major axis and the eccentricity oscillate throughout the orbit. In Paper II we compare the orbit-averaged time-evolution equations for the osculating and contact orbital elements assuming either adiabaticity or delta-function mass-loss at periastron.

## 5. Anisotropic wind mass-loss

Stellar mass-loss is time-dependent and depends in general on the latitude and longitude along the surface of the star. For example, there is much observational evidence that the equatorial regions around rapidly rotating hot-stars generally have an enhanced wind density, perhaps even (e.g., as in Be stars) a circumstellar disk (e.g., Owocki et al., 1998). More recent studies assume that the Be-star wind has a disk component in the equatorial plane plus a weak spherical wind above the poles (e.g., Bogomazov, 2005).

The formation of bipolar jets is strongly linked with the accretion onto protostars. Since the accretion process can persist for a long time after the star is born, jets usually accompany a disc (e.g., Blandford & Payne, 1982)

In this section, we do not study the case of a circumstellar disk. However, we assume that the mass-flux rate as well as the wind velocity is different at each point along the spherical surface of the star. We refer here to this case as anisotropic wind mass-loss. The results of this section are also applicable to the case of anisotropic planet evaporation (including jets) around the host star in a planet-star binary system (e.g., Boué et al., 2012; Veras et al., 2013).

Veras et al. (2013) were the first to study the case of anisotropic wind mass-loss from the star in a star-planet binary system. However, as we mentioned in Section 3, in this work there are some inconsistencies arising from a mistreated transformation from one reference system to the other. Both the importance of the choice of reference frame and the importance of the use of the right set of equations for the reference system chosen were pointed out in Section 3 as well. For the case of the anisotropic wind mass-loss we describe here, we discuss in Paper II the effect of the corrected time-evolution equations we derived in Section 3 [see equations and ] on the secular time-evolution equations of the orbital elements.

If the mass-losing body of mass has a mass-loss flux per solid angle , then we can express the mass-loss rate as

 ˙M1=−14π∫π0∫2π0J(ϕ,θ,t)sinθdϕdθ (75)

where and are the longitude and co-latitude respectively on the surface of the mass-losing body (donor).

If the mass at a specific point is lost with a velocity relative to the center-of-mass and along the radius of the donor, then assuming the spin axis of the donor to be fixed with respect to the orbit, the perturbing acceleration arising from the anisotropic wind mass-loss is given by

 Aaniso=−FanisoM1(t) (76)

where is the perturbing force arising from the anisotropic wind mass-loss. In the inertial reference system mentioned in Section 3 and depicted in Figure 3, the components of this perturbing force are given by

 Faniso,x =14π∫π0∫2π0J(ϕ,θ,t)u(ϕ,θ,t)dϕsinθcosϕdθ (77) Faniso,y =14π∫π0∫2π0J(ϕ,θ,t)u(ϕ,θ,t)dϕsinθsinϕdθ (78) Faniso,z =14π∫π0∫2π0J(ϕ,θ,t)u(ϕ,θ,t)dϕcosθdθ. (79)

It is important to mention here that both the velocity and the mass-flux rate do not depend on the relative position between the two bodies in the binary system. Thus, based on equations one sees that the perturbing force is not phase-dependent. However, both the mass-loss rate and the mass of the donor and thus both the perturbing force and acceleration remain time-dependent. In a similar way to the isotropic wind mass-loss case [see Section 4] where the mass-loss rate and the total mass was a function of time this time dependence induces implicitly a perturbation even though this perturbation is not explicitly phase-dependent. In Paper II we derive the secular time-evolution equations for the osculating orbital elements in the case of anisotropic wind mass-loss and for the phase-independent perturbing acceleration described by equations .

Picking a different reference system can lead to equations that are more easy to qualitatively interpret. Given the fact that the wind structure depends on the rotation of the star as well as its magnetic field we usually describe the form of the wind in the inertial reference frame . However, given the form of the wind in this system we can calculate the inertial perturbing force components and we can then transform for considerable simplifications to the reference system using equation or to reference system using equation . In this new systems the perturbing force components are given respectively by

 ⎛⎜⎝Faniso,pFaniso,qFaniso,l⎞⎟⎠=⎛⎜⎝Q11Q21Q31Q12Q22Q32Q13Q23Q33⎞⎟⎠⎛⎜⎝Faniso,xFaniso,yFaniso,z⎞⎟⎠ (80)

and

 ⎛⎜⎝Faniso,rFaniso,τFaniso,n⎞⎟⎠=⎛⎜⎝cosνsinν0−sinνcosν0001⎞⎟⎠⎛⎜⎝Faniso,pFaniso,qFaniso,l⎞⎟⎠. (81)

In Paper II we examine qualitatively some simple structures of anisotropic wind mass-loss including jets. We initially begin by describing the wind structure in the inertial reference frame and then by transforming to different reference systems that simplify our equations (reference systems and ) we derive the secular time-evolution equations of the orbital elements in these systems in the adiabatic regime.

## 6. Non-isotropic ejection/accretion in mass-transfer (reaction-forces)

Wind mass-loss is only one of the ways two bodies in a binary system can exchange mass. Another possibility is mass ejection or accretion from specific points and with specific velocities. Due to linear momentum conservation, this ejection and accretion induces kicks/reaction-forces to the mass-losing and the mass-accreting star respectively. To distinguish this mass-transfer process from the case of anisotropic wind mass-loss, we choose to refer to this process as non-isotropic ejection/accretion.

The most characteristic example of anisotropic ejection/accretion is Roche-Lobe-Overflow (RLOF). In this case, the mass is ejected from the Lagrangian point . This mass can escape from the system be re-accreted from the donor (self-accretion) hit the radius of the accretor (direct impact) hit the radius of an accretion disk that has been formed around the companion (e.g., Sepinsky et al., 2007, 2009, 2010). We refer to the case when all the mass ejected is accreted as conservative mass-transfer while to the case when some mass is lost from the system as non-conservative.

In the next two sections, we study the cases of non-isotropic ejection and accretion in both conservative and non-conservative manner, keeping the ejection/accretion points and velocities as general parameters of the problem.

### 6.1. Conservative case (˙M=˙Jorb=0)

We assume an eccentric binary system consisting of two masses and . We refer to these bodies as the donor and the accretor, respectively.

Binary component loses mass at a rate , with absolute velocity and from point relative to the center-of-mass of the body. Binary component accretes mass at a rate , absolute velocity and to point relative to the center-of-mass of the body. The mass-transfer stream is assumed to remain always confined to the orbital plane.

Following the work of Hadjidemetriou (1963) and Sepinsky et al. (2007), the equation of motion of the two-body system can be written in an inertial frame as follows

 ¨r =−G(M1+M2)r3r+˙M2M2(w2+ωorb×rA2) −˙M1M1(w1+ωorb×rA1)−˙M2M2v2+˙M1M1v1 (82) +¨M2M2rA2−¨M1M1rA1 =−G(M1+M2)r3r+˙M2M2(w2−v2+ωorb×rA2) −˙M1M1(w1−v1+ωorb×rA1) +¨M2M2rA2−¨M1M1rA1 (83)

where and are the orbital velocities of the two bodies with respect to the inertial frame and is the orbital angular frequency. Equation takes into account the motion of the center-of-mass of the bodies and ignores the forces on the bodies from the mass-transfer stream under the assumption that the stream’s mass is negligible. We note here that when writing the linear momentum conservation equations, we assume the absolute velocities and to be in the opposite direction to the direction of motion of the bodies.

From equation we can see that anisotropic accretion and ejection introduces the reaction forces

 ˙M2M2(w2−v2+ωorb×rA2) (84) − ˙M1M1(w1−v1+ωorb×rA1), (85)

for accretion and ejection respectively. The terms with the second order time derivative in equation are related to the acceleration of the centers-of-mass of the two stars relative to their positions in the unperturbed no mass-loss case [see Section 2 in Sepinsky et al. (2007) for a detailed description of the relevant terms].

The reaction forces and emerge from applying the linear momentum conservation law in the accretor and the donor respectively. As a reminder, we mention here that in the isotropic mass ejection/accretion case we have

 ∑i{(w1−v1)i+ωorb×rA1i}=0 (86) ∑i{(w2−v2)i+ωorb×rA2i}=0. (87)

These terms are zero due to isotropy.

If all the mass lost from the donor is accreted by the companion, and any orbital angular momentum transported by the transferred mass is immediately returned to the orbit, then as we mentioned before the mass-transfer is conservative. In this case, both the total system mass and orbital angular momentum are conserved and we can then write . Equations and in the case of conservative mass-transfer yield

 ¨r =−G(M1+M2)r3r