The Evolution of Binaries in a Gaseous Medium: Three-Dimensional Simulations of Binary Bondi-Hoyle-Lyttleton Accretion

TL;DR

This study uses 3D hydrodynamic simulations to explore how binary stars interact with surrounding gas, revealing how gas-driven forces influence their orbital evolution and potential coalescence.

Contribution

It provides new insights into the timescales and mechanisms of binary orbital decay in gaseous environments, challenging previous analytic models.

Findings

Binary separation influences bow-shock structure formation.

Center-of-mass motion is slowed faster than orbital inspiral.

Orbital inspiral timescale scales with semi-major axis to the 0.19 power.

Abstract

Binary stars are common. While only those with small separations may exchange gas with one another, even the widest binaries interact with their gaseous surroundings. Drag forces and accretion rates dictate how these systems are transformed by these interactions. We perform three-dimensional hydrodynamic simulations of Bondi-Hoyle-Lyttleton flows, in which a binary moves supersonically relative to a homogeneous medium, using the adaptive mesh refinement code FLASH. We simulate a range of values of the initial semi-major axis of the orbit relative to the gravitational focusing impact parameter of the pair. When the binary separation is less than the gravitational focusing impact parameter, the pair orbits within a shared bow shock. When the pair is wider, each object has an individual bow-shock structure. The long-term evolution of the binary is determined by the timescales for accretion, slowing of the center of mass, and orbital inspiral. We find a clear hierarchy of these timescales; a binary's center-of-mass motion is slowed over a shorter timescale than the pair inspirals or accretes. In contrast to previous analytic predictions, which assume an unperturbed background medium, we find that the timescale for orbital inspiral is proportional to the semi-major axis to the $0.19 \pm 0.01$ power. This positive scaling indicates that gaseous drag forces can drive binaries either to coalescence or to the critical separation at which gravitational radiation dominates their further evolution. We discuss the implications of our results for binaries embedded in the interstellar medium, active galactic nuclei disks, and common envelope phases.

Tables

Table 1. Table 1: a / R BH 𝑎 subscript 𝑅 BH a/R_{\rm BH} in a Variety of Encounters.

Binary Type	Mass $\left[M_{\odot }\right]$	$a$	$a/R_{\mathrm{BH}}$
			Warm ISM		Hot ISM
			Motion in Plane of Disk	Motion Perpendicular to Disk
Stellar Binary	1 + 1	2 $R_{\odot }$	$5\times {10}^{-3}$	0.03 - 0.24	0.05 - 0.24
		1 AU	0.05	2.8 - 25	5 - 26
		0.05 pc	531	${10}^4-{10}^5$	${10}^4-{10}^5$
Massive Stellar Binary	10 + 10	30 $R_{\odot }$	$7\times {10}^{-4}$	0.04 - 0.35	0.07 - 0.4
		1 AU	0.005	0.3 - 2.5	0.5 - 2.6
Binary Black Holes	30 + 30	10 $R_{\odot }$	$8\times {10}^{-5}$	0.004 - 0.04	0.008 - 0.04
		100 $R_{\odot }$	$8\times {10}^{-4}$	0.04 - 0.4	0.08 - 0.4
SMBHs	${10}^6$ - ${10}^8$	1 pc	${10}^{-4}-{10}^{-2}$	0.01 - 10	0.02 - 11
Double Planet	$3.7M_{\mathrm{Jup}}+3.7M_{\mathrm{Jup}}$	3.9 AU	$56$	${10}^3-{10}^4$	${10}^3-{10}^4$

Table 2. Table 2: Orbit-averaged Simulation Data

$a_0$	$⟨{\tau }_{\mathrm{insp}}⟩$	$⟨{\tau }_{\mathrm{stop}}⟩$	$⟨{\tau }_{\mathrm{accr}}⟩$	$⟨\dot{M}⟩$
$\left[R_a\right]$	$\left[{\tau }_{\mathrm{stop},\mathrm{BH}}\right]$	$\left[{\tau }_{\mathrm{stop},\mathrm{BH}}\right]$	$\left[{\tau }_{\mathrm{stop},\mathrm{BH}}\right]$	$\left[{\dot{M}}_{\mathrm{BH}}\right]$
0.0	—	${0.142}_{-0.000}^{+0.000}$	${1.265}_{-0.007}^{+0.004}$	${0.988}_{-0.002}^{+0.007}$
0.16	${0.450}_{-0.016}^{+0.026}$	${0.143}_{-0.003}^{+0.002}$	${1.418}_{-0.070}^{+0.070}$	${0.883}_{-0.006}^{+0.008}$
0.3	${0.503}_{-0.035}^{+0.049}$	${0.145}_{-0.004}^{+0.003}$	${1.595}_{-0.085}^{+0.089}$	${0.786}_{-0.026}^{+0.022}$
0.41	${0.544}_{-0.096}^{+0.088}$	${0.145}_{-0.002}^{+0.002}$	${1.883}_{-0.132}^{+0.145}$	${0.663}_{-0.036}^{+0.043}$
0.55	${0.545}_{-0.088}^{+0.134}$	${0.147}_{-0.004}^{+0.003}$	${2.166}_{-0.181}^{+0.226}$	${0.579}_{-0.049}^{+0.041}$
0.74	${0.628}_{-0.109}^{+0.121}$	${0.148}_{-0.009}^{+0.003}$	${2.132}_{-0.139}^{+0.153}$	${0.588}_{-0.035}^{+0.029}$
1.0	${0.603}_{-0.219}^{+0.521}$	${0.148}_{-0.014}^{+0.006}$	${2.040}_{-0.243}^{+0.151}$	${0.613}_{-0.041}^{+0.079}$

Equations

$$R_{\mathrm{a}}=\frac{2GM}{v_{\infty }^2},$$

$$R_B=\frac{GM}{c_{s,\infty }^2}=\frac{1}{2}{\mathcal{M}}_{\infty }^2{R}_{\mathrm{a}}$$

$${\dot{M}}_{\mathrm{HL}}=\pi R_{\mathrm{a}}^2{\rho }_{\infty }{v}_{\infty }=4\pi G^2{M}^2{\rho }_{\infty }{v}_{\infty }^{-3}$$

$$F_{\mathrm{HL}}={\dot{M}}_{\mathrm{HL}}{v}_{\infty }=4\pi G^2{M}^2{\rho }_{\infty }{v}_{\infty }^{-2},$$

$${\dot{E}}_{\mathrm{HL}}=\frac{1}{2}{\dot{M}}_{\mathrm{HL}}{v}_{\infty }^2=2\pi G^2{M}^2{\rho }_{\infty }{v}_{\infty }^{-1}.$$

$${\tau }_{\mathrm{stop},\mathrm{HL}}\equiv \frac{\frac{1}{2}M{v}_{\infty }^2}{{\dot{E}}_{\mathrm{HL}}}=\frac{v_{\infty }^3}{4\pi G^{2}M{\rho }_{\infty }}.$$

$${\tau }_{\mathrm{accr},\mathrm{HL}}=\frac{M}{{\dot{M}}_{\mathrm{HL}}}=\frac{\frac{1}{2}M{v}_{\infty }^2}{\frac{1}{2}{\dot{M}}_{\mathrm{HL}}{v}_{\infty }^2}=\frac{\frac{1}{2}M{v}_{\infty }^2}{{\dot{E}}_{\mathrm{HL}}}={\tau }_{\mathrm{stop},\mathrm{HL}}.$$

$$\dot{M}_{\rm BH}=\frac{4\pi G^{2}M^{2}\rho_{\infty}}{\big{(}v_{\infty}^{2}+c_{s,\infty}^{2}\big{)}^{3/2}}=\dot{M}_{\rm HL}\bigg{(}\frac{\mathcal{M}_{\infty}^{2}}{1+\mathcal{M}_{\infty}^{2}}\bigg{)}^{3/2}$$

$${\dot{M}}_{\mathrm{BH}}=\pi R_{\mathrm{BH}}^2{\rho }_{\infty }(v_{\infty }^2+c_{s,\infty }^2{)}^{1/2}$$

$$R_{\rm BH}=\frac{2GM}{v_{\infty}^{2}+c_{s,\infty}^{2}}=R_{\rm a}\bigg{(}\frac{\mathcal{M}_{\infty}^{2}}{1+\mathcal{M}_{\infty}^{2}}\bigg{)}.$$

$$F_{\rm BH}=\dot{M}_{\rm BH}(v_{\infty}^{2}+c_{s,\infty}^{2})^{1/2}=F_{\rm HL}\bigg{(}\frac{\mathcal{M}_{\infty}^{2}}{1+\mathcal{M}_{\infty}^{2}}\bigg{)}.$$

$$\dot{E}_{\rm BH}=\text@frac{1}{2}\dot{M}_{\rm BH}(v_{\infty}^{2}+c_{s,\infty}^{2})=\dot{E}_{\rm HL}\bigg{(}\frac{\mathcal{M}_{\infty}^{2}}{1+\mathcal{M}_{\infty}^{2}}\bigg{)}^{1/2}.$$

$$\tau_{\rm stop,BH}\equiv\frac{\frac{1}{2}Mv_{\infty}^{2}}{\dot{E}_{\rm BH}}=\tau_{\rm stop,HL}\bigg{(}\frac{1+\mathcal{M}_{\infty}^{2}}{\mathcal{M}_{\infty}^{2}}\bigg{)}^{1/2}.$$

$$\tau_{\rm accr,BH}=\frac{M}{\dot{M}_{\rm BH}}=\tau_{\rm stop,BH}\bigg{(}\frac{1+\mathcal{M}_{\infty}^{2}}{\mathcal{M}_{\infty}^{2}}\bigg{)}.$$

$$\frac{a}{R_{\mathrm{BH}}}=\frac{GM/v_{\mathrm{orb}}^2}{2GM/(v_{\infty }^2+c_{s,\infty }^2)}=\frac{v_{\infty }^2+c_{s,\infty }^2}{2v_{\mathrm{orb}}^2}=\frac{v_{\infty }^2(1+{\mathcal{M}}_{\infty }^{-2})}{2v_{\mathrm{orb}}^2},$$

$${\mathbf{R}}_{\mathrm{CM}}=\frac{m_1}{M}{\mathbf{r}}_1+\frac{m_2}{M}{\mathbf{r}}_2.$$

$${\mathbf{V}}_{\mathrm{CM}}=\frac{m_1}{M}{\mathbf{v}}_1+\frac{m_2}{M}{\mathbf{v}}_2$$

$${\mathbf{R}}_i$$

$${\mathbf{V}}_i$$

$$a_1=\frac{m_2}{M}a$$

$$v_{\mathrm{orb}}$$

$$P=(\gamma -1)\rho ϵ.$$

$$P_{\infty }=\frac{{\rho }_{\infty }{c}_{s,\infty }^2}{\gamma }=\frac{{\rho }_{\infty }}{\gamma {\mathcal{M}}_{\infty }^2}.$$

$$\dot{m}_{i}=\frac{1}{\Delta t}\int_{\rm sink_{i}}\big{(}\rho-\rho_{\rm sink}\big{)}dV.$$

$$\dot{{\bf p}}_{i}=\frac{1}{\Delta t}\int_{\rm sink_{i}}{\bf v}\big{(}\rho-\rho_{\rm sink}\big{)}dV.$$

$${\mathbf{F}}_{\dot{p}}={\mathbf{F}}_{{\dot{p}}_1}+{\mathbf{F}}_{{\dot{p}}_2}$$

$$d{\bf F}_{\rm DF_{i}}=-\frac{Gm_{i}\rho({\bf r}^{\prime})dV}{|{\bf r}_{i}-{\bf r}^{\prime}|^{3}}\big{(}{\bf r}_{i}-{\bf r}^{\prime}\big{)}$$

$${\bf F}_{\rm DF_{i}}=-\int_{\rm domain}\frac{Gm_{i}\rho({\bf r}^{\prime})dV}{|{\bf r}_{i}-{\bf r}^{\prime}|^{3}}\big{(}{\bf r}_{i}-{\bf r}^{\prime}\big{)}.$$

$${\mathbf{F}}_{\mathrm{DF}}={\mathbf{F}}_{\mathrm{D}{\mathrm{F}}_1}+{\mathbf{F}}_{\mathrm{D}{\mathrm{F}}_2}.$$

$${\dot{M}}_{\mathrm{BH},1}=\frac{4\pi G^2{m}_1^2{\rho }_{\infty }}{(v_{\infty }^2+c_{s,\infty }^2{)}^{3/2}}=\frac{4\pi G^2(M/2{)}^2{\rho }_{\infty }}{(v_{\infty }^2+c_{s,\infty }^2{)}^{3/2}}=\frac{1}{4}{\dot{M}}_{\mathrm{BH}}.$$

Full text

accretion, accretion disks - binaries: general - hydrodynamics - methods: numerical

††software: FLASH (Fryxell et al., 2000); yt (Turk et al., 2011); astropy (Astropy Collaboration et al., 2013).

1 Introduction

Many objects in the universe are found in binary systems. The components of these binaries range from substellar mass objects to supermassive black holes (SMBHs) in the centers of galaxies (Cowley, 1992; Duchêne & Kraus, 2013; Cherepashchuk, 2014; Postnov & Yungelson, 2014). In fact, roughly one-third of all main-sequence (MS) stars are members of binary systems, while $\sim$ 10% of solar-type MS stars and an overwhelming 73% of O-type MS stars have 2 or more stellar companions (Moe & Di Stefano, 2017).

