The evolution of the mass ratio of accreting binaries: the role of gas temperature

Abstract

1 INTRODUCTION

The accretion of gas on to binary systems is of widespread astrophysical importance. For supermassive black hole binaries in galactic nuclei, the details of gas accretion is important in setting the observational signatures of such binaries (Artymowicz & Lubow 1996; MacFadyen & Milosavljević 2008) and in determining the alignment of black hole spins (Bogdanović, Reynolds & Miller 2007; Lodato & Gerosa 2013).

Another important application is to the formation of stellar mass binaries. This is a central problem for understanding star formation since a large fraction of stars reside in binary systems (Duquennoy & Mayor 1991; Fischer & Marcy 1992; Mason et al. 1998; Preibisch et al. 1999; Close et al. 2003; Basri & Reiners 2006). It is well known that the initial mass of a protobinary forming from a collapsing cloud is a small fraction of the mass of the parent cloud (except for the widest binaries). The initial protobinary forms from the first material to collapse, which is located near the clouds centre and thus has low specific angular momentum. The protobinary then acquires most of its mass from the remaining infalling material, which has greater specific angular momentum. Because of its greater specific angular momentum, this material is expected to first form a circumbinary disc. The binary then accretes material through inflowing streams linked to this circumbinary disc. The ultimate mass ratio of the binary (q = M₂/M₁, where M₂ and M₁ are the masses of the secondary and primary components, respectively) is thus set by the relative accretion rates from the circumbinary disc on to the primary and secondary components, not by the mass ratio of the initial protobinary fragment.

A number of numerical studies have focused on understanding accretion on to isolated binaries. These studies have ranged from ballistic calculations (Artymowicz 1983; Bate 1997) to hydrodynamical calculations using smoothed particle hydrodynamics (SPH) in 3D (Bate & Bonnell 1997; Bate 2000) and grid codes in 2D (Ochi, Sugimoto & Hanawa 2005; Hanawa, Ochi & Ando 2010; de Val-Borro et al. 2011; Takakuwa et al. 2014). Simulations of entire clusters have also been used to predict ensemble properties of binaries (Goodwin, Whitworth & Ward-Thompson 2004a, 2004b; Delgado-Donate & Clarke 2008; Bate 2009; Offner, Hansen & Krumholz 2009). The resolution of the simulations of entire clusters is inevitably poorer than in simulations dedicated to a single isolated binary.

SPH (and also ballistic) simulations of circumbinary accretion have always found a strong tendency for the majority of the accreted mass to flow on to the secondary. This is because material inflowing from the circumbinary disc encounters the secondary first, since it orbits further from the centre of mass of the binary than the primary. Therefore, SPH studies of binaries at both the individual binary and cluster scale exhibit a preference for creating binaries with high q values, since material falling preferentially on the secondary implies |$\dot{q}>0$|⁠. This preference is in mild conflict with the observed statistics of solar mass binaries, which show a flat q distribution (Raghavan et al. 2010).

A re-examination of the accretion of gas on to binaries using 2D grid-based techniques by Ochi et al. (2005) made the opposite prediction, finding more mass accreted by the primary. Although these simulations also showed that material encounters the secondary first, they found a significant fraction of it continued to flow past the secondary and became bound to the primary. Unless q ≪ 1 this results in a decreasing mass ratio (i.e. |$\dot{q}<0$|⁠). Higher resolution grid studies (Hanawa et al. 2010) later confirmed that this effect was not a result of insufficient resolution. Ochi et al. (2005) hypothesized that the discrepancy with SPH simulations was due to excessive numerical dissipation in the SPH code causing accreting particles to spiral into the secondary instead of transitioning to the primary.

However, as well as using different numerical techniques, the SPH and grid-based studies of accretion on to isolated binaries also differ in the temperature used for the accreting gas. Because the dynamical importance of pressure to the flow of material on to the binary is a function of temperature, it is possible for gas temperature to qualitatively alter the accretion dynamics. It has been suggested that the discrepancy between SPH and grid-based simulations is due not to numerical technique but to this difference in temperature (Clarke 2008).

