Relativistic Dynamics and Mass Exchange in Binary Black Hole Mini-disks

Our goal is to study the hydrodynamics of mini disks when the BBH separation is in the relativistic regime. As a first step, we limit the problem to examining 2D inviscid hydrodynamics in the orbital plane. To approach a stationary state as rapidly as possible, we choose initial conditions that are nearly hydrostatic and in which the mini-disks extend close to (either near or somewhat outside) the tidal truncation radius estimated from Newtonian orbital mechanics. In addition, to disentangle purely gravitational effects from fluid effects, we will complement these hydrodynamic simulations with test-particle simulations. The spacetime itself is described by a high-order PN scheme in which gravity is due entirely to the two BHs of the binary. In plausible SMBBH circumstances, the gas mass is a tiny fraction of the binary mass, so its contribution to gravity should always be negligible.

Throughout the paper, unless otherwise noted, we use geometrized units in which $G=c=1$. When used as tensorial indices, we reserve Greek letters (e.g., $\alpha ,\beta ,\gamma ,\ldots$) for spacetime indices and Roman letters (e.g., $i,j,k,\ldots$) as indices spanning spatial dimensions.

We use Harm3d (Noble et al. 2009) to solve the equations of general relativistic hydrodynamics on a background spacetime in flux-conservative form. These equations amount to conservation of baryon number density and conservation of stress-energy (see Noble et al. 2009 for more details). Taken together, they can be written as

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{U}}({\boldsymbol{P}})=-{\partial }_{i}{{\boldsymbol{F}}}^{i}({\boldsymbol{P}})+{\boldsymbol{S}}({\boldsymbol{P}}),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{P}}$ are the "primitive" variables, ${\boldsymbol{U}}$ the "conserved" variables, ${{\boldsymbol{F}}}^{i}$ are the fluxes, and ${\boldsymbol{S}}$ the source terms. In terms of the primitive variables and metric functions, they can be expressed as

$$\begin{eqnarray}&&{\boldsymbol{U}}({\boldsymbol{P}})=\sqrt{-g}{[\rho {u}^{t},{{T}^{t}}_{t}+\rho {u}^{t},{{T}^{t}}_{j}]}^{T},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\boldsymbol{F}}}^{i}({\boldsymbol{P}})=\sqrt{-g}{[\rho {u}^{i},{{T}^{i}}_{t}+\rho {u}^{i},{{T}^{i}}_{j}]}^{T},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\boldsymbol{S}}({\boldsymbol{P}})=\sqrt{-g}{[0,{{T}^{\kappa }}_{\lambda }{{{\rm{\Gamma }}}^{\lambda }}_{t\kappa }-{{ \mathcal F }}_{t},{{T}^{\kappa }}_{\lambda }{{{\rm{\Gamma }}}^{\lambda }}_{j\kappa }-{{ \mathcal F }}_{j}]}^{T},\end{eqnarray} \tag{ 4 }$$

where g is the determinant of the metric, ${{{\rm{\Gamma }}}^{\lambda }}_{\alpha \beta }$ are the Christoffel symbols, and ${u}^{\alpha }$ are the components of the fluid's four-velocity. The stress-energy tensor can be written as

$$\begin{eqnarray}&&{T}_{\alpha \beta }=\rho {{hu}}_{\alpha }{u}_{\beta }+{{pg}}_{\alpha \beta },\end{eqnarray} \tag{ 5 }$$

where $h=1+\epsilon +p/\rho$ is the specific enthalpy, is the specific internal energy, p is the gas pressure, and ρ is the rest-mass density.

The gas's thermodynamics are governed by an adiabatic equation of state with index ${\rm{\Gamma }}=5/3$ and local cooling. The cooling is introduced via an explicit source term in the stress-energy conservation equation: ${{\rm{\nabla }}}_{\lambda }{{T}^{\lambda }}_{\beta }=-{{ \mathcal L }}_{c}{u}_{\beta }.$ The fluid rest-frame cooling rate per unit volume ${{ \mathcal L }}_{c}$ is determined via the prescription of Noble et al. (2012) in which the gas is cooled at a rate of

$$\begin{eqnarray}&&{{ \mathcal L }}_{c}=\displaystyle \frac{\rho \epsilon }{{t}_{\mathrm{cool}}}\left(\displaystyle \frac{{\rm{\Delta }}S}{{S}_{0}}+\left|\displaystyle \frac{{\rm{\Delta }}S}{{S}_{0}}\right|\right),\end{eqnarray} \tag{ 6 }$$

where ${\rm{\Delta }}S\equiv S\mbox{--}{S}_{0}$ and ${t}_{\mathrm{cool}}$ is the cooling timescale. This prescription cools the gas to the initial entropy S₀ in regions of increased entropy on a timescale ${t}_{\mathrm{cool}}$ usually set to the orbital period of the fluid element. The cooling time is set differently in each of the four distinct regions; one for the circumbinary region, one for each mini-disk, and one for the cavity between the mini-disk and circumbinary regions. We define the circumbinary region as the full azimuthal extent of an annulus extending in radius from $r=1.5a$ out to the end of our numerical domain, where r is the PN Harmonic (PNH) radial coordinate (Blanchet 2014) with its origin at the binary center-of-mass and a is the binary separation. In this region the cooling timescale is the local Keplerian orbital period about the total mass, ${t}_{\mathrm{cool}}=2\pi {(r+M)}^{3/2}/\sqrt{M}$, where M is the total BH mass. The mini-disk regions are defined by the areas within which ${r}_{i}\leqslant 0.45a$, where r_i is the PNH radial distance from the ith BH. In these regions, ${t}_{\mathrm{cool}}$ is again the local orbital period, but now it is calculated in terms of the local Boyer–Lindquist (BL) coordinates with respect to the closest BH (see Section 2.4 for how these are calculated), i.e., ${t}_{\mathrm{cool}}=2\pi {r}_{\mathrm{BL}}^{3/2}/\sqrt{{m}_{i}}$. Finally, in the cavity region between the mini-disks and circumbinary region where there are no quasi-stable circular orbits, ${t}_{\mathrm{cool}}$ is set to be the (constant) period of a local Keplerian orbit at $r=1.5a$.

An important part of our prescription is the use of an approximate, analytic BBH spacetime (Mundim et al. 2014). Because it is analytic, numerically integrating the Einstein Field Equations forward in time is unnecessary; avoiding that computational load allows us to follow inspiraling SMBBHs over the timescales on which gas accumulates (hundreds of binary orbits).

The spacetime is broken into four regions: (1) Inner Zone for BH1 (IZ1), (2) Inner Zone for BH2 (IZ2), (3) Near Zone (NZ), and (4) Far Zone (FZ). For a schematic representation of the spacetime see Figure 1.

Zoom In Zoom Out Reset image size

The IZ is defined as the region where the distance to an individual BH is much less than the binary separation. In these regions, the metric is approximately Schwarzschild, but boosted into the binary center-of-mass frame and perturbed by the gravity of the other BH. The velocity of the boost is set to be the instantaneous velocity of the BH, as found by the PN equations of motion under the quasi-circular approximation (Blanchet 2014). The equations of motion are accurate to 3.5 PN order and include gravitational radiation losses, making the binary's inspiral rate consistent with Einstein's equations for the separations explored in this paper. We note that the spacetime construction is consistent in that the same PN equations of motion are used in the NZ spacetime calculation. The metric in the IZ ensures that our BHs have true horizons by being constructed via perturbation theory around the Schwarzschild BH solution; it is written in horizon-penetrating Cook–Scheel (CS) harmonic coordinates to avoid coordinate singularities (Cook & Scheel 1997). Farther from the BHs, the metric is described by the NZ metric. This metric encompasses the domain where the separation from an individual BH is much greater than the BH's mass but still much less than the gravitational wavelength, $2\pi c/{{\rm{\Omega }}}_{\mathrm{bin}}$ in the general case, $\pi c/{{\rm{\Omega }}}_{\mathrm{bin}}$ when q = 1. In this domain, the metric can be described by PN theory (Blanchet 2014) in PNH coordinates. The mini-disks are contained entirely within the NZ and IZ domains. No problems are created at the zone boundaries because each zone defining the boundary, by construction, shares a region of common validity with its partners (labeled BZ in the figure). In these regions, the metrics are stitched together through asymptotic matching. The resultant global spacetime has been shown to satisfy the Einstein field equations to the expected level of accuracy (Mundim et al. 2014). Finally, the BH trajectories are updated according to equations of motion accurate to 3.5PN-order; in our largest separation run (a = 100M), the separation barely changes over the course of the simulation, while in our smallest separation case (a = 20M), the orbit shrinks by almost 20% in only 14 orbital periods (see Figure 2).