Multiple-star systems are found throughout the galaxy, traveling through gaseous environments with widely varying densities, temperatures, and compositions. The range of gaseous environments include the interstellar medium (ISM; Chanamé & Gould, 2004), accretion flows in active galactic nuclei (AGN; Artymowicz et al., 1993), the dense envelopes of evolved companion stars (Sabach & Soker, 2015; Di Stefano, 2018), and the intracluster medium in young star clusters (Naiman et al., 2011, 2018; Abbate et al., 2018; Roupas & Kazanas, 2019; Kaaz et al., 2019).

A gravitating object moving through gas captures mass and momentum from its surroundings (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944). The drag force it feels due to its gravitationally induced wake robs it of linear momentum (Chandrasekhar, 1943; Ostriker, 1999). For binary or multiple objects, drag forces could lead to a metamorphosis of the size and shape of their shared orbit.

The effects of accretion and drag can have profound consequences. For example, a metal poor star moving through the ISM could accrete, over its history, enough material to significantly enhance its birth metallicity (Shen et al., 2017). More dramatically, the orbit of a stellar-mass binary black hole (BBH) embedded in an AGN disk may be tightened by gaseous interaction to the point that the binary merges under the influence of gravitational wave (GW) radiation (Stone et al., 2017). Whether such binaries form within AGN disks or are dynamically captured (Bartos et al., 2017), they may be the only case of BBH mergers that are accompanied by an electromagnetic signature (due to the surrounding gas; Stone et al., 2017).

Our goal in this paper is to examine the nature of accretion and drag experienced by binaries moving through gaseous environments and to learn the extent to which we can derive physical intuition from the classical Bondi-Hoyle-Lyttleton (BHL) formalism (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952). To that end we perform 3D hydrodynamical studies of binary systems moving supersonically through a gaseous medium.

In the next section, we introduce the BHL formalism for a single accretor, discuss how this formalism can be extended to binary systems, and consider the scales for astrophysically typical encounters. In Section 3, we describe our numerical method and present the results of our hydrodynamic simulations in Sections 4 and 5. In Section 6, we discuss the implications of our results for astrophysical systems. Finally, in Section 7, we conclude.

2 Scales and Typical Encounters in Binary BHL Flows

In this section we define the characteristic scales of BHL and their extension to binary systems. Finally, we briefly discuss typical astrophysical environments for binary BHL flows.

2.1 Characteristic Scales: Single-object BHL

A massive object moving through an initially uniform background gas accelerates matter in its direction as it moves. If the motion of the mass is supersonic relative to the background, a bow shock forms, and shocked gas piles up in a stagnation region downstream of the moving object. If it is absorbing, the object accretes gravitationally bound material from the wake, while the rest is eventually advected away. Although the material in this wake is, therefore, transient, the overall structure of the wake persists, and it exerts a gravitational drag force on the point mass until there ceases to be relative motion between the object and the gas (Chandrasekhar, 1943).

The BHL formalism provides analytical estimates for these interactions. Let $M$ be the mass of the object, $v_{\infty}$ be the speed of the object relative to the gas, and $c_{s,\infty}$ be the sound speed of the gas. The Mach number of the flow is $\mathcal{M}_{\infty}=v_{\infty}/c_{s,\infty}$. The object’s influence on the gas is characterized by its accretion radius. If the object moves highly supersonically, then we can use the Hoyle-Lyttleton accretion radius

[TABLE]

where $G$ is the gravitational constant (Hoyle & Lyttleton, 1939). For very subsonic motion, the accretion radius is given by the Bondi radius:

[TABLE]

(Bondi, 1952).

2.1.1 Hoyle-Lyttleton (HL) Scalings

In general, BHL predicts an accretion rate that is given by the flux of matter through a circular cross section far upstream from the object, where the radius of the cross section is the accretion radius of the object. In the highly supersonic regime, this gives

[TABLE]

(Hoyle & Lyttleton, 1939). The drag force on the object due to the wake is

[TABLE]

and the rate at which energy is dissipated is

[TABLE]

Finally, the approximate timescale for the object to stop is

[TABLE]

In the HL formalism, the accretion timescale is identical to the stopping timescale because

[TABLE]

2.1.2 Bondi-Hoyle (BH) Scalings

For mild Mach numbers, the work of Bondi & Hoyle (1944) and Bondi (1952) led to an interpolation formula for the accretion rate

[TABLE]

where the factor of four was found by Shima et al. (1985). It is helpful to define quantities analogous to those of the HL formalism. If we factor $\dot{M}_{\rm BH}$ similarly to $\dot{M}_{\rm HL}$ in Equation (3), then

[TABLE]

where the characteristic velocity is $(v_{\infty}^{2}+c_{s,\infty}^{2})^{1/2}$ and the BH accretion radius is

[TABLE]

In keeping with Equation (4), we define the BH drag force as an accretion rate times a characteristic velocity, or

[TABLE]

For $\mathcal{M}_{\infty}\ll 1$, $F_{\rm BH}\approx\mathcal{M}_{\infty}^{2}F_{\rm HL}$, whereas the piecewise analytical dynamical friction (DF) formula for a subsonic, nonaccreting perturber goes as $\mathcal{M}_{\infty}^{3}F_{\rm HL}$ (Ostriker, 1999). For a highly supersonic object ($\mathcal{M}_{\infty}\gg 1$), $F_{\rm BH}\to F_{\rm HL}$.

Similarly, we can define the characteristic energy dissipation rate to be

[TABLE]

The translational kinetic energy of the object is $\text@frac{1}{2}Mv_{\infty}^{2}$, which gives a BH stopping time of

[TABLE]

The BH accretion timescale is

[TABLE]

While the above equations are only approximately correct, coefficients for drag and accretion for the above canonical case have been obtained numerically by Shima et al. (1985), Ruffert & Arnett (1994), Ruffert (1994, 1995), and Blondin & Raymer (2012). See Edgar (2005) for a review of the BHL accretion formalism, including derivations of the above equations, and Ostriker (1999) for the first analytical treatment of the gaseous DF drag approximated by the expressions for $F_{\rm HL}$ and $F_{\rm BH}$.

2.2 Characteristic Scales: Binary BHL

Now consider what happens if we replace the single object of mass $M$ with two objects, $m_{1}$ and $m_{2}$, whose total mass is $M$, and whose center of mass (CM) moves with velocity ${\bf v}_{\infty}$ relative to the gas. Free parameters now include the binary mass ratio, the initial eccentricity of the orbit, the inclination of the orbital plane relative to ${\bf v}_{\infty}$, and the semi-major axis of the orbit, $a$. We will examine the simplifying case of an equal-mass binary in an initially circular orbit with inclination of $\pi/2$ (that is, the incoming gas intercepts the binary “edge-on” to the orbit).

Now there are two key length scales in the problem: the gravitational focusing scale, $R_{\rm BH}$, and the binary semi-major axis, $a$. Here, $R_{\rm BH}$, defined by Equation (10), is the BH accretion radius of the CM. We can re-express this ratio of length scales in terms of other parameters of the binary:

[TABLE]

where we have used the orbital speed, $v_{\rm orb}=\sqrt{GM/a}$, for a binary in a circular orbit. Therefore, in the limit of highly supersonic motion (large $\mathcal{M}_{\infty}$), $a/R_{\rm BH}\approx v_{\infty}^{2}/2v_{\rm orb}^{2}$.

We can imagine two limiting cases for the ratio of $a/R_{\rm BH}$. When $a\ll R_{\rm BH}$, the binary appears essentially single as gas is focused toward the system’s CM from the background medium. In this case, we might imagine that the derivations of the previous section provide a reasonable description of flow around, accretion toward, and drag forces on the system’s CM. In the opposite limit, when $a\gg R_{\rm BH}$, the binary is much wider than the gravitational focusing scale. In this case, we might expect that each object undergoes an essentially isolated BHL interaction with the gaseous background. That is, the system evolves according to the BHL prescription for two independent objects (each of mass $\frac{1}{2}M$). In this limit, we would expect the accretion rate of the system and the total DF force on the binary to be $\propto(m_{1}^{2}+m_{2}^{2})=M^{2}/2$.

If we allow $a\sim R_{\rm BH}$, the two length scales are on the same order. Gas is focused toward both the individual objects and the system’s CM. Orbital motion projects phase-dependent velocities for the individual objects with magnitude similar to that of the CM motion. Without a clear expectation for the flow structures in this intermediate regime, we turn to hydrodynamical simulations to provide a solution for the small- and large-scale structures that contribute to the evolution of the orbital and CM motion of the binary. The purpose of the remainder of this paper is to explore binary BHL flow with varying $a/R_{\rm BH}$ and to characterize the evolution of a binary as a function of this ratio.

2.3 Typical Encounters in Binary BHL

Astrophysical binary systems span many orders of magnitude in mass and semimajor axis and these systems interact with a wide variety of ambient gases. In Table 1, we calculate $a/R_{\rm BH}$ for a sampling of binary systems moving through the ISM. For each gaseous medium in Table 1, we assume an isothermal ideal gas and calculate the sound speed as $c_{s,\infty}=\sqrt{k_{B}T/m_{p}}$. We see from the table that a stellar-mass binary with a separation of 1 AU moving through the warm ISM may have $a/R_{\rm BH}$ values from several hundredths to a few times 10. Increasing the mass of the binary by a factor of 10 gives results of $0.3\lesssim a/R_{\rm BH}\lesssim 3$ when crossing the disk.

We discuss ambient gases of much greater density in Section 6. For example, a binary embedded in the model AGN disk presented in Section 6.4 may encounter number densities, $n_{\infty}$, of $10^{8}-10^{14}$cm*-3*, depending on its radial location in the disk. A stellar-mass binary with a separation of 1 AU located at a radius of 2 pc in this disk will have $a/R_{\rm BH}\approx 0.6$. The same binary embedded at a radius of 0.5 pc will have $a/R_{\rm BH}\approx 2.5$. As another example, we discuss in Section 6.5 a common envelope (CE) interaction in the triple system PSR J0337+1715 . We find that the models of Sabach & Soker (2015) imply envelope densities of $10^{16}\lesssim n_{\infty}\lesssim 10^{19}$ and the embedded binary would have $0.8<a/R_{\rm BH}<2.7$. The scaling with $\rho_{\infty}$ in Equation (13) indicates that the binary in one of these environments will transform 10 orders of magnitude more quickly than in the ISM.

3 Numerical Approach

We perform simulations of supersonic BHL flows around an embedded binary system. The traditional approach to numerical studies of BHL is to perform simulations in the frame of a point mass by constructing a wind tunnel that sends gas flowing past a single gravitating sink particle that is fixed at the coordinate origin. The sink particle accretes gas from the domain, and the forces experienced by the sink particle due to the gas can be recorded at each time step. We modify this approach by instead placing two sink particles in orbit within a wind tunnel. The sink particles are free to move about the domain, and we start the simulations in the initial frame of the CM of the binary. The sink particles accrete gas from the domain and actively respond to the forces they experience due to the gas. In the following subsections, we give a brief history of numerical studies of BHL that are most relevant to the current work, we describe our numerical method, and we outline the parameters chosen for the suite of simulations that we run.

3.1 Previous Numerical Approaches to BHL

BHL flows have a long history of hydrodynamical study, beginning with Hunt (1971), which has been reviewed by Edgar (2004) and Foglizzo et al. (2005). Of particular importance to our work are the studies of Ruffert & Arnett (1994) and Ruffert (1995, 1996, 1999), who performed Eulerian hydrodynamic simulations of BHL. They model the central, accreting boundary condition with a low-pressure sink that effectively absorbs all incoming material and removes it from the domain. These papers examine the effects of gas adiabatic index, Mach number, sink size, resolution, and the size of the computational domain on the rates of accretion and drag. They showed the qualitative features of BHL flows in 3D: a symmetric bow shock sets up in front of the accretor with an opening angle and standoff distance that are set by $\gamma$ and $\mathcal{M}_{\infty}$. High-density material collects in a downstream wake from which material is accreted or is advected downstream. Quantitatively, they found higher accretion rates with smaller adiabatic index and with lower Mach number.

Ruffert’s results, along with the more recent, higher resolution simulations of Blondin & Pope (2009) and Blondin & Raymer (2012), have shown that accretion rates are also sensitive to the size of the computational domain, the size of the sink, and the resolution of the grid. We therefore studied those three aspects of our numerical setup before proceeding with our production runs. The results of those studies are presented in Appendix B.

In this paper, we use a very similar numerical approach to extend the efforts of Ruffert & Arnett (1994), Ruffert (1995, 1996, 1999), Blondin & Raymer (2012), MacLeod & Ramirez-Ruiz (2015), MacLeod et al. (2017), and Murguia-Berthier et al. (2017) to binary systems by simulating two sink particles orbiting within a supersonic BHL flow. While methodologically distinct, we note that a similar scenario has been simulated in 3D by Farris et al. (2010) in the ultrarelativistic regime in which a pair of black holes decays rapidly to merger while traversing through a surrounding gaseous medium.

3.2 Notation for Two-body Orbits

Before we describe our numerical method in further detail, it is necessary to define the notation we will use to describe the orbit of the binary. We follow the notation of Hilditch (2001, Chapter 2) with the exception that we use $M$ to represent the total mass of the binary. Let us summarize the quantities relevant to this paper. If ${\bf r}_{1}$ and ${\bf r}_{2}$ are the positions of $m_{1}$ and $m_{2}$, respectively, relative to the coordinate origin, then the CM is located at

