Bondi-Hoyle-Lyttleton accretion by binary stars
MMNRAS , 1–18 (2019) Preprint 30th October 2019 Compiled using MNRAS L A TEX style file v3.0
Bondi-Hoyle-Lyttleton accretion by binary stars
T.A.F. Comerford (cid:63) , R.G. Izzard , , R.A. Booth and G. Rosotti , Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, United Kingdom. Astrophyhsics Research Group, Faculty of Engineering and Physics, University of Surrey, Guildford, GU2 7XH, United Kingdom. Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands.
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Binary stars often move through an ambient medium from which they accrete materialand angular momentum, as in triple-star systems, star-forming clouds, young globu-lar clusters and in the centres of galaxies. A binary form of Bondi-Hoyle-Lyttletonaccretion results whereby the accretion rate depends on the binary properties: thestellar masses and separation, and the relative wind speed. We present the results ofsimulations performed with the hydrodynamic code gandalf , to determine the massaccretion rates over a range of binary separations, inclinations and mass ratios. Whenthe binary separation is short, the binary system accretes like a single star, while accre-tion onto stars in wide binaries is barely affected by their companion. We investigateintermediate-separation systems in some detail, finding that as the binary separationis increased, accretion rates smoothly decrease from the rate equal to that of a singlestar to the rate expected from two isolated stars. The form of this decrease depends onthe relative centre-of-mass velocity of the binary and the gas, with faster-moving bin-aries showing a shallower decrease. Accretion rates vary little with orbital inclination,except when the orbit is side-on and the stars pass through each others’ wakes. Thespecific angular momentum accretion rate also depends on the inclination but is neversufficient to prevent the binary orbit from contracting. Our results may be applied toaccretion onto protostars, pollution of stars in globular and nuclear clusters, and windmass-transfer in multiple stellar systems.
Key words: binaries: general – accretion – hydrodynamics – methods: numerical
Bondi–Hoyle–Lyttleton accretion (BHLA) describes how amassive object moving relative to a uniform fluid gravit-ationally focusses and accretes material. Based on the ori-ginal ballistic model of Hoyle & Lyttleton (1939), the BHLArate accounts for gas pressure and shock formation, andsmoothly interpolates between stationary and supersonic ac-cretors (Bondi & Hoyle 1944; Bondi 1952). For a detailedreview of BHLA see Edgar (2004). We The BHLA rate isapplied to a variety of situations, ranging in scale from windmass transfer in binary stars (Boffin & Jorissen 1988) andthe accretion of gas by stars in stellar clusters (Thoul et al.2002), to accretion by entire stellar clusters (Lin & Mur-ray 2007) and galaxies (Sakelliou 2000). In many such cases,the accretion rates predicted by theory closely match theobserved consequences.It is possible that the accretor in any of these contextsis a gravitationally bound binary, rather than single star sys- (cid:63)
E-mail: [email protected] tem. However, until recently, there has been little investiga-tion of BHLA onto binary star accretors. This is a problemfor accurate modelling of mass transfer in triple star systems,in which a wind from a giant star encounters a binary in or-bit around the giant. In particularly close binaries, this couldalso pose issues if a similar triple system ever undergoes com-mon envelope evolution in which all three stars orbit withina single envelope. Modifications to accretion rates causedby binary accretors could also serve to alter predictions ofstellar pollution in much the same way as wind accretionproduces stars like barium and extrinsic carbon stars.The Galactic halo object CS 22964-161 is a double-linedspectroscopic binary in which both stars are carbon andbarium rich dwarfs (Thompson et al. 2008). The system islikely a triple containing the relatively binary we see now,with a more distant white dwarf companion that transferredcarbon and barium to the binary in its wind when it ascen-ded the asymptotic giant branch (AGB; Abate et al. 2016).Soker (2004) seeks to explain the shapes of irregularplanetary nebulae as a result of BHLA onto a close binaryin a triple system with an AGB star. In this situation, they © a r X i v : . [ a s t r o - ph . S R ] O c t T.A.F. Comerford et al. find that it is possible for the stars in the binary to accretesufficient angular momentum to launch jets which shape thesurrounding AGB wind, however only limited considerationis given to the mass accretion, its effect on the binary andits subsequent evolution.We wish to explore the effects of binary BHLA in situ-ations not studied in that work by relaxing the assumptionthat the accretion only occurs within a narrow column, andthat the binary separation is much smaller than the Bondi-Hoyle radius, r a = GM v , (1)where M is the combined mass of the binary, and v is itscentre-of-mass velocity relative to the gas. For comparison,it is also convenient to define the Bondi radius, r B = GMc s , (2)where c s is the sound speed in the ambient gas.Lin & Murray (2007) consider the case of gas accretiononto a stellar cluster, modelled by a relatively shallow Plum-mer potential. They find that under the condition that thecluster’s internal velocity dispersion is slower than both itsmotion and the sound speed of ambient gas, bulk accretionby the cluster far outweighs the combined accretion rate ofits stars. Most of this gas is retained inside the cluster andis not, at least initially, accreted onto stars.Dimensional arguments show that accretion that isdominated by gravity has a rate, (cid:219) M ∝ M , where M isthe mass of the accreting object. If a binary star has twicethe mass of two otherwise identical single stars, this naivelyyields an extra factor of four in accretion rate. The accretedmass is then split between two stars and each accretes attwice the rate compared to if it were isolated. Naturally,this reasoning cannot hold in all binary systems, as somewill have separations wider than the characteristic lengthscale of BHLA flow, but it is reasonable to assume that theaccretion enhancement factor approaches four as the binaryseparation becomes arbitrarily short.Because this problem is analytically intractable, wehave carried out a suite of smoothed particle hydrodynamic(SPH) simulations using the code GANDALF (Hubber et al.2018). Our simulations probe a range of binary separations,velocities and inclinations, with the aim of determining anumerical prefactor to the analytic accretion rate which de-pends on the binary properties. The original description of Hoyle & Lyttleton (1939) con-siders a star, modelled as a point mass, moving relativeto a uniform ideal fluid at a constant, supersonic velocity.The gravity of the star focuses passing material into a one-dimensional shock in its wake. As gas decelerates into theshock, it loses speed in the direction perpendicular to thewake. This deceleration causes some material to becomegravitationally bound to the star, i.e. it will eventually be ac-creted. The initial impact parameter at which material willbecome marginally bound is the Hoyle-Lyttleton radius, r a .By considering the rate at which incoming gas with impact Figure 1.
A schematic of a typical binary system consideredin this paper. v orb represents the binary’s orbital velocity, while v CoM is the centre-of-mass velocity relative to the surrounding gas.Note that this figure shows an equal-mass binary, while we alsoconsider binaries with unequal masses (see § parameter less than r a approaches the star, the accretionrate is, (cid:219) M HL = π r ρ v = πρ G M v , (3)termed the Hoyle-Lyttleton accretion rate. Here, ρ is theambient gas density, M is the mass of the star, v is the rel-ative velocity of the star and gas, and G is the gravitationalconstant.This description suffers from several drawbacks, themost serious being that the accretion column has multival-ued velocity: material must simultaneously flow towards andaway from the star. It also treats the fluid as collisionless,except in the shock.A more realistic model was proposed in Bondi & Hoyle(1944) which introduces a finite shock width. A three-dimensional accretion column avoids the problem of mul-tivalued velocity, although it is still assumed that the fluidmoves ballistically outside the shock. A convenient propertyof this model is that the accretion rate at low relative velo-cities approaches the form of the Bondi rate, (cid:219) M B ≈ πρ G M c . (4)This results in an accretion rate which interpolates betweenthe Bondi rate and the supersonic Hoyle-Lyttleton rate, (cid:219) M BHL = × πρ G M (cid:16) v + c (cid:17) / , (5)known as the Bondi-Hoyle-Littleton (BHL) rate. While 4 im-plies that the numerical constant should be 2, it is not partic-ularly well constrained by the treatment in Bondi (1952). So,following the numerical simulations in Shima et al. (1985),we set it to 4 to match the hypersonic, i.e. large v , Hoyle-Lyttleton rate (Eq. 3).In systems with binary accretors we expect three re-gimes of accretion depending on the ratio of binary orbitalseparation, a , and the typical accretion radius, r a . If a (cid:29) r a ,the system behaves similarly to two single stars and accre-tion flows should be of the standard BHLA form onto each MNRAS000
A schematic of a typical binary system consideredin this paper. v orb represents the binary’s orbital velocity, while v CoM is the centre-of-mass velocity relative to the surrounding gas.Note that this figure shows an equal-mass binary, while we alsoconsider binaries with unequal masses (see § parameter less than r a approaches the star, the accretionrate is, (cid:219) M HL = π r ρ v = πρ G M v , (3)termed the Hoyle-Lyttleton accretion rate. Here, ρ is theambient gas density, M is the mass of the star, v is the rel-ative velocity of the star and gas, and G is the gravitationalconstant.This description suffers from several drawbacks, themost serious being that the accretion column has multival-ued velocity: material must simultaneously flow towards andaway from the star. It also treats the fluid as collisionless,except in the shock.A more realistic model was proposed in Bondi & Hoyle(1944) which introduces a finite shock width. A three-dimensional accretion column avoids the problem of mul-tivalued velocity, although it is still assumed that the fluidmoves ballistically outside the shock. A convenient propertyof this model is that the accretion rate at low relative velo-cities approaches the form of the Bondi rate, (cid:219) M B ≈ πρ G M c . (4)This results in an accretion rate which interpolates betweenthe Bondi rate and the supersonic Hoyle-Lyttleton rate, (cid:219) M BHL = × πρ G M (cid:16) v + c (cid:17) / , (5)known as the Bondi-Hoyle-Littleton (BHL) rate. While 4 im-plies that the numerical constant should be 2, it is not partic-ularly well constrained by the treatment in Bondi (1952). So,following the numerical simulations in Shima et al. (1985),we set it to 4 to match the hypersonic, i.e. large v , Hoyle-Lyttleton rate (Eq. 3).In systems with binary accretors we expect three re-gimes of accretion depending on the ratio of binary orbitalseparation, a , and the typical accretion radius, r a . If a (cid:29) r a ,the system behaves similarly to two single stars and accre-tion flows should be of the standard BHLA form onto each MNRAS000 , 1–18 (2019) inary BHLA star. In this case, one would expect the accretion rate to be, (cid:219) M = (cid:219) M + (cid:219) M = (cid:20) q ( + q ) + ( + q ) (cid:21) (cid:219) M BHL , (6)where (cid:219) M BHL is the BHL accretion rate for a star with a massequal to M + M (which we term the ‘single-star rate’).When a (cid:28) r a , the accretion flow should be similar toa point mass containing the total mass of the binary at theBHL rate, with deviations only occurring close the the stars.In this case, one might expect that accreted mass is dividedbetween the stars in the same ratio as Bondi accretion, withthe total accretion rate set by the point-mass BHL rate ( (cid:219) M in Eq. 5). Alternatively, if the accreted material has sufficientangular momentum, circumstellar or circumbinary discs mayform (a general investigation of BHLA onto a disc appearsin Moeckel & Throop 2009).In the intermediate case, where a ∼ r a , it is unclearwhich regime should dominate, as the two characteristic ve-locity scales, c s and the binary’s orbital velocity, v orb = (cid:114) GMa , (7)are approximately equal; that is to say, r a a ≈ v (cid:113) v + c s . (8)As a result, we have chosen to perform the majority of ourhydrodynamic simulations on systems in the range of orbitalseparation where the ratios in Eq. 8 are of order unity.To first order, we expect the intermediate-case accretionrate to exceed that of wide binaries because of the proximity,and added gravitational pull, of the companion. However,closer binaries have higher orbital velocities, which may dis-rupt the accretion flow and subsequently reduce the accre-tion rate. This could be counteracted if the binary retainsgas that is not bound to either star. In this case, the gasdensity near the individual stars will be increased, poten-tially offsetting the decrease due to orbital velocity.The stability of BHLA has is the topic of some debate,with conclusions depending on dimensionality and assumedsymmetries, as well as physical quantities such as Mach num-ber, M , and accretor size (Foglizzo et al. 2005). A variety ofinstability mechanisms and types are have been proposed,generally sharing the property that at low Mach numbersthe flow is stable, only becoming unstable for M (cid:38) oraccretors with sizes much less than the accretion radius. Be-cause none of our simulations fall in this range, we do notobserve nor further address instability of the flow in thispaper. We simulate hydrodynamics with the GANDALF ∇ h smoothed particle hydrodynamics (SPH) code (Hubber et al.2018) using an M4 cubic spline kernel (Monaghan & Lattan-zio 1985) and a leapfrog kick-drift-kick integration scheme.To accurately simulate energy dissipation in shocks we solvethe energy equation for an ideal gas, using an adiabatic in-dex of 5/3, and without artificial conductivity. Because we expect strong shocks, we use the Monaghan (1997) alpha-viscosity, with a viscosity switch as described in Morris &Monaghan (1997). This prevents particles jumping acrossa shock in a single simulation timestep and ensures thatthe shocks are accurately simulated. We do not include theeffect of self-gravity of the gas, nor the subsequent dynam-ical friction on the stars. This reduces the computationalrequirements, as well as keeping our results independent ofthe gas density, and is in accord with the standard BHLAprescriptions in which the density of the gas is assumed tobe negligible.Barai et al. (2011) show that spherical Bondi accretioncan be accurately modelled using the SPH code gadget-3 ,and at similar resolution to our simulations. They note twounphysical effects that occur as a result of the numericalmethod and its finite resolution. The first is excess heatingof the gas near the accretor as a result of artificial viscosity.This effect should be diminished in our simulations becauseof the use of a viscosity switch, which reduces the artificialviscosity in regions away from shocks. The switch allowsthe artificial viscosity to decay towards a minimum specifiedvalue in regions where the flow’s velocity has positive orzero divergence. For a radial inflow, if the velocity is locallydescribed by a power law v ∝ r − α , the condition for negativedivergence is α > . Our simulations show α is typicallyclose to 1 in the vicinity of the accretor (Appendix A), sothe viscosity should still be its minimum value here. That isnot to say that there will be no heating from viscosity, butthat the amount should be less than observed in Barai et al.(2011).The second unphysical effect produced in their simula-tions is a flow of material across the simulation boundarywhich occurs as a result of their initial setup and boundaryconditions. We avoid this problem by initially filling the en-tire domain of the simulation with uniform fluid and usingperiodic boundary conditions.Our simulations begin with a cuboid of fluid at uni-form density and temperature, and SPH particles in a cubiclattice. The particle separation is linked to the smoothinglength by the parameter η , where, h i = η (cid:18) m i ρ i (cid:19) / , (9)and h i , m i , and ρ i are the smoothing length, mass, and dens-ity of particle i . In all our simulations, we keep η at its defaultvalue of 1.2, meaning that, initially, the smoothing length ofthe particles is 1.2 times their separation. Most simulationsuse approximately 2, 4 or 8 million particles with a particlemass of × − in simulation units as defined in § MNRAS , 1–18 (2019)
T.A.F. Comerford et al. due to its lower numerical viscosity. We perform both sets ofsimulations with the same particle spacing, since this meansthat the SPH and MFV simulations are similarly resolved.
Our simulations span a range of binary separations, inclin-ations, mass ratios and velocities. We do not treat eccentricbinaries in this work because this would require two addi-tional parameters, the eccentricity and longitude of periap-sis. Two main effects limit the ranges of parameters used inthe simulations. The first is that if the binary separation istoo wide, the size of the simulation box needed to preventinteraction with the boundary becomes impractically large.The nature of our simulation setup also precludes simula-tions of very high-velocity binaries, because these will gen-erally leave the simulation domain before reaching a steadystate.The parameter space is also limited by the finite resol-ution of our simulations. As shown in Ruffert (1994), whenthe size of the accretor is large, the bow shock becomes at-tached to the accretor, producing accretion rates higher thanin better-resolved simulations. Overall, this effect results inthe constraint that the accretor size (and so, the smoothinglength) must be smaller than r a by a factor that dependson the accretor’s velocity. Appendix A shows the effects ofvarying resolution on our accretion rates. This imposes anupper limit on the centre-of-mass velocity, and lower lim-its on the binary separation and mass ratio. There are nosuch constraints on the inclination, so we perform simula-tions with inclinations of 0 (face on), 45, and ◦ (side on),to capture the range of behaviour without redundant simu-lations. In addition to the binary simulations, we also carryout a series of single-star runs. These serve as a comparisonto the analytic rates (Eqs. 3–5) and results in other works,and provide a baseline to compare the effect of binary ac-cretors. We do not vary the total mass of the system in oursimulations, so all binary accretors have a combined massequal to that of the single stars.Because stars accrete SPH particles with finite mass,our measured accretion rates are also subject to shot noise.Assuming a Poisson distribution, the expected relative errorin accretion rate is, σ shot = √ N = (cid:118)(cid:117)(cid:117)(cid:116) (cid:16) v + c s (cid:17) / h G M δ t , (10)where N is the number of particles accreted during time δ t over which the accretion rate is measured. For the accretionrate to not be significantly affected by shot noise, we imposethe condition, h (cid:28) r (cid:113) v + c s δ t , (11)which, when δ t is similar to the time taken for an accretorto cross its own BHL radius, reverts to the condition h (cid:28) r a ,which is similar to the condition that ensures that the accre-tion is well-resolved. Since σ shot ∝ ( δ t ) − / , shot noise affectsthe orbit-averaged accretion rates less than the derived in-stantaneous rates.The dimensions of the parameter space are reflected inthe numbering of the simulations in Figures 2 and 3, and in Table B1. The first digit refers to the mass ratio of thesimulated binary, the second to its inclination, the third toits separation, and the fourth digit to the centre-of-massvelocity. Note that the labelling of velocities is not unique,as more than ten different velocities were used in total. From each of our simulations we extract the time-varyingaccretion rate of each accreting star. To calculate the meanaccretion rate we average over a suitable time interval. Insimulations with constant accretion rates (single stars andface-on systems), the time interval over which we averagestarts when the accretion rates reach their steady state, andends at the end of the simulation, before the binary en-counters the simulation boundary. In simulations with time-varying accretion rates, the time interval is truncated at theend so that it is equal to an integral number of binary or-bital periods. This is done to avoid any offset in the averagerate due to variation with orbital phase. Typically, accretionrates reach a steady state after the accretor has travelled adistance equal to its BHL radius, r a , or the sound crossingtime of the BHL radius, in the case of stationary accret-ors with v = . In face-on ( i = ◦ ) systems, the accretionrate remains constant for the duration of the simulations,as one would expect from the symmetry of the system. In-clined systems have accretion rates that vary with time. Themost extreme cases are side-on ( i = ◦ ) orbits for whichthe instantaneous accretion rate varies by up to a factor oftwo during an orbit. Figure 3 shows the accretion rates asa function of orbital phase for a variety of different binaryconfigurations. Our simulations use dimensionless variables and our resultsare converted to physical quantities by multiplying by appro-priate scale factors. A physical quantity, x , can be expressedas the product X x (cid:48) , where X is a scale factor with the sameunits as x , while x (cid:48) is the dimensionless quantity used in thesimulations.We define our simulation units such that the simulationmass unit M = M (cid:12) , and the isothermal sound speed is unity, c s = (cid:112) kT / µ = , where µ is the mass of an H molecule. Asa result, the scale factors depend on the temperature of thegas, T .The resulting scale factors for velocity, length, time andaccretion rate are, V = (cid:18) T K (cid:19) / m s − , (12a) R = . × (cid:18) T K (cid:19) − m , (12b) T = . × (cid:18) T K (cid:19) − / s , (12c) (cid:219) M = . × × (cid:18) T K (cid:19) / kg s − , (12d)respectively. When T = K, the scale factors for length,time and mass accretion rate are approximately
AU, yr, and . × − M (cid:12) / yr, respectively. MNRAS000
AU, yr, and . × − M (cid:12) / yr, respectively. MNRAS000 , 1–18 (2019) inary BHLA In this section we present the results of our simulations byconsidering the initial parameter space one dimension at atime. We start with a general overview of the phenomenaobserved in the simulations ( § § § § § § In each simulation, we observe the typical features expectedof BHLA: a bow shock is formed ahead of the star(s), witha region of higher density gas behind (Figure 2). In thisregion, fluid particle trajectories curve toward the accretor,such that the flow is almost radial near the accretor. Becausethe simulations are carried out in a finite domain, the bowshock interacts with the simulation boundary downstream ofthe accretor(s). This is visible in Figure 2 as the high-densityregions at the top and bottom of each plot.Among supersonic accretors, the interaction with thesimulation boundary has no effect on the accretion rate; thiswas verified by performing identical simulations in widerboxes. However, a slight change in accretion rate was ob-served in the subsonic accretors, even though simulationboxes are always at least twice as large as the character-istic accretion radius. The varying accretion rate appears tobe the result of a periodic change in density, varying witha period equal to the simulation’s sound-crossing time. Thevariation is most extreme in simulations with single accret-ors (i.e. standard Bondi accretion), where the instantaneousaccretion rate varies by about 5 per cent. We find that thevariation is not present in the larger simulation boxes, al-though the mass resolution is subsequently lower, resultingin increased shot noise in the accretion rate.Our higher-resolution simulations, in which accretion isbetter resolved, resulted in an accretion rate around 10 percent higher than our original stationary simulation. Becausethis factor does not appear to depend on binary separation,we apply it as a correction to the accretion rates of all sta-tionary accretors.Our SPH simulations also exhibit a wake of low density,high temperature gas. Over time, the wake collapses intoa series of beads, visually similar to the Plateau-Rayleighinstability (Papageorgiou 1995). The existence of the wakeappears to have little effect on the accretion, except whenone star encounters the other’s wake in side-on systems. Inthese cases, the instantaneous accretion rate drops by up to50 per cent before quickly returning to normal (subplots ( c ) and ( d ) in Fig. 3). However, because the duration of wakeinteraction is only up to 10 per cent of the orbital period,the time-averaged accretion rates are changed only by a fewpercent, which is of similar order to the shot noise in themeasured accretion rates. We performed a small number of simulations using the MFVfluid scheme to determine whether our results are affected by SPH artefacts, such as surface tension (Springel 2010).Because of the more computationally intensive nature of theMFV simulations, we restricted this to the more complexsimulations, namely those in which one star interacts withthe other’s wake.The MFV simulations display very similar flow struc-ture to the SPH, although there are some slight differencesin the density structure between the two simulations (Figure4). The post-shock region in the SPH simulations is approx-imately 10 per cent denser than in MFV, while the low-density wake observed in our SPH simulations is much lesspronounced with MFV, and none of the beading occurs. Weinterpret the beading in the wake as being due to the sur-face tension between SPH particles with differing specificentropy, which artificially prevents mixing. The low-densitywake has little effect on the accretion rates, except when i = ◦ (side-on). In this case, the stars pass through thewake, creating dips in the accretion rate ( § Our simplest set of simulations concern accretion onto asingle sink particle as a test of the standard BHLA form-alism. The relative velocity is ranged between and insimulation units, probing stationary ( v = ), and supersonic( v > ) accretion, as well as the transition region between.The resulting accretion rates (Fig. 5) follow the analyticBHLA rate (Eq. 5) until v = , where the accretion be-comes unresolved. Faster than this, the simulation resolutionis insufficient to resolve the bow shock, and the resulting ac-cretion rate is simply proportional to the accretor velocity.For this reason, we restrict our parameter space to v ≤ in subsequent simulations. From here on, we refer to theseaccretion rates as the ‘single star rates’. The second dimension of the parameter space is the binarysemi-major axis, a , which we vary between and , while theinclination is held at (accretion face-on to the orbit) andthe mass ratio at (equal mass stars, each being half themass of the stars in § at small separations and to / at wide separations, wherethe binary can be treated as a single accretor, and two non-interacting accretors, respectively. However, this does nothold for binaries with v = ; these show a slight increaseover the equivalent single-star rate for close binaries, andthe ratio drops to 1 when a / r a ≥ / . This appears to bethe result of unresolved accretion, where, as (cid:219) M ∝ v , thereincrease accretion in close binaries owing to their high or-bital velocities, which are non-negligible compared to totalvelocity of the binary.The separation dependence of the accretion rate de-pends on the binary’s velocity relative to the gas; in slowerbinaries, the accretion rate drops to half the single-star rate MNRAS , 1–18 (2019)
T.A.F. Comerford et al.
Figure 2.
Sample snapshots from four of our simulation runs, labelled with the simulation parameters (also available in Table B1, towhich the simulation numbers refer). The accretors and their previous trajectories are marked with circles and dashed lines, respectively;the colours correspond to the line colours in Fig 3. The colour shows the gas density, integrated along the z -axis of the simulationdomain and normalised such that the ambient density is . The snapshots are taken at different times in each simulation, in order tobest demonstrate the observed phenomena; see Fig 3 for the orbital phase of each snapshot. at shorter binary separations. This can be explained by con-sidering the total velocity of each accretor. Because, in face-on binaries, the centre-of-mass velocity and orbital velocityof the individual accretors are perpendicular, they are addedin quadrature to obtain the magnitude of the total velocityof each star relative to the gas. As a result, at a given separ-ation, the difference between single and binary accretors ismost pronounced at slow centre-of-mass velocities, and theeffect of duplicity is more pronounced than in fast accretors.In Figure 6, this is indicated by the dotted vertical lines.In the region to the left of the dotted line, the velocity ofeach individual accretor is dominated by the binary orbitalvelocity; to the right, the centre-of-mass velocity dominates.There is also an increase in accretion rate observed invery close ( a (cid:46) . ) binaries, giving an accretion rate upto 20 per cent higher than in our single star simulations.There are two possible explanations for this phenomenon.First, that this is a genuine physical phenomenon and thatthe higher rate could be justified by considering that theaccretors are sweeping out a larger volume per unit time.In effect, the accretors’ motion effectively expands the ac-cretion radius, r a , giving a higher mean accretion rate. Thesecond explanation, which we believe to be more likely, isthat the increase results from unresolved accretion. In less massive, faster accretors, the BHL rate decreases, but thelimit at which accretion is unresolved increases. The com-bination of these two effects enables the unresolved accretionto noticeably affect the accretion rate. The magnitude of thiseffect is difficult to estimate, however, because the motionsof the fluid and accretors are much more complex than thesingle-star case. When comparing our accretion rates at arbitrary inclina-tions to our face-on rates, the dependence of accretion rateon inclination is not particularly strong; the greatest devi-ation from the face-on case is 37 per cent. While our meanaccretion rates are similar at inclinations of ◦ and ◦ , thegenerally larger bars in the right-hand plot of Figure 7 indic-ate that the side-on accretion rates vary more over an orbitalperiod. This is expected because stars in side-on ( i = ◦ )binaries encounter the other star’s bow shock more often,and have a total velocity which varies more extremely thanwhen face-on or at ◦ .From Figure 7, accretion rates are most enhanced inslow, but not stationary, wide binaries. This is perhaps tobe expected, because in an inclined binary stars spend con- MNRAS000
Sample snapshots from four of our simulation runs, labelled with the simulation parameters (also available in Table B1, towhich the simulation numbers refer). The accretors and their previous trajectories are marked with circles and dashed lines, respectively;the colours correspond to the line colours in Fig 3. The colour shows the gas density, integrated along the z -axis of the simulationdomain and normalised such that the ambient density is . The snapshots are taken at different times in each simulation, in order tobest demonstrate the observed phenomena; see Fig 3 for the orbital phase of each snapshot. at shorter binary separations. This can be explained by con-sidering the total velocity of each accretor. Because, in face-on binaries, the centre-of-mass velocity and orbital velocityof the individual accretors are perpendicular, they are addedin quadrature to obtain the magnitude of the total velocityof each star relative to the gas. As a result, at a given separ-ation, the difference between single and binary accretors ismost pronounced at slow centre-of-mass velocities, and theeffect of duplicity is more pronounced than in fast accretors.In Figure 6, this is indicated by the dotted vertical lines.In the region to the left of the dotted line, the velocity ofeach individual accretor is dominated by the binary orbitalvelocity; to the right, the centre-of-mass velocity dominates.There is also an increase in accretion rate observed invery close ( a (cid:46) . ) binaries, giving an accretion rate upto 20 per cent higher than in our single star simulations.There are two possible explanations for this phenomenon.First, that this is a genuine physical phenomenon and thatthe higher rate could be justified by considering that theaccretors are sweeping out a larger volume per unit time.In effect, the accretors’ motion effectively expands the ac-cretion radius, r a , giving a higher mean accretion rate. Thesecond explanation, which we believe to be more likely, isthat the increase results from unresolved accretion. In less massive, faster accretors, the BHL rate decreases, but thelimit at which accretion is unresolved increases. The com-bination of these two effects enables the unresolved accretionto noticeably affect the accretion rate. The magnitude of thiseffect is difficult to estimate, however, because the motionsof the fluid and accretors are much more complex than thesingle-star case. When comparing our accretion rates at arbitrary inclina-tions to our face-on rates, the dependence of accretion rateon inclination is not particularly strong; the greatest devi-ation from the face-on case is 37 per cent. While our meanaccretion rates are similar at inclinations of ◦ and ◦ , thegenerally larger bars in the right-hand plot of Figure 7 indic-ate that the side-on accretion rates vary more over an orbitalperiod. This is expected because stars in side-on ( i = ◦ )binaries encounter the other star’s bow shock more often,and have a total velocity which varies more extremely thanwhen face-on or at ◦ .From Figure 7, accretion rates are most enhanced inslow, but not stationary, wide binaries. This is perhaps tobe expected, because in an inclined binary stars spend con- MNRAS000 , 1–18 (2019) inary BHLA Figure 3.
Accretion rates as a function of orbital phase in thefour simulations shown in Fig 2. Accretion rates are normalised tothe mean single star rate (subplot a ), and the accretion rates forboth stars are plotted individually where applicable. The verticalline shows the orbital phase at which the snapshots in Fig. 2were taken. The shaded segments in plots (c) and (d) indicatethe period when one star is passing through the wake of the otherstar. The colour of the shading indicates which star is currentlypassing through a wake, corresponding to the colours of the linein this figure and Figure 4, and the dots in Figure 2. siderable time moving through each other’s high-densityshocked region. In slow accretors, this high-density regionhas a greater extent because of the large opening angle ofthe bow shock, while in wide binaries, the density enhance-ment due to proximity of the companion is weakest.Compared to the face-on simulations with v = , therealso appears to be an increase in accretion rate in both in-clined binaries. This is due to the spuriously low rates meas-ured for simulations with a = . and a = (Figure 6); whenthe absolute accretion rates are considered instead, theseresults are consistent with the rates at other binary separa-tions. The analytic prescriptions in § q tends to zero.The simulation results support this prediction; we findthat the accretion rates monotonically increase with decreas-ing q , and that the results when q = / are indistinguish-able from the single-star rates. Binaries with q = / haveaccretion rates between 0 and 20 per cent higher than theequal-mass ( q = ) rates; the rates in binaries with q = / are enhanced between 0 and 50 per cent. In the limit ofwide binaries, Eq. 6 predicts increases of 11 and 36 per cent,respectively. The observed increase in accretion above theequal-mass rate was most extreme at slow centre-of-massvelocities and is weakly positively correlated with binaryseparation. The correlation between accretion rate and sep-aration is unexpected, and partially counters the relations in § (cid:219) q ∝ q q − + q . (13)Since (cid:219) q is always negative, this implies that the accretionshould be oligarchic, and that the initially more massivecomponent will eventually dominate the system. We cantest this prediction by measuring the ratio (cid:219) M / (cid:219) M in oursimulations; according to equation 6, this should be equal to q . Larger ratios correspond to more equal division of gasaccreted by the binary, while lower values imply that the ac-cretion is dominated by the larger component. We find thatin all of our simulations, (cid:219) M / (cid:219) M > q ; i.e. that the meas-ured accretion rates are more equal than predicted by theanalytic theory. However, the rate of accretion onto the lessmassive star is never sufficient to make (cid:219) q positive. The evolution of the binary orbit depends on both the massand angular momentum accreted by the binary. While so farour analysis has focussed on mass accretion, the angular mo-mentum accretion rates can similarly be extracted from oursimulations. Because our simulations do not include dynam-ical friction, all changes in the orbital elements are due toaccretion of mass and angular momentum by the stars. Like-wise, we can also measure the change in the orbital elements a and e . While the simulations are too short to determinelong-term evolution, we can compute instantaneous rates ofchange.The angular momentum ( J ) of a circular orbit is givenby, J ∝ (cid:112) M a , (14)which can be manipulated to give the rate of change of or-bital separation, (cid:219) aa = (cid:219) JJ − (cid:219) MM . (15)As a result, the binary orbit either expands or contractsdepending on the angular momentum and mass accretionrates.The quantity d ln (cid:219) J / d ln (cid:219) M is the ratio of specific angu-lar momenta of the accreted material and the binary; whenthis ratio is less than 1.5, the orbit shrinks. Excluding sim-ulations with insufficient resolution, we find typical valuesof this ratio to be in the range -0.1 to 0.2, with a slightnegative correlation with the binary centre-of-mass velocity.Because this ratio is ratio is less than 1.5, the binary orbitshrinks in all our simulations in which the accretion is suffi-ciently resolved. On timescales much longer than the orbitalperiod, the decrease in binary separation enhances the effectof BBHLA (Figure 6).All the binaries in our simulations started on circularorbits, however they gained eccentricity up to ∼ − in thetime before the accretion rates reached a steady state. Once MNRAS , 1–18 (2019)
T.A.F. Comerford et al.
Figure 4.
Comparison between our SPH and MFV simulations. SPH data are shown in the left-hand column, MFV in the centre, andthe difference between SPH and MFV on the right. The top row shows accretion rates and the bottom row shows density snapshots. Allunits and symbols are as in Figures 2 and 3. accretion is steady, we find that the eccentricity varies ap-proximately sinusoidally with a period equal to the binaryorbit. Interestingly, this variation occurs in both face-on andside-on binaries, although the initial eccentricity increaseis smaller in face-on binaries. The time-averaged (cid:219) e in thesteady state is zero, although the duration of our simula-tions is too short to fully determine the secular variation ineccentricity. While we did not perform simulations of ini-tially eccentric binaries, it is perhaps reasonable to expectthat BBHLA should reduce the orbital eccentricity, by ana-logy with BHLA within a binary system in which (cid:219) e / e isgenerally negative (Theuns et al. 1996; Karakas et al. 2000). In this section, we discuss the implications and applicationsof this work, and compare to existing literature.
The primary pre-existing study of BBHLA is Soker (2004)(hereafter S04), in which the author determines that whena binary orbits inside the wind from a tertiary companion,the accreted material forms a disk with sufficient angular momentum to launch jets. It is difficult to draw direct com-parisons between that work, which focussed on angular mo-mentum accretion onto a binary orbiting a companion, andour work, which concerns mass accretion and orbital evol-ution in idealised circumstances. We can, however test theassumptions to determine the extent to which the simpleanalytic model is realistic.The analysis in S04 is predicated on the assumed order-ing of length scales, R (cid:28) w (cid:46) a (cid:28) r a (cid:28) a trip , (16)where R is the radius of one of the stars in the binary, w isthe width of the Bondi-Hoyle accretion column, a is the bin-ary separation ( a in S04), r a is the BHLA radius ( R acc inS04), and a trip the distance to the tertiary companion (whichmay be treated as infinite in our analysis). Our simulationsverify that when a (cid:28) r a , the binary produces a single com-bined bow shock, and that the form of the flow is similarto the single star case outside the region within a distanceabout a of the binary.However, our simulations do not display the narrowcolumn described in Hoyle & Lyttleton (1939); the gas ap-proaches the binary from a wider region. This explains ourresulting eccentricity variation with orbital phase. Duringhalf the orbit the star encounters material moving in the op- MNRAS000
The primary pre-existing study of BBHLA is Soker (2004)(hereafter S04), in which the author determines that whena binary orbits inside the wind from a tertiary companion,the accreted material forms a disk with sufficient angular momentum to launch jets. It is difficult to draw direct com-parisons between that work, which focussed on angular mo-mentum accretion onto a binary orbiting a companion, andour work, which concerns mass accretion and orbital evol-ution in idealised circumstances. We can, however test theassumptions to determine the extent to which the simpleanalytic model is realistic.The analysis in S04 is predicated on the assumed order-ing of length scales, R (cid:28) w (cid:46) a (cid:28) r a (cid:28) a trip , (16)where R is the radius of one of the stars in the binary, w isthe width of the Bondi-Hoyle accretion column, a is the bin-ary separation ( a in S04), r a is the BHLA radius ( R acc inS04), and a trip the distance to the tertiary companion (whichmay be treated as infinite in our analysis). Our simulationsverify that when a (cid:28) r a , the binary produces a single com-bined bow shock, and that the form of the flow is similarto the single star case outside the region within a distanceabout a of the binary.However, our simulations do not display the narrowcolumn described in Hoyle & Lyttleton (1939); the gas ap-proaches the binary from a wider region. This explains ourresulting eccentricity variation with orbital phase. Duringhalf the orbit the star encounters material moving in the op- MNRAS000 , 1–18 (2019) inary BHLA Accretor velocity − − A cc r e t i o n r a t e / B o nd i r a t e Analytic BHL rateResolution limit
Figure 5.
Accretion rates in our single stars versus velocity (insimulation units), normalised to the stationary Bondi accretionrate. The dashed orange line is the analytic BHLA rate (5), andthe green dot-dashed line is the expression for the unresolved ac-cretion rate derived in Ruffert (1994). The error bars show thestandard deviation in our computed accretion rate, but are smal-ler than the points in most of our simulations. See Appendix Afor a more detailed study on the impact of varying simulationresolution on our results. posite direction to its orbital motion, slowing the star andincreasing eccentricity, while during the other half of the or-bit, the accreted material moves in the same direction as thestar, reducing the eccentricity. If, instead, the material ar-rives in a narrow stream, we would observe strongly varyingaccretion rates and rapidly increasing eccentricity.Our work is similar to Lin & Murray (2007) who ex-amine the accretion of gas onto stellar clusters encounteringa slab of gas. They find that the main effect depends onthe ratio of cluster centre-of-mass velocity to the velocitydispersion within the cluster. When the centre-of-mass ve-locity exceeds the velocity dispersion, the stars accrete indi-vidually, and the cluster as a whole retains little gas. Whenthe dispersion is higher, gas is accreted onto the cluster,but the stars move too quickly to accrete much individually.The gas is retained in the cluster, and accreted over time bythe stars. They consider accretion during a single encounterwith a slab of gas; it is less clear what the effect of continuousaccretion would be.The ratio of cluster velocity to velocity dispersion isanalogous to the ratio of binary centre-of-mass velocity andorbital velocity in our work. The first case, where the centre-of-mass velocity is high, corresponds to our wide binaries,where the stars form separate bow shocks. The second case,where the cluster is found to retain gas, does not appearto have an equivalent in our simulations; we have no caseswhere gas is retained by the binary without being quicklyaccreted by either star. There are several possible reasonswhy we may not have seen this: the first is that many pointmasses may be required to retain gas, while we only havetwo. Another possibility is that our simulations lack suf-ficient resolution to produce circumbinary or circumstellardisks. In this case, formation of a circumbinary disk wouldbe equivalent to the cluster retaining gas after the encounter
Figure 6.
Accretion rates in all our face-on twin simulationsat four different velocities. Each rate is normalised to the corres-ponding single-star ( a = ) rate. The horizontal dashed lines showan accretion rate ratio of . , which is the expected value at largeseparation. The vertical dotted lines show the binary separationat which the orbital velocity equals the centre-of-mass velocity.To the left, orbital velocity dominates the total velocity of eachaccretor; to the right, centre-of-mass velocity dominates. Figure 7.
Accretion rates in all our inclined twin simulations:each rate is normalised to the face-on accretion rate. The leftsubplot includes all simulations with inclinations of ◦ , and theright, simulations with i = ◦ . The x -axis plots the binary separ-ation, while the colour shows the centre-of-mass velocity in sim-ulation units. As before, the bars show the standard deviation inaccretion rate over one orbit. with the slab, without gas becoming associated with the in-dividual cluster stars.A pair of recent papers, Kaaz et al. (2019) and Ant-oni et al. (2019), consider BHLA onto clusters and binarystars, respectively. They reach the same broad conclusionsas us, namely that in systems with a (cid:29) r a the stars accreteseparately, while at closer separations the system accretesgas more like a single body. While we have used SPH andMFV, they use an adaptive mesh refinement (AMR) fluidsimulation, so the fact that we reach similar conclusions im-plies that our results are largely free of spurious simulationartefacts. The simulations in Antoni et al. (2019) are bothperformed at a higher resolution and run for a longer time MNRAS , 1–18 (2019) T.A.F. Comerford et al.
Figure 8.
Derivatives of properties of the binary orbit (mass m , semi-major axis a , eccentricity e , and angular momentum J ), plottedover one orbital period, for a face-on binary (left) and a side-on binary (right). The points and bars to the right of each plot show theorbital averages of each derivative, as well as its standard deviation. than our simulations, allowing the authors to investigate thelong term evolution of the binary orbit. In all binaries thatthey simulate, they find that (cid:219) v / v < (cid:219) a / a < (cid:219) m / m ; that thebinary should come to a halt relative to the gas before theorbit substantially shrinks or the total mass substantiallyincreases. Our simulations agree with the ordering and ap-proximate ratio of (cid:219) a / a and (cid:219) m / m , however, we do not calcu-late (cid:219) v / v , due to the shorter duration of our simulations, andlack of gas drag on the stars. In the astrophysical contextsto which we plan to apply these results, there are additionalsources of momentum which can re-accelerate the binaryafter its centre-of-mass velocity has decreased. In a cluster,this acceleration could be due to a multi-body interactionwhich does not disrupt the binary; in a triple star system,the lost momentum can be regained from the orbit of thebinary about the tertiary companion. We have shown that by being associated in a binary, twostars may accrete up to twice as much material as if theywere isolated. This has implications for several astrophysicalsituations. In young globular clusters, stars may accrete gasleft over from star formation, or ejected by intermediate-mass stars during their AGB phase. This gas is a possiblesource of α -rich elements observed in the spectra of “secondgeneration” stars (Bastian & Lardo 2018). Since binary stel-lar systems are, on average, more massive than single stars,mass segregation results in binaries accumulating towardsthe centre of the cluster, in close proximity to the massiveAGB stars and residual gas. In these regions, the rate ofthree-body encounters is highest, so it is unlikely that thebinary system remains bound for long (Ivanova et al. 2005).Whether the woulds last long enough to substantially af-fect the accretion of material onto its constituent stars isdependent on the properties of the cluster.A similar case can be made for stars near the Galacticcentre, with similar caveats regarding the rate of multiple- body encounters which could unbind the binary. BHLA ontobinaries orbiting within a gas disk affects its orbit aroundthe Galactic centre (Baruteau et al. 2011), and accretiononto black hole binaries may substantially alter their massand prompt a gravitational-wave-producing merger (Yi et al.2018).The final case in which we consider the role of binaryBHLA is in wind mass transfer in hierarchical triple stars,as described by Soker (2004). In this case, a companion toa binary produces a wind which encounters the binary. Theassumptions we make in this paper are most valid in the casewhere the system is very hierarchical, i.e. that the binaryorbital separation is much less than the distance to the thirdstar. We also assume that incoming gas is uniform, withno velocity or density gradients. This assumption not onlyrequires the large separation mentioned above, but also thatthe wind is supersonic and only weakly focussed through thecentral Lagrangian point. While our simulations are sufficient to determine generaltrends, this study has been primarily limited by the timeand computing power requirements. We elected to study awide region of parameter space with relatively low-resolutionsimulations, meaning that there are cases where the phenom-ena we were looking for became unresolved, and that thesimulations could not run for long enough to determine thelong-term evolution. For systems which appeared to be onlymarginally resolved, an increase in simulation resolution byaround a factor of two would be sufficient to fully resolvethe bow shock, although this increases the computationalcomplexity by a factor of at least eight. One important phe-nomenon we have neglected is the formation of disks, botharound the binary and around the individual stars. Our sim-ulations lacked the resolution to accurately form accretiondisks, although their formation is expected due to the an-gular momentum of the infalling material. As a result, we
MNRAS000
MNRAS000 , 1–18 (2019) inary BHLA have been unable to draw definitive conclusions regardingthe angular momentum accretion onto the individual stars.We have also not accounted for radiative cooling, mo-lecular dissociation or ionisation. These processes are sinksof thermal energy and alter the properties of the accretingmaterial, but complicate the modelling due to their depend-ence on the gas density and properties of the star.Similarly, a more complete treatment of this problemcould include the effects of a magnetic field. Unsurprisingly,the addition of magnetic fields has the potential to com-plicate the simulations, whether the field originates in theinterstellar medium, from one or both of the accreting stars,or the donor star in a triple system. In Lee et al. (2014), theauthors perform simulations of BHLA of molecular gas witha magnetic field, finding that in the case of Bondi accretion,the accretion rate is significantly reduced when the mag-netic pressure exceeds one per cent of the gas pressure. Formoving accretors, the critical magnetic pressure increasesroughly proportionally to the accretor’s Mach number, somagnetic fields will have less of an effect on fast accretors.Applying the criteria in Lee et al. (2014), the conditionon the ambient magnetic field to not disrupt the accretionis, B (cid:46) . × − T (cid:18) T
100 K (cid:19) / (cid:18) ρ − kg m − (cid:19) / ( + M ) / , (17)where M is the accretor’s Mach number.To determine whether BHLA can be magnetically in-hibited in a multiple stellar system, we can examine the ex-ample case of a companion accreting from the wind of a giantstar with a radius of R = R (cid:12) and a magnetic field witha surface strength of .
01 T corresponding to the strongestmagnetic fields measured in G/K giants by Auri`ere et al.(2015).The field strength a distance r from the giant is, assum-ing a dipole field geometry, approximately, B = .
01 T (cid:18) R (cid:12) r (cid:19) =
10 T R (cid:12) r , (18)hence at a minimum distance to the accreting star of r = R = R (cid:12) then B = / × − T = . × − T . Adopting T = and ρ = − kg m − in the vicinity of the ac-cretor, Equation 17 gives a threshold magnetic field strengthof at least . ( + M ) T . As a result, we would not expectBHLA to be strongly affected by the presence of the mag-netic field in a typical red-giant wind accretion scenario, evenfor low Mach numbers.In realistic physical scenarios, the accreted material willhave some angular momentum, so formation of a disk aroundeither or both of the accretors. In these cases, the interactionof the disk with magnetic fields can reduce the accretion rateand produce outflows (Mo´scibrodzka & Proga 2009).There are many directions in which one could expandon this work. In addition to wider and finer sampling ofthe parameter space, and higher resolution simulations, onecould also introduce more varied physical context. One pos-sible route would be to abandon the assumption of a homo-geneous medium, and allow the incoming material to have variations in density and temperature (in directions bothparallel and perpendicular to the relative velocity; MacLeod& Ramirez-Ruiz 2015). Introducing a longitudinal densitygradient and angular momentum would produce a modelthat better represents accretion in a triple stellar system,while a transverse gradient would represent the model inSoker (2016), in which an asymptotic giant branch star swal-lows a binary that merges inside the envelope of the largerstar. Using the code gandalf , we have performed a suite of sim-ulations of Bondi-Hoyle-Lyttleton accretion (BHLA) ontosingle and binary accretors, with a range of separations, ve-locities, inclinations and mass ratios, using SPH and MFVschemes. We recover the analytic form of the BHLA ratein single stars, and find the expected result that in veryclose binaries the accretion rate approaches the single-starrate. In wide binaries the mass accretion rate drops to halfthe single-star rate, where in slow-moving binaries, the ratereaches the asymptote at lower binary separations.We find that the inclination of the binary relative tothe centre-of-mass velocity vector has comparatively littleeffect on the mass accretion rates but is important whenconsidering the orbital evolution of the binary because it al-ters the angular momentum accretion rate. In binaries wherethe stars have unequal masses, the accretion rate of the lessmassive star is enhanced by its proximity to the companion,and the more massive star accretes at a disproportionatelyhigher rate, leading to a decrease in the mass ratio, q . If thistrend continues over timescales in which the stars’ massesincrease substantially, it implies that binaries should evolveaway from equal mass ratios toward systems that are dom-inated by a single component.While this initial study was limited by available resolu-tion and computation time, we have demonstrated the gen-eral trends in accretion rates as a function of binary separa-tion, velocity, inclination, and mass ratio. Simulations withhigher spatial resolution would demonstrate the behaviourof the gas in the immediate vicinity of the stars, includingthe potential formation of circumstellar and circumbinarydiscs. Simulating many binary orbits would also constrainthe long-term orbital evolution. Future work could involvedetermining the role of eccentricity in binary BHLA, as wellas more physically realistic simulations of the applicationsin globular clusters, dense stellar environments and multiplestellar systems. ACKNOWLEDGEMENTS
MNRAS , 1–18 (2019) T.A.F. Comerford et al.
The Data Centric System was funded by a BIS Na-tional E-infrastructure capital grant ST/K00042X/1, STFCcapital grant ST/K00087X/1, DiRAC Operations grantST/K003267/1 and Durham University. DiRAC is partof the National e-Infrastructure. RGI thanks the STFCfor Rutherford fellowship funding under grant numberST/L003910/1. Hello to Jason Isaacs. G.R. acknowledgessupport from the Netherlands Organisation for Scientific Re-search (NWO, program number 016.Veni.192.233).Figures were produced using Matplotlib (Hunter 2007),SPLASH (Price 2007), and Inkscape ( https://inkscape.org ). References
Abate C., Stancliffe R. J., Liu Z.-W., 2016, A&A, 587, A50Antoni A., MacLeod M., Ramirez-Ruiz E., 2019, ApJ (inpress), arXiv 1901.07572,Auri`ere M., et al., 2015, A&A, 574, A90Barai P., Proga D., Nagamine K., 2011, MNRAS, 418, 591Baruteau C., Cuadra J., Lin D. N. C., 2011, ApJ, 726, 28Bastian N., Lardo C., 2018, ARA&A, 56, 83Boffin H. M. J., Jorissen A., 1988, A&A, 205, 155Bondi H., 1952, MNRAS, 112, 195Bondi H., Hoyle F., 1944, MNRAS, 104, 273Edgar R., 2004, New Astron. Rev., 48, 843Foglizzo T., Galletti P., Ruffert M., 2005, A&A, 435, 397Hopkins P. F., 2015, MNRAS, 450, 53Hoyle F., Lyttleton R. A., 1939, Proceedings of the CambridgePhilosophical Society, 35, 405Hubber D. A., Rosotti G. P., Booth R. A., 2018, MNRAS, 473,1603Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Ivanova N., Belczynski K., Fregeau J. M., Rasio F. A., 2005,MNRAS, 358, 572Kaaz N., Antoni A., Ramirez-Ruiz E., 2019, arXiv 1901.03649,Karakas A. I., Tout C. A., Lattanzio J. C., 2000, MNRAS, 316,689Lanson N., Vila J.-P., 2008, Siam. J. Numer. Anal., 46, 1912Lee A. T., Cunningham A. J., McKee C. F., Klein R. I., 2014,ApJ, 783, 50Lin D. N. C., Murray S. D., 2007, ApJ, 661, 779MacLeod M., Ramirez-Ruiz E., 2015, ApJ, 803, 41Moeckel N., Throop H. B., 2009, ApJ, 707, 268Monaghan J. J., 1997, Journal of Computational Physics, 136,298Monaghan J. J., Lattanzio J. C., 1985, A&A, 149, 135Morris J. P., Monaghan J. J., 1997, Journal of ComputationalPhysics, 136, 41Mo´scibrodzka M., Proga D., 2009, MNRAS, 397, 2087Papageorgiou D. T., 1995, Physics of Fluids, 7, 1529Price D. J., 2007, Publ. Astron. Soc. Australia, 24, 159Ruffert M., 1994, ApJ, 427, 342Sakelliou I., 2000, Monthly Notices of the Royal AstronomicalSociety, 318, 1164Shima E., Matsuda T., Takeda H., Sawada K., 1985, MNRAS,217, 367Soker N., 2004, MNRAS, 350, 1366Soker N., 2016, MNRAS, 455, 1584Springel V., 2010, ARA&A, 48, 391Theuns T., Boffin H. M. J., Jorissen A., 1996, MNRAS, 280, 1264Thompson I. B., et al., 2008, The Astrophysical Journal, 677, 556Thoul A., Jorissen A., Goriely S., Jehin E., Magain P., Noels A.,Parmentier G., 2002, A&A, 383, 491Yi S.-X., Cheng K. S., Taam R. E., 2018, ApJ, 859, L25 MNRAS000
Abate C., Stancliffe R. J., Liu Z.-W., 2016, A&A, 587, A50Antoni A., MacLeod M., Ramirez-Ruiz E., 2019, ApJ (inpress), arXiv 1901.07572,Auri`ere M., et al., 2015, A&A, 574, A90Barai P., Proga D., Nagamine K., 2011, MNRAS, 418, 591Baruteau C., Cuadra J., Lin D. N. C., 2011, ApJ, 726, 28Bastian N., Lardo C., 2018, ARA&A, 56, 83Boffin H. M. J., Jorissen A., 1988, A&A, 205, 155Bondi H., 1952, MNRAS, 112, 195Bondi H., Hoyle F., 1944, MNRAS, 104, 273Edgar R., 2004, New Astron. Rev., 48, 843Foglizzo T., Galletti P., Ruffert M., 2005, A&A, 435, 397Hopkins P. F., 2015, MNRAS, 450, 53Hoyle F., Lyttleton R. A., 1939, Proceedings of the CambridgePhilosophical Society, 35, 405Hubber D. A., Rosotti G. P., Booth R. A., 2018, MNRAS, 473,1603Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Ivanova N., Belczynski K., Fregeau J. M., Rasio F. A., 2005,MNRAS, 358, 572Kaaz N., Antoni A., Ramirez-Ruiz E., 2019, arXiv 1901.03649,Karakas A. I., Tout C. A., Lattanzio J. C., 2000, MNRAS, 316,689Lanson N., Vila J.-P., 2008, Siam. J. Numer. Anal., 46, 1912Lee A. T., Cunningham A. J., McKee C. F., Klein R. I., 2014,ApJ, 783, 50Lin D. N. C., Murray S. D., 2007, ApJ, 661, 779MacLeod M., Ramirez-Ruiz E., 2015, ApJ, 803, 41Moeckel N., Throop H. B., 2009, ApJ, 707, 268Monaghan J. J., 1997, Journal of Computational Physics, 136,298Monaghan J. J., Lattanzio J. C., 1985, A&A, 149, 135Morris J. P., Monaghan J. J., 1997, Journal of ComputationalPhysics, 136, 41Mo´scibrodzka M., Proga D., 2009, MNRAS, 397, 2087Papageorgiou D. T., 1995, Physics of Fluids, 7, 1529Price D. J., 2007, Publ. Astron. Soc. Australia, 24, 159Ruffert M., 1994, ApJ, 427, 342Sakelliou I., 2000, Monthly Notices of the Royal AstronomicalSociety, 318, 1164Shima E., Matsuda T., Takeda H., Sawada K., 1985, MNRAS,217, 367Soker N., 2004, MNRAS, 350, 1366Soker N., 2016, MNRAS, 455, 1584Springel V., 2010, ARA&A, 48, 391Theuns T., Boffin H. M. J., Jorissen A., 1996, MNRAS, 280, 1264Thompson I. B., et al., 2008, The Astrophysical Journal, 677, 556Thoul A., Jorissen A., Goriely S., Jehin E., Magain P., Noels A.,Parmentier G., 2002, A&A, 383, 491Yi S.-X., Cheng K. S., Taam R. E., 2018, ApJ, 859, L25 MNRAS000 , 1–18 (2019) inary BHLA APPENDIX A: RESOLUTION STUDY
To ensure that our simulations are adequately resolved, weperformed a resolution study consisting of a series of simula-tions with varying resolution and accretor velocity, with theaim of comparing our results with those of Ruffert (1994,hereafter R94).In R94 and its subsequent papers, the authors derive anexpression for the BHL accretion rate (their Eq. 15) whichinterpolates between cases when the flow is well- and poorly-resolved. They then go on to test their analytic rates usinga set of adaptive mesh refinement (AMR) simulations, andfind them to be in good agreement.Qualitatively, they find that when the accretion is well-resolved, the bow shock is separated from the accretor, andthere is little to no resolution-dependence in the accretionrates. However for larger accretors, the shock cone becomesattached to the sink particle, yielding unreliable accretionrates that increase proportionally with the square of the ac-cretor radius. The transition between these regimes dependson the velocity of the accretor velocity, with faster-movingaccretors becoming unresolved at smaller physical scales.If the accretor size is held at about . Bondi radii,the default in our main simulations, the results of R94 implythat accretion is unresolved at velocities of 6 or higher, whichagrees well with our results (Fig. 5).The simulations in this resolution study span a rangeof accretor velocities between Mach 1 and 4, and accretorsizes between . and . stationary Bondi radii. At higherresolutions than this, the computational cost increases rap-idly, while at lower resolutions, our simulation fail to con-verge. The results of our simulations match those of R94well (Fig. A1); at the three different accretor velocities weobserve both the flat, well-resolved region and the slope upto the unresolved rate for large accretors, with the transitionpoint at the size predicted in R94.Our higher-resolution accretion rates generally fall be-low the analytic predictions by between 10 and 20 percent,which is consistent with other simulations of BHLA (Edgar2004).We also compare the results of our simulation 0009 tothe expected density and velocity profiles of spherical Bondiaccretion (Figure A2). The analytic treatment in Bondi(1952) shows that for an adiabatic index of / , the fluiddensity and velocity increase without bound near the ac-cretor, but the flow never becomes supersonic. Our simula-tions do show a sonic surface within one to two smoothinglengths of the sink, which deviates from the theory, but doesensure that the conditions immediately adjacent do not af-fect the extended flow. The cause of the extra accelerationis the artificial pressure gradient created when particles areremoved from the simulation during accretion. APPENDIX B: LIST OF SIMULATIONPARAMETERS
Figure A1.
Our version of Fig. 6 in Ruffert (1994), showing thedependence of the accretion rates on the physical size of the ac-cretor (normalised to the Bondi rate and radius respectively). Thethree colours correspond to different accretor velocities, and thedotted lines are Eq. 15 from R94. Points marked with stars indic-ate the standard resolution used in our single-star simulations.MNRAS , 1–18 (2019) T.A.F. Comerford et al. − − Radius D e n s i t y − − Radius . . . . . . V e l o c i t y Figure A2.
Density and velocity profiles of spherical Bondi accretion for simulation 0009. The grey region to the left represents theradius of the sink particle; analytic profiles from Bondi (1952) are also plotted in grey. Density and velocity are normalised to the ambientdensity and sound speed; the radius is in simulation units, in which the sink particle has a radius of 0.039 and the Bondi radius (2) is0.6. MNRAS000
Density and velocity profiles of spherical Bondi accretion for simulation 0009. The grey region to the left represents theradius of the sink particle; analytic profiles from Bondi (1952) are also plotted in grey. Density and velocity are normalised to the ambientdensity and sound speed; the radius is in simulation units, in which the sink particle has a radius of 0.039 and the Bondi radius (2) is0.6. MNRAS000 , 1–18 (2019) inary BHLA Table B1: List of our SPH simulation runs. ‘Resolution’ refers to thesize of the simulation domain in units of the SPH particle spacing while‘dimension’ is the absolute size in simulation units. x , y , and z specifythe quantity along the corresponding coordinate axes; the y and z valuesare equal in all cases. The next four columns list the binary mass ratio( q ), inclination ( i ), semi-major axis ( a ), and centre-of-mass velocity ( v ).Generally, the four digits of the simulation ID correspond to the valuesof q , i , a , and v , respectively. There are some exceptions to this, as morethan 10 different velocities were used. The final three columns containthe total binary accretion rate, normalised by the ambient gas density( (cid:219) M / ρ ); ratio of accretion rates of each of the stars ( (cid:219) M / (cid:219) M , where M is the more massive star); and the fractional change of semi-major axis,normalised by gas density ( (cid:219) a /( ρ a ) ). Simulations marked with a dagger( † ) were repeated using an MFV simulation.ID Resolution Dimension / q i ( ◦ ) a v (cid:219) M / ρ (cid:219) M / (cid:219) M (cid:219) a /( ρ a ) x y , z x y , z †
512 128 16 4 - - - 2 0.356 - -0002 512 128 16 4 - - - 4 0.0683 - -0004 256 128 8 4 - - - 0.5 0.979 - -0005 512 128 16 4 - - - 3 0.118 - -0006 512 128 16 4 - - - 5 0.0389 - -0007 512 96 16 3 - - - 7 0.0354 - -0008 512 96 16 3 - - - 6 0.0316 - -0009 128 128 4 4 - - - 0 0.919 - -0009b 203 203 8 8 - - - 0 1.008 - -0009c 322 322 16 16 - - - 0 1.090 - -0010 256 128 8 4 1 90 1 1 0.569 0.999 -1.140011 512 128 16 4 1 90 1 2 0.167 1.000 -0.310012 512 128 16 4 1 90 1 4 0.0666 0.997 -0.220014 512 128 16 4 1 90 1 2.5 0.108 1.000 -0.30015 512 128 16 4 1 90 1 1.5 0.369 1.000 -0.520019 128 128 4 4 1 - 1 0 0.546 1.000 -0.970020 256 128 8 4 1 90 1/2 1 0.616 1.000 -2.550021 512 128 16 4 1 90 1/2 2 0.248 1.000 -0.780022 512 128 16 4 1 90 1/2 4 0.0654 1.000 -0.450029 128 128 4 4 1 - 1/2 0 0.580 1.000 -2.190030 256 128 8 4 1 90 1/4 1 0.685 1.000 -6.570031 512 128 16 4 1 90 1/4 2 0.335 1.000 -3.770032 512 128 16 4 1 90 1/4 4 0.0674 1.010 -0.920039 128 128 4 4 1 - 1/4 0 0.629 1.000 -5.510040 256 128 8 4 1 90 1/8 1 0.786 1.000 -17.810041 512 128 16 4 1 90 1/8 2 0.359 1.000 -8.470042 512 128 16 4 1 90 1/8 4 0.0729 1.000 -1.920049 128 128 4 4 1 - 1/8 0 0.759 1.000 -16.290050 256 128 8 4 1 90 1/16 1 0.953 1.000 -49.140051 512 128 16 4 1 90 1/16 2 0.401 1.000 -22.680052 512 128 16 4 1 90 1/16 4 0.082 1.000 7.960059 128 128 4 4 1 - 1/16 0 0.987 1.000 -50.280070 256 128 8 4 1 90 0.7 1 0.604 1.000 -1.850071 512 128 16 4 1 90 0.7 2 0.194 0.999 -0.470072 512 128 16 4 1 90 0.7 4 0.068 1.000 -0.330079 128 128 4 4 1 - 0.7 0 0.567 1.000 -1.440110 256 128 8 4 1 0 1 1 0.532 0.911 -0.790111 512 128 16 4 1 0 1 2 0.230 0.681 -0.510112 512 128 16 4 1 0 1 4 0.0631 1.000 -0.230114 512 128 16 4 1 0 1 2.5 0.135 0.718 -0.390115 512 128 16 4 1 0 1 1.5 0.375 0.967 -0.450120 512 128 16 4 1 0 1/2 1 0.566 0.963 -2.560121 †
512 128 16 4 1 0 1/2 2 0.240 0.944 -0.950122 512 128 16 4 1 0 1/2 4 0.0691 1.030 -0.46
MNRAS , 1–18 (2019) T.A.F. Comerford et al. †
512 128 16 4 4 0 1/2 2 0.300 5.920 -1.132122 512 128 16 4 4 0 1/2 4 0.0759 2.140 -0.722130 256 128 8 4 4 0 1/4 1 0.746 4.530 -7.082131 512 128 16 4 4 0 1/4 2 0.328 4.690 -3.122132 512 128 16 4 4 0 1/4 4 0.0792 2.050 -1.62210 256 128 8 4 4 45 1 1 0.664 6.730 -1.22
MNRAS000
MNRAS000 , 1–18 (2019) inary BHLA MNRAS , 1–18 (2019) T.A.F. Comerford et al.
This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000