Here, we explore the hypothesis that |$\dot{q}$| depends on the temperature of the accreting gas, using 2D SPH simulations at a variety of resolutions and compare our results with SPH and grid-based results in the literature. We also aim to discover the resolution required in order to accurately determine the mass flow rates on to the primary and secondary components of a binary (of value to cluster simulations). Section 2 describes the simulations performed, Section 2 reports on the results and Section 4 discusses the application to astrophysical systems. Section 5 presents our Conclusions.

2 BINARY MODEL

We model a binary on an anticlockwise circular orbit, with mass ratio q = M₂/M₁, total mass M = M₁ + M₂ and separation a. We use a modified version of the SPH code gadget2 (Springel 2005) in 2D. We have chosen 2D SPH in order to facilitate the comparison with the 2D grid simulations of Ochi et al. (2005), Hanawa et al. (2010) and since numerical convergence is faster in 2D than 3D (for N particles, the SPH smoothing length-scales as N^−1/2 in 2D and N^−1/3 in 3D). In a real circumbinary disc, accretion is driven by the effective viscosity in the disc (Balbus & Hawley 1991; Boss 2000; Lodato 2008; Turner et al. 2014) which feeds the binary with material that has similar specific angular momentum to that of the binary itself, even if the material falling into the disc at large radius has a significantly higher specific angular momentum. In order to avoid the computational expense of modelling this situation (with numerical or some prescribed ‘physical’ viscosity driving the secular evolution of the circumbinary disc) we instead follow previous authors and inject gas at large radius R_infa with a ratio (j_inf) of specific angular momentum to that of the binary (⁠|$\sqrt{GMa}$|⁠) of order unity (see Table 1). We set the initial inward radial velocity of the injected particles so that these particles are initially marginally gravitationally bound to the binary. The stars are treated as sink particles so that particles are accreted if they stray within the accretion radius of each star (0.01a). We set the SPH smoothing length h using |$h=1.2\sqrt{\frac{m}{\Sigma }}$|⁠, where m is the particle mass and Σ is the surface density (Price 2012). This choice gives roughly 18 neighbours per SPH particle in 2D. Details of how to access the exact version of the code used, along with parameter files and initial conditions can be found in Section 6.

Note that the number of injected particles is higher for the ‘hot’ simulations (c = 0.25: 5–8 in Table 1), since a significant fraction of the marginally bound injected particles are pushed out by pressure gradients. These particles are deleted from the simulation soon thereafter. The ‘cold’ simulations (with c = 0.05: 1–4 in Table 1) are each matched with a corresponding hot simulation (5–8) so that they have roughly the same number of particles actively contributing to the simulation at any time (see Table 1 for the corresponding values of |$\dot{N}$|⁠).

For each simulation, we also measure h/H (SPH smoothing length over disc scaleheight) at the outer edge of the circumsecondary disc. This quantity measures how well resolved the secondary's disc is at the point where it interacts with infalling material. As expected, h/H scales as N^−1/2, except for the lowest resolution simulations where it increase more rapidly with declining particle number. This rapid increase occurs when the role of numerical viscosity is enhanced to the point that the accretion time-scale of the circumsecondary disc becomes comparable with the infall time-scale.

Figure 1.

Equipotential surfaces of Φ, the modified potential (gravitational potential plus term to include centrifugal force of non-inertial frame) for the values of Φ obtained at the three Lagrange points L1,L2 and L3. The colour map shows Φ elsewhere. The stars are marked in red (secondary) and green (primary) and the centre of mass of the system is marked in blue. The green and blue lines are the locations where the x and y components of the force resulting from this potential (i.e. −∇Φ) are zero.

Open in new tabDownload slide

3 RESULTS

3.1 Dependence of mass ratio evolution on temperature

Figure 2.

Accretion rates and corresponding qΓ (equation 3) for cold (left) and hot (right) simulations. Red lines represent qΓ, blue is accretion on to the primary and green is accretion on to the secondary (both per unit accreted mass). The solid, dashed and dot–dashed lines represent h/H = 0.04, 0.1 and 0.25 (⁠|$\dot{N} = 1000,200$| and 100), respectively, for the cold simulation and h/H = 0.2, 0.45 and 0.75 (⁠|$\dot{N}=2000,500$| and 200), respectively, for the hot simulation. The first ∼30 orbits are not shown as the simulation has yet to settle into a quasi-steady state.