[TABLE]

Similarly, the velocity of the CM is

[TABLE]

where ${\bf v}_{1}$ and ${\bf v}_{2}$ are the velocities of $m_{1}$ and $m_{2}$ relative to the origin. Consistent with Hilditch (2001), we use capital letters with subscripts for the barycentric quantities. The barycentric position and velocity of particle $m_{i}$ are, respectively,

[TABLE]

and

[TABLE]

The position of $m_{1}$ relative to $m_{2}$ is ${\bf r}={\bf r}_{1}-{\bf r}_{2}$ and the semi-major axis of this “relative” orbit is denoted by $a$. The semi-major axis of the barycentric orbit of $m_{1}$ is

[TABLE]

and similarly for $m_{2}$. For a circular orbit, the orbital speed is

[TABLE]

3.3 Hydrodynamical Setup

We model BHL accretion onto a binary system moving relative to a uniform gas using FLASH (Fryxell et al., 2000). FLASH is an Eulerian hydrodynamic code with adaptive mesh refinement (AMR). We use FLASH’s directionally split Piecewise Parabolic Method Riemann solver (Colella & Woodward, 1984). The gas pressure is defined by an ideal gas, gamma-law equation of state,

[TABLE]

where $\rho$, $\epsilon$, and $\gamma$ are the density, internal energy, and adiabatic index of the gas, respectively. We initialize a 3D Cartesian grid with the coordinate origin at the center of the grid. The size of the computational domain in $x$, $y$, and $z$ is $L$, giving a domain volume of $L^{3}$.

We use the active sink particles (ASP) unit in FLASH (Federrath et al., 2010) to model the binary system. We represent the binary with two equal point masses, $m_{1}$ and $m_{2}$, placed at positions ${\bf r}_{1}$ and ${\bf r}_{2}$ relative to the grid. We begin the simulations with the CM of the binary at rest (${\bf V}_{\rm CM}={\bf 0}$) at the origin of the coordinate system (${\bf R}_{\rm CM}={\bf 0}$). The two masses are initialized in a circular orbit of semi-major axis $a_{0}$ in the $x-y$ plane with the orbital angular momentum vector pointing in the $+z$ direction. For comparison to single-particle BHL, we also include the option to place a single particle of mass $M$ on the grid instead of a binary.

To model the supersonic motion of the binary relative to a uniform gas, we fill the domain with a wind of density $\rho_{\infty}$ that flows in the $+x$ direction at speed $v_{\infty}$. The wind speed is related to the sound speed, $c_{s,\infty}$, through the Mach number $\mathcal{M}_{\infty}=v_{\infty}/c_{s,\infty}$. The initial gas pressure is set according to

[TABLE]

An inflow boundary condition at the $-x$ boundary of the computational domain replenishes this homogeneous wind throughout the simulation. All other boundaries of our computational domain have outflow conditions, which extrapolate the thermodynamic quantities and allow gas to freely leave and enter the domain. In our code units, $v_{\infty}=1$. With this choice and $M=(2G)^{-1}$, $R_{a}=2GM/v_{\infty}^{2}=1$. The result is a dimensionless setup with length unit $R_{a}$ and time unit $\,R_{\rm a}/v_{\infty}$.

3.4 ASP Dynamics and Diagnostics

We utilize many features built in to the ASP unit of FLASH with some additional diagnostics. See Federrath et al. (2010) for a complete description; we briefly describe our implementation here. The two massive particles that represent our binary system are gravitationally softened point masses surrounded by a spherical inflow boundary of radius $R_{\rm s}$. The softening radius of each point mass is $0.4R_{\rm s}$. Once initialized, each particle moves in response to gravitational accelerations from the other particle and the gas in the domain. The inflow boundary, or sink, that surrounds each particle removes gas from the domain. Once material falls into the sink, it is added to the sink particle properties and then it is deleted. This boundary condition represents perfectly efficient accretion without feedback. In our implementation of the ASP, we turn off the creation of new sink particles and we turn off self-gravity of the gas.

3.4.1 Accreted Quantities

Gas absorbed by sink particles carries mass and momentum, which are acquired by the particle. The accretion step is performed by integrating the quantity over all cells within the volume of the sink and by adding the summed values to the particle’s properties. To compute an accretion rate, we divide each accreted quantity by the time step. The accretion rates are recorded, and then the accreted gas is deleted by setting the density within the sink cells to $\rho_{\rm sink}=10^{-2}\rho_{\infty}$, and setting all components of the gas velocity to zero. The mass accretion rate for particle $m_{i}$ is calculated as a volume integral over the sink

[TABLE]

The total mass accretion rate for the binary is $\dot{M}\equiv\dot{m}_{1}+\dot{m}_{2}$.

The accretion of linear momentum is integrated in each Cartesian coordinate. For particle $i$, this is given by

[TABLE]

The accretion of linear momentum represents a force on each particle. We will denote these “momentum transport forces” on $m_{1}$ and $m_{2}$, respectively, as ${\bf F}_{\dot{p}_{1}}$ and ${\bf F}_{\dot{p}_{2}}$. Finally,

[TABLE]

is the total momentum transport force on the binary CM.

3.4.2 Dynamical Friction

Each sink particle exerts a gravitational force on the gas, which restructures the gas in the domain. In turn, the redistributed gas exerts a gravitational force on each particle. This gravitational drag force, or gaseous DF, will occur whenever the gas distribution is not spherically symmetric about the particle. The total (vector) DF force acting on each particle is found by summing up the gravitational force on the particle that is due to each cell in the domain. Recalling that ${\bf r}_{i}$ is the location of $m_{i}$ in the grid and letting ${\bf r}^{\prime}$ denote the location of a particular gas cell, we find the force on $m_{i}$ due to the cell is

[TABLE]

and the total DF force on $m_{i}$ is the integral over the domain

[TABLE]

The total DF force on the CM is

[TABLE]

3.4.3 Particle Advancement

The CM position and CM velocity of each particle are updated to account for the accreted gas. Before the gas is evolved further, each particle’s motion is advanced using ASP’s leapfrog integrator. The advancement step applies the momentum transport, DF, and sink-on-sink accelerations to each sink particle. At each time step, we record the accretion rates and forces as well as the total mass, position, and velocity of each particle.

3.5 Simulation Parameters

To summarize, the free physical parameters of our numerical setup are the background density $\rho_{\infty}$, the flow Mach number $\mathcal{M}_{\infty}$, the adiabatic index $\gamma$, and the initial semi-major axis of the binary $a_{0}$. In order to study the effects of changing binary separation relative to the gravitational focusing scale, we vary $a_{0}$ while holding $\mathcal{M}_{\infty}$ and $\gamma$ constant. We simulate seven values of $a_{0}/R_{a}$: 0.0, 0.16, 0.3, 0.41, 0.55, 0.74, and 1.0. Each of these simulations adopts $\rho_{\infty}=10$, $\gamma=5/3$, and $\mathcal{M}_{\infty}=2$.

We have studied the effects of domain size, spatial resolution, sink radius, and background density on our results. We present these validations of our method in Appendix B. As a result of these studies, we adopt the following in our production runs. We use a domain of volume $L^{3}=(32R_{a})^{3}$, with the CM of the binary initialized at the center. The base-level mesh consists of four blocks of $8^{3}$ cells in the $x$, $y$, and $z$ directions. We adaptively apply between one and seven levels of refinement. The finest grid cells are $\delta_{\rm min}=0.016R_{a}$ wide in $x$, $y$, and $z$. The ASP unit maximally refines the region of each sink particle. We further instruct the AMR unit to refine on the second derivative of the internal energy of the gas. The radius of each sink particle is set to $R_{\rm s}=0.05R_{\rm a}$ so that each sink has $R_{\rm s}/\delta_{\rm min}=3.125$ grid cells across its radius. We run the simulations to $t=50\,R_{\rm a}/v_{\infty}$, which is 18$\,R_{\rm a}/v_{\infty}$ beyond the box-crossing time of the wind. We note that the $a_{0}=1.0R_{\rm a}$ simulation falls in the $a\gtrsim R_{\rm BH}$ regime because $1.0R_{a}=1.25R_{\rm BH}$.

4 Numerical Results

In this section, we describe the results of our gas dynamical simulations of binary BHL flows. We analyze the morphology of these flows, the resultant accretion onto the individual members, and the forces on the binary CM and orbit.

4.1 Flow

The presence of a binary significantly alters the morphology of BHL accretion flows. At a given time, the manifestation of the binary’s presence depends strongly on $a_{0}$ as compared to $R_{\rm a}$, and on orbital phase. We explore these two aspects of BHL flows in the following subsections and in Figures 1 - 4. These figures show slices through the orbital ($z=0$) plane.

4.1.1 Binary Separation and Binary BHL Flows

Figure 1 illustrates the effect of semi-major axis on flow morphology. Here we plot density (in units of $\rho_{\infty}$) for our single-particle simulation (upper left) and for binary simulations with five different values of semi-major axis. The single particle, with the same total mass as the binaries, represents the limit of zero binary separation ($a_{0}\rightarrow 0.0R_{\rm a}$).

The $a_{0}=0.0R_{\rm a}$ snapshot of Figure 1 illustrates key features of canonical BHL flow. Supersonic motion relative to the background gas sets up an upstream bow shock front that funnels material into a downstream wake. The shock front, which has an opening angle proportional to $c_{s,\infty}/v_{\infty}$, is smooth and axisymmetric about the line of motion of the point mass.

With the introduction of a binary, rather than a single, central object, some aspects of the flow remain similar while others diverge substantially. With $a_{0}=0.16R_{\rm a}$, the size and shape of the primary shock remain very similar to that of $a_{0}=0.0R_{\rm a}$. As $a_{0}$ increases, the effect of the binary on the overall bow shock structure becomes more apparent. By $a_{0}=1.0R_{\rm a}$, each particle sets up its own primary shock (during parts of the orbit), and the overall structure is the superposition of the two shock cones.

Figure 1 also shows that the postshock density and velocity structures are also modified by the presence of a binary. In the single-object case, the highest densities surround the accreting object. Flow converges in the wake and compresses as it falls inward toward the accretor. With increasing binary separation relative to $R_{\rm a}$, we see that this density enhancement divides into separate convergences near $m_{1}$ and $m_{2}$. This is particularly apparent in the $a_{0}=1.0R_{\rm a}$ snapshot. We can understand this transition by considering that, with the division of the central object into a binary, the formally divergent potential at the origin (with $a_{0}=0.0R_{\rm a}$) has been replaced by a potential that is softened on the scale of the binary separation. This restructuring of material in the postshock region has implications for the rates of accretion and drag, which we discuss in detail in Sections 4.2 and 4.3.

Trailing the bow shock, flow is relatively smooth in the case of a single accretor. With the introduction of a binary, orbital motion imparts new features on the postshock flow. For compact binary separations (e.g. $a_{0}=0.16R_{\rm a}$), some of the primary features are spiral shocks launched because the orbital motion of the binary is supersonic relative to the surrounding gas. These features are most apparent in Figure 2, in which we plot the velocity divergence, $\nabla\cdot{\bf V}_{\rm gas}$. Shocked material is revealed in purple in the figure, where $\nabla\cdot{\bf V}_{\rm gas}<0$ indicates abrupt convergence of gas at a shock front. Especially for the binaries with $a_{0}=0.16R_{\rm a}$ and $a_{0}=0.3R_{\rm a}$, spiral shocks propagate outward away from the binary and are advected downstream in the wake.

As the binary separation increases, the relative velocity between the binary object’s orbital motion and the gas decreases. The particle moving through $-y$ positions transitions to subsonic relative motion first, because the particle in the lower half of the orbit moves in the same direction as the gas. This transition to subsonic relative motion is visible in the bottom row of Figure 2 through the disappearance of the leading shock front and in Figure 1 through the symmetric density distribution around the object. In the $a_{0}=0.41R_{\rm a}$ slice of Figure 2, the secondary shock ahead of the lower particle is barely discernible. By $a_{0}=1.0R_{\rm a}$, both secondary shocks are absent, and the only inflection of the binary envelope is from the superposition of two primary shock fronts. Here, all orbital motion of the binary is subsonic relative to the postshock sound speed.

4.1.2 Time-dependent Properties of the Flow

In the simulations involving a binary central object, orbital motion imparts time variability and orbital phase dependence to the flow structure. We examine those features in detail here. Figures 3 and 4 show the time series of the flow for $a_{0}=1.0R_{a}$ and $a_{0}=0.41R_{a}$, respectively. In both of these time series, we show the binary over half of an orbit only. The second half of each orbit is qualitatively similar to the first, except that the roles of $m_{1}$ and $m_{2}$ are reversed.

Figure 3 shows $a_{0}=1.0R_{a}$ over the half-orbit that begins at $t=42.1\,R_{\rm a}/v_{\infty}$. In the first image, $m_{1}$ is fully upstream from $m_{2}$. Although at this instant the particles lie along a line parallel to the wind direction (at $\infty$), the bow shock is not symmetric about this line but carries distortions that highlight the preceding orbital motion of the binary. We will trace the origin of this distortion through the remaining snapshots of Figure 3.

Situated downstream, $m_{2}$ interacts with material that has already been influenced by $m_{1}$. This dense, lower-velocity material piles up around $m_{2}$, resulting in a higher-than-average, and almost spherical, density distribution about $m_{2}$. This density enhancement remains with $m_{2}$ as time advances across the upper row of Figure 3, approaching the primary bow shock by the first image of the second row.