Zoom In Zoom Out Reset image size

The solution for an individual mini-disk in hydrostatic equilibrium in a binary spacetime is not known. However, solutions for hydrostationary torii orbiting a central mass are known (Chakrabarti 1985; De Villiers & Hawley 2003). These solutions suppose that the metric has killing symmetries, ${({\partial }_{t})}^{\alpha }$ and ${({\partial }_{\phi })}^{\alpha }$, and that the metric can be expressed in a spherical coordinate system whose only non-zero diagonal metric component is ${g}_{t\phi }$. Although this is only a crude approximation to the actual spacetime, and tidal forces do create significant departures from stationary behavior, it is a feasible means of constructing an initial state. Following Noble et al. (2012), we construct the solution for an initially isentropic accretion disk in hydrostatic equilibrium around an individual BH as a function of local BL coordinates $({r}_{\mathrm{BL}},{\theta }_{\mathrm{BL}})$, neglecting the presence of the binary companion. We specify the radial location of the inner edge of the disk (${r}_{\mathrm{in}}$) and pressure maximum (${r}_{p}$) under the constraint that $(H/r)=0.1$ at the pressure maximum. For one run ($20M$), ${r}_{\mathrm{in}}$ is inside the innermost stable circular orbit (ISCO); in this case, we set the initial density and pressure to the floor value for $r\lt {r}_{\mathrm{ISCO}}$. Our 2D simulation then evolves the data in the orbital plane, ignoring the material above and below. As we show later, the lack of exact hydrostatic balance leads to a transient, but it decays in one to three orbits.

We relate the local BL coordinates in which the mini-disk initial data are constructed to the PNH coordinates describing the NZ through an intermediate CS harmonic coordinate (Cook & Scheel 1997). For any location specified in Cartesian PNH coordinates, we calculate the corresponding CS coordinates (X_CS) and Jacobian $\tfrac{\partial {X}_{\mathrm{PNH}}}{\partial {X}_{\mathrm{CS}}}$ as in Appendix B of Mundim et al. (2014). Going the other way, we can express the CS coordinates and Jacobian $\tfrac{\partial {X}_{\mathrm{CS}}}{\partial {X}_{\mathrm{BL}}}$ in terms of the local BL coordinate system (Gallouin et al. 2012) through

$$\begin{eqnarray}T & = & {t}_{\mathrm{BL}}+\displaystyle \frac{{r}_{+}^{2}+{\chi }^{2}}{{r}_{+}-{r}_{-}}\mathrm{ln}\left|\displaystyle \frac{{r}_{\mathrm{BL}}-{r}_{+}}{{r}_{\mathrm{BL}}-{r}_{-}}\right|\\ X+{iY} & = & ({r}_{\mathrm{BL}}-m+i\chi ){e}^{i{\phi }_{\mathrm{IK}}}\sin {\theta }_{\mathrm{BL}}\\ Z & = & ({r}_{\mathrm{BL}}-m)\cos {\theta }_{\mathrm{BL}}\\ {\phi }_{\mathrm{IK}} & = & {\phi }_{\mathrm{BL}}+\displaystyle \frac{\chi }{{r}_{+}-{r}_{-}}\mathrm{ln}\left|\displaystyle \frac{{r}_{\mathrm{BL}}-{r}_{+}}{{r}_{\mathrm{BL}}-{r}_{-}}\right|\end{eqnarray} \tag{ 7 }$$

where χ is the dimensional spin parameter (zero for this paper), ${r}_{\pm }=m\pm \sqrt{{m}^{2}-{\chi }^{2}}$, m is the mass of the individual BH, and the initial BH orbital phase (${\phi }_{\mathrm{IK}}$) is assumed to be zero. Combining these two coordinate transformations completes the rule for transforming quantities between the BL and PNH systems.

The actual simulation is done in a different coordinate system we call "warped coordinates" derived from a spherical coordinate system (see Section 2.5). The last stage in initial condition preparation is therefore to transform the data, originally prepared in BL coordinates, from Cartesian PNH coordinates (as described above) to spherical PNH coordinates, and finally to the warped coordinates, i.e.,

$$\begin{eqnarray}&&{X}_{\mathrm{BL}}\to {X}_{\mathrm{CS}}\to {X}_{\mathrm{PNHC}}\to {X}_{\mathrm{PNHS}}\to {X}_{\mathrm{WARP}}\end{eqnarray} \tag{ 8 }$$

where X_PNHC, X_PNHS, and X_WARP are the Cartesian, spherical, and warped representation of PNH coordinates, respectively. Finally, the solution is bi-linearly interpolated onto the numerical grid.

We performed seven different hydrodynamic evolutions, differing primarily by initial binary separation. We tabulate their parameters in Table 1. For the "small," 30M, and 20M runs, the mini-disks' initial outer radii were approximately the Newtonian truncation radius $(0.3a);$ for the "large" and 40M runs, the initial radii were larger $(0.4a)$. This distinction permitted us to explore whether the quasi-steady state obtained is independent of the initial mini-disk size. The binary separations range from 100M (found to be quasi-Newtonian by Zilhão et al. 2015) to 20M, at which relativistic effects are substantial. In Figure 3, we plot the initial density contours around BH1 (the BH initially on the positive x axis) for the 20M binary separation and "small" 50M and 100M binary separation runs. For reference, we mark the black hole horizon, ISCO, the inner edge of the IZ-NZ BZ (the outer edge of the IZ-NZ BZ lies well outside the mini-disk), and the Newtonian estimate for the tidal truncation radius. The majority of the disk mass is located within the IZ-NZ BZ, well outside the ISCO. For the 50M and 100M runs, we put in place a circumbinary disk following the prescription of Noble et al. (2012) with ${r}_{\mathrm{in}}=3{a}_{0}$ and ${r}_{p}=5{a}_{0}$, with only floor values of density and pressure between the circumbinary disk and the mini-disks. We did not observe any substantial inflow from the circumbinary into the mini-disks in these runs because, in the absence of MHD accretion stresses, there is no mechanism to drive accretion into the central cavity aside from numerical diffusion, and that is kept small by our spherical grid. Having seen this behavior in the larger separation runs, in the 20M, 30M, and 40M runs, we set the initial gas density and pressure everywhere outside the mini-disks and in the cavity to floor values, eliminating any circumbinary disk. While astrophysical mini-disks will be influenced by accretion streams from the circumbinary disk, removing such streams allows us to focus on purely gravitational effects of GR on the mini-disk structure.

Zoom In Zoom Out Reset image size

Table 1.  Initial Data Parameters

Run Name	Initial Separation	${r}_{\mathrm{in}}$	${r}_{p}$
100M Large	100M	10M	$18.5M$
100M Small	100M	10M	15M
50M Large	50M	6M	11M
50M Small	50M	6M	9M
40M	40M	5M	9M
30M	30M	3M	$5.5M$
20M	20M	$2.5M$	$4.2M$

Note. Initial data parameters used for the hydrodynamic runs. Radial coordinates are in the local BL system centered on the individual BH. M is the total mass of the binary in all entries.

Download table as:  ASCIITypeset image

Our hydrodynamic simulations are performed in the equatorial plane of a dynamic, double fish-eye (warped) spherical coordinate system whose origin is at the binary center-of-mass (Zilhão & Noble 2014). Cells within this coordinate system are spaced uniformly in numerical spatial coordinates $\{{x}^{i}\}$. By this means, we are able to focus resolution in the vicinity of the BHs and rarefy resolution in the dynamically less interesting portions of the cavity. Near the BHs, the grid is approximately Cartesian, while, farther out, the grid is spherical.

Given physical PNHS coordinates $(r,\phi )$, Zilhão & Noble (2014) considered three regions of warping: one for each BH and one for the region between the BHs. The coordinate transformations from PNHS to WARP take a different form in each region and are smoothly interpolated (analytically) to match each other at the regions' boundaries. The transformations are designed so that ${\rm{\Delta }}r(r)$ has a local minimum at $r=a/2$ and ${\rm{\Delta }}\phi (\phi )$ has a local minimum at the instantaneous azimuthal positions of the BHs. Approximately 32 cells span each black hole horizon in each dimension, a resolution chosen to ensure that the near-horizon spacetime is well resolved. Parameters control various aspects of the transformations. In particular, it is important to achieve a smooth transition from nearly Cartesian cells near the BH horizons to nearly spherical cells far from the BHs. The values we used are stated in Table 2. Quantitative expressions for our transformations and the the parameter definitions may be found in Equations (29)–(32) of Zilhão & Noble (2014). For convenience, we also include in Table 2 the number of cells within the Newtonian tidal truncation radius for each binary separation.

Table 2.  Warped Grid Parameters

Initial Separation	100M	50M	40M	30M	20M
${\delta }_{x1}$	0.1	0.1	0.1	0.1	0.1
${\delta }_{x2}$	0.1	0.1	0.1	0.1	0.1
${\delta }_{x3}$	0.25	0.18	0.21	0.21	0.18
${\delta }_{x4}$	0.25	0.18	0.21	0.21	0.18
${\delta }_{y3}$	0.25	0.18	0.21	0.21	0.1
${\delta }_{y4}$	0.25	0.18	0.21	0.21	0.1
${\delta }_{z}$	0.4	0.4	0.4	0.4	0.4
ax1	1.7	1.67	1.64	1.60	1.50
ax2	1.7	1.67	1.64	1.60	1.50
az	0	0	0	0	0
hx1	20	20	20	20	20
hx2	20	20	20	20	20
hx3	20	20	20	20	20
hx4	20	20	20	20	20
hy3	20	20	20	20	10
hy4	20	20	20	20	10
hz	20	20	20	20	20
s1	4	4	0.01	0.01	0.01
s2	4	4	0.01	0.01	0.01
s3	4	4	0.01	0.01	0.01
b1	10	10	10	6.5	6.5
b2	10	10	10	6.5	6.5
b3	40	40	40	40	6.5
${R}_{\mathrm{out}}$	500M	750M	120M	90M	60M
Cell Count	400 × 400	400 × 400	300 × 320	300 × 320	600 × 640
Cells per Mini-Disk	25942	19258	14782	14798	48532

Note. Parameters of the warped grid used for each initial binary separation. Please see Equations (29)–(32) of Zilhão & Noble (2014) for the explicit expressions defining the warped system and the significance of these parameters. For all runs, the inner radial cutout is set to 1M.

Download table as:  ASCIITypeset image

The warped coordinates are derived from spherical coordinates, and therefore possess a coordinate singularity at the origin. For this reason, we excise a sphere of radius M from the computational domain. In Figure 4, we plot the grid for the inner $1.5a$ of the 20M, 50M, and 100M runs. Throughout the length of our simulations, $\lesssim 1 \%$ of the initial gas is lost through the origin cutout.

Zoom In Zoom Out Reset image size

We consistently impose outflow boundary conditions at the radial x¹ boundaries, requiring u^r to be oriented out of the domain there. If at any time the velocity points inward, the value is set to zero, and we recalculate the primitive velocity. To enforce our restriction to 2D, we specify all x² ghost zones as copies of the physical cell in the equatorial plane. This effectively amounts to considering all primitive variables as vertically averaged quantities that are constant w.r.t. θ. Finally, we apply periodic boundary conditions in the azimuthal x³ direction and cover all $2\pi$.

To differentiate between hydrodynamical and purely gravitational effects, we also simulated ensembles of particles orbiting a single BH in the BBH system. We performed these calculations in Bothros (Noble et al. 2007, 2009, 2011), which solves the geodesic equations of motion

$$\begin{eqnarray}&&{\partial }_{\lambda }{x}^{\alpha }={N}^{\alpha }\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\partial }_{\lambda }{N}_{\alpha }={{{\rm{\Gamma }}}^{\kappa }}_{\alpha \eta }{N}_{\kappa }{N}^{\eta }.\end{eqnarray} \tag{ 10 }$$

Here ${N}^{\alpha }$ is the particle's four-velocity, and λ is the affine parameter, or in the case of time-like paths, the proper time for the particle. We terminate the geodesic if it exits the central cavity or falls within an ISCO. Bothros was originally written for stationary spacetimes, so support for handling dynamic spacetimes was added for this investigation. Spatial and temporal derivatives needed for evaluation of the Christoffel symbols are computed by centered fourth-order finite differences, using a resolution of ${10}^{-6}M$ in space and time for all runs. The Christoffel symbols are calculated for every sub-step of the multiple sub-step Bulirsch–Stoer procedure (Press et al. 1992) used to integrate the geodesics in time.

For a given particle orbiting the non-spinning BH, we specify the initial four-velocity in local BL coordinates as

$$\begin{eqnarray}&&{u}^{\alpha }=\gamma (1,0,0,\eta {{\rm{\Omega }}}_{K})\end{eqnarray} \tag{ 11 }$$

where ${{\rm{\Omega }}}_{K}$ is the Keplerian orbital frequency and η is a parameter selected with uniform probability density within $[0.9,1.1]$ in order to sample local tidal effects. This four-velocity is then transformed to the X_PNHC coordinate system and evolved in Bothros using the spacetime of Section 2.3. We perform three runs of 5376 test particles each at 100M, 50M, and 20M binary separations. The geodesics are launched from 16 × 16 uniformly spaced points in $({r}_{\mathrm{BL}},{\phi }_{\mathrm{BL}})$ space. We use different ranges of r_BL depending on the initial separation: ${r}_{\mathrm{BL}}\in [15,30]M$ for 100M, ${r}_{\mathrm{BL}}\in [7.5,15]M$ for 50M, and ${r}_{\mathrm{BL}}\in [3.00,6.75]M$ for 20M binary separations. The full ${\phi }_{\mathrm{BL}}\in [0,2\pi ]$ extent is used for all runs. At each location, 21 geodesics are each launched with a different value of η.

Our simulations span a range of binary separations from $a\approx 100M$, where one might expect physics to be quasi-Newtonian, to $a\approx 20M$, where relativistic effects, including binary inspiral, become very important. Their general character during the first few binary orbits is illustrated by the snapshots of the mini disk around BH1 shown in Figure 5. A casual look suggests there is little variation, either as a function of time for fixed separation or as a function of separation; a closer look reveals both significant variability and important trends with separation.

Zoom In Zoom Out Reset image size

The first snapshot, at 0.5 orbits, shows the disks at their greatest extent, as the disks undergo an expansion due to the departures from hydrostatic balance in our initial conditions caused by the omission of tidal forces in our approximate hydrostatic balance equations. However, one or two more orbits suffice for this transient to decay, letting the mini disks achieve their approximate long-term structure. Once equilibration has completed, we find that both the hydrodynamic and test particle mini-disks have settled down to similar stable configurations.

Similarly, although all three separations exhibited similar initial transients, they also show a dependence on separation. In particular, note the greater gas density on the side toward BH2 (the left side) at smaller separations. Although not visible in these snapshots, gas is readily shared back and forth between the two mini disks in a bar-like region centered on the L1 point. The fraction of all available gas finding itself in this region grows sharply with decreasing separation, rising more than an order of magnitude by a = 20M (see Section 3.4 for further details).

Finally, another effect better seen in other figures (see below) is the development of spiral shocks within the mini-disks. Qualitatively similar spiral features have been observed in simulations of accretion disks in cataclysmic variables (Ju et al. 2016) and BBHs (Ryan & MacFadyen 2016). However, we find that relativistic effects drive the formation of an m = 1 mode in addition to the m = 2 mode predicted by Newtonian gravity (see Section 3.5 for further details).

As discussed in the Introduction, our principal goal is to explore how mini-disk dynamics may change as relativistic effects become more important. The simplest way to highlight how they enter is to study the lowest-order PN corrections to the metric. At this level of approximation, one can isolate the gravitational potential Φ in a frame corotating with the binary through the relation

$$\begin{eqnarray}&&{g}_{{tt}}=-(1+2{\rm{\Phi }}).\end{eqnarray} \tag{ 12 }$$

This potential combines what in Newtonian language would be called the genuine gravitational potential with the centrifugal contribution of the corotating frame. In Figure 6, we plot Φ, scaled by the binary separation, along the line connecting the BHs for all separations we simulated. In Newtonian gravity, all the rescaled binary potentials would be identical. The fact that they differ, and in a fashion that is monotonic with binary separation, illustrates the relativistic modifications to Newtonian expectations. The net effect of the PN corrections is to create shallower potential wells at closer binary separations, particularly in the secondary's Roche lobe. There, this effect is, in relative terms, stronger for mass ratios farther from unity (see Figure 7). For example, for a separation of 20M, when q = 1, the PN potential difference between the L1 point and the nearest edge of the mini-disk is $0.75\times$ the Newtonian potential difference; the corresponding ratios (for the secondary) when q = 0.33 and q = 0.1 are 0.53 and 0.38, respectively. In the primary's Roche lobe, the trend with mass ratio is rather weaker, with the ratio of PN to Newtonian potential difference ≃0.75–1 and varying little with q (see Figure 8). These PN changes in the potential will be a recurring theme in our description of the simulations' behavior.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of our primary goals is to see whether relativistic effects alter the structure of mini-disks. To accomplish this, the first step is to demonstrate that the structures examined are not influenced by transients due to our initial conditions. We do so in two ways. The first is by contrasting simulations in which the initial conditions differ, and confirming that after the decay of transients they reach similar states. The second makes use of time-averaging over a period later than the $\simeq 2{t}_{\mathrm{bin}}$ required for transient decay.

3.3.1. Contrasting Initial Conditions

To test our results' sensitivity to initial conditions, we ran two versions of each of the larger separation cases (100M and 50M). In the "small" runs, the disks were initially filled with matter to a point just inside our analytic estimate of the tidal truncation radius, ${r}_{t}\simeq 0.3a$, an estimate appropriate to circular-orbit equal-mass Newtonian binaries (Paczynski 1977; Papaloizou & Pringle 1977; Artymowicz & Lubow 1994). In the "large" runs, in the initial state the disks extended to ≃(0.4–0.45)a.

Although it hides departures from axisymmetry, we have chosen the distribution of mass enclosed within a given radius as a diagnostic of internal structure. In terms of the conserved rest mass, this quantity is

$$\begin{eqnarray}&&M(\lt {r}_{\mathrm{BL}})={\int }_{{r}_{\min }}^{r}\,{{dr}}^{\prime }d{\phi }_{\mathrm{BL}}d\theta \rho {u}_{\mathrm{BL}}^{t}\sqrt{-{g}_{\mathrm{BL}}},\end{eqnarray} \tag{ 13 }$$

where ${r}_{\min }=2{M}_{1}$. Because the code data are defined with respect to warped PNH coordinates, to obtain the data necessary for the integral, we first transformed the code data into the BL system, resolving four-vector components with respect to the local BL basis, and then interpolated onto a regular grid in BL coordinates.

In Figure 9, we show how well the initial conditions relax to essentially the same structure by contrasting the large and small initial mass-enclosed distributions for 100M and 50M separation with the distribution averaged over the interval $[2{t}_{\mathrm{bin}},3{t}_{\mathrm{bin}}]$. Although the 100M Small case differs somewhat from 100M Large in the inner quarter of mass even in the later state of the disk, the remainder of the two 100M mass-enclosed distributions are very close, as are the entire distributions for the Large and Small 50M cases. Because our prescription for the initial condition omits tidal forces, and tidal forces partially cancel the nearby BH's gravity, the Small cases are over-pressured; as a result, they expand. However, tidal forces also prevent matter from staying near the disk when it lies at distances beyond r_t. For this reason, the Small cases cease expanding when they reach the tidal truncation limit, while in Large cases gas beyond r_t does not stay attached to the disk.

Zoom In Zoom Out Reset image size

On the basis of this success in quickly reaching a state almost entirely independent of initial conditions, we ran only one case for the a = 20M, 30M, and 40M separations. For a = 40M, the initial disk extended to $\simeq 0.4a;$ for the two smaller separations, it was filled to only $\simeq 0.3a$.

3.3.2. Time Averaging

Before continuing on to other measures of disk structure, it is necessary to elaborate on our methods of time-averaging. We chose to begin time-averaging for all runs at $2{t}_{\mathrm{bin}}$ because the results shown in Figure 5 demonstrate that initial transients have almost entirely decayed by this time. We ended it at $3{t}_{\mathrm{bin}}$ for all cases except a = 20M because the 100M simulation stopped at about this time, and we wished to enforce a consistent procedure on all our cases. The smallest separation case demanded special treatment because in it the binary's separation shrinks appreciably over a single binary orbital period. In this case, we therefore average successive intervals of duration ${t}_{\mathrm{bin}}$ beginning at roughly 2, 4, 6, 8, 10, 12, 14, and 15.5 binary orbits. These correspond to the times at which the binary separation is $19.5M$, $19.0M$, $18.5M$, $18.0M$, $17.5M$, $17.0M$, $16.5M$, and $16.0M$.

In Figure 10, we plot the time-averaged density contours both on a linear and a logarithmic scale, overlaid with contours of the binary potential evaluated in the frame corotating with the binary. Several features can be clearly seen in these two representations. In a linear scale it is apparent that, rather than being axisymmetric around its BH, the density on the side of each mini-disk nearer the L1 point is higher than on the opposite side. Interestingly, the orientation of the density gradient is rotated by $\simeq \pi /5$ relative to the line between the two BHs when a = 100M, but this rotation diminishes with smaller binary separation. On a logarithmic scale, we see that as the binary separation shrinks, the relative amount of mass in the zone between the disks rises, hinting that the density cutoff at the edge of the disk facing the center-of-mass is becoming less sharp.

Zoom In Zoom Out Reset image size

For a more quantitative view of the time-averaged structure, we study the surface density,

$$\begin{eqnarray}&&{\rm{\Sigma }}({r}_{\mathrm{BL}},{\phi }_{\mathrm{BL}})=\displaystyle \frac{\int \,d\theta \rho \sqrt{-{g}_{\mathrm{BL}}}}{\sqrt{{g}_{{\phi }_{\mathrm{BL}}{\phi }_{\mathrm{BL}}}(\theta =\pi /2)}},\end{eqnarray} \tag{ 14 }$$

in local BL coordinates. The surface density profile of all the BH1 disks is displayed in Figure 11. To be more precise, this figure shows the surface density averaged within a pair of wedges, each centered on BH1 and stretching 45° on either side of the line between the two BHs. Several qualitative facts are apparent in this figure. First, the location of the near-side surface density peak varies little between a = 100M and a = 30M, but it moves sharply outward when a = 20M. On the other hand, the far side peak position changes little with separation, but the net sense is to move inward as the separation becomes smaller. In addition, although at a = 100M the two surface density peaks, near-side and far-side, are very nearly equal in height, the near-side peak becomes increasingly dominant as a decreases.

Zoom In Zoom Out Reset image size

3.3.3. Quantitative Measures of Tidal Truncation

The time-averaged surface density profile can be used to determine quantitatively the edge of the mini-disks in two different ways. In the first, we define it as the point where the radial gradient of the time- and azimuthally averaged surface density has the greatest magnitude, both on the side nearest ($0.75\pi \leqslant {\phi }_{\mathrm{BL}}\leqslant 1.25\pi$) to and farthest ($-0.25\pi \leqslant {\phi }_{\mathrm{BL}}\,\leqslant 0.25\pi$) from the binary companion (see Figure 12). In the second, we search for an outer limit on bound streamlines of the fluid element's velocity field in snapshots taken during the quasi-steady state (see Figure 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Table 3, we tabulate our estimates for the edge of the mini-disk around BH1 using both methods on both the sides nearest $(\phi =\pi )$ and farthest from the binary companion $(\phi =0)$. Echoing the behavior of the near-side surface density peak, for separations down to $a\gtrsim 30M$, the gradient definition yields near-side truncation radii approximately $0.23a$, while for $a\simeq 20M$, it moves to ≃(0.3–0.35)a. On the far-side, the radii fluctuate between 0.23a and 0.30awith little trend as a function of separation. The streamline truncation radii tell a different story. They vary more irregularly with binary separation, ranging between $0.3a\mbox{--}0.4a$ (near-side) and $0.19a\mbox{--}0.24a$ (far-side). These contrasting results may be due, in part, to a larger measurement error in determining the truncation radius from streamlines: contrasting BH1 and BH2 by this measure at the same time, the radius can differ by as much as $0.05a$. In a fashion qualitatively similar to Newtonian gravity, the disks appear to be significantly larger on the near-side than on the far-side, particularly when defined by streamlines, but also (for $a\lesssim 20M$) when defined by surface density gradients; in other words, the truncation "surface" is asymmetric, resembling a Roche-lobe shape more than a sphere. In addition, particularly as judged by the surface density gradient, the disks extend farther on the near-side as the binary separation shrinks. For larger separations, then, these results are in qualitative agreement with the Newtonian ${r}_{t}\simeq 0.3a$ (Paczynski 1977; Papaloizou & Pringle 1977; Artymowicz & Lubow 1994), but even in that regime it is worth noting the non-circularity of the disks.

Table 3.  Truncation Measurements

$a[M]$	${grad}({\rm{\Sigma }}(r,\phi ))$	Streamlines
100	0.22–0.30	0.30–0.24
50	0.24–0.23	0.33–0.17
40	0.24–0.28	0.34–0.24
30	0.22–0.24	0.40–0.19
19.5	0.31–0.27	0.40–0.21
19.0	0.32–0.26	0.33–0.23
18.5	0.31–0.29	0.30–0.22
18.0	0.31–0.29	0.33–0.22
17.5	0.34–0.28	0.39–0.22
17.0	0.33–0.28	0.35–0.19
16.5	0.35–0.29	0.32–0.21
16.0	0.29–0.27	0.34–0.23

Note. Measured radial extent of the mini-disks on the NEAR $({\phi }_{\mathrm{BL}}=\pi )$-FAR $({\phi }_{\mathrm{BL}}=0)$ sides of the disk as fractions of the binary separation for the surface density and streamline method.

Download table as:  ASCIITypeset image

As previously remarked in Section 3.1, matter sloshes back and forth within a bar-like region centered on the binary center-of-mass. Its mass during early times of our simulations is largely due to a transient whose origin lies in our only approximately hydrostatic initial condition for the mini disks. That this should be so is demonstrated most dramatically by the 50M and 100M binary separation runs, in which the ratio of sloshing mass to total mass quickly drops from its initial peak as gas settles back onto the mini-disks (see Figure 14). However, if we consider the mass in this region only after decay of the transient, we find a remarkable increase as the separation diminishes below $\simeq 50M$. To quantify this statement, we must first define this region more precisely: for our purposes, it is a rectangle in the corotating frame with dimensions of $(\tilde{x}\times \tilde{y})=(0.4a(t)\times 0.6a(t))$, centered on the center-of-mass. Because the separation a in the initially 20M separation run decreases with time, we adjust this box size as a function of time so that its dimensions are always the same fraction of a(t). The motivation for choosing these dimensions can be seen by inspection of Figure 15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Employing this definition, we can calculate the ratio ${M}_{\mathrm{slosh}}/{M}_{\mathrm{cav}}$ as a function of time measured in binary orbital periods for the 20M, 30M, 40M, 50M Large, and 100M Large hydrodynamic runs. Here ${M}_{\mathrm{slosh}}$ is the mass within the box just defined, while ${M}_{\mathrm{cav}}$ is the mass within a circle of radius $1.5a(t)$ from the center-of-mass. This quantity is plotted in Figure 14. All the different separation runs begin with similar values of this ratio, $\simeq (5\mbox{--}10)\times {10}^{-3}$. In the quasi-Newtonian 100M and 50M separation runs, this ratio drops to ≃(2–8) × 10⁻⁴ within three binary orbits. In sharp contrast to this behavior, as the simulations become progressively more relativistic, the mean fractional gas content of the sloshing region increases monotonically. In the 20M separation run, this ratio maintains its initial value—$\simeq 5\times {10}^{-3}$—albeit with a factor of several fluctuations, for the entire 14-orbit duration of the simulation. In other words, a decrease in separation by about a factor of two leads to an increase in ${M}_{\mathrm{slosh}}/{M}_{\mathrm{cav}}$ of an order of magnitude. The most natural explanation for this dramatic change is also the progressively shallower gravitational potential found as the system becomes increasingly relativistic (see Figure 6). As a result, gas can more easily flow out of the potential well of a single BH into the center-of-mass region. However, it cannot remain there in a stationary state because the L1 region is dynamically unstable. Once in the sloshing region, gas can only move back and forth between the relatively stable regions near the mini-disks. Upon reaching the truncation radius of a disk, the streams shock.

To gain further insight into the properties of the sloshing region, in Figure 15, we plot density contours at four selected times. This figure demonstrates that the back-and-forth motions are concentrated into discrete streams, but sometimes there are two streams at once, while at other times there is only one. Due to the q = 1 symmetry of our system, the internal motions of the dual stream features are equal in magnitude and opposite one another. We speculate that the number of streams is related to interaction with the spiral structure (see Section 3.5), which exhibits both m = 1 and m = 2 components.

The stream pairs carrying mass between the mini-disks appear to form in a quasi-periodic manner. To characterize this behavior quantitatively, we first pre-whiten the a = 20M data by removing any secular linear trends in ${\tilde{M}}_{\mathrm{slosh}}(t)\,\equiv {M}_{\mathrm{slosh}}(t)/{M}_{\mathrm{cav}}(t)$. We calculate this pre-whitened function (ψ) as

$$\begin{eqnarray}&&\psi (t)=\displaystyle \frac{{\tilde{M}}_{\mathrm{slosh}}(t)-{\tilde{M}}_{\mathrm{fit}}(t)}{| {M}_{0}| }\end{eqnarray} \tag{ 15 }$$

where M₀ is the largest value of ${\tilde{M}}_{\mathrm{slosh}}$, ${\tilde{M}}_{\mathrm{fit}}(t)$ is a linear fit for ${\tilde{M}}_{\mathrm{slosh}}(t)$, and we use only data for $2000\leqslant t/M\leqslant 6000$. This corresponds roughly to 3.5–11 binary orbital periods.

We then compute Ψ, the Fourier transform of ψ. In Figure 16, we plot $| {\rm{\Psi }}{| }^{2}$, the Fourier power density, as a function of frequency. It possesses two distinct peaks, with angular frequencies roughly 2 and 2.75 times the mean orbital frequency of the binary $({\bar{\omega }}_{\mathrm{bin}})$. The width of these peaks is comparable to the range of binary orbital frequencies during the period included in the Fourier analysis, from $\simeq 1.05\times$ the initial frequency to $\simeq 1.23\times$ that frequency.

Zoom In Zoom Out Reset image size

Within the main body of the mini-disks, we observe the formation of spiral density waves. They can be seen in density plots (Figures 5, 10, and 15) and stand out even more clearly in a plot showing the local sound speed (Figure 17) because a strong shock raises the local temperature by a much larger factor than it increases the density. These spiral density waves are of special interest, in part, because they have the potential to alter the locally emitted spectrum, but also because they may contribute significantly to accretion stresses. In fact, spiral density waves were put forward very early on as a candidate mechanism to explain all accretion stresses (Lynden-Bell & Pringle 1974), though subsequent studies determined that they were unlikely to be a general answer to the problem (see (Shu & Lubow 1981; Papaloizou & Lin 1995; Boffin 2001 and references therein). In regions where the spiral density wave pattern speed is less than the local fluid orbital frequency of the disk, the spiral density wave orbits backward through the disk, thereby carrying negative angular momentum with respect to the fluid in the disk. Once the spiral density wave amplitude becomes nonlinear, it can steepen into shocks, allowing it to couple to the disk fluid via dissipation (Papaloizou & Lin 1995; Goodman & Rafikov 2001; Heinemann & Papaloizou 2012; Rafikov 2016). The resulting loss in the matter's angular momentum introduces a radial component to the flow and transports angular momentum outward in the disk, thereby facilitating accretion onto the BH.

Zoom In Zoom Out Reset image size

In accretion disks in binary systems, these spiral density waves are the result of the perturbing non-axisymmetric gravitational potential of the binary companion (see, e.g., Savonije et al. 1994 and references therein). The density waves in the disk due to the perturbing gravitational potential are launched at Lindblad resonances (Goldreich & Tremaine 1979; Papaloizou & Lin 1984). In Newtonian gravity, the perturbing potential (in a coordinate system centered on the central mass surrounded by the disk) due to the binary companion on a circular orbit is given by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{pert}}=-\displaystyle \frac{{M}_{2}}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{2}| }+\displaystyle \frac{{M}_{2}({\boldsymbol{r}}\cdot {{\boldsymbol{r}}}_{2})}{{r}_{2}^{3}},\end{eqnarray} \tag{ 16 }$$

where M₂ is the mass of the perturber, ${\boldsymbol{r}}$ is the position vector, and ${{\boldsymbol{r}}}_{2}$ is the position vector of the secondary (see, e.g., (Binney & Tremaine 1987; Savonije et al. 1994)). The first term in Equation (16) is the Newtonian potential of the secondary; the second term accounts for centrifugal forces due to the fact that the corotating frame is non-inertial.

To analyze the effects of the perturbing potential, it is useful to expand it as a Fourier series in the azimuthal angle. The m-th component of the perturbing potential excites a spiral density wave with m arms. If the perturber's orbit is not too close to the outer edge of the disk, the dominant azimuthal Fourier component of the Newtonian perturbing potential is the m = 2 mode (Savonije et al. 1994). Consequently, in Newtonian simulations of these systems, a two-armed spiral density wave is often found to develop in the disk (see, e.g., the simulations of Savonije et al. 1994; Makita et al. 2000; Ju et al. 2016 or the early evolution of the disks studied in Kley et al. 2008).

Because the tidal perturbation remains constant in the frame corotating with the binary, the induced spiral density waves' pattern speed is the binary orbital frequency, which is slower than the orbital velocity of the disk fluid. In the fluid frame, a mode of order m creates a spiral wave whose phase speed is m times the binary orbital frequency, as the fluid encounters m spiral arms in each orbit.

Although Newtonian simulations find the development of clearly visible two-armed spiral density waves as the result of the dominant m = 2 mode in the perturbing potential, our simulations show a more complicated flow morphology in the mini-disk for all binary separations studied. To study the structure of the spiral waves in our relativistic simulations, we compute the amplitude of different m modes in the disk surface density as a function of time (see, e.g., Zurek & Benz 1986; Heemskerk et al. 1992):

$$\begin{eqnarray}&&{D}_{m}(t)={\int }_{{r}_{\min }}^{{r}_{\max }}\,{{dr}}_{\mathrm{BL}}d{\phi }_{\mathrm{BL}}{\rm{\Sigma }}({r}_{\mathrm{BL}},{\phi }_{\mathrm{BL}};t){e}^{-{im}{\phi }_{\mathrm{BL}}}\sqrt{-{g}_{\mathrm{BL}}},\end{eqnarray} \tag{ 17 }$$

where ${r}_{\min }=0.13a$ and ${r}_{\max }=0.38a$. The results are insensitive to the selection of ${r}_{\min }.$ As a check, we performed this azimuthal mode analysis in two frames, a frame co-moving with BH1, and a frame rotating at the binary orbital period. Both calculations give the same mode strengths, as expected.

The time-dependence of these mode amplitudes for initial binary separations of 20M, 50M, and 100M are shown in Figure 18. The m = 1 amplitude is larger than the m = 2 amplitude in all cases, though ${D}_{1}/{D}_{2}$ appears to diminish over time for the larger separations. At late times in the a = 50M run, ${D}_{1}/{D}_{2}\simeq 3$. When the separation is as small as 20M, ${D}_{1}/{D}_{2}$ oscillates around a constant mean value of $\simeq 4$. Although the 100M case did not run long enough to reach "late times," ${D}_{1}/{D}_{2}$ appears to be roughly ∼2, suggesting a slow increase in this ratio with decreasing separation.

Zoom In Zoom Out Reset image size

The origin of this departure from Newtonian experience appears to lie in the different structure of the tidal potential in the PN regime. To demonstrate this, we computed the amplitudes with respect to m in the potential in much the same fashion as we did for the surface density; we call them S_m. In order to calculate them, we used the same PN approximation for the perturbing potential as we did when discussing the tidal potential itself, i.e., ${\rm{\Phi }}=-(1+{g}_{{tt}})/2$. However, before computing the mode integrals, we subtracted off the isotropic Newtonian potential due to the individual BH. The Newtonian ratio ${S}_{1}/{S}_{2}=0.143$, irrespective of binary separation, because Newtonian gravity is scale-free. However, relativistic gravity is scale-dependent (see Table 4). Even at a = 100M, ${S}_{1}/{S}_{2}$ is close to triple the Newtonian value, and it grows monotonically with decreasing binary separation, exceeding unity for separations slightly less than 30M. Relativistic contributions to m = 1 modes can be significant at separations as large as 100M because they first appear in the PN expansion in the lowest-order non-Newtonian terms (see, e.g., the explicit expression for the PN Lagrangian governing a test-particle in a binary system in Ratkovic et al. 2005). In addition, the same low-order terms also create new contributions to m = 2 modulation, altering the Newtonian component in both magnitude and phase.

Table 4.  Azimuthal Potential Modes

$a[M]$	${S}_{1}/{S}_{2}$
Newtonian	0.143
100	0.412
50	0.665
40	0.825
30	0.990
19.5	1.209
19.0	1.212
18.5	1.239
18.0	1.276
17.5	1.298
17.0	1.335
16.5	1.359
16.0	1.382

Note. Strength of the m = 1/m = 2 modes of the perturbing potential as a function of separation.

Download table as:  ASCIITypeset image

As a final point, we note that we have performed a test simulation of a BBH system in which initially only BH1 was surrounded by a mini-disk, in order to test if shocks induced by the sloshing might produce spiral density waves in the mini-disks. The morphology of the fluid flow and mode strengths were very similar to those from simulations with two mini-disks. We therefore conclude that the main driving factor of the spiral waves in our mini-disks is indeed the tidal field of the binary companion.

Previous discussions of the shapes of disks in binary systems have been almost exclusively posed in terms of Newtonian gravity. In addition, when describing their size, it has generally been customary to specify only a mean radius from the mass at the center of such a disk, often without a precise operational definition of what exactly makes that radius the "size" of the disk. In the course of our investigation of PN effects on the tidal truncation of disks in binaries, we have found that quantitative description demands such a definition: these disks are, at the $\sim 20 \%$ level, non-circular even in the Newtonian limit, and their sizes (when more precisely defined) do begin to change as gravity departs from Newtonian.

In principle, there could be many definitions of the "size" of a disk in a binary system. Here we explored the implications of two choices. One, a definition focusing on the surface mass density profile, is suited to answering questions such as "if an external fluid stream strikes the disk, where does it encounter the full inertia of the disk?" The other, a definition focusing on the shapes of disk material's orbits, is designed to discriminate between fluid elements that repeatedly orbit the central mass of a disk from those that may pass near the disk, but then swing away from it.

Even in the Newtonian limit (exemplified by our 100M separation simulation), the disks appear to be asymmetric in terms of both these measures. The far side (away from the other mass in the binary system) extent defined by the surface density gradient agrees with the near side extent defined by streamlines, and both are consistent with the usual Newtonian estimate for equal-mass binaries ($\simeq 0.30a$), but the near-side surface density definition and the far side streamline definition both give extents that are $\sim 20 \%$ smaller. In other words, by the surface density definition, disks are cut off more tightly on the near side than the far, while by the stream line definition they extend farther from the central mass on the near side. Thus, at this level of precision, it is clear that a single "truncation radius" fails to be an adequate description: the disks are neither round nor admit a single definition appropriate to all intuitive meanings of the "disk edge."

Relativistic effects also begin to alter mini-disk shapes once a is a few tens of M or less, particularly with respect to the near-side surface density definition. On the near-side, the disk extends to a progressively larger and larger fraction of a as the binary separation moves farther and farther into the relativistic regime. When $a\simeq 20M$, the (surface density) near edge is $\simeq 35 \%$ farther from the central mass than in the Newtonian limit. As we have already remarked, this effect can be attributed directly to the progressively shallower potential gradient in the L1 region produced by increasingly relativistic gravity (as illustrated in Figure 6). When the binary mass ratio is less than unity, the relativistic effects are, if anything, stronger (Figure 7). A similar result, also in the PN formulation but posed in terms of the volume occupied by the secondary's Roche lobe, was found by Ratkovic et al. (2005). It is worth noting, however, that for their parameters (separations of 5M and 10M), the Roche lobe volume decreases with greater relativistic effects when $q\lt 0.7$.

On the basis of our results, we might speculate that at still smaller separations the mini-disks continue to exist until the truncation radius (defined in any of these senses) becomes only somewhat larger than the absolute minimum scale set by the ISCO. However, once the inspiral becomes more rapid than the internal accretion rate within such a disk, matter from the outer edge, which now extends beyond the truncation radius, must lose its binding to an individual BH, and instead travel through both its original Roche lobe and the partner's. Such effects are closely related to the topic of the next subsection, the "sloshing region."

Gas extending across the L1 region in a system with binary disks has been previously observed in stellar systems (Mayama et al. 2010); it has also appeared in Newtonian hydrodynamics simulations of binary stellar systems (de Val-Borro et al. 2011; Nelson & Marzari 2016) and SMBBHs (Farris et al. 2014; D'Orazio et al. 2016), although without comment in the latter work. However, we have found a number of new properties of this gas. It is in a constant state of back-and-forth motion, and, as shown in Figure 16, the variations in this motion are strongly modulated at two frequencies, one very close to $2{\bar{\omega }}_{\mathrm{bin}}$, the other at $2.75{\bar{\omega }}_{\mathrm{bin}}$. Second, and most surprisingly, the fraction of disk mass in this sloshing region rises steeply as relativistic effects become important.

The quantity of gas in the sloshing region is determined through a complex mechanism. When sloshing gas shocks against the edge of a mini-disk, it "spalls" off gas parcels, giving them enough energy to travel across the L1 region. The amount required is comparatively small. Although Figures 6–8 show the effective gravitational potential in the rotating frame of the binary, they do not include the contribution to the effective potential associated with the orbital motion of a fluid element around one of the BHs. Even in the Newtonian limit, this contribution causes the total effective potential to be rather shallower than the curves shown in this figure; when PN effects further decrease the change in effective potential from the edge of the mini-disk to the L1 point, the energy delivery threshold for liberation of gas becomes only a small fraction of GM/a. In agreement with this potential picture, the mean fraction of available gas in the sloshing region (${M}_{\mathrm{slosh}}/{M}_{\mathrm{cavity}}$) increases monotonically with decreasing binary separation, with a particularly sharp increase for $a\lesssim 30M$. In addition, animations of our simulations show correlations between fluctuations in spiral wave structure and ejection of matter into the sloshing region. It is possible that when accretion streams from a circumbinary disk enter the picture, their impact may also influence gas injection into the sloshing region.

These new properties may lead to a significant new observable. Every time a sloshing stream shocks against the edge of a mini-disk, the associated energy dissipation is available for photon radiation. Moreover, the periodic character of these motions means that if the cooling time post-shock is shorter than the dynamical period, the radiation may be strongly modulated. Because the characteristic speed of this motion is comparable to the orbital speed, and therefore increases $\propto {a}^{-1/2}$, and the amount of mass involved increases much more rapidly with decreasing binary separation when $a\lesssim 30M$, the energy available for this potentially periodic radiation signal should increase sharply as the binary inspirals through the last few tens of M in separation.

More quantitatively, for BBH systems with mass ratios of the order of unity, we may estimate the time-averaged heat release in the sloshing region by

$$\begin{eqnarray}&&{L}_{\mathrm{slosh}}\sim ({M}_{\mathrm{slosh}}/{M}_{\mathrm{cav}})(\dot{M}{t}_{\mathrm{in}}){v}_{\mathrm{orb}}^{2}{\omega }_{\mathrm{slosh}},\end{eqnarray} \tag{ 18 }$$

where the total mass of the mini-disks ${M}_{\mathrm{cav}}$ is determined by the accretion rate $\dot{M}$ and the inflow time through an individual mini disk ${t}_{\mathrm{in}}$. We further estimate the speed of the sloshing to be comparable to the binary orbital speed ${v}_{\mathrm{orb}}$. The repetition rate of the sloshing cycle is ${\omega }_{\mathrm{slosh}}$, which we have already determined to be $\simeq 2{\omega }_{\mathrm{bin}}$. For greater insight, this luminosity estimate may be rewritten as

$$\begin{eqnarray}&&{L}_{\mathrm{slosh}}\sim ({M}_{\mathrm{slosh}}/{M}_{\mathrm{cav}})({r}_{g}/a){[{(h/r)}^{2}{\alpha }^{\prime }]}^{-1}{L}_{\mathrm{acc}},\end{eqnarray} \tag{ 19 }$$

in which ${L}_{\mathrm{acc}}$ is the luminosity released by accretion onto the BHs, h/r is the aspect ratio of the mini-disks, and ${\alpha }^{\prime }$ is the ratio between vertically integrated and time-averaged accretion stresses and the similarly vertically integrated and time-averaged disk pressure. We write it as ${\alpha }^{\prime }$ rather than the conventional α as a reminder that stresses associated with spiral waves can add to the usual correlated MHD turbulence. We have already found that when a = 20M, ${M}_{\mathrm{slosh}}/{M}_{\mathrm{cav}}\simeq 3\,\times {10}^{-3}$. Although the disk aspect ratio must certainly depend on such parameters as the accretion rate, for the time being we simply note that if the accretion rate is not far below Eddington, when the system is in the PN regime, $h/r\sim {10}^{-2}$ and ${\alpha }^{\prime }\sim {10}^{-1}$ should be reasonably conservative estimates. If so, ${L}_{\mathrm{slosh}}$ might be somewhere in the neighborhood of $\sim {10}^{-1}{L}_{\mathrm{acc}}$. Further work will be required to make sensible predictions about its spectrum.

The optical depth in the sloshing region can be estimated in similar fashion:

$$\begin{eqnarray}&&{\tau }_{\mathrm{slosh}}\sim 4({M}_{\mathrm{slosh}}/{M}_{\mathrm{cav}}){({r}_{g}/a)}^{1/2}{[{(h/r)}^{2}{\alpha }^{\prime }]}^{-1}(\dot{m}/\eta ),\end{eqnarray} \tag{ 20 }$$

where $\dot{m}$ is the ratio of the accretion rate to the Eddington rate and η is the rest-mass efficiency of energy release by accretion. Here, to be consistent with the definition of ${M}_{\mathrm{slosh}}$ found in Section 3.4, we estimate the area occupied by the sloshing region as $\simeq {a}^{2}/4$. Repeating the estimate we just made for the mini-disk internal inflow rate, this expression suggests that the sloshing region should generally be optically thin to Thomson scattering because ${({r}_{g}/a)}^{1/2}\lesssim 1$, while $\dot{m}/\eta$ is likely to be bounded above by ∼10, but could be considerably less. Thus, in many instances, the $\sim 10 \%$ addition to the time-averaged bolometric luminosity should, in fact, be modulated with the period of the sloshing, half the binary orbital period.

The sloshing may also lead to a wholly new kind of mass transfer. In our simulation, with its unity mass ratio, the sloshing is, on average, perfectly symmetric. However, the general case for binaries is a mass ratio different from unity; when this is true, the sloshing is likely to lose its symmetry, so that there is a net mass flow from one mini disk to the other. As a result, mass that leaves the inner edge of the circumbinary disk and arrives at the outer edge of the mini-disk around, for example, the secondary, may, through the sloshing mechanism, find its way to the mini-disk around the primary. This will be an interesting phenomenon for future simulations to explore.

Perhaps somewhat surprisingly, although spiral waves in disks within binaries were first studied in the context of cataclysmic variables (Lynden-Bell & Pringle 1974), disks in SMBBHs whose separations are small enough to be in the relativistic regime may be the environment in which their effects are strongest. The reason is that the angular momentum transport that can be accomplished by spiral shocks increases with disk temperature, in part because a higher ratio of disk sound speed to orbital speed makes the spiral opening angle larger, leading to stronger shocks (Savonije et al. 1994; Ju et al. 2016); disks in SMBBHs are especially hot in this sense. In such disks, if the accretion rate in Eddington units is $\dot{m}\gtrsim 0.01$, the local pressure is dominated by radiation out to $\sim 100M$, and the ratio of the effective sound speed (including radiation pressure) to orbital speed is $\gtrsim 0.1$ for all radii within $\sim 15(\dot{m}/0.1)M$ (Shakura & Sunyaev 1973). Thus, it is possible that spiral shocks, driven by some combination of tidal forces and accretion stream shocks, may be of special interest in relativistic SMBBHs. We caution, however, that there are three-dimensional effects that might complicate the situation (see Section 4.4 for further discussion of this point). We also remark that our simulations described disk thermodynamics in an extremely simplistic fashion, so the details of spiral shock behavior in them should not be taken as predictive of real systems.

The spiral shocks in relativistic SMBBHs may be enhanced in another way as well. As the binary peels streams of gas off the inner edge of the circumbinary disk, a portion of the streams falls into the central cavity and shocks against the outer edges of the mini-disks. These shocks could also contribute to spiral shocks in two ways. First, the accretion stream shocks may heat the mini-disk and therefore enhance the effectiveness of angular momentum transport in spiral shocks. The degree to which the temperature throughout the mini-disk is raised by the accretion shocks will depend on how rapidly the heat they generate is radiated. Some estimates (Roedig et al. 2014) suggest that this reradiation is quite rapid, but the immediate post-shock temperature is so high that the ultimate post-shock temperature might still exceed the temperature expected from ordinary disk dissipation processes. Second, other simulations (Noble et al. 2012; Shi et al. 2012; D'Orazio et al. 2013; Farris et al. 2014) have shown that the rate at which matter is delivered from a circumbinary disk to mini-disks exhibits strong periodic modulation, with a period comparable to the binary orbital period. This periodicity, which does not affect ordinary mass transfer through Roche lobe overflow, could drive the formation of additional spiral shocks.

It is also of interest that spiral shocks in relativistic SMBBHs have a different symmetry from those in Newtonian systems. Previous Newtonian simulations (see, e.g., Savonije et al. 1994; Makita et al. 2000; Ju et al. 2016), as well as analytic theory (Spruit et al. 1987; Savonije et al. 1994; Rafikov 2016), predict the dominant spiral shock mode should be m = 2. So, too, does the recent work of Ryan & MacFadyen (2016), in which mini-disk dynamics driven by the BH at the center of the disk were treated in full general relativity, but the tidal field due to a BH companion was added as a perturbation. By contrast, our work found that the amplitude of surface density modulation associated with m = 1 modes was always a few times greater than that due to m = 2 modes. This discrepancy is almost certainly explained by the fact that in all these previous efforts the companion's tidal field was described in terms of Newtonian gravity, as in Equation (16). In this context, it should be noted that the Newtonian approximation was actually appropriate to the work of Ryan & MacFadyen (2016) because the binary separation in their simulations was fixed at 1000M.

In fact, not only do the m = 1 modes dominate m = 2 for separations as large as 100M, but the ratio of m = 1 to m = 2 grows as the binary separation decreases and general relativistic effects become even more prominent. The ratio might increase still further at very small separation because inspiral, whose timescale is proportional to a⁴, is another source of m = 1 behavior. The inspiral rate due to the loss of energy to GW emission depends on the mass ratio of the binary, and is highest for equal mass binaries. Therefore, the two BHs in the binary approach one another faster than their respective mini-disks; on their own, the mini-disks would inspiral much more slowly due to less efficient GW emission. As a result, at any given time, the BHs are slightly offset from the centers of their own disks, with the magnitude of this offset comparable to the orbital shrinkage during a mini-disk fluid dynamical time. Relative to the size of the mini disk, this offset can become sizable in the late stages of inspiral:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}x}{{r}_{t}}\sim 0.1\displaystyle \frac{{q}^{1/2}}{{(1+q)}^{3/2}}{({r}_{t}/0.3a)}^{1/2}{(a/10M)}^{-5/2}.\end{eqnarray} \tag{ 21 }$$

Analogs to this effect have been seen in other kinds of systems. For example, an inverse version can be observed in self-gravitating disks around single masses: in this case, the central mass is forced onto a spiral trajectory by the gravitational potential of an m = 1 mode in the disk (Adams et al. 1989; Heemskerk et al. 1992; Korobkin et al. 2011; Mewes et al. 2016a, 2016b). Similarly, m = 1 spiral shocks can be created in accretion disks around black holes recoiling after a merger (Corrales et al. 2010; Ponce et al. 2012).

We close this section by pointing out a curious link between spiral shocks, accretion onto the BHs, and sloshing. As shown in Figure 16, the mass in the sloshing region is modulated at two frequencies, $2{\omega }_{\mathrm{bin}}$ and $2.75{\omega }_{\mathrm{bin}}$. Surprisingly, we find (see Figure 19) that the accretion rate onto BH1 (and presumably its twin, BH2) is also modulated, but exclusively at a frequency of $2{\omega }_{\mathrm{bin}}$. We emphasize that this is not the modulation associated with accretion from the inner edge of the circumbinary disk (Noble et al. 2012; Shi et al. 2012; D'Orazio et al. 2013; Farris et al. 2014); it is the accretion rate inside a mini disk wholly isolated from any external circumbinary disk.

Zoom In Zoom Out Reset image size

In our purely hydrodynamic simulations, the only way matter can accrete through a mini-disk is by stresses associated with spiral shocks. However, to the extent that these spiral shocks are well described by their time-averaged structure, there is no reason for them to cause any time-dependence in the accretion rate because, on average, they are stationary in the corotating frame. It is possible, though, that their strength, and therefore the accretion rate, might be modulated due to interaction between the dominant m = 1 and the weaker m = 2 mode. We further speculate, prompted by the $2{\omega }_{\mathrm{bin}}$ modulation of the sloshing mass, that this periodic variation of the spiral shock amplitude is related to periodic behavior in the sloshing region. Some of the angular momentum carried outward by the spiral shocks may be given to sloshing mass; conversely, periodic motions in the sloshing region may drive corresponding changes in amplitude in the spiral shocks.

All our calculations were confined to 2D dynamics in the disk equatorial plane. Real disks are, of course, 3D objects. Although many effects we investigated (e.g., the location of the disks' outer edges) are unlikely to be substantially affected by a change in dimensionality, others are.

Two in particular are worth comment. The first is that the thickness of the accretion streams from the circumbinary disk to the mini-disks relative to the thickness of the mini-disks themselves could be an important parameter regulating how the material in these streams joins the disk and the character of waves launched into the disk by stream impact. For example, if the streams are thicker than the disks, how far inward do the shocks wrap around the disks? When the streams shock against a disk edge, do they rise in temperature sufficiently that they are always thicker than the disk? Similar questions might also apply to the sloshing streams.

The second is that there can be significant 3D effects in the propagation of the spiral waves within the mini-disks. Lubow & Ogilvie (1998) and Ogilvie & Lubow (1999) studied this issue for the case of waves driven by tidal gravity in a binary. When the disk temperature decreases vertically away from the midplane, they found that the spiral waves are channeled toward the surface, but also that a more gradual drop in gas density at the surface can limit the concentration of the waves there. The degree of wave focusing can, of course, have significant implications for where the waves steepen into shocks and dissipate their energy.

Proper resolution of these questions must, of course, await fully 3D calculations.