Mass Transfer in Binary Stars using SPH. I. Numerical Method
aa r X i v : . [ a s t r o - ph . S R ] N ov Draft version August 15, 2018
Preprint typeset using L A TEX style emulateapj v. 03/07/07
MASS TRANSFER IN BINARY STARS USING SPH. I. NUMERICAL METHOD
Charles-Philippe Lajoie & Alison Sills
Department of Physics and Astronomy, McMaster University,Hamilton, ON L8S 4M1, Canada
Draft version August 15, 2018
ABSTRACTClose interactions and mass transfer in binary stars can lead to the formation of many differentexotic stellar populations, but detailed modeling of mass transfer is a computationally challengingproblem. Here, we present an alternate Smoothed Particle Hydrodynamics approach to the modelingof mass transfer in binary systems that allows a better resolution of the flow of matter betweenmain-sequence stars. Our approach consists of modeling only the outermost layers of the stars usingappropriate boundary conditions and ghost particles. We arbitrarily set the radius of the boundaryand find that our boundary treatment behaves physically and conserves energy well. In particular,when used with our binary relaxation procedure, our treatment of boundary conditions is also shownto evolve circular binaries properly for many orbits. The results of our first simulation of mass transferare also discussed and used to assess the strengths and limitations of our method. We conclude thatit is well suited for the modeling of interacting binary stars. The method presented here representsa convenient alternative to previous hydrodynamical techniques aimed at modeling mass transfer inbinary systems since it can be used to model both the donor and the accretor while maintaining thedensity profiles taken from realistic stellar models.
Subject headings: binaries: close — stars: evolution — hydrodynamics — methods: numerical INTRODUCTIONGravitation rules. It is what forms dark matter halosand giant molecular clouds. It is also what compressesthese clouds of gas to such an extent that stars formout of them. And these newly-formed stars are oftenbound to each other by the means of this force, formingbinary or multiple stellar systems. Although the detailsof binary star formation are still not fully understood(e.g. Tohline 2002), it is now acknowledged that stellarmultiplicity is more the rule than the exception. Obser-vations suggest that over 50% of the stars in our Galaxyare part of double, triple, quadruple, or even sextuplesystems (Abt & Levy 1976; Duquennoy & Mayor 1991;Fischer & Marcy 1992; Halbwachs et al. 2003). Becausestars grow in size considerably as they evolve, it is es-timated that those binaries with a period of less than ∼ days will inevitably interact at some point of theirlife (Paczy´nski 1971). When such close interactions oc-cur, material is transferred from one star to the other andthe course of evolution of the stars is irreversibly altered.Likewise, such close interactions can also be triggeredby stellar encounters in dense stellar environments (e.g.captures and exchanges; Pooley & Hut 2006).The first realization of the importance of binary in-teractions may have been by Crawford (1955), who sug-gested an interesting solution to the paradox of the Al-gol system, in which the more evolved star is also theleast massive. Crawford suggested that the close prox-imity of the two components must have led to significantmass transfer from the initially more massive star to theleast massive one until the mass ratio was reversed. Thisdiscovery opened the way to many more types of stars,such as cataclysmic variables and X-ray binaries, heliumwhite dwarfs, and blue stragglers, whose very existence Electronic address: [email protected],[email protected] could now be understood in terms of close binary evo-lution. Even for stars that are not transferring mass, aclose companion can have all sorts of effects on their ob-servable properties, such as increased chromospheric andmagnetic activity (e.g. RS Canum Venaticorum stars;Rodono 1992).Since the work of Crawford (1955), much effort hasbeen put into better understanding binary evolution andits by-products (Morton 1960; Paczy´nski 1965, 1971;Paczy´nski & Sienkiewicz 1972; Iben 1991). However, theusual Roche lobe formalism used to study binary starsapplies only to the ideal case of circular and synchronizedorbits. The addition of simple physics such as radiationpressure, for example, is enough to significantly mod-ify the equipotentials of binary systems (Dermine et al.2009) and estimates of fundamental parameters suchas the Roche lobe radius R L consequently become un-certain. Since the evolution of close binaries depends,among other things, on the rate at which mass transferproceeds and R L , analytical prescriptions for the masstransfer rate also become uncertain. The determinationand characterization of the mass transfer rate in binarystars therefore represents a key issue that needs to be ad-dressed, especially given that surveys of binary stars haveshown that a non-negligible fraction ( ∼
20% in the sam-ple of Petrova & Orlov 1999; see also Raguzova & Popov2005) of semi-detached or contact systems have eccentricorbits.Recent analytical work by Sepinsky et al. (2007a,b,2009), who investigated the secular evolution of eccen-tric binaries under episodic mass transfer, have shownthat, indeed, eccentric binaries can evolve quite differ-ently from circular ones. However, because mass transfercan occur on short dynamical timescales, it is necessaryto use other techniques to characterize it fully. In partic-ular, hydrodynamics has shown to be useful for studying LAJOIE & SILLStransient phenomena and episodes of stable mass trans-fer. Simple ballistic models (e.g. Warner & Peters 1972;Lubow & Shu 1975; Flannery 1975) and two-dimensionalhydrodynamical simulations of semi-detached binaries(e.g. Sawada et al. 1986; Blondin et al. 1995) have gen-erally been used in the past to study the general char-acteristics of the flow between two stars and the proper-ties of accretion disks. Later three-dimensional modelswith higher resolution allowed for more realistic stud-ies of coalescing binaries (e.g. Rasio & Shapiro 1994,1995) and accretion disks (e.g. Bisikalo et al. 1998) insemi-detached binaries, all mainly focusing on the sec-ular and hydrodynamical stability of binaries and onthe structure of the mass transfer flow. More recently,Motl, Tohline, & Frank (2002) and D’Souza et al. (2006)used grid-based hydrodynamics to simulate the coales-cence of n = 3 / . ⊙ ). They also investi-gated the onset of dynamical and secular instabilities inclose binaries and were able to get a good agreement withtheoretical expectations.Characterization of the accretion process and the be-haviour of the accreted material onto the secondary starrequires that the secondary be modeled realistically. Inmost simulations to date, the secondary has often beenmodeled using point masses or simple boundary condi-tions to approximate the surface accretion. Moreover, aspointed out by Sills & Lombardi (1997), the use of poly-tropes, instead of realistic models, may lead to signifi-cantly different internal structures for collision products,which may arguably be applicable to mass-transferringbinaries. Therefore, more work remains to be done inorder to better understand how mass transfer operatesin binary systems.Here, we present our alternate hydrodynamical ap-proach to the modeling of mass transfer in binaries. Wefirst discuss the usual assumptions made when study-ing binary stars §
2. Our computational method, alongwith our innovative treatment of boundary conditions,are introduced and tested in §
3. We discuss initial con-ditions for binary stars in SPH and the applicability ofour technique to binary stars in §
4. Our first mass trans-fer simulation is analyzed in §
5. The work presentedhere is aimed at better understanding the process of masstransfer and will be applied to specific binary systems inanother paper (Lajoie & Sills 2010; hereafter, Paper II). BRIEF THEORY OF BINARY SYSTEMSOur understanding of interacting binary stars isprompted by observations of systems whose existencemust be a result of mass transfer. These observationsand simple physical considerations have driven the theo-retical framework we use to model these systems.2.1.
The Roche approximation
Under the assumptions of point masses (M and M ),circular and synchronous orbits, the gravitational poten-tial in a rotating reference frame gives rise to equipo-tentials with particular shapes, as described in Eggleton(2006). The equipotential that surrounds both stars andintersects at a point between the two stars is the Rochelobe and forms a surface within which a star’s potentialdominates over its companion’s. Eggleton (1983) approx- imates the equivalent Roche lobe radius (R L ) to an ac-curacy of 1% over the whole range 0 < q < ∞ , by R L a ≈ . q / . q / + ln(1 + q / ) , (1)where a is the semi-major axis. Analytical studies ofbinary stars usually rely on this definition of the Rochelobe radius, but as discussed in §
1, such estimates for R L may not always be reliable.2.2. Rate of mass transfer
Once a star overfills its Roche lobe, mass is transferredto the companion’s potential well. The rate at whichmaterial is transferred is, in general, a strong function ofthe degree of overflow of the donor, ∆ R = R ∗ − R L , where R ∗ is the radius of the star. In general, mass is assumedto be transferred through the L point and, using simplephysical assumptions such as an isothermal and inviscidflow and an ideal gas pressure law, Ritter (1988) showsthat the mass transfer rate can be expressed as˙ M = − ˙ M exp (cid:16) R ph − R L H p (cid:17) . (2)Here, R ph is the photospheric radius, H P is the pres-sure scale height of the donor star, and ˙ M is the masstransfer rate of a star exactly filling its Roche lobe. Thismass transfer rate depends rather strongly on the de-gree of overflow ∆ R and goes to zero exponentially ifthe star is within its Roche lobe. Ritter (1988) pro-vides estimates for ˙ M of about 10 − M ⊙ yr − and H P /R ≃ − for low-mass main-sequence stars. Thismodel of mass transfer has been successfully applied tocataclysmic variables where the photosphere of the donoris located about one to a few H P inside the Roche lobe.For cases where the mass transfer rates are much larger,Paczy´nski & Sienkiewicz (1972) (see also Eggleton 2006and Gokhale et al. 2007) derives, in a similar way, themass transfer rate for donors that can be approximatedby polytropes of index n . In such cases, the mass transferrate is ˙ M = − ˙ M (cid:16) R ∗ − R L R ∗ (cid:17) n +3 / (3)where ˙ M is a canonical mass transfer rate which dependson M , M , and a . The dependence of the mass transferrate on (∆ R/R ) is again found, although somewhat dif-ferent than that of Ritter (1988). This rate is also zerowhen ∆ R ≤ Orbital evolution due to mass transfer
If the mass of the stars changes, then both the periodand separation ought to readjust. This behaviour can beASS TRANSFER IN BINARIES. I. 3shown by taking the time derivative of the total angularmomentum of a system of two point mass orbiting eachother with an eccentricity e :˙ L tot L = ˙ M M + ˙ M M −
12 ˙ MM + 12 ˙ aa − e ˙ e − e . (4)Equation 4 shows that as the masses of the stars changeand as mass and angular momentum are being lost fromthe system, both the orbital separation and the eccen-tricity change. The exact behaviour of these quantitiesdepends of course on the degree of conservation of bothtotal mass and angular momentum. For the usual as-sumptions of circular orbits and conservative mass trans-fer, we can further impose that e = 0, ˙ M = − ˙ M and˙ J tot = 0, therefore reducing Equation 4 to˙ aa = 2 ˙ M M (cid:16) − M M (cid:17) . (5)Assuming M to be the donor and more massive (i.e.˙ M < M > M ) we find that the separation a de-creases until the mass ratio is reversed, at which point theseparation starts increasing again. For main-sequence bi-naries, where the most massive star is expected to overfillits Roche lobe first, the separation is therefore expectedto decrease upon mass transfer.2.4. Limitations
The theoretical framework derived in this section hasgenerally been applied to the study of close binaries.However, strictly speaking, it is not valid in most in-stances. Close binaries are not all circular and synchro-nized, and the Roche lobe formalism therefore does notapply. This, in turn, makes estimates of mass trans-fer rates rather uncertain. Moreover, conservative masstransfer is more an ideal study case than a realistic oneand the secular evolution of binary system becomes acomplex problem. To circumvent these difficulties, ap-proximations to the mass transfer and accretion rates aswell as to the degree of mass loss have to be made. How-ever, to better constrain these approximation or avoidover-simplifications, one can use hydrodynamics, whichis well suited for modeling and characterizing episodes ofmass transfer. We therefore discuss our hydrodynamicstechnique and show how it can be used to better con-strain mass transfer rates in binary systems. COMPUTATIONAL METHOD3.1.
Smoothed Particle Hydrodynamics
Smoothed Particle Hydrodynamics (SPH) was intro-duced by Lucy (1977) and Gingold & Monaghan (1977)in the context of stellar astrophysics. Its relatively sim-ple construction and versatility have allowed for themodeling of many different physical problems such asstar formation (Price & Bate 2009; Bate et al. 1995),accretion disks (e.g. Mayer et al. 2007), stellar col-lisions (Lombardi et al. 1995; Sills et al. 1997, 2001),galaxy formation and cosmological simulations (e.g.Mashchenko, Couchman, & Wadsley 2006; Stinson et al.2009; Governato et al. 2009). Our code derives from thatof Bate (1995), which is based on the earlier version ofBenz (1990) and Benz et al. (1990). Here, we only em-phasize on the main constituents of our code. The reader is referred to these early works for complementary de-tails.SPH relies on the basic assumption that the value ofany smooth function at any point in space can be ob-tained by averaging over the known values of the func-tion around this point. This averaging is done using aso-called ‘smoothing kernel’ to determine the contribu-tion from neighbouring particles. The smoothing ker-nel can take many forms (see e.g. Price 2005); here weuse the compact and spherically symmetric kernel firstsuggested by Monaghan & Lattanzio (1985). To preventthe interpenetration of particles in shocks and allow forthe dissipation of kinetic energy into heat, we includean artificial viscosity term in the momentum and energyequations. The artificial viscosity can also take variousforms (e.g. Lombardi et al. 1999); we use the form givenby Monaghan (1989) with α = 1 and β = 2. We allow forthe smoothing length to change both in time and space,and we use individual timesteps for the evolution of allthe required quantities. In this work, we assume an equa-tion of state for ideal gases of the form P = ( γ − ρu ,where γ = 5 / ∼
24 CPUs for simulations of ∼ particles.3.2. Boundary conditions
As discussed by Deupree & Karakas (2005), the innerparts of stars in close binaries generally remain unaf-fected by the presence of a companion, and only thestructure of the outermost layers is modified by closetidal interactions. This result prompted us to modelonly the outer parts of the stars with appropriate bound-ary conditions. Such an approach effectively reduces thetotal number of SPH particles in our simulations with-out decreasing the spatial resolution. Conversely, for thesame amount of CPU time, modeling only the outermostlayers of stars allows for the use of more particles, there-fore enhancing the spatial and mass resolutions. More-over, CPU time is spent solely on particles actually tak-ing part in the mass transfer or being affected by thecompanion’s tidal field.SPH codes calculate hydrodynamical quantities by av-eraging over a sufficiently large number of neighbours.For particles located close to an edge or a boundary, twothings happen. First, since there are no particles on oneside of the boundary, a pressure gradient exists and theparticles tend to be pushed further out of the domainof interest. Second, if the number of neighbours for eachparticle is kept fixed by requirements, as it is in our code,then the smoothing length is changed until enough neigh-bours are enclosed by the particle’s smoothed volume.This lack of neighbours therefore effectively decreases thespatial resolution at the boundary and underestimatesthe particle’s density. In such circumstances, the imple-mentation of boundary conditions is required.3.2.1.
Ghost particles
Boundary conditions have often been implemented us-ing the so-called ghost particles, first introduced byTakeda, Miyama, & Sekiya (1994) (see also Monaghan1994). Ghost particles, like SPH particles, contributeto the density of SPH particles and provide a pressuregradient which prevents the latter from approaching or LAJOIE & SILLS
Fig. 1.—
Illustration showing how boundary conditions are im-plemented using fixed ghost particles. SPH and ghost particles arerepresented by solid and open dots respectively and the boundaryis represented by the thick solid line. The density calculation ofSPH particles includes the contributions from all particles locatedwithin two smoothing lengths of them. penetrating the boundary. Ghost particles can be cre-ated dynamically every time an SPH particle gets withintwo smoothing lengths of the boundary. When this oc-curs, the position of each ghost is mirrored across theboundary from that of its parent SPH particle (alongwith its mass and density). Therefore, the need forghosts (and a boundary) occurs only when a particlecomes within reach of the boundary. Here, however,we use a slightly different approach, based on the workof Morris, Fox, & Zhu (1997) and Cummins & Rudman(1999). The approach of these authors differs from themirrored-ghost technique in that the ghosts are createdonce, at the beginning of the simulation, and their rela-tive position remains fixed in time during the simulation.We further improve upon this technique in order to modelthe outer parts of self-gravitating objects. Starting fromour relaxed configurations, we identify any particles asghosts if they are located within three smoothing lengths inside of the boundary, which, at this point, is arbitrar-ily determined. Particles located above the boundaryare tagged as SPH particles, whereas the remaining onesare erased and replaced by a central point mass whosetotal mass accounts for both the particles removed andthose tagged as ghosts. Point masses interact with eachother and
SPH particles via the gravitational force only.Point masses are also used when modeling massive or gi-ant branch stars whose steep density profile in the coreis hard to resolve with SPH particles. Figure 1 illus-trates how our boundary conditions are treated in ourcode. Here, we use three smoothing lengths of ghostsas a first safety check in order to prevent SPH particlesfrom penetrating the boundary. We also enforce that noparticle goes further than one smoothing length insidethe boundary by repositioning any such particle abovethe boundary.We further ensure conservation of momentum by im-parting an equal and opposite acceleration to the central point mass (since the point mass and the ghosts movetogether), which we write as a jpm = − X i m i m pm a i (6)where a i is the hydrodynamical acceleration impartedto particle i from ghost j . This term is added to theusual gravitational acceleration of the central point mass.Ghosts are moved with the central point mass and aregiven a fixed angular velocity. Note that at this point, a i has already been calculated by our code so that this cal-culation requires no extra CPU time. Ghost particles arealso included in the viscosity calculations to realisticallymimic the interface.3.2.2. Applications to single stars
We now show that our new boundary condition treat-ment is well suited for the modeling of stars in hydro-static equilibrium. We first relax a 0.8-M ⊙ star withrotation ( ω = 0 .
10; solar units). The relaxation of a starrequires a fine balance between the hydrodynamical andgravitational forces, and therefore allows to assess the ac-curacy of our code. We model our stars using theoreticaldensity profiles as given by the Yale Rotational EvolutionCode (YREC; Guenther et al. 1992). SPH particles arefirst spaced equally on a hexagonal close-packed latticeextending out to the radius of the star. The theoreticaldensity profile is then matched by iteratively assigning amass to each particle. Typically, particles at the centre ofthe stars are more massive than those located in the outerregions, by a factor depending on the steepness of thedensity profile. As discussed by Lombardi et al. (1995),this initial hexagonal configuration is stable against per-turbations and also tends to arise naturally during therelaxation of particles. Stars are relaxed in our codefor a few dynamical times to allow for the configura-tion to redistribute some of its thermal energy and settledown. Once the star has reached equilibrium, we removethe particles in the central regions and implement ourboundary conditions. We give the star, the ghosts andthe central point mass translational and angular veloci-ties. The final relaxed configuration of our star, with theboundary set to ∼
75% of the star’s radius, is shown inFigure 2 along with the density and pressure profiles inFigure 3. By using our initial configuration for the setupof ghosts, we ensure that the ghosts’ position, internalenergy, and mass are scaled to the right values and thatthe ensuing pressure gradient maintains the global hy-drostatic equilibrium. The total energy (which includesthe gravitational, kinetic, and thermal energies) for themodel of Figure 2 evolved under translation and rotationis conserved to better than 1% over the course of one fullrotation. These results suggest that our treatment ofboundaries is adequate for isolated stars in translationand rotation. APPLICATIONS TO BINARY STARSWe now discuss the modeling of binary stars with ournew boundary conditions. In particular, we develop aself-consistent technique for relaxing binary stars andshow that our code, along with our new boundary con-ditions, can accurately follow and maintain two stars onrelatively tight orbit for many tens of dynamical times.ASS TRANSFER IN BINARIES. I. 5
Fig. 2.—
Example of a star modeled with ghost particles (smallred dots), SPH particles (black dots), which are found beyond theboundary (green line), and a central point mass (big red dot).
Fig. 3.—
Density and pressure profiles of the ghosts (red) andSPH particles (black) showing the gradients at the boundary, asexpected for realistic models of stars. In all panels, only the par-ticles located within 2 h of the equatorial plane are plotted. Theprofiles of the ghost particles exactly match that of the theoreti-cal profiles, showing that the overall profiles are very close the thetheoretical ones. We emphasize that the location of the boundary is, atthis point, arbitrary. We will discuss this issue in moredetails in § Binary star relaxation
As discussed in § Fig. 4.—
Logarithm of density in the orbital plane for a de-tached relaxed binary with stars of 0 .
80 M ⊙ fully modeled (left)and modeled with boundary conditions (right). binary systems for SPH simulations.To account for the different hydrostatic equilibriumconfiguration of binary stars, we use the fact that thestars should be at rest in a reference frame that is cen-tered at the centre of mass of the binary and that rotateswith the same angular velocity as the stars. A centrifugalterm is added to the acceleration of the SPH particles,which we find by requiring that it cancels the net gravita-tional acceleration of the stars (e.g. Rosswog et al. 2004;Gaburov et al. 2010). By assuming an initial orbital sep-aration, we find, at each timestep, the necessary angularvelocity Ω that cancels the stars’ net gravitational accel-eration. Note that we do not account for the Coriolisforce since the system is assumed to be at rest in the ro-tating frame. The net acceleration of the centre of massof star i is calculated in the following way: a i cm = P m j ( a hyd j + a grav j ) M i (7)where M i , a hyd and a grav are respectively the total massand the hydrodynamical and gravitational accelerationsof star i , and the summation is done over particles j that are bound to star i . We therefore get the followingcondition for the angular velocity:Ω i = P m j ( a hyd j + a grav j ) P m j r j . (8)where r j is the distance of particle j to the axis of ro-tation of the binary. We find the angular velocity forboth stars and take the average value and add it to theacceleration of each SPH particle. To ensure that theorbital separation is kept constant during the relaxation,we also reset the stars to their initial positions after eachtimestep by a simple translation. An example of sucha detached relaxed binary is shown in the left panel ofFigure 4. 4.2. Circular orbits
After the initial relaxation, the stars are put in an in-ertial reference frame by using the value of the angularvelocity obtained from the binary relaxation procedure.To assess the physicality and numerical integration ca-pabilities of our hydrodynamical code, we have evolved LAJOIE & SILLS
Fig. 5.—
Normalized energies and orbital separation as a functionof time for a 0.6+0.6 M ⊙ circular binary modeled with ∼ § different wide binaries on circular orbits for a small num-ber of orbits. This is important for simulations of masstransfer as any changes in orbital separation must bedriven by the mass transfer itself and not the initial con-ditions and/or numerical errors inherent to our method.Figure 5 shows the normalized orbital separation anddifferent energies as a function of time for a 0 .
60 + 0 . ⊙ binary with each star fully modeled (i.e. without aboundary). Each star contains ∼ ,
000 SPH particlesand the initial separation is 3 .
25 R ⊙ . Our relaxation pro-cedure yields an angular velocity of Ω = 0 . ∼ .
25% for over sixorbits. The slight (anti-symmetric) variations observedin the orbital separation and the kinetic energy are anindication that the system is on a slightly eccentric or-bit. At this level, however, we estimate that our codecan properly (and physically) evolve two stars orbitingaround each other.Figure 6 shows the comparison between the normal-ized energies and orbital separation of a 0 .
80 + 0 .
80 M ⊙ binary modeled with a different number of particles (thehigh resolution simulation is also shown in Figure 4).The solid and dotted lines correspond, respectively, tothe systems modeled without and with boundary con-ditions. In the simulations shown here, the boundaryis set at 75% of the star’s radius. Figure 6 shows thatintroducing the boundary conditions smoothes the os-cillations seen for the fully modeled binaries. This is Fig. 6.—
Normalized total energy and orbital separation for a0.80+0.80 M ⊙ circular binary modeled with (dotted line) and with-out (solid line) boundary conditions. The upper two panels are fora low-resolution simulation containing ∼ ,
000 (full star) parti-cles and the lower two panels are for a high-resolution simulationcontaining ∼ ,
000 particles (full star). explained by the fact that the calculation of the gravita-tional force on the central point mass, and hence mostof the star’s mass, is done using a direct summation (in-stead of a binary tree), therefore improving the accuracyand reducing the oscillations in the orbital separation.Figure 6 also shows that increasing the spatial resolutionincreases the quality of the evolution of the orbits. Forcomparison, the orbital separation of our high-resolutionsimulation varies by less than 0 . ∼ ∼ .
25% over the wholeduration of the simulations and the gravitational andthermal energies (not shown) remain constant to betterthan ∼ . .
80 + 0 . ⊙ binary system over four orbits and a relatively lownumber of particles. Note that only the primary is mod-eled with the use of our boundary conditions for rea-sons that are discussed in §
6. Again, the orbital sepa-ration remains fairly close to its initial value, to within2% at the end of four orbits. The total energy is con-served to better than 0 .
5% and only the kinetic en-ergy oscillates significantly. We note that, when com-paring our results of circular orbits with those of otherauthors, our binary relaxation procedure yields quan-titatively comparable orbital behaviours. The resultsof Benz et al. (1990) show oscillations of ∼ −
2% inthe orbital separation over three orbits whereas the sim-ASS TRANSFER IN BINARIES. I. 7
Fig. 7.—
Comparison of the normalized energies and orbitalseparation as a function of time for a 0.80+0.48 M ⊙ circular bi-nary modeled with (dotted line) and without (solid line) boundaryconditions. ulations from Dan, Rosswog, & Br¨uggen (2008) of twounequal-mass binaries showed a constant orbital sep-aration for many tens of orbits with an accuracy of ∼ ∼ . ∼
50% when modeled with boundary conditions.This depends of course on the location of the boundary,which at this point is arbitrary but chosen wherever itmakes the most sense for the problem at hand. Note,however, that the boundary should be placed at least afew smoothing lengths from the surface of the star.4.3.
Eccentric Orbits
The case of an eccentric orbit is interesting since tidalforces are time-dependent and can vary greatly depend-ing on the orbital phase. Here, we show that despitethe large tidal forces at periastron, our boundary con-ditions are well suited for the modeling of such systems.The system we model consists of two main-sequence starswith masses of 1 .
40 and 1 .
50 M ⊙ with an eccentricity of e = 0 .
15 and evolved for over four orbits ( P orb ≃
44 codeunits). The total number of particles nears ∼ , ∼
75% of the stars’radius, which, as shown in Figure 8, is deep inside the
Fig. 8.—
Radii enclosing different fractions of the total boundmass (in SPH particles) to the primary star as a function of timefor the 1 .
40 + 1 .
50 M ⊙ binary with e = 0 .
15. The dotted linerepresents the location of the boundary. star so that the effects of tidal force are negligible. In-deed, Figure 8 shows the radius of the primary enclosingdifferent fractions of the total bound mass (in SPH par-ticles) as a function of time for our binary system. Forexample, the radii containing 60% to 90% of the totalbound mass are shown to not change significantly dur-ing the whole duration of this simulation. In fact, onlythe outer radius of the star, containing over 95% of thebound mass, oscillates during each orbit. Therefore, inthis case, the choice of the location of the boundary (dot-ted line) is well justified and Figure 8 shows that the useof our method for eccentric binaries is adequate. Replac-ing the core of a star with a central point mass and aboundary remains a valid approximation as long as theboundary is deep enough inside the envelope of the star. SIMULATION OF MASS TRANSFERWe now present the results from the simulation pre-sented in § § Determination of bound mass
To estimate the mass transfer rates, we have to de-termine the component to which every SPH particle isbound. We use a total-energy (per unit mass) criterion,as presented by Lombardi et al. (2006), and determine LAJOIE & SILLS
Fig. 9.—
Logarithm of the density in the orbital plane for the1 .
40 + 1 .
50 M ⊙ binary with e = 0 .
15. The orbital period is about44 time units and the central point masses are not shown. whether a particle is bound to the primary, the sec-ondary, or the binary as a whole. In particular, giventhat most of the mass of the two stars is contained inthe point masses, we use the latter as the main compo-nents to which particles are bound. For a particle to bebound to any of the components, we require its total en-ergy relative to the component under consideration to benegative. The total energy of any particle with respectto both the point masses and the binary’s centre of massis defined as E ij = 12 v ij + u i − G ( M j − m i ) d ij (9)where v ij and d ij are the relative velocity and separa-tion, respectively, between particle i and component j .Moreover, we require the separation d ij to be less thanthe current separation of the two centres of mass of thestars (in this case, the point masses). For particles thatsatisfy both of these criteria for both stellar components,we assign them to the stellar component for which the to-tal energy is most negative. If only the energy conditionis satisfied for the stellar components, or the energy withrespect to the binary is negative, the particle is assignedto the binary component. Finally, if the total relative en- Fig. 10.—
Change in bound mass onto the primary as a functionof time for the 1 . .
50 M ⊙ binary with e = 0 .
15. The number ofparticles transferred during each episode is ∼ ergy is positive, the particle is unbound and is assignedto the ejecta.5.2. Estimates of mass transfer rates
Using the total mass bound to the stellar componentsas a function of time, we can determine the mass transferrates for the system modeled. We use a simple approachto determine the instantaneous mass transfer rates basedon the difference of the total mass bound of each compo-nent between two successive timesteps, i.e.˙ M = M it − M it − ∆ t , (10)where M refers to the bound mass and the indices referto component i and timesteps t and t − − M ⊙ yr − . The number of particles trans-ferred during each episodes is 175 −
200 which, given themasses of the SPH particles, may limit our ability to re-solve lower mass transfer rates. Indeed, SPH requires aminimum number of neighbouring particles to calculatethe density and in cases where only a handful of particlesare transferred, the SPH treatment may not be adequate.Given the masses of the particles, low numbers of parti-cles therefore set lower limits to the mass transfer ratesthat our simulations can resolve. Using more particles ofsmaller masses would definitely allow for the resolution oflower mass transfer rates, as discussed in §
6. Apart fromthe main episodes of mass transfer, we also observe sec-ondary peaks occurring before the main episodes of masstransfer. Our simulation shows that some of the mate-rial lost both by the primary and the secondary falls backonto both components and this is what is observed here. DISCUSSIONAnalytical approaches used to study the evolution ofbinaries usually rely on prescriptions or approximationswhen dealing with mass transfer rates. To relax some ofthese approximations, hydrodynamical modeling can beused, although it remains a hard task for many physi-cal and numerical reasons. To circumvent some of thesedifficulties, we introduced an alternate technique usingboundary conditions and ghost particles to model onlyASS TRANSFER IN BINARIES. I. 9the outermost parts of both the donor and the accre-tor stars. The location of the boundary is arbitrar-ily set. Since only the surface material is involved inmass transfer, our approach allows for better spatial andmass resolutions in the stream of matter. Moreover, ourmethod allows for the modeling of both the accretor andthe donor simultaneously while using less CPU time andmaintaining realistic density profiles, taken from stellarmodels.Our code was shown to work particularly well for starsof equal mass and stars that are centrally condensed. In-deed, replacing the dense core of a massive star by a pointmass is a good approximation. This is generally true forstars with masses & . ⊙ , although the evolutionarystage of the star may modify its density profile. Low-mass stars have shown to be more difficult to properlymodel with our approach (see Figure 7) and we suggestthat limiting our new boundary conditions to centrallycondensed stars will in general yield better results.We also discussed the setup of proper initial conditionsfor modeling binary stars, which we have tested on bothequal- and unequal-mass binaries. We demonstrated thatour relaxation procedure is consistently implemented inour code and that it allows for the evolution over manyorbits of equal-mass detached circular binaries and main-tain their orbital separation to within 1 − & − M ⊙ yr − ) are consistent with the esti-mates of Chen & Han (2008) who investigated the forma-tion of blue stragglers through episodes of mass transferonto main-sequence stars. We therefore suggest that ourmethod consisting of modeling both stars simultaneouslywith appropriate boundary conditions can be applied tothe problem of mass transfer in main-sequence binariesand help clarify the origin of blue stragglers.Using particles of lower mass allows for the modeling oflower mass transfer rates, although this requires the useof more particles and CPU time. Pushing the boundaryfurther out or using a point mass to model the secondary(as in Church et al. 2009) can also allow for the use ofmore particles and a better mass transfer rate resolu-tion. Using a point mass as the secondary would howevercounter the benefit of our method to be able to modeltwo interacting stars simultaneously.The physicality of our simulations depends of courseon the physical ingredients we put in our code. As such,we do not include the effects of radiation pressure andenergy transport mechanisms by radiation or convec-tion. These effects may have significance especially whenstudying the long-term evolution of mass-transferring bi-naries, where radiative cooling in the outer layers of thestars and envelope might be more important. Moreover,like any numerical technique, our method has some of itsown limitations, and we discuss them now.6.1. Solid boundary
By construction, our boundary is “semi-impermeable”,in that it does not allow particles to go through it. We setthree smoothing lengths of ghosts and enforce that par-ticles be artificially repositioned if they happen to crossthe boundary. However, it is possible that these con-ditions fail when dealing with large mass transfer rates and if some particles find themselves inside the bound-ary, around the central point mass, our code has to bestopped. This particle penetration limits our ability toadequately model episodes of extremely (and unrealisti-cally) large mass transfer rates ( & − M ⊙ yr − ; seePaper II). However, we do not think that our boundaryshould be so particle-tight since we do expect some mix-ing in the envelope of the secondaries. Although movingthe boundary to a smaller radius could fix the issue ofparticle penetration, it would counter the use and ben-efits of our approach. Instead, we suggest using sinkparticles at the centre of the stars in order to account fordeep mixing. Sink particles are like point masses but,in addition, their mass and momentum are allowed toincrease as SPH particles get accreted.Also, as discussed in § − years) than the duration of our simulations.6.2. Relaxation of unequal-mass binaries
Our relaxation procedure for binary stars has proven tobe especially efficient for equal-mass binaries. However,for unequal-mass binaries, we do not quite achieve thesame level of accuracy for the evolution of orbital sep-aration. We think the reason for this difference comesfrom the fact that the equal-mass systems we model areperfectly symmetric, i.e. the two stars are exact replicasof each other, whereas in the case of unequal-mass bina-ries, symmetry is broken. For equal-mass binaries, thegravitational acceleration calculations are exactly equaland opposite and the two stars are evolved identically.But given the adaptive nature of our code, stars (andparticles) of different masses may be evolved on differ-ent timesteps and care should be taken if the two stars(and particles) are to be evolved consistently. We per-formed test runs during which we forced our code to usea smaller common timestep, but this approach does notimprove the results. On the other hand, using a more di-rect summation approach for the calculation of the grav-itational force (by using a smaller binary tree openingangle) is found to improve the results. Indeed, resultsfrom test runs using such an approach show much im-provements in evolving stars on circular orbits. However,doing so also makes our simulations significantly longerto run. Therefore, we think the observed behaviour ofour unequal-mass binaries is the result of our calcula-tion of gravitational acceleration through a binary tree,0 LAJOIE & SILLSalthough more work remains to be done.6.3.