At $t=43.7\,R_{\rm a}/v_{\infty}$, a second shock front is forming ahead of $m_{2}$ as the particle approaches its greatest velocity relative to the gas. Between this snapshot and the next, $m_{2}$ crosses the leading boundary of the bow shock of $m_{1}$. The overdensity that has been carried along by $m_{2}$ remains with the particle during the crossing. In the next two images ($t=44.3$ and $44.8$ $\,R_{\rm a}/v_{\infty}$), this dense material is impacted by the incoming wind. By $t=45.4\,R_{\rm a}/v_{\infty}$ (first panel of third row), the over-density is largely detached from $m_{2}$ and is being pulled into the approaching potential well of $m_{1}$. In the same panel, a bridge of material connects the two particles and, for the first time in these snapshots, the density enhancement about $m_{1}$ is greater than that of $m_{2}$.

Two flow properties transition in the last row of Figure 3. First, the dense material that is attached to $m_{2}$ at $t=45.4\,R_{\rm a}/v_{\infty}$ becomes fully detached from the particle as the downstream region is increasingly dominated by the gravitational pull of $m_{1}$. Second, as $m_{2}$ moves down and to the left, the shock set up by the particle rotates with it, leading the $+y$ portion of the bow shock to protrude farther upstream than the $-y$ portion.

Figure 4 shows the $a_{0}=0.41R_{a}$ binary for the half-orbit beginning at $t=40.1\,R_{\rm a}/v_{\infty}$. We additionally annotate several features, which will be referenced here. As discussed in the context of Figure 2, the particle traveling upwind through the post-bow-shock material launches a secondary spiral shock wave behind the primary shock.

The interaction of the binary with spiral shocks sets up the fan-like structures emanating from the binary in Figure 4. Material crossing the primary shock interacts with each particle twice, resulting in four fan structures per orbit. In the first snapshot ($t=40.1\,R_{\rm a}/v_{\infty}$), the particles are located at $x\approx 0$. There is a large standoff distance between the primary shock and the binary. Behind the primary shock, a secondary shock is set up in front of $m_{1}$. We label these structures “1” and “2,” respectively, in all of the snapshots shown, so that one can follow each structure over time. From $t=40.2$ to $40.3\,R_{\rm a}/v_{\infty}$, shock 1 is deflected by the gravitational pull of $m_{1}$. By $t=40.4\,R_{\rm a}/v_{\infty}$, shock 2 is pushed past shock 1 by $m_{1}$ as the particle advances in the $-x$ direction. From $t=40.6$ to $t=41.0\,R_{\rm a}/v_{\infty}$, $m_{2}$ advances toward $+y$ values, and shock 1 is defected by the gravitational pull of $m_{2}$. At $t=41.2\,R_{\rm a}/v_{\infty}$, a new secondary shock (set up by $m_{2}$ and labeled “3” in Figure 4) crosses shock 1 as shock 1 continues to advect downstream. In the final image, shock 3 is advancing to meet shock 2 and will cross shock 2 as the orbit continues.

In addition to the fan structures, the density very close to the binary is of interest. In particular, it appears that the particle at $+y$ positions exhibits a higher local density than its counterpart at $-y$ values. A persistent bridge of dense material always connects the two particles, giving a slightly higher density enclosed by the orbit than in the $a_{0}=1.0R_{\rm a}$ simulation.

4.2 Accretion

In this section, we examine instantaneous and orbit-integrated accretion rates onto the objects.

4.2.1 Instantaneous Accretion Rates

In Figure 5, we plot the instantaneous accretion rates versus time for $a_{0}=0.16,0.41$, and $1.0R_{a}$ for simulation times $\geq 30\,R_{\rm a}/v_{\infty}$. As expected from prior studies of canonical BHL accretion (with $R_{\rm s}=0.05R_{a}$, $\mathcal{M}_{\infty}\sim 2$, and at our resolution), $\dot{M}(a_{0}=0)\approx 0.7\dot{M}_{\rm HL}\approx 1.0\dot{M}_{\rm BH}$ (Ruffert, 1994; Blondin & Raymer, 2012).

Although $\dot{M}(a_{0}=0)$ reaches a steady state, the accretion rates for the members of the binary systems are periodic. For the tightest binary, $a_{0}=0.16R_{a}$, and $\dot{m}_{1}$ and $\dot{m}_{2}$ are exactly out of phase, each with an accretion rate that follows the orbital period. As such, $\dot{M}$ is almost constant in time and approaches $\dot{M}(a_{0}=0)$. For more widely separated binaries, the individual accretion rates of the two particles no longer sum to a constant value. In the case of $a_{0}=1.0R_{a}$, $\dot{m}_{1}$ shows a sharp spike once per orbit before declining to a minimum value, and then the rate climbs to a local maximum before the next sharp spike in accretion. While the curve for $\dot{m}_{2}$ looks very similar, the variation in $\dot{M}$ reveals that the individual accretion rates are not exactly out of phase. In the case of $a_{0}=0.41R_{a}$, the individual accretion rates reveal one peak and one valley per orbit with transitions between the two that are variable from one orbit to the next. In addition, $\dot{m}_{1}$, for example, reveals a slight overall increase and decrease over many orbits. The superposition of these trends in $\dot{m}_{1}$ and $\dot{m}_{2}$ yields an $\dot{M}$ that has irregularly spaced peaks and valleys with variation from orbit to orbit.

Let us now connect features of the accretion rates depicted in Figure 5 to the flow morphology discussed in Section 4.1.2. In Figure 6, we show $\dot{m}_{1}$, $\dot{m}_{2}$, and $\dot{M}$ for the $a_{0}=1.0R_{a}$ simulation for simulation times $t=39.9-48.7$ $\,R_{\rm a}/v_{\infty}$, as in the panels of Figure 3. In the first row of Figure 3, $m_{2}$ moves from a position directly downstream from $m_{1}$ at $t=42.1\,R_{\rm a}/v_{\infty}$ into a lower-density pocket behind the primary shock by $t=43.2\,R_{\rm a}/v_{\infty}$. Figure 6 reveals a decreasing accretion rate for $m_{2}$ during the same period of time. Although $m_{2}$ has carried along higher-density material during this time interval, the material is able to expand into the low-density region behind the primary shock, giving rise to the local minimum in $\dot{m}_{2}$ at $t\approx 43.2\,R_{\rm a}/v_{\infty}$.

Between $t=43.2$ and $43.7\,R_{\rm a}/v_{\infty}$, $m_{2}$ has turned into the oncoming wind and approaches the primary shock front; the particle nears its highest velocity relative to the gas at $t=43.7\,R_{\rm a}/v_{\infty}$ and a second shock front forms in front of $m_{2}$. Meanwhile, the high-density concentration that has followed $m_{2}$ is now stopped by this second shock, and a stream of material flows from $m_{1}$ to $m_{2}$. This confluence of flow gives $\dot{m}_{2}$ a sharp peak at $t=43.7\,R_{\rm a}/v_{\infty}$.

By $t=44.3\,R_{\rm a}/v_{\infty}$, $m_{2}$ has crossed the primary bow shock of $m_{1}$, and now $\dot{m}_{2}$ declines. From $t=44.3$ to $44.8\,R_{\rm a}/v_{\infty}$, $m_{2}$ begins to lose the high-density clump that has followed the particle since its time downstream from $m_{1}$. Then $\dot{m}_{2}$ bottoms out as the clump of material finally detaches from $m_{2}$ and is swept downstream ($t=45.4$ to $46.5\,R_{\rm a}/v_{\infty}$). At the same time, $m_{1}$ is moving into a position downstream from $m_{2}$. There is an increase in $\dot{m}_{1}$ from $t=44.3-45.9\,R_{\rm a}/v_{\infty}$ as $m_{1}$ moves with its slowest relative velocity compared to the surrounding gas and interacts with dense material stripped from the vicinity of $m_{2}$. From $45.9$ to $46.5\,R_{\rm a}/v_{\infty}$, $\dot{m}_{1}$ declines slightly as the clump of material that has detached from $m_{2}$ crosses outside of the orbit of $m_{1}$ as it is advected downstream and out of reach.

Figure 7 shows the accretion rates for the $a_{0}=0.41R_{\rm a}$ simulation for the time period $t=39.6$ to $t=41.9\,R_{\rm a}/v_{\infty}$, also shown in the snapshots of Figure 4. From $t=40.1$ to $40.8\,R_{\rm a}/v_{\infty}$, $\dot{m}_{2}$ steadily increases as $m_{2}$ moves in the same direction as the flow. At $40.8\,R_{\rm a}/v_{\infty}$, $m_{2}$ turns upstream. There is a steepening of $\dot{m}_{2}$ to a maximum at $41.0\,R_{\rm a}/v_{\infty}$ as $m_{2}$ deflects shock 1 while also stripping material from the region surrounding $m_{1}$. Between $41.0$ and $41.2$ $\,R_{\rm a}/v_{\infty}$, shock 3 forms in front of $m_{2}$. Once this material crosses shock 1 (at $41.2\,R_{\rm a}/v_{\infty}$), it expands into the lower-density region behind shock 2, reducing the rate of accretion onto $m_{2}$.

Meanwhile, $m_{1}$ experiences a decreasing accretion rate from $t=40.1\,R_{\rm a}/v_{\infty}$ to $t=40.8\,R_{\rm a}/v_{\infty}$ as the particle takes up a position in front of $m_{2}$. There is a flattening of $\dot{m}_{1}$ from $t=40.2$ to $40.3\,R_{\rm a}/v_{\infty}$ when shock 1 is pulled into $m_{1}$ and shock 2 crosses shock 1. The overall decrease in $\dot{m}_{1}$ continues through $t=40.8\,R_{\rm a}/v_{\infty}$ while $m_{1}$ is in the upstream position and $m_{2}$ pulls material away from $m_{1}$. Just before $t=40.8\,R_{\rm a}/v_{\infty}$, the slope of $\dot{m}_{1}$ turns positive as $m_{1}$ moves to $-y$ values. There, $m_{1}$ moves in the same direction as the wind, reducing the relative motion between the particle and the gas. From there the cycle depicted in Figure 4 begins again, but this time with the roles of $m_{1}$ and $m_{2}$ reversed.

4.2.2 Orbit-averaged Accretion Rates

Having examined the origins of time-variable accretion rates in binary BHL flows, we now compare time-integrated quantities for binaries of varying separation.

It is apparent from Figure 5 that $\dot{m}_{1}$, $\dot{m}_{2}$, and $\dot{M}$ for $a_{0}=0.41$ and $1.0R_{\rm a}$, are lower than the corresponding accretion rates in the $a_{0}=0.16R_{a}$ case. To look for trends in the accretion rate as a function of $a_{0}$, let us compare the time-averaged total accretion rate for each pair. We use angled brackets, $\langle~{}\rangle$, to indicate quantities averaged over an integer number of orbits. Each time average begins at $t=30\,R_{\rm a}/v_{\infty}$ and terminates at the end of the last complete orbit before the end of the simulation ($t=50\,R_{\rm a}/v_{\infty}$). The $a_{0}=1.0R_{a}$ binary (with $P_{\rm orb}\approx 8.89\,R_{\rm a}/v_{\infty}$) completes two full orbits during this time. The $a_{0}=0.16R_{a}$ binary completes 33 full orbits.

In Figure 8, we plot the orbit-averaged total accretion rate versus $a_{0}$. As the binary separation decreases, $\langle{\dot{M}}\rangle$ approaches $\dot{M}_{\rm BH}\propto M^{2}\propto(m_{1}+m_{2})^{2}$. That is, binaries of smaller separations accrete like a single particle of mass $m_{1}+m_{2}$. Wider binaries, on the other hand, have a lower accretion rate and fall closer to $\langle{\dot{M}}\rangle=0.5\dot{M}_{\rm BH}$, which assumes that $m_{1}$ and $m_{2}$ accrete independently of one another. To see where the factor of $0.5$ comes from, let us write the BH accretion rate for $m_{1}=M/2$:

[TABLE]

Because $m_{2}=M/2=m_{1}$, $\dot{M}_{\rm BH,2}=\dot{M}_{\rm BH,1}$. Then the total accretion rate for the pair (treating them as independent particles) is $\dot{M}_{\rm BH,independent}=0.25\dot{M}_{\rm BH}+0.25\dot{M}_{\rm BH}=0.5\dot{M}_{\rm BH}$. For more general $m_{1}$ and $m_{2}$, we would say that because $\dot{M}_{\rm BH,1}\propto m_{1}^{2}$ and $\dot{M}_{\rm BH,2}\propto m_{2}^{2}$, $\dot{M}_{\rm BH,independent}\propto(m_{1}^{2}+m^{2}_{2})$.

4.3 Forces

When either a single object or binary moves through a gaseous medium, it restructures the distribution of the surrounding gas. The restructured gas exerts forces on the CM and the components of the binary. In this section, we describe how these forces are calculated from our numerical simulations, and we show an example calculation for our $a_{0}=0.41R_{\rm a}$ simulation.

4.3.1 Forces on the Center of Mass

First, let us consider how the forces change the translational motion of the binary CM. At each time step, we calculate the DF and momentum transport forces acting on each particle. The DF forces on $m_{1}$ and $m_{2}$ are calculated according to Equation (28), while the momentum transport forces are calculated according to Equation (25). We can ask how these forces influence the motion of the CM by considering the component of each force along the CM velocity relative to the gas.