Open in new tabDownload slide

It is clear that significantly more material is accreted by the primary in the hot simulation than in the cold. Averaging across the steady state, this translates into a higher value of Γ for the cold simulation (where Γ ≈ 5) compared with the hot simulation (Γ ≈ 3.5). We note that for the cold simulation, this value is very similar to that found in the previous 3D SPH simulation of Bate & Bonnell (1997) for the same parameters, implying that the dimensionality of the problem is not critical here.

Our result for the hot simulation is close to, but somewhat larger than, the value (Γ ∼ 2.5) obtained by Ochi et al. (2005) for their corresponding model (2-12). However, the value reported by Ochi et al. is obtained by averaging over a short time period (t ∼ 30 in the units of Fig. 2) compared to the period of the variability (t ∼ 90). If we average our results on a similarly short time period, Γ can take values between 2.5 and 4 depending on when the average is taken. Furthermore, the total duration of the Ochi et al. simulations is shorter than the period we have excluded from our analysis on account of it not being in a steady state. As such, the value of Γ reported by Ochi et al. may still be partially measuring transient accretion driven by their initial conditions. Finally, Ochi et al. use a different definition of accretion (counting mass within the entire disc) which we found to produce greater variability. Given these differences, our results do not suggest a large discrepancy between grid-based and SPH results.

Although we have limited our analysis to simulations that have already settled to a quasi-steady state, there is still some time variability in the accretion rates (particularly in the hot simulations) as has been identified by previous authors (Ochi et al. 2005; MacFadyen & Milosavljević 2008; Hanawa et al. 2010). We find very little change in the accretion rates or average Γ with resolution for the models detailed in Table 2, indicating that our results are numerically converged (see Section 3.3 for an exploration of lower resolution models). Note also that the discs become better resolved as a function of time (as more particles accumulate in our simulation domain) and so the lack of any secular changes also indicates that numerical convergence has been reached.

3.2 Dependence of flow morphology on temperature

To better understand the flow morphology, we averaged the flow velocity on a grid over the final 60 binary orbits (the first 30 orbits were ignored as they represent a transient settling phase). Figs 3 and 4 show streamlines constructed from this averaged flow for the highest resolution cold and hot simulations, respectively. The colour scheme indicates the average density of the flow, normalized to the average surface density of the circumbinary disc.

Figure 3.

Average density and streamlines of the velocity field for the highest resolution cold simulation (simulation 4). The density and velocity field were obtained by averaging the velocity and density on a grid across roughly 60 orbits. The colour scheme shows the average surface density, normalized to the surface density in the circumbinary disc. Contours of the Roche potential and the Lagrange points are shown as in Fig. 1.

Open in new tabDownload slide

Figure 4.

The same as Fig. 3 but for the hot simulation (simulation 8).

Open in new tabDownload slide

In both hot and cold simulations, there are two spiral arms which join the circumbinary disc to the circumstellar discs via the L3/L2 Lagrange points. In both cases, the spiral arm via L2 is much more pronounced than that through L3. This is consistent with the lower value of the Roche potential at L2 (Ochi et al. 2005) which provides an energetically favourable route for material inflowing from the circumbinary disc. The morphology of the flow in the vicinity of L2 is however quite different in the two cases: in the cold simulation the flow is narrow and the trajectory of its peak surface density is offset from L2 by around 60° (measured with respect to the secondary). The flow traverses the Roche lobe (solid line in Fig. 3) more or less radially and impacts the disc edge well inside the Roche lobe. At this point the material becomes entrained in the disc flow.

In the hot simulation, we find that material of similar normalized surface density (coded orange) occupies a broad swathe of space between the secondary and the circumbinary disc and covers nearly half the surface of the secondary's Roche lobe. This much broader structure is simply a result of the much larger pressure scalelength in the hot simulations. Moreover, in contrast to the nearly radial flow seen in the cold simulations, material entering the secondary's Roche lobe describes more or less tangential trajectories, skirting the boundary of the Roche lobe ‘below’ the secondary (at negative y).¹ At the lower right of the secondary (at small negative x), the flow divides into a portion that is retained within the secondary's Roche lobe and one which flows out through the saddle region in the vicinity of L1. The latter flow of material then enters the primary's Roche lobe to the lower left of the primary. The pronounced orange figure of eight structure which wraps around the Roche lobes represents weakly bound material which in some cases orbits the binary several times before eventually being accreted into the disc of one of the stars. As noted above, the main feeding sources for this flow are the two spiral arms from the circumbinary disc, especially that in the region of L2.

