The evolution of hierarchical triple star-systems

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 behavior 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 modeling 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.


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 ∼ 10% of the low-mass stars are in triples (Moe and Di Stefano, 2016;Raghavan et al., 2010;Tokovinin, 2008Tokovinin, , 2014b a fraction that gradually increases (Duchêne and Kraus, 2013) to ∼ 50% for spectral type B stars (Moe and Di Stefano, 2016;Remage Evans, 2011;Sana et al., 2014). 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 Fabrycky and Tremaine (2007); Ford et al. (2000); Liu et al. (2015a); Naoz and Fabrycky (2014); Naoz et al. (2013), 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 (Hamers et al., 2013;Kratter and Perets, 2012;Michaely and Perets, 2014;Naoz et al., 2016;Perets and Kratter, 2012;Shappee and Thompson, 2013), 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. 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.

Parameters
Stellar Orbital Single star m -Binary m 1 , m 2 a, e Triple m 1 , m 2 , m 3 i mutual , a in , e in , g in , h in aout, eout, gout, hout 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.

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.

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 overshooting [1] .
[1] Overshooting refers to a chemically mixed region beyond the boundary of the convective core (e.g. Maeder and Meynet, 1991;Massevitch et al., 1979;Stothers, 1963), 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).

Timescales
Fundamental timescales of stellar evolution are the dynamical (τ dyn ), thermal (τ th ), and nuclear timescale (τ nucl ). 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 R and m 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 L 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, c is the speed of light, and m nucl is the amount of mass available as fuel. For core hydrogen burning, = 0.007 and M nucl ≈ 0.1M . Assuming a mass-luminosity relation of L ∝ M α , with empirically α ≈ 3 − 4 (e.g. Eker et al., 2015;Salaris and Cassisi, 2005), it follows that massive stars live shorter and evolve faster than low-mass stars. For the Sun, τ dyn ≈ 30 min, τ th ≈ 30 Myr, and τ nucl ≈ 10 Gyr. Typically, τ dyn < τ th < τ nucl , which allows us to quantitatively predict the structure and evolution of stars in broad terms.

Hertzsprung-Russell diagram
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 Figure 1 Hertzsprung-Russell diagram Evolutionary tracks for seven stars in the HR-diagram with masses 1, 1.5, 2.5, 4, 6.5, 10, and 15M 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. 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.2M 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 13M reach effective temperatures as low as (10 3.7 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 10 8 K, helium core burning commences and the red giant phase has come to an end. For stars less massive than 2M , helium ignites degenerately in a helium flash. For stars more massive than 13M , 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 5 · 10 8 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 Evolution of stellar radius Radius as a function of stellar age for two stars with masses 4 and 6.5M 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.

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 stellarwind 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 ( 10 −7 − 10 −4 M yr −1 ) and low terminal velocities (5-30 kms −1 ). For stars up to 8M , 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 (10 −7 − 10 −4 M yr −1 ) and high velocities (several 100-1000 kms −1 ) (e.g. Kudritzki and Puls, 2000). For stars of more than ∼30M , 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.

Stellar remnants
The evolution of a star of less than ∼6.5M 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 ∼11M experience an additional nuclear burning phase. Carbon burning leads to the formation of a degenerate oxygen-neon (ONe) core. Stars up to ∼8M follow a similar evolutionary path discussed above, but they end their lives as oxygen-neon white dwarfs. In the mass range ∼8-11M , the oxygen-neon core reaches the Chandrasekhar mass, and collapses to a neutron star (NS).
More massive stars than ∼11M 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. 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.32M (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.

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. Janka, 2012;Lai, 2004), can give rise to a natal-kick v k to the star. Neutron stars are expected to receive a kick at birth of about 400 kms −1 (e.g. Cordes et al., 1993;Hobbs et al., 2005;Lyne and Lorimer, 1994), however smaller kick velocities in the range of 50 kms −1 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 and Nelemans, 2015;Repetto et al., 2012;Wong et al., 2014;Zuo, 2015).

Binary evolution
The evolution of a binary can be described by the masses of the stars m 1 and m 2 , the semi-major axis a, and the eccentricity e. 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) where the mass ratio q = m 1 /m 2 . 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.

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(Huang, , 1963.
If none of the wind-matter is accreted, the wind causes the orbit to widen. From angular momentum conservation, it follows as: where m 1 and m 2 are the masses of the stars,ṁ 1 is the mass lost in the wind of the star with mass m 1 (ṁ 1 0), a is the semi-major axis of the orbit, anḋ a wind,no−acc the change in the orbital separation with no wind accretion. Eq. 5 can be rewritten to: where a f and a i are the semi-major axis of the orbit before and after the wind mass loss, and ∆m wind is the amount of matter lost in the wind (∆m wind 0). 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 m 2 accretes at a rate oḟ m 2 = −βṁ 1 . Note that Eq. 7 reduces to Eq. 5 for complete non-conservative mass transfer i.e. β = 0. 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 (Chaty, 2011;Tauris and van den Heuvel, 2006) 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 kms −1 (e.g. Höfner, 2015), which is comparable to the velocity of stars in a binary of a ≈ 10 3 R . Hydrodynamical simulations of such binaries suggest that the wind of the donor star is gravitationally confined to the Roche lobe of the donor star (de Val-Borro et al., 2009;Podsiadlowski, 2007, 2011). The wind can be focused towards the orbital plane and in particular towards the companion star. This scenario (often called wind Rochelobe 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 (Blind et al., 2011;Karovska et al., 2005) 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 (10 −7 − 10 −4 M yr −1 ) for orbital separations approximately larger then 10 6 R (P ≈ 10 6 yr where P is the orbital period).
Lastly, anisotropic mass-loss might occur in fastrotating 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 √ a. 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.

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: anḋ whereȧ gr andė gr 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 (Huang, 1966;Schatzman, 1962;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 lowmass 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.2M , 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).

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 J > 1 3 J b , with J b the orbital angular momentum and J = IΩ, where I 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 m 2 raises tides on a star with mass m 1 , the change of binary parameters due to tidal friction is: whereq ≡ m 2 /m 1 , and Ω b = 2π/P is the mean orbital angular velocity. The star with mass m 1 has an apsidal motion constant k am , gyration radius k, and spin angular frequency Ω. τ TF represents the typical timescale on which significant changes in the orbit take place through tidal evolution. The parameters f n (e 2 ) are polynomial expressions given by (Hut, 1981): f 1 (e 2 ) = 1 + 31 2 e 2 + 255 8 e 4 + 185 16 e 6 + 25 64 e 8 f 2 (e 2 ) = 1 + 15 2 e 2 + 45 8 e 4 + 5 16 e 6 f 3 (e 2 ) = 1 + 15 4 e 2 + 15 8 e 4 + 5 64 e 6 f 4 (e 2 ) = 1 + 3 2 e 2 + 1 8 e 4 f 5 (e 2 ) = 1 + 3e 2 + 3 8 e 4 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 T ( Eq. 10-12) is subject to debate due to quantitative uncertainties in tidal dissipation mechanisms (Meibom and Mathieu, 2005;Willems, 2003;Witte and Savonije, 1999). 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 attributed [2] to turbulent friction in the convective regions of the star (Goldman and Mazeh, 1991;Zahn, 1977Zahn, , 1989. 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(Zahn, , 1977. Despite the uncertainties in tidal dissipation mechanisms, it is generally assumed that circularization and synchronization is achieved before RLOF in a binary.

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. Hjellming and Webbink, 1987;Pols and Marinus, 1994; [2] For alternative mechanisms, see (Goodman and Dickson, 1998;Savonije and Witte, 2002;, 2000). Soberman et al., 1997;Webbink, 1985). 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 CEevolution upon RLOF [3] . 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).

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 [3] Unless q 0.7, such that the orbit and the Roche lobe expand significantly upon mass transfer (e.g. Eqs. [18][19] 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 (Nandez et al., 2015;Passy et al., 2012b;Ricker and Taam, 2012). 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 E gr is the binding energy of the envelope, E orb,i and E orb,f are the orbital energy of the pre-and postmass 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 E gr 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 R is the radius of the donor star, M d,env is the envelope mass of the donor and λ ce depends on the structure of the donor (de Kool et al., 1987;Dewi and Tauris, 2000;Loveridge et al., 2011;Xu and Li, 2010). The parameters λ ce and α are often combined in one parameter αλ ce . According to the alternative γ-formalism (Nelemans et al., 2000), angular momentum is used to expel the envelope of the donor star, according to: where J b,i and J b,f are the orbital angular momentum of the pre-and post-mass transfer binary respectively. The parameters m d and m a represent the mass of the donor and accretor star, respectively, and ∆m d 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 (αλ ce ≈ 0.25, see Camacho et al., 2014;Portegies Zwart, 2013;Toonen and Nelemans, 2013;Zorotovic et al., 2010). However, for the formation of the second WD in double WDs, the orbit only shrinks moderately (αλ ce ≈ 2, see Nelemans et al., 2000Nelemans et al., , 2001van 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(Nelemans et al., , 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 nonconservatively 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.

Stable mass transfer
Whereas the duration of the CE-phase is likely of the order of 10 3 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 i and f denote the pre-and postmass 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 β ≡ −ṁ a /ṁ d of the transferred matter, and the matter that is lost carries specific angular momentum h 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 β = 1 orṁ a =ṁ d . 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 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 cluster [4] . 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.

Supernova explosions in binaries
If the collapsing star is part of a binary or triple, natal kick v k 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;Pijloo et al., 2012;Tauris and Takens, 1998) through: where a i and a f are the semi-major axis of the pre-SN and post-SN orbit, ∆m is the mass lost by the collapsing star, m t,i is the total mass of the system pre-SN, r i is the pre-SN distance between the two stars, v i is the pre-SN relative velocity of the collapsing star relative to the companion, and [4] Stable 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.
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 6.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 v i ·v k 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.

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 m 1 and m 2 , and the mass of the outer star in the outer orbit m 3 ; -the semi-major axis a, the eccentricity e, the argument of pericenter g of both the inner and outer orbits. Parameters for the inner and outer orbit are denoted with a subscript 'in' and 'out', respectively; -the mutual inclination i r between the two orbits. The longitudes of ascending nodes h 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 h in − h out = π for a reference frame with the z-axis aligned along the total angular momentum vector .
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.

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 (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 a out a in < a out a in | crit 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, 2004;Aarseth and Mardling, 2001). 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 (Freire et al., 2011;Iben and Tutukov, 1999;Kiseleva et al., 1994;Portegies Zwart et al., 2011). Regarding wind mass losses from the inner binary exclusively, the fractional mass losses |ṁ|/(m 1 + m 2 ) > |ṁ|/(m 1 + m 2 + m 3 ). Therefore, the fractional orbital increasesȧ in /a in >ȧ out /a out , 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.

Lidov-Kozai mechanism
Secular dynamics can play a mayor role in the evolution of triple systems. The key effect is the Lidov-Kozai mechanism (Kozai, 1962;Lidov, 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).
When the three-body Hamiltonian is expanded to quadrupole order in a in /a out , the timescale for the Lidov-Kozai cycles is (Kinoshita and Nakai, 1999): where P in and P out 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 e max is a function of the initial mu-tual inclination i i as (Innanen et al., 1997): in the test-particle approximation , i.e. nearly circular orbits (e in = 0, e out = 0) with one of the inner two bodies a massless test particle (m 1 m 0 , m 2 ) and the inner argument of pericenter g in = 90 • . 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 a in /a out i.e. the octupole level of approximation, even richer dynamical behaviour is expected than for the quadrupole approximation (e.g. Blaes et al., 2002;Ford et al., 2000;Lithwick and Naoz, 2011;Naoz et al., 2013;Shappee and Thompson, 2013;Teyssandier et al., 2013). The octupole term is nonzero 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 oct measure the importance of the octupole term compared to the quadropole term, and is defined by: Generally, when | oct | 0.01, the eccentric Lidov-Kozai mechanism can be of importance 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).

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 | oct | 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) (Michaely and Perets, 2014;Shappee and Thompson, 2013). 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.

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 m d is the mass of the distorted star that instigates the precession and m a the companion star in the two-body orbit. The distorted star has a classical apsidal motion constant k am , radius R, 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 timescales [5] 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 (Liu et al., 2015a;Naoz et al., 2012;Petrovich, 2015).

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 (Kiseleva et al., 1998;Mazeh and Shaham, 1979). 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 (Kiseleva et al., 1998;Mazeh and Shaham, 1979). 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 tides [6] . The significance of LKCTF is illustrated by Fabrycky and Tremaine (2007), who show that MS binaries with [5] See for example Eq. 8 and 9. The timescales for GW emission are shorter for binaries with small separations and large eccentricities.
[6] A balance is also possible between the eccentricity excitations of the Lidov-Kozai mechanism and other sources of precession, see Sect.2.3.4). orbital periods of 0.1-10d are produced from binaries with much longer periods up to 10 5 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 (Eggleton and Kiseleva-Eggleton, 2001;Fabrycky and Tremaine, 2007;Hamers et al., 2013;Kisseleva-Eggleton and Eggleton, 2010;Mazeh and Shaham, 1979), 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 r p = a(1 − e). The Roche lobe is smaller for stars that rotate supersynchronously 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 (Church et al., 2009;Lajoie and Sills, 2011;Layton et al., 1998;Regös et al., 2005;van der Helm et al., 2016) as well as analytical approaches (Davis et al., 2013;Dosopoulou and Kalogera, 2016a,b;Sepinsky et al., 2007bSepinsky et al., , 2009Sepinsky et al., , 2010. 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 (Dosopoulou and Kalogera, 2016a,b;Sepinsky et al., 2007bSepinsky et al., , 2009Sepinsky et al., , 2010. 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 a out /a in 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 nonconservatively, 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 commonenvelop 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 a out /a in 3, v esc,in > v esc,out unless m 3 f (m 1 + m 2 ). The factor f is of the order of one, e.g. for circular orbits f = 2. In the catalogue of Tokovinin (2014b) it is uncommon that m 3 (m 1 + m 2 ). Out of 199 systems, there are 3 systems with (m 1 + m 2 ) < m 3 < 2(m 1 + m 2 ) and none with m 3 2(m 1 + m 2 ). 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 , where a ring is the radius of the circumbinary ring (Soberman et al., 1997). For example, for a binary with q = 2, 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., , 2015Martin 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 spiralin 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 10 4 R 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. Passy et al., 2012b;Ricker and Taam, 2012), 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 CEmatter 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.

Mass loss from the outer star
In about 20% of the multiples in the Tokovinin catalogue of multiple star systems in the Solar neighbour-hood, 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 (Sabach and Soker, 2015;Tauris and van den Heuvel, 2014). 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 αλ ce 3 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 point [7] of the outer orbit. As a consequence [7] Here 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 mass and angular momentum loss, the outer orbit shrinks with [8] η ≈ 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 R acc,column is the width of the accretion column at the binary location, and R acc,B−H 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 q 3.3. In this case the accretion is in a steady-state and mainly towards the most massive star of the inner binary.

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, 2000;Bond and Livio, 1990;De Marco et al., 2015;Zijlstra, 2015) or multiple companions (e.g. Bear and Soker, 2016;Bond et al., 2002;Exter et al., 2010;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).
[8] There 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 ∆J J b , where ∆J is the amount of angular momentum that is lost from the system, and J b the orbital angular momentum. It does not represent ∆J Ja where Ja is the angular momentum of the accretor star.
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, 1994;Soker et al., 1992). 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. 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 semimajor axis of the outer orbit is affected by the SN in the inner binary as:

Supernova explosions in triples
where m t,i is the total mass of the pre-SN triple, r i and r f are the pre-SN and post-SN distance between the star in the outer orbit and the center of mass of the inner binary, v i is the pre-SN relative velocity of the center of mass of the inner binary relative to the outer star, v sys is the systemic velocity the inner binary due to the SN, and v c ≡ Gm t,i a i the orbital velocity in a circular orbit. Note that in comparison with the circular velocity in binaries (Eq. 22), here m t,i refers to m 1 +m 2 +m 3 and a = a out . 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 6.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 6.1).

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) [9] . 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). 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 .

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, threebody 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, 2013;Portegies Zwart et al., 2009). This is a component library with a homogeneous interface structure based on Python. AMUSE can be downloaded for free at amusecode.org and github.com/amusecode/amuse. 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).

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 dt trial and evolve the stars as single stars for this timestep (step 1). The trial timestep is estimated with: dt trial = min(dt star , dt wind , dt R , f prev dt prev ), (35) where dt star , dt wind , dt R and f prev dt prev 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 dt trial , where the minimum is adopted as a trial timestep in TrES . The timestep dt star 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 dt wind = f wind m/ṁ wind , where f wind = 0.01 anḋ m wind is the wind mass loss rate given by the stellar evolution code. The numerical factor f wind 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 dt R = f R f R · R/Ṙ, where f R and f R are numerical factors. We take f R = 0.005 and whereṘ andṘ prev 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 dt trial < f prev dt prev , where dt prev is the previous successfully accomplished timestep and f prev 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 dt star . 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 dt trial . This process is done iteratively until the conditions are met or until the timestep is sufficiently small, dt min = 10 −3 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 0.9dt trial represents 90% of the previous timestep and dt trial (Ṙ ) is a newly calculated timestep according to Eq. 35 for which the time derivative of the radiuṡ R from the last trial timestep is used. During mass transfer, the timestep is estimated by: whereṁ MT is the mass loss rate from mass transfer (Sect. 3.4.5), and the numerical factors f MT = 0.01 and f MT,prev = 5. If a star starts filling its Roche lobe in a timestep, dt trial = 0.5dt trial . If dt trial < dt trial,MT , 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 dt MT , 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.

Stellar evolution
Single stellar evolution [10] is included through the fast stellar evolution code SeBa (Nelemans et al., 2001;Portegies Zwart and Verbunt, 1996;Toonen and Nelemans, 2013;Toonen et al., 2012). 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, 1971Eggleton, , 1972) that solves the stellar structure equations.