The velocity of the CM relative to the grid is ${\bf V}_{\rm CM}$ and the velocity of the wind relative to the grid is ${\bf v}_{\infty}=v_{\infty}\hat{x}$, so the velocity of the CM relative to the gas is ${\bf V}_{\rm rel}={\bf V}_{\rm CM}-{\bf v}_{\infty}$. Letting $\hat{V}_{\rm rel}$ be the unit vector in the direction of ${\bf V}_{\rm rel}$, we find the scalar component of the DF and momentum transport forces along ${\bf V}_{\rm rel}$ are, respectively,

[TABLE]

and the momentum transport force on the CM is

[TABLE]

Notice that $F_{\rm DF,CM}$ and $F_{\dot{p}\rm,CM}$ may be positive or negative. A positive value indicates a net thrust, while a negative value exerts a net drag. The net force on the CM

[TABLE]

is the total drag (or thrust) that results from the sum of the two projections.

Figure 9 shows the forces on the CM of the $a_{0}=0.41R_{a}$ binary from $t=30\,R_{\rm a}/v_{\infty}$ until the end of the simulation. The top panel shows the DF forces on the CM (Equation (31)). The contribution due to $m_{1}$ is ${\bf F}_{\rm DF_{1}}\cdot\hat{V}_{\rm rel}$, while the contribution from $m_{2}$ is ${\bf F}_{\rm DF_{2}}\cdot\hat{V}_{\rm rel}$. The periodic modulation of the blue and gray curves is due to the orbital motion of each particle. Even though the density structure at large scales is relatively stable, the position and velocity of each particle relative to the gas distribution change as the particles complete each orbit. This gives rise to a changing angle between each force vector and $\hat{V}_{\rm rel}$ as well as a changing magnitude of the force vectors themselves (${\bf F}_{\rm DF_{1}}$ and ${\bf F}_{\rm DF_{2}}$). In the same panel, $F_{\rm DF,CM}$ (in black) also shows variation due to the orbital motion about a nearly constant value of -4.2 $F_{\rm BH}$, which is less than the single-particle value of -5.1 $F_{\rm BH}$. Overall, the negative value of $F_{\rm DF,CM}$ indicates that the large-scale structure always exerts a drag on the CM motion of the pair.

The center panel of Figure 9 shows the momentum transport forces acting along the CM motion of the binary, which are calculated according to Equation (32). The $m_{1}$ and $m_{2}$ curves show ${\bf F}_{\dot{p}_{1}}\cdot\hat{V}_{\rm rel}$ and ${\bf F}_{\dot{p}_{2}}\cdot\hat{V}_{\rm rel}$, respectively. The blue and gray curves dip to negative values when accreted material is moving opposite the CM motion of the particle. Their sum is always positive, though, so the momentum transport force always exerts a net thrust on the CM. This is similar to the single-particle case, although the median value of $|F_{\dot{p},\rm CM}|\approx$ 0.18 $F_{\rm BH}$ is $74\%$ less than the magnitude of the single-particle value of 0.69 $F_{\rm BH}$. Finally, the bottom panel of Figure 9 shows the net force on the CM, given by Equation (33). The combination of DF and momentum transport forces always exerts a drag on the CM. The orbit-averaged net force on the CM for the $a_{0}=0.41R_{a}$ binary is -4.06 $F_{\rm BH}$, a value that is only $7\%$ lower than in the single-particle case.

4.3.2 Forces on the Orbital Motion

We now consider the forces affecting the motion of the binary about the CM. To do this, we find the projections of the forces along the barycentric velocity vectors of the two particles. The barycentric velocities of $m_{1}$ and $m_{2}$ are ${\bf V}_{1}$ and ${\bf V}_{2}$, respectively (Equation (19)). The unit vectors in the directions of these velocities are then $\hat{V}_{1}={\bf V}_{1}/V_{1}$ and $\hat{V}_{2}={\bf V}_{2}/V_{2}$. The DF force on the orbital motion of the binary is

[TABLE]

and the momentum transport force on the orbit is

[TABLE]

which gives a net “orbital” force of

[TABLE]

which is the force projected along the objects’ respective directions of orbital motion.

In Figure 10, we show these orbital forces calculated for $a_{0}=0.41R_{a}$. The DF forces, ${\bf F}_{\rm DF_{1}}\cdot\hat{V}_{1}$ and ${\bf F}_{\rm DF_{2}}\cdot\hat{V}_{2}$, shown in the upper panel of Figure 10, show a sinusoidal fluctuation between large negative and positive values. Because the overdense wake located at positive $x$ values remains the primary driver of DF forces, this periodic change in the sign of $F_{\rm DF,orb}$ is dominated by the changing direction of the barycentric velocities over each single orbit. The net force on the orbit due to DF (black curve in the upper panel) is small and negative, always exerting a small drag on the orbit.

The momentum accretion forces, ${\bf F}_{\dot{p}_{1}}\cdot\hat{V}_{1}$ and ${\bf F}_{\dot{p}_{2}}\cdot\hat{V}_{2}$, shown in the center panel of Figure 10, exhibit more complex behavior throughout the orbit. Consider the blue curve for $m_{1}$. The momentum accretion often causes a drag on the orbital motion of $m_{1}$. The most negative value (largest drag force) occurs while $m_{1}$ is moving in the same direction as the wind (for example, the bottom three panels of Figure 4). At some instances, the blue curve turns positive. This occurs when $m_{1}$ is traveling head on into the wind. At this time, $m_{1}$ is moving supersonically relative to the material behind the primary shock. The sharp (positive) peaks in the blue curve occur when a secondary shock forms ahead of $m_{1}$ (as in the first snapshot in Figure 4). In these cases, $m_{1}$ shocks and focuses material locally, accreting some of this material from downstream relative to its instantaneous motion. Because this material is moving in the same direction as $m_{1}$ as it is accreted, it imparts a thrust.

Although the momentum transport forces on $m_{1}$ and $m_{2}$ are sometimes positive, these forces always sum to a negative value (black curve in the middle panel of Figure 10). This reflects the fact that, in the case of $a_{0}=0.41R_{a}$, the highest accretion rates happen during the subsonic portions of the orbit of each particle. Material accreted from upstream exerts a drag force on each particle’s orbital motion, and the overall result of momentum accretion is a net drag on the orbit.

The net force on the orbit due to both momentum accretion and DF (bottom panel of Figure 10) is always negative (a drag). Momentum accretion contributes nearly half of the drag force on the orbit. We note that momentum accretion is a local phenomenon that depends on the details of the flow in the region of the orbiting binary. An understanding of the long-term evolution of the orbit of an astrophysical binary, therefore, requires an understanding of the dynamics and microphysics at small scales. This conclusion can be contrasted to the CM forces shown in Figure 9, in which the forces are dominated by the DF force, which results from the large-scale structure of the gas.

5 Analysis: Timescales for Binary Transformation

A binary interacting with an ambient gas will transform as a result of drag and accretion. The mass of the binary grows by accretion over a timescale $M/\dot{M}$. Forces on the CM reduce the CM kinetic energy, $E_{\rm trans}=M{\bf V}_{\rm CM}^{2}/2$, over a timescale $E_{\rm trans}/\dot{E}_{\rm trans}$. Forces on the binary orbit tighten the separation over a timescale $a/\dot{a}$. In this section we explore these timescales and their hierarchy. We first derive analytic results from BHL theory as context and then we discuss the numerical results from our simulation models.

5.1 Analytical Timescales

Before turning to our numerical simulations, let us see what the BHL formalism tells us about these timescales. The BHL stopping timescales are $\tau_{\rm stop,HL}$ and $\tau_{\rm stop,BH}$ (equations (6) and (13), respectively) with $M=m_{1}+m_{2}$. Similarly, the accretion timescales are $\tau_{\rm accr,HL}$ and $\tau_{\rm accr,BH}$, given by equations (7) and (14), respectively, with $M=m_{1}+m_{2}$.

We derive expressions for the BHL inspiral timescales in Appendix A. For large Mach numbers, the HL characteristic scalings give an orbital inspiral time that can be expressed in terms of the stopping time:

[TABLE]

In the regime of mild Mach numbers, the BH formalism gives an orbital inspiral time of

[TABLE]

Because $\tau_{\rm stop,HL}\propto\rho_{\infty}^{-1}$, the inspiral timescales are inversely proportional to the ambient density. If the density of the medium is increased, then the binary evolves more quickly.

A less intuitive feature of the two inspiral timescales is the leading-order dependence on $a^{-3/2}$. This dependence arises from the fact that the accretion radius of each particle is inversely proportional to the square of its speed relative to the gas. As the binary separation decreases, the orbital speed of each object increases, so the BHL force decreases with decreasing $a$. The lengthening of $\tau_{\rm insp,BH}$ with decreasing $a$, for fixed $\rho_{\infty}$, implies that the binary will never merge.

A natural question to ask is whether the stopping time is longer than the inspiral timescale for a given binary. In the BH regime, the ratio of the timescales is

[TABLE]

For $\mathcal{M}_{\infty}=2$, this implies that $\tau_{\rm stop,BH}>\tau_{\rm insp,BH}$ where $a\gtrsim 0.21R_{a}$ (this transition would occur at $a=0.198R_{a}$ if the HL timescales are used instead).

5.2 Numerical Timescales and Their Hierarchy

Having established the analytic predictions for timescales of accretion, stopping, and inspiral, we now proceed to measure these same quantities in our numerical simulations. For our simulated binaries, we measure orbit-averaged timescales of mass accretion

[TABLE]

slowdown of the CM

[TABLE]

and orbital inspiral

[TABLE]

by taking the median over an integer number of orbits from $t=30\,R_{\rm a}/v_{\infty}$ until the end of each simulation ($t=50\,R_{a}/v_{\infty}$). For $a_{0}=0.0R_{\rm a}$, the median and percentile ranges of $\langle{\tau_{\rm accr}}\rangle$ and $\langle{\tau_{\rm stop}}\rangle$ are computed over $30\leq t\leq 50\,R_{\rm a}/v_{\infty}$.

Numerically derived orbit-averaged timescales are shown in Figure 11. The difference in $\langle{\tau_{\rm accr}}\rangle$ and $\langle{\tau_{\rm stop}}\rangle$ at $a_{0}=0.0R_{\rm a}$ conveys the difference in efficiency of accretion and drag. The fact that $\langle{\tau_{\rm stop}}\rangle\ll\langle{\tau_{\rm accr}}\rangle$ implies that drag coefficients, compared to the BH predictions, are generally larger than unity, while accretion coefficients are of order unity or lower. This hierarchy of timescales has been consistently observed in previous simulations (e.g. Blondin & Raymer, 2012; Thun et al., 2016) and arises from large-scale gravitational interactions contributing to the DF drag (Ostriker, 1999), while accretion is a local property.

As we introduce binaries of increasing separation, the timescales change, but their basic hierarchy stays the same. The lengthening accretion timescale with larger binary separation arises from a drop in $\langle{\dot{M}}\rangle$ with larger $a_{0}$, as shown in Figure 8. We also see that $\langle{\tau_{\rm stop}}\rangle$ changes by only 5% as the binary separation increases. This reflects the fact that the slowdown of the CM is dominated by the DF force that is due to the downstream wake. The wider binaries show larger time-variability because the primary shock changes shape as the particles move in and out of each other’s wake (see Figure 1), but the overall structure downstream from the binary is similar in each simulation. Critically, the overall size of the wake and the standoff distance of the primary shock are similar in all simulations, implying the Coulomb logarithm, $\ln(r_{\rm max}/r_{\rm min})$ in the DF force is essentially unchanged (Ostriker, 1999; Thun et al., 2016, Section 5.3).

Figure 11 shows that the timescales fall into a clear hierarchy with $\langle{\tau_{\rm accr}}\rangle>\langle{\tau_{\rm insp}}\rangle>\langle{\tau_{\rm stop}}\rangle$. The most striking aspect of this hierarchy is that $\langle{\tau_{\rm insp}}\rangle>\langle{\tau_{\rm stop}}\rangle$ for all of the binary systems we simulated. In Section 5.1, we found that a simple application of BHL predicts that $\langle{\tau_{\rm stop}}\rangle$ should be longer than $\langle{\tau_{\rm insp}}\rangle$ for all binaries with $a_{0}\geq 0.21R_{a}$. Instead, the stopping time is always shorter, which means that the CM motion will come to a stop before the orbit decays.

Another prediction of the BHL formalism is the shape of $\langle{\tau_{\rm insp}}\rangle$ as a function of $a_{0}$. We plot $\tau_{\rm insp,BH}$ versus $a_{0}$ in Figure 12 (dashed line). The plot reveals the $a^{-3/2}$ leading-order shape of $\tau_{\rm insp,BH}$. We include in Figure 12 the calculated values of $\langle{\tau_{\rm insp}}\rangle$ of Figure 11. Instead of following $a^{-3/2}$, $\langle{\tau_{\rm insp}}\rangle$ decreases with decreasing $a_{0}$.

To obtain the $a_{0}$-dependence of $\langle{\tau_{\rm insp}}\rangle$, we perform a least-squares fit to a power law of the form

[TABLE]

We fit to the median values (the black points in Figure 12) and use one standard deviation (approximately the average of the upper and lower error bars) as the uncertainty on each point when performing the fit. We find $A=0.64\pm 0.01$ and $B=0.19\pm 0.01$, giving the functional form

[TABLE]

This fit is shown in Figure 12.

5.3 Physical Interpretation of $\langle{\tau_{\rm insp}}\rangle$