To further explore the modes of binary accretion demonstrated in Figs 3 and 4, we computed the mass flux across the boundary of the Roche lobe in our simulations. In 2D the flux is |$F= \Sigma {\boldsymbol v} \bullet {\rm d}{\boldsymbol S}$|⁠, which we compute using SPH interpolation to estimate |${\boldsymbol v}$| (Price 2012, 2007) and calculating the surface normal for the critical Roche equipotential. This flux was then averaged over the final 60 orbits at each point on the boundary of the Roche lobe and plotted in Fig. 5 for the highest resolution hot and cold simulations.

Figure 5.

Flux across Roche lobe boundary for the highest resolution cold (left) and hot (right) simulations. Negative flux implies an inflow. The flux is normalized by the particle injection rate, |$\dot{N}$|⁠. See text for the details of how the flux is calculated.

Open in new tabDownload slide

Fig. 5 (left) shows that in the cold simulation the net flux of particles is into the secondary at all points but is strongly concentrated in the region of the high-density feeding stream (Fig. 3) near L2. For the primary, a significant fraction of the inflow occurs in the vicinity of (0.0, 0.5), which can be seen in Fig. 3 to connect to L3 via a spiral arm. The remainder of the primary's Roche lobe is divided into alternating regions of mild inflow and mild outflow, of which the most significant is the outflow in the region of L1. The outflow near L1 is partly re-accreted by the primary and partly forms a minor input to the adjacent secondary lobe. Overall, the transfer of material between the Roche lobes is small in the cold simulation, with only 6(18) per cent of the material finally accreted by the primary(secondary) having spent some time in the Roche lobe of the other star.

The picture is quite different in the hot simulations. As noted above, the flow into the secondary's Roche lobe in the vicinity of L2 is much broader and the flow of material once it has entered the Roche lobe is nearly parallel to the Roche lobe boundary. Although much of this material is retained by the secondary, a significant minority exits the secondary's Roche lobe on the approach to the L1 point and transitions to the primary's Roche lobe. In the hot simulation, this provides the major accretion stream into the primary's Roche lobe (i.e. at around (−0.6,−0.1)). Since this flow is nearly parallel to the critical equipotential, the remainder of the primary's Roche lobe is again crossed by alternating bands of mild inflow and outflow; the outflow at small positive y in the vicinity of L1 is linked to a corresponding mild inflow to the adjacent portion of the secondary's lobe.

This circulation of material between the lobes does not greatly affect the dominant accretion mode of the secondary. The fraction of material accreting on to the secondary that has ever been in the primary's Roche lobe is only 21 per cent in the hot simulation, a figure scarcely different from that in the cold simulation. However, 65 per cent of the material accreted by the primary in the hot simulation has spent time in the secondary's Roche lobe (compared to only 6 per cent in the cold simulations).

We thus find that the hotter conditions previously modelled by Ochi et al. (2005) and Hanawa et al. (2010) indeed favour an increased net flow from secondary to primary.² This dominant flow boosts the ratio of accretion on to the primary to that on to the secondary when compared with the cold simulation (see Fig. 2).

As noted above, our hot simulation is broadly consistent with model 2-12 of Ochi et al. (2005), a finding which undermines the argument that previous ‘cold’ results reported in the literature were corrupted by numerical issues associated with SPH. Nevertheless, we next discuss whether our results are likely to be corrupted by excessive numerical dissipation.

3.3 Numerical effects

The models presented here represent some of the highest resolution simulations of this problem conducted to date. This statement may seem surprising given the relatively modest numbers of particles involved (see Table 2) but is based on an evaluation of the ratio of a resolution element to the local pressure scaleheight in the vicinity of the Roche lobes (h/H). For our highest resolution cold simulation, the ratio of SPH smoothing length to pressure scaleheight in this region is 0.2–0.5 while the corresponding value for the previous (3D) cold SPH simulations of Bate & Bonnell (1997) is of order unity. For our highest resolution hot simulation this ratio is 0.04–0.1, somewhat less than in the grid-based simulations of Hanawa et al. (2010) and comparable with the simulation of Ochi et al. (2005).³ We note that grid-based codes have not previously been used to model ‘cold’ conditions, because of the difficulty of resolving the inner regions of the circumstellar discs in this case.

