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 open-source 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 phase-dependent 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: Roche-lobe overflow

The instantaneous Roche-lobe 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 Roche-lobe 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 Roche-lobe radius of the accretor can be obtained by inverting the mass ratio in Eq. (A.9), thus replacing q_d with q_a.

Mass-loss is implemented in MESA following the prescription of Ritter (1988) and Kolb & Ritter (1990). This formalism differs between mass-loss from the optically thin region (R_d ≤ R_L,d), and that from the optically thick region (R_d>R_L,d). Mass-loss 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)Ṁ₀ is the mass loss rate when the donor star just fills its Roche-lobe: (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 & Meyer-Hofmeister 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 Roche-lobe, the Roche-lobe will lie in the optically thick region of the donor. The mass loss is then calculated via: (A.18)where Ṁ₀ 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₃ is a function that depends on the adiabatic exponent Γ₁, and is given by: (A.19)The mass loss equations depend on the eccentric anomaly through the Roche-lobe 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 re-emission). 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 = γ²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 Roche-lobe radius as the torque in a tidal friction model, and there is an expected saturation when co-rotation 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⁴, 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 Bondi-Hoyle 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 angular-momentum loss due to wind-mass loss, assuming a spherical-symmetric 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: phase-dependent 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 phase-dependent wind-mass loss through tidal interactions, and phase-dependent 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 phase-dependent 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 mass-loss, 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 tidally-induced large-scale 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 weak-friction 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₂₋₅ 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 time-scale (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 time-scale 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₂ 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 non-resonant interactions have been described by Goldreich & Tremaine (1979) and Artymowicz & Lubow (1994), by using a linear-perturbation 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 half-angular-momentum radius of the disk). The second assumption is that the non-axisymmetric 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⁵ 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 disk-surface-distribution 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α¹) 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 Stefan-Boltzmann 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 post-AGB or post-RGB binary. Post-AGB and post-RGB 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 disk-binary 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₄ = − 1.3652 and c₅ =−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 mass-loss 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