Issue 
A&A
Volume 550, February 2013



Article Number  A100  
Number of page(s)  13  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201220327  
Published online  01 February 2013 
BINSTAR: a new binary stellar evolution code
Tidal interactions
^{1} Institut d’Astronomie et d’Astrophysique, Université Libre de Bruxelles, ULB CP226, 1050 Brussels, Belgium
email: siess@astro.ulb.ac.be
^{2} Argelanderinstitut für Astronomie, University of Bonn, Germany
^{3} MaxPlanckInstitut für Astrophysik, Karl Schwarzschild Str. 1, 85741 Garching, Germany
Received: 3 September 2012
Accepted: 6 December 2012
We provide a detailed description of a new stellar evolution code, BINSTAR, which has been developed to study interacting binaries. Based on the stellar evolution code STAREVOL, it is specifically designed to study low and intermediatemass binaries. We describe the stateoftheart input physics, which includes treatments of tidal interactions, mass transfer and angular momentum exchange within the system. A generalised Henyey method is used to solve simultaneously the stellar structure equations of each component as well as the separation and eccentricity of the orbit. Test simulations for cases A and B mass transfer are presented and compared with available models. The results of the evolution of Algol systems are in remarkable agreement with the calculations of the Vrije Universiteit Brussel (VUB) group, thus validating our code. We also computed a large grid of models for various masses (2 ≤ M/M_{⊙} ≤ 20) and seven metallicities (Z = 0.0001, 0.001, 0.004, 0.008, 0.01, 0.02, 0.03) to provide a useful analytical parameterisation of the tidal torque constant E_{2}, which allows the determination of the circularisation and synchronisation timescales for stars with a radiative envelope and convective core. The evolution of E_{2} during the main sequence shows noticeable differences compared to available models. In particular, our new calculations indicate that the circularisation timescale is constant during core hydrogen burning. We also show that E_{2} weakly depends on core overshooting but is substantially increased when the metallicity becomes lower.
Key words: binaries: general / stars: evolution / stars: interiors / accretion, accretion disks
© ESO, 2013
1. Introduction
The study of binary stars is fundamental to many areas of astronomy. Binary stars allow a precise determination of stellar masses and a number of astrophysical phenomena (Type Ia supernovae, novae, Xray bursts and possibly gammaray bursts, etc.) can only be explained as the result of binary evolution. However, comprehensive and detailed stellar models of various physical processes involved in binarystar evolution are lacking and suffer many theoretical uncertainties associated with mass transfer, chemical mixing and orbital evolution, especially in the domain of low and intermediatemass stars (Pols 1994; Nelson & Eggleton 2001).
The stateoftheart in binary modelling consists of two approaches. On one hand, synthetic binary models (Nelemans et al. 2001; Hurley et al. 2002; Izzard et al. 2006; Belczynski et al. 2008) merge the results of detailed singlestar modelling, but fitted to formulae or tabulated for fast evaluation, with binaryevolution algorithms. While these are excellent tools to explore the binarystar parameter space in statistical studies, they still assume each star in the binary evolves like an isolated, single star. They also cannot tell us much about the interior structure of the stars and do not properly follow mixing of chemicals or angular momentum transport. On the other hand, binary stellar evolution codes have been developed to follow in detail the evolution of each component of the system. Most of these simulations consider relatively massive stars (e.g. De Loore & Doom 1992; Pols 1994; Wellstein et al. 2001; Nelson & Eggleton 2001; Petrovic et al. 2005) and focused on mass transfer prior to the supernova explosion (e.g. Podsiadlowski et al. 1992), on the formation of gammaray bursts (Cantiello et al. 2007) or more specifically on helium white dwarfs (Benvenuto & De Vito 2003).
Models of the evolution of low and intermediatemass binaries for various initial masses and metallicities are missing and BINSTAR has been designed to this end. Low and intermediatemass binary systems are of great interest because they often display unusual chemical compositions. Good examples are the Ba and the CH stars, rich in barium and carbon respectively, which are thought to have been polluted by material from a companion star which is now a white dwarf (e.g. Jorissen 1999; Jorissen et al. 1998). Similarly, the carbonenhanced metalpoor (CEMP) stars probably result from binary mass transfer involving a carbonrich companion star (Masseron et al. 2010), which are important in studies of stellar archaeology and nearfield cosmology. Postasymptotic giant branch (postAGB) stars also represent a challenge to modellers, as their evolution is thought to involve circumbinary discorbit interactions (e.g. Dermine et al. 2013).
The presence of a companion star is also essential for the understanding of nova explosions which are an important source of some rare isotopes such as ^{15}N and ^{17}O (e.g. José& Hernanz 1998) or Xray binaries. Many outstanding questions related to the effect of binarity on lightelement production (e.g. in ^{3}Herich stars), on the pulsation properties of stars (e.g. Kamath et al. 2010) or whether mass transfer is conservative or not (e.g. van Rensbergen et al. 2011) need to be addressed with a set of detailed low and intermediatemass (M ≲ 8 M_{⊙}) binarystar models. BINSTAR will be the tool for that study.
The paper is organised as follows: in the next two technical sections, a detailed presentation of the stellar and binary physics is presented. Section 4 discusses some numerical aspects of the code, benchmark applications as well as a study of the tidal torque constant are presented in Sect. 5 before concluding.
2. Input stellar physics
BINSTAR is an extension of the 1D stellar evolution code STAREVOL, that handles the simultaneous calculation of the binary orbital parameters (separation and eccentricity) as well as the two stellar components. The stellar input physics is the same as described in Siess (2006). In short, the equation of state derives from Pols et al. (1995) which is based on a minimisation of the Helmholtz free energy and that takes into account the effects of pressure ionisation and Coulomb shielding.
The radiative opacity tables are supplied by Iglesias & Rogers (1996) above 8000 K, supplemented at lower temperature by Ferguson et al. (2005). Conductive opacities are computed from Hubbard & Lampe (1969); Iben (1975); Raikh & Iakovlev (1982); Itoh et al. (1983); Mitake et al. (1984). The effect of carbon enrichment in the low temperature molecular opacities in carbon rich stars is treated analytically following Marigo (2002).
Concerning the nuclear physics, the neutrino energy loss rates are given by Itoh et al. (1996) and take into account the effects of plasma, pair, bremsstrahlung, recombination, and photo neutrino emission. Screening factors are calculated using the formulation of Mitler (1977). Our nuclear network includes N_{sp} = 53 species (from neutrons to ^{37}Cl) coupled to a network of 185 nuclear reactions taking into account all relevant (n, p, αcapture), weak (electron capture, βdecay) and electromagnetic (photodisintegration) interactions. This network ensures accurate nuclear energy production up to neon burning.
Neutrons are followed approximately by use of a fake neutron sink particle (Jorissen & Arnould 1989). If available, charged particle nuclear reactions rates are preferentially taken from Iliadis et al. (2010), then Angulo et al. (1999) else Caughlan & Fowler (1988). Neutron capture rates are from Bao et al. (2000) and experimental beta decay rates are taken from Horiguchi et al. (1996). Nuclear burning and chemical transport are solved simultaneously after the structure has converged and the evolution of the composition is governed by a system of N_{sp} coupled differential equations: (1)for k = 1,N_{sp}, where Y_{i,j,k,l} are the molar mass fractions of species i, j, k and l, ⟨ σ v ⟩ _{kl} and ⟨ σ v ⟩ _{ij} is the average nuclear cross section between the particles k,l and i,j respectively, and D is the diffusion constant.
The treatment of convection is based on the mixing length theory (MLT) as described in Cox & Giuli (1968) with the parameter α_{mlt} = 1.75 calibrated from our solar fit model. Nonstandard mixing processes can also be taken into account, in particular the treatment of thermohaline mixing as described in Siess (2009), diffusive overshooting (Siess et al. 2002) and semiconvection (Langer et al. 1985). In case of convective mixing, the diffusion coefficient in Eq. (1) is given by (2)where H_{p} is the local pressure scale height and v_{conv} the convective velocity as given by the MLT.
3. Input binary physics
In the following subsections, we describe the key input binary physics implemented in BINSTAR. The evolution of the orbital angular momentum is presented in Sect. 3.1, while the treatment of mass transfer via Rochelobe overflow (RLOF) and the resulting torques applied to each star are given in Sect. 3.2. A corresponding prescription for mass transfer via stellar winds is described in Sect. 3.3. The calculation of tides is given in Sect. 3.4, and Sect. 3.5 provides a summary of the different mass transfer processes.
3.1. Angular momentum evolution
The binary system consists of a donor (subscript “d”) and gainer (subscript “g”) of mass M_{d} and M_{g}, radii R_{d} and R_{g}, spins Ω_{d} and Ω_{g}. The total angular momentum (AM) of the binary star system (J_{Σ}) is the sum of the orbital (J_{orb}) and stellar angular momentum of the gainer (J_{g}) and donor (J_{d}) (3)where (4)with a the binary separation and e the eccentricity. The time derivative of Eq. (4) gives the evolution of the binary separation (5)To compute the separation at the next timestep (a(t + Δt) = a(t) + ȧΔt), one has to determine the various mass transfer rates (Ṁ_{d} and Ṁ_{g}), the rate of change of the orbital angular momentum () as well as ė. The quantity is calculated by imposing AM conservation i.e. (6)where are the torques applied to the star and the rate of change of the system angular momentum as should be accounted for in nonconservative evolution.
3.2. Mass and angular momentum transfer via Rochelobe overflow
3.2.1. Mass loss rate from the donor star
The Rochelobe radius for the donor star, R_{L1,d}, is determined using the formalism described by Eggleton (1983), i.e. (7)where the mass ratio (8)The Rochelobe radius of the accretor can be found by inverting the value of q in Eq. (7).
We treat mass transfer via RLOF according to the prescription described in Kolb & Ritter (1990). In this formalism, the authors make the distinction between 2 regimes depending whether the donor’s Roche radius lies within the optically thick or thin region of the star. Mass loss from the optically thin region of the donor star (where the optical depth τ < 2/3) is calculated using the formalism given by Ritter (1988), i.e. (9)where is the pressure scale height at radius R_{L1,d}. Ṁ_{0} is the mass transfer rate if the donor star exactly fills its Roche lobe, and is given by (10)Here, all quantities refer to the donor: T_{eff,d} is the effective temperature, ρ_{ph,d} the photospheric density, μ_{ph,d} the mean molecular mass at the photosphere and ℛ is the ideal gas constant. F(q) is a function of the mass ratio and is given by (11)In turn, g(q) depends on the distance from the centre of mass of the donor star to the L_{1} point x_{L}, in units of the binary separation a, and is given by (12)Finally, the pressure scale height at the location of the L_{1} point, is calculated from (13)where H_{P,ph} is the pressure scale height at the donor’s photosphere and γ(q) a function of the mass ratio q(14)However, Eq. (9) is only valid if the Roche lobe radius and the donor radius are (approximately) equal, and if the two radii move in step, i.e. Ṙ_{d} = Ṙ_{L1,d}. Indeed, some studies calculate mass transfer iteratively until the condition R_{L1} = R_{d} is satisfied (e.g. Benvenuto & De Vito 2003; Kalogera et al. 2004). However, care must be taken when applying this procedure. This condition holds if the donor star has a “sharp” surface, i.e. H_{P,ph}/R_{d} ≪ 1. This is the case for main sequence stars, but not for giants or AGB stars due to their extended atmospheres. Pastetter & Ritter (1989) showed that applying Ṙ_{d} = Ṙ_{L1,d} to the study of mass loss from thermallypulsating AGB stars gives incorrect results. The other condition (i.e. that the donor and Roche lobe radius move in step) implies that the mass transfer rate is stationary, i.e. . This is not the case during turnon and turnoff of mass transfer (D’Antona et al. 1989).
If either of these conditions are not satisfied, the stellar radius significantly overfills its Roche lobe and the Roche lobe will lie in the optically thick region of the donor star, where τ > 2/3. In this situation, we calculate the mass transfer rate via (15)where P_{ph} and P_{L1} is the pressure at the photosphere and at the L_{1} point respectively, T is the temperature and μ is the mean molecular weight and θ is a function of the adiabatic exponent, Γ_{1} = (d lnP/d lnρ)_{ad}, and is given by (16)
3.2.2. Torque on the donor star
Material lost from the donor at the L_{1} point has the specific AM of the stellar surface where we can reasonably assume R_{d} = R_{L1,d} (if the radius significantly exceeds the Roche radius then other physics will take over, such as commonenvelope evolution or the formation of a circumbinary disc). The torque applied on the donor resulting from RLOF mass transfer is then (17)
3.2.3. Gainer star
In conservative evolution, all the mass leaving the donor is accreted by the gainer. However this may not be always the case as various hydrodynamical simulations of close (e.g. Sytov et al. 2007) or wide (Theuns & Jorissen 1993; Podsiadlowski & Mohamed 2007) binaries show. The presence of a hot spot can also contribute to the ejection of mass (van Rensbergen et al. 2008) and the observation of circumbinary discs in wide postAGB binaries (van Winckel 2003; de Ruyter et al. 2006) may also attest of this nonconservatism. To account for this process, we use the standard 0 ≤ β ≤ 1 prescription where the mass accretion rate of the gainer is given by (18)Material lost from the donor can either impact directly onto the surface of the gainer or distribute into an accretion disc. Assuming ballistic trajectories, we know that a particle, at its closest approach, reaches a radius R_{min} = 0.0425a(q[1 + q])^{1/4} (Ulrich & Burger 1976).
If R_{g} < R_{min} the stream goes around the star and after one orbit collides with itself. After multiple collisions the matter accumulates in a narrow ring about the circularisation radius (R_{circ}) and later spreads to form an accretion disc, assumed to be Keplerian (for a review see Frank et al. 2002) of extent R_{disc} ≈ 2 R_{circ} (Ulrich & Burger 1976). At the surface of the gainer, the rate of angular momentum deposition, i.e. the applied torque, is given by (19)where (20)On the other hand, if R_{g} > R_{disc} the stream leaving the L_{1} point impacts the star with a specific AM (21)where the radius vector R = (x,y) and the velocity v = (ẋ,ẏ) at impact are estimated by computing the ballistic trajectory of a test particle in the potential of the gainer (Flannery 1975). In Eq. (21), x_{g} = M_{g}/(M_{d} + M_{g}) is the distance from the centre of mass of the gainer to the centre of mass of the binary system, and Ω_{orb} is the orbital angular speed. Distances are expressed in units of the orbital separation and time in units of the orbital period. The resulting torque is (22)
3.3. Mass and angular momentum transfer via stellar winds
3.3.1. Wind loss
Mass is assumed to be lost in spherical shells by the wind leading to an AM loss rate (23)where is the wind mass loss rate from star i = 1,2 as prescribed by e.g. de Jager et al. (1988), Reimers (1975) or Vassiliadis & Wood (1993). The spherical wind carries away the specific AM of the surface layers as well as some fraction of the orbital AM. This contribution, although not explicit in the previous equation, is taken into account in our global AM conservation equation ( in Sect. 3.5)
3.3.2. Wind accretion
This is the opposite of wind loss: material is accreted in shells such that a positive surface torque is applied to the stellar surface. In the classical BondiHoyle scheme, the AM accreted could well be zero (Blondin & Raymer 2012), but instead we use the specific spin AM of the companion and accrete a fraction f_{Jacc} onto the considered star (24)where f_{Jacc} is a free parameter usually set to 1 and the BondiHoyle rate is given by (25)where , , and . The free parameter α_{BH} is of the order of 1.5 and β_{w} ≈ 1/8 (Hurley et al. 2002). The accretion rate depends on the companion’s mass loss rate and should only be used in the limit of fast winds (v^{wind} ≫ v_{orb}).
In case of slow winds, such as those occurring in red giant branch (RGB) or AGB stars, this prescription is not valid, although detailed hydrodynamical simulations (e.g. Theuns et al. 1996) suggest that values of α_{BH} ten times less the standard value could be used for slow winds. However, a better approach would consist of incorporating a more physically realistic prescription based on the results from hydrodynamical simulations. This remains to be implemented (Abate et al. 2011).
3.4. Tides
The gravitational forces acting on stars in a binary system induce a deformation of their structure and the appearance of tidal bulges. Because the stellar material is viscous, the bulges make an angle α with the line connecting the centre of the stars. As a result of that tidal lag, the gravitational attraction generates a torque on the bulges and force the synchronisation of the stellar rotation rates with the orbital motion and, on a longer timescale, to circularise the orbit. The orbital parameters change because of the dissipation (conversion) of the star’s rotational energy into heat. In the weak friction model, it is simply assumed that the angle α is proportional to Ω − ω (e.g. Zahn 2008), where ω is the orbital angular speed of the binary.
In convective stars, the kinetic energy of the large scale currents raised by the tides cascades down to ever smaller scales until it dissipates into heat due to viscous friction of the convective environment. In radiative stars, the main dissipation mechanism is radiative damping of gravity modes that are excited by the tidal potential (see Zahn 2005, for a review).
We separate the tides into the convective equilibrium tide and the radiative dynamical tide. There are two general approaches, those by Hut (1981) and Zahn (1975, 1977, 1989). The prescription of Hut is valid for any e, and relies on calculation of timescales from the prescription of Zahn. An improved prescription for the equilibrium tides which is valid for any e was subsequently proposed by Zahn (1989).
The Hut (1981) formalism, which is developed in the framework of the weak friction model, describes the evolution of the orbital parameters for any eccentricity e, stellar spin Ω and separation (a), the latter simply deriving from AM conservation. We have (26)and (27)where f_{2 − 5} are polynomials in e (see Appendix A), the momentum of inertia of the star, and M its total mass. The convectiveequilibrium and radiativedynamic tidal timescales are estimated below.
3.4.1. Tidal timescale for radiative stars
In radiative stars, Zahn (1977, Eq. (5.9)) gives (for ω = Ω) (28)where E_{2} is a constant that depends on the stellar structure (see Appendix B) and is the ratio of the mass of the companion to the mass of the considered star. Zahn (1977) also gives the synchronisation timescale for radiative stars (his Eq. (5.7)), (29)
3.4.2. Tidal timescale for stars with a convective envelope
Zahn (1989) provides an improved prescription of equilibrium tides to accommodate situations in which the orbital period becomes shorter than the convective turnover timescale, as may be the case for close binaries. In this context, the weak friction approximation breaks down and one has to account for the multiple Fourier component of the tidal potential. With this formalism, the synchronisation and circularisation timescales (Eqs. (20) and (21) of Zahn 1989) are given by (30)and (31)The factors λ_{lm} (see Appendix C) are similar in magnitude to the apsidal motion constant k_{2} (which was used previously in their place) and the friction time (32)Fortunately, Zahn gives a prescription for their calculation in his paper which has been implemented in BINSTAR (see Appendix C below).
3.4.3. Summary
The computation of the tidal timescales allows us to estimate the torque applied onto each star i(33)as well as the contribution of star i to the rate of change of eccentricity ė_{i}, where I_{i} is the moment of inertia. The rate of change of the orbital eccentricity is then (34)The change in ȧ due to the tides is accounted for by the conservation of angular momentum (Sect. 3.5).
Fig. 1 Different modes of mass exchange between stars in a binary system; see text for details. 
3.5. Mass and angular momentum conservation
Figure 1 summarises the complex mass transfers taking place in a binary system. The mass lost by the system during Δt is (35)or, expressed in terms of the mass transfer rate, (36)The rate of angular momentum loss from the system is (37)where (38)is the specific angular momentum and f_{Σ} is a free parameter characterising the specific angular momentum of the ejected material.
When solving for the evolution of the separation (Eq. (5)), the unknown is computed using Eq. (6) where the rate of change of stellar angular momenta (i = 1,2) are given by (39)
4. Numerical considerations
BINSTAR uses a staggered grid to provide spatial centring of the stellar structure equations (e.g. Bodenheimer et al. 2007). The dependent variables are the natural logarithm of the stellar radius (lnr) and of the temperature (lnT), the variable lnf from which the EOS depends (see Pols et al. 1995) and the normalised luminosity (L_{m}/L_{norm}) where L_{m} is the stellar luminosity at mass coordinate m and L_{norm} the maximum value of  L_{m}  over the discretised structure. The radius r and luminosity L_{m} are defined at the cell boundary while T, lnf and the composition are defined at the centre of the cell. In our calculations, the evolution of the chemical abundances is done after the stellar structures have converged, i.e. the structure and nucleosynthesis equations are decoupled. This approximation is valid provided the timestep Δt is adequately chosen.
We use various criterion to constrain Δt based on the structural, nuclear, orbital and mass transfer timescales. Formally the structural timestep is determined from the time evolution of relevant variables A = { r, lnf, T, ρ,L } throughout all stellar grid points (shells k) during the previous timestep Δt^{n}(40)where ℱ_{A} ≈ 0.05 represents the amount of change which is acceptable for the quantity A. The nuclear timestep is calculated as (41)where ℱ_{i} (a few percent) has the same meaning as the parameter ℱ_{A}, Y_{i} is the molar abundance of species i where i = H, He, C or Ne and is the rate of destruction of Y_{i} and is estimated using the nuclear reaction rates of the slowest reactions of the H, He, C and Ne burning chains.
The timestep is also constrained by the requirement that the orbital elements (e and a) do not vary too rapidly (42)Outside the hydrogen and helium main sequence phases, the timestep is not allowed to exceed a fraction ~ 30% of the KelvinHelmholtz timescale (Δt_{KH} = GM^{2}/RL). The timestep is further constrained by the requirement that mass is not lost (accreted) by the donor (gainer) in too few timesteps. The maximum mass transfer timescale, Δt_{Ṁ} is constrained by (43)where ℱ_{acc} ≈ 0.01 and ℱ_{lost} ≈ 10^{4}. Finally, the new timestep is given by the minimum for the 2 stars of (44)Before we start the iteration process, the spatial resolution of each stellar component is adapted according to the structural changes. In short, shells are added whenever the relative difference between 2 consecutive shells δV = ln  V_{k + 1}/V_{k}  of one of the variables V = { r,  L  , lnf, T, P, ρ } exceeds a given threshold Δ_{max}. On the other hand, if at a given shell the δV for all the variables are less than Δ_{min} the shell is removed. Typically, Δ_{max} ≈ 1.2 meaning that we allow a difference of at most 20% between 2 shells and Δ_{min} = 0.02.
Once the timestep has been selected and the rezoning performed, a generalised Henyey method is used to solve simultaneously the stellar structure equations of the two stars as well as the evolution of the orbital parameters (Eqs. (5) and (34)). At each iteration during the convergence process, the mass transfer rates are updated according to the changes in the stellar structure and orbital parameters. Mass is added or removed from a mass coordinate M_{0} up to the surface. The value of M_{0} arbitrarily corresponds to the mass coordinate where the temperature first reaches 10^{6} K. In the layers where the mass has been modified, the Lagrangian time derivative of a variable A is calculated according to Neo et al. (1977), (45)where M is the stellar mass and q is the pseudoLagrangian coordinate, (46)The first and second term in the RHS of Eq. (45) are usually referred to as the nonhomologous and homologous terms, respectively. We also ensure that when mass is accreted, the composition of the surface layers corresponds to that of the accreted material, i.e. the matter is not mixed. The integration of the stellar structure equations is done in one shot from the centre to the surface where specific boundary conditions are applied. We have two options: either stop the integration at an optical depth τ = 2/3 and impose a blackbody photospheric boundary condition or use a grey atmosphere approximation. In the first case the surface boundary conditions at the outermost shell N are where L_{edd} is the Eddington limit, and the other symbols have their usual meaning. In the grey approximation, the boundary conditions write where a is the Stefan constant. We also emphasise that all the derivatives are evaluated analytically, except for the opacity and the dependence of ε_{nuc} on temperature. Finally, the code has been largely parallelised using OpenMP instructions, specifically for the calculation of the equation of state, stellar opacity and nucleosynthesis. This decreases CPU time by more than a factor of two.
5. Applications
To test our new code, we have computed the evolution of two close binaries undergoing case A (hydrogen core burning donor) and case B (hydrogen shell burning donor) mass transfer (for an exhaustive description of the various modes of mass transfer, see Eggleton 2006). Models were selected according to literature cases (and availability of the numerical data) for a comparative study. We also investigate the behaviour of the tidal torque constant E_{2} as a function of mass and metallicity. This parameter determines the circularisation and synchronisation timescales of stars with convective cores and radiative envelopes.
Fig. 2 Case A mass transfer. The initial binary parameters are M_{d} = 12 M_{⊙} (VUB solid black and ULB red dotted lines), M_{d} = 7.2 M_{⊙} (VUB short dashed green and ULB longdashed blue line) and P = 4d. The figure shows the evolution as a function of stellar mass of the stellar radius a), the mass transfer rate (b), left axis) and period (b), right axis), the evolutionary path in the HRD of the binary system c) and the time evolution of the mass transfer rate d). In panel d) the transfer time t_{transfer} is defined as t − 18.4 Myr for the VUB model and t − 29.3 Myr in our case. The final masses are 2.66 M_{⊙} and 16.54 M_{⊙}. 
5.1. Algol systems
In the following test cases we assume, for the sake of comparison, that the orbit is circular and that the stars are not rotating, implying that tidal interactions are not considered. Furthermore, the evolution is assumed to be fully conservative, i.e. and Ṁ_{d} + Ṁ_{g} = 0.
5.1.1. Case A mass transfer
In this calculation, we consider a 12 + 7.2 M_{⊙} system with an initial period of 4 days. Core overshooting is taken into account by artificially mixing the core composition with that of the shells located within a distance of α = 0.3 pressure scale height, H_{p}, beyond the Schwarzschild limit. We assume instantaneous mixing of that region and do not include wind mass loss, which has a negligible impact. The global properties of the evolution are illustrated in Fig. 2 and compared with a similar simulation performed by the Vrije Universiteit Brussel (VUB) group (van Rensbergen et al. 2008). The agreement between these independent calculations is remarkable and all the main features of the evolution are present.
Mass transfer begins while the donor is on the main sequence. As we initially have q > 1, then conservative evolution causes the orbit to shrink, ȧ/a < 0 (Eq. (5)), along with the orbital period and the Roche lobe radius (Eq. (7)). Mass loss is therefore initially unstable because even though the donor star loses mass, its radius decreases more slowly than the (shrinking) Roche radius. The mass transfer rate rises very rapidly as (R_{d} − R_{L1}) becomes larger and reaches in our simulation a maximum value of 5 × 10^{4} M_{⊙} yr^{1}. This phase of rapid growth lasts until the donor star restores thermal equilibrium, i.e after ~2 × 10^{5} yr which corresponds to the KelvinHelmholtz timescale of the primary at the beginning of mass transfer.
During this period of rapid mass transfer, the mass ratio is reversed (q < 1) after which the orbital period lengthens, reaching 30 days by the end of the simulation. The subsequent evolution is on a nuclear timescale, τ_{nuc}, and mass transfer is now driven by the radial expansion of the star, caused by core hydrogen burning, and Ṁ ≈ M/τ_{nuc}.
Once hydrogen is depleted in the core, central convection switches off, the donor contracts and for a short period of time the system becomes detached, ceasing mass transfer. This is marked in the HertzsprungRussell diagram (HRD) by the presence of a hook in the track near log (T_{eff}/K) = 4.2 and log (L/L_{⊙}) = 4 and by the discontinuity in the M − R relation around M ~ 5 M_{⊙}. However, as hydrogen shell burning switches on, the primary reexpands and a second mass transfer episode follows, termed case AB mass transfer (Ziółkowski 1970). During this subsequent mass transfer episode, the period of the system increases significantly.
Fig. 3 Similar to Fig. 2 but for Case B mass transfer. The initial binary parameters are M_{d} = 6 M_{⊙} (VUB solid black and ULB red dotted lines), M_{d} = 3.6 M_{⊙} (VUB short dashed green and ULB longdashed blue line) and P = 3.5 d. In panel d), t_{transfer} is defined as t − 75.1 Myr for the VUB model and t − 65.2 Myr in our case. The final masses are 0.86 M_{⊙} and 8.73 M_{⊙}. 
When helium ignites in the donor’s core, the structure contracts as the star reaches the helium main sequence and mass transfer stops. The evolution was terminated at neon ignition, before the donor went supernova (we note that the VUB model did not reach that evolutionary stage).
For the gainer star, the mass accretion timescale, M_{g}/  Ṁ_{g}  is shorter than its KelvinHelmholtz timescale, , and so the star is unable to maintain thermal equilibrium. The gainer therefore ascends the HRD with a higher luminosity compared to a main sequence star of the same mass.
At the end of the calculation, the final stellar masses are almost identical between the two model sequences. Compared with the VUB models, we obtain a slightly longer duration of the slow phase and a longer detached phase preceding case AB mass transfer. These minor differences are the consequence of the slightly higher mass transfer rate in the VUB model which is attributed to their somewhat larger donor radius, and hence a larger Roche lobe overfilling factor. Concerning the gainer, we note that near the peak mass transfer rate, the radius of the VUB gainer is larger and shifted towards higher mass. Again these effects are due to the higher mass accretion rate which drives the gainer even further out of thermal equilibrium and imposes a higher luminosity than in a relaxed configuration.
5.1.2. Case B mass transfer
The initial system configuration for this model is a 6.0+3.6 M_{⊙} with a period of 3.5 days. An overshooting of α = 0.23H_{p} is applied at the edge of the convective core and mass loss by stellar winds is neglected.
Following Kippenhahn & Weigert (1967), case B mass transfer starts once the more massive star has left its main sequence and crosses the Hertzsprung gap. Similarly to case A, it proceeds in two steps: a rapid phase during which the mass loss rate rises up to about 10^{4} M_{⊙} yr^{1} and the mass ratio is reversed, followed by a slow phase during which little mass is accreted but most of the orbital changes take place. During the phase of rapid massexchange, our mass transfer rate (Fig. 3) also shows the typical double peak feature (van der Linden 1987) where the second bump in Ṁ_{RLOF} is associated with the development of a thin surface convection zone. The maximum mass transfer rate which is approximately given by Ṁ ~ M/τ_{KH} (Paczyński 1971), is higher compared to a case A because when the donor star fills its Roche lobe, it has a larger radius and hence shorter τ_{KH}. Consequently the duration of the RLOF phase is also significantly shorter.
Another difference compared to case A is the maintenance of a relatively high mass transfer rate after the fast phase. This is also related to the structure of the star, where hydrogenshell burning is driving the expansion of the donor’s radius. In contrast to the corehydrogen burning donor in case A, case B evolution occurs on the nuclear timescale of hydrogenburning shell, τ_{HBS}, (Paczyński 1971) which is much shorter than τ_{nuc}. Because Ṁ ≈ M/τ_{HBS}, the mass loss rate is larger. Mass transfer ceases when helium eventually ignites: the donor then contracts and the star relaxes towards its helium main sequence location in the HRD.
Our simulation compares again very well with that of the VUB group. The mass transfer rate evolution is very similar; both the peak value and the duration of each phase is well reproduced and the final masses are again almost equal. In the MR diagram, we note that the VUB model experiences a short contact phase which is only just avoided in our simulation. The slightly larger initial stellar radius of the secondary VUB model is responsible for this feature and is due to overshooting (see below). One also notices that the stellar age at the onset of mass transfer (t_{transfer} in the figures) is always longer in our simulation. The reason for this difference is related to the treatment of core overshooting which produces a slightly larger convective core in the VUB model, thus extending the main sequence lifetime and the slow radius expansion. In our simulation, the donor star experiences a case BB RLOF during the short period of time that lasts between the end of core helium burning and carbon ignition (De Greve & De Loore 1977; Delgado & Thomas 1981).
5.2. Tidal timescale in main sequence stars with radiative envelopes
As presented in Sect. 3.4.1 for stars with a radiative envelope, tidal interaction operates through the dissipation of energy by gravity modes generated by the periodic tidal potential. In this theory, the characteristic tidal timescales τ_{circ} and τ_{sync} are inversely proportional to the quantity E_{2} (Eqs. (28) and (29)). The only available modern calculations that provide E_{2} values are the models computed by Claret (2004, 2005, 2006, 2007) which include core overshooting. Comparison between the 5 M_{⊙}, Z = 0.02 model of Claret (2004) and our 5 M_{⊙} with similar physics, reveals noticeable discrepancies. As shown in Fig. 4, our E_{2} values, at that metallicity and mass, are systematically lower by a factor varying from 10 up to 100 at the end of corehydrogen burning. In our models, the tidal torque constant decreases progressively as the convective core shrinks while in Claret (2004) it varies little until the end of the main sequence where it drops abruptly.
It is very difficult to understand the disagreement between these models considering that the evolutionary path in the HRD, core masses and main sequence lifetime are almost identical (compare the blue, longdashed and orange, short longdashed lines in Fig. 4). The differences are persistent over the considered mass range (2 M_{⊙} ≤ M ≤ 20 M_{⊙}) and tend to increase with stellar masses where the effect of wind loss on the structure become important. As mentioned by Claret (2004), the determination of E_{2} requires dedicated numerical attention because it involves the calculation of the derivative of the BruntVäisälä frequency at the edge of the convective core. To check if the differences were ascribed to some inadequate numerical resolution, we ran additional tests with increased resolution but we obtained the same results.
Fig. 4 Dependence of the tidal torque constant E_{2} on the input physics for a 5 M_{⊙} model. The solid (black) line corresponds to our standard Z = 0.02 model, the reddotted, greendashed and bluelongdashed curves refer to the overshooting models with parameters α = 0.05, 0.1 and 0.2, respectively. For illustration, we also plot the 5 M_{⊙}, Z = 0.0001 (cyandot–shortdashed line), Z = 0.004 (magentadotlongdashed line) models with no core overshooting as well as the result of Claret (2004) for M = 5 M_{⊙}, Z = 0.02 (orange, shortlongdashed). The mid and bottom panel display the evolution of radius and mass of the convective core for these models. 
In the top panel of Fig. 5, we display the evolution of E_{2}R^{9} which is inversely proportional to the circularisation time. The important point to mention is that, in contrast to the results of Khaliullin & Khaliullina (2010) who used models from Claret (2004), our circularisation timescale is not a monotonically decreasing function of time but instead remains relatively constant over the entire main sequence lifetime. This implies that tidal effects are constantly acting during corehydrogen burning with the same strength.
Fig. 5 Evolution of the tidal torque constant E_{2} (bottom panel), radius of the convective core (middle) and the quantity for solar composition models of mass M = 2, 2.5, 3.0, 4.0, 5.0, 7.0 10, 15 and 20 M_{⊙} as a function of the normalised main sequence lifetime. In each panel, the highest mass track corresponds to the top curve. In the lower panel, the solid lines on the upper right corner indicate the (constant) values of E_{2} for each mass as provided by Eq. (51). 
As emphasised in Claret (2004), the calculation of E_{2} (see Appendix B) is not straightforward and requires the knowledge of the stellar structure. To include the effects of tides in population synthesis calculations, Hurley et al. (2002) adopt the following fit (using the data provided by Zahn 1975) in terms of the stellar mass M(51)The short lines in the upper right corner at the bottom of Fig. 5 shows the values of E_{2} given by this formula for the considered masses. This parametrisation always overestimates E_{2} by at least a factor of ten and the discrepancy increases with time as E_{2} drops near the end of the main sequence with the reduction of the convective core size. Our average (52)with t_{MS} the main sequence lifetime, is a factor about 15 − 25 smaller than the value provided by Eq. (51). The Hurley et al. (2002) study largely underestimates the circularisation and synchronisation timescale and the assumption that at the time of RLOF the orbits will be circularised and the stellar components synchronised may not necessarily hold. In Appendix B we provide a new fit for this variable as a function of mass and metallicity that can easily be implemented in binary population synthesis codes.
Our models also allow us to study the dependence of E_{2} on stellar parameters (mass, age and core overshooting). For a given mass, E_{2} decreases with increasing metallicity. This is related to the structure of the convective core and in particular to its radial extent. Stars with a lower metal content have more compact core because, in order to sustain their higher luminosity, they have to contract more to increase the central temperature and the nuclear energy production. As shown in Fig. 4, including overshooting at the border of the convective core increases its mass but marginally its radius. As a consequence, because E_{2} is mostly dependent on the radius, core overshooting does not lead to a large increase in E_{2}, a conclusion also reached by Claret & Gimenez (1998).
6. Summary
This paper provides a detailed description of the binary physics implemented in the new stellar evolution code BINSTAR. One of the strength of this code is the simultaneous resolution of the stellar structure equations of each star along with the evolution of the orbital separation and eccentricity. The mass transfer rate is calculated at each iteration during the convergence process following the optically thin and thick case as described in Kolb & Ritter (1990). The prescriptions for the equilibrium and dynamical tides have been implemented following Zahn (1975, 1977, 1989) and our numerical scheme has been designed to ensure angular momentum conservation at all times.
Test simulations for case A and base B mass transfer have been computed and show remarkable agreement with the published models of the VUB group (van Rensbergen et al. 2008). The rate of mass transfer, the duration of RLOF phases, the final masses and the evolution in the HRD are well reproduced, providing a successful validation of our binary stellar evolution code.
We also investigated the effect of tides on main sequence stars with radiative envelope, i.e. for initial masses M_{init} ≥ 2 M_{⊙}. We find that our tidal torque constant E_{2} is about a factor 10 to 100 smaller than previous estimates by Claret (2004). We could not find the source of this discrepancy considering the apparent similarities in the stellar structures. Contrarily to Khaliullin & Khaliullina (2010), our selfconsistent models indicate that the circularisation timescales are relatively constant during the entire main sequence lifetime. This also implies the circularisation timescales computed by BINSTAR are longer than previously anticipated.
This code, with input physics inherited from STAREVOL, is mostly dedicated to the study of low and intermediate mass binaries. In the future, we will use this code to study mass transfer in eccentric binaries.
Acknowledgments
We are most grateful to Walter van Rensbergen and Jean Pierre de Greve for sharing their great expertise and for many fruitful discussions. L.S. is an FNRS Researcher and thanks the MaxPlanck Institut for Astrophysik in Garching for its hospitality during final elaboration of this work. P.J.D. and R.D. acknowledge support from the Communauté française de Belgique – Actions de Recherche Concertées as well as from FNRS fundings. This research has received funding from the Seventh Framework Programme of the European Community under grant agreement 220440.
References
 Abate, C., Pols, O. R., Izzard, R. G., Mohamed, S., & de Mink, S. E. 2011, in Evolution of Compact Binaries, eds. L. Schmidtobreick, M. R. Schreiber, & C. Tappert, ASP Conf. Ser., 447, 81 [Google Scholar]
 Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Bao, Z. Y., Beer, H., Käppeler, F., et al. 2000, Atomic Data and Nuclear Data Tables, 76, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Benvenuto, O. G., & De Vito, M. A. 2003, MNRAS, 342, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., & Raymer, E. 2012, ApJ, 752, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Bodenheimer, P., Laughlin, G. P., Rózyczka, M., & Yorke, H. W., eds. 2007, Numerical Methods in Astrophysics: an Introduction (Taylor and Francis Group) [Google Scholar]
 Cantiello, M., Yoon, S.C., Langer, N., & Livio, M. 2007, A&A, 465, L29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Caughlan, G. R., & Fowler, W. A. 1988, Atomic Data and Nuclear Data Tables, 40, 283 [Google Scholar]
 Claret, A. 2004, A&A, 424, 919 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A. 2005, A&A, 440, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A. 2006, A&A, 453, 769 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A. 2007, A&A, 467, 1389 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A., & Cunha, N. C. S. 1997, A&A, 318, 187 [NASA ADS] [Google Scholar]
 Claret, A., & Gimenez, A. 1998, A&AS, 133, 123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cox, J. P., & Giuli, R. T. 1968, Principles of stellar structure (New York: Gordon and Breach) [Google Scholar]
 D’Antona, F., Mazzitelli, I., & Ritter, H. 1989, A&A, 225, 391 [NASA ADS] [Google Scholar]
 De Greve, J.P., & De Loore, C. 1977, Ap&SS, 50, 75 [NASA ADS] [CrossRef] [Google Scholar]
 de Jager, C., Nieuwenhuijzen, H., & van der Hucht, K. A. 1988, A&AS, 72, 259 [NASA ADS] [Google Scholar]
 De Loore, C. W. H., & Doom, C. 1992, Structure and Evolution of Single and Binary Stars (Kluwer academic Publishers) [Google Scholar]
 de Ruyter, S., van Winckel, H., Maas, T., et al. 2006, A&A, 448, 641 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delgado, A. J., & Thomas, H.C. 1981, A&A, 96, 142 [NASA ADS] [Google Scholar]
 Dermine, T., Izzard, R. G., Jorissen, A., & Van Winckel, H. 2013, A&A, in press, DOI: 10.1051/00046361/201219430 [Google Scholar]
 Eggleton, P. P. 1983, ApJ, 268, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Eggleton, P. P. 2006, Evolutionary Processes in Binary and Multiple Stars, ed. P. Eggleton (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Flannery, B. P. 1975, MNRAS, 170, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics, eds. J. Frank, A. King, & D. Raine (Cambridge, UK: Cambridge University Press), 398 [Google Scholar]
 Horiguchi, T., Tachibana, T., & Katakura, J. 1996, in Japanese Nuclear Data Committee and Japan Atomic Energy Research Institute Nuclear Data Center [Google Scholar]
 Hubbard, W. B., & Lampe, M. 1969, ApJS, 18, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
 Iben, Jr., I. 1975, ApJ, 196, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Iliadis, C., Longland, R., Champagne, A. E., Coc, A., & Fitzgerald, R. 2010, Nucl. Phys. A, 841, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Itoh, N., Mitake, S., Iyetomi, H., & Ichimaru, S. 1983, ApJ, 273, 774 [NASA ADS] [CrossRef] [Google Scholar]
 Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Izzard, R. G., Dray, L. M., Karakas, A. I., Lugaro, M., & Tout, C. A. 2006, A&A, 460, 565 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jorissen, A. 1999, in Asymptotic Giant Branch Stars, eds. T. Le Bertre, A. Lebre, & C. Waelkens, IAU Symp. 191, 437 [Google Scholar]
 Jorissen, A., & Arnould, M. 1989, A&A, 221, 161 [NASA ADS] [Google Scholar]
 Jorissen, A., Van Eck, S., Mayor, M., & Udry, S. 1998, A&A, 332, 877 [NASA ADS] [Google Scholar]
 José, J., & Hernanz, M. 1998, ApJ, 494, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Kalogera, V., Henninger, M., Ivanova, N., & King, A. R. 2004, ApJ, 603, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Kamath, D., Wood, P. R., Soszyński, I., & Lebzelter, T. 2010, MNRAS, 408, 522 [NASA ADS] [CrossRef] [Google Scholar]
 Khaliullin, K. F., & Khaliullina, A. I. 2010, MNRAS, 401, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., & Weigert, A. 1967, ZAp, 65, 251 [Google Scholar]
 Kippenhahn, R., & Weigert, A. 1994 (New York: SpringerVerlag Berlin Heidelberg), also Astronomy and Astrophysics Library, Stellar Structure and Evolution, XVI, 468, 192 [Google Scholar]
 Kolb, U., & Ritter, H. 1990, A&A, 236, 385 [NASA ADS] [Google Scholar]
 Langer, N., El Eid, M. F., & Fricke, K. J. 1985, A&A, 145, 179 [NASA ADS] [Google Scholar]
 Marigo, P. 2002, A&A, 387, 507 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masseron, T., Johnson, J. A., Plez, B., et al. 2010, A&A, 509, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mitake, S., Ichimaru, S., & Itoh, N. 1984, ApJ, 277, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Mitler, H. E. 1977, ApJ, 212, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Nelemans, G., Yungelson, L. R., Portegies Zwart, S. F., & Verbunt, F. 2001, A&A, 365, 491 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nelson, C. A., & Eggleton, P. P. 2001, ApJ, 552, 664 [NASA ADS] [CrossRef] [Google Scholar]
 Neo, S., Miyaji, S., Nomoto, K., & Sugimoto, D. 1977, PASJ, 29, 249 [NASA ADS] [Google Scholar]
 Paczyński, B. 1971, ARA&A, 9, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Pastetter, L., & Ritter, H. 1989, A&A, 214, 186 [NASA ADS] [Google Scholar]
 Petrovic, J., Langer, N., & van der Hucht, K. A. 2005, A&A, 435, 1013 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Podsiadlowski, P., & Mohamed, S. 2007, Baltic Astron., 16, 26 [Google Scholar]
 Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246 [NASA ADS] [CrossRef] [Google Scholar]
 Pols, O. R. 1994, A&A, 290, 119 [NASA ADS] [Google Scholar]
 Pols, O. R., Tout, C. A., Eggleton, P. P., & Han, Z. 1995, MNRAS, 274, 964 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in C. The art of scientific computing, 2nd edn. (Cambridge: University Press) [Google Scholar]
 Raikh, M. E., & Iakovlev, D. G. 1982, Ap&SS, 87, 193 [Google Scholar]
 Reimers, D. 1975, Circumstellar envelopes and mass loss of red giant stars (New York: SpringerVerlag), 229 [Google Scholar]
 Ritter, H. 1988, A&A, 202, 93 [NASA ADS] [Google Scholar]
 Siess, L. 2006, A&A, 448, 717 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Siess, L. 2009, A&A, 497, 463 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Siess, L., Livio, M., & Lattanzio, J. 2002, ApJ, 570, 329 [NASA ADS] [CrossRef] [Google Scholar]
 Sytov, A. Y., Kaigorodov, P. V., Bisikalo, D. V., Kuznetsov, O. A., & Boyarchuk, A. A. 2007, Astron. Rep., 51, 836 [Google Scholar]
 Theuns, T., & Jorissen, A. 1993, MNRAS, 265, 946 [NASA ADS] [CrossRef] [Google Scholar]
 Theuns, T., Boffin, H. M. J., & Jorissen, A. 1996, MNRAS, 280, 1264 [NASA ADS] [CrossRef] [Google Scholar]
 Ulrich, R. K., & Burger, H. L. 1976, ApJ, 206, 509 [NASA ADS] [CrossRef] [Google Scholar]
 van der Linden, T. J. 1987, A&A, 178, 170 [NASA ADS] [Google Scholar]
 van Rensbergen, W., De Greve, J. P., De Loore, C., & Mennekens, N. 2008, A&A, 487, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Rensbergen, W., de Greve, J. P., Mennekens, N., Jansen, K., & de Loore, C. 2011, A&A, 528, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Winckel, H. 2003, ARA&A, 41, 391 [NASA ADS] [CrossRef] [Google Scholar]
 Vassiliadis, E., & Wood, P. R. 1993, ApJ, 413, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Wellstein, S., Langer, N., & Braun, H. 2001, A&A, 369, 939 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zahn, J.P. 1970, A&A, 4, 452 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1975, A&A, 41, 329 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1977, A&A, 57, 383 [NASA ADS] [Google Scholar]
 Zahn, J.P. 1989, A&A, 220, 112 [NASA ADS] [Google Scholar]
 Zahn, J.P. 2005, in Tidal Evolution and Oscillations in Binary Stars, eds. A. Claret, A. Giménez, & J.P. Zahn, ASP Conf. Ser., 333, 4 [Google Scholar]
 Zahn, J.P. 2008, in EAS Publ. Ser. 29, eds. M.J. Goupil, & J.P. Zahn, 67 [Google Scholar]
 Ziółkowski, J. 1970, Acta Astron., 20, 213 [NASA ADS] [Google Scholar]
Appendix A: Hut’s polynomials
The polynomials f_{2−5} used in Eqs. (26) and (27) are given by (Hut 1981) and (A.5)
Appendix B: Calculation of E_{2}
The equations derived by Zahn (1977) are summarised by Claret & Cunha (1997). The parameter E_{n} is given by (B.1)where f labels the edge of the convection zone, s the surface, x = r/R, g the gravity, R the stellar radius, M the stellar mass and γ_{n} is a numerical constant and (B.2)is the usual gamma function. The function B can be calculated from (known) stellar variables (B.3)where we have neglected the ϕ∇_{μ}/δ term given in Kippenhahn & Weigert (1994). It remains to calculate H_{n} which is defined by (B.4)where, for n = 2, X and Y are the solution of the differential equations (B.5)and (B.6)
B.1. Boundary conditions
The variable X is defined in Zahn (1970) Eqs. (36) and (37), (B.7)At x = 0 from Zahn’s Eq. (2b) we have (B.8)and (B.9)hence also dχ/dx = 0. Just before Eq. (36) Zahn (1970) gives the derivative (B.10)which is zero because χ = 0 at x = 0. The variable Y is defined by Eq. (33) of Zahn (1970) by (B.11)where m(x) is the mass coordinate. The function y_{1}(x) is introduced before Eq. (29) and has the property y_{1}x^{ − (n + 1)} → 1 as x → 0. Hence y_{1} → x^{3} because n = 2. This implies that at x = 0 we have Y = M x^{6}/m(x) = 0 because m(x) → 0 more slowly than x^{6} (approximately m(x) ~ x^{3} from integrating dm/dx = 4πx^{2}ρR^{3}). We can then differentiate Y to find (B.12)which is zero by a similar argument.
B.2. The X equation
We can rewrite (B.13)as x^{2}X′′ + Ax^{2}X′ − 6X = 0, where A = −dlnρ/dx. Now using a series solution (B.14)our differential equation becomes (B.15)Collecting like terms of x^{n + s} (n = 0... ), the x^{s} terms give a_{0}s(s − 1) − 6a_{0} = 0 i.e. s(s − 1) = 6 hence s = −2 or 3. If we use the a_{n} series with s = 3 and isolate the coefficients of the polynomials we obtain the following recursion relation (B.16)If we now use the b_{n} series with s = −2, X = b_{0}x^{2} + b_{1}x^{1} + b_{2} + b_{3}x + b_{4}x^{2} + ... the relation between the coefficients is the same as for a_{n} (Eq. (B.16)). The boundary condition X(0) = 0 implies that b_{0} = b_{1} = b_{2} = 0 and by recursion this implies that b_{n} = 0 ∀n , so the solution with s = −2 is simply X = 0 and is thus discarded. The first coefficient a_{0} is undefined but cancels in the final expression for H_{n} in Eq. (B.4). Note that y ~ x^{3} for small x (neglecting higher powers) as suggested by Zahn (1970). This gives us X′ ~ x^{2} and X′′ ~ x hence the differential equation can be integrated for small x and then numerically for all for 0 < x < 1.
B.3. The Y equation
We can rewrite the equation (B.17)as x^{2}Y′′ + AxY′ + BY = 0, where and . This can be solved exactly – it is an EulerCauchy equation. Let x = e^{t} then dx = e^{t}dt and hence (B.18)which gives us a new equation (B.19)This is a standard secondorder differential equation with constant coefficients. The auxiliary equation is (B.20)which has roots (B.21)and then Y = a_{1}e^{D1t} + a_{2}e^{D2t} hence (B.22)At x = 0 we have hence A = 0, B = −6 and D = 3 or − 2. The x^{2} solution must have zero coefficient, because Y = 0 at x = 0 hence Y(x = 0) = a_{1}x^{3}, which agrees nicely with the definition of Y in terms of y_{1} given above. Again the a_{1} is unknown but it cancels in the final expression for H_{n}, Eq. (B.4). Near x = 0 we have Y′ = 3a_{1}x^{2} and Y′′ = 6a_{1}x from which the differential equation can be integrated numerically for 0 < x < 1.
B.4. E_{2} fits
A LevenbergMarquardt algorithm (Press et al. 1992) was used to fit the values of E_{2} as a function of time t along the main sequence and stellar mass M for a given metallicity. The functional form of the fit is given by (B.23)where and with M in solar mass. The coefficients a_{i} are given in Table B.1 for each computed metallicity based on models computed without core overshooting. Our fits are accurate to within 3% of the computed value for M ∈ [1.6,20] . We also derived a fit to the main sequence lifetime as a function of mass and metallicity (B.24)where . This duration thus represents a lower limit of the main sequence lifetime. The fit is accurate to within 1.5% for Z ∈ [10^{4},0.03] and M ∈ [1.6,20] . We would like to stress that Eq. (B.24) was derived from models which did not include core overshooting.
Appendix C: λ_{lm} parameters of Zahn formalism
Zahn (1989) describes an alternative method for calculating the tidal torques in stars with convective envelopes. This involves dropping the k_{2} calculation in favour of the parameters λ_{lm} where l and m are integers which describe a given harmonic used in the expansion of the tidal potential and are related to the tidal timescale by Π_{lm} = 2π/lω − mΩ. According to Zahn (1989), we have (C.1)where x_{b} = r_{conv}/R is the base of the convective envelope. So in general two integrals are required, although in practice there is often no root, so no x_{a}, and just one integral is required (the second one, with limits x_{a} to 1). Note that unless the tidal timescale and orbital period are similar, all the λ_{lm} are essentially identical and the theory should be similar to that of the older 1970s papers which use k_{2} instead. The value of x_{a} is given by (C.2)where α′ = 0.762α_{MLT} and t_{f} (Eq. (32)) is the convective turnover timescale. When the tidal timescale is large compared to the orbital period the RHS of Eq. (C.2) is also large and there is no root to the equation, as would be expected. However, when the tidal timescale is small the RHS is also small so to determine if there is a root to the equation we simply go through the shells in the convective envelope looking for LHS − RHS < 0 and choose the shell for which the difference LHS − RHS is smallest, hence we have r = r_{a} at that shell, hence x_{a} = r_{a}/R. To calculate the parameter E, we use the fact that at the base of the convective envelope, the density is given by (C.3)The variable h is related to the pressure scale height by the relation H_{p} = Rhx^{2}/q where x = r_{conv}/R and q = m_{conv}/M where r_{conv} and m_{conv} are the radius and mass coordinate at the base of the convective envelope. Using , we obtain (C.4)Zahn suggests this has a value of 45.48 in a fully convective star. With our premain sequence models, which are fully convective, we obtained a value between 45 and 47.
All Tables
All Figures
Fig. 1 Different modes of mass exchange between stars in a binary system; see text for details. 

In the text 
Fig. 2 Case A mass transfer. The initial binary parameters are M_{d} = 12 M_{⊙} (VUB solid black and ULB red dotted lines), M_{d} = 7.2 M_{⊙} (VUB short dashed green and ULB longdashed blue line) and P = 4d. The figure shows the evolution as a function of stellar mass of the stellar radius a), the mass transfer rate (b), left axis) and period (b), right axis), the evolutionary path in the HRD of the binary system c) and the time evolution of the mass transfer rate d). In panel d) the transfer time t_{transfer} is defined as t − 18.4 Myr for the VUB model and t − 29.3 Myr in our case. The final masses are 2.66 M_{⊙} and 16.54 M_{⊙}. 

In the text 
Fig. 3 Similar to Fig. 2 but for Case B mass transfer. The initial binary parameters are M_{d} = 6 M_{⊙} (VUB solid black and ULB red dotted lines), M_{d} = 3.6 M_{⊙} (VUB short dashed green and ULB longdashed blue line) and P = 3.5 d. In panel d), t_{transfer} is defined as t − 75.1 Myr for the VUB model and t − 65.2 Myr in our case. The final masses are 0.86 M_{⊙} and 8.73 M_{⊙}. 

In the text 
Fig. 4 Dependence of the tidal torque constant E_{2} on the input physics for a 5 M_{⊙} model. The solid (black) line corresponds to our standard Z = 0.02 model, the reddotted, greendashed and bluelongdashed curves refer to the overshooting models with parameters α = 0.05, 0.1 and 0.2, respectively. For illustration, we also plot the 5 M_{⊙}, Z = 0.0001 (cyandot–shortdashed line), Z = 0.004 (magentadotlongdashed line) models with no core overshooting as well as the result of Claret (2004) for M = 5 M_{⊙}, Z = 0.02 (orange, shortlongdashed). The mid and bottom panel display the evolution of radius and mass of the convective core for these models. 

In the text 
Fig. 5 Evolution of the tidal torque constant E_{2} (bottom panel), radius of the convective core (middle) and the quantity for solar composition models of mass M = 2, 2.5, 3.0, 4.0, 5.0, 7.0 10, 15 and 20 M_{⊙} as a function of the normalised main sequence lifetime. In each panel, the highest mass track corresponds to the top curve. In the lower panel, the solid lines on the upper right corner indicate the (constant) values of E_{2} for each mass as provided by Eq. (51). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.