CIRCUMBINARY MAGNETOHYDRODYNAMIC ACCRETION INTO INSPIRALING BINARY BLACK HOLES

The PN approximation is based on a perturbative expansion of all fields, assuming slow motion v/c ≪ 1 and weak fields GM/(rc²) ≪ 1.⁴ These assumptions allow us to search for solutions that can be expressed as a divergent asymptotic series about a flat Minkowski background spacetime. These perturbations obey differential equations determined by the PN-expanded Einstein field equations. One then solves such equations perturbatively and iteratively to construct an approximate solution.

The two-body problem in the slow-motion/weak-field limit is better understood by classifying the spacetime into different regions, where different assumptions hold and different approximations can be used (see, e.g., Thorne 1980; Alvi 2000, 2003; Yunes et al. 2006; Yunes & Tichy 2006; Yunes 2007; Johnson-McDaniel et al. 2009 for a review). Here we concentrate on the near zone, which is the region sufficiently far from the horizons that the weak-field approximation of PN is valid, but less than a reduced gravitational wave wavelength $\lambdabar$ away from the center of mass of the system, so that retardation effects can be treated perturbatively. We note that in the far zone, i.e., the radiation zone where retardation effects can no longer be treated perturbatively, a multipolar post-Minkowskian expansion can be used rather than a PN one. Very close to each black hole (i.e., in the inner zone), perturbed Schwarzschild solutions are used (which can be extended to include spin by using perturbed Kerr solutions). In this paper, the BBH metric will be approximated with only the near zone solution.

The PN approximation, of course, has its limits: binary systems eventually become so closely separated that a slow-motion/weak-field description is inappropriate. For example, consider the simple case of a test particle spiraling into a non-spinning (Schwarzschild) black hole in a quasi-circular orbit. Eventually, the particle will reach the innermost stable circular orbit, at which point its orbital velocity v/c ∼ 0.41. Clearly, such a velocity is not much less than unity, and thus, the PN approximation need not be an accurate description of the relativistic orbital dynamics. Similarly, when comparable-mass BBHs inspiral, they eventually reach a separation at which the PN approximation is a bad predictor of the dynamics, since the small-velocity/weak-field assumptions are violated.

The determination of the formal region of validity of the PN approximation is crucial, but it can only be assessed when one possesses a more accurate, perhaps numerical, description of the orbital dynamics. This is indeed the case when considering EMRIs, consisting of a stellar-mass compact object spiraling into a supermassive black hole. When considering EMRIs, one can model the spacetime through black hole perturbation theory, i.e., by decomposing the full metric as that of the supermassive black hole plus a perturbation induced by the stellar-mass compact object, without assuming slow motions (see, e.g., Section 4 of Hughes 2009 for a recent review, Mino et al. 1997; Sasaki & Tagoshi 2003; Barack 2009; Poisson et al. 2011 for related topics). To leading order, the orbital dynamics are then described by geodesics of the small object in the spacetime of the supermassive black hole. The orbital motions are slowly perturbed by the radiation reaction due to the emission of gravitational waves.

The comparison of black hole perturbation theory and PN theory predictions has allowed for the construction of different measures to estimate the PN region of validity. When considering the reduction in signal-to-noise ratio induced by filtering an "exact" black hole perturbation theory gravitational wave with a 4PN⁵ order filter, Poisson (1995) found that PN theory is sufficiently accurate provided v ≲ 0.2. When considering the 5.5 and 4PN order predictions for the loss of the binary's binding energy for non-spinning and spinning background black holes, respectively, relative to an "exact" black hole perturbation theory prediction, Yunes & Berti (2008) and Zhang et al. (2011) found that the former is accurate provided v ≲ 0.29, which corresponds to an orbital separation of a ≳ 11 M. The difference between these estimates is due to the different measures used and the different order of the PN approximation employed.⁶

The region of validity of the PN approximation for comparable-mass binaries has not been as well studied. This is because black hole perturbation theory is not applicable here, and one must rely on full numerical relativistic simulations. Currently, state-of-the-art simulations can only model the last few tens of orbits prior to merger, while the determination of the formal region of validity would require knowledge of at least the last thousand orbits. Nonetheless, analytical arguments exist, suggesting that the region of validity in the comparable-mass case is larger than in the EMRI case (i.e., the PN expansion is valid for even larger velocities) (Simone et al. 1997; Blanchet 2003; Mora & Will 2004). Moreover, the NINJA (Numerical INJection Analysis)-2 project (Ajith et al. 2012), a collaboration between numerical relativists and gravitational wave data analysts, has established that certain 3PN order gravitational waveforms are sufficiently accurate for use as templates provided v ≲ 0.33 (a ≳ 8 M).

The PN order of a given term is determined by the exponent of the perturbation parameter contained in that term. In the near zone, the PN expansion of a BBH spacetime metric is a series expansion in the orbital velocity v/c ≪ 1 and the field strength GM_A/(r_Ac²), where M_A and r_A are the masses of the Ath particle and the distances from the Ath particle to a field point, respectively. Here, we may consider 1/c as the PN parameter which goes to zero in the Newtonian limit c → ∞. Note that by the virial theorem $v^2/c^2 = {\cal {O}}[{\it GM}/(a \, c^2)]$. A term proportional to (1/c)ⁿ beyond the Newtonian (leading-order) expression is said to be of the (n/2)th PN order.

The near zone metric will be described here by a resummed PN expression. One begins with the 2.5PN expansion of the metric for non-spinning point particles in a quasi-circular orbit in harmonic coordinates, given for example in Blanchet et al. (1998). Such a metric, however, describes black holes as point particles, which is why one then applies a "background resummation," as in Yunes & Tichy (2006), Yunes (2007), and Johnson-McDaniel et al. (2009). This resummation is intended to improve the strong-field behavior of the metric close to each point particle, i.e., it recovers the horizon of each individual black hole. The metric can then be formally written as

$$\begin{equation} g_{\mu \nu }(t,\,\textbf {\em x}) = g_{\mu \nu }[\textbf {\em x};\,\textbf {\em y}_{A}(t),\,\textbf {\em v}_{A}(t)], \end{equation} \tag{ 1 }$$

where $\textbf {\em x}$ is a spatial vector from the binary's center of mass to a field point, while $\textbf {\em y}_{A}(t)$ and $\textbf {\em v}_{A}(t)$ are the particle's spatial location and 3 velocity with A = (1, 2). This metric depends on the mass of each individual black hole, but also on the binary orbital evolution $\lbrace \textbf {\em y}_{A}(t),\,\textbf {\em v}_{A}(t)\rbrace$ that must be prescribed separately. We will use Greek letters (e.g., μ, ν, λ, κ) to represent spacetime indices [0, 1, 2, 3], and Roman letters (e.g., i, j, k, l) to represent spatial indices [1, 2, 3].

The orbital evolution can also be prescribed within the PN approximation. We simplify the analysis by considering only quasi-circular orbits. This simplification is justified because gravitational wave emission tends to circularize binaries very efficiently, as demonstrated in the weak (Peters 1964) and strong field regimes (Sperhake et al. 2008; Hinder et al. 2008, 2010). We will here use a 3.5PN expansion for the orbital phase evolution ϕ(t), as given for example by Equation (234) of Blanchet (2006), which depends on the quantity Θ = ν(t_c − t)/(5M), where t is time, t_c is the time of coalescence, M is the total mass, and ν = M₁M₂/M² is the symmetric mass ratio. The PN orbital frequency is calculated from Ω_bin = dϕ/dt in this paper.

Although waveforms can be fully characterized by harmonics of the orbital phase, the latter depend on the orbital trajectories. One may model these in harmonic coordinates via

$$\begin{eqnarray} y_{1}^{i}(t) &=& \frac{M_{2}}{M} a(t) [\cos {\phi (t)},\,\sin {\phi (t)},\, 0], \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} y_{2}^{i}(t) &=& {-}\frac{M_{1}}{M} a(t) [\cos {\phi (t)},\,\sin {\phi (t)},\, 0], \end{eqnarray} \tag{ 3 }$$

where $a(t)=|\textbf {\em y}_{1}(t)-\textbf {\em y}_{2}(t)|$ is the orbital separation as a function of time. This separation can be calculated via the balance law (see, e.g., Blanchet 2006), which states that the local rate of change of the binary's orbital binding energy is exactly balanced by the gravitational wave luminosity carried out to future null infinity and into black hole horizons, namely,

$$\begin{equation} { \frac{dE_{\rm Orb}}{dt} = - \cal {L}.} \end{equation} \tag{ 4 }$$

The rate of change of the orbital separation is then given by

$$\begin{equation} { \frac{da}{dt} = - \left(\frac{dE_{\rm Orb}}{da}\right)^{-1} {\cal L}. } \end{equation} \tag{ 5 }$$

Assuming the initial condition a(t = 0) = a₀, we then find

$$\begin{equation} { t = t_c - \int _0^{a} da^{\prime } \left(\frac{dE_{\rm Orb}}{da^{\prime }}\right) \; {\cal L}^{-1}, } \end{equation} \tag{ 6 }$$

where the integrand is expanded in a Taylor PN series and the coalescence time t_c is defined by

$$\begin{equation} { t_c = \int _0^{a_0} da \left(\frac{dE_{\rm Orb}}{da}\right) \; {\cal L}^{-1}. } \end{equation} \tag{ 7 }$$

In this paper, we only use the part of ${\cal {L}}$ that is carried to future null infinity and, although Equations (6) and (7) are typically Taylor expanded to evaluate t(a), here they are inverted through the Newton–Raphson method to yield a(t).

Given the above analysis, we can now calculate the time of coalescence and the number of orbits in each of the simulations carried out. As already discussed, RunSE is artificially kept at a fixed semimajor axis, so a(t) = a₀ for all times, and thus, formally t_c = ∞. On the other hand, RunIn keeps a(t) fixed to a(t) = a₀ = 20 M for t < t_shrink ≡ 40, 000 M, after which it is allowed to evolve according to the PN equations of motion (EOMs). The time of coalescence for this run can be found by evaluating Equation (7) to obtain t_c ∼ 14, 000 M, although to leading order it is approximately described by Peters (1964)

$$\begin{equation} { t_c \sim \frac{5}{256 \nu }\left(\frac{a_0}{M}\right)^4 M = \frac{5}{256}\left(\frac{a_0}{M}\right)^4 \frac{(1+q)^2}{q} M, } \end{equation} \tag{ 8 }$$

where q = M₂/M₁ is the binary mass ratio. A BBH clearly takes longer to merge for systems that start at larger initial separations and that possess extreme mass ratios. Obviously, t_c is always defined as the length of time to coalesce, when the binary is allowed to inspiral. The total simulation time of RunIn is then t_c + t_shrink ∼ 54, 000 M.

We have plotted a few diagnostics to get a sense of the evolution of the binary system in RunIn. Figure 1 shows the orbital evolution of the binary in the x–y plane, after it is allowed to inspiral. Figure 2 plots the orbital separation as a function of time. Observe that initially the semimajor axis is artificially kept fixed, while after t > t_shrink it is allowed to decrease due to gravitational radiation reaction. Figure 3 plots the number of orbits N_orbits traced by the binary system as a function of time, which is given by

$$\begin{eqnarray} &&N_{\rm orbits} = \left\lbrace \begin{array}{@{}ll @{}}\displaystyle \frac{1}{2\pi } \,\Omega _{\rm bin,0}\,t \ ({\rm if} \ t < t_{\rm shrink}), \\ [2ex]{\displaystyle \frac{1}{2\pi }} \, [ \Omega _{\rm bin,0}\,t_{\rm shrink} + \phi (t-t_{\rm shrink})] \ ({\rm if} \ t \ge t_{\rm shrink}). \end{array} \right. \end{eqnarray} \tag{ 9 }$$

We have here defined Ω_{bin, 0} = Ω_bin(r = a₀) to be the (constant) PN orbital frequency at a fixed semimajor axis, and we have set ϕ(0) = 0 for the PN phase evolution. The number of orbits is obviously a piece-wise function since when t < t_shrink, N_orbits increases linearly, as the binary is artificially kept at fixed a₀, while when t > t_shrink, N_orbits can be approximated to leading order by

$$\begin{equation} { N_{\rm orbits} \sim \frac{1}{64 \pi \nu }\left(\frac{a_0}{M}\right)^{5/2} = \frac{1}{64 \pi }\left(\frac{a_0}{M}\right)^{5/2} \frac{(1+q)^2}{q}. } \end{equation} \tag{ 10 }$$

Therefore, for RunIn there are approximately 100 total orbits, while for RunSE there are approximately 127 orbits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since we assume that the gas does not self-gravitate and alter the spacetime dynamics, we need only solve the GRMHD equations on a specified background spacetime, g_μν(x^λ), where {x^λ} represents a set of general spacetime coordinates. The EOM originate from the local conservation of baryon number density, the local conservation of 4-momentum, and the induction equations from Maxwell's equations (please see Noble et al. 2009 for more details). They take the form of a set of conservation laws:

$$\begin{equation} {\partial }_t {\bf U}\left({{\mathbf {P}}}\right) = -{\partial }_i {\bf F}^i\left({{\mathbf {P}}}\right) + \mathbf {S}\left({{\mathbf {P}}}\right), \end{equation} \tag{ 11 }$$

where U is a vector of "conserved" variables, Fⁱ are the fluxes, and $\mathbf {S}$ is a vector of source terms. Explicitly, these are

$$\begin{equation} {\bf U}\,({{\mathbf {P}}}) = \sqrt{-g} [ \rho u^t,\, {T^t}_t + \rho u^t,\, {T^t}_j,\, B^k ]^T \end{equation} \tag{ 12 }$$

$$\begin{equation} {\bf F}^i({{\mathbf {P}}}) = \sqrt{-g} [ \rho u^i,\, {T^i}_t + \rho u^i,\, {T^i}_j,\, (b^i u^k - b^k u^i)]^T \end{equation} \tag{ 13 }$$

$$\begin{equation} \mathbf {S}\,({{\mathbf {P}}}) = \sqrt{-g} [ 0,\, {T^\kappa }_\lambda {\Gamma ^\lambda }_{t \kappa } - \mathcal {F}_t,\, {T^\kappa }_\lambda {\Gamma ^\lambda }_{j \kappa } - \mathcal {F}_j,\, 0 ]^T, \end{equation} \tag{ 14 }$$

where g is the determinant of the metric, Γ^λ_μκ is the metric's affine connection, $B^i = {^{^*}\!\!F}^{it}/\sqrt{4\pi }$ is our magnetic field (proportional to the field measured by observers traveling normal to the spacelike hypersurface), ${^{^*}\!\!F}^{\mu \nu }$ is the Maxwell tensor, u^μ is the fluid's 4-velocity, b^μ = (1/u^t)(δ^μ_ν + u^μu_ν)B^ν is the magnetic 4-vector or the magnetic field projected into the fluid's co-moving frame, and $W = u^t / \sqrt{-g^{tt}}$ is the fluid's Lorentz function. The MHD stress-energy tensor, T_μν, is defined as

$$\begin{equation} T_{\mu \nu } = \left(\rho h + 2 p_{m} \right) u_\mu u_\nu + \left(p + p_{m}\right) g_{\mu \nu } - b_\mu b_\nu , \end{equation} \tag{ 15 }$$

where p_m = b^μb_μ/2 is the magnetic pressure, p is the gas pressure, ρ is the rest-mass density, h = 1 + + p/ρ is the specific enthalpy, and is the specific internal energy, We evolve the quantity (ρu^t + T^t_t) instead of T^t_t in order to reduce the magnitude of the internal energy's numerical error (Gammie et al. 2003). Note that the terms proportional to Γ^λ_tκ and Γ^λ_ϕκ in the source no longer vanish as the metric is now dependent on time and azimuthal coordinate, ϕ. Also, note that we add a negative source term $(-\mathcal {F}_\mu )$ to the local energy conservation equation to model energy/momentum loss from radiative cooling; please see Section 3.4 for more details.

The MHD evolution is facilitated by calculating and using so-called primitive variables: the rest-mass density (ρ), the internal energy density (u = ρ), the velocities relative to the observer moving normal to the spacelike hypersurface, $\tilde{u}^i = u^i - u^t g^{ti} / g^{tt}$. The magnetic field Bⁱ is considered both a primitive and a conserved variable. We employ piecewise parabolic reconstruction of the primitive variables for calculating the local Lax-Friedrichs flux at each cell interface (Gammie et al. 2003). We use a three-dimensional version of the FluxCT to impose the solenoidal constraint, $\partial _i \sqrt{-g} B^i = 0$ (Tóth 2000). The electromotive forces are calculated midway along each cell edge using piecewise parabolic interpolation of the fluxes from the induction equation. A second-order accurate Runge–Kutta method is used to integrate the EOM using the method of lines once the numerical fluxes are found. The primitive variables are found from the conserved variables using the "2D" (two-dimensional) scheme of Noble et al. (2006). A conservation equation for the entropy density is evolved and used to replace the total energy equation of the 2D method whenever the plasma becomes too magnetically dominated or—specifically—when ρ < 0.02p_m; this procedure helps us avoid numerical instabilities and negative pressures from developing. Please see Noble et al. (2009) for more details.

The MHD evolution is performed in the same way as in single black hole cases except that the metric is evaluated⁸ at the present sub-step's time before the MHD field components are updated. The metric is required in many facets of the update procedure. For example, it is used in calculating the 4-velocity from the primitive velocities, source terms, and geometric factors in the EOM, and for deriving the primitive variables from the conserved variables. The affine connection is calculated via finite differencing the PN metric to evaluate

$$\begin{equation} {\Gamma ^\mu }_{\nu \kappa } = \frac{1}{2} g^{\mu \sigma } (\partial _\nu g_{\kappa \sigma } + \partial _\kappa g_{\nu \sigma } - \partial _\sigma g_{\nu \kappa } ). \end{equation} \tag{ 16 }$$

The spatial finite differences use fourth-order centered stencils away from the physical boundaries and backward/forward stencils adjacent to the physical boundaries. Since the metric is evaluated and stored at the cell centers and faces, but the connection is only evaluated at the centers, fourth-order stencils require only three cells' worth of data to compute. The time derivatives are second-order accurate, but use a time spacing 10⁻³ times that used in the MHD integration. This means that additional evaluations of the metric are made at advanced and retarded times at the cell centers to calculate the time derivatives for each connection evaluation. We have verified that the truncation error from the time derivatives is smaller than that from the spatial derivatives for the mesh spacings used in these simulations. Also, the connection's spatial finite differencing is one order more accurate than that of the MHD procedure, implying it is not the primary source of error in the calculation. Please see Appendix B for a discussion on our resolution tests.

In this project, we avoid evolving the gas in the neighborhood of the black holes, choosing instead to focus on establishing reasonable prior conditions for the gas that ultimately feeds the BBH. We therefore excise a spherical domain, which includes the binary, from our calculation.

It is common to begin with initial conditions devoid of large transient artifacts. This is often done by starting from a torus of material in equilibrium (via pressure and rotational support) about the central gravitating source (e.g., a black hole) (Chakrabarti 1985; De Villiers et al. 2003). Unfortunately, such tori will not be near equilibrium in our spacetime as it is (t, ϕ) dependent. Furthermore, the equations describing their structure assume that the metric has the same form (i.e., shares the same zero-valued elements) as the Kerr metric in Boyer–Lindquist coordinates. We resolve these issues in the following way. First, since we hold the binary at fixed separation for several orbits, the spacetime initially has a helical Killing symmetry, with Killing vector $\mathcal {K}^a = \Omega _\mathrm{bin} (\partial _\phi)^a + (\partial _t)^a$. In other words, the spacetime is invariant in a frame rotating with the binary, while the separation is held constant. Since the torus will lie a few a₀ away from the binary, its dynamical response time—comparable to its orbital period—will be longer than the binary period, implying that a torus near equilibrium in this helically symmetric spacetime will also be near equilibrium in its time average. Due to its helical symmetry, its time average is also its azimuthal average. We therefore start with a torus in equilibrium in a (t, ϕ)-independent spacetime, $\hat{g}_{\mu \nu }$, found by averaging over ϕ:

$$\begin{equation} \hat{g}_{\mu \nu } = \frac{ \int g_{\mu \nu } \sqrt{g_{\phi \phi }} \, d\phi }{ \int \sqrt{g_{\phi \phi }} \, d\phi } . \end{equation} \tag{ 17 }$$

We have verified that the same components that are zero-valued in Boyer–Lindquist coordinates are consistent with zero to within our PN-order accuracy in the $\hat{g}_{\mu \nu }$ metric. This means we can employ a similar torus solution method as described in Chakrabarti (1985). A description of our modifications to the procedure—including the generalization to our ϕ-averaged spacetime—is provided in Appendix A. Note that we now ensure that the equilibrium solution is found iteratively to greater precision, instead of the approximate method described in De Villiers et al. (2003) which has been used in prior work of the authors in single black hole disk evolutions (Noble et al. 2009, 2010) and by others studying the hydrodynamic circumbinary case (Farris et al. 2011). We find that our procedure produces initial tori that are much closer to equilibrium than the approximate scheme. Please see Appendix A for more details.

Previous studies have shown that a gap develops near 2.5a for equal mass binaries (MacFadyen & Milosavljević 2008; Shi et al. 2012). We aim to study how this gap develops, so we choose to start material outside this radius. We therefore set up a disk with inner edge located at r_in = 3a₀ and pressure maximum located at r_p = 5a₀; from prior experience, r_p is approximately the radius at which the disk transitions from accreting to decreting since matter must shed its angular momentum to fluid elements further out in order to accrete. These outer elements gain angular momentum and form a time-averaged decretion flow away from the central potential. We will therefore focus on r < r_p = 5a₀ in our analyses. The initial disk extends to r_out ≃ 12a₀, is isentropic with $p/\rho ^\Gamma = 0.01$, and is tuned to have an aspect ratio of H/r = 0.1 at r = r_p, where H is the density scale height defined as the first moment of the rest-mass density with respect to distance from the midplane:

$$\begin{equation} H \equiv \frac{\langle \rho \sqrt{g_{\theta \theta }} \, |\theta - \pi /2| \rangle }{\langle \rho \rangle }, \end{equation} \tag{ 18 }$$

and where the 〈X〉 denotes the average over origin-centered spheres:

$$\begin{equation} \langle X \rangle \equiv \frac{\int X \, \sqrt{-g} \, d\theta d\phi }{\int \sqrt{-g} \, d\theta d\phi }. \end{equation} \tag{ 19 }$$

More information about the initial torus and its solution method is given in Appendix A. In the disk, we add random, cell-scale noise to the internal energy, u, in order to hasten the development of turbulence; the random noise is evenly distributed over the range ±5 × 10⁻³.

Once the torus is in place on the grid, a surrounding nonmagnetized atmosphere is added as our numerical scheme requires us to maintain positive values of ρ and p. The atmosphere is initially static, uⁱ = 0, and in approximate pressure equilibrium: ρ_atm = 1 × 10⁻⁷ρ_max (r/M)^−3/2, u_atm = 3.3 × 10⁻⁶ u_max (r/M)^−5/2, where ρ_max and u_max are—respectively—the initial maxima of ρ and u. We note that when either ρ or u are found to go below, respectively, ρ_atm or u_atm, they are set to those atmosphere values without any modification to the magnetic field or fluid velocity; this happens very rarely once the disk's turbulence saturates.

The magnetic field is initialized as a set of dipolar loops that follow density contours in the disk's interior. We set the azimuthal component of the vector potential and differentiate it to yield Bⁱ; B^ϕ(t = 0) = 0 in our configuration. The vector potential component is

$$\begin{equation} A^\phi = A^\phi _0 \, \max \left[ \left(\rho - \frac{1}{4}\rho _\mathrm{max} \right), 0 \right]. \end{equation} \tag{ 20 }$$

The magnitude of the field, A^ϕ₀, is set such that the ratio of the disk's total internal energy to its total magnetic energy is 100.

The domain on which the MHD EOM are solved is a uniformly discretized space of spatial coordinates {x⁽ⁱ⁾} that are isomorphic to spherical coordinates {r, θ, ϕ}:

$$\begin{equation} r(x^{(1)}) = M e^{x^{(1)}}, \end{equation} \tag{ 21 }$$

$$\begin{eqnarray} &&\theta (x^{(2)}) = \frac{\pi }{2} \bigg[\! 1 + (1\,{-}\,\xi ) (2 x^{(2)}\,{-}\, 1 ) +\left(\xi \,{-}\, \frac{2 \theta _c}{\pi } \right) (2 x^{(2)}\,{-}\, 1 )^n \!\bigg], \nonumber\\ \end{eqnarray} \tag{ 22 }$$

and ϕ = x⁽³⁾. We set n = 9, ξ = 0.87, and θ_c = 0.2. The logarithmic radial coordinates are such that the radial cell extents are smaller at smaller radii in order to resolve smaller scale features of the accretion flow there. The x⁽²⁾↔θ mapping concentrates more cells near the plane of the disk and the binary's orbit, the equator of our coordinate system. Let each grid cell in our numerical domain be labeled by three spatial indices that each cover [0, N⁽ⁿ⁾ − 1], where {N⁽ⁿ⁾} are the number of cell divisions along each dimension. A cell with indices (i, j, k) is located at (x⁽¹⁾_i, x⁽²⁾_j, x⁽³⁾_k), where x⁽ⁿ⁾_j = x⁽ⁿ⁾_b + (j + (1/2))Δx⁽ⁿ⁾. The grid we used is completely specified by: x⁽¹⁾_b = ln (r_min/M), Δx⁽¹⁾ = ln (r_max/r_min)/N⁽¹⁾, r_min = 15 M, r_max = 260 M, N⁽¹⁾ = 300, x⁽²⁾_b = 0, Δx⁽²⁾ = 1/N⁽²⁾, N⁽²⁾ = 160, x⁽³⁾_b = 0, Δx⁽³⁾ = 2π/N⁽³⁾, and N⁽³⁾ = 400.

We chose our resolution and grid extent based upon a number of criteria. First, our θ and ϕ resolutions were set in order to adequately resolve the MRI based on guidelines of Hawley et al. (2011) and Sorathia et al. (2012). The radial resolution was chosen to resolve the spiral density waves—generated by the binary's time-varying tidal field—by several radial zones. We find that our grid adequately resolves the MRI, as measured by the criteria of Hawley et al. (2011) and Sorathia et al. (2012), throughout the domain of interest (i.e., r < 5a₀) ∀t. Please see Appendix B for a quantitative description of these resolution criteria and for a demonstration of how well we resolve the MRI.

The radial extent of the grid was inspired by Shi et al. (2012) and the limits of our near-zone PN metric. The inner edge of the computational domain is placed outside the orbit of the BBH , in order to avoid simulating matter in the immediate vicinity of the black holes. This relaxes the Courant condition on the MHD evolution, allowing us to take larger time steps, and makes the simulations affordable with our available computational resources. The gain in computational efficiency, however, comes with a cost: it excludes any possibility of studying jets in our simulations as jets are launched when large-scale magnetic field is attached to black hole horizons (e.g., Blandford & Znajek 1977; McKinney & Gammie 2004; Hawley & Krolik 2006; Palenzuela et al. 2010b). There is no physical reason why we cannot extend the method described herein to the entire spacetime; such a calculation is left for future work. Since the near-zone PN metric that we use is only valid at distances more than 10M_i from the black hole with mass M_i (Yunes & Tichy 2006; Yunes et al. 2006; Johnson-McDaniel et al. 2009), then—in the equal mass case considered here—we have r_min ⩾ 10(M/2) + max (a(t))/2 = 5M + a₀/2 = 15 M. Shi et al. (2012) found that an inner radial boundary located at r_min ≲ 1.1a was sufficiently far away from the gap and deep in the potential as to not significantly alter the development and evolution of the surface density peak at the edge of the gap. These two constraints justify our choice of r_min = 15M = 0.75a₀ and suggest that our inner boundary condition may begin affecting the gap's evolution when a(t) ≲ r_min/1.1 ≃ 13.6M ≃ 0.68a₀, which occurs after approximately t = 51235 M in RunIn. We set r_max = 260M = 13a₀ to encompass the initial torus.

All cells are advanced in time with the same time increment (Δx⁰ = Δt), which itself changes in time; Δx⁰ is set to 0.45Δt_min, where Δt_min is the shortest cell crossing time of any MHD wave over the entire domain.

Boundary conditions were imposed through assignment of primitive variables and Bⁱ in ghost zones. Outflow boundary conditions are imposed at r = r_min and r = r_max which amounts to extrapolating the primitive variables at zeroth-order into the ghost zones. Additionally, u^r is set to zero—and $\tilde{u}^i$ recalculated—whenever it points into the domain at r = r_max. We note that we attempted to implement a similar condition on u^r at r = r_min, but found it to be unstable during the earliest part of the simulation. Even though it was successfully used in MacFadyen & Milosavljević (2008) and Shi et al. (2012), we found that this condition was inconsistent with the tendency of negative radial pressure gradients developing ahead of each black hole as it moved around its orbit. This pressure gradient moves small amounts of material onto the grid, elevating the density just above the floor ahead of the black holes. Even without the special condition on u^r at r_min and just using zeroth-order extrapolation of the primitive variables there, we observe insignificant amounts (≪1% of the total) of positive mass flux into the domain there and a nearly flat $\dot{M}(r,t)$ profile over r_min < r < 2a with no noticeable artifacts near r_min (e.g., Figure 7).

Depending on internal properties of the gas (e.g., density, accretion rate), the disk may or may not be optically thin, geometrically thin, or have a constant aspect ratio H/r. As these disk characteristics are sensitive to the assumed initial conditions and thermodynamics of the system, we therefore must select which kind of disk to model a priori and verify the consistency of our assumptions a posteriori. We chose to model a disk with intermediate thickness (H/r = 0.1) in order to both address the fact that binary's torque will likely heat the gas efficiently and our expectation that the disk will be dense, optically thick and radiating efficiently.

We chose the ideal-gas "Γ-law" equation of state to close the MHD EOM: P = (Γ − 1)ρ. We set Γ = 5/3, which describes reasonably well the behavior of a plasma whose specific thermal energy is smaller than an electron's rest-mass energy (i.e., it is not relativistically hot). The gas is cooled to a target entropy, the initial entropy of the disk, at a rate equal to

$$\begin{equation} \mathcal {L}_c = \frac{\rho \epsilon }{T_\mathrm{cool}} \left(\frac{\Delta S}{S_0} + \left|\frac{\Delta S}{S_0}\right| \right), \end{equation} \tag{ 23 }$$

where ΔS ≡ S − S₀ and T_cool = 2π(r/M)^3/2 is the cooling time, which we approximate as the Newtonian period of a circular equatorial orbit at radius r. Our procedure is similar to those used by Noble et al. (2009) and Penna et al. (2010). The term in the parentheses acts as a switch ensuring that $\mathcal {L}_c \ge 0$ always, and is zero when the local entropy, $S=p/\rho ^\Gamma$, is below the target entropy, S₀ = 0.01, which is the constant value used in the initial data's torus. Hence, the cooling function should release any heat generated through dissipation since the initial state. We do not cool unbound material—i.e., fluid elements that satisfy (ρh + 2p_m)u_t < −ρ—since we do not want to include cooling that results from application of density or pressure floors. Since $\mathcal {L}_c$ is the cooling rate in the local fluid frame, its implementation in the EOM must be expressed in the Eulerian frame:

$$\begin{equation} \mathcal {F}_\mu = \mathcal {L}_c u_\mu. \end{equation} \tag{ 24 }$$

Another advantage of the cooling function is that it provides us with a proxy for bolometric emissivity that is consistent with the disk's thermodynamics—unlike a posteriori estimates of synchrotron and/or bremsstrahlung luminosity that have typically been made in numerical relativity simulations (e.g., Bode et al. 2010; Farris et al. 2010, 2011). We will use $\mathcal {L}_c$ to make predictions of the total luminosity from circumbinary disks. These predictions are made by integrating $\mathcal {L}_c$ over the domain in the Eulerian frame; we expect to verify their accuracy using full GR ray-tracing in future work.

At the beginning of both simulations, orbital shear transforms part of the radial component of the magnetic field to toroidal, creating a laminar Maxwell stress. Meanwhile, in the same region, the MRI grows, its amplitude exponentially growing on the local dynamical timescale, ≃ 500 M at the initial inner edge of the disk, r = 60 M. The turbulence in the inner disk reaches nonlinear saturation at t ≃ 10,000 M. Under the combined influence of the initial laminar and later turbulent Maxwell stress, matter flows inward (see Figure 4).

4.1.1. Surface Density

Soon after t ≃ 10, 000 M, the inward flow begins to pile up at r ≃ 50 M, between two and three times the binary separation (the dashed line in both panels of Figure 4 marks the location of 2a(t) in order to guide the eye). We define the surface density Σ as

$$\begin{equation} \Sigma (r,\phi) \equiv \int \, d\theta \sqrt{-g} \rho \big/ \sqrt{g_{\phi \phi }(\theta =\pi /2)}; \end{equation} \tag{ 25 }$$

when we quote it as Σ(r) that denotes an azimuthal average of Equation (25). In later discussion, we will sometimes normalize the surface density to Σ₀, the maximum surface density in the initial condition; in code-units Σ₀ = 0.0956. In RunSE, Σ(r ∼ 2a) grows steadily for the duration of the simulation, but after t ≃ 20,000 M, the logarithmic rate of growth (i.e., dln Σ(r)/dt) gradually becomes slower and slower. Because a number of azimuthally averaged properties like Σ(r) all begin to evolve more slowly after t = 40000 M, we call the period from then until the end of RunSE the "secularly evolving epoch." For the same reason, we began the binary orbital evolution of RunIn at that time.

Zoom In Zoom Out Reset image size

Once this secularly evolving state is reached, Σ(r) rises sharply from the inner boundary at r = 16 M to r ≃ 50 M, initially ∝r^2.5, but at late times in RunSE, ∝exp (3r/a) (Figure 5). At first, the azimuthally averaged surface density profile forms a relatively flat plateau at radii greater than 50M ≃ 2.5a, but by t = 30000  M, a distinct local maximum appears at r ≃ 50 M and persists for the remainder of the simulation. This maximum is noticeably asymmetric in the sense that |dΣ/dr| is always considerably smaller in the disk body (i.e., r > 2.5a) than in the gap region inside r = 2a (Figure 5). This behavior resembles closely what has been seen previously in the Newtonian regime (e.g., MacFadyen & Milosavljević 2008; Shi et al. 2012).

Zoom In Zoom Out Reset image size

By construction, the behavior of RunIn is identical to that of RunSE up to t = 40,000 M, when the binary inspiral was begun. In fact, at large radius, the behavior of the surface density profile in RunIn continues to be very similar to that of RunSE even after the binary begins to shrink. Near the surface density peak and at smaller radii, however, things change. In RunIn, the location of the peak moves inward as the binary becomes smaller, and the slope of the disk's inner edge becomes noticeably shallower as the inspiral accelerates. Comparing the curve of the dashed line in the RunIn panel of Figure 4 to the curving edge of the colors denoting higher surface density, one can see that the location of the disk's inner edge follows the evolution of the binary until shortly before the end of the simulation.

However, speaking in terms of azimuthally averaged surface density obscures an important aspect of circumbinary disks: near and inside their inner edges, their structure is generically far from axisymmetric. In Figure 6, we show Σ(r, ϕ) at t = 40000 M. As mentioned previously, at radii smaller than ≃ 2a ≃ 40 M, there is relatively little matter. The reason this gap forms is that, unlike a time-steady, axisymmetric potential, the time-dependent quadrupolar potential of the binary does not conserve either the energy or the angular momentum of test particles. Consequently, closed orbits do not exist, and torques driven by the binary can rapidly expel some matter to the outside, while matter on other trajectories can be forced inward (Shi et al. 2012). As a result, even though the rate at which matter enters the gap is comparable to the outer-disk accretion rate, at any given time, relatively little matter can be found in the region within ≃ 2a, and the matter that is present follows trajectories with little resemblance to stationary circular orbits. Instead, a pair of streams leave the inner edge of the disk and curve inward toward each member of the binary. Part of their flow gains enough angular momentum to return to the disk, but part crosses the inner simulation boundary, traveling toward the domain of the binary.

Zoom In Zoom Out Reset image size

In addition, several tens of binary orbits after matter begins to pile up at r ≃ 2.5a, a distinct "lump" (as Shi et al. 2012 called it) forms in the region of the surface density peak. The density contrast between this lump and adjacent regions grows steadily in time. As shown by Shi et al. (2012), the lump's growth is driven by the mechanics of the streams that fall inward from the edge of the circumbinary disk and are then redirected outward after absorbing additional angular momentum from the binary torques. Thus, despite the relatively slow variation of azimuthally averaged disk properties during this period, the "lump" continues to evolve secularly. Although RunSE was not continued long enough to see this effect, Shi et al. (2012) found that the eccentricity of the lump's orbit also grows slowly.

4.1.2. Accretion Rate, Internal Stresses, and Angular Momentum Budget

It is also useful to characterize the global dynamics of circumbinary disks in terms of the radial dependence of the net mass flow, i.e., the accretion rate as a function of radius. We show in Figure 7 how this quantity changed gradually during the secularly evolving epoch of RunSE by dividing the time from 30, 000 M until the end of the simulation at 76, 000 M into four segments and averaging over each one separately. The accretion rate is constant as a function of radius only inside the gap region, at most times increasing gradually outside r ≃ 2a. During the first part of this period, the accretion rate rises steadily to radii beyond 5a, but after t ≃ 50, 000 M, the accretion rate in the outer disk gradually falls. At the end of the simulation, $\dot{M}(r)$ is actually about a factor of two greater at r ≃ 3a than anywhere else. Averaging over the entire secularly evolving epoch, the rate at which mass passes through the inner boundary is a bit less than half the accretion rate at r = 5a. Although the first analytic theories of circumbinary disks (Pringle 1991) assumed that no accretion would pass the inner edge of such a disk, Newtonian simulations, both purely hydrodynamic (MacFadyen & Milosavljević 2008) and MHD (Shi et al. 2012), have generally seen leakage fractions of a few tens of percent; our fraction is thus only somewhat greater than that previously found.

Zoom In Zoom Out Reset image size

As mentioned earlier, Maxwell stresses due to correlations induced in MHD turbulence by orbital shear dominate angular momentum transport within accretion disks. Because the ratio of Maxwell stress to magnetic pressure, 〈b^(r)b_(ϕ)〉/〈p_m〉, is fixed (Hawley et al. 2011) at ≃ 0.3–0.4 in a point-mass potential (here the notation X^(μ) denotes the magnitude of the μ-component of four-vector X projected into the fluid frame), the stress is linearly proportional to the magnetic pressure. A useful measure of the strength of magnetic effects is therefore the plasma β ≡ 〈p〉/〈p_m〉. In most previous accretion disk simulations, this quantity is ∼100 in the midplane and drops to ∼O(1) a few scale heights out of the plane. We show its dependence on position in the poloidal plane at several times during the secularly evolving epoch of RunSE in Figure 8; to be more precise, we show the ratio of the time- and azimuthally averaged gas pressure to the similarly averaged magnetic pressure. As that figure illustrates, the level of magnetization is rather larger than usual (i.e., β is smaller than usual), but gradually diminishes over time. At t = 30, 000 M, β ≃ 1 in the midplane at r ∼ 3–5a and ≃ 3 in the region of the surface density peak (r ∼ 2–3a); by the end of the simulation, it is ≳ 10 in the disk body for the whole range 2a < r < 5a and reaches as much as ≃ 30 in the lump.

Zoom In Zoom Out Reset image size

Ever since the work of Shakura & Sunyaev (1976), it has been popular to measure the vertically integrated, azimuthally averaged, and time-averaged internal disk stress in units of the similarly integrated and averaged pressure. In order to avoid unphysical pressures found in the unbound regions, we computed the stresses and pressures only in bound material. Outside r ∼ 4a, where the disk resembles an ordinary accretion disk, we find that the Maxwell stress alone has magnitude ≃ 0.3–0.5 in these units. This is roughly three to five times larger than the stress levels found in general relativistic simulations of MHD flows in the Kerr metric (Krolik et al. 2005). In the gap region, the ratio of Maxwell stress to pressure rises about a factor of two, while the Reynolds stress in the gap rises dramatically (as also found by Shi et al. 2012). These large Reynolds stresses are entirely due to the strong binary torque, which pushes part of the inflowing streams back out to the disk with additional angular momentum.

An overview of angular momentum flow in the system can be gleaned from Figure 9, in which we show the radial derivatives of the time-averaged angular momentum fluxes integrated on shells, i.e., the time-averaged torques due to the several mechanisms acting. Several important points stand out in this figure. The first is that the binary torques are delivered primarily in the gap region a ≲ r ≲ 2a. The torque density dT/dr peaks at r ≃ 1.45a, and the region surrounding that peak dominates the integral over all radii. Moreover, all these torques are positive in net, but they are locally negative both at small radii (r ≲ a) and at large ones (r ≳ 1.9a). Thus, most of the angular momentum the binary gives the disk is delivered in the gap, where the gas density is very much lower than in the disk proper. This point has previously been emphasized by Shi et al. (2012). Second that angular momentum is conveyed to the disk proper by fluid flows, i.e., Reynolds stresses. That is why the Reynolds stress is large and positive from r ≃ 1.8a to r ≃ 2.5a. Outside those regions, Maxwell stress, which always acts so as to remove angular momentum from the gas and carry it outward, dominates the internal stresses. Finally, the net angular momentum change at any given radius is generally positive in the inner disk because matter continues to pile up between r ≃ 2a and r ≃ 5a throughout the simulation.

Zoom In Zoom Out Reset image size

4.1.3. Disk Thickness

At t = 40, 000 M in RunIn, we began to evolve the binary orbit, letting it compress as gravitational radiation removes its orbital energy. The rate of orbital evolution is extremely sensitive to separation: $\dot{a}/a \propto a^{-4}$ when a/r_g ≫ 1. Consequently, even at the relatively small initial separation assumed here (a₀ = 20 M), orbital evolution is comparatively slow at first. However, it accelerates dramatically after t ≃ 50, 000 M. By the end of RunIn (t = 54, 000 M), $\dot{a}/a$ is quite rapid, and a ≃ 8 M, small enough to make our PN expansion problematic.

While the binary orbit changes relatively slowly, the inner edge of the disk moves inward in pace with the change in the binary separation, staying close to ≃ 2a(t) (as shown in Figures 4 and 10) until t ≃ 50, 000 M. However, as the orbital evolution becomes more rapid (after t ≃ 50, 000 M), although the inner edge of the disk continues to move inward in terms of absolute distance (Figure 4), it begins to recede in terms of r/a(t) (Figure 10). At the end of the simulation, the disk edge has moved in to ≃ 20 M, but that is ≃ 2.5a. Simultaneous with this evolution, the slope of the inner edge also becomes gentler (Figure 5). In other words, the contrast between the surface density in the disk body and in the gap weakens, particularly when considering the outer part of the gap. As shown by the RunSE panel in Figure 10, none of this adjustment (in r/a(t) terms) occurs without binary evolution.

Zoom In Zoom Out Reset image size

Another view of this process may be seen in Figure 11. In that figure, we see the way matter accumulates in the inner disk over time, at first during the secularly evolving epoch and later during the binary orbital evolution of RunIn. The left-hand panel shows what happens when referred to an absolute radius scale. When the binary begins to shrink, the quantity of matter found at small radii grows abruptly, particularly in the original gap region: the amount of mass inside r = 40 M almost doubles, and the mass inside r = 20 M increases by a factor of five during the period of binary orbital evolution. The right-hand panel shows the same events from a different point of view. In this figure, we see that the mass enclosed within small multiples of a(t) declines rapidly as the binary's shrinkage accelerates. For larger multiples (e.g., 3a and 4a), the mass enclosed continues to rise for a while after binary orbital evolution, but eventually drops once the compression becomes rapid. In particular, the mass within the gap region (i.e., r < 2a(t)) falls by roughly a factor of 40 during the period of orbital evolution, although this ratio is in fact a bit ill defined because 2a(t) is almost at the simulation's inner boundary by the end of the simulation.

Zoom In Zoom Out Reset image size

The accretion rate behaves differently. It falls (see Figure 12) from ≃ 30, 000 M to 40, 000 M, even before the binary begins to compress. Without binary orbital evolution (RunSE), it levels out from ≃ 40, 000 M to 50, 000 M, before declining more gradually from ≃ 50, 000 M until the end of RunSE at ≃ 76, 000 M. In RunIn, the onset of binary evolution at t = 40, 000 M leads to a continuing decrease in the rate at which mass flows through the inner boundary that levels out only after ≃ 50, 000 M. Although the accretion rates in the simulations with and without binary orbital evolution decline at different times and at different rates, the final accretion rate in RunIn, when the binary separation has shrunk to 8 M, is only 20%–30% less than at the same time in RunSE.

Zoom In Zoom Out Reset image size

Another consequence of the changing relationship between disk material and the binary is a diminution in the integrated torque when the binary compresses (Figure 13). During the initial slow stages of energy loss due to gravitational wave emission, the binary continues to exert nearly as much torque on the disk as in RunSE, in which the binary orbit does not change at all. However, once the orbital shrinkage begins to accelerate, the torque plummets; at the end of RunIn, it has fallen to ≃ 1/5 of the value at that time in RunSE. The greater part of this diminution in torque is due to the fact that at this stage in the binary's evolution, its separation diminishes so rapidly that the region between a and 2a, where most of the torque is expressed, moves inward faster than the matter can follow. There is consequently much less matter on which these torques can be exerted. The connection between available matter and torque is shown clearly in the right-hand panel of Figure 13, in which one can easily see that for nearly the entire inspiral the torque density at the location of its maximum (r = 1.45a(t)) is almost exactly proportional to the surface density there. However, there is also a smaller part due to an artifact of the simulation. Its inner boundary lies at r_min = 15 M. As soon as a(t) becomes smaller than 15 M, part of the region in which the binary torque is applied is no longer in the problem volume, so we cannot calculate any torque occurring there. As shown in the right-hand panel of Figure 13, this effect becomes significant at t ≃ 5.2 × 10⁴ M, when a(t) ≃ 13 M. By the end of RunIn, a ≃ 8 M, so that nearly the entire region where the torque is exerted (a ≲ r ≲ 2a) has left the problem volume. At that point, even if there were significant matter there, our calculation can say neither what its mass is nor what torque it feels.

Zoom In Zoom Out Reset image size

We define the (Eulerian frame) cooling rate per unit radius of the disk by

$$\begin{equation} \frac{dL}{dr} = \int \, \sqrt{-g} \, d\theta \, d\phi \mathcal {L}_c u_t. \end{equation} \tag{ 26 }$$

During the approximate stationary state, it is best described in terms of two separate regimes. As shown in Figure 14, at large radius (r ≳ 2a), it is very well described by a power law, dL/d(r/a₀) ≃ 5 × 10⁻⁴(r/a₀)⁻²Σ₀a₀. At around r ≃ 2a, the cooling rate per unit radius reaches a local maximum and declines inward. This distinction neatly corresponds to two different mechanisms for generating the requisite heat: the dissipation of MHD turbulence associated with mass accretion (at large radius) and the dissipation of fluid kinetic energy given to the relatively small amount of gas in the gap by the binary torques (at small radius). In fact, this identification is confirmed semi-quantitatively. In time-steady accretion, the luminosity per unit radius is $(3/2)\dot{M} c^4/[(r/r_g)^2 {\it GM}]$ at radii where the local orbital angular momentum per unit mass is large compared to the net angular momentum flux per unit mass. Our disk is never in inflow equilibrium, and this expression is not exact when $\dot{M}$ is a function of radius. Nonetheless, taking it as an estimator, it predicts

$$\begin{equation} \frac{dL}{dr/a_0} = 4 \times 10^{-4} (\dot{M}/0.01) (r/a_0)^{-2} \Sigma _0 a_0. \end{equation} \tag{ 27 }$$

As Figure 7 shows, the mean accretion rate in code units at r = 2a in RunSE was ≃ 0.01, while $\dot{M}$ at larger radii is typically similar or perhaps a factor of two greater. Thus, this prediction of the luminosity profile on the basis of the time-averaged accretion rate and expectations derived from time-steady accretion onto a solitary mass quite accurately matches the actual luminosity profile seen in the simulation.

Zoom In Zoom Out Reset image size

Integrated over radius, the total luminosity reaches a peak $\hat{L} \simeq 5.5 \times 10^{-3}$ at t ≃ 33, 000 M (Figure 15), where $\hat{L}$ is the integrated luminosity in units of GMΣ₀c. After reaching this peak, $\hat{L}$ falls slowly, reaching ≃ 3 × 10⁻³ at t ≃ 76, 000 M in RunSE; averaged over the entire secularly evolving period in this simulation, it is 3.8 × 10⁻³.

Zoom In Zoom Out Reset image size

The light output from RunIn remains very close to that in RunSE until the binary orbital evolution becomes rapid at t ≃ 50, 000 M. After that time, it falls more sharply, so that by the time at which RunIn stops, $\hat{L} \simeq 2.7 \times 10^{-3}$; this is, however, still two-thirds the luminosity in RunSE at the same time. As the binary shrinks, the radial distribution of the luminosity changes in parallel, with the peak in surface brightness moving inward. We attribute the gradual decline in luminosity to the gradual decline in accretion rate. The sharp drop in the final stages of binary orbital shrinkage is due to the interaction of a boundary effect with genuine dynamics. As shown by Shi et al. (2012), gas streams flow inward from the inner edge of a secularly evolving circumbinary disk to radii ≃ 1.2a, where they can be strongly torqued and some of their material flung back outward toward the disk. The outward-moving matter shocks against the disk proper at a radius near that of the surface density peak, and the heat dissipated in these shocks contributes significantly to the luminosity. When the binary shrinks, this mechanism is weakened for two reasons. The inner boundary of our simulation (r = 0.8a₀) eventually becomes larger than 1.2a(t); when it does, matter is no longer thrown outward by binary torques. At the same time, however, it is possible that the retreat of the disk's inner edge when measured in terms of a(t) might also lead to weaker inward streams.

The fact that the energy deposited by binary torques is ultimately radiated in the disk proper leads to a method of estimating the relative contributions to the total luminosity coming from accretion and binary torques. For that reason, and also because the accretion rate diminishes as the region of the surface density peak is approached from larger radius, it is a reasonable approximation to suppose that most of the luminosity from the region of the surface density peak inward has its source in the binary torques. We can therefore estimate the work done by the torques by bounding it between L(r < 2a₀) and L(r < 3a₀). On this basis, accretion would account for ≃ 1/2–3/4 of the total (i.e., $\hat{L} \simeq 1.8$–2.9 × 10⁻³) and the binary torque for ≃ 1/4–1/2 ($\hat{L}\simeq 0.9$–2 × 10⁻³).

The rest-mass efficiency of this luminosity is comparable to the rest-mass efficiency due to accretion that goes all the way to the black hole. Measured in terms of the time-dependent luminosity relative to the time-averaged accretion rate through the inner boundary, the efficiency in RunSE falls from a peak ≃ 0.06 achieved for 20, 000 M ≲ t ≲ 45, 000 M to ≃ 0.03 at the end of this simulation. There are several reasons that this efficiency is so great even though the potential at r = 50 M is an order of magnitude shallower than the potential at the innermost stable circular orbit (the "ISCO"). One is that the accretion rate in the circumbinary disk is roughly twice the accretion rate through the inner boundary, so the local accretion dissipation in the disk is boosted by that same factor of two relative to the rate at which mass passes the inner boundary. Another is that in a conventional disk around a single black hole the dissipation rate in the region just outside the ISCO is depressed relative to larger radii because some of the potential energy released is transported outward by the inter-ring stresses. In the Novikov–Thorne model (in which the stresses are assumed to vanish at the ISCO), almost 40% of the total luminosity is released outside r = 40 M when the black hole has no spin. This fraction is smaller when the spin is greater, and may be further reduced to the degree that the net angular momentum flux is smaller (Krolik et al. 2005). Last, of course, additional energy is deposited in the disk by the work done by the binary torques.

Translating the peak cooling rate into physical units gives

$$\begin{equation} L_{\rm disk} \simeq 2.4 \times 10^{40} (\hat{L}/10^{-3}) M_6 \tau _0{\rm \,erg\,s^{-1}}. \end{equation} \tag{ 28 }$$

Here, τ₀ is the Thomson optical depth through a disk of surface density Σ₀, $\hat{L}$ is the luminosity in code units, i.e., 3–5 × 10⁻³, and M₆ ≡ M/(10⁶M_☉). In Eddington units, this becomes $L_{\rm disk}/L_E \simeq 1.7 \times 10^{-4} (\hat{L}/10^{-3}) \tau _0$. Thus, for such a system to be readily observable at cosmological distances, it will be necessary both for the disk to be optically thick to Thomson scattering and for the mass of the binary to be relatively large. As a gauge of what might reasonably be expected, we note that in a steady-state accretion disk around a solitary black hole, the optical depth of the disk at r/r_g = 20 would be ${\sim} 2 \times 10^3 (\alpha /0.1)^{-1} (\eta /\dot{m})$, where η is the usual rest-mass efficiency and $\dot{m}$ is the accretion rate in Eddington units. With this disk surface density, the luminosity would approach that of a typical AGN when M₆ is at least ∼1.

If this light were radiated thermally, the corresponding effective temperature would be

$$\begin{equation} T_{\rm eff} \simeq 4 \times 10^4 (\hat{L}/10^{-3})^{1/4} M_6^{-1/4} \tau _0^{1/4}\hbox{ K}, \end{equation} \tag{ 29 }$$

where we have assumed that the radiating area is 2π(2a)². Thus, it would emerge primarily in the ultraviolet for fiducial values of black hole mass and optical depth.

The luminosity (assumed to pass freely through the disk) exhibits a noticeable modulation as a function of time, with peak-to-trough contrast of ≃ 5%. Its Fourier power spectrum shows a strong, sharp peak at a frequency 1.47Ω_bin (see Figure 16) and a weaker peak at 0.26Ω_bin. The latter is the orbital frequency at the radius of the surface density maximum, ≃ 2.4a; because the lump is located at this radius, we call this frequency Ω_lump. The former we identify with the rate at which the lump approaches the orbital phase of a member of the binary, 2(Ω_bin − Ω_lump) = 1.46Ω_bin. The peak at this beat frequency has high statistical significance: it is approximately 54 times the level of the fluctuation power at nearby frequencies. In addition, its stable periodicity is confirmed by the fact that it persists over at least ∼25 oscillation periods. It also has a natural physical interpretation. When the lump draws near one of the black holes, a new stream forms, falls inward, and is split into two pieces, one of which gains angular momentum, sweeps back out to the disk, and ultimately shocks against the disk gas. It is this process, whose frequency is 1.46Ω_bin, that modulates the light curve. If the binary mass ratio were far from unity, we expect that the modulation frequency would fall to ≃ Ω_bin − Ω_lump.

Zoom In Zoom Out Reset image size

In many respects, the behavior we found in this PN regime resembles what was previously found in the Newtonian limit (Shi et al. 2012). There is very good agreement in the shapes of their azimuthally averaged surface density profiles, with any contrasts attributable to their somewhat different initial conditions. In both cases, during the secularly evolving epoch the surface density at the disk's inner edge rises ∝exp (3r/a), reaches a maximum at r ≃ 2.5a, and then declines to larger radii.

At early times in both, there is a pair of streams leading from the disk edge to the inner boundary, which they typically reach at an orbital phase slightly ahead of the nearest member of the binary. Both also develop a strong m = 1 asymmetry (a "lump") in the surface density at r ≃ 2.5a at late times. This asymmetry ultimately causes, in both the Newtonian and post-Newtonian simulations, a single stream in the gap to become dominant. Almost the only contrast in this regard is that the orbit of the lump developed a growing eccentricity in the Newtonian case, but not in RunSE.

The level of magnetization is likewise qualitatively similar: the mean plasma β in the Newtonian case fell from ∼1 at ≃ 6a to ≃ 0.3 at r ≃ 2a, while the value (averaged over the secularly evolving epoch in RunSE) in our simulations was ≃ 1.5 at r = 6a, grew to ≃ 2.5 at the surface density peak, and then decreased inward. The magnetic stress-to-pressure ratio α in the disk body follows the same pattern of close resemblance. It was ≃ 0.3 in the disk body in the Newtonian case and ≃ 0.2 in RunSE. In the gap, the similarity was more qualitative than quantitative: in both cases, it rose steeply into the gap, but reached only ≃ 0.7 at r ≃ a in the PN simulation, whereas it climbed to ≃ 10 in the Newtonian one.

Most strikingly, the luminosity estimated by Shi et al. (2012) scales extremely well to the PN case. Shi et al. (2012) could not directly compute the luminosity because they assumed an isothermal equation of state. However, they argued that the work done by the binary torques would be delivered to the disk and ultimately dissipated there into heat. Rewriting in our units their value for the rate at which the torques did work on the gas gives a luminosity of 0.018GMΣ_pc(a/r_g)^−1/2, where Σ_p is the surface density at the maximum; for our separation (a = 20M) and our surface density at the maximum (≃ 0.55 averaged over the secularly evolving epoch in RunSE) that becomes 2.2 × 10⁻³GMΣ₀c. This prediction agrees well with the upper end of our estimated range for the binary torque share of the luminosity.

Milosavljević & Phinney (2005) predicted that at some point well before the merger, the BBH should begin compressing so fast by gravitational radiation that internal stresses within the disk would not allow it to move inward rapidly enough to stay near the binary. At the order of magnitude level, this breakaway point would be expected to come when the gravitational radiation time

$$\begin{equation} t_{\rm gr} = \frac{5}{64}\left(\frac{a}{M}\right)^4 \frac{(1+q)^2}{q} M \end{equation} \tag{ 30 }$$

becomes shorter than the characteristic disk inflow time

$$\begin{eqnarray} t_{\rm in} &=& \alpha ^{-1}(H/r)^{-2}(d\ln \Sigma /d\ln r)^{-1}\Omega ^{-1} \nonumber \\ &=& \alpha ^{-1}(H/r)^{-2}(d\ln \Sigma /d\ln r)^{-1}(r/r_g)^{3/2}M. \end{eqnarray} \tag{ 31 }$$

In these equations Ω is the local disk orbital frequency. The logarithmic derivative of the surface density enters because spreading of the inner edge is more rapid when it is especially sharp.

With a typical estimate of the stress level, α ∼ 0.01, the binary separation at which t_gr and t_in would match, and the disk and binary might decouple, is

$$\begin{equation} a_{\rm dec} = 70 (d\ln \Sigma /d\ln r)^{-2/5}\left(\frac{H/r}{0.15}\right)^{-4/5} M, \end{equation} \tag{ 32 }$$

where the fiducial radius at which the inflow time was computed is r_* = 2a, and we set q = 1. Indeed, it was this sort of estimate that led us to choose the initial conditions for our simulation.

However, scaling to the actual parameters of our simulation leads to a considerably smaller predicted value,

$$\begin{eqnarray} &&a_{\rm dec} \simeq 10 [(d\ln \Sigma /d\ln r)/6]^{-2/5}\left(\frac{\alpha }{0.2}\right)^{-2/5}\left(\frac{H/r}{0.15}\right)^{-4/5}M, \nonumber \\ \end{eqnarray} \tag{ 33 }$$

which is much closer to what is found in RunIn. Thus, the substantially stronger magnetic stresses than predicted by usual α-based estimates lead to decoupling at a much smaller binary separation. Nonetheless, in the end the inward motion of the disk is limited by angular momentum transport, so once the magnitude of those stresses are known, a_dec can be estimated quite well by this means.

On the other hand, the meaning of the term "binary runaway" should also be made more nuanced. As we have seen, the accretion rate through the simulation inner boundary decreases as the binary shrinks, but almost the same decrease in accretion rate occurs when the binary does not shrink. Moreover, the continuing advance of the disk during the period of orbital evolution brought matter rapidly inward. The amount of matter within 30 M rose by about a factor of four while the binary shrank, so that almost as much matter could be found within that radius as had been within 40 M at the beginning of the binary orbital evolution. Thus, acceleration of binary orbital evolution does lead to a state in which the surface density at r < 3a is smaller than would be expected if the orbital evolution were slower, and this diminution in the mass close to the binary does lead to consequences such as sharply diminished torque (as discussed in Section 4.2) and luminosity (see Section 4.3). On the other hand, neither the torque nor the luminosity falls by as much as an order of magnitude because the decoupling of binary and disk matter is not complete. Most notably, accretion continues at a rate only tens of percent lower than in the absence of inspiral.

Continuing accretion into the binary orbital region has a particularly interesting consequence. The material in those streams should, just as happens when the binary orbital evolution is slower, be captured into orbit around one or the other of the members of the binary. It will then settle into two smaller disks, one around each black hole. In the conditions of our simulation, the inflow time in the individual black hole disks might be only slightly shorter than the merger time because decoupling occurs when the binary separation is not a great deal larger than the ISCO, even if both black holes spin rapidly. If the circumbinary disk were cooler than in our simulation, so that a_dec ≫ M, the inflow time in the smaller disks would be shorter than that in the circumbinary disk by a sizable ratio: ∼3 × 10⁻⁴ at decoupling if the scale heights in the inner disks and the outer disk are the same, and even smaller at later times. In either case, there could be interesting hydrodynamic interaction between the two smaller disks as the binary compresses.

In terms of an ultimate comparison to observations, a larger question is posed by the disk's thermal state. As just shown, the binary separation at decoupling scales as (H/r)^−4/5, so the disk's internal pressure (H∝c_s∝(p/ρ)^1/2) can influence a_dec. For reasons of numerical convenience, we chose parameters yielding a relatively thick disk. Although the factors controlling the saturation of magneto-rotational turbulence are still not well understood, a scaling with disk pressure remains plausible. If the effective sound speed of the gas were lower, decoupling might occur at rather larger binary separation, well outside the domain of relativistic orbits.

Several factors can influence the actual equation of state of the disk. In ordinary AGNs, local heating due to accretion can make the disk radiation dominated inside r ≃ 100 M when $\dot{m} \gtrsim 0.3$ and the central mass M > 10⁶ M_☉ (Krolik 1999). Larger central masses lead to radiation dominance even when $\dot{m}$ is smaller. When the disk surrounds a binary, the local heating should be similar at radii r ≳ 2a, as shown by Figure 14; the diminished accretion inside ∼4a is compensated by dissipation of the work done by the binary torques and delivered to the disk. In those circumstances, the disk scale height for r ≳ 2a is independent of radius, giving an aspect ratio $H/r \simeq (3/2) (\dot{m}/\eta) (r/r_g)^{-1}$, where $\dot{m}$ is now the local accretion rate in Eddington units. Thus, our aspect ratio of ≃ 0.15 would correspond to a nominal $\dot{m} \simeq 0.4(r/40r_g)$; that is, this $\dot{m}(r)$ is the mass accretion rate at r that would, if it reached the black hole in a flow with η = 0.1, produce that fraction of an Eddington luminosity.

Given the relatively large leakage fraction through the inner edge of the circumbinary disk, the accretion rate onto the two black holes will in general be smaller than the accretion rate in the disk by only a factor of a few. They might therefore generate a sizable luminosity with a spectrum not too different from that of a generic AGN. Because the density in the gap is considerably smaller than in the disk, this luminosity may irradiate the disk, particularly if it is relatively thick. The inner edge of the disk could be heated by this means, as well.

Thus, considerable uncertainty remains in the thickness profile of a circumbinary disk surrounding a relatively compact BBH. The particular thickness we have simulated lies within that range of uncertainty, but may not be generic.

Given a sufficient external mass supply rate, the luminosity from the circumbinary disk alone could be great enough to be detected, even from a cosmological distance (Section 4.3). However, until the binary separation becomes as small as that considered here, the luminosity from matter accreting onto the two individual black holes would dominate the circumbinary luminosity by a large ratio. In many respects, a BBH with a ≫ 10 M should strongly resemble a conventional AGN. Although truncated at its inner edge, the bulk of the circumbinary disk is very much like an ordinary accretion disk, and there may exist smaller conventional accretion disks around each of the black holes. However, there are also several significant contrasts between this situation and that of an ordinary AGN. Near r ≃ 2a, there is additional heating in the circumbinary disk due to the binary torques. On the other hand, in the gap region, from r ≃ 2a inward to the outer edges of the individual black hole disks (≃ a/3 from each black hole), the dissipation rate is likely to be much smaller than in a conventional disk because the flow is laminar, not turbulent; moreover, the small surface density in that region means that whatever light does issue from that region is unlikely to be effectively thermalized, and would therefore emerge at considerably shorter wavelengths. Within the small disks around each black hole, one might expect behavior qualitatively, but not quantitatively (Beloborodov & Illarionov 2001), similar to the innermost portion of a conventional disk. For the parameters of this simulation (a = 20 M, equal-mass Schwarzschild black holes), however, the fraction of the total light emitted that close in is no more than 10% of the total. Nonetheless, whatever light is emitted in those individual disks would sample the hottest temperatures found in the flow. Thus, one might expect an overall spectrum roughly similar to that of an ordinary AGN, but the portions intermediate in wavelength between those characteristic of the innermost part of the disk and those characteristic of r ≃ 2a would be significantly suppressed, with the equivalent energy appearing at much shorter wavelengths.

The periodic modulation in the heating rate of the circumbinary disk that we have found might make its emission easier to isolate. The key question governing that "might" is how effectively optical depth in the disk blurs the modulation. Our cooling function, which operates at a characteristic rate ∼Ω(r), filters out variations on timescales ≪Ω⁻¹, but optical depth would impose a rather more severe upper bound on the maximum effective frequency of variation. According to Kylafis & Klimis (1987), a constant-density sphere of radius R and optical depth τ suppresses the amplitude of a periodic signal of angular frequency ω injected at its center by a factor ≃ 3c/(Rτω) when ωRτ/c ≫ 1. To apply this estimator, we suppose that a stratified disk segment with scale height H can be approximated by a homogeneous sphere of radius R = H. For binary separation a, the relevant orbital radius is ≃ 2.5a and the signal frequency ω = 1.46Ω_bin ≃ 1.5(GM/a³)^1/2 for total binary mass M. As shown in Figure 6, the surface density in the lump grows to be ≃ Σ₀, so we take τ = τ₀. The suppression factor is then ≃ 0.024(τ₀/1000)⁻¹(a/20r_g)^1/2(H/0.15r)⁻¹. In other words, when the relevant region of the disk is optically thick, the luminosity in the modulated component is independent of surface density, so that its fractional modulation decreases as the luminosity increases.

As we have discussed in Section 5.3, there is also considerable uncertainty in H/r, even given a disk surface density. If the disk thickness and optical depth were determined by the considerations of conventional time-steady accretion flows around single black holes, the characteristic cooling time τH/c can be identified with (αΩ)⁻¹. The fluctuation suppression factor could then be estimated by ≃ 3αΩ/ω. Here, the relevant orbital frequency is Ω(2.4a) = 0.26Ω_bin, while 0.74Ω_bin ⩽ ω ⩽ 1.46Ω_bin (the upper limit applies in the equal-mass case, the lower limit when the masses are very unequal). Thus, the suppression factor estimated in this way is ≃ 0.1(α/0.2) for equal black hole masses, rising to double that in the limit of very unequal masses. However, this estimate is made uncertain by the fact that the disk in the vicinity of the surface density peak is certainly not in a state of inflow equilibrium. In addition, just as for the other estimate, the association of the periodic modulation with the lump means that any estimate based on assumptions of axisymmetry likely underestimates the local optical depth.

Both estimates suggest that the modulation will be suppressed by at least a factor of several, but both are also subject to considerable uncertainty, making the actual outcome unclear. It is worth pointing out that in the event the modulation is detectable, the period of the modulation would allow an estimate of the binary orbital period. When the binary mass ratio is unity, the binary orbital frequency is 0.68 times the frequency of the modulation; as the mass ratio departs from unity, the binary orbital frequency should rise toward ≃ 1.36 times the modulation frequency.

Finally, we remark that our predictions of EM signals from circumbinary disks around merging black holes are complementary to those previously made (Bode et al. 2012) in two ways. In the previous work, the period of EM emission began when the binary separation shrank to 8 M; that is when our calculation ends. In addition that effort expressly excluded the disk proper, which they defined as r ⩾ 16 M; in our work that is the location of the overwhelming majority of the emission. Our effort also differs from previous work in this area in that we explicitly include radiation losses in the gas's energy equation and also tie the rate of radiation directly to the instantaneous local thermodynamic state of the gas (albeit in only a formal way). In addition, our discussion of the observability of periodic modulation in the light curve takes into account possible suppression of the variation due to optical depth in the source.

Here, we describe a method for calculating axisymmetric non-magnetized gas distributions supported by pressure gradients and rotation within an axisymmetric spacetime. We will assume a general spacetime with Killing vectors (∂/∂ϕ)^a and (∂/∂t)^a that can be expressed in the simple form (in coordinates similar to spherical Boyer–Lindquist coordinates):

$$\begin{equation} g_{\mu \nu } = \left[\begin{array}{@{} cccc @{}}g_{t t} & 0 &0 &g_{t \phi } \\ [0.25cm] 0 & g_{r r} & 0 & 0 \\ [0.25cm] 0 & 0 & g_{\theta \theta } & 0 \\ [0.25cm] g_{t \phi } & 0 & 0 & g_{\phi \phi } \\ [0.25cm] \end{array} \right], \end{equation} \tag{ A1 }$$

which means that the inverse metric is

$$\begin{equation} g^{\mu \nu } = \left[\begin{array}{@{} cccc @{}}-\frac{g_{\phi \phi }}{A} & 0 &0 & \frac{g_{t \phi }}{A} \\ [0.25cm] 0 & \frac{1}{g_{r r}} & 0 & 0 \\ [0.25cm] 0 & 0 & \frac{1}{g_{\theta \theta }} & 0 \\ [0.25cm] \frac{g_{t \phi }}{A} & 0 & 0 & -\frac{g_{t t}}{A} \\ [0.25cm] \end{array} \right], \end{equation} \tag{ A2 }$$

where A = g²_tϕ − g_tt g_ϕϕ. We have verified that the ϕ-average of our PN spacetime, $\hat{g}_{\mu \nu }$, has this form to within the accuracy of our PN procedure.

The initial state of the simulation consists of matter in axisymmetric hydrostatic equilibrium with a specific angular momentum profile, ℓ. We start from the discussion of De Villiers et al. (2003), which is based on Chakrabarti (1985) and other citations mentioned therein. The disk is centered about the equator of the black hole's spin; we will eventually assume that it is initially isentropic. The time independent and axisymmetric Euler–Lagrange equations reduce, essentially, to

$$\begin{equation} \frac{\partial _i h}{h} + \frac{1}{2} u^2_t \partial _i u_t^{-2} - \frac{\Omega }{1 - \ell \Omega } \partial _i \ell = 0 , \end{equation} \tag{ A3 }$$

where the angular frequency—Ω = u^ϕ/u^t—is not a simple function of the specific angular momentum—ℓ = −u_ϕ/u_t. The 4-velocity, u^μ, in our symmetry has zero components: u^r = u^θ = u_r = u_θ = 0. One can show, from the normalization condition u_μu^μ = −1, that

$$\begin{equation} u_t = {-} [{-}g^{tt} + 2 \ell \,g^{t\phi } - \ell ^2 g^{\phi \phi }]^{-1/2}, \end{equation} \tag{ A4 }$$

and

$$\begin{equation} \Omega = \frac{ g^{t\phi } - \ell \,g^{\phi \phi } }{ g^{tt} - \ell \,g^{t\phi } }. \end{equation} \tag{ A5 }$$

The solutions assume that

$$\begin{equation} \Omega = \eta \,\lambda ^{-q}, \end{equation} \tag{ A6 }$$

where η and q are yet to be determined parameters, and λ is defined by

$$\begin{eqnarray} &&\lambda ^2 = \frac{\ell }{\Omega } = \ell \,\frac{g^{tt} - \ell \,g^{t\phi }}{g^{t\phi } - \ell \,g^{\phi \phi }}. \end{eqnarray} \tag{ A7 }$$

We can eliminate Ω from this system by combining Equations (A6) and (A7) to yield a nonlinear algebraic equation for ℓ = ℓ(r, θ) in terms of the metric:

$$\begin{equation} R(\ell) = g^{t\phi } [ \ell ^2 + \lambda ^2(\ell) ] - g^{tt} \ell - g^{\phi \phi } \ell \,\lambda ^2(\ell) = 0 , \end{equation} \tag{ A8 }$$

where

$$\begin{eqnarray} &&\ell = \Omega \,\lambda ^2 = \eta \,\lambda ^{2-q}, \end{eqnarray} \tag{ A9 }$$

or

$$\begin{equation} \lambda = \left(\frac{\ell }{\eta } \right)^{1/(2-q)}. \end{equation} \tag{ A10 }$$

Also, we can show that

$$\begin{eqnarray} \Omega = \eta ^{-2/(q-2)} \ \ell ^{q/(q-2)} \equiv & k \,\ell ^\zeta, \end{eqnarray} \tag{ A11 }$$

where k = η^{−2/(q − 2)} and ζ = q/(q − 2).

Typically, one "solves" Equation (A8) by approximating λ² with its Schwarzschild value: λ² ≃ −g^tt/g^ϕϕ (De Villiers et al. 2003; Noble et al. 2009; Farris et al. 2011). For our ϕ-averaged spacetimes, the Schwarzschild approximation is not so good (see Figure 17).⁹ Therefore, we need a better solution. We can solve this equation to round off precision by using a Newton–Raphson scheme. To do so, we will need to know ∂R(ℓ)/∂ℓ:

$$\begin{equation} \frac{ \partial R }{ \partial \ell } = g^{t\phi } \left[ 2 \,\ell + \frac{ \partial \lambda ^2 }{ \partial \ell } \right] - g^{tt} - g^{\phi \phi } \left[ \lambda ^2 + \ell \,\frac{ \partial \lambda ^2 }{ \partial \ell } \right], \end{equation} \tag{ A12 }$$

where ∂λ²/∂ℓ = 2λ²/[(2 − q)ℓ]. Let us come back to Equation (A3). A solution to this equation yields our disk solution. The solution process involves integrating it from the inner disk's edge—located at r_in—to a point within the disk:

$$\begin{equation} \int ^h_{h_\mathrm{in}} \frac{d h}{h} = - \frac{1}{2} \int ^{u_t}_{u_{t \mathrm{in}}} \frac{ d (u_t)^{-2} }{ (u_t)^{-2} } + \int ^{\ell }_{\ell _\mathrm{in}} \frac{k \,\ell ^\zeta }{ 1 - k \,\ell ^{\zeta + 1}} d\ell. \end{equation} \tag{ A13 }$$

With the boundary condition, h_in = h(r_in, θ = π/2) = 1, one can solve this integral equation for h = h(r, θ):

$$\begin{equation} h = \frac{ u_{t \mathrm{in}} f(\ell _\mathrm{in}) }{ u_{t}(r,\theta) f(\ell (r,\theta)) } , \end{equation} \tag{ A14 }$$

where u_t is given by Equation (A4), ℓ is found using Newton–Raphson on Equation (A8), ℓ_in = ℓ(r_in, π/2) is a boundary value, and f(ℓ) ≡ |1 − k ℓ^{ζ + 1}|^{1/(ζ + 1)}.

Zoom In Zoom Out Reset image size

We want a distribution that resembles a torus, which has a pressure maximum at some radius r_p, and has finite extent. The parameters {ℓ_in, q, η} determine whether we get such a solution. We would like to replace one of the degrees of freedom with r_p, however, there is not a closed-form solution for r_p in terms of any of the original parameters. We know that the fluid attains the Keplerian angular momentum at the pressure maximum (r_p) as the pressure gradient must be zero there. It means that ℓ_in should be super-Keplerian at the inner edge (r_in), and ℓ_out should be sub-Keplerian at the outer edge (r_out). We therefore know that ℓ_in > ℓ_K(r_in), ℓ_out < ℓ_K(r_out), and ℓ_p = ℓ_K(r_p), where ℓ_K is the Keplerian specific angular momentum (see the next section below). However, we only have two free parameters, and now have three constraints (if we specify all ℓ_in, ℓ_out, and ℓ_p). It may be possible to change Equation (A9) to look like

$$\begin{equation} \ell = \eta \left(\lambda - \lambda _0 \right)^{2-q}, \end{equation} \tag{ A15 }$$

and then find λ₀ with this third constraint. Using these three constraints, however, does not yield a closed-form solution for q, η, and λ₀. Therefore, another Newton–Raphson procedure would be required. Hence, we relax the constraint on ℓ_out < ℓ_K(r_out), and let r_out be a result of our procedure. Using λ₀ = 0, we find that

$$\begin{eqnarray} q &=& 2 - \frac{ \log {(\ell _\mathrm{in}/\ell _p)} }{ \log {(\lambda _\mathrm{in} / \lambda _p)} }, \end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray} \eta &=& \frac{ \ell _p }{ \lambda _p^{2-q} }, \end{eqnarray} \tag{ A17 }$$

where λ_p = λ(ℓ_p, r_p) and λ_in = λ(ℓ_in, r_in) given by Equation (A7).

We follow the solution process based on one described in Chakrabarti (1985) and is the following.

• 1.  
  Choose values of (r_in, r_p, ℓ_in), and derive q and η using Equations (A16) and (A17).
• 2.  
  Calculate λ (and then ℓ) at a new (r, θ) via Newton–Raphson (using Equations (A8) and (A10)) close to the previous location, so that the old location's value can be used as a seed to the Newton iteration to successfully find a solution.
• 3.  
  Calculate u_t via Equation (A4) then h at (r, θ) via Equation (A14).

The solution process progresses throughout (r, θ) space until the boundary of the disk is found, where h(r, θ) = 1. The path we take starts at (r, θ) = (r_in, π/2), moves in increasing r along the θ = π/2 line, and we test to see if h is increasing. If h is not increasing at first, then we stop and try a different set of parameters. If h is increasing at first, we then proceed until the outer edge of the disk is reached; this radius is denoted as r_out. Then, ∀r ∈ [r_in, r_out], we start from the θ = π/2 solution and proceed backward and forward in θ along constant r to find h(r, θ).

The initial torus solution used herein is parameterized by ℓ_in = 8.743, η = 1.961, q = 1.642, r_in = 3a₀ = 60 M, r_p = 5a₀ = 100 M, r_out = 11.75a₀ = 235 M. These parameters give a total mass of 7290 in code units. Less than 1% of this mass is accreted during either RunSE or RunIn.

A.1. General Keplerian Velocity

The stationary torus solution described in Section A requires the Keplerian, or circular equatorial geodesic orbits, of the spacetime. Since we do not calculate $\hat{g}_{\mu \nu }$ in closed form, we require equations for these orbits based on a generalized metric of the form equation (A1). Here, we state the equations governing the Keplerian orbits.

Keplerian orbits in our spacetime have 4-velocity, u^μ = [u^t, 0, 0, u^ϕ] = u^t[1, 0, 0, Ω_K]. Ω_K is found from the r-component of the geodesic equation, du^r/dτ = −Γ^r_μνu^μu^ν = 0, which ultimately yields

$$\begin{eqnarray} \Omega _{K \, \pm } &=& \frac{1}{ {\Gamma ^r}_{\phi \phi } } [ - {\Gamma ^r}_{t \phi } \pm \sqrt{ ({\Gamma ^r}_{t \phi })^2 - {\Gamma ^r}_{\phi \phi } {\Gamma ^r}_{t t} }] \nonumber \\ &=& {-} \frac{1}{ \partial _r g_{\phi \phi } } [ \partial _r g_{t \phi } \pm \sqrt{ (\partial _r g_{t \phi })^2 -(\partial _r g_{\phi \phi }) (\partial _r g_{t t}) }]. \nonumber \\ \end{eqnarray} \tag{ A18 }$$

Ω_K − and Ω_K + are the prograde angular velocity and retrograde angular velocity, respectively. The 4-velocity components are found by the normalization condition:

$$\begin{eqnarray} u^t &=& [ - g_{t t} - 2 \,\Omega \,g_{t \phi } - \Omega ^2 g_{\phi \phi } ]^{-1/2}, \end{eqnarray} \tag{ A19 }$$

$$\begin{eqnarray} u^\phi &=& \Omega \,u^t. \end{eqnarray} \tag{ A20 }$$

The Keplerian specific angular momentum, ℓ_K, is found by the relation between ℓ and Ω:

$$\begin{equation} \ell = - \frac{ g_{t \phi } + \Omega \,g_{\phi \phi } }{ g_{t t} + \Omega \,g_{t \phi } }. \end{equation} \tag{ A21 }$$

We begin this section with the angular momentum equation because of its import to accretion physics. We follow the notation and derivation procedure outlined in Farris et al. (2011).

An extensive quantity, J, is the integral over the spatial volume of the time component of the its associated current, j^μ: $J = \int j^t \, \sqrt{-g} \, dV$, where dV is the spatial volume component in the spacelike hypersurface (e.g., drdθdϕ). We are interested in the azimuthal component of the momentum as that is the dominant component of the gas's and the binary's momenta. We therefore recognize that j^μ = T^μ_νϕ^ν, and ϕ^ν = (∂_ϕ)^ν = ∂x^ν/∂ϕ = [0, 0, 0, 1] in spherical coordinates, which is what we use. We wish to calculate d²J/dtdr. If J is locally conserved perfectly, ∇_μj^μ = 0. In our case it will not be conserved exactly, and exploring the radial gradient of its volume integral will help us understand how MHD stresses and the binary's gravitational torque compete over the run of the flow. This quantity is

$$\begin{eqnarray} \frac{ d }{ d r } \int (\nabla _\mu j^\mu ) \, \sqrt{-g} \, dr d\theta d\phi & = & \partial _r \int (\partial _\mu \sqrt{-g} \, j^\mu ) \, dr d\theta d\phi \nonumber \\ & = & \! \partial _r \partial _t J + \partial _r \partial _i \int j^i \, \sqrt{-g}\, dr d\theta d\phi \nonumber \\ & = & \! \partial _r \partial _t J + \partial _r \int {T^r}_\phi \,\sqrt{-g} \, d\theta d\phi , \nonumber \\ \end{eqnarray} \tag{ C1 }$$

where the last equality results from the fact that j^θ is zero on the axis, and j^ϕ(ϕ = 0) = j^ϕ(ϕ = 2π). On the other hand, we know from the stress-energy EOM—$\nabla _\mu {T^\mu }_\nu = -\mathcal {F}_\nu$—that

$$\begin{eqnarray} \frac{ d }{ d r } \int (\nabla _\mu j^\mu) \, \sqrt{-g} \, dr d\theta d\phi & = & \int (\nabla _\mu j^\mu ) \, \sqrt{-g} \, d\theta d\phi \nonumber \\ & = & \frac{ d T }{ d r } - \int \mathcal {F}_\phi \, \sqrt{-g} \, d\theta d\phi , \nonumber \\ \end{eqnarray} \tag{ C2 }$$

where the torque density, dT/dr, can be expressed as

$$\begin{eqnarray} &&\frac{ d T }{ d r } =\! \int {T^\mu }_\nu {\Gamma ^\nu }_{\mu \phi } \, \sqrt{-g} \, d\theta d\phi \ =\! \ \frac{1}{2} \int T^{\mu \nu } \partial _\phi g_{\mu \nu } \, \sqrt{-g} \, d\theta d\phi. \nonumber \\ \end{eqnarray} \tag{ C3 }$$

We remind the reader that $\mathcal {F}_\nu$ is the radiative cooling flux (see Section 3 for details).

Therefore, equating the two Equations (C1) and (C2), we have

$$\begin{eqnarray} \partial _r \partial _t J & = &\frac{ d T }{ d r } - \lbrace \mathcal {F}_\phi \rbrace - \partial _r \lbrace {T^r}_\phi \rbrace \nonumber\\ \ms & = &\frac{ d T }{ d r } - \lbrace \mathcal {F}_\phi \rbrace - \partial _r \lbrace {M^r}_\phi \rbrace - \partial _r \lbrace {R^r}_\phi \rbrace - \partial _r \lbrace {A^r}_\phi \rbrace, \qquad \end{eqnarray} \tag{ C4 }$$

where we have used here the shorthand

$$\begin{equation} \left\lbrace X \right\rbrace \equiv \int \sqrt{-g} \, X d\theta \, d\phi \ = \langle X \rangle \int \sqrt{-g} \, d\theta \, d\phi. \end{equation} \tag{ C5 }$$

Also, M^r_ϕ, R^r_ϕ, and A^r_ϕ are—respectively—the Maxwell (MHD) stress, Reynolds stress, and advected flux of angular momentum. We note that M^μ_ν = 2p_mu^μu_ν + p_mδ^μ_ν − b^μb_ν is the EM part of T^μ_ν, while (R^μ_ν + A^μ_ν) = T_H^μ_ν = ρhu^μu_ν + pδ^μ_ν is the hydrodynamic part. The Reynolds stress alone is more complicated to calculate as we have to find the perturbation from the mean flow:

$$\begin{equation} {R^r}_\phi = \rho h \, \delta u^r \, \delta u_\phi , \end{equation} \tag{ C6 }$$

where

$$\begin{equation} \delta u^\mu \equiv u^\mu - \lbrace \rho u^\mu \rbrace / \lbrace \rho \rbrace. \end{equation} \tag{ C7 }$$

We note that we include the enthalpy as it technically contributes to the stress; its contribution is insignificant, however, for our relatively cool flow. The quantities {R^μ_ν} and {A^μ_ν} are not calculated during the simulation, but found approximately from other shell-integrated quantities we do calculate; {M^μ_ν}, $\lbrace \mathcal {F}_\mu \rbrace$, and dT/dr are calculated as stated above during the run. One can easily show from Equations (C6) and (C7) that

$$\begin{equation} \lbrace {R^r}_\phi \rbrace \ = \ \lbrace \rho h \, \delta u^r \, \delta u_\phi \rbrace \ \simeq \ \lbrace {{T_H}^r}_\phi \rbrace - \lbrace {A^r}_\phi \rbrace, \end{equation} \tag{ C8 }$$

where T_H^r_ϕ is the hydrodynamic part of T^r_ϕ, and {A^r_ϕ} is calculated approximately as

$$\begin{equation} \lbrace {A^r}_\phi \rbrace \simeq \frac{\lbrace \rho \ell \rbrace \lbrace \rho h u^r \rbrace }{\lbrace \rho \rbrace }. \end{equation} \tag{ C9 }$$

Here, ℓ = −u_ϕ/u_t as it is defined in Appendix A. The approximations used to find Equations (C8–C9) include: (1) h ≃ 1 and (2) u_t ≃ −1. We have demonstrated that these assumptions are valid to the few percent level in the bound portion of the flow for our simulations described in this paper.

Torques and stresses do work on the gas, transporting angular momentum. This work can be dissipated in the disk, changing its internal energy, which is eventually radiated away in part. Here, we calculate the partitions in which the energy can move into; this calculation is nearly identical to that for d²J/dtdr in Appendix C.1. The current associated with E is e^μ = T^μ_νt^ν, where t^μ is the 4-vector along time coordinate, t^μ = [1, 0, 0, 0]. They are related by $E = \int e^t \sqrt{-g} \,dV$. Just as with j^μ, the divergence of e^μ is not exactly zero, because of the time-dependent spacetime. Using a similar analysis as before, we get

$$\begin{equation} \partial _r \partial _t E \ = \ dW/dr \ - \ \lbrace \mathcal {F}_t \rbrace \ - \ \partial _r \lbrace {T^r}_t \rbrace, \end{equation} \tag{ C10 }$$

where dW/dr = {(1/2)T^μν∂_tg_μν} is the work done by the spacetime on the matter.

The current j^μ = ρu^μ is associated with the conserved quantity M, so we have $M = \int \rho u^t \sqrt{-g} \,dV$, and

$$\begin{equation} \frac{ d M }{ d t } \ = \ - \int \rho u^r \sqrt{-g} \,d\theta d\phi \end{equation} \tag{ C11 }$$

by using a similar technique to obtain Equation (C1).