Given our agreement with lower resolution results in the literature, it is unsurprising that we have demonstrated that our results are numerically converged. The SPH smoothing length varies by around a factor 5 between models 1/4 and 5/8, yet the difference in fraction accretion rates (shown in Fig. 3) is at the 10 per cent level. This suggests that resolution is not a key factor in determining which component of the binary accretes material. This is difficult to reconcile with the conjecture of Ochi et al. (2005, which ascribed the results from SPH simulations to excessive numerical dissipation in SPH) since such dissipation is a function of resolution.

We find that for all particles, J only undergoes significant changes on time-scales shorter than the dynamical time at shock-fronts. An example is shown in Fig. 6, where a particle in the near radial accretion stream in the cold simulation impacts the secondary's circumstellar disc. Away from obvious shock-fronts (such as where the particle in Fig. 6 strikes the disc), we find J to be conserved to within a few per cent on a dynamical time-scale. This means that the material that is skirting the lower edge of the secondary's Roche lobe in Fig. 4 should change its orbital radius due to numerical viscosity by no more than a few per cent in its first half-revolution around the secondary. Evidently, such a tiny change would have very little effect on which streamlines should remain confined within the secondary's Roche lobe and thus on the resulting flow from secondary to primary.

Figure 6.

The Jacobi constant for a single particle at different points along its trajectory (cold simulation).

Open in new tabDownload slide

We therefore find no evidence that our results are driven by unphysically large numerical viscosity effects. In order to compare grid-based and SPH results in more detail it will be necessary to compare numerically converged results for the same parameters and for simulations that have been run to a quasi-steady state. Such a comparison is currently not possible since the higher resolution simulations reported by Hanawa et al. (2010) are subject to large secular variations in accretion rates throughout the duration of the experiment.

Finally, we explore the effects of degrading the resolution further, motivated by a desire to understand the potential systematic effects of underresolving binary gas flow in large-scale cluster simulations. In Fig. 7, we plot the relative accretion rates on to the primary and secondary for hot and cold simulations, where we have parametrized the resolution in terms of the ratio of the SPH smoothing length to pressure scalelength at the edge of the circumsecondary disc. The lowest resolution simulations contain only a few hundred particles in the circumsecondary disc. We include error bars on each simulation to show the range of variation in Γ in the steady state.

Figure 7.

Average Γ for hot (green triangles) and cold (blue circles) simulations, averaged over the final 20 binary orbits, for a range of resolutions. The error bars show the variation in Γ (1 standard deviation) over the same period. Resolution is reported as h/H (SPH smoothing length to disc scaleheight) at the edge of the circumsecondary disc. The dashed line shows the Γ for purely ballistic particles (Bate 1997).

Open in new tabDownload slide

We find that provided h/H < 1 at the edge of the secondary's disc, the differential accretion pattern is only modestly affected. For the cold simulations, the effect of underresolving the accretion flow is to overestimate accretion on to the primary. The magnitude of this overestimation remains less than 5 per cent so long as h/H < 1, but once h/H > 1 Γ rapidly transitions to the value found for the accretion of purely ballistic particles (Bate 1997). The hot simulations are subject to larger temporal variability in Γ (see Fig. 2 and the error bars in Fig. 7), but the average Γ still only changes by roughly 10 per cent when h/H < 1. There is also some evidence for a step change in Γ similar to the cold simulations when h/H exceeds 1. We caution that the fact that the well-resolved hot simulation have Γ similar to that of ballistic simulations makes detecting changes in Γ at low resolution potentially difficult to detect.

We thus conclude that poor resolution in cluster-scale simulations (where the normalized sound speeds tend to be close to our cold case) is unlikely to severely compromise the mass ratio of binaries formed in such simulations.

4 DISCUSSION