The sharp difference in scaling between $\langle{\tau_{\rm insp}}\rangle\propto a_{0}^{0.19\pm 0.01}$, computed from the hydrodynamic simulations, and $\tau_{\rm insp,BH}\propto a^{-3/2}$, from BHL theory, is apparent in Figure 12. While BHL theory predicts that the binary will take an infinite time to merge, $\tau_{\rm insp,fit}$ grows shorter as the orbit decays.

By assuming a background density of $\rho_{\infty}$, $\tau_{\rm insp,BH}$ fails to capture the enhancement of the local density that is due to the gravitational pull of the binary CM. If the CM of the binary were at rest relative to the gas, we would expect the density as a function of distance from the CM to be similar to a Bondi profile (Bondi, 1952, albeit truncated inside the orbit of the binary since the potential is smoothed to a distance of approximately $a/2$). We use a Bondi density profile to express the density of the gas at a distance of $a/2$ from the CM of the binary:

[TABLE]

(for $\gamma=5/3$ and $a/2\leq R_{B}$; Shapiro & Teukolsky, 1983). Replacing $\rho_{\infty}$ with $\tilde{\rho}$ in Equation (38) cancels the leading $a^{-3/2}$ dependence of $\tau_{\rm insp,BH}$. The modified timescale is then

[TABLE]

where the coefficient $\eta$ is an overall normalization factor. We choose the normalization such that $\widetilde{\tau}_{\rm insp,BH}=\langle{\tau_{\rm insp}}\rangle$ at $a_{0}=0.16R_{\rm a}$, which yields $\eta=4.96$.

We plot $\widetilde{\tau}_{\rm insp,BH}(\eta=4.96)$ in Figure 12. With this normalization, $\widetilde{\tau}_{\rm insp,BH}$ nearly replicates the power-law fit and the numerical data. The normalization factor of $\eta\approx 5$ indicates that Equation (45) overestimates the density by a factor of approximately 5. This may be because the softening of the potential to the size of the orbit and the pressure support provided by the orbital motion of the binary both work to reduce the local density that can be realized versus the pure Bondi case. In addition, the incoming wind sweeps material from the region of the binary, constantly modifying the density structure of the local region.

5.4 Comparison to Previous Work

In previous analytical and semianalytical work, there has been considerable disagreement of how the inspiral timescale and forces scale with binary separation. Here we briefly discuss these results in light of our semianalytical interpretation of the inspiral timescale of Section 5.3, noting the different physical scenarios and assumptions that led previous authors to varying results.

To simplify the comparison to previous studies, we first write our inspiral timescales in terms of the barycentric orbital Mach number

[TABLE]

which says

[TABLE]

Substituting the last expression into equations (38) and (46) yields the forms

[TABLE]

and

[TABLE]

respectively. We plot equations (47) and (48) along with our simulation data, $\langle{\tau_{\rm insp}}\rangle$, in Figure 13. In these units, the fit to our simulation data, Equation (44), becomes

[TABLE]

implying

[TABLE]

The gaseous DF force acting on a point mass in straight-line motion relative to a background gas was first calculated analytically by Ostriker (1999). As in the case of a collisionless gas, the integral nature of the DF force gives rise to a Coulomb logarithm $\ln(r_{\rm max}/r_{\rm min})$. The radii $r_{\rm max}$ and $r_{\rm min}$ represent the range of size scales spanned by the asymmetric wake trailing the point mass. Kim & Kim (2007) applied the methods of Ostriker (1999) to a perturber moving in a circular orbit relative to an initially uniform gas, finding that the orbital motion truncates the maximum coherence length of the wake to $r_{\rm max}\approx a$.

Kim et al. (2008) extended these calculations to a symmetric binary pair of perturbers in a circular orbit with CM at rest ($\mathcal{M}_{\infty}=0$). Kim & Kim (2007) and Kim et al. (2008) report their results with piecewise functions of $\mathcal{M}_{\rm orb}$111In their notation, the barycentric Mach number is $\mathcal{M}\equiv V_{p}/c_{s}$, where $V_{p}$ is the orbital speed of their perturber and $c_{s}$ is the unperturbed sound speed of the background gas. Their quantity $\mathcal{M}$ is identical to our $\mathcal{M}_{\rm orb}$. The difference in notation is to avoid confusion with $\mathcal{M}_{\infty}$, which is the CM Mach number of the binary. We note for convenience that they label the mass and orbit radius of the perturber as $M_{p}$ and $R_{p}$, respectively, which are $m_{1}=M/2$ and $a/2$ in our notation, respectively. The linear momentum of the perturber is $M_{p}V_{p}$, which is $Ma/4$ in our notation.. To compare their results to ours, we construct an inspiral timescale from their azimuthal drag forces as follows. The azimuthal DF force on one perturber due to its own wake is $\mathcal{F}\mathcal{I}_{\phi,1}$. The azimuthal DF force on that same perturber due to the companion’s wake is $\mathcal{F}\mathcal{I}_{\phi,2}$. Here, $\mathcal{F}$ and $\mathcal{I}_{\phi,1}$ are given by equations 11 and 14 of Kim & Kim (2007), respectively, and $\mathcal{I}_{\phi,2}$ by equation 5 of Kim et al. (2008). When computing $\mathcal{I}_{\phi,1}$ ,we take $r_{\rm min}=0.05R_{\rm a}$, which is the size of the sink in our simulations and, as such, is the minimum distance to which we integrate the DF force. The inspiral timescale is the linear momentum of the perturber divided by the total azimuthal force on the perturber, which simplifies to

[TABLE]

We plot $\tau_{\rm insp,K08}$ as a function of $\mathcal{M}_{\rm orb}$ in Figure 13.

Figure 13 shows that for $\mathcal{M}_{\rm orb}\lesssim 1.1$, $\tau_{\rm insp,K08}$ increases much more rapidly than $\langle{\tau_{\rm insp}}\rangle$. For $\mathcal{M}_{\rm orb}\gtrsim 1.1$, $\tau_{\rm insp,K08}$ traces $\tau_{\rm insp,BH}\propto\mathcal{M}_{\rm orb}^{3}$, which similarly assumes a uniform background density.

Stahler (2010) also considered angular momentum loss by a binary orbiting in a uniform gas with CM at rest in the linear regime, but came to a very different conclusion. They computed the angular momentum and energy carried by waves induced by the time-varying quadrupole moment of the binary potential. Their equation 39 gives the total torque, $\Gamma$, on a binary with moment of inertia $I$ and orbital frequency $\omega$. For comparison,222Stahler (2010) denotes the total mass and total semi-major axis of the binary as $M_{\rm tot}$ and $a_{\rm tot}$, respectively, which are $M$ and $a$ in our notation, respectively. In their equation 39, $\rho_{0}$ and $c_{s}$ are the unperturbed density and sound speed, respectively. In constructing $J/\Gamma$, we use the fact that, in our notation, $I\omega^{2}=Ma^{2}/4=M\mathcal{M}_{\rm orb}^{2}c_{s,\infty}^{2}$ for an equal-mass binary., we construct an inspiral timescale by dividing the total angular momentum of the binary, $J=I\omega$, by their expression for $\Gamma$. The resulting timescale

[TABLE]

is shown in green in Figure 13. It is apparent that the scaling with $\mathcal{M}_{\rm orb}$ for $\mathcal{M}_{\rm orb}\gtrsim 1.1$ is opposite in sign from both $\tau_{\rm insp,K08}$ and $\tau_{\rm insp,BH}$. Sánchez-Salcedo & Chametla (2014) report on this discrepancy and demonstrate that it can be traced to the multipole expansion that underlies the solution from Stahler (2010).

Similar to Kim & Kim (2007), Sánchez-Salcedo & Chametla (2014) calculated the torque on a single circular-orbit perturber in the linear regime, but they allowed the orbit center to move at ${\bf V}_{\rm CM}$ relative to the initially uniform background gas. They considered the cases of face-on orbits, in which ${\bf V}_{\rm CM}$ is parallel to the orbital angular momentum vector, as well as edge-on cases (as in the present work), in which the orbital angular momentum vector is perpendicular to ${\bf V}_{\rm CM}$. Their Figures 9 (face-on cases) and 12 (edge-on cases) present computed values of the azimuthal DF force on a single perturber due to its own wake, denoted $F_{\rm DF,\phi}$, for different combinations of $\mathcal{M}_{\infty}$ and $\mathcal{M}_{\rm orb}$. Although their method is slightly different, the DF torque on their face-on perturber with subsonic CM motion ($\mathcal{M}_{\infty}=0.5$) tends to the $\mathcal{M}_{\infty}=0$ torque calculated by Kim & Kim (2007), with $\tau_{\rm insp}\sim\tau_{\rm insp,BH}$ for $\mathcal{M}_{\rm orb}\gtrsim 1.1$.

By allowing for CM motion, the calculations of Sánchez-Salcedo & Chametla (2014) are most applicable to our particular study of $\tau_{\rm insp}$. To apply their values of $F_{\rm DF,\phi}$ to a perturber in an equal-mass binary system, the contribution of the companion’s wake to the azimuthal force must be included. The total azimuthal force on the perturber (due to both wakes) is $F^{(1)}_{\rm DF,\phi}=\mathcal{R}_{\phi}F_{\rm DF,\phi}$ where $0.1\lesssim\mathcal{R}_{\phi}\lesssim 1$ (see their §5.1). Only a few values of $\mathcal{R}_{\phi}$ are given, but Kim et al. (2008) showed that $\mathcal{R}_{\phi}$ is a function of $\mathcal{M}_{\rm orb}$. Without an understanding of how $\mathcal{R}_{\phi}$ scales with $\mathcal{M}_{\rm orb}$, we are unable to include results from Sánchez-Salcedo & Chametla (2014) in Figure 13, although such a comparison is highly desirable.

Our numerical result exhibits a key difference relative to each of these previous works, which all assume a constant-density background. We can improve the analytic representation to bring it into agreement with the numerical findings by considering the enhanced density that arises from flow convergence onto the binary pair. Modifying $\tau_{\rm insp,BH}$ by adopting a local density, $\tilde{\rho}\propto r^{-n}$ (where $n=1.5$ in the BHL case; Section 5.3), changes the scaling of the analytic predictions and can reconcile the results of the analytic models with our simulations, as shown in Figure 12.

6 Discussion

In this section, we discuss some implications of our simulation results for astrophysical binaries.

6.1 Implications of Timescale Hierarchy and Scaling

In our models, we focused on the ratio of the initial separation of the binary, $a_{0}$, relative to the accretion radius of its center of mass, $R_{a}$, as an essential parameter in describing the flow. This ratio informs the rate of mass accretion, and it determines the nature of the drag forces acting on both the orbital and translational motion of the pair. The timescales for accretion, $\tau_{\rm accr}$, orbital inspiral, $\tau_{\rm insp}$, and slowdown of the center of mass, $\tau_{\rm stop}$ inherit this dependence on $a_{0}/R_{a}$. The timescales follow a hierarchy, with $\tau_{\rm accr}>\tau_{\rm insp}>\tau_{\rm stop}$ for all $a_{0}/R_{\rm a}$ that we simulated. In contrast to the expectation from BHL, $\tau_{\rm insp}$ is a monotonically increasing function of $a_{0}$.

This hierarchy of timescales implies that astrophysical binaries moving through gaseous surroundings will slow relative to the surrounding gas before they either inspiral or accrete significant quantities of gas. Binaries will be slowed when their gravitational cross section has intercepted a gas mass of about 15% of the binary mass (because $\tau_{\rm stop}\approx 0.15\tau_{\rm stop,BH}$). To be driven toward merger, they would need to traverse a larger column of approximately 50% their mass, and to double in mass, they need to intercept nearly twice their own mass.333We note, however, that the numerically derived accretion rate (and therefore timescale) is sensitive to the size of the absorbing sink boundary condition, such that smaller objects would be expected to accrete less material (Appendix B). Further, this mass-doubling timescale assumes that the object is able to accrete efficiently from any gas that falls through the sink boundary without, for example, forming a disk or feeding back substantially into its surroundings. This suggests that a binary in motion through a uniform gaseous environment would be slowed such that it would enter into a phase of binary Bondi accretion in which the binary CM was approximately at rest relative to the gas.

The scaling of $\tau_{\rm insp}$ with binary separation is also of great importance. Because this timescale has a positive slope, binaries with shrinking $a_{0}/R_{\rm a}$ continue to be driven to merger on finite timescales. The ratio $a_{0}/R_{\rm a}$ can shrink either because the binary separation is decreasing or as $R_{\rm a}$ increases when the binary CM motion is slowed. Therefore, whereas the naive scaling $\tau_{\rm insp,BH}$ indicates that pairs of objects can rarely be driven fully to merger by gas – instead being halted at finite separation by the increasing timescale as their orbits shrink – our updated results indicate that a tightening orbit yields a trend toward ever-shorter $\tau_{\rm insp}$.

6.2 Critical Separation between Gas-dominated and GW-dominated Inspiral

The fact that gaseous forces drive binaries to tighter separations leads to the conclusion that GWs may become important at some separation, especially if the binary is composed of compact objects. For a circular-orbit, equal-mass binary, Peters (1964) showed that the inspiral timescale is

[TABLE]

The inspiral timescale that we found numerically, $\tau_{\rm insp,fit}$ given in Equation (44), says that $\tau_{\rm insp,gas}\propto a_{0}^{0.19}$. This mild, positive power law allows the gas to be the dominant driver of orbital inspiral at large separations. At smaller separations, $\tau_{\rm insp,GW}<\tau_{\rm insp,gas}$ and GW radiation dominates the orbital decay until the binary finally merges. Below we find the boundary between gas-dominated and GW-dominated inspiral, by computing the critical separation, $a_{0,\rm crit}$, at which $\tau_{\rm insp,gas}=\tau_{\rm insp,GW}$.

Our expression for $\tau_{\rm insp,fit}$ in Equation (44) is given in units of $\tau_{\rm stop,BH}$. To simplify the algebra, let us write our fit in terms of $\tau_{\rm stop,HL}$, instead. Let us define

[TABLE]

where, in these units, $A=0.71$ instead of $0.64$. The power $B=0.19$ is unchanged. In astrophysical units,

[TABLE]

Equating $\tau_{\rm insp,gas}$ and $\tau_{\rm insp,GW}$ gives the critical separation

[TABLE]

with $A=0.71$ and $B=0.19$. In astrophysical units,

[TABLE]

In Figure 14, we plot $a_{0,\rm crit}$ versus number density, $n_{\infty}$, for three sample binary systems (black lines), where we have converted mass density to number density through $\rho_{\infty}=m_{p}n_{\infty}$. The vertical shaded regions give the typical densities encountered in several astrophysical environments. The ISM densities use values from Kim & Ostriker (2017), the AGN conditions use the disk models of Bartos et al. (2017), and the CE density range is modeled after the red giant proposed by Sabach & Soker (2015) to have initiated a CE phase in the hierarchical triple system PSR J0337+1715 . The gray hatched region shows where the enclosed gas mass becomes similar to the binary mass ($M_{\rm encl}>2M_{\sun}$). For higher-mass binaries, the lower boundary of the hatched region moves up and to the right.

Each black line in the figure plots Equation (53) for a different combination of $M$ and $v_{\infty}$. Above each line, $\tau_{\rm insp,gas}$ is faster than $\tau_{\rm insp,GW}$ and we say that the binary is in the gas-dominated regime. Below each line, $\tau_{\rm insp,GW}<\tau_{\rm insp,gas}$, and the binary is in the GW-dominated regime.

Below, we consider some example systems in the context of Figure 14. For convenience, we also plot Equation (52) for two different choices of $M$ and $v_{\infty}$ in Figure 15.

6.3 Binaries in the ISM

The multiphase ISM spans a large range of density and temperature conditions. Here, for the sake of a concrete example, we focus on the warm ISM, which occupies about one-half of the interstellar volume, has a number density of approximately 1 cm*-3*, and has a temperature of the order of $10^{4}$K. A binary traveling at $10$ km s*-1* relative to the gas has a Mach number of 1.1. At this speed, a $2M_{\odot}$ binary with a separation of 1 au has $\tau_{\rm insp,gas}<\tau_{\rm insp,GW}$. For this system, our numerical fits of $\langle{\tau_{\rm insp}}\rangle$ and $\langle{\tau_{\rm stop}}\rangle$ (equations (44) and (C2), respectively) give a stopping timescale of $2\times 10^{14}$ yr and an orbital inspiral timescale of $5\times 10^{13}$ yr (see Figure 15). These timescales indicate that binaries traversing the warm ISM will suffer little dynamical effect from this interaction over the age of the universe.

Nonetheless, this effect might potentially be measurable in the case of pulsar binaries with extremely precise timing. Using $2\dot{P}_{\rm orb}/{P}_{\rm orb}=3\dot{a}/{a}=3/\tau_{\rm insp}$, the inspiral time gives a period derivative of $\dot{P}_{\rm orb}\approx 10^{-13}$.444In this estimate, we have neglected the pulsar’s influence on its surroundings, which may not be justified because many pulsars are surrounded by hot, low-density wind nebulae. In cases where the standoff distance of the cavity is similar to the binary separation or bow shock size, drag forces on the binary orbit would likely be substantially reduced. By comparison, for the binary radio pulsar PSR J0751+1807, Nice et al. (2008) found an orbital period derivative of $(-3.1\pm 0.5)\times 10^{-14}$, a precision that is two orders of magnitude better than our sample case presented here.

6.4 Binaries Interacting with AGN Disks

The central regions of galaxies hosting AGNs experience an inflow of gas that settles into a cold accretion disk surrounding the central SMBH. Stars and binaries coexist with the SMBH and accretion disk in these central regions, forming a dense cluster around the black hole. Dynamical interactions between the stars and disk can trap stars within the disk plane (e.g. Syer et al., 1991; Artymowicz et al., 1993; Ivanov et al., 1998; Šubr & Karas, 2005; Miralda-Escudé & Kollmeier, 2005; Karas & Šubr, 2007; Perets et al., 2007; Baruteau et al., 2011; McKernan et al., 2011, 2012; Just et al., 2012; del Valle & Escala, 2014; Kennedy et al., 2016; Panamarev et al., 2018).

Our result that $\tau_{\rm stop}<\tau_{\rm insp}$ for all of the binaries we simulated demonstrates that trapping by the AGN disk is not reserved for single stars. Instead, the CM motion of a binary can be slowed by the gaseous disk without causing a merger of the binary (we expect this conclusion to hold as long as $a_{0}<H$, the scale height of the disk). That is, binaries formed outside of the disk can remain intact after undergoing capture into the disk plane. When this process involves BBHs, the gaseous drag can aid the trend toward merger until the separation is reduced sufficiently that GWs take over, as discussed in Section 6.2. These mergers are of additional interest because BBH mergers in dense, gaseous environments may be accompanied by detectable electromagnetic emission (e.g. Farris et al., 2015; Stone et al., 2017; Kelly et al., 2017).

To give a concrete example, we use the disk model of Bartos et al. (2017, their Equations (1)-(4)) with a central black hole mass of $M_{\bullet}=10^{6}M_{\sun}$. The central black hole dominates the gravitational force within a sphere of influence of $\approx 1$ pc, so we assume that binaries orbiting the black hole in this region do so in Keplerian orbits with a CM speed of $v_{\infty}=\sqrt{GM_{\bullet}/R_{\rm disk}}$. At $R_{\rm disk}=0.1$ pc, $n_{\infty}=3\times 10^{11}$cm*-3*, and the Keplerian speed is $v_{\infty}=207$ km s*-1*. This speed gives a CM Mach number of 4.6 (neglecting the effect of rotation of the disk). A $20M_{\sun}$ BBH with $a_{0}=2$ AU falls in the gas-dominated regime (see Figure 14). If inclined such that it is continuously within the disk plane, we find $\tau_{\rm insp}=2\times 10^{5}$ yr (see Figure 15), $\tau_{\rm stop}=4\times 10^{4}$ yr, and $\tau_{\rm accr}=5\times 10^{5}$ yr. This implies that the binary is entrained within the disk flow (and any eccentricity is damped) within approximately 14 orbits, and the BBH inspirals within approximately 70 orbits. Over this time, the BBH might accrete up to 40% its own mass, growing to approximately 30$M_{\odot}$. The brevity of these timescales relative to a typical AGN activity cycle, on the order of $10^{6}$ yr, indicates that under certain conditions AGN disks very likely play a role in tightening BBHs.

6.5 Embedded Binaries in the CE

Triple-star systems are common. While about 10% of solar-mass stars are in triple or higher-order systems, nearly 60% of O-type stars are in similar hierarchical configurations (e.g. Moe & Di Stefano, 2017, Figure 39). For long-term dynamical stability, these systems must be hierarchical (Naoz & Fabrycky, 2014). Many such triple systems containing compact objects must have experienced significant mass-transfer events as their stellar progenitors evolved (Tauris & van den Heuvel, 2014; Sabach & Soker, 2015; Toonen et al., 2016; Di Stefano, 2018). In some of these hierarchical triple systems, the outer star evolves to transfer mass toward an inner, more compact binary. Depending on the stability of this mass transfer, the donor star’s mass transfer could either proceed at a low rate or lead to runaway orbital tightening and the engulfment of the more compact pair within the envelope of the donor star. What happens when a close binary becomes embedded in a CE has been a topic of previous uncertainty (e.g. Sabach & Soker, 2015). In particular, would the binary spiral inward within the stellar envelope, or be driven to merger?

Unraveling the outcomes of these triple-star CE encounters requires an understanding of gas accretion and orbital evolution as studied in this paper. If one assumes a simple formalism, such as the timescales derived from BHL, it would appear that a binary embedded in a CE should merge faster than its CM spirals into the donor star’s envelope. We have shown here that this is not necessarily the case. Instead, we have found that the slowdown of the binary’s CM (which is the driver of the CE inspiral) occurs more quickly than the embedded binary’s own separation decays. Therefore, the embedded binary can spiral into the envelope gas, and that process alone will not cause the binary to merge.

Since the binary may be able to complete part of the CE inspiral intact, there are many possible outcomes for the interaction. The binary may survive the CE if the envelope can be ejected early in the encounter, the binary may be tidally disrupted before envelope ejection, or, in the case of unsuccessful envelope ejections, one or both of the binary components could merge with the core of the donor star. The revised formalism presented in this paper should be used when examining the various outcomes expected from these complex interactions. For example, Sabach & Soker (2015) considered a formation scenario for PSR J0337+1715 in which the outer binary system remains intact as its CM spirals into the envelope of the donor star, reaching high-enough densities to be tidally disrupted. Our results indicate this channel is indeed plausible, as the timescale for the drag-induced inspiral of the CM of the embedded binary is shorter than the timescale for the embedded binary’s finite semi-major axis to decay.

7 Conclusions

This paper has studied the evolution of binaries moving through gaseous media as a function of the binary system’s initial semi-major axis, $a_{0}$, relative to its gravitational focusing radius, $R_{a}$. We extended the analytical and numerical formalisms of single-particle BHL accretion to the case of a binary system moving perpendicular to the orbital angular momentum vector of the pair. We performed 3D hydrodynamical simulations to study the effects of accretion and drag on the translational and orbital motion of the system. Below we summarize the main conclusions of this work:

The introduction of a finite binary separation softens the gravitational potential of the pair relative to a single object and introduces time-dependent structures into the flow (as described in Section 4.1 and seen in Figures 1-4). We find that the ratio of the binary separation to the gravitational focusing scale is crucial in describing the flow structure. When $a/R_{\rm BH}<1$, the binary orbits within a bow shock that is otherwise similar to that from a single mass. In these cases, the orbital motion is supersonic relative to the postshock gas and the pair launches spiral shocks into its surroundings (Figure 2). By contrast, when $a/R_{\rm BH}>1$, the binary members are each trailed by their own shock structure. 2. 2.

The orbiting binary’s motion introduces periodic variation in the accretion rate onto the binary components (see Section 4.2 and Figure 5). The overall accretion rate onto the binary is again modulated by the ratio of the binary separation to the gravitational focusing scale. When $a/R_{\rm BH}\ll 1$, the accretion rate onto the binary approaches that of a single particle, with scaling $\dot{M}\propto M^{2}$. When $a/R_{\rm BH}>1$, the accretion rate approaches the sum of two independent particles with scaling $\dot{M}\propto m_{1}^{2}+m_{2}^{2}$ (Figure 8). 3. 3.

Gaseous DF and the accretion of momentum give rise to forces that modify the CM and orbital motion of the pair. In Section 4.3, we discuss the projection of these forces onto the CM and orbital velocities. CM deceleration is dominated by DF forces from the large-scale wake that trails the binary (Figure 9). Gas structures and flow on the scale of the binary orbit give rise to forces that drive orbital tightening (e.g., Kim et al., 2008). In this case, momentum accretion and DF become similarly important (Figure 10). 4. 4.

We measure orbit-averaged timescales for accretion, $\langle{\tau_{\rm accr}}\rangle$, slowdown of the CM, $\langle{\tau_{\rm stop}}\rangle$, and orbital inspiral of the binary, $\langle{\tau_{\rm insp}}\rangle$, with Equations (40)-(42), respectively. For all values of $a_{0}$ that we simulated, the timescales exhibit a clear hierarchy with $\tau_{\rm stop}<\tau_{\rm insp}<\tau_{\rm accr}$ (Figure 11). Therefore, a binary moving through a gaseous medium slows before being driven to merger or gaining significant mass (Section 5). 5. 5.

BHL theory predicts that $\tau_{\rm insp}\propto a_{0}^{-3/2}$, which implies that as the orbit tightens, the time to coalescence lengthens. By contrast, in our simulations, we find a mild, positive scaling, $\langle{\tau_{\rm insp}}\rangle\propto a^{0.19}$ (given by Equation (44) and shown in Figure 12). Previous semianalytic work predicting the sign and scaling of this relation has arrived at a variety of conclusions, as described in Section 5.4. We show that the simulation results can be reproduced by considering the enhanced density in which the binary orbits that is due to the convergence of gravitationally focused gas (Equation (46)), a key effect that was neglected in previous analytic studies. 6. 6.

The magnitude of the forces we derive implies that binaries in environments from the ISM to star-forming clouds, AGN disks, and stellar envelopes during CE phases all may experience transformation due to interaction with their gaseous surroundings. We discuss particular systems and the implications of the observed hierarchy of timescales, and we derive the transition between systems whose orbital evolution is dominated by gaseous interaction (at large separations) and those whose orbital evolution is driven by gravitational radiation (at small separations) in Section 6.

There are many aspects of this problem left to explore. In this work, we have simulated equal-mass binary systems in circular orbits inclined such that the incoming wind hits the orbital plane edge-on. A caveat in the extension of our results to astrophysical binaries is the assumption of a uniform background medium, which is valid when velocity and density scale heights are larger than both $a$ and $R_{a}$ of the binary and when conditions change over timescales longer than the orbital period.

Realistic binary systems are often subject to nonuniform gases with gradients in density, pressure, and Mach number. Previous work has shown that breaking the symmetry of canonical BHL has significant effects on the flow. For example, MacLeod & Ramirez-Ruiz (2015) and MacLeod et al. (2017) found that in single-object BHL, the presence of density gradients typical of CE encounters impart a net rotation and reduce accretion from postshock flow. Recently, Xu & Stone (2019) have examined the criteria that define how the inhomogeneity of the background gas affects the stability of postshock flow. Krumholz et al. (2006) simulated BHL flows in cold, turbulent media of relevance to star-forming clouds and found that the vorticity imparted by turbulence can dominate the flow structure and lead to highly variable accretion rates.

Future simulation work could relax our current assumptions of equal mass, circular orbits, edge-on inclination, and a homogeneous background medium to extend the applicability of these results to a broader range of astrophysical systems and environments.

We gratefully acknowledge the anonymous referee and I. El Mellah for helpful suggestions on an earlier draft of this manuscript. We thank Z. Haiman, M. Hance, D. Lee, D.N.C. Lin, J. Naiman, E. Ostriker, A. Rosen, and S. Schrøder for helpful discussions. The calculations for this research were carried out in part on the UCSC supercomputer Hyades, which is supported by the National Science Foundation (award number AST-1229745) and UCSC. A.A. is grateful for support from the NSF REU program Lamat at UCSC and the California Space Grant Consortium. M.M. is grateful for support for this work provided by NASA through Einstein Postdoctoral Fellowship grant No. PF6-170155 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. E.R.-R. acknowledges support from NASA ATP grant NNX14AH37G, NSF grant AST-1615881, the Julie Packard Foundation, the Radcliffe Institute for Advanced Study, the Heising-Simons Foundation, and the Danish National Research Foundation (DNRF132).

Appendix A Derivation of Analytical Timescales for Orbital Inspiral

A principal objective of this paper is to learn how the predictions of BHL must be modified to describe a binary system. In this section, we derive the analytic timescales for orbital inspiral given by equations (38) and (37).

To associate a timescale with the decay of the orbit, we imagine the CM at rest and calculate the Bondi-Hoyle energy dissipation rate for $m_{1}$ and $m_{2}$ as they orbit the CM. For an equal-mass binary in a circular orbit, $m_{1}=m_{2}=M/2$, and the orbital speed is $v_{\rm orb}^{2}=GM/a$. If $V_{1}$ and $V_{2}$ are the barycentric speeds of $m_{1}$ and $m_{2}$, respectively, then $V_{1}=V_{2}=v_{\rm orb}/2$. The BH accretion radius, Equation (10), of $m_{1}$ is

[TABLE]

with an energy dissipation rate, Equation (12), of

[TABLE]

The total energy dissipation rate for $m_{1}$ and $m_{2}$ is just twice $\dot{E}_{\rm BH,1}$, that is,

[TABLE]

Replacing $m_{1}=M/2$ and $V_{1}=v_{\rm orb}/2$ and factoring out a 1/4 from the denominator give

[TABLE]

The total orbital energy of the binary is $E_{\rm orb}=-Gm_{1}m_{2}/2a=-GM^{2}/8a$, so the timescale for orbital decay is

[TABLE]

If we factor $v_{\rm orb}^{2}$ out of the parentheses and substitute $v_{\rm orb}=(GM/a)^{1/2}$ then

[TABLE]

Now the term $c_{s,\infty}^{2}/GM$ can be expressed as

[TABLE]

so we can write

[TABLE]

It is convenient to introduce a factor of $R_{\rm a}$ in order to obtain the quantity $(a/R_{\rm a})$, giving

[TABLE]

The first term in brackets is simply $\tau_{\rm stop,HL}$, so we have arrived at Equation (38):

[TABLE]

In the Hoyle-Lyttleton formalism, the sound speed is omitted from the characteristic velocity in the expressions (A1) and (A2), resulting in a total orbital energy dissipation rate of $\dot{E}_{\rm HL,orb}=2\pi G^{2}M^{2}\rho_{\infty}v_{\rm orb}^{-1}$ and a decay timescale of

[TABLE]

which is Equation (37). This quantity is the leading coefficient of $\tau_{\rm insp,BH}$. That is,

[TABLE]

Appendix B Validation of Numerical Method

In this section, we characterize the sensitivity of our results to background density, $\rho_{\infty}$, the length of the computational domain in each Cartesian direction, $L$, the radius of the absorbing sinks used to represent the binary, $R_{\rm s}$, and the length of each side of the finest grid cells, $\delta_{\rm min}$. The studies presented in this section adopt $\mathcal{M}_{\infty}=2$, $\gamma=5/3$, and $a_{0}=1.0R_{a}$. In the panels of Figure 16, we plot $\langle{\tau_{\rm accr}}\rangle$, $\langle{\tau_{\rm insp}}\rangle$, and $\langle{\tau_{\rm stop}}\rangle$, normalized to $\tau_{\rm insp,BH}$, for each series of studies.

B.1 Grid Refinement

We maximally refine the region of our two sink particles so that sink cells have length of $\delta_{\rm min}$ in each Cartesian direction. We test sensitivity to $\delta_{\rm min}$ by running simulations with $R_{\rm s}/\delta_{\rm min}=$ 3 and 6 with sink sizes of $R_{\rm s}=0.05R_{\rm a}$ and $R_{\rm s}=0.1R_{\rm a}$. Comparing points with the same $R_{\rm s}$ in panel (C) of Figure 16 reveals that there is very little difference in the orbit-average timescales when the resolution is doubled from our fiducial value of $R_{\rm s}/\delta_{\rm min}=$ 3. When $R_{\rm s}=0.05R_{\rm a}$, $\langle{\tau_{\rm stop}}\rangle$ and $\langle{\tau_{\rm accr}}\rangle$ change by less than 1% and $\langle{\tau_{\rm insp}}\rangle$ increases by 2.5%. Accretion rates are larger when $R_{\rm s}=0.1R_{\rm a}$ and, in this case, we find that a resolution of twice our fiducial value increases $\langle{\tau_{\rm accr}}\rangle$ by $2.6\%$ and the other timescales change by less than 1%. We adopt $R_{\rm s}/\delta_{\rm min}=3$ in our production runs.

B.2 Sink Radius

The influence of sink radius, $R_{\rm s}$, on accretion rates in single-object BHL flows has been well studied (e.g. Ruffert & Arnett, 1994; Blondin & Raymer, 2012; MacLeod & Ramirez-Ruiz, 2015; MacLeod et al., 2017). For example, Ruffert & Arnett (1994) found that doubling the size of the sink increased accretion rates by $\approx 15\%-25\%$ in their simulations with $\mathcal{M}_{\infty}=3$ and $\gamma=5/3$. We test $R_{\rm s}$ in our setup by running simulations with half and twice our fiducial value of $R_{s}=0.05R_{\rm a}$. These simulations adopt $R_{\rm s}/\delta_{\rm min}=3$, $\rho_{\infty}=10$ and $L=32R_{a}$. When $R_{\rm s}/R_{\rm a}$ is doubled from 0.05 to 0.1, $\langle{\tau_{\rm stop}}\rangle$ decreases by 1% and $\langle{\tau_{\rm insp}}\rangle$ and $\langle{\tau_{\rm accr}}\rangle$ each decrease by 25%. A doubling of $R_{\rm s}/R_{\rm a}$ from 0.025 to 0.05 similarly decreases $\langle{\tau_{\rm stop}}\rangle$ by 1%, while $\langle{\tau_{\rm insp}}\rangle$ decreases by 17% and $\langle{\tau_{\rm accr}}\rangle$ decreases by 27%.

The filled points in panel C of Figure 16 show that the hierarchy of timescales $\langle{\tau_{\rm accr}}\rangle\geq\langle{\tau_{\rm insp}}\rangle\geq\langle{\tau_{\rm stop}}\rangle$ is preserved across values of $R_{\rm s}$. However, because sink size influences rates of mass accretion, the accretion timescale is sensitive to this choice. For an astrophysical accretor of size less than $R_{\rm s}$, the mass accretion rate we derive therefore represents the upper limit defined by the flux of material into the object’s vicinity, rather than the precise amount that is able to accrete. Because sink size affects the momentum accreted, we note that the exact scaling of $\tau_{\rm insp,fit}$ as a function of $a_{0}$ will be influenced by the choice of sink radius. We adopt $R_{\rm s}=0.05R_{\rm a}$ in our production runs.

B.3 Background Density

We performed tests of background density, $\rho_{\infty}$, to ensure that the timescales we derive could be scaled to any astrophysical density. We ran simulations with $\rho_{\infty}=$ 10, 100, and 1000. These simulations were run with $R_{\rm s}=0.05R_{a}$, $L=16R_{a}$, and $R_{\rm s}/\delta_{\rm min}=6$. Panel (A) of Figure 16, shows that the background density scales out of the orbit-average timescales so that the values are the same when normalized to $\tau_{\rm stop,BH}$. In all cases, $\langle{\tau_{\rm accr}}\rangle\gg\langle{\tau_{\rm insp}}\rangle\gg\langle{\tau_{\rm stop}}\rangle$. Increasing our fiducial value by a factor of 100 (to $\rho_{\infty}=1000$) only changes $\langle{\tau_{\rm accr}}\rangle/\tau_{\rm stop,BH}$ and $\langle{\tau_{\rm insp}}\rangle/\tau_{\rm stop,BH}$ by $7.6\%$ each and decreases $\langle{\tau_{\rm stop}}\rangle/\tau_{\rm stop,BH}$ by $4\%$. We adopt $\rho_{\infty}=10$ for our production runs.

B.4 Extent of the Computational Domain

The DF drag forces exerted on an object depend on the size of the wake, with approximate scaling dictated by a Coulomb logarithm $\ln(r_{\rm max}/r_{\rm min})$, where $r_{\rm min}$ is the size scale outside of which a net force is contributed and $r_{\rm max}$ is the system size (Ostriker, 1999). While we are forced to choose a single system size in our numerical simulations, we measure the approximate scaling here.

In our setup, the CM of the binary is located at the origin and we place all outer boundaries at a distance $L/2$ from the origin. We run simulations with half and double a fiducial value of $L=16R_{\rm a}$. Panel (B) of Figure 16 shows that $\langle{\tau_{\rm accr}}\rangle$ and $\langle{\tau_{\rm insp}}\rangle$ do not vary significantly with changing $L$. Doubling $L$ from 16 to 32 reduces $\langle{\tau_{\rm accr}}\rangle$ by 6.4% and reduces $\langle{\tau_{\rm insp}}\rangle$ by 2.6%.

On the other hand, the same doubling of $L$ reduces $\langle{\tau_{\rm stop}}\rangle$ by 43%, which is expected because this timescale is dominated by DF due to gas at larger scales. In otherwise similar simulations of single-object BHL, Thun et al. (2016) and MacLeod et al. (2017) found $r_{\rm min}\approx 0.5R_{\rm a}$ and $r_{\rm max}\approx L/2$. A fit of the DF force for the simulations presented here reveals that $r_{\rm min}\sim 3R_{\rm a}$ for the larger box sizes of $L=16R_{\rm a}$ and $L=32R_{\rm a}$ (adopting $r_{\rm max}\approx L/2$). Since $\tau_{\rm stop}\propto 1/F_{\rm DF,CM}$, our measured stopping time can be approximately scaled to larger astrophysical systems by dividing by $\ln[(L/2)/(3R_{\rm a})]$.

Appendix C Fitting Formulae to Simulation Results

In this section, we provide formulae for $\langle{\tau_{\rm stop}}\rangle$, $\langle{\tau_{\rm accr}}\rangle$, and $\langle{\dot{M}}\rangle$ as a function of $a_{0}$. The fit for $\langle{\tau_{\rm insp}}\rangle$ was given in Section 5.2. The remaining quantities $\langle{\tau_{\rm stop}}\rangle$, $\langle{\tau_{\rm accr}}\rangle$, and $\langle{\dot{M}}\rangle$ do not follow a single power law, so we fit piecewise functions for $a_{0}\leq 0.55$ and $a_{0}>0.55$ of the functional form

[TABLE]

which goes to a constant value because the two particles become independent as $a_{0}\to\infty$. We perform least-squares fits of the median values given in Table 2 to the model, $f(a_{0})$, weighting the points according to one standard deviation of their time variability. For$\langle{\tau_{\rm stop}}\rangle$, we find $A=0.009$ and $B=0.142$, giving

[TABLE]

in units of $\tau_{\rm stop,BH}$. For $\langle{\tau_{\rm accr}}\rangle$, we obtain $A=1.379$ and $B=1.265$, or

[TABLE]

in units of $\tau_{\rm stop,BH}$. For the accretion rate, $\langle{\dot{M}}\rangle$, in units of $\dot{M}_{\rm BH}$, we find $A=-0.724$ and $B=0.992$, that is,

[TABLE]

In panels (A), (B), and (C) of Figure 17, we plot $\langle{\tau_{\rm stop}}\rangle$, $\langle{\tau_{\rm accr}}\rangle$, and $\langle{\dot{M}}\rangle$, respectively, with points and error bars. The associated fit of each quantity is shown with a gray line in each panel.