The evolution of hierarchical triple star-systems

The evolution of hierarchical triple star-systems

[    [    [ \orgnameLeiden Observatory, Leiden University, \postcodePO Box 9513 \cityLeiden, \cnyThe Netherlands

Field stars are frequently formed in pairs, and many of these binaries are part of triples or even higher-order systems. Even though, the principles of single stellar evolution and binary evolution, have been accepted for a long time, the long-term evolution of stellar triples is poorly understood. The presence of a third star in an orbit around a binary system can significantly alter the evolution of those stars and the binary system. The rich dynamical behaviour in three-body systems can give rise to Lidov-Kozai cycles, in which the eccentricity of the inner orbit and the inclination between the inner and outer orbit vary periodically. In turn, this can lead to an enhancement of tidal effects (tidal friction), gravitational-wave emission and stellar interactions such as mass transfer and collisions. The lack of a self-consistent treatment of triple evolution, including both three-body dynamics as well as stellar evolution, hinders the systematic study and general understanding of the long-term evolution of triple systems. In this paper, we aim to address some of these hiatus, by discussing the dominant physical processes of hierarchical triple evolution, and presenting heuristic recipes for these processes. To improve our understanding on hierarchical stellar triples, these descriptions are implemented in a public source code TrES  which combines three-body dynamics (based on the secular approach) with stellar evolution and their mutual influences. Note that modelling through a phase of stable mass transfer in an eccentric orbit is currently not implemented in TrES , but can be implemented with the appropriate methodology at a later stage.



addressref=aff1, corref=aff1, ]\initsS\fnmS. \snmToonen addressref=aff1, ]\initsA\fnmA. \snmHamers addressref=aff1, ]\initsS\fnmS. \snmPortegies Zwart


binaries (including multiple): close \kwdstars: evolution

1 Introduction

The majority of stars are members of multiple systems. These include binaries, triples, and higher order hierarchies. The evolution of single stars and binaries have been studied extensively and there is general consensus over the dominant physical processes (Postnov and Yungelson, 2014; Toonen et al., 2014). Many exotic systems, however, cannot easily be explained by binary evolution, and these have often been attributed to the evolution of triples, for example low-mass X-ray binaries (Eggleton and Verbunt, 1986) and blue stragglers (Perets and Fabrycky, 2009). Our lack of a clear understanding of triple evolution hinders the systematic exploration of these curious objects. At the same time triples are fairly common; Our nearest neighbour Cen is a triple star system (Tokovinin, 2014a), but more importantly % of the low-mass stars are in triples (Tokovinin, 2008; Raghavan et al., 2010; Tokovinin, 2014b; Moe and Di Stefano, 2016) a fraction that gradually increases (Duchêne and Kraus, 2013) to % for spectral type B stars (Remage Evans, 2011; Sana et al., 2014; Moe and Di Stefano, 2016).

The theoretical studies of triples can classically be divided into three-body dynamics and stellar evolution, which both are often discussed separately. Three-body dynamics is generally governed by the gravitational orbital evolution, whereas the stellar evolution is governed by the internal nuclear burning processes in the individual stars and their mutual influence.

Typical examples of studies that focused on the three-body dynamics include Ford et al. (2000); Fabrycky and Tremaine (2007); Naoz et al. (2013); Naoz and Fabrycky (2014); Liu et al. (2015a), and stellar evolution studies include Eggleton and Kiseleva (1996); Iben and Tutukov (1999); Kuranov et al. (2001). Interdisciplinary studies, in which the mutual interaction between the dynamics and stellar aspects are taken into account are rare (Kratter and Perets, 2012; Perets and Kratter, 2012; Hamers et al., 2013; Shappee and Thompson, 2013; Michaely and Perets, 2014; Naoz et al., 2016), but demonstrate the richness of the interacting regime. The lack of a self consistent treatment hinders a systematic study of triple systems. This makes it hard to judge the importance of this interacting regime, or how many curious evolutionary products can be attributed to triple evolution. Here we discuss triple evolution in a broader context in order to address some of these hiatus.

In this paper we discuss the principle complexities of triple evolution in a broader context (Sect. 2). We start by presenting an overview of the evolution of single stars and binaries, and how to extend these to triple evolution. In the second part of this paper we present heuristic recipes for simulating their evolution (Sect. 3). These recipes combine three-body dynamics with stellar evolution and their mutual influences, such as tidal interactions and mass transfer. These descriptions are summarized in a public source code TrES  with which triple evolution can be studied.

2 Background

We will give a brief overview of isolated binary evolution (Sect. 2.2) and isolated triple evolution (Sect. 2.3). We discuss in particular under what circumstances triple evolution differs from binary evolution and what the consequences are of these differences. We start with a brief summary of single star evolution with a focus on those aspects that are relevant for binary and triple evolution.

Parameters Stellar Orbital
Single star -
Binary , ,
Triple , , , , , ,
, , ,
Table 1: Necessary parameters to describe a single star system, a binary and a triple. For stellar parameters, age and metallicity of each star can be added. The table shows that as the multiplicity of a stellar system increases from one to three, the problem becomes significantly more complicated.

2.1 Single stellar evolution

Hydrostatic and thermal equilibrium in a star give rise to temperatures and pressures that allow for nuclear burning, and consequently the emission of the starlight that we observe. Cycles of nuclear burning and exhaustion of fuel regulate the evolution of a star, and sets the various phases during the stellar lifetime.

The evolution of a star is predominantly determined by a single parameter, namely the stellar mass (Tbl. 1). It depends only slightly on the initial chemical composition or the amount of core overshooting111 Overshooting refers to a chemically mixed region beyond the boundary of the convective core (e.g. Stothers, 1963; Massevitch et al., 1979; Maeder and Meynet, 1991), as predicted by basic stellar evolutionary theory, i.e. the Schwarzschild criterion. A possible mechanism is convection carrying material beyond the boundary due to residual velocity. For the effects of overshooting on stellar evolution, see e.g. Bressan et al. (2015)..

2.1.1 Timescales

Fundamental timescales of stellar evolution are the dynamical (), thermal (), and nuclear timescale (). The dynamical timescale is the characteristic time that it would take for a star to collapse under its own gravitational attraction without the presence of internal pressure:


where and are the radius and mass of the star. It is a measure of the timescale on which a star would expand or contract if the hydrostatic equilibrium of the star is disturbed. This can happen for example because of sudden mass loss.

A related timescale is the time required for the Sun to radiate all its thermal energy content at its current luminosity:


where is the luminosity of the star. In other words, when the thermal equilibrium of a star is disturbed, the star will move to a new equilibrium on a thermal (or Kelvin–Helmholtz) timescale

Finally, the nuclear timescale represents the time required for the star to exhaust its supply of nuclear fuel at its current luminosity:


where is the efficiency of nuclear energy production, is the speed of light, and is the amount of mass available as fuel. For core hydrogen burning, and . Assuming a mass-luminosity relation of , with empirically (e.g. Salaris and Cassisi, 2005; Eker et al., 2015), it follows that massive stars live shorter and evolve faster than low-mass stars.

For the Sun, 30 min, 30 Myr, and 10 Gyr. Typically, , which allows us to quantitatively predict the structure and evolution of stars in broad terms.

2.1.2 Hertzsprung-Russell diagram

Figure 1: \csentenceHertzsprung-Russell diagram Evolutionary tracks for seven stars in the HR-diagram with masses 1, 1.5, 2.5, 4, 6.5, 10, and 15 at solar metallicity. Specific moments in the evolution of the stars are noted by blue circles as explained in the text. The tracks are calculated with SeBa (Portegies Zwart and Verbunt, 1996; Toonen et al., 2012). The dashed lines show lines of constant radii by means of the Stefan–Boltzmann law.

The Hertzsprung-Russell (HR) diagram in Fig.1 shows seven evolutionary tracks for stars of different masses. The longest phase of stellar evolution is known as the main-sequence (MS), in which nuclear burning takes place of hydrogen in the stellar core. The MS occupies the region in the HR-diagram between the stellar birth on the zero-age MS (ZAMS, blue circles in Fig. 1) and the end of the MS-phase (terminal-age MS (TAMS), blue circles in Fig. 1). Stars more massive than 1.2 contract slightly at the end of the MS when the stellar core runs out of hydrogen. This can be seen in Fig. 1 as the hook in the tracks leading up to the TAMS.

After the TAMS, hydrogen ignites in a shell around the core. Subsequently the outer layers of the star expand rapidly. This expansion at roughly constant luminosity results in a lower effective temperature and a shift to the right in the HR-diagram. Stars of less than 13 reach effective temperatures as low as (K) 5000K before helium ignition. At this point (denoted by a blue circle in Fig. 1) they start to ascend the red giant branch (RGB) which goes hand in hand with a strong increase in luminosity and radius. On the right of the RGB in the HR-diagram, lies the forbidden region where hydrostatic equilibrium cannot be achieved. Any star in this region will rapidly move towards the RGB. The red giant star consists of a dense core and an extended envelope up to hundreds of solar radii. When the temperature in the core reaches K, helium core burning commences and the red giant phase has come to an end. For stars less massive than 2, helium ignites degenerately in a helium flash. For stars more massive than 13, helium ignites before their effective temperature has decreased to a few thousand Kelvin; the shift to the right in the HR-diagram is truncated when helium ignites.

During helium burning the stellar tracks make a loop in the HR-diagram, also known as the horizontal branch. This branch is marked in Fig. 1 by a blue circle at its maximum effective temperature. The loop goes hand in hand with a decrease and increase of the stellar radius. As the burning front moves from the core to a shell surrounding the core, the outer layers of the star expand again and the evolutionary track bends back to right in the HR-diagram.

As the core of the star reaches temperatures of K, carbon ignites in the star (denoted by a blue circle in Fig. 1).

As the core of the star becomes depleted of helium, helium burning continues in a shell surrounding the inert carbon-oxygen core. The star has now has reached the supergiant-phases of its life. The star ascents the asymptotic giant branch (AGB) reaching its maximum size of about a thousand solar radii.