Our isothermal simulations have demonstrated the critical role played by the normalized sound speed in determining the relative accretion rates on to the primary and secondary. In real (non-isothermal) discs, we would expect that it is the gas temperature in the vicinity of the Roche lobes that most critically affects the accretion pattern (since it is here that we have seen the decisive role of pressure gradients in affecting the flow morphology). This allows us to make some crude estimates of the values of normalized sound speed to be expected in astrophysical systems.

In the case of protobinaries, a 1 au binary with primary mass in the regime typical for T Tauri stars will have orbital velocity ∼30 km s⁻¹. For a typical disc, heated mainly by irradiation from the central star (Chiang & Goldreich 1997; D'Alessio, Calvet & Hartmann 2001), we expect T ∼ 100 K for gas on an au scale. These typical numbers yield c ∼ 0.02. Since the temperature of a disc heated by irradiation scales as R^−1/2, the sound speed at the relevant distance (a) scales as a^−1/4. By contrast, the orbital velocity of the binary scales as a^−1/2. Thus, the normalized sound speed scales as a^1/4. Eventually (at ≈100 au), the temperature attains a ‘floor’ value of the order of 10 K and the normalized sound speed rises more steeply with binary separation (∝ a^0.5).⁴

Taken together this suggests that c lies in the range of ∼0.01–0.1 for protobinaries of roughly solar mass. These values would be a factor of a few larger or smaller for earlier or later type pre-main-sequence stars. We thus conclude that c values similar to our cold simulation are characteristic of young binaries although higher values are possible in higher mass, wide systems and/or where heating by accretion significantly boosts the gas temperature.

In the case of supermassive black hole binaries, the aspect ratio of the disc (which is similar to the normalized sound speed in the vicinity of the binary) is a very weak function of system parameters and is of the order of 0.01 (Natarajan & Pringle 1998). Thus, the cold regime is the relevant one in this case also.

5 CONCLUSION

We have investigated the quantitative and qualitative discrepancies previously reported in the literature between SPH and grid-based studies of accretion on to protobinaries. Our results are broadly in line with results in the literature obtained by both methods and demonstrate that the differences reported are in fact due to the different gas temperatures previously used in SPH and grid-based calculations. The fact that (for a given temperature) our (2D SPH) calculations agree with previous 3D SPH calculations and 2D grid-based calculations implies that neither the dimensionality of the simulation nor the numerical technique are critical in determining the modes of accretion on to binaries.

We find that as the ratio of the sound speed in the accreting gas to the orbital speed of the binary is increased (i.e. the temperature is increased), there is an increasingly dominant flow between the secondary and primary Roche lobes. In our cold and hot simulations (where this ratio is, respectively, 0.05 and 0.25), the fraction of particles ultimately accreted by the primary that have previously spent time in the secondary's Roche lobe changes from 6 to 65 per cent. The existence of this cross Roche lobe flow accounts for the slower increase of q in the hotter simulations.

It would be premature to lay out the implications for protobinary evolution from this pilot study at a single value of q. We have however demonstrated that it is possible to obtain numerically converged quasi-steady flow solutions where the role of numerical dissipation in the critical region around the binary Roche lobes is demonstrably small away from shocks. We emphasize that quantitative statements should not be based on the initial transient stage of the simulation. We have found that a quasi-steady state is only attained after ∼30 binary orbits, and we base our accretion rates on long time-scale averaging once this initial phase is concluded.

Finally, we have explored the effect of very poor resolution. We find that provided the ratio of the resolution length to the pressure scaleheight at the edge of the secondary's disc (h/H) remains below 1, the relative accretion rates are only mildly compromised. Our results suggest that even in the poorly resolved regime expected in cluster scale simulations, the evolution of binary mass ratios should be captured correctly.

6 MATERIALS AND METHODS

In the interests of reproducibility and transparency, all the code needed to reproduce this work has been made freely available online at https://bitbucket.org/constantAmateur/binaryaccretion. See the readme file in this repository for further details.

All figures in this paper were generated using the python package matplotilb (Hunter 2007).

Matthew Young gratefully acknowledges the support of a Poynton Cambridge Australia Scholarship. This work has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG.

We are indebted to Eduardo Delgado-Donate who was working on this problem at the time of his untimely death in 2007.

We thank the referee for comments that improved the clarity of the paper.

REFERENCES