Orbital evolution
TrES solves the orbital evolution through a system of first-order ordinary differential equations (ODE): [10] Note that SeBa is not used to model binary evolution in TrES .
in =ȧ in,GR +ȧ in,TF +ȧ in,wind +ȧ in,MṪ a out =ȧ out,GR +ȧ out,TF +ȧ out,wind + a out,MṪ e in =ė in,3b +ė in,GR +ė in,TḞ e out =ė out,3b +ė out,GR +ė out,TḞ g in =ġ in,3b +ġ in,GR +ġ in,tides +ġ in,rotatė g out =ġ out,3b +ġ out,GR +ġ out,tides + g out,rotatė where θ ≡ cos(i), J b,in and J b,out are the orbital angular momentum of the inner and outer orbit, and Ω 1 , Ω 2 and Ω 3 the spin frequency of the star with mass m 1 , m 2 and m 3 , respectively, and I the moment of inertia of the corresponding star.ẋ represents the time derivative of parameter x. 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 dt trial . Within the ODE solver, parameters that are not given in Eq. 39 (e.g. gyration radius), are assumed to be constant during dt trial . An exception to this is the stellar radius, mass and moment of inertia. Even though dt trial 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Ω I 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 where m c and R c are the mass and radius of the core, k 2 = 0.1 and k 3 = 0.21 (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 m c ≡ 0M , and compact objects for which m − m c ≡ 0M . 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 quadrupoleorder level of approximation assuming J b,out J b,in , 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 .

The secular approach
The components in Eq. 39 (ė 3b ,ġ 3b , andḣ 3b ) that describe the secular three-body dynamics, including the regular Lidov-Kozai cycles, are derived using the orbital-averaging technique (e.g. Luo et al., 2016;Michaely and Perets, 2014;Naoz, 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 a in /a out (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. | oct | < 10 −4 . For example, this includes systems with m 1 ≈ m 2 . 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. Antognini et al., 2014;Antonini and Perets, 2012;Antonini et al., 2014;Bode and Wegg, 2014;Katz and Dong, 2012;Luo et al., 2016;Naoz et al., 2013). 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 (| oct | 0.05).

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 β 1→2 of the wind mass lost from m 1 at a rateṁ 1 can be accreted by m 2 , and β 2→1 for the mass flowing in the other direction. Following Eq. 7, the effect on the inner orbit is then: a in,wind =ȧ wind (ṁ 1 , β 1→2 ) +ȧ wind (ṁ 2 , β 2→1 ). (44) The wind mass lost from the inner binary isṁ in = (1 − β 1→2 )ṁ 1 + (1 − β 2→1 )ṁ 2 , of which the outer star can accrete a fraction β in→3 . We do not allow the inner binary to accrete mass from the outer star. The winds widen the orbit according to: a out,wind =ȧ wind (ṁ in , β in→3 )+ȧ wind,no−acc (ṁ 3 ), (45) see Eq. 5 and Eq. 7.
The wind matter carries away an amount of angular momentum which affects the spin of the star. Under the assumption that the wind mass decouples from the star as a spherical shell: If wind matter is accreted by a star, we assume the accretor star spins up i.e. the stellar spin angular momentum increases with the specific angular momentum of the wind matter. For example for m 1 , the total change in the spin due to winds is:

Tides and precession
Tidal friction is included in TrES as described in Eq. 10-12. The dominant tidal dissipation mechanism is linked with the type of energy transport in the outer zones of the star. We follow Hurley et al. (2002), and distinguish three types: damping in stars with convective envelopes, radiative envelopes (i.e. dynamical tide), and degenerate stars. The quantity k am /τ TF of Eq. 10-12 is given for these three types of stars in their Eq.30, 42 [11] and 47, respectively. We assume that radiative damping takes place in MS stars with M > 1.2M , in helium-MS stars and horizontal branch stars. Excluding compact objects, all other stars are assumed to have convective envelopes. For the mass and radius of the convective part of the stellar envelope, we follow Hurley et al. (2000) (their Sec.7.2) and Hurley et al. (2002) (their Eq.36-40), respectively, with the modification that MS stars to have convective envelopes in the mass range 0.3-1.2M . Regarding the gyration radius k, for stars with convective or radiative envelopes we assume k = k 2 , for compact objects k = k 3 (see Eq. 41).

Stability of mass transfer initiated in the inner
binary When one of the inner stars fills its Roche lobe, we test for the stability of the mass transfer: • Tidal instability; Tidal friction can lead to an instability in the binary system and subsequent orbital decay (see Sect. 2.2.3). The tidal instability takes place in compact binaries with extreme mass ratios. It occurs when there is insufficient angular momentum [11] Note that there is an error in Eq.42 of Hurley et al. (2002). The factor M R 2 /a 5 should be raised to the power 1/2, which means that kam/τTF ∝ R instead of kam/τTF ∝ R 2 .
to keep the star in synchronization i.e. J > 1 3 J b , with J = IΩ. When RLOF occurs due to a tidal instability, we assume that a CE develops around the inner binary. This will lead further orbital decay, and finally either a merger or ejection of the envelope.
• RLOF instability; The stability of the mass transfer depends on the response of the radius and the Roche lobe to the imposed mass loss. In the fundamental work of Hjellming and Webbink (1987), theoretical stability criteria are derived for polytropes. Stability criteria have been improved with the use of more realistic stellar models, see (e.g. Ge et al., 2010Ge et al., , 2015Passy et al., 2012a;Woods et al., 2010Woods et al., , 2012. Our incomplete understanding of the stability of mass transfer leads to differences between synthetic binary populations (Toonen et al., 2014). The response of the Roche lobe is strongly dependent on the envelope of the donor star and the mass ratio of the system [12] . Therefore the stability of mass transfer is often described by a critical mass ratio q crit < q ≡ m d /m a for different types of stars. For unevolved stars with radiative envelopes, mass transfer can proceed in a stable manner for relatively large mass ratios. We assume q crit = 3, unless the star is on the MS for which we take q crit = 1.6 de Mink et al., 2007b). Stars with convective envelopes are typically unstable to mass transfer, unless the donor is considerably less massive than the companion. For giants, we adopt q crit = 0.362 + [3(1 − M c /M )] −1 , where M c is the core mass of the donor star (Hjellming and Webbink, 1987). For naked helium giants, low-mass MS stars (M < 0.7M ), and white dwarfs, we follow Hurley et al. (2002) and adopt q crit = 0.784, q crit = 0.695 and q crit = 0.628, respectively.

Common-envelope evolution in the inner binary
As CE-evolution is a fast, hydrodynamic process, the ODE solver routine is disabled during the modelling of the CE-phase. If the donor star is a star without a clear distinction of the core and the envelope (i.e. MS stars, helium MS stars and remnants), we assume the phase of unstable mass transfer leads to a merger. For other types of donor stars, the treatment of the CE-phase consists of three different models that are based on combinations of the formalisms [12] and to a lesser degree also the accretion efficiency and the corresponding angular momentum loss mode (e.g. Soberman et al., 1997;Toonen et al., 2014). described in Sect. 2.2.5. In model 1 and model 2, the αformalism (Eq. 14) and γ-formalism (Eq. 17) are used to determine the outcome of the CE-phase, respectively. When two giants are involved, the double-CE is applied (Eq. 15). In the standard model, the γformalism is applied unless the CE is triggered by a tidal instability or the binary contains a remnant star. This is based on modelling the evolution of double white dwarfs (Nelemans et al., 2000(Nelemans et al., , 2001Toonen et al., 2012). The standard values of αλ ce and γ are taken to be 2.0 and 1.75 (Nelemans et al., 2000). The companion star in the inner binary is probably not able to accrete from the overflowing material of the CE-phase, because of its relatively long thermal timescale compared to the short timescale on which the CE is expected to place. Therefore, we assume that the CE occurs completely non-conservatively.
The effect of a CE-phase on the outer star of a triple is poorly studied or constrained (Sect. 2.3.8). For stable, hierarchical systems, if the CE-material is energetic enough to escape from the inner binary, the matter is likely energetic enough to escape from the triple as well (Sect.2.3.7). We assume that the escaping CE-matter has expanded and diluted sufficiently to avoid a second CE-phase with the outer star, as we only consider stable triples with a out /a in 3. We allow the matter to escape as a fast wind in a nonconservative way i.e. according to Eq. 45 withṁ 3 = 0 and β in→3 = 0, i.e. a out,wind =ȧ wind,no−acc (ṁ in ).
There may be friction between the CE-matter and the outer star, if the CE-matter is primarily expelled in the orbital plane and the inner and outer orbital planes are parallel. Therefore we multiply Eq. 48 with a factor f fric = min 1, where i crit is a minimum inclination necessary for the friction to take place.

Stable mass transfer in a circular inner binary
We assume mass transfer in the inner binary takes place on either the thermal or nuclear timescale of the donor star. Mass transfer driven by angular momentum loss is currently not implemented in the code. The mass transfer rate is estimated byṁ MT = m/τ MT , where τ MT is the timescale of mass transfer. The thermal timescale of a star is given by Eq. 2, where R and L are given by the stellar evolution code used in TrES , SeBa. The nuclear timescale of a MS or helium-MS star is estimated by Eq. 3. For other stars we take τ nucl = R/Ṙ, whereṘ is the time derivative of the radius, calculated from the current and previous timestep. If the star is shrinking, which is possible for horizontal branch stars or evolved AGB stars, we estimate the nuclear timescale by 10% of the stellar age. Rejuvenation of the accretor star, and the opposite process for the donor star are taken into account by SeBa. Their method is explained in Appendix A.2.1 of Toonen et al. (2012). The orbital evolution of the inner binary is approximated with Eq. 19 where β and η are taken to be constants. If the companion star in the inner binary fills its Roche lobe during the mass transfer phase in response to the accretion, a contact binary is formed. We allow the inner binary to go through a CE-phase as described in Sect. 3.4.4, which likely leads to a merger of the system.
The degree of conservativeness β of the mass transfer is one of the major uncertainties in binary evolutionary calculations. The accretor star is expected to spinup due to the accretion. Even if the companion only accretes a few percent of its own mass, the accretor is spin-up to critical rotation (Packet, 1981). This has been invoked to limit the amount of accretion that can take place. However, as some binaries have managed to experience a phase of fairly conservative mass transfer (e.g. φ Per Pols, 2007), additional mechanisms of angular momentum loss must play a role during mass transfer (de Mink et al., 2007a).
The lack of synchronisation can affect the size of the Roche lobe significantly. For example, for a star that is rotating 100 times faster than synchronization, the Roche lobe is only 5-10% of that of the classical Roche lobe (based on Sepinsky et al., 2007a). For simplicity, we make the common assumption that any circularized system entering RLOF, is and will remain synchronized during the mass transfer phase.
For stable, hierarchical triples systems, the matter lost by the inner binary is likely energetic enough to escape from the triple (Sect. 2.3.7), and we model this as a fast wind i.e. according to Eq. 48. To incorporate the effect of friction between the matter and the outer star, we multiply Eq. 48 with a factor 0 < f crit < 1 (Eq. 49), similar to the case of a CE-phase in the inner binary (Sect. 3.4.4. During mass transfer, the effect of wind mass losses on the inner and outer orbit (Eqs. 44 and 45) are taken into account simultaneously.

Mass transfer initiated by the outer star
Mass transfer from an outer star onto a binary is an intriguing new evolutionary pathway opened up by stellar triples. Even though it is relatively common for triples (about 1%), it is a complex process that has not been studied in much detail. The study of (de Vries et al., 2014) focuses on two triples undergoing mass transfer initiated in the outer star, however, the study is limited in parameter space (as they study two triples), and in time (as the hydrodynamical simulations they perform are expensive). For this reason, we have not implemented the process of outer mass transfer, and currently the code is stopped when the outer star fills its Roche lobe.

Supernova explosions in TrES
During a SN event, the star collapses on a dynamical timescale, for which the secular approach is not valid. The ODE solver routine is therefore disabled, and the orbital evolution due to the SN event is solved for in a separate function as detailed below.
The amount of mass-loss in the SN-event, and the type or remnant that is left behind are determined by the stellar evolution code SeBa. The effect of the SN ejecta on the companion stars (e.g. compositions and velocities) is usually small (e.g. Hirai et al., 2014;Kalogera, 1996;Liu et al., 2015b;Rimoldi et al., 2015), unless the pre-supernova separation between the stars is smaller than a few solar radii. For this reason, we assume the dynamics of the companion stars are not affected by the expanding shell of material, and the companion stars neither accrete nor are stripped of mass.
We make the common assumption that the SN takes place instantaneously. As a result, the positions of the stars just before and after the SN are not changed. As TrES is based on orbit-averaged techniques, we do not follow the position of the stars along the orbit as a function of time. In order to obtain the position at the moment of the SN, we randomly sample the mean anomaly from a uniform distribution. The natal kick is randomly drawn from either of three distributions (Paczynski (1990), Hansen and Phinney (1997) or Hobbs et al. (2005)) in a random direction. Our method simply consists of two coordinate transformations (thus we do not use Eq. 21,33,63,70,73,nor Eq.78 directly). We convert from our standard orbital parameters of i and a, e, g, h for the inner and outer orbit to orbital vectors i.e. eccentricityê and angular momentum vectorĴ b for both orbits. After the mass of the dying star is reduced and the natal kick is added to it, we convert back to the orbital elements. The reason for performing two coordinate transformation, to orbital vectors and back, is that the orbital elements in the code are defined with respect to the 'invariable' plane, i.e. in a frame defined by the total angular momentum. In the case of a SN, however, the total orbital angular momentum vector is not generally conserved, which implies that the coordinate frame changes after the SN. In contrast, the orbital vectors are defined with respect to an arbitrary inertial frame that is not affected by the SN. The post-SN orbital vectors are transformed to the orbital elements in the new 'invariable' plane, i.e. defined with respect to the new total angular momentum vector. An additional advantage of the double coordinate transformation is that the pre-supernova orbit can be circular as well as have an arbitrary eccentricity.
If the post-supernova eccentricity of an orbit is larger than one, the orbit is unbound. We distinguish four situations: • Both the inner as the outer orbit remain bound, and the system remains a triple. The simulation of the evolution of the triple is continued. • When the inner orbit remains bound, and the outer orbit becomes unbound, the outer star and inner binary remain as separated systems. We assume the outer star does not dynamically affect the inner binary. With the default options in TrES the simulation is stopped here unless the user specifies otherwise. • When both the inner as the outer orbit become unbound, the stars evolve further as isolated stars. As in the previous scenario, by default the simulation is stopped unless the user specifies otherwise. • The inner orbit becomes unbound, but at the moment just after the SN the outer star remains bound the inner system. In this case, TrES cannot simulate the evolution of this system further. The evolution of these systems should be followed up with an N-body code.

Examples
In Sect. 2.2 and Sect. 2.3, we discussed several physical processes and how they affect the long-term evolution of inner binaries and triple systems. Here we illustrate those processes by simulating the evolution of a few realistic triple star systems. For example, the evolution of Gliese 667 displays the Lidov-Kozai cycles, and the evolution of Eta Carinae illustrates the effect of precession and stellar winds. The evolutionary pathways of the triple systems are simulated with the new triple code TrES , such that the examples below also demonstrate the capabilities of TrES .

Gliese 667
Gliese 667 is a nearby triple system in the constellation of Scorpius. The orbital parameters of the system are described in Tbl. 3 based on (Tokovinin, 2008). The outer star is in an orbit of a out > 230AU, but for simplicity, we will assume a out = 250AU in the following. The outer star is also a planetary host-star; up to five planets have been claimed, of which two have been confirmed so far (Feroz and Hobson, 2014). The orbit  -Escudé et al., 2012). In the following, we will neglect the dynamical effect of the presence of planets on the evolution of the triple. Gliese 667 is a prime example of a triple system undergoing  show the evolution of the inner eccentricity and mutual inclination for the first 3Myr after the birth of the system. under the assumption of e out = 0.5, i = 90 • , g in = 0.1 and g out = 0.5. For different values for the outer eccentricity, arguments of pericenters, and mutual inclination the general behaviour of Figs. 3-5 remains the same, but the timescale and amplitude of the Lidov-Kozai cycles varies (to the point where the cycles are not notable). Figs.3 and 4 show the cyclic behaviour of eccentricity and inclination in Gliese 667. When the eccentricity is at its maximum, the inclination between the orbits is minimal. The timescale of the oscillations is a few 0.1Myr, which is consistent with the order of magnitude approximation of 0.4Myr of Eq. 24. The octupole parameter oct < 0.001, which indicates that the eccentric Lidov-Kozai mechanism is not of much importance here.
For the same timescale as Fig. 3 and 4, Fig. 5 shows the evolution of the inner-orbital semi-major axis a in . The change in the inner semi-major axis of Gliese 667 is negligibly small, however, the figure illustrates the effect of Lidov-Kozai cycles with tidal friction or LKCTF. When the inner eccentricity is at its maximum, and the inner stars are at their closest approach during pericenter passage, the inner semi-major axis decreases due to tidal forces. In this way, LKCTF could lead to RLOF in or a merger of the inner system.
Regarding the long-term evolution of this system, it is analogous to that discussed for the first 3 Gyr, After 10 Gyr (approximately the age of the Galactic thin disk, e.g. del Peloso et al., 2005;Oswalt et al., 1996;Salaris, 2009), the system is still detached. The orbital separation has decreased by only ∼7km. The stellar masses are sufficiently low, that the stars do not evolve off the MS within 10 Gyr, and as such do not experience a significant growth in radius that could lead to RLOF. Even taking into account the low metallicity of Gliese A (Cayrel de Strobel et al., 2001), and the corresponding speed-up of the evolutionary timescales, a 0.73M star is not massive enough to evolve of the MS within 10 Gyr. Furthermore, the stars are not massive enough to lose a considerable amount of matter in

Eta Carinae
Eta Carinae is a binary system with two massive stars (m 1 ∼ 90M and m 2 ∼ 30M ) in a highly eccentric orbit (e = 0.9) with a period of 5.5yr (Damineli et al., 1997). Both stars are expected to explode as supernovae at the end of their stellar lives. Eta Carinae is infamous for its 'Great Eruption'. From 1837 to 1857, it brightened considerably, and in 1843 it even became the second brightest star in the sky (de Vaucouleurs and Eggen, 1952). The system is surrounded by the Homunculus Nebulae, that was formed during the Great Eruption, and heavily obscures the binary stars (Humphreys and Davidson, 1999). The kinetic energy of the Homunculus Nebulae is large i.e. 10 49.7 erg (Smith et al., 2003) and comes close to that of normal supernovae. However, as both stars have survived the Great Eruption, Eta Carinae is often referred to as a 'supernova imposter'. The cause of the Great Eruption remains unexplained. A massive outflow, as during the Great Eruption, can be driven by a strong interaction between two stars (e.g. Harpaz and Soker, 2009;Smith, 2011) or a merger of two stars (e.g. Soker and Tylenda, 2003). Recently, Portegies Zwart and van den Heuvel (2016) tested the hypothesis that the Eta Carinae system is formed from the merger of a massive inner binary of a triple system. According to their model, the merger was triggered by the gravitational interaction with a massive third companion star, which is the current ∼30M companion star in Eta Carinae. Here, we simulate the evolution of their favourite model with the initial conditions as given by Tbl. 3. Furthermore, Portegies Zwart and van den Heuvel (2016) assume that the argument of periastron does not affect the tidal evolution, and therefore we arbitrarily set g in = 0.1 and g out = 0.5. During the early evolution of the triple, the system experiences Lidov-Kozai cycles with a timescale of a few kyr, see Fig. 6. The octupole parameter oct = 0.068 > 0.01, which indicates that the system is in the eccentric Lidov-Kozai regime with a timescale of order t oct ∼ t kozai / oct ∼few tens of kyr.
The evolution of the semi-major axis shows two characteristics in Fig. 7. Firstly, as for Gliese 667, the system is affected by LKCTF i.e. the semi-major axis shrinks periodically, due to strong tides at pericenter when the inner eccentricity is at its maximum. Secondly, as the primary star is very massive, strong winds remove large amounts of mass while the star is still on the MS (Sect. 2.1.3). The dynamical effect of such a fast wind is that the inner and outer orbits expand (Sect. 2.2.1). Initially the inner orbit expands faster than the outer orbit, as expected for a strong wind from the inner orbit (Sect. 2.3.1). However, due to the combination of stellar winds with LKCTF for the Eta Carinae progenitor, its outer orbit expands faster, and the triple becomes more dynamically stable.
On a longer timescale, the triple moves from the eccentric Lidov-Kozai regime to the regular regime (| oct | 0.01) as the inner binary loses matter and angular momentum in the stellar winds (Fig. 8). After 3Myr, the octupole parameter has decreased from oct = 0.068 initially, to oct = 0.001. The long-term evolution of the progenitor candidate of Eta Carinae shows another interesting feature in Fig. 9 and 10, i.e. the Lidov-Kozai cycles are quenched. Here the precession due to the distortion and rotation of the stars dominates over the precession caused by the Lidov-Kozai mechanism. As a result, the amplitudes of the cycles in inner eccentricity and mutual inclination are reduced. After approximately 1.5Myr, the evolution of the system is completely dominated by tides, i.e. the system circularizes and the inner semimajor decreases accordingly (Fig. 11). After circularization of the inner binary has been achieved ( 2Gyr), the inner semi-major axis increases again due to the stellar winds from the inner binary. The evolution of the Eta Carinae progenitor (Fig. 6-11) illustrates that both three-body dynamics and stellar evolution matter, and neither can be neglected.
After about 3Myr, the primary star fills its Roche lobe and initiates a mass transfer phase (Fig. 8). The inner semi-major axis is about 0.5AU (109R ), and the massive primary has increased in size to 51R . Even though, the inner semi-major axis (and Roche lobe) of the primary are ∼ 10% smaller around 2Gyr, there is no RLOF yet as the radius of the primary star is ∼50% smaller. At RLOF, the masses of the inner stars have reduced from the initial 110M and 30M to 75M and 29M , and both stars are still on the MS. The mass transfer phase proceeds in an unstable manner (Sect. 3.4.3), and a common-envelope develops that leads to a merger of the inner stars (Sect. 3.4.4). A new star is formed, that is still on the MS, with a mass of 104M . We assume this merger proceeds conservatively, and therefore the outer orbit is not affected, such that a out ≈ 32AU, e out ≈ 0.2 and P ≈ 15.7yr. Prior to the merger, the outer semi-major axis has increased from the initial value of 25 AU to 32AU due to the stellar winds.
The resulting binary is similar to the current Eta Carinae system in mass and orbital period. It is not an exact match, as the evolution of this specific triple is shown for illustrative purposes, and has not been fitted to match the currently observed system. A progenitor study of Eta Carinae to improve the match is beyond the scope of this paper.
We note that the current eccentricity of our remaining binary is low i.e. e ≈ 0.2 compared to the observed e = 0.9. In our simulations, the post-merger eccentricity is equal to the pre-merger outer eccentricity. During the evolution of the Eta Carinae progenitor, the outer eccentricity has remained roughly equal to its initial value of e out = 0.2. The outer eccentricity is not affected strongly by stellar evolution or Lidov-Kozai cycles. If we study the evolution of an alternative progenitor similar to the favourite model of Portegies Zwart and van den Heuvel (2016), but with e out = 0.9, the system is dynamically unstable at birth. In order for the triple to be dynamically stable e out 0.81 for the standard i = 60 • , or up to e out 0.84 for i = 90 • . The evolutions of the dynamically stable systems with e out 0.7 show similar behaviour as our initial system (Tbl. 3), and the merger leads to a similar binary as in the case of our initial system. For dynamically stable systems with higher outer eccentricities, the merger time decreases strongly, and the inner system does not reach circularization before the merger takes place. The merger product is more massive, as less mass is lost in stellar winds. In the simulation of (Portegies Zwart and van den Heuvel, 2016), the outer orbit has an eccentricity e out = 0.2 initially, but becomes highly eccentric in the merger phase due to asymmetric mass loss.

MIEK-mechanism
In this section, we illustrate the dynamical effect of mass loss on a triple system from the point of view of transitions in dynamical regimes, e.g. the regime without Lidov-Kozai cycles, with regular or with eccentric Lidov-Kozai behaviour. Here, we focus on the transition from the regular Lidov-Kozai regime to the eccentric regime, i.e. where the octupole term is significant. This transition has been labelled 'mass-loss induced eccentric Kozai' or MIEK (Sect. 2.3.3). The canonical example of MIEK-evolution is a triple with the initial conditions as given by Tbl. 3 (Michaely and Perets, 2014;Shappee and Thompson, 2013). To reproduce the experiment of Shappee and Thompson (2013), we simulate the evolution of this triple with  Figure 11 Inner semi-major axis evolution The evolution of the inner semimajor-axis on the same timescale as Fig. 9 and 10. In the first Myr, the evolution of the system is dominated by the Lidov-Kozai mechanism and the inner semimajor-axis remains more or less constant. In the following Myr, the system circularizes and the semimajor-axis decreases by a factor 2. After synchronisation and circularisation has been reached, the inner semimajor-axis increases due to the ejection of stellar winds. TrES including three-body dynamics and wind mass losses, however, without stellar evolution in radius, luminosity, or stellar core mass etc. Starting from the birth of the triple system, its orbit is susceptible to Lidov-Kozai cycles (Fig. 13 and 14). The timescale of the cycles is approximately 0.1Myr.
The cycles are in the regular regime, i.e. oct = 0.002. As time passes, the stars evolve. The primary star evolves off the MS at 49Myr, and after 55Myr it reaches the AGB with a mass of 6.9M (Fig. 12). Subsequently, it quickly loses a few solar masses in stellar winds, before it becomes an oxygen-neon white dwarf of 1.3M at 56Myr. The outer orbit widens to about a out ∼ 350AU due to the wind mass losses. Figure 12 Radius and mass evolution The evolution of the radii (dashed line) and mass (dash-dotted line) as a function of time for the stars in a triple that transitions from a region with regular to eccentric Lidov-Kozai behaviour i.e. MIEK. The initial conditions of the triple are given in Tbl. 3, based on Shappee and Thompson (2013) and Michaely and Perets (2014). The primary star is shown in blue, the secondary in green and the tertiary in red. The Roche lobe of the primary and secondary are overplotted (blue solid line and green solid line respectively). The Roche lobe of the tertiary is 6500-8000R . After about 55.5Myr, the primary star fills its Roche lobe. The wind mass loss allows the triple to transition to the eccentric Lidov-Kozai regime at about 56Myr, i.e. oct = 0.045 at this time. The system is driven into extremely high eccentricities, and also the amplitude of the Lidov-Kozai cycle in inclination increases. The evolution of the system as shown in Figs. 13 and 14 is qualitatively similar to that found by Shappee and Thompson (2013) based on N-body calculations and Michaely and Perets (2014) based on the secular approach. In these studies, stellar winds are implemented ad-hoc with a constant mass loss rate for a fixed time interval starting at a fixed time. Moreover, the system is followed for multiple Myrs after the mass loss event in both papers, such that the inclination rises above 90 • , and the inner and outer orbit become retrograde to each other. In our case the simulation is stopped before such a flip in inclination develops, as RLOF is initiated in the inner binary when the inner eccentricity is high. Figure 13 Inner eccentricity evolution The evolution of the inner eccentricity e in as a function of time for a triple that transitions from a region with regular to eccentric Lidov-Kozai behaviour i.e. MIEK. The initial conditions of the triple are given in Tbl. 3, based on Shappee and Thompson (2013) and Michaely and Perets (2014). For this figure, the stars are not allowed to evolve in TrES , except for wind mass losses. If stellar evolution is taken into account fully, RLOF initiates at 55.5Myr, before the transition to MIEK can develop. However, if we fully include stellar evolution, as in the standard version of TrES , the triple is not driven into the octupole regime. On the AGB, the radius of a 7M -star can reach values as large as ∼1000R (Fig. 12), and therefore RLOF initiates before the MIEK-mechanism takes place. Even if the inner binary would be an isolated binary, RLOF would occur for initial separations of a < 15AU. For triples, RLOF can occur for larger initial (inner) separations, as the Lidov-Kozai cycles can drive the inner eccentricity to higher values. For wider inner binaries i.e. a in > 16AU, the MIEK-mechanism does not occur either, as the triple is dynamically unstable. This example indicates that the parameter space for the MIEK-mechanism to occur is smaller than previously thought, and so it may occur less frequently. Moreover, this example demonstrates the importance of taking into account stellar evolution when studying the evolution of triples.
For the canonical triple with a in = 10AU and a out = 250AU, RLOF occurs at 55.5Myr, just a few 0.1Myr after the primary star arrives on the AGB. In that time, the radius of the primary increased by a factor ∼3, and tides can no longer be neglected. The tidal forces act to circularize and synchronize the inner system, such that e in = 0 at RLOF. The eccentric Kozai-mechanism does not play a role at this point, i.e. oct = 0.0008. The mass of the core has not had enough time to grow to the same size as in the example without RLOF, i.e. the core mass is 1.25M instead of 1.3M . Stellar winds have reduced the mass of the primary star to 6.8M . As the primary has a convective envelope and is more massive than the secondary, a CE-phase develops. We envision three scenarios based on the different models for CE-evolution (Sect. 2.2.5 and 3.4.4). First, the CE-phase leads to a merger of the inner binary, when the inner orbit shrinks strongly, as for the αmodel of CE-evolution with αλ ce = 0.25 (Sect. 2.2.5). Second, the CE-phase leads to strong shrinkage of the orbit, but not enough for the inner stars to merge. In this scenario, the envelope of the donor star is completely removed from the system, and the outer orbit widens to about 350AU, under the assumption that the mass removal affects the outer orbit as a fast wind. Assuming αλ ce = 2 (Eq. 14, Sect. 2.2.5), a in ∼ 0.33AU and oct = 0.0009 after the CE-phase. In this scenario, the triples does not enter the octupole regime, and the MIEK-mechanism does not manifest. Lastly, the CEphase does not lead to a strong shrinkage of the inner orbit, as for the γ-model of CE-evolution with γ = 1.75 (Eq. 17, Sect. 2.2.5). The inner semimajor-axis even increases from 6.0 to 7.3AU. In this scenario, oct = 0.02, such that the perturbations from the octupole level become significant. In this last scenario, the triple undergoes the MIEK mechanism, despite and because of the mass transfer phase.

Discussion and conclusion
In this paper, we discuss the principle complexities of the evolution of hierarchical triple star systems. Hierarchical triples are fairly common and potentially longlived, which allows for their evolution to be affected by (secular) three-body dynamics, stellar evolution and their mutual influences. We present an overview of single star evolution and binary evolution with a focus on those aspects that are relevant for triple evolution. Subsequently, we describe the processes that are unique to systems with multiplicities of higher order than for binaries.
In some cases, the evolution of a hierarchical triple can be adequately described by the evolution of the inner and outer binary separately. In other cases, the presence of the outer star significantly alters the evolution of the inner binary. Several examples of the latter are given in detail. These examples also show the richness of the regime in which both three-body dynamics and stellar evolution play a role simultaneously. Moreover, the examples demonstrate the importance of coupling three-body dynamics with stellar evolution.
Additionally, we present heuristic recipes for the principle processes of triple evolution. These descriptions are incorporated in a public source code TrES for simulating the evolution of hierarchical, coeval, dynamically stable stellar triples. We discuss the underlying (sometimes simplifying) assumptions of the heuristic recipes. Some recipes are exact or adequate (e.g. gravitational wave emission, wind mass loss or Lidov-Kozai cycles), and others are admittedly crude (e.g. mass transfer). The recipes are based on simple assumptions and should be seen as a starting point for discussion and further study. When more sophisticated models become available of processes that influence triple evolution, these can be included in TrES , and subsequently the effect on the triple populations can be studied. For now, the accuracy levels of the heuristic recipes are sufficient to initiate the systematic exploration of triple evolution (e.g. populations, evolutionary pathways), while taking into account three-body dynamics and stellar evolution consistently. We note that simulating through a phase of stable mass transfer in an eccentric inner orbit is currently beyond the scope of the project. However, appropriate methodology for eccentric mass transfer (e.g. Dosopoulou and Kalogera, 2016b;Sepinsky et al., 2007bSepinsky et al., , 2009) has been developed that we aim to implement at a later stage.
The triple evolution code TrES is based on the secular approach to solve for the dynamics of the triple system. It has been shown that this approach is in good agreement with N-body simulations of systems in which the secular approximations are valid (Hamers et al., 2013;Michaely and Perets, 2014;Naoz et al., 2013). The advantage of the secular approach is that the computational time is orders of magnitudes shorter than for an N-body simulation. The secular approach, however, is not valid when the evolutionary processes occur on timescales shorter than the dynamical timescale of the system. In these cases, we either stop the simulation (e.g. during a dynamical instability) or simulate the process as an instantaneous event (such as a common-envelope phase). Lastly, the secular approximation becomes inaccurate when the triple hierarchy is weaker (e.g. Antognini et al., 2014;Antonini and Perets, 2012;Bode and Wegg, 2014;Katz and Dong, 2012;Luo et al., 2016). In this case, the timescale of the perturbation from the outer star onto the inner binary during its periastron passage, is comparable to the dynamical timescale of the inner binary. This can result in extremely high eccentricities and collisions between the stars in the inner binary. With the secular approach, as in TrES , these occurrences are probably underestimated in systems with moderate hierarchies (see also Naoz et al., 2016).
TrES is written in the Astrophysics Multipurpose Software Environment, or AMUSE (Portegies Zwart, 2013;Portegies Zwart et al., 2009), which is based on Python. AMUSE including TrES can be downloaded for free at amusecode.org and github.com/amusecode/amuse. Due to the nature of AMUSE, the triple code can be easily extended to include a detailed stellar evolution code or a direct N-body code. Regarding the latter, this is interesting in the context of triples with moderate hierarchies where the orbit-averaged technique breaks down (as discussed above). Furthermore, it is relevant for triples that become dynamically unstable during and as a consequence of their evolution. For example, Perets and Kratter (2012) show triples that become dynamically unstable due to their internal wind mass losses, are responsible for the majority of stellar collisions in the Galactic field. Consequently, the majority of stellar collisions do not take place between two MS stars, but involve an evolved star of giant-dimensions. Another interesting prospect is the inclusion of triples in simulations of cluster evolution, where triples are often not taken into account e.g. in the initial population, through dynamical formation nor a consistent treatment of the evolution of triple star systems. However, dynamical encounters involving triples are common, reaching or even exceeding the encounter rate involving solely single or binary stars, in particular in low-to moderate-density star clusters (Leigh and Sills, 2011;Leigh and Geller, 2013). Therefore, the evolution of triples might not only be important for the formation and destruction of compact or exotic binaries, but also for the dynamical evolution of clusters in general.

Competing interests
The authors declare that they have no competing interests.
Author's contributions ST wrote the draft of this paper and the public source code TrES . AH derived the equations for the orbital evolution of a triple during a supernova explosion as given in Appendix 6.1. AH also composed the ODE solver routine from updating the routine presented in Hamers et al. (2013). SPZ assisted with the construction of the heuristic recipes for triple evolution. All of the authors contributed corrections and improvements on the draft of the manuscript. All authors read and approved the final manuscript.