Figure 2: \csentenceEvolution of stellar radius Radius as a function of stellar age for two stars with masses 4 and 6.5 at solar metallicity. Specific moments in the evolution of the stars are noted by blue circles as for Fig. 1. The radius evolution is calculated with SeBa (Portegies Zwart and Verbunt, 1996; Toonen et al., 2012). The figure also shows that high-mass stars evolve faster and live shorter than lower-mass stars.

Fig. 2 shows the variation of the outer radius as the star evolves in its lifetime. It illustrates the dramatic increases in radius during the RGB- and AGB-phases as previously discussed. Shrinkage of star occur after helium ignition, and to a lesser degree at the end of the MS. The radial evolution is of particular interest for binaries and triples, as a star is more likely to initiate mass transfer (i.e. fill its Roche lobe) when its envelope is extended e.g. on the RGB or AGB.

2.1.3 Stellar winds

During the lifetime of a star, a major fraction of the star’s mass is lost by means of stellar winds (Lamers and Cassinelli, 1999; Owocki, 2013). The winds deposit enriched material back into the ISM and can even collide with previously ejected matter to form stellar-wind bubbles and planetary nebulae.

Stellar winds develop for almost all stars, but the mass losses increases dramatically for more evolved stars and for more massive stars. The winds of AGB stars (for a review Höfner, 2015) are characterized by extremely high mass-loss rates ( ) and low terminal velocities (5-30). For stars up to 8, these ’superwinds’ remove the entire stellar envelope. AGB-winds are driven by radiation pressure onto molecules and dust grains in the cold outer atmosphere. The winds are further enhanced by the stellar pulsations that increase the gas density in the extended stellar atmosphere where the dust grains form.

For massive O and B-type stars, strong winds already occur on the MS. These winds (e.g. Puls et al., 2008; Vink, 2015) are driven by another mechanism, i.e. radiation pressure in the continuum and absorption lines of heavy elements. The winds are characterized with high mass-loss rates () and high velocities (several 100-1000) (e.g. Kudritzki and Puls, 2000). For stars of more than 30, the mass-loss rate is sufficiently large that the evolution of the star is significantly affected, as the timescale for mass loss is smaller than the nuclear timescale. In turn the uncertainties in our knowledge of the stellar wind mechanisms, introduces considerable uncertainties in the evolution of massive stars.

2.1.4 Stellar remnants

The evolution of a star of less than 6.5 comes to an end as helium burning halts at the end of the AGB. Strong winds strip the core of the remaining envelope and this material forms a planetary nebula surrounding the core. The core cools and contracts to form a white dwarf (WD) consisting of carbon and oxygen (CO).

Slightly more massive stars up to 11 experience an additional nuclear burning phase. Carbon burning leads to the formation of a degenerate oxygen-neon (ONe) core. Stars up to 8 follow a similar evolutionary path discussed above, but they end their lives as oxygen-neon white dwarfs. In the mass range 8-11, the oxygen-neon core reaches the Chandrasekhar mass, and collapses to a neutron star (NS).

More massive stars than 11 go through a rapid succession of nuclear burning stages and subsequent fuel exhaustion. The nuclear burning stages are sufficiently short, that the stellar envelope hardly has time to adjust to the hydrodynamical and thermal changes in the core. The position of the star in the HR-diagram remains roughly unchanged. The stellar evolution continues until a iron core is formed after which nuclear burning cannot release further energy. The star then collapses to form a NS or a black hole (BH). An overview of the initial mass ranges and the corresponding remnants are given in Tbl. 2.

Initial mass () Remnant type Remnant mass ()
1–6.5 CO WD 0.5-1.1
6.5–8 ONe WD 1.1–1.44
8–23 NS 1.1–2
23 BH 5
Table 2: Initial stellar mass range and the corresponding remnant type and mass. Note that the given masses represent approximate values. They are dependent on the metallicity and on stellar evolutionary processes that are not understood well, such as stellar winds and core overshooting (footnote 1).

When a star is part of a compact stellar system, its evolution can be terminated prematurely when the star looses its envelope in a mass-transfer phase. After the envelope is lost, the star may form a remnant directly. If on the other hand, the conditions to sustain nuclear burning are fulfilled, the star can evolve further as a hydrogen-poor helium rich star i.e. helium MS star or helium giant star.

Due to the mass loss, the initial mass ranges given in Tbl. 2 can be somewhat larger. Furthermore, if a star with a helium core of less than 0.32(e.g. Han et al., 2002) looses its envelope as a result of mass transfer before helium ignition, the core contracts to form a white dwarf made of helium instead of CO or ONe.

2.1.5 Supernova explosions

When a high-mass star reaches the end of its life and its core collapses to a NS or BH, the outer layers of the star explode in a core-collapse supernova (SN) event. The matter that is blown off the newly formed remnant, enriches the ISM with heavy elements. Any asymmetry in the SN, such as in the mass or neutrino loss (e.g. Lai, 2004; Janka, 2012), can give rise to a natal-kick to the star. Neutron stars are expected to receive a kick at birth of about 400 (e.g. Cordes et al., 1993; Lyne and Lorimer, 1994; Hobbs et al., 2005), however smaller kick velocities in the range of have been deduced for neutron stars in high-mass X-ray binaries (Pfahl et al., 2002). Also, whether or not black holes that are formed in core-collapse supernova receive a kick is still under debate (e.g. Gualandris et al., 2005; Repetto et al., 2012; Wong et al., 2014; Repetto and Nelemans, 2015; Zuo, 2015).

2.2 Binary evolution

The evolution of a binary can be described by the masses of the stars and , the semi-major axis , and the eccentricity . A useful picture for binaries is the Roche model, which describes the effective gravitational potential of the binary. It is generally based on three assumptions: 1) the binary orbit is circular 2) the rotation of the stellar components are synchronized with the orbit 3) the stellar components are small compared to the distance between them. The first two assumptions are expected to hold for binaries that are close to mass transfer because of tidal forces (Sect. 2.2.3). Under the three assumptions given above, the stars are static in a corotating frame of reference. The equipotential surface around a star in which material is gravitationally bound to that star is called the Roche lobe. The Roche radius is defined as the radius of a sphere with the same volume of the nearly spherical Roche lobe, and is often approximated (Eggleton, 1983) by:


where the mass ratio . If one of the stars in the binary overflows its Roche lobe, matter from the outer layers of the star can freely move through the first Lagrangian point L1 to the companion star. Binaries with initial periods less than several years (depending on the stellar masses) will experience at least one phase of mass transfer, if the stars have enough time to evolve.

If the stars do not get close to Roche lobe overflow (RLOF), the stars in a binary evolve effectively as single stars, slowly decreasing in mass and increasing in radius and luminosity until the remnant stage. The binary orbit can be affected by stellar winds, tides and angular momentum losses such as gravitational wave emission and magnetic braking. These processes are discussed in the following three sections. In the last three sections of this chapter we describe how RLOF affects a binary.

2.2.1 Stellar winds in binaries

Wind mass loss affects a binary orbit through mass and angular momentum loss. Often the assumption is made that the wind is spherically symmetric and fast with respect to the orbit. In this approximation, the wind does not interact with the binary orbit directly, such that the process is adiabatic. Furthermore, the orbital eccentricity remains constant (Huang, 1956, 1963).

If none of the wind-matter is accreted, the wind causes the orbit to widen. From angular momentum conservation, it follows as:


where and are the masses of the stars, is the mass lost in the wind of the star with mass (), is the semi-major axis of the orbit, and the change in the orbital separation with no wind accretion. Eq. 5 can be rewritten to:


where and are the semi-major axis of the orbit before and after the wind mass loss, and is the amount of matter lost in the wind ().

While the two stars in the binary are in orbit around each other, the stars can accrete some of the wind material of the other star. Including wind accretion, the orbit changes as:


where the star with mass accretes at a rate of . Note that Eq. 7 reduces to Eq. 5 for complete non-conservative mass transfer i.e. . Wind accretion is often modelled by Bondi-Hoyle accretion (Bondi and Hoyle, 1944; Livio and Warner, 1984). This model considers a spherical accretion onto a point mass that moves through a uniform medium. Wind accretion is an important process known to operate in high-mass X-ray binaries (Tauris and van den Heuvel, 2006; Chaty, 2011) and symbiotic stars (Mikolajewska, 2002; Sokoloski, 2003).

The assumptions of a fast and spherically symmetric wind are not always valid. The former is not strictly true for all binary stars i.e. an evolved AGB-star has a wind of 5-30 (e.g. Höfner, 2015), which is comparable to the velocity of stars in a binary of . Hydrodynamical simulations of such binaries suggest that the wind of the donor star is gravitationally confined to the Roche lobe of the donor star (Mohamed and Podsiadlowski, 2007; de Val-Borro et al., 2009; Mohamed and Podsiadlowski, 2011). The wind can be focused towards the orbital plane and in particular towards the companion star. This scenario (often called wind Roche-lobe overflow (wRLOF) or gravitational focusing) allows for an accretion efficiency of up to 50%, which is significantly higher than for Bondi-Hoyle accretion. A requirement for wRLOF to work is that the Roche lobe of the donor star is comparable or smaller than the radius where the wind is accelerated beyond the escape velocity. wRLOF is supported by observations of detached binaries with very efficient mass transfer (Karovska et al., 2005; Blind et al., 2011)

Furthermore, the assumption of adiabatic mass loss is inconsistent with binaries in which the orbital timescale is longer than the mass-loss timescale. The effects of instantaneous mass loss has been studied in the context of SN explosions, and can even lead to the disruption of the binary system (see also Sect. 2.2.7). However, also wind mass-loss can have a non-adiabatic effect on the binary orbit (e.g. Hadjidemetriou, 1966; Rahoma et al., 2009; Veras et al., 2011) if the mass-loss rate is high and the orbit is wide. Under the assumption that mass-loss proceeds isotropically, the wind causes the orbit to widen, as in the case for adiabatic mass loss. However, the eccentricity may decrease or increase, and may even lead to the disruption of the binary system (see e.g. Veras et al., 2011, for a detailed analysis of the effects of winds on sub-stellar binaries in which an exoplanet orbits a host star). Toonen, Hollands, Gaensicke, & Boekholt, in prep. show that also (intermediate-mass) stellar binaries can be disrupted during the AGB-phases when the mass loss rates are high () for orbital separations approximately larger then (yr where is the orbital period).

Lastly, anisotropic mass-loss might occur in fast-rotating stars or systems that harbour bipolar outflows. Rotation modifies the structure and evolution of a star, and as such the surface properties of the star where the wind originates (see Maeder and Meynet, 2012, for a review). For an increasing rate of rotation until critical rotation, the stellar winds increasingly depart from a spherical distribution (see e.g. Georgy et al., 2011). Additionally, the bipolar outflows or jets are associated with protostars, evolved post-AGB stars and binaries containing compact objects. Their origin is most likely linked to the central object or the accretion disk (e.g. O’Brien, 1991).

The effect of anisotropic mass loss on the orbit of a binary system is important primarily for wide binaries (e.g. Parriott and Alcock, 1998; Veras et al., 2013). Specifically, Veras et al. (2013) show that the relative contribution of the anisotropic terms to the overall motion scale as . If the mass loss is symmetric about the stellar equator, the mass loss does not affect the orbital motion in another way than for the isotropic case. Veras et al. (2013) conclude that the isotropic mass-loss approximation can be used safely to model the orbital evolution of a planet around a host star until orbital separations of hundreds of AU. For a fixed total mass of the system, the effects of anisotropic mass loss are further diminished with decreasing mass ratio (i.e. for systems with more equal masses), such that the assumption of isotropic mass-loss is robust until even larger orbital separations for stellar binaries.

2.2.2 Angular momentum losses

Angular momentum loss from gravitational waves (GW) and magnetic braking act to shrink the binary orbit (e.g. Peters, 1964; Verbunt and Zwaan, 1981). Ultimately this can lead to RLOF of one or both components and drive mass transfer.

The strength of GW emission depends strongly on the semi-major axis, and to lesser degree on the eccentricity. It affects the orbits as:




where and are the change in orbital separation and eccentricity averaged over a full orbit (Peters, 1964). Accordingly, GW emission affects most strongly the compact binaries. These binaries are a very interesting and the only known source of GWs for GW interferometers such as LIGO, VIRGO and eLISA.

Magnetic braking extracts angular momentum from a rotating magnetic star by means of an ionized stellar wind (Schatzman, 1962; Huang, 1966; Skumanich, 1972). Even when little mass is lost from the star, the wind matter can exert a significant spin-down torque on the star. This happens when the wind matter is forced to co-rotate with the magnetic field. If the star is in a compact binary and forced to co-rotate with the orbit due to tidal forces, angular momentum is essentially removed from the binary orbit as well (Verbunt and Zwaan, 1981). This drain of angular momentum results in a contraction of the orbit.

Magnetic braking plays an important role in the orbital evolution of interacting binaries with low-mass donor stars, such as cataclysmic variables and low-mass X-ray binaries (Knigge et al., 2011; Tauris and van den Heuvel, 2006). For magnetic braking to take place, the donor star is expected to have a mass between 0.2-1.2, such that the star has a radiative core and convective envelope to sustain the magnetic field. The strength of magnetic braking is still under debate and several prescriptions exist (see Knigge et al., 2011, for a review).

2.2.3 Tides

The presence of a companion star introduces tidal forces in the binary system that act on the surface of the star and lead to tidal deformation of the star. If the stellar rotation is not synchronized or aligned with the binary orbit, the tidal bulges are misaligned with the line connecting the centres of mass of the two stars. This produces a tidal torque that allows for the transfer of angular momentum between the stars and the orbit. Additionally, energy is dissipated in the tides, which drains energy from the orbit and rotation. Tidal interaction drives the binary to a configuration of lowest energy e.g. it strives to circularize the orbit, synchronize the rotation of the stars with the orbital period and align the stellar spin with respect to the orbital spin. See Zahn (2008) and Zahn (2013) for recent reviews.

For binaries with extreme mass ratios, a stable solution does not exist (Darwin, 1879; Hut, 1980). In this scenario a star is unable to extract sufficient angular momentum from the orbit to remain in synchronized rotation. Tidal forces will cause the orbit to decay and the companion to spiral into the envelope of the donor star. This tidal instability occurs when the angular momentum of the star , with the orbital angular momentum and , where is the moment of inertia and the spin angular frequency.

Hut (1981) derives a general qualitative picture of tidal evolution and its effect on the orbital evolution of a binary system. Hut (1981) considers a model in which the tides assume their equilibrium shape, and with very small deviations in position and amplitude with respect to the equipotential surfaces of the stars. If a companion star with mass raises tides on a star with mass , the change of binary parameters due to tidal friction is:


where , and is the mean orbital angular velocity. The star with mass has an apsidal motion constant , gyration radius , and spin angular frequency . represents the typical timescale on which significant changes in the orbit take place through tidal evolution. The parameters are polynomial expressions given by (Hut, 1981):


The degree of tidal interaction strongly increases with the ratio of the stellar radii to the semi-major axis of the orbit (Eq. 10 11 and 12). Therefore, tidal interaction mostly affect the orbits of relatively close binaries, unless the eccentricities are high and/or the stellar radii are large.

The tidal timescale ( Eq. 10-12) is subject to debate due to quantitative uncertainties in tidal dissipation mechanisms (Witte and Savonije, 1999; Willems, 2003; Meibom and Mathieu, 2005). Tidal dissipation causes the misalignment of the tidal bulges with the line connecting the centres of mass of the two stars. For stars (or planets) with an outer convection zone, the dissipation is often attributed222For alternative mechanisms, see (Goodman and Dickson, 1998; Tas00, 2000; Savonije and Witte, 2002). to turbulent friction in the convective regions of the star (Zahn, 1977, 1989; Goldman and Mazeh, 1991). For stars with an outer radiation zone, the dominant dissipation mechanism identified so far is radiative damping of stellar oscillations that are exited by the tidal field i.e. dynamical tides (Zahn, 1975, 1977). Despite the uncertainties in tidal dissipation mechanisms, it is generally assumed that circularization and synchronization is achieved before RLOF in a binary.

2.2.4 Mass transfer

Whether or not mass transfer is stable depends on the response of the donor star upon mass loss, and the reaction of the Roche lobe upon the re-arrangement of mass and angular momentum within the binary (e.g. Webbink, 1985; Hjellming and Webbink, 1987; Pols and Marinus, 1994; Soberman et al., 1997). If the donor star stays approximately within its Roche lobe, mass transfer is dynamically stable. When this is not the case, the donor star will overflow its Roche lobe even further as mass is removed. This leads to a runaway situation that progresses into a common-envelope (CE, Paczynski, 1976). During the CE-phase, the envelope of the donor star engulfs both stars causing them to spiral inwards until both stars merge or the CE is expelled.

Due to the mass loss, the donor star falls out of hydrostatic and thermal equilibrium. The radius of the star changes as the star settles to a new hydrostatic equilibrium, and subsequently thermal equilibrium. The stellar response upon mass loss depends critically on the structure of the stellar envelope i.e. the thermal gradient and entropy of the envelope. In response to mass loss, stars with a deep surface convective zone tend to expand, whereas stars with a radiative envelope tend to shrink rapidly. Therefore, giant donor stars with convective envelopes favour CE-evolution upon RLOF 333Unless , such that the orbit and the Roche lobe expand significantly upon mass transfer (e.g. Eqs. 18-19). As giants have radii of several hundreds to thousands of Solar radii, the orbit at the onset of mass transfer is of the same order of magnitude. On the other hand, donor stars on the MS with radiative envelope often lead to dynamically stable mass transfer in binaries with short orbital periods (e.g. Toonen et al., 2014).

2.2.5 Common-envelope evolution

During the CE-phase, the core of the donor star and the companion are contained within a CE. Friction between these objects and the slow-rotating envelope is expected to cause the objects to spiral-in. If this process does not release enough energy and angular momentum to drive off the entire envelope, the binary coalesces. On the other hand if a merger can be avoided, a close binary remains in which one or both stars have lost their envelopes. The evolution of such a star is significantly shortened, or even terminated prematurely if it directly evolves to a remnant star.

The systems that avoid a merger lose a significant amount of mass and angular momentum during the CE-phase. The orbital separation of these systems generally decreases by two orders of magnitude, which affects the further evolution of the binary drastically. The CE-phase plays an essential role in the formation of short-period systems with compact objects, such as X-ray binaries, and cataclysmic variables. In these systems the current orbital separation is much smaller than the size of the progenitor of the donor star, which had giant-like dimensions at the onset of the CE-phase.

Despite of the importance of the CE-phase and the enormous efforts of the community, the CE-phase is not understood in detail (see Ivanova et al., 2013, for a review). The CE-phase involves a complex mix of physical processes, such as energy dissipation, angular momentum transport, and tides, over a large range in time- and length-scales. A complete simulation of the CE-phase is still beyond our reach, but great progress has been made with hydrodynamical simulations in the last few years (Ricker and Taam, 2012; Passy et al., 2012b; Nandez et al., 2015). The uncertainty in the CE-phase is one of the aspects of the theory of binary evolution that affects our understanding of the evolutionary history of a specific binary population most (e.g. Toonen and Nelemans, 2013; Toonen et al., 2014).

The classical way to treat the orbital evolution due to the CE-phase, is the -formalism. This formalism considers the energy budget of the initial and final configuration (Tutukov and Yungelson, 1979);


where is the binding energy of the envelope, and are the orbital energy of the pre- and post-mass transfer binary. The -parameter describes the efficiency with which orbital energy is consumed to unbind the CE. When both stars have loosely bound envelopes, such as for giants, both envelopes can be lost simultaneously (hereafter double-CE, see Brown, 1995; Nelemans et al., 2001). In Eq. 14 is then replaced by the sum of the binding energy of each envelope to its host star:


The binding energy of the envelope of the donor star in Eq. 14 and 15 is given by:


where is the radius of the donor star, is the envelope mass of the donor and depends on the structure of the donor (de Kool et al., 1987; Dewi and Tauris, 2000; Xu and Li, 2010; Loveridge et al., 2011). The parameters and are often combined in one parameter .

According to the alternative -formalism (Nelemans et al., 2000), angular momentum is used to expel the envelope of the donor star, according to:


where and are the orbital angular momentum of the pre- and post-mass transfer binary respectively. The parameters and represent the mass of the donor and accretor star, respectively, and is the mass lost by the donor star. The -parameter describes the efficiency with which orbital angular momentum is used to blow away the CE.

Valuable constraints on CE-evolution have come from evolutionary reconstruction studies of observed samples of close binaries and from comparing those samples with the results of binary population synthesis studies. The emerging picture is that for binaries with low mass ratios, the CE-phase leads to a shrinkage of the orbit. For the formation of compact WD-MS binaries with low-mass MS companions, the orbit shrinks strongly (, see Zorotovic et al., 2010; Toonen and Nelemans, 2013; Portegies Zwart, 2013; Camacho et al., 2014). However, for the formation of the second WD in double WDs, the orbit only shrinks moderately (, see Nelemans et al., 2000, 2001; van der Sluys et al., 2006). When binaries with approximately equal masses come in contact, mass transfer leads to a modest widening of the orbit, alike the -formalism (Nelemans et al., 2000, 2001). The last result is based on a study of the first phase of mass transfer for double WDs, in which the first WD is formed. Woods et al. (2012) suggested that this mass transfer episode can occur stably and non-conservatively even with donor star (early) on the red giant branch. Further research is needed to see if this evolutionary path suffices to create a significant number of double WDs.

2.2.6 Stable mass transfer

Whereas the duration of the CE-phase is likely of the order of yr (i.e. the thermal timescale of the envelope), stable mass transfer occurs on much longer timescales. Several driving mechanisms exist for stable mass transfer with their own characteristic mass transfer timescales. The donor star can drive Roche lobe overflow due to its nuclear evolution or due to the thermal readjustment from the mass loss. Stable mass transfer can also be driven by the contraction of the Roche lobe due to angular momentum losses in the system caused by gravitational wave radiation or magnetic braking.

When mass transfer proceeds conservatively the change in the orbit is regulated by the masses of the stellar components. For circular orbits,


where the subscript  and  denote the pre- and post-mass transfer values. In general, the donor star will be the more massive component in the binary and the binary orbit will initially shrink in response to mass transfer. After the mass ratio is approximately reversed, the orbit widens. In comparison with the pre-mass transfer orbit, the post-mass transfer orbit is usually wider with a factor of a few (Toonen et al., 2014).

If the accretor star is not capable of accreting the matter conservatively, mass and angular momentum are lost from the system. The evolution of the system is then dictated by how much mass and angular momentum is carried away. Assuming angular momentum conservation and neglecting the stellar rotational angular momentum compared to the orbital angular momentum, the orbit evolves as (e.g. Massevitch and Yungelson, 1975; Pols and Marinus, 1994; Postnov and Yungelson, 2014):


where the accretor star captures a fraction of the transferred matter, and the matter that is lost carries specific angular momentum equal to a multiple of the specific orbital angular momentum of the binary:


Different modes of angular momentum loss exist which can lead to a relative expansion or contraction of the orbit compared to the case of conservative mass transfer (Soberman et al., 1997; Toonen et al., 2014). For example, the generic description of orbital evolution of Eq. 19 reduces to that of conservative mass transfer (Eq. 18) for or . Also, Eq. 19 reduces to Eq. 7 describing the effect of stellar winds on the binary orbit, under the assumption of specific angular momentum loss equal to that of the donor star ( or ). Depending on which mode of angular momentum loss is applicable, the further orbital evolution and stability of the system varies.

Stable mass transfer influences the stellar evolution of the donor star and possibly that of the companion star. The donor star is affected by the mass loss, which leads to a change in the radius on long timescales compared to a situation without mass loss (Hurley et al., 2000). Stable mass transfer tends to terminate when the donor star has lost most of its envelope, and contracts to form a remnant star or to a hydrogen-poor helium rich star. In the latter case the evolution of the donor star is significantly shortened, and in the former it is stopped prematurely, similar to what was discussed previously for the CE-phase.

If the companion star accretes a fraction or all of the transferred mass, evolution of this star is affected as well. Firstly, if due to accretion, the core grows and fresh fuel from the outer layers is mixed into the nuclear-burning zone, the star is ’rejuvenated’ (see e.g. Vanbeveren and De Loore, 1994). These stars can appear significantly younger than their co-eval neighbouring stars in a cluster444Stable mass transfer is one of the proposed evolutionary pathways for the formation of blue stragglers (for a review Davies, 2015). These are MS stars in open and globular clusters that are more luminous and bluer than other MS stars in the same cluster.. Secondly, the accretor star adjusts its structure to a new equilibrium. If the timescale of the mass transfer is shorter than the thermal timescale of the accretor, the star will temporarily fall out of thermal equilibrium. The radial response of the accretor star will depend on the structure of the envelope (as discussed for donor stars in Sect. 2.2.4). A star with a radiative envelope is expected to expand upon mass accretion, whereas a star with a convective envelope shrinks. In the former case, the accretor may swell up sufficiently to fill its Roche lobe, leading to the formation of a contact binary.

2.2.7 Supernova explosions in binaries

If the collapsing star is part of a binary or triple, natal kick alters the orbit and it can even unbind the system. Under the assumption that the SN is instantaneous and the SN-shell does not impact the companion star(s), the binary orbit is affected by the mass loss and velocity kick (Hills, 1983; Kalogera, 1996; Tauris and Takens, 1998; Pijloo et al., 2012) through:


where and are the semi-major axis of the pre-SN and post-SN orbit, is the mass lost by the collapsing star, is the total mass of the system pre-SN, is the pre-SN distance between the two stars, is the pre-SN relative velocity of the collapsing star relative to the companion, and


is the orbital velocity in a circular orbit. A full derivation of this equation and that for the post-SN eccentricity is given in Appendix A.1. Note that the equation for the post-SN eccentricity of Pijloo et al. (their Eq.8a 2012) is incomplete.

Eq. 21 shows that with a negligible natal kick, a binary survives the supernova explosion if less than half of the mass is lost. Furthermore, the binary is more likely to survive if the SN occurs at apo-astron. With substantial natal kicks compared to the pre-SN orbital velocity, survival of the binary depends on the magnitude ratio and angle between the two (through in Eq. 21). Furthermore, the range of angles that lead to survival is larger at peri-astron than apo-astron (Hills, 1983). If the direction of the natal kick is opposite to the orbital motion of the collapsing star, the binary is more likely to survive the SN explosion.

2.3 Triple evolution

The structures of observed triples tend to be hierarchical, i.e. the triples consist of an inner binary and a distant star (hereafter outer star) that orbits the centre of mass of the inner binary (Hut and Bahcall, 1983). To define a triple star system, no less than 10 parameters are required (Tbl. 1):
- the masses of the stars in the inner orbit and , and the mass of the outer star in the outer orbit ;
- the semi-major axis , the eccentricity , the argument of pericenter of both the inner and outer orbits. Parameters for the inner and outer orbit are denoted with a subscript ’’ and ’’, respectively;
- the mutual inclination between the two orbits.
The longitudes of ascending nodes specify the orientation of the triple on the sky, and not the relative orientation. Therefore, they do not affect the intrinsic dynamical evolution. From total angular momentum conservation for a reference frame with the z-axis aligned along the total angular momentum vector (Naoz et al., 2013).

In some cases, the presence of the outer star has no significant effect on the evolution of the inner binary, such that the evolution of the inner and outer binary can be described separately by the processes described in Sect. 2.1 and 2.2. In other cases, there is an interaction between the three stars that is unique to systems with multiplicities of higher orders than binaries. In this way, many new evolutionary pathways open up compared to binary evolution. The additional processes are described in the following sections, such as the dynamical instability and Lidov-Kozai cycles.

2.3.1 Stability of triples

The long-term behaviour of triple systems has fascinated scientists for centuries. Not only stellar triples have been investigated, but also systems with planetary masses, such as the Earth-Moon-Sun system by none other than Isaac Newton. It was soon realised that the three-body problem does not have closed-form solutions as in the case for two-body systems. Unstable systems dissolve to a lower order systems on dynamical timescales (van den Berk et al., 2007).

It is hard to define the boundary between stable and unstable systems, as stability can occur on a range of timescales. Therefore, many stability criteria exist (Mardling, 2001; Georgakarakos, 2008), that can be divided in three categories: analytical, numerical integration and chaotic criteria. The commonly used criterion of Mardling and Aarseth (1999):


where systems are unstable if and . This criterion is based on the concept of chaos and the consequence of overlapping resonances. The criterion is conservative, as the presence of chaos in some cases is not necessarily the same as an instability. By comparison with numerical integration studies, it was shown that Eq. 23 works well for a wide range of parameters (Aarseth and Mardling, 2001; Aarseth, 2004).

Most observed triples have hierarchical structures, because democratic triples tend to be unstable and short-lived (van den Berk et al., 2007). Hierarchical triples that are born in a stable configuration can become unstable as they evolve. Eq. 23 shows that when the ratio of the semi-major axes of the outer and inner orbit decreases sufficiently, the system enters the instability regime. Physical mechanisms that can lead to such an event, are stellar winds from the inner binary and stable mass transfer in the inner binary (Kiseleva et al., 1994; Iben and Tutukov, 1999; Freire et al., 2011; Portegies Zwart et al., 2011). Regarding wind mass losses from the inner binary exclusively, the fractional mass losses . Therefore, the fractional orbital increases , following Eq. 5. Perets and Kratter (2012) shows that such a triple evolution dynamical instability (TEDI) lead to close encounters, collisions, and exchanges between the stellar components. They find that the TEDI evolutionary channel caused by stellar winds is responsible for the majority of stellar collisions in the Galactic field.

2.3.2 Lidov-Kozai mechanism

Secular dynamics can play a mayor role in the evolution of triple systems. The key effect is the Lidov-Kozai mechanism (Lidov, 1962; Kozai, 1962), see Sect. 4.1 for an example of a triple undergoing Lidov-Kozai cycles. Due to a mutual torque between the inner and outer binary orbit, angular momentum is exchanged between the orbits. The orbital energy is conserved, and therefore the semi-major axes are conserved as well (e.g. Mardling and Aarseth, 2001). As a consequence, the orbital inner eccentricity and mutual inclination vary periodically. The maximum eccentricity of the inner binary is reached when the inclination between the two orbits is minimized. Additionally, the argument of pericenter may rotate periodically (also known as precession or apsidal motion) or librate. For a comprehensive review of the Lidov-Kozai effect, see Naoz (2016).

The Lidov-Kozai mechanism is of great importance in several astrophysical phenomena. For example, it can play a mayor role in the eccentricity and obliquity of exoplanets (e.g. Holman et al., 1997; Veras and Ford, 2010; Naoz et al., 2011) including high-eccentricity migration to form hot Jupiters (e.g. Wu and Murray, 2003; Correia et al., 2011; Petrovich, 2015), and for accretion onto black holes in the context of tidal disruption events (e.g. Chen et al., 2009; Wegg and Nate Bode, 2011) or mergers of (stellar and super-massive) black hole binaries (e.g. Blaes et al., 2002; Miller and Hamilton, 2002; Antonini et al., 2014). In particular for the evolution of close binaries, the Lidov-Kozai oscillations may play a key role (e.g. Harrington, 1969; Mazeh and Shaham, 1979; Kiseleva et al., 1998; Fabrycky and Tremaine, 2007; Naoz and Fabrycky, 2014), e.g. for black hole X-ray binaries (Ivanova et al., 2010), blue stragglers (Perets and Fabrycky, 2009), and supernova type Ia progenitors (Thompson, 2011; Hamers et al., 2013).

When the three-body Hamiltonian is expanded to quadrupole order in , the timescale for the Lidov-Kozai cycles is (Kinoshita and Nakai, 1999):


where and are the periods of the inner and outer orbit, respectively. The dimensionless quantity depends weakly on the mutual inclination, and on the eccentricity and argument of periastron of the inner binary, and is of order unity (Antognini, 2015). The timescales are typically much longer than the periods of the inner and outer binary.

Within the quadrupole approximation, the maximum eccentricity is a function of the initial mutual inclination as (Innanen et al., 1997):


in the test-particle approximation (Naoz et al., 2013), i.e. nearly circular orbits (, ) with one of the inner two bodies a massless test particle () and the inner argument of pericenter . In this case, the (regular) Lidov-Kozai cycles only take place when the initial inclination is between 39.2-140.8. For larger inner eccentricities, the range of initial inclinations expands.

For higher orders of i.e. the octupole level of approximation, even richer dynamical behaviour is expected than for the quadrupole approximation (e.g. Ford et al., 2000; Blaes et al., 2002; Lithwick and Naoz, 2011; Naoz et al., 2013; Shappee and Thompson, 2013; Teyssandier et al., 2013). The octupole term is non-zero when the outer orbit is eccentric or if the stars in the inner binary have unequal masses. Therefore it is often deemed the ’eccentric Lidov-Kozai mechanism’. In this case the z-component of the angular momentum of the inner binary is no longer conserved. It allows for a flip in the inclination such that the inner orbit flips from prograde to retrograde or vice versa (hereafter ’orbital flip’). Another consequence of the eccentric Lidov-Kozai mechanism is that the eccentricity of the inner binary can be excited very close to unity. The octupole parameter measure the importance of the octupole term compared to the quadropole term, and is defined by:


Generally, when , the eccentric Lidov-Kozai mechanism can be of importance (Naoz et al., 2011; Shappee and Thompson, 2013).

The dynamical behaviour of a system undergoing regular or eccentric Lidov-Kozai cycles can lead to extreme situations. For example, as the eccentricity of the inner orbit increases, the corresponding pericenter distance decreases. The Lidov-Kozai mechanism is therefore linked to a possible enhanced rate of grazing interactions, physical collisions, and tidal disruptions events of one of the stellar components (Ford et al., 2000; Thompson, 2011), and to the formation of eccentric semi-detached binaries (Sect. 2.3.6).

2.3.3 Lidov-Kozai mechanism with mass loss

Eqs. 24 and 26 show that the relevance of the Lidov-Kozai mechanism for a specific triple strongly depends on the masses and mass ratios of the stellar components. If one of the components loses mass, the triple can change from one type of dynamical behaviour to another type. For example, mass loss from one of the stars in the inner binary, can increase significantly. As a result the triple can transfer from a regime with regular Lidov-Kozai cycles to a regime where the eccentric Lidov-Kozai mechanism is active. This behaviour is known as mass-loss induced eccentric Kozai (MIEK) (Shappee and Thompson, 2013; Michaely and Perets, 2014). See also Sect. 4.3 for an example of this evolutionary pathway.

The inverse process (inverse-MIEK), when a triple changes state from the octupole to the quadrupole regime, can also occur. Eq. 26 shows this is the case when mass loss in the inner binary happens to create an fairly equal mass binary, or when the semi-major axis of the outer orbit increases. This latter is possible when the outer star loses mass in a stellar wind (Sect. 2.2.1).

Another example comes from Michaely and Perets (2014), who studied the secular freeze-out (SEFO). In this scenario mass is lost from the inner binary such that the Lidov-Kozai timescale increases (Eq. 24). This induces a regime change from the quadrupole regime, to a state where secular evolution is either quenched or operates on excessively long time-scales.

The three examples given above illustrate that the dynamical evolution of a triple system is intertwined with the stellar evolution of its components. Thus, in order to gain a clear picture of triple evolution, both three-body dynamics and stellar evolution need to be taken into account simultaneously.

2.3.4 Precession

Besides precession caused by the Lidov-Kozai mechanism, other sources of precession exist in stellar triples. These include general relativistic effects (Blaes et al., 2002):


Furthermore, orbital precession can be caused by the distortions of the individual stars by tides (Smeyers and Willems, 2001):


and by intrinsic stellar rotation (Fabrycky and Tremaine, 2007):


where is the mass of the distorted star that instigates the precession and the companion star in the two-body orbit. The distorted star has a classical apsidal motion constant , radius , and a spin frequency . The precession rates in Eq. 27, as well as Eq. 28 and Eq. 29, are always positive. This implies that relativistic effects, tides and stellar rotation mutually stimulate precession in one direction. Note, that precession due to these processes also take place in binaries, which affects the binary orientation, but not the evolution of the system.

If the timescales555See for example Eq. 8 and 9. The timescales for GW emission are shorter for binaries with small separations and large eccentricities. for these processes become comparable or smaller than the Lidov-Kozai timescales, the Lidov-Kozai cycles are suppressed. Because Lidov-Kozai cycles are driven by tidal forces between the outer and inner orbit, the additional precession tends to destroy the resonance (Liu et al., 2015a). As a result of the suppression of the cycles, the growth of the eccentricity is limited, and orbital flips are limited to smaller ranges of the mutual inclination (Naoz et al., 2012; Petrovich, 2015; Liu et al., 2015a).

2.3.5 Tides and gravitational waves

As mentioned earlier, the Lidov-Kozai mechanism can lead to very high eccentricities that drives the stars of the inner binary close together during pericenter passage. During these passages, tides and GW emission can effectively alter the orbit (Mazeh and Shaham, 1979; Kiseleva et al., 1998). Both processes are dissipative, and act to circularize the orbit and shrink the orbital separation (Sect. 2.2.2 and 2.2.3). The combination of Lidov-Kozai cycles with tides or GW emission can then lead to an enhanced rate of mergers and RLOF. For GW sources, the merger time of a close binary can be significantly reduced, if an outer star is present that gives rise to Lidov-Kozai cycles on a short timescale (Thompson, 2011). This is important in the context of supernova type Ia and gamma-ray bursts.

The combination of Lidov-Kozai cycles with tidal friction (hereafter LKCTF) can also lead to an enhanced formation of close binaries (Mazeh and Shaham, 1979; Kiseleva et al., 1998). This occurs when a balance can be reached between the eccentricity excitations of the (regular or eccentric) Lidov-Kozai mechanism and the circularisation due to tides666A balance is also possible between the eccentricity excitations of the Lidov-Kozai mechanism and other sources of precession, see Sect.2.3.4). . The significance of LKCTF is illustrated by Fabrycky and Tremaine (2007), who show that MS binaries with orbital periods of 0.1–10d are produced from binaries with much longer periods up to d. Observationally, 96% of close low-mass MS binaries are indeed part of a triple system (Tokovinin et al., 2006). Several studies of LKCTF for low-mass MS stars exist (Mazeh and Shaham, 1979; Eggleton and Kiseleva-Eggleton, 2001; Fabrycky and Tremaine, 2007; Kisseleva-Eggleton and Eggleton, 2010; Hamers et al., 2013), however, a study of the effectiveness of LKCTF for high-mass MS triples or triples with more evolved components is currently lacking. Due to the radiative envelopes of high-mass stars, LKCTF is likely less effective compared to the low-mass MS case. However, evolved stars develop convective envelopes during the giant phases for which tidal friction is expected to be effective. Hence, in order to understand the full significance of LKCTF for triple evolution, it is necessary to model three-body dynamics and stellar evolution consistently.

2.3.6 Mass transfer initiated in the inner binary…

In Sect. 2.2.4, we described the effect of mass transfer on a circularized and synchronized binary. However, as Lidov-Kozai cycles can lead effectively to RLOF in eccentric inner binaries, the simple picture of synchronization and circularisation before RLOF, is no longer generally valid for triples. In an eccentric binary, there does not exist a frame in which all the material is corotating, and the binary potential becomes time-dependent. Studies of the Roche lobe for eccentric and/or asynchronous binaries, show that the Roche lobe can be substantially altered (Plavec, 1958; Regös et al., 2005; Sepinsky et al., 2007a). In an eccentric orbit, the Roche lobe of a star at periastron may be significantly smaller than that in a binary that is circularized at the same distance . The Roche lobe is smaller for stars that rotate super-synchronously at periastron compared to the classical Roche lobe (Eq. 4), and larger for sub-synchronous stars. It is even possible that the Roche lobe around the accretor star opens up. When mass is transferred from the donor star through L1, it is not necessarily captured by the accretor star, and mass and angular momentum may be lost from the binary system.

The modification of the Roche lobe affects the evolution of the mass transfer phase, e.g. the duration and the mass loss rate. Mass transfer in eccentric orbits of isolated binaries has been studied in recent years with SPH techniques (Layton et al., 1998; Regös et al., 2005; Church et al., 2009; Lajoie and Sills, 2011; van der Helm et al., 2016) as well as analytical approaches (Sepinsky et al., 2007b, 2009, 2010; Davis et al., 2013; Dosopoulou and Kalogera, 2016a, b). These studies have shown that (initially) the mass transfer is episodic. The mass transfer rate peaks just after periastron, and its evolution during the orbit shows a Gaussian-like shape with a FWHM of about 10% of the orbital period.

The long-term evolution of eccentric binaries undergoing mass transfer can be quite different compared to circular binaries. The long-term evolution has been studied with analytics adopting a delta-function for the mass transfer centred at periastron (Sepinsky et al., 2007b, 2009, 2010; Dosopoulou and Kalogera, 2016a, b). Under these assumptions, the semi-major axis and eccentricity can increase as well as decrease depending on the properties of the binary at the onset of mass transfer. In other words, the secular effects of mass transfer can enhance and compete with the orbital effects from tides. Therefore, rapid circularization of eccentric binaries during the early stages of mass transfer is not generally justified.

The current theory of mass transfer in eccentric binaries predicts that some binaries can remain eccentric for long periods of time. The possibility of mass transfer in eccentric binaries is supported by observations. For example, the catalogue of eccentric binaries of Petrova and Orlov (1999) contains 19 semi-detached and 6 contact systems out of 128 systems. That circularisation is not always achieved before RLOF commences, is supported by observations of some detached, but close binaries e.g. ellipsoidal variables (Nicholls and Wood, 2012) and Be X-ray binaries, which are a subclass of high-mass X-ray binaries (Raguzova and Popov, 2005).

2.3.7 … and its effect on the outer binary

Mass transfer in the inner binary can affect the triple as a whole. The most simple case is during conservative stable mass transfer, when the outer orbit remains unchanged. If the inner orbit is circularized and synchronised, mass transfer generally leads to an increase in the inner semi-major axis by a factor of a few (Eq. 18). When the ratio decreases, the triple approaches and possibly crosses into the dynamically unstable regime (Sect. 2.3.1 and Eq. 23).

If the mass transfer in the inner binary occurs non-conservatively, the effect on the outer binary is completely determined by the details of the mass loss from the inner binary. We conceive three scenarios for this to take place. First, if during stable mass transfer to a hydrogen-rich star, matter escapes from the inner binary, it is likely the matter will escape from L2 in the direction of the orbital plane of the inner binary. Second, during mass transfer to a compact object, a bipolar outflow or jet may develop. Thirdly, matter may be lost from the inner binary as a result of a common-envelop phase, which we will discuss in Sect. 2.3.8.

If matter escapes the inner binary, its velocity must exceed the escape velocity:


and analogously, to escape from the outer binary, and the triple as a whole:


For stable triples, such that , unless . The factor is of the order of one, e.g. for circular orbits . In the catalogue of Tokovinin (2014b) it is uncommon that . Out of 199 systems, there are 3 systems with and none with . Therefore, if the inner binary matter is energetic enough to escape from the inner binary, it is likely to escape from the triple as a whole as well.

For isolated binary evolution, it is unclear if the matter that leaves both Roche lobes is energetic enough to become unbound from the system, e.g. when mass is lost through L2. Instead a circumbinary disk may form that gives rise to a tidal torque between the disk and the binary. This torque can efficiently extract angular momentum from the binary i.e. Eq. 19 with , where is the radius of the circumbinary ring (Soberman et al., 1997). For example, for a binary with , angular momentum is extracted more than 10 faster if the escaping matter forms a circumbinary disk compared to a fast stellar wind (Soberman et al., 1997; Toonen et al., 2014). Hence, the formation of a circumbinary disk leads to a stronger reduction of the binary orbit, and possibly a merger.

If a circumbinary disk forms around the inner binary of a triple, we envision two scenarios. Firstly, the outer star may interact with the disk directly if its orbit crosses the disk. Secondly, the disk gives rise to two additional tidal torques, with the inner binary and the outer star. It has been shown that the presence of a fourth body can lead to a suppression of the Lidov-Kozai cycles, however, bodies less massive than a Jupiter mass have a low chance of shielding (Hamers et al., 2015, 2016; Martin et al., 2015; Muñoz and Lai, 2015).

2.3.8 The effect of common-envelope on the outer binary

Another scenario exists in which material is lost from the inner binary; in stead of a stable mass transfer phase, mass is expelled through a CE-phase. For isolated binaries the CE-phase and its effect on the orbit is an unsolved problem, and the situation becomes even more complicated for triples. In the inner binary the friction between the stars and the material is expected to cause a spiral-in. If the outer star in the triple is sufficiently close to the inner binary, the matter may interact with the outer orbit such that a second spiral-in takes place. If on the other hand, the outer star is in a wide orbit, and the CE-matter is lost in a fast and isotropic manner, the effect on the outer orbit would be like a stellar wind (Sect. 2.2.1). Veras and Tout (2012) study the effect of a CE in a binary with a planet on a wider orbit. Assuming that the CE affects the planetary orbit as an isotropic wind, they find that planetary orbits of a few  are readily dissolved. In this scenario the CE-phase operates virtually as a instantaneous mass loss event, and therefore the maximum orbital separation for the outer orbit to remain bound is strongly dependent on the uncertain timescale of the CE-event (Sect. 2.2.1). Disruption of a triple due to a CE-event may also apply to stellar triples, however, the effect is likely less dramatic, as the relative mass lost in the CE-event to the total system mass is lower.

Note that most hydrodynamical simulations of common-envelope evolution show that matter is predominantly lost in the orbital plane of the inner binary (e.g. Ricker and Taam, 2012; Passy et al., 2012b), however, these simulations are not been able to unbind the majority of the envelope. In contrast, in the recent work of Nandez et al. (2015), the envelope is expelled successfully due to the inclusion of recombination energy in the equation of state. These simulations show a more spherical mass loss. Roughly 60% of the envelope mass is ejected during the spiral-in phase in the orbital plane, while the rest of the mass is ejected after the spiral-in phase in a closely spherical way (priv. comm. Jose Nandez).

The first scenario of friction onto the outer orbit has been proposed to explain the formation of two low-mass X-ray binaries with triple components (4U 2129+47 (V1727 Cyg) (Portegies Zwart et al., 2011) and PSR J0337+1715 (Tauris and van den Heuvel, 2014)), however, it could not be ruled out that the desired decrease in the outer orbital period did not happen during the SN explosion in which the compact object was formed. Currently, it is unclear if the CE-matter is dense enough at the orbit of the outer star to cause significant spiral-in. Sabach and Soker (2015) suggests that if there is enough matter to bring the outer star closer, the CE-phase would lead to a merger in the inner binary.

2.3.9 Mass loss from the outer star

In about 20% of the multiples in the Tokovinin catalogue of multiple star systems in the Solar neighbourhood, the outer star is more massive than the inner two stars. For these systems the outer star evolves faster than the other stars (Sect. 2.1). In about 1% of the triples in the Tokovinin catalogue, the outer orbit is significantly small that the outer star is expected to fill its Roche lobe at some point in its evolution (de Vries et al., 2014).

What happens next has not been studied to great extent, i.e. the long-term evolution of a triple system with a mass-transferring outer star. It is an inherently complicated problem where the dynamics of the orbits, the hydrodynamics of the accretion stream and the stellar evolution of the donor star and its companion stars need to be taken into account consistently.

Such a phase of mass transfer has been invoked to explain the triple system PSR J0337+1715, consisting of a millisecond pulsar with two WD companions (Tauris and van den Heuvel, 2014; Sabach and Soker, 2015). Tauris and van den Heuvel (2014) note one of the major uncertainties in their modelling of the evolution of PSR J0337+1715 comes from the lack of understanding of the accretion onto the inner binary system and the poorly known specific orbital angular momentum of the ejected mass during the outer mass transfer phase. Sabach and Soker (2015) proposes that if the inner binary spirals-in to the envelope of the expanding outer star, the binary can break apart from tidal interactions.

To the best of our knowledge, only de Vries et al. (2014) have performed detailed simulations of mass transfer initiated by the outer star in a triple. They use the same software framework (AMUSE, Sect. 3) as we use for our code TrES . de Vries et al. (2014) simulate the mass transfer phase initiated by the outer star for two triples in the Tokovinin catalogue, Tau and HD97131. For both systems, they find that the matter lost by the outer star does not form an accretion disk or circumbinary disk, but instead the accretion stream intersects with the orbit of the inner binary. The transferred matter forms a gaseous cloud-like structure and interacts with the inner binary, similar to a CE-phase. The majority of the matter is ejected from the inner binary, and the inner binary shrinks moderately to weakly with depending on the mutual inclination of the system. In the case of HD 97131, this contraction leads to RLOF in the inner binary. The vast majority of the mass lost by the donor star is funnelled through L1, and eventually ejected from the system by the inner binary through the L3 Lagrangian point777Here the L3 Lagrangian point is located behind the inner binary on the line connecting the centres of mass of the outer star and inner binary. of the outer orbit. As a consequence of the mass and angular momentum loss, the outer orbit shrinks with888There is an inconsistency in the meaning of between Eqs.4 and 5 in de Vries et al. (2014), denoted as in their equations. In their fits represents the ratio , where is the amount of angular momentum that is lost from the system, and the orbital angular momentum. It does not represent where is the angular momentum of the accretor star. 3–4 in Eq. 19. During the small number of outer periods that are modelled, the inner and outer orbits approaches contraction at the same fractional rate. Therefore the systems remain dynamically stable.

Systems that are sufficiently wide that the outer star does not fill its Roche lobe, might still be affected by mass loss from the outer star in the form of stellar winds. Soker (2004) has studied this scenario for systems where the outer star is on the AGB, such that the wind mass loss rates are high. Assuming Bondi-Hoyle-Littleton accretion and


where is the width of the accretion column at the binary location, and the Bondi-Hoyle accretion radius, Soker (2004) finds that a large fraction of triples the stars in the inner binary may accrete from an accretion disk around the stars. The formation of an accretion disk depends strongly on the orientation of the inner and outer orbit. When the inner and outer orbits are parallel to each other, no accretion disk forms. On the other hand, when the inner orbit is orientated perpendicular to the outer orbit, an accretion disk forms if . In this case the accretion is in a steady-state and mainly towards the most massive star of the inner binary.

2.3.10 Triples and planetary nebulae

Interesting to mention in the context of triple evolution are planetary nebulae (PNe), in particular those with non-spherical structures. The formation of these PNe is not well understood, but maybe attributed to interactions between an AGB-star and a companion (e.g. Bond and Livio, 1990; Bond, 2000; De Marco et al., 2015; Zijlstra, 2015) or multiple companions (e.g. Bond et al., 2002; Exter et al., 2010; Soker, 2016; Bear and Soker, 2016). Where a binary companion can impose a non-spherical symmetry on the resulting PN, and even a non-axisymmetry (see e.g. Soker and Rappaport, 2001, for eccentric binaries), triple evolution can impose structures that are not axisymmetric, mirrorsymmetric, nor pointsymmetric (Bond et al., 2002; Exter et al., 2010; Soker, 2016).

Since the centers of many elliptical and bipolar PNe host close binaries, the systems are expected to have undergone a CE-phase. In the context of triples, PNe formation channels have been proposed that concern outer stars on the AGB whose envelope matter just reaches or completely engulfs a tight binary system, e.g. the PN SuWt 2 (Bond et al., 2002; Exter et al., 2010). Another proposed channel involves systems with a very wide outer orbit of tens to thousands of AU and in which the outer star interacts with the material lost by the progenitor-star of the PN (Soker et al., 1992; Soker, 1994). For a detailed review of such evolutionary channels, see Soker (2016). Under the assumptions that PN from triple evolutionary channels give rise to irregular PNe, Soker (2016) and Bear and Soker (2016) find that about 1 in 6-8 PNe might have been shaped by an interaction with an outer companion in a triple system.

2.3.11 Supernova explosions in triples

Pijloo et al. (2012) study the effect of a supernova explosion in a triple star system under the same assumptions as Hills (1983). The authors show that for a hierarchical triple in which the outer star collapses in a SN event, the inner binary is not affected, and the effect on the outer orbit can be approximated by that of an isolated binary. For a SN taking place in the inner binary, the inner binary itself is modified similar to an isolated binary (Sect. 2.2.7). The effect on the outer binary can be viewed as that of an isolated binary in which the inner binary is replaced by an effective star at the center of mass of the inner binary. The effective star changes mass and position (as the center of mass changes) instantaneously in the SN event. The semi-major axis of the outer orbit is affected by the SN in the inner binary as:


where is the total mass of the pre-SN triple, and are the pre-SN and post-SN distance between the star in the outer orbit and the center of mass of the inner binary, is the pre-SN relative velocity of the center of mass of the inner binary relative to the outer star, is the systemic velocity the inner binary due to the SN, and


the orbital velocity in a circular orbit. Note that in comparison with the circular velocity in binaries (Eq. 22), here refers to and .

A full derivation of the change in semi-major axis (Eq. 33) and eccentricity (Eq. 78) of the outer orbit due to a SN in the inner orbit, are given in Appendix A.1. Note that the equation for the post-SN eccentricity of the outer orbit in Pijloo et al. (2012, their Eq.27) is incomplete (see Appendix A.1).

2.4 Quadruples and higher-order hierarchical systems

Although quadruples star systems are less common than triple systems, hierarchical quadruples still comprise about 1% of F/G dwarf systems in the field (Tokovinin, 2014a, b). While for triples, there is one type of hierarchy that is stable on long-terms (compared to stellar lifetimes), quadruples can be arranged in two distinct long-term stable configurations: the ‘2+2’ or ‘binary-binary’ configuration, and the ‘3+1’ or ‘triple-single’ configuration. In the first case, two binaries orbit each others barycentre, and in the latter case a hierarchical triple is orbited by a fourth body. In the sample of F/G dwarfs of Tokovinin (2014a, b), the ‘2+2’ systems comprise about 2/3 of quadruples, and the ‘3+1’ about 1/3.

The secular dynamics of the ‘2+2’ systems were investigated by Pejcha et al. (2013) using N-body methods. Pejcha et al. (2013) find that in these systems, orbital flips and associated high eccentricities are more likely compared to equivalent triple systems (i.e. with the companion binary viewed as a point mass). The ‘3+1’ configuration was studied by Hamers et al. (2015) 999Hamers et al. (2015) derive the orbit-averaged Hamiltonian expressions for the ‘3+1’ as well as the ‘2+2’ configuration.. For highly hierarchical systems, i.e. in which the three binaries are widely separated, the global dynamics can be qualitatively described in terms of the (initial) ratio of the Lidov-Kozai time-scales of the two inner most binaries compared to that of the outer two binaries. This was applied to the ‘3+1’ F/G systems of Tokovinin (2014a, b), and most (90%) of these systems were found to be in a regime in which the inner three stars are effectively an isolated triple, i.e. the fourth body does not affect the secular dynamical evolution of the inner triple.

We note that in the case of ‘3+1’ quadruples and with a low-mass third body (in particular, a planet), the third body can affect the Lidov-Kozai cycles that would otherwise have been induced by the fourth body. In particular, under specific conditions the third body can ‘shield’ the inner binary from the Lidov-Kozai oscillations, possibly preventing the inner binary from shrinking due to tidal dissipation, and explaining the currently observed lack of circumbinary planets around short-period binaries (Martin et al., 2015; Muñoz and Lai, 2015; Hamers et al., 2016). A similar process could apply to more massive third bodies, e.g. low-mass MS stars.

It is currently largely unexplored how non-secular effects such as stellar evolution, tidal evolution and mass transfer affect the evolution of hierarchical quadruple systems, or, more generally, in higher-order multiple systems. The secular dynamics of the latter could be efficiently modelled using the recent formalism of Hamers and Portegies Zwart (2016).

3 Methods

In the previous section, we gave an overview of the most important ingredients of the evolution of stars in single systems, binaries and triples. For example, nuclear evolution of a star leads to wind mass loss, that affects the dynamics of binaries and triples, and can even lead to a dynamical instability in multiple systems. Three-body dynamics can give rise to oscillations in the eccentricity of the inner binary system of the triple, which can lead to an amplified tidal effect and an enhanced rate of stellar mass transfer, collisions, and mergers. Additionally, a triple system can transition from one to another dynamical regime (i.e. without Lidov-Kozai cycles, regular and eccentric Lidov-Kozai cycles) due to stellar evolution, e.g. wind mass loss or an enhancement of tides as the stellar radius increases in time. These examples illustrate that for the evolution of triple stars, stellar evolution and dynamics are intertwined. Therefore, in order to study the evolution of triple star systems consistently, three-body dynamics and stellar evolution need to be taken into account simultaneously.

In this paper, we present a public source code TrES  to simulate the evolution of wide and close, interacting and non-interacting triples consistently. The code is designed for the study of coeval, dynamically stable, hierarchical, stellar triples. The code is based on heuristic recipes that combine three-body dynamics with stellar evolution and their mutual influences. These recipes are described here.

The code can be used to evaluate the distinct evolutionary channels of a specific population of triples or the importance of different physical processes in triple evolution. As an example, it can be used to assess the occurrence rate of stable and unstable mass transfer initiated in circular and eccentric inner orbits of triple systems (Toonen, Hamers, & Portegies Zwart in prep.). We stress that modelling though a phase of stable mass transfer in an eccentric orbit is currently not implemented in TrES , but we aim to add this to the capabilities of TrES  in a later version of the code.

The code TrES  is based on the secular approach to solve the dynamics (Sect. 3.3) and stellar evolution is included in a parametrized way through the fast stellar evolution code SeBa (Sect. 3.2). TrES  is written in the Astrophysics Multipurpose Software Environment, or AMUSE (Portegies Zwart et al., 2009; Portegies Zwart, 2013). This is a component library with a homogeneous interface structure based on Python. AMUSE can be downloaded for free at and In the AMUSE framework new and existing code from different domains (e.g. stellar dynamics, stellar evolution, hydrodynamics and radiative transfer) can be easily used and coupled. As a result of the easy coupling, the triple code can be easily extended to include a detailed stellar evolution code (i.e. that solve the stellar structure equations) or a direct N-body code to solve the dynamics of triples that are unstable or in the semi-secular regime (Sect. 3.3.1).

3.1 Structure of TrES

The code consist of three parts:
step 1) stellar evolution,
step 2) stellar interaction,
step 3) orbital evolution.
At the beginning of each timestep we estimate an appropriate timestep and evolve the stars as single stars for this timestep (step 1). The trial timestep is estimated with:


where , , and are the minimum timesteps due to stellar evolution, stellar wind mass losses, stellar radius changes and the previous timestep. Each star gives rise to a single value for , where the minimum is adopted as a trial timestep in TrES . The timestep is determined internally by the stellar evolution code (SeBa, Sect.3.2). It is the maximum attainable timestep for the next iteration of this code and is mainly chosen such that the stellar masses that evolve due to winds, are not significantly affected by the timesteps. Furthermore, when a star changes its stellar type (e.g. from a horizontal branch star to an AGB star), the timestep is minimized to ensure a smooth transition. For TrES , we require a more strict constraint on the wind mass losses, such that , where and is the wind mass loss rate given by the stellar evolution code. The numerical factor establishes a maximum average of 1% mass loss from stellar winds per timestep. Furthermore, we ensure that the stellar radii change by less then a percent per timestep through , where and are numerical factors. We take and


where and represent the time derivative of the radius of the current and previous timestep, respectively. This limit is particularly important since the degree of tidal interaction strongly depends on the stellar radius. Lastly, we require that , where is the previous successfully accomplished timestep and a numerical factor with a value of 100.

The trial timestep is accepted and the code continues to step 2, only if:
case a) no star has started to fill its Roche lobe,
case b) stellar radii have changed by less than 1%,
case c) stellar masses have changed by less than 5%
within the trial timestep. We have tested that these percentages give accurate results with respect to the orbital evolution. Condition b is not applied at moments when the stellar radius changes discontinuously, such as during the helium flash or at white dwarf formation.

Conditions b & c are not applied, when a massive star collapses to a neutron star or black hole. Note that when a star undergoes such a supernova explosion, the timestep is minimized through . Additionally step 2 and 3 are skipped and the triple is adjusted according to Sect. 3.4.7.

If conditions a, b and c are not met, the timestep is reverted and step 1 is tried again with a smaller timestep . This process is done iteratively until the conditions are met or until the timestep is sufficiently small, yr. If the change in the stellar parameters is too large (i.e. case b & c), the new trail timestep is taken to be:


where represents 90% of the previous timestep and is a newly calculated timestep according to Eq. 35 for which the time derivative of the radius from the last trial timestep is used.

During mass transfer, the timestep is estimated by:


where is the mass loss rate from mass transfer (Sect. 3.4.5), and the numerical factors and . If a star starts filling its Roche lobe in a timestep, . If , mass transfer is allowed to commence.

Step 2 in our procedure regards the modelling of the stellar interactions such as stable mass transfer, contact evolution and common-envelope evolution (Sect. 3.4).

The last step involves the simulation of the orbital evolution of the system by solving a system of differential equations (Sect. 3.3). If the evolution leads to the initiation of RLOF during the trial timestep, both the orbit and stellar evolution are reverted to the beginning of the timestep. If the time until RLOF is shorter than 1% of , the latter is taken to be the new trial timestep and mass transfer is allowed to commence. If not, the timestep is taken to be the time until RLOF that was found during the last trial timestep. If during the orbital evolution the system becomes dynamically unstable, the simulation is terminated. The stability criterion of Mardling and Aarseth (2001) is used (Eq.23). In all other cases, the trial timestep is accepted, and the next iteration begins.

3.2 Stellar evolution

Single stellar evolution101010Note that SeBa is not used to model binary evolution in TrES . is included through the fast stellar evolution code SeBa (Portegies Zwart and Verbunt, 1996; Nelemans et al., 2001; Toonen et al., 2012; Toonen and Nelemans, 2013). SeBa is a parametrized stellar evolution code providing parameters such as radius, luminosity and core mass as a function of initial mass and time. SeBa is based on the stellar evolution tracks from Hurley et al. (2000). These tracks are fitted to the results of a detailed stellar evolution code (based on Eggleton, 1971, 1972) that solves the stellar structure equations.

3.3 Orbital evolution

TrES  solves the orbital evolution through a system of first-order ordinary differential equations (ODE):


where , and are the orbital angular momentum of the inner and outer orbit, and ,  and  the spin frequency of the star with mass ,  and , respectively, and the moment of inertia of the corresponding star. represents the time derivative of parameter . Eq. 39 includes secular three-body dynamics (with subscript 3b), general relativistic effects (GR), tidal friction (TF), precession, stellar wind effects and mass transfer (MT). The quadrupole terms of the three-body dynamics are based on Harrington (1968), and the octupole terms on Ford et al. (2000) with the modification of Naoz et al. (2013). Gravitational wave emission is included as in Eq. 8 and 9 (Peters, 1964). Our treatment of tidal friction and precession is explained in Sect. 3.4.2. Magnetic braking is currently not included. The treatment of stellar winds and mass transfer is described in Sect. 3.4.1-3.4.5.

The ODE solver routine uses adaptive timesteps to simulate the desired timestep . Within the ODE solver, parameters that are not given in Eq. 39 (e.g. gyration radius), are assumed to be constant during . An exception to this is the stellar radius, mass and moment of inertia. Even though is chosen such that the parameters do not change significantly within this timestep (Sect. 3.1), there is a cumulative effect that can violate angular momentum conservation on longer timescales if is not taken into account. As a non-interacting star evolves and the mass, radius and moment of inertia change, the spin frequency of the star evolves accordingly due to conservation of spin angular momentum. The change in the spin frequency is:




where and are the mass and radius of the core, and (Hurley et al., 2000). Thus we approximate the moment of inertia with a component for the core and for the envelope of the star. This method works well for evolved stars that have developed dense cores, as well as for MS stars with , and compact objects for which .

The initial spin periods of the stellar components of the triple are assumed to be similar to that of ZAMS stars. Based on observed rotational velocities of MS stars from Lang (1992), Hurley et al. (2000) proposed the fit:


As in Hamers et al. (2013), we make the simplifying assumption that the stellar spin axes are aligned with the orbital axis of the corresponding star. For the vast majority of stellar triples the magnitude of the spin angular momenta are small compared to that of the orbital angular momenta. A consequence of this assumption, is that tidal friction from a spin-orbit misalignment is absent. The change in the mutual inclination in Eq. 39 is based on the conservation of total angular momentum (Hamers et al., 2013).

The set of orbital equations of Eq. 39 are solved by a routine based on the ODE solver routine presented in Hamers et al. (2013). This routine uses the CVODE library, which is designed to integrate stiff ODEs (Cohen et al., 1996). It has been verified by comparing integrations with example systems presented in Ford et al. (2000), Blaes et al. (2002) and Naoz et al. (2011), and comparing with analytical solutions at the quadrupole-order level of approximation assuming , given by Kinoshita and Nakai (1999).

The ODE consists of a combination of prescriptions for the main physical processes for triple evolution (e.g. three-body dynamics and tides) which are described in Sect. 2. In order to set up the ODE, we have assumed that the physical processes are independent of one another, such that their analytical treatments can be added linearly. Michaely and Perets (2014) show that the dynamics of a hierarchical triple including mass loss and transfer can be well modelled with this approach. They find that the secular approach shows excellent agreement with full N-body simulations. Additionally, we note that the ODE of Eq. 39 is valid as long as the included processes occur on timescales longer than the dynamical timescale. In the next section we discuss when this criterion is violated, and describe the alternative treatments in TrES .

3.3.1 The secular approach

The components in Eq. 39 (, , and ) that describe the secular three-body dynamics, including the regular Lidov-Kozai cycles, are derived using the orbital-averaging technique (e.g. Michaely and Perets, 2014; Naoz, 2016; Luo et al., 2016). With this method, the masses of the components of the triple system are distributed over the inner and outer orbit. The three-body dynamics is then approximated by the interaction between these ellipses. Furthermore, the Hamiltonian is expanded up to third order in (i.e. the octopole term). The quadrupole level of approximation refers to the second-order expansion, which is the lowest non-trivial expansion order. The quadrupole level is sufficient for systems in which the octupole parameter (Eq. 26) is sufficiently low, i.e. . For example, this includes systems with . To accurately model the dynamics of a population of systems with a wide range of mass ratios, eccentricities and orbital separations, one needs to include the octupole term as well.

However, whereas the orbit-averaged equations of motion are valid for strongly hierarchical systems, they become inaccurate for triples with weaker hierarchies (e.g. Antonini and Perets, 2012; Katz and Dong, 2012; Naoz et al., 2013; Antonini et al., 2014; Antognini et al., 2014; Bode and Wegg, 2014; Luo et al., 2016). For these systems, the timescale of the perturbations due to the outer star can become comparable to or shorter than the dynamical timescales of the system. This is problematic as the orbit-average treatment neglects any modulations on short orbital timescales per definition. For moderately hierarchical systems, the outer star can significantly change the angular momentum of the inner binary between two successive pericenter passages causing rapid oscillations in the corresponding eccentricity. The orbit-average treatment is valid when:


as derived by Antonini et al. (2014).

In the non-secular regime, the inner binary can be driven to much higher eccentricities than the secular approximation predicts, and subsequently lead to more collisions of e.g. black holes (Antonini et al., 2014), neutron stars (Seto, 2013) or white dwarfs (Katz and Dong, 2012). These are interesting in the context of gravitational wave emission and type Ia supernovae. Due to the very short timescale of the eccentricity oscillations, and therefore rapid changes of the periapse distance, tidal or general relativistic effects do not play a role (Katz and Dong, 2012). With the secular approach, as in TrES , the maximum inner eccentricity and therefore the number of collisions is probably underestimated for moderately hierarchical systems (see also Naoz et al., 2016).

Furthermore, Luo et al. (2016) showed recently that the rapid oscillations accumulate over time and alter the long-term evolution of the triple systems (e.g. whether or not an orbital flip occurs). The non-secular behaviour discussed by Luo et al. (2016) occurs in systems in which the mass of the outer star is comparable or larger than that of the inner binary, and in which the octupole term is important ().

3.4 Stellar interaction

3.4.1 Stellar winds in TrES

The mass loss rate of each star depends on the evolutionary stage and is determined by the stellar evolution code. We make the common assumption that the wind is fast and spherically symmetric with respect to the orbit, as discussed in Sect. 2.2.1. For the inner binary, we assume a fraction of the wind mass lost from at a rate can be accreted by , and for the mass flowing in the other direction. Following Eq. 7, the effect on the inner orbit is then: