Online material
Appendix A: Binary physics in MESA
In this appendix we will list the physics used in the binary evolution module of MESA. For the exact software implementation we refer the reader to the opensource code, and future instrument papers.
Appendix A.1: Evolution of the orbital parameters
The orbital angular momentum of a binary system is: (A.1)where a is the binary separation, and e is the eccentricity of the orbit. M_{d} and M_{a} are respectively the mass of the donor and the accretor, and G is the gravitational constant. The evolution of the binary separation is obtained from the time derivative of Eq. (A.1): (A.2)In MESA, the mass loss and accretion rates, the change in orbital angular momentum, and the change in eccentricity are computed, and then used to update the orbital separation.
The change in mass of both components depends on the mass lost in stellar winds, the mass lost through RLOF, and the fraction of mass that is accreted by the companion. For the donor star the net mass loss is: (A.3)while for the accretor: (A.4)where ϵ_{wind} is the fraction of the wind mass loss that is accreted by the companion and ϵ_{rlof} is the fraction of the RLOF mass loss accreted by the companion. The amount of mass lost to infinity is then: (A.5)The change in orbital angular momentum has three main components: mass lost to infinity in stellar winds, mass lost to infinity during RLOF and the change in angular momentum due to resonances between the binary and a potential CB disk: (A.6)In the wide binary models considered here, the strongest contributions to will come from mass loss during RLOF and the interaction between the binary and a CB disk.
The evolution of the eccentricity is governed by the tidal interactions between the two stars, the effect of phasedependent mass loss/transfer which can pump or decrease the eccentricity based on the mass ratio and the interaction between the binary and a CB disk: (A.7)The relation between the orbital period and the separation is: (A.8)All physical mechanisms relevant to the terms given in Eqs. (A.3)–(A.7) are explained in the following subsections.
Appendix A.2: Rochelobe overflow
The instantaneous Rochelobe radius of the donor star (R_{L,d}) at a specific moment in the orbit is approximated by adapting the fitting formula for the effective Rochelobe radius of Eggleton (1983) by allowing the binary separation to depend on the orbital phase (θ): (A.9)where a is the separation, e the orbital eccentricity. The mass ratios q_{d} and q_{a} are defined as: (A.10)The Rochelobe radius of the accretor can be obtained by inverting the mass ratio in Eq. (A.9), thus replacing q_{d} with q_{a}.
Massloss is implemented in MESA following the prescription of Ritter (1988) and Kolb & Ritter (1990). This formalism differs between massloss from the optically thin region (R_{d} ≤ R_{L,d}), and that from the optically thick region (R_{d}>R_{L,d}). Massloss from the former region is given by: (A.11)where R_{d} is the donor radius and H_{P,L1} is the pressure scale height at the inner Lagrange point (L1), which can be linked to the pressure scale height at the photosphere of the donor H_{P,ph} as: (A.12)Here P_{ph} is the pressure at the donor photosphere and ρ_{ph} is the photospheric density of the donor. γ(q) is a factor that depends on the mass ratio as: (A.13)Ṁ_{0} is the mass loss rate when the donor star just fills its Rochelobe: (A.14)μ_{ph,d} and ρ_{ph,d} are the mean molecular weight and density at the donor photosphere, T_{eff,d} is the effective temperature of the donor and ℛ is the ideal gas constant. F(q) is a function of the mass ratio depending on the Roche geometry, given by (see also Meyer & MeyerHofmeister 1983): (A.15)where g(q) is: (A.16)An approximation of the distance between the centre of mass of the donor and the L1 point in units of the binary separation (x_{L}) is given by Frank et al. (2002): (A.17)In the case that the stellar radius significantly overfills its Rochelobe, the Rochelobe will lie in the optically thick region of the donor. The mass loss is then calculated via: (A.18)where Ṁ_{0} is given in Eq. (A.14), P_{L1} and P_{ph} are respectively the pressure at L1 and at the stellar photosphere. Both the temperature T and the mean molecular weight μ depend on the pressure. F_{3} is a function that depends on the adiabatic exponent Γ_{1}, and is given by: (A.19)The mass loss equations depend on the eccentric anomaly through the Rochelobe radius. The time steps that MESA uses in the evolution are much larger than one orbital period, thus to obtain the average mass loss rate to be used in a time step, the above equations are integrated over the orbit.
RLOF is not necessarily conservative. To allow for mass to be lost from the system the formalism of Tauris & van den Heuvel (2006) and Soberman et al. (1997) is used. This system describes three fractions (α,β,δ) to lose mass from the system:

α: mass lost from the vicinity of the donor as a fast wind (Jeans mode). This is modelled as a spherically symmetric outflow from the donor star in the form of a fast wind. The mass lost in this way carries the specific angular momentum of the donor star.

β: mass lost from the vicinity of the accretor as a fast wind (Isotropic reemission). A flow in which matter is transported from the donor to the vicinity of the accretor, where it is ejected as a fast isotropic wind. Mass lost in this way carries the specific angular momentum of the accretor.

δ: mass lost from a CB coplanar toroid. The radius of the coplanar toroid is determined by γ as R_{toroid} = γ^{2}a.
The accretion efficiency of RLOF is then given by: (A.20)When the mode of mass loss is known, the change in angular momentum can be calculated. Mass accreted by the companion will not change the total angular momentum of the system, thus only fractions α,β and δ will have an influence on . The total effect of mass loss through RLOF on the change in angular momentum is given by Tauris & van den Heuvel (2006): (A.21)
Appendix A.3: Wind mass loss
The amount of mass lost through stellar winds is determined by any of the wind loss prescriptions defined in the stellar part of MESA. The binary module offers the possibility of boosting this mass loss using the CRAP mechanism (Tout & Eggleton 1988), also known as tidally enhanced mass loss. In this model it is assumed that tidal interactions or magnetic activity are responsible for the enhancement of the wind. The enhancement is expected to have a similar dependence on radius over Rochelobe radius as the torque in a tidal friction model, and there is an expected saturation when corotation is reached (in their model when R = R_{L}/ 2). (A.22)The factor B_{wind} is estimated by Tout & Eggleton (1988) to be of the order B_{wind} ≈ 10^{4}, but can vary significantly depending on which system needs to be explained.
Part of the mass lost due to stellar winds is accreted by the companion. In the case of fast winds (v_{wind} ≫ v_{orb}) the accretion fraction is given by the BondiHoyle mechanism (Hurley et al. 2002). Here we give the equations for the donor, those of the accretor can be obtained by switching the d and a subscripts: (A.23)The velocities are: The wind velocity used here is set proportional to the escape velocity from the stellar surface. Based on observed wind velocities in cool supergiants (v_{wind} = 5−35 km s^{1}) (Kučinskas 1998), β_{W} is taken to be 1/8. The free parameter α_{BH} is set to 3/2 based on Boffin & Jorissen (1988).
The angularmomentum loss due to windmass loss, assuming a sphericalsymmetric wind is: (A.27)For the change in total angular momentum of the binary due to stellar winds, the contribution of the accretor need to be added and .
Appendix A.4: phasedependent mass loss
When mass loss, weather due to stellar winds, RLOF or other reasons, is not constant during the binary orbit, it will have an effect on the eccentricity of the system. There are multiple causes of the periodicity of mass loss and accretion. For example, the mass loss might be caused by stellar pulsations, which can transfer more mass at the periastra that coincide with a maximum stellar radius. The methods implemented in MESA are phasedependent windmass loss through tidal interactions, and phasedependent RLOF on eccentric orbits. The latter will boost mass loss near the periastron passage while during apastron the star is completely inside its Roche lobe. The effect of phasedependent mass loss on the orbital eccentricity was studied by Soker (2000) based on the theoretical work of Eggleton (2006, Ch. 6.5).
In calculating the effect on the eccentricity we have to distinguish between mass lost to infinity, and mass that is accreted by the companion. Assuming isotropic mass loss, the change in eccentricity for mass lost from the system is given by: (A.28)where θ is the true anomaly, and Ṁ_{∞}(θ) is the mass lost from the system at a specific phase θ during the orbit. From this equation one can see that constant mass loss will not change the eccentricity, if however the fraction of mass lost near periastron is significantly larger than during the remainder of the orbit, Eq. (A.28) predicts a positive ė.
The effect on ė of mass accreted by the companion, again under the assumption of isotropic massloss, is given by: (A.29)where Ṁ_{acc}(θ) is the mass accreted by the companion at a specific phase θ during the orbit. As mass transfer is expected near periastron, the eccentricity will increase if M_{d} < M_{a}, which is required in the case of stable mass transfer.
To obtain the total change in eccentricity in one orbit, Eq. (A.28) and (A.29) are integrated over the orbit: (A.30)
Appendix A.5: Tidal forces
The gravitational forces acting on the components in binary systems induce a deformation of their structure, and creates tidal bulges. These bulges are lagging behind, thus the gravitational attraction generates a torque on those bulges, which forces the synchronisation and circularisation of the stellar and orbital rotation. The orbital parameters change because the star’s rotational energy is dissipated into heat. There are two cases one has to consider, stars with convective envelopes and stars with radiative envelopes. In the former, the kinetic energy of tidallyinduced largescale currents is dissipated into heat by viscous friction of the convective environment. In radiative stars it is mainly radiative damping of gravity modes that functions as dissipation process. See Zahn (2005) for a review on tidal dissipation.
We use the formalism developed by Hut (1981) to calculate the effect of tides on circularisation and synchronisation. This formalism, developed based on the weakfriction model, results in (where we again give only the equations relevant for the donor star, and those of the accretor are obtained by switching the d and a subscripts): (A.31)for the circularisation, while the synchronisation is given by: (A.32)Here f_{2−5} are polynomials in e: and r_{g} is the radius of giration. Furthermore k/T depends on the star being convective or radiative. In the case of a convective envelope, k/T is given by Hurley et al. (2002, Eqs. (30)–(33)), based on Rasio et al. (1996): (A.37)where τ_{conv} is the eddy turnover timescale (the time scale on which the largest convective cells turn over): (A.38)and f_{conv} is a numerical factor depending on the tidal pumping timescale P_{tid}:
M_{env} and R_{env} are respectively the mass in the outer stellar convective zone, and the stellar radius at the base of that zone. L is the luminosity.
In the case of a radiative envelope, the damping is caused by a range of oscillations that are driven by the tidal field. k/T is then (Hurley et al. 2002, Eqs. (42)–(43) based on Zahn 1975, 1977): (A.41)where E_{2} is a second order tidal coefficient which can be fitted with: (A.42)The total effect on the change in orbital eccentricity is then obtained by combining the tidal forces on both components: (A.43)
Appendix A.6: Circumbinary disks
Circumbinary disks can form around binaries during the RLOF phase, if part of the mass can leave the system through the outer Lagrange points and form a Keplerian disk around the binary. The CB disk – binary resonant and nonresonant interactions have been described by Goldreich & Tremaine (1979) and Artymowicz & Lubow (1994), by using a linearperturbation theory. There are two main assumptions in these models; the first is that the disk is thin (0.01 < H/R < 0.1, where H and R are respectively the thickness and the halfangularmomentum radius of the disk). The second assumption is that the nonaxisymmetric potential perturbations are small around the average binary potential.
The effect of the disk – binary resonances on the orbital parameters has been the subject of many studies. In MESA we follow the approach of Artymowicz & Lubow (1994) and Lubow & Artymowicz (1996). The model of Lubow & Artymowicz (1996) for small and moderate eccentricities (e < 0.2) is based on the result of SPH simulations, which show that for a disk in which the viscosity is independent of the density, the torque is independent of resonance strength and width. The variation in the orbital separation due to the inner and outer Lindblad resonances is (Lubow & Artymowicz 1996): (A.44)where l and m are integer numbers indicating which resonance has the strongest contribution. J_{D} and J_{B} are respectively the angular momentum in the disk and in the binary, and τ_{v} is the viscous evolution time scale, which is the time scale on which matter diffuses through the disk under the effect of the viscous torque. It is given by: (A.45)where α_{D} is the viscosity parameter of the disk and Ω_{b} is the orbital angular frequency. With a disk viscosity of α_{D} = 0.1, the viscous time scale is typically of the order of 10^{5} yr.
To determine the angular momentum of a Keplerian disk, one need to know the surface mass distribution σ in the disk. J_{D} is then given by: (A.46)where the Keplerian velocity of an element in the disk at distance r from the centre of mass, when neglecting the mass of the disk is: . R_{in} and R_{out} are the inner and outer boundaries of the disk. We assume the surface distribution in the disk to depend on the radius to the power of a free parameter δ: (A.47)The distribution constant D_{c} depends on the total disk mass and δ. With this distribution the disk angular momentum can be calculated: (A.48)which depending on parameter δ has the following solutions: (A.49)As shown by the SPH simulations of Lubow & Artymowicz (1996), only the inner part of the disk plays a role in the disk – binary interactions. Thus in the previous equations the angular momentum is only calculated from the region between R_{in} and six times the binary separation (R_{out} = 6a).
To calculate the distribution constant in Eq. (A.47), one has to know the total mass in the disk (M_{disk}). The total disk mass can then be related to the distribution constant as: (A.50)which depending on δ has two solutions: (A.51)The main difference between the implementation here and that of Dermine et al. (2013) is that the disksurfacedistribution parameter δ in our model is consistent throughout the calculation of the disk angular momentum and the distribution constant D_{c}, while Dermine et al. (2013) used δ = 2.5 in the former, while δ = 1 was used in the latter.
We consider two processes to determine the inner radius of the CB disk. The interaction between the binary and the disk will clear the inner part of the disk. SPH simulations performed by Artymowicz & Lubow (1994) show that the inner radius depends on the Reynolds number of the disk gas (ℜ = (H/R)^{2}α^{1}) and the eccentricity. A fitting formula of their results was derived by Dermine et al. (2013): (A.52)The second limiting factor on the inner radius is the dust condensation temperature. Based on Dullemond & Monnier (2010, Eqs. (1)–(12)), and adding an offset for the binary separation we obtain: (A.53)where L_{d} and L_{a} are respectively the luminosity of the donor and accretor, and σ_{bol} is the StefanBoltzmann constant. T_{cond} is the dust condensation temperature which we take at 1500 K. The last factor in Eq. (A.53), is the average binary separation weighted by the luminosity of both components. The actual inner radius of the disk is the maximum of that determined from SPH simulations and the dust condensation radius: (A.54)There are no direct observations of the outer radius of a CB disk around a postAGB or postRGB binary. PostAGB and postRGB disks are likely similar to protoplanetary disks, thus it can be assumed that the surface mass distribution of a CB disk does not follow the same behaviour along the whole radius. The inner part of the disk follows σ(r) ~ r^{− δ},δ = 1 − 2, while the outer part has an exponential decline in σ. For the CB diskbinary interaction only the inner part of the disk is important. The outer disk radius here then represents the radius where the exponential drop in mass distribution starts, and the maximum disk mass represents the mass in this inner part of the disk. We assume an outer disk radius of 250 AU. See for example Bujarrabal et al. (2007, 2013, 2015). Keep in mind that this outer radius is not the same as the six time the separation radius that is used in calculating the effective disk angular momentum.
The change in eccentricity due to Lindblad resonances can be given in function of ȧ/a, and depends on the eccentricity. For small eccentricities this can be calculated analytically. In the range the m = l resonance dominates, while in the region the m = 2, l = 1 resonance is dominant. For small e the analytic form of ė is (see also Lubow & Artymowicz 2000): (A.55)For average eccentricities the above equation is extrapolated according to an efficiency decreasing with 1 /e: (A.56)where the constants are c_{4} = − 1.3652 and c_{5} =−1.9504 so that Eq. (A.56) smoothly connects to Eq. (A.55), while going to 0 at e = 0.7. For larger eccentricities resonances that damp the eccentricity start to dominate, thus ė_{disk} = 0 in the range e> 0.7 (e.g. Roedig et al. 2011).
The change in angular momentum due to disk – binary interactions is: (A.57)where ȧ/a is given by Eq. (A.44).
The mass feeding the CB disk is the fraction of the mass lost from the binary system through the outer Lagrange point during RLOF. The CB disk massloss rate is determined by the maximum mass in the disk M_{disk,max}, and the life time of the disk τ_{disk}. Both are input parameters in our model. The rate in which the disk loses/gains mass is then: (A.58)In this implementation the disk will continue existing for a period τ_{disk} after RLOF stops.
© ESO, 2015