Binary Black Hole Mergers in Gaseous Environments: "Binary Bondi" and "Binary Bondi-Hoyle-Lyttleton" Accretion
aa r X i v : . [ a s t r o - ph . H E ] D ec Binary Black Hole Mergers in Gaseous Environments:“Binary Bondi” and “Binary Bondi-Hoyle-Lyttleton” Accretion
Brian D. Farris, Yuk Tung Liu, and Stuart L. Shapiro ∗ Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801
Merging supermassive black hole-black hole (BHBH) binaries produced in galaxy mergers arepromising sources of detectable gravitational waves. If such a merger takes place in a gaseous en-vironment, there is a possibility of a simultaneous detection of electromagnetic and gravitationalradiation, as the stirring, shock heating and accretion of the gas may produce variability and en-hancements in the electromagnetic flux. Such a simultaneous detection can provide a wealth ofopportunities to study gravitational physics, accretion physics, and cosmology. We investigate thisscenario by performing fully general relativistic, hydrodynamic simulations of merging, equal-mass,nonspinning BHBH binaries embedded in gas clouds. We evolve the metric using the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation with standard moving puncture gauge conditionsand handle the hydrodynamics via a high-resolution shock-capturing (HRSC) scheme. We considerboth “binary Bondi accretion” in which the binary is at rest relative to the ambient gas cloud,as well as “binary Bondi-Hoyle-Lyttleton accretion” in which the binary moves relative to the gascloud. The gas cloud is assumed to be homogeneous far from the binary and governed by a Γ-lawequation of state. We vary Γ between 4 / /
3. For each simulation, we compute the gas flow andaccretion rate and estimate the electromagnetic luminosity due to bremsstrahlung and synchrotronemission. We find evidence for significant enhancements in both the accretion rate and luminosityover values for a single black hole of the same mass as the binary. We estimate that this luminos-ity enhancement should be detectable by LSST for a 10 M ⊙ binary in a hot gas cloud of density n ∼
10 cm − and temperature T ∼ K at z = 1, reaching a maximum of L ∼ × erg s − ,with the emission peaking in the visible band. PACS numbers: 04.25.D-, 04.25.dg, 47.75.+f
I. INTRODUCTION
All bulge galaxies (including the Milky Way) arebelieved to contain a central supermassive black hole(SMBH) with a mass M between 10 M ⊙ and 10 M ⊙ [1, 2, 3]. It is also believed that galaxy mergers commonlylead to the formation of a massive black hole-black hole(BHBH) binary system in the merged remnants [4, 5]. Inthe standard picture, the BHBH binary separation de-creases first through dynamical friction due to distantstellar encounters, then through gravitational slingshotinteractions in which nearby stars are ejected at velocitiescomparable to the binary’s orbital velocity, and finallythrough the emission of gravitational radiation, leadingto coalescence of the binary [6]. These low-frequencygravitational waves will be detectable by LISA and willcontain a wealth of information about the inspiral. How-ever, it has been argued that the gaseous accretion flowwhich forms around the binary can be a source of elec-tromagnetic radiation as well [7]. This raises the excitingpossibility of a simultaneous detection of electromagneticand gravitational waves from BHBH mergers.This picture is supported by a number of observedAGNs that may be harboring BHBH binaries. VLBI ob-servations of the elliptical galaxy 0402+379 have discov- ∗ Also at Department of Astronomy & NCSA, University of Illinoisat Urbana-Champaign, Urbana, IL 61801 ered two radio sources at a projected separation of only7 pc. The existence of jets, as well as variability asso-ciated with BH activity, indicate that the system maybe a BHBH binary [8, 9, 10]. Another candidate is OJ287, a BL Lac object whose light curve shows variabil-ity with a period of ∼
12 yr. It is believed that thismay be a massive BHBH binary in which the smallerBH orbits with a period of 12 yr, penetrating the diskof the primary and giving rise to the observed variabil-ity [11, 12, 13]. The quasar SDSS 092712.65+294344is not believed to be a binary system, but rather a re-coiling BH which is the product of a binary merger.This suggestion is supported by a systematic shift of2650 km s − in its emission lines [14, 15]. Another candi-date is quasar SDSS J153636.22+044127.0, in which twobroad line emission systems are observed, separated invelocity by 3500 km s − . This observation has been in-terpreted as a BHBH binary system in which each objecthas its own emission system [16].Information from a simultaneous detection of elec-tromagnetic and gravitational waves may be useful forstudying fundamental aspects of gravitational physics.For example, in some modified gravity scenarios, thepropagation velocity for gravitons may differ from thatof photons [17, 18]. Additionally, the measurement of theluminosity distance from the gravitational wave signal atan accuracy of 1 − ∼ M ≈ M ⊙ accreting at 10% of the Eddingtonrate, the subsequent evolution of this disk, assumed op-tically thick, gives rise to a source which initially peaksin the UV band and hardens to EUV and soft X -rayemission at late times [7, 25]. Additionally, there is asudden change in the mass of the binary during the finalstage of the merger, as gravitational waves carry awaya few percent of the mass. The abrupt change in thepotential creates waves in the disk which can grow intoshocks and increase the luminosity of the disk in a char-acteristic way [26, 27, 28], giving rise to a potentiallydetectable prompt X -ray signal. Another possibility isthat the merged BH remnant experiences a recoil veloc-ity which may, in principle, be as high as several thou-sand km s − [29], although it is likely to be much lower( <
200 km / s) in most galaxy mergers [30]. This recoilingBH may “shake” or penetrate the disk, creating shockswhich could give rise to a transient signal.Various methods have been used to model plausiblesources of electromagnetic emission from BH mergers. Inone approach, the dynamics of the inspiral is ignored, fo-cusing on the effect of BH kicks and/or BH mass loss onthe hydrodynamical flow [26, 27, 28, 31, 32, 33, 34, 35].In another approach, the behavior of the gas is modeledby following the motion of collisionless “particle tracers”on geodesics [36]. Only recently has a fully relativistic,dynamical simulation of a BHBH binary merger in a hy-drodynamic setting been performed [37].In this paper we study BHBH binary mergers in thepresence of ambient gas. Modeling such systems re-quires fully general relativistic dynamical simulations, in-cluding relativistic hydrodynamics. The development ofstable algorithms to integrate Einstein’s field equationsof general relativity numerically in 3 + 1 dimensions,such as the BSSN formalism [38, 39] and the general-ized harmonic approach [40], together with the identifi-cation of suitable gauge conditions, has enabled severalpioneering simulations that demonstrate how to trackthe late inspiral, plunge and merger of a BHBH binaryin vacuum [41, 42, 43]. More refined follow-up simula-tions of these strong-field, late phases, joined onto an-alytic, post-Newtonian calculations of the early inspiralepoch [44], are now capable of producing accurate gravi-tational waveforms for merging BHBH binaries with com-panions spanning a range of mass ratios and spins. With the problem of gravitational wave emission froma vacuum BHBH binary inspiral well in hand, it is nowimportant to turn to the problem of electromagneticemission from BHBH binary coalescence in an astrophys-ically realistic gaseous environment.When gas accretes onto a SMBH in the center of agalaxy the specific angular momentum of the gas ˜ L may,in many cases, be much larger than that of a circularorbit near the horizon, ˜ L ∼ M c [45] leading to a disk-like accretion flow. For this reason, simulations to datehave focused on disk-like accretion. However, it has beenargued that for hot flows, in which the gas is near thegalaxy virial temperature, the gas is supported by pres-sure against free infall and the flow may well be describedby the spherical “Bondi” accretion model [6]. In thisso-called “cooling-flow model of quasar fueling”, inter-galactic gas during the early stages of galaxy formationis accelerated toward the center of dark matter halos andis shock-heated to the virial temperature. This gas ac-cretes then onto the growing SMBH at nearly the Edding-ton rate. For a 10 M ⊙ BH, the gas is expected to have adensity n ∼
10 cm − and temperature T ∼ K − K[6, 46].Classical “Bondi accretion” refers to spherically sym-metric, steady-state accretion of an adiabatic gas ontoa stationary star. The gas is assumed to be homoge-neous and at rest far from the star and flow adiabaticallywith a Γ-law equation of state (EOS). The problem wasfirst studied in Newtonian gravitation for a point massby Bondi [47], and later extended to accretion onto aBH in full general-relativity [48, 49, 50]. Accretion ontoa star moving with constant velocity through a cloudwhich is asymptotically uniform and at rest was firststudied qualitatively in Newtonian gravity by Hoyle andLyttleton [51] and later extended by Bondi and Hoyle[52]. The general-relativistic version for accretion ontoa single BH has been studied via numerical simulations[53, 54, 55, 56]. We refer to this scenario as the “Bondi-Hoyle-Lyttleton accretion” (BHL) problem.This paper is the first in a series of papers in whichwe will explore hydrodynamic gas flows around mergingBHBH binaries. Here we focus on hot flows with zero netangular momentum and restrict our attention to simula-tions in which the binary is placed in a gas cloud whichhas constant density and temperature at infinity. Wetreat two cases: one in which the gas is asymptoticallyat rest (binary Bondi accretion), the other where thegas is moving with constant velocity with respect to theBHBH center of mass (binary BHL accretion). Bondi andBHL accretion onto single BHs represent classic prob-lems which are well understood. We seek to use thisunderstanding as a foundation for tackling the problemof accretion onto a merging BHBH binary. The initialconditions of the gas during the merger of two supermas-sive BHs in realistic astrophysical environments is stillan open question and subject to uncertainties in cosmo-logical structure and galaxy formation scenarios, and inour understanding of the formation history and role ofmassive BHs in this process. The binary Bondi and BHLproblems provide excellent settings in which to begin arigorous probe of the binary accretion problem.The structure of the paper is as follows. In Sec. II wediscuss the unique computational challenge posed by thevast dynamical range characterizing our problem and weintroduce our technique to tackle it. In Secs. III and IV,we briefly outline the basic gravitational field and matterevolution equations and their specific implementation inour relativistic, hydrodynamic numerical scheme. Herewe also provide an overview of our initial data, gaugeconditions, and diagnostics. In Sec. IV D, we reviewcode tests which were performed to validate our numer-ical scheme. In Sec. V we review the analytic Bondisolution for a single, stationary BH and demonstrate ourability to reproduce those results. We also compare withprevious relativistic simulations of BHL accretion onto asingle BH. In Sec. VI, we summarize the results of ourbinary BHBH merger simulations. In Sec. VII we sum-marize our findings and comment on future directions.We adopt geometrized units and set G = c = 1. II. COMPUTATIONAL CHALLENGE
Simulating realistic accretion flows is extremely chal-lenging due to the enormous dynamic range in the rele-vant length scales defining the problem. One importantlength scale is set by the ADM mass M of the centralgravitating source, R = M . This is the length scale atwhich relativistic effects become significant. Another im-portant length scale is the transonic radius, R s ∼ M/a ∞ ,where a ∞ is the asymptotic sound speed in the gas cloud.This corresponds to the radius inside of which the accre-tion flow onto a stationary central mass M becomes su-personic. When the central gravitating object is movingsupersonically with velocity V ∞ relative to the gas cloud,we follow [51] and define the characteristic length scale, R BHL = MV ∞ . (1)Gas approaching the central object with impact param-eter less than ∼ R BHL will be captured. The two lengthscales R s and R BHL can be combined into a single char-acteristic radius [47], R a ≡ M ( a ∞ + V ∞ ) , (2)within which gas is bound to M and will be accreted.Another important length scale for the BHBH problem isthe binaray separation, d . When d ≪ R a , the accretionat radius r ≫ d is basically BHL flow onto a centralgravitating object of mass equal to the total ADM massof the binary. When d ≫ R a , the accretion in regionsnear each BH is again a BHL flow but for a source ofmass equal to the mass of a single BH.The typical densities and temperatures for interstellargas in a hot Bondi-like accretion flow are n ∞ ∼
10 cm − and T ∼ K respectively [6]. Assuming that V ∞ < ∼ a ∞ ,we find that R a ∼ M . With current computationalresources, it is impossible to perform a relativistic sim-ulation that follows the complete binary inspiral fromseparation d ∼ R a to d ∼ M , assuming a realistic asymp-totic gas temperature. We address this issue by first per-forming “prototype” simulations in which we artificiallyincrease the asymptotic temperatures in order to makethe accretion radius R a much smaller. We then use theseresults and scaling to perform “realistic” simulations inwhich we treat realistic temperatures, but restrict ourgrid to domains much smaller than R a . The details ofthese approaches are given in Section VI. III. BASIC EQUATIONS
Throughout this paper, we use Latin indices to denotespatial components (1-3) and Greek indices to denotespacetime components (0-3).The basic gravitational field and relativistic hydrody-namics equations are discussed in [39, 57] where theirnumerical implementation is described and code tests aresummarized. Here, we briefly sketch these equations andtheir implementation.
A. Evolution of gravitational fields
We write the spacetime metric in the standard 3 + 1form: ds = − α dt + γ ij ( dx i + β i dt )( dx j + β j dt ) . (3)where α , β i , and γ ij are the lapse, shift, and spatial met-ric, respectively. The extrinsic curvature K ij is given by( ∂ t − L β ) γ ij = − αK ij , (4)where L β is the Lie derivative with respect to β i . Weevolve γ ij and K ij using the BSSN formulation [38, 39].The fundamental variables for BSSN evolution are φ ≡
112 ln[det( γ ij )] , (5)˜ γ ij ≡ e − φ γ ij , (6) K ≡ γ ij K ij , (7)˜ A ij ≡ e − φ ( K ij − γ ij K ) , (8)˜Γ i ≡ − ˜ γ ij,j . (9)The evolution and constraint equations for these fieldsare summarized in [38, 39]. We assume in this paperthat the mass of the gas is negligible compared to themass of the BHs, thus we do not include matter sourceterms in our metric evolution equations.Adding Kreiss-Oliger dissipation to the BSSN evolu-tion equations outside the BH can reduce high-frequencynumerical noise associated with AMR refinement inter-faces [43]. We use this technique here and have foundit useful in reducing Hamiltonian and momentum con-straint violations.We adopt the standard puncture gauge conditions:an advective “1+log” slicing condition for lapse and a“Gamma-freezing” condition for shift [58]. Thus we have ∂ α = − αK, (10) ∂ β i = (3 / B i , (11) ∂ B i = ∂ ˜Γ i − ηB i , (12)where ∂ = ∂ t − β j ∂ j . The η parameter is set to 0 . /M in all simulations. B. Evolution of hydrodynamic fields
The fundamental matter variables are the rest-massdensity ρ ≡ nm B , specific internal energy ǫ , pressure P , and four-velocity u µ . The stress-energy tensor for anideal gas is given by T µν = ρ hu µ u ν + P g µν , where h = 1 + ǫ + P/ρ is the specific enthalpy. Weevolve the “conservative” variables ρ ∗ , ˜ S i , and ˜ τ . Theyare defined as ρ ∗ = −√ γρ n µ u µ (13)˜ S i = √ γT µν n µ γ νi (14)˜ τ = √ γT µν n µ n ν − ρ ∗ . (15)Here n α = ( α − , − α − β i ) is the timelike unit vector nor-mal to the t = constant time slices. Evolution equationsare given by Eqs. (34), (36), and (38) of [57]. ∂ t ρ ∗ + ∂ j ( ρ ∗ v j ) = 0 , (16) ∂ t ˜ S i + ∂ j ( α √ γT ji ) = 12 α √ γT αβ ∂ i g αβ , (17) ∂ t ˜ τ + ∂ i ( α √ γT i − ρ ∗ v i ) = s , (18)where γ ≡ det( γ ij ) = e φ and v i ≡ u i /u is the fluid’s3-velocity. The energy source term s is given by s = − α √ γT µν ∇ ν n µ = α √ γ [( T β i β j + 2 T i β j + T ij ) K ij − ( T β i + T i ) ∂ i α ] . (19) C. Equation of State
To complete the system of equations, we must specifyan equation of state (EOS). While our code can handleany EOS of the form P = P ( ρ , ǫ ), we adopt a standardΓ-law EOS, P = (Γ − ρ ǫ. (20) We perform simulations with Γ = 4 /
3, 5 /
3, and 13 / / kT > ∼ m B c ∼ K ). Thechoice of Γ = 5 / kT < ∼ m e c ∼ K ). The choice of Γ = 13 / m e c ∼ K < ∼ kT < ∼ m B c ∼ K )[49]. While it is not expected that tem-peratures in a realistic accretion flow will be high enoughto make the matter behave as a Γ = 4 / / X be thefractional abundance by number of hydrogen ions. Thethermal energy density can be expressed as ρ ǫ = X (3 + 3 / nkT + (1 − X )(6 + 3 / nkT . (21)Here, each relativistic electron contributes a factor of3 nkT , and each nonrelativistic nucleus contributes a fac-tor of (3 / nkT , where n ≡ ρ /m B is the baryon numberdensity of the fluid. The pressure can be expressed as P = X nkT + (1 − X )3 nkT . (22)Combining Eq. (20), Eq. (21) and Eq. (22) we find thatΓ = (13 / X + (7 / − X ) X + (5 / − X ) . (23)Adopting cosmological abundances [59], we set X =0 .
92, which gives Γ = 1 . / . X = 1, we henceforthset Γ = 13 / / / P = 2 nkT , (24)appropriate for pure ionized hydrogen. IV. NUMERICAL METHODSA. Initial Data
For each simulation discussed in this paper, we con-sider the gas to be adiabatic (apart from shocks) withuniform density and pressure at infinity. We take the gaseither be asymptotically at rest (binary Bondi accretion),or moving uniformly (binary BHL accretion). We use the
TWOPUNCTURES code [60] to construct the initial data forthe binary BHBH metric. We choose the bare mass andmomentum of the punctures according to [61] in order toensure that the BHBH binary orbit is initially close toquasicircular.We restrict our analysis to equal-mass, nonspinningBHs in this paper. However, we note that upon merger,the remnant settles down to a spinning black hole withADM mass M f = 0 . M . Here M f denotes the finalADM mass, while M denotes the initial ADM mass of thebinary. We also measure a final spin parameter a/M f =0 .
69, in good agreement with the estimate of [62]. Forthis case, there is no recoil velocity.For all of our runs we use the analytic relativistic Bondisolution for accretion onto a stationary Schwarzschild BHas initial data. This solution is reviewed in Sec. V. Toimplement this data, we treat the binary as a single grav-itating object of mass M , where M is the ADM mass ofthe binary, and apply the analytic Bondi solution every-where outside a radius r . We follow the method of [63]and adjust the fluid parameters within the radius r ac-cording to ρ ( r ) = ρ ( r ) + dρ dr (cid:12)(cid:12)(cid:12)(cid:12) r r − r r . (25)This recipe ensures that the density and its first deriva-tive are continuous across r . Eq. (25) is, of course, notthe correct initial data for a quasistationary flow inside r , but we find that the system settles into a quasista-tionary flow within several δt ∼ r /a ( r = r ), where a is the local sound speed. For all our BHBH runs, we set r /M = 5 .
5, though we find that our results are not sen-sitive to this parameter. The fluid 4-velocity inside r isset to be radially inward, with the magnitude set accord-ing to Eq. (34). The fluid pressure is set according to apolytropic EOS, P = Kρ Γ0 , with K = P ∞ /ρ Γ0 , ∞ , where ρ , ∞ and P ∞ are the rest-mass density and pressure atinfinity, respectively.For our code tests in which we evolve only a single,stationary BH, our hydrodynamic initial data is con-structed in isotropic coordinates, which become singu-lar on the horizon at r = M/
2. For these cases, we set r /M = 0 .
6, and verify the finding of [63] that the fluidevolution quickly relaxes to the equilibrium Bondi solu-tion (see Sec. V B).
B. Evolution of metric and matter
We evolve the BSSN field equations with fourth-orderaccurate, centered, finite-difference stencils, except onshift advection terms, where we use fourth-order accurateupwind stencils. We apply Sommerfeld outgoing waveboundary conditions to all BSSN fields. Our code is em-bedded in the Cactus parallelization framework [64], and our fourth-order Runge-Kutta timestepping is managedby the
MoL (Method of Lines) thorn, with a Courant-Friedrichs-Lewy (CFL) factor set to 0.5 in all BHBHsimulations. We use the Carpet [65] infrastructure toimplement the moving-box adaptive mesh refinement. Inall AMR simulations presented here, we use second-ordertemporal prolongation, coupled with fifth-order spatialprolongation. The apparent horizon (AH) of the BH iscomputed with the
AHFinderDirect
Cactus thorn [66].We write the general relativistic hydrodynamics equa-tions in conservative form. They are evolved by ahigh-resolution shock-capturing (HRSC) technique [57]that employs the monotonized central (MC) reconstruc-tion scheme [67] coupled to the Harten, Lax, and vanLeer (HLL) approximate Riemann solver [68]. Theadopted hydrodynamic scheme is second-order accuratefor smooth flows, and first-order accurate when disconti-nuities (e.g. shocks) arise. Throughout the evolution, weimpose limits on the pressure to prevent spurious heatingand negative values of the internal energy ǫ . Specifically,we require P min ≤ P ≤ P max , where P max = 10 Kρ Γ0 and P min = Kρ Γ0 /
2. Whenever P exceeds P max or dropsbelow P min , we reset P to P max or P min , respectively.We check that this fix is applied only inside the apparenthorizon, which is causally disconnected from the rest ofthe grid.At each timestep, we need to recover the “primitivevariables” ρ , P , and v i from the “conservative” vari-ables ρ ∗ , ˜ τ , and ˜ S i . We perform the inversion as speci-fied in Eqs. (57)–(62) of [57], but with a slightly modifiedanalytic quartic solver from the GNU Scientific Librarythat outputs only the real roots. We use the same tech-nique as in [69] to ensure that the values of ˜ S i and ˜ τ yieldphysically valid primitive variables, except we reset ˜ τ to10 − ˜ τ , max (where ˜ τ , max is the maximum value of ˜ τ ini-tially) when either ˜ S i or ˜ τ is unphysical [i.e., violate oneof the inequalities (34) or (35) in [69]]. The restrictionsusually apply only to the region near the puncture insidethe horizon.For all of our “prototype” calculations, we set our outerboundary at 328 M and use 9 AMR refinement levels.The maximum resolution near each puncture is δx/M =0 . M and use 8 AMR refinement levels. Forthese cases, the highest resolution near each puncture isalso δx/M = 0 . perturbation , and neglectits influence on the hydrodynamic flow as well as anydeviation from adiabaticity. We also assume that theself-gravity of the gas can be ignored, and we do notinclude matter source terms in the evolution equationsfor the gravitational fields. C. Diagnostics
1. Rest-mass accretion rate
An important diagnostic is the rest-mass flux throughthe horizon of each BH. In order to compute this quan-tity we begin by defining a function f ( t, x, y, z ) such that f ( t, x, y, z ) = 0 on the 3d hypersurface which correspondsto the world tube of a BH horizon. Thus, we can write f = p ( x − x h ( t )) + ( y − y h ( t )) + ( z − z h ( t )) − R ( t, θ, φ ) . (26)Here { x h ( t ) , y h ( t ) , z h ( t ) } represents the coordinate posi-tion of the BH center at coordinate time t and R ( t, θ, φ )represents the coordinate distance from the BH center tohorizon in the ( θ, φ ) direction. Here θ and φ are sphericalpolar coordinates with the origin set at the BH center.The rest-mass flux is then given by,˙ M = − Z α √ γρ u µ ∂ µ f Jdθdφ , (27)where J is the Jacobian J = (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( t, f, θ, φ ) ∂ ( t, x, y, z ) (cid:12)(cid:12)(cid:12)(cid:12) − = (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( f, θ, φ ) ∂ ( x, y, z ) (cid:12)(cid:12)(cid:12)(cid:12) − , (28)and where the integral is over each BH (apparent) hori-zon.The subtlety in this calculation arises from the factthat the position and shape of each horizon changes intime for a dynamical spacetime, so we cannot neglect thetime derivatives of f in Eq. (27). For a full derivation ofEq. (27), see Appendix A.We note that the ˙ M defined above is inherently gaugedependent. This is true even in a spacetime with atimelike Killing vector field. In our adopted gauge thesum of mass fluxes across the BH horizons is equal tothe total mass flux through a large sphere measuredby a distant observer in the cases in which the flow is(quasi)stationary. Quasi-stationary flow is realized dur-ing the BHBH inspiral phase, as well as after the binarymerger once the remnant settles down. See Appendix Afor a detail discussion on various issues involving ˙ M andgauge choices.
2. Luminosity
In order to compute the observed electromagnetic ra-diation exactly, it would be necessary to employ a fullyrelativistic, radiative-transfer integrator suitable for adynamical spacetime in 3 + 1 dimensions. While ad-vancements toward constructing such a scheme have beenmade in various approximations (e.g. [70, 71, 72, 73, 74,75, 76, 77, 78, 79, 80, 81, 82]), this is an extremely chal-lenging task which has not yet been fully accomplished. In this paper, we instead adopt a crude method for esti-mating luminosities from our accreting gas, assuming it isan optically thin medium. Results should be interpretedas order of magnitude estimates only. Moreover, since weare primarily interested in enhancements and variationsof the luminosity which may correlate with detectablegravitational wave signals, even crude estimates can pro-vide valuable information.We assume that the radiation propagates through anoptically thin gas and we neglect the roles of radiationpressure and radiative cooling on the hydrodynamic evo-lution. Our solution is therefore valid only when theluminosity is below the Eddington limit and the radia-tive flux is below the thermal gas flux. To calibrate ourcrude approximation scheme we compare our results tothe more exact calculation of optically thin emission fromaccretion onto a single, static Schwarzschild BH at restin a gas cloud [49]. We first review this solution below,then describe the approximations we make in this paperto simplify the calculation.The total luminosity L ∞ received at infinity from gasaccreting steadily onto a stationary Schwarzschild BH isgiven by [49], and accounts for Doppler and gravitationalredshifts as well as photon capture from the BH. Theresult is L ∞ = Z ∞ L ν (cid:18) dν dν (cid:19) dν . (29)Here, dν /dν accounts for Doppler and gravitational red-shifts, and is given by dν dν = (1 − v ) / − v cos Θ ′ (1 − M/r ) / , (30)where v is the proper velocity of the fluid as measuredby a stationary observer, the specific luminosity L ν isgiven by L ν = 8 π Z ∞ M r dr Z cos Θ ∗ − j ν (1 − v )(1 − v cos Θ ′ ) d (cosΘ ′ ) , (31)and where | cosΘ ∗ | = " (cid:18) Mr (cid:19) (cid:18) Mr − (cid:19) + 1 / . (32)Here the emissivity j ν is the specific emissivity (the en-ergy emitted isotropically per time per volume per fre-quency interval) in the comoving frame of the fluid.For dynamical systems containing inspiraling BHBHbinaries, simple analytic expressions for the luminositysimilar to Eq. (31) do not exist. In our simulations weemploy a crude approximation for the luminosity. Wesimply compute L ( t ) ≈ Z V j dV , (33)where dV is the proper 3-volume element in the fluid.We exclude from our integration points less than ξr h from each horizon. Here r h is the apparent horizon ra-dius, ξ is a constant chosen so that for a single, isolatedpuncture, r = ξr h corresponds to the surface of constantSchwarzschild radius ˆ r = 3 M . We have made this choicebecause within this radius, 50% of the emitted radia-tion by a stationary source will be captured (see Eq. (24)in [49]). For realistic temperatures (see Sec. VI C), wefind that the contribution to the luminosity is dominatedby emission from near the binary. Our measured lumi-nosity is thus insensive to the outer limit of integrationand we chose to integrate to the outer boundary of ourintegration domain for definiteness. For our high tem-perature prototype runs, however the gas has a nonneg-ligible contribution to the luminosity even far from thebinary, thus the measurement is sensitive to the outerlimit of the integration. We have chosen to integrate to r out = 44 . M , which corresponds to r out ≈ R a for ourprototype BHBH Bondi runs. We have verified that fora single, stationary puncture, our luminosity estimate iswithin a factor of ∼ P = βP M ≡ βB / (8 π ). Wethus assume that the magnetic pressure is some fraction1 /β of its equipartition value. Simulations of magne-tized accretion flows have demonstrated that the mag-netic fields do not typically reach their full equipartitionvalue [83]. We have chosen β = 10 to account for this. Weprovide further details of the bremsstrahlung and syn-chrotron emissivities which we use in Appendix B. D. Code Tests
Our general relativistic hydrodynamic code hasbeen thoroughly tested by passing a robust suite oftests. These tests include maintaining stable rotatingstars in stationary equilibrium, reproducing the exactOppenheimer-Snyder solution for collapse to a BH, andreproducing analytic solutions for relativistic shocks andspherical Bondi accretion onto isolated BHs [57]. Ourcode has also been used to simulate the collapse ofvery massive, rotating stars to black holes [84]; mergingBHBH binaries [85], BHNS binaries [69, 86], and rela-tivistic hydrodynamic matter in the presence of punctureblack holes [63]. Recently, our code has been generalizedto incorporate (optically thick) radiation transport andits feedback on fluids in dynamical spacetimes [82]. All of the above tests and simulations were performedon grids with uniform spacing. In some of the simula-tions, we utilized the multiple-transition fisheye transfor-mation [87] so that a uniform computational grid spacingcorresponds to physical coordinates with spatially vary-ing resolution. Recently, we have modified our code sothat we can use the moving-box AMR infrastructure pro-vided by Carpet [65]. To test our new code, we have per-formed shock-tube tests and 3+1 simulations of lineargravitational waves, single stationary and boosted punc-ture BHs, puncture BHBH binaries, and rapidly and dif-ferentially rotating relativistic stars. Our AMR code hasalso been used to perform simulations of BHNS mergers[86].All of our 3+1 AMR code tests were performed as-suming equatorial symmetry (i.e., symmetry about the z = 0 orbital plane), which we assume in all evolutionspresented in this paper. We have also checked that ourAMR code is able to accurately reproduce the analyticBondi solution, and we have shown good agreement withthe results of [53] for BHL accretion. Results of the lattertwo tests are summarized in Sec. V B. V. ACCRETION ONTO A SINGLE BHA. Analytic Relativistic Bondi Solution
Steady state, adiabatic, spherically symmetric accre-tion onto point masses was first considered by Bondifor Newtonian gravitation [47]. This work was later ex-tended to handle accretion onto BHs in general relativ-ity [48, 49, 88, 89]. A thorough discussion of the relativis-tic solution may be found in Appendix G of [50]. Here webriefly summarize this solution. We consider sphericallysymmetric, steady-state accretion onto a SchwarzschildBH of mass M . We assume that the BH is at rest in aninfinite gas cloud which has rest-mass density ρ ∞ , pres-sure P ∞ , and fluid 4-velocity u i = 0 at infinity. We as-sume that the gas is adiabatic with an adiabatic index Γ.We can solve the equations of relativistic hydrodynamicsto derive an exact solution for the accretion flow.
1. Key Equations
For simplicity, we will derive the solution inSchwarzschild coordinates, then transform to isotropiccoordinates. We begin by recasting the steady-state rel-ativistic continuity and Euler equations into an integralform, 4 πρ ˆ u ˆ r ≡ ˙ M = const . (34) h (cid:18) − M ˆ r + ˆ u (cid:19) ≡ h ∞ = const . (35)Here, ˆ r and ˆ u denote the radius and radial fluid 4-velocityin Schwarzschild coordinates.For adiabatic flow, the Γ-law EOS implies the poly-tropic relation P = Kρ Γ0 , with K =const. Thus P = Kρ γ = (Γ − ρ h − . (36)For an adiabatic gas, the speed of sound is, a ≡ h / (cid:18) dPdρ (cid:19) / = (cid:18) Γ Pρ h (cid:19) / . (37)This leads to the following relations between the en-thalpy, the speed of sound, and the temperature a = (Γ − h − h (38) kT ≡ P n = m B −
1Γ ( h − . (39)Any solution satisfying Eqs. (34) and (35) that maintainsthe causality constraint a < u s and the speed of sound aregiven by ˆ u s = M R s , (40) a s = M/ R s − M/ R s . (41)The transonic radius sets the length scale at which thefluid parameters begin to deviate from their asymptoticvalues during inward flow. We derive the transonic radiusby plugging Eq. (41) into Eq. (G.30) from Appendix Gof [50]. This leads to a polynomial equation which wesolve using a numerical root solver. In the Newtonianlimit, in which a ∞ ≪ M/R s ≪
1, we recover theexpressions, R s = − Ma ∞ ≤ Γ < / Ma ∞ Γ = 5 / . (42)Knowing the transonic radius, the accretion rate isgiven by ˙ M = 4 πρ s ˆ u s ˆ R s = 4 πλM ρ ∞ a − ∞ , (43)where λ = λ (Γ , a ∞ ) is a parameter of order unity, whichin the Newtonian limit is given by λ = (cid:18) (cid:19) (Γ+1) / − (cid:18) − (cid:19) − (5 − / − . (44)In the relativistic domain, we must solve a cubic equationfor to find λ (see [50] for details). Given values for ρ , ∞ and T ∞ , we can determine ˙ M and h ∞ using Eqs. (38),(39), and (43). We can then use Eqs. (34) and (35), alongwith the EOS, Eq. (36), to obtain the complete analyticsolution.
2. Transformation to isotropic coordinates
In order to provide initial data for our simulations, wetransform the solution to isotropic coordinates, the coor-dinate system adopted for our metric (puncture) initialdata. The relation between Schwarzschild radius ˆ r andisotropic radius r outside the horizon is given byˆ r = r (cid:18) M r (cid:19) (45) r = ˆ r − M + p ˆ r (ˆ r − M )2 . (46)Since the flow is purely radial, we have u = − u r = − ˆ u ˆ r drd ˆ r = ˆ u (cid:18) d ˆ rdr (cid:19) − = ˆ u (1 − M/ r )(1 + M/ r ) , (47)where u r is the radial component of the 4-velocity. Thetime component of the 4-velocity, u , remains unchangedby the transformation. We find the 3-velocity v = v r = u r /u from the normalization condition u α u α = −
1, or − α ( u ) + e φ ( u ) v = − r > r . Inside r we set theinitial data according to Eq. (25). We choose r = 0 . M ,which is outside the horizon, to avoid the initial coordi-nate singularity of u at the horizon. There is no singu-larity during the evolution as the adopted 1+log slicingcondition is horizon penetrating.
3. Moving Bondi solutions
We also wish to consider the problem of adiabatic,steady-state accretion onto a BH which moves with con-stant velocity V ∞ through a gas which is uniform atinfinity. This problem was first studied in Newtoniangravitation in [51] and [52], and is often referred to as“Bondi-Hoyle-Lyttleton (BHL) accretion”. The problemhas also been studied numerically in Newtonian theory by[54, 55], and for BHs in full general relativity by [53, 56].For cases in which the gas is not at rest asymptotically,we construct initial data following the method describedin [63], in which we “boost” the stationary Bondi solu-tion by taking the initial density at any coordinate pointto be approximately the stationary isotropic Bondi solu-tion, and computing the initial velocity field employing aLorentz transformation for the velocities. Denoting thestationary Bondi flow velocity by v i , the moving Bondiflow by v i ′ , and the “boost” velocity by V ∞ , we set v x ′ c l = V ∞ + v x c l V ∞ v x c l (49) FIG. 1: Semi-analytic profiles of density, temperature andfluid 4-velocity for spherical, adiabatic Bondi accretion ontoa Schwarzschild BH. The solid line gives ρ /ρ , ∞ , the dot-ted line shows T /T ∞ , and the dashed line shows the fluid4-velocity in Schwarzschild coordinates, ˆ u . The adiabatic in-dex is given by Γ = Γ ∗ according to Eq. (53). The asymptotictemperature is T ∞ = 10 K . The horizon is labeled R h . v y ′ c l = v y c l γ b (1 + V ∞ v x c l ) (50) v z ′ c l = v z c l γ b (1 + V ∞ v x c l ) , (51)where γ b ≡ (1 − V ∞ ) − / and c l ≡ α/ψ . Here α is thelapse, and ψ ≡ e φ , with φ defined in Eq. (5). For anyhydrodynamic quantity g ( x, y, z ), we compute g ( x ′ , y ′ , z ′ ) = g ( γ b x, y, z ) , (52)which accounts for the coordinate transformation be-tween the x α and x α ′ frame on the t ′ = 0 timeslice. Theabove method of computing initial data is only strictlyvalid asymptotically. This is because we use the Bowen-York initial data for the moving BH metric, which are notobtained by boosting the stationary BH metric. How-ever, we find that deviations propagate away quickly,leaving a stationary flow. In Sec. V B we describe a sim-ulation of Bondi accretion which we perform in a movingreference frame. In this case, we use the same techniquedescribed above to construct hydrodynamic initial data.
4. Effective adiabatic index
In some cases, we consider a gas cloud which the elec-trons have nonrelativistic temperatures at infinity, butachieve relativistic temperatures near the binary. In thiscase we follow [49] and define an “effective adiabatic in-dex” Γ ∗ according toΓ ∗ = (cid:26) / kT /m e c ≤ / / kT /m e c > / , (53)Thereby replacing the actual continuous transition by asimpler discrete transition. This transition occurs at aSchwarzschild radius [50] r ∗ M ≈ m p m e ≈ . (54)We solve for r ∗ using Eqs. (34), (35), (39), and the factthat kT ( r ∗ ) = (2 / m e c . Using the continuity of ρ and P at r = r ∗ , we obtain the full equilibrium flowsolution. In practice, r ∗ is outside the outer boundaryof the computational grid in our simulations, so we stillimplement a constant Γ. However, taking into accountthis transition alters our initial data significantly, sincethe outer Γ = 5 / a < ∼ u up to r = r ∗ ,increasing the gas temperature near the black hole whencompared to gas in which Γ = 13 / B. Relativistic Bondi test
We test our code’s ability to accurately simulate hy-drodynamic accretion onto a moving puncture. This isparticularly important in an AMR code in which mattercrosses moving refinement zone boundaries. To test ourcode, we simulated spherical Bondi accretion in a framein which both the BH and gas cloud are moving with thesame velocity. We consider a BH moving with velocity V BH ≡ P x /M = 0 .
37, initially situated at the origin.Here P x is the momentum of the BH as measured by astationary coordinate observer at infinity. We consider agas with adiabatic index Γ = 13 / a ∞ = 0 .
022 also moving with velocity v = V BH atinfinity. With this choice the transonic radius is givenby areal radius R a = 45 . M in the comoving frame. Weevolve for a duration t = 713 . M , which is approximatelyequal to two free-fall times at the accretion radius. Bythis time, the BH has moved to a coordinate location of x = 247 M .To assess the agreement with the analytic solution, wecompare invariant quantities [63]. Two such invariantsare the fluid rest-mass density, ρ , and the the rate ofchange of the fluid rest-mass density, as measured by acomoving observer, dρ /dτ . We choose a set of 6 arealradii in the comoving frame. At each radius, we com-pute ρ and dρ /dτ analytically. By construction, con-tours of ρ and dρ /dτ are spherical and coincide in0 FIG. 2: Snapshots showing spherical Bondi accretion onto asingle BH simulated in a frame in which the BH and gas movewith velocity V BH = 0 .
37. Contours of constant density ρ (solid black lines) and constant dρ /dτ (dotted red lines) areshown. Their overlap indicates agreement with the analyticBondi solution. The adiabatic index is set to Γ = 13 / a ∞ = 0 . this frame. However, because ρ and dρ /dτ are bothcoordinate-independent quantities, their contours mustcontinue to coincide in any coordinate frame. We use thisfact to check that our numerical simulation, performedin a frame in which the BH and gas are moving, matchesthe analytic solution (see Fig. 2). C. Relativistic Bondi-Hoyle-Lyttleton test
We will simulate cases of accretion onto binaries inwhich the BHBH center of mass moves relative to theasymptotic gas cloud. Here we verify that we can accu-rately simulate “BHL” accretion of a single puncture BHmoving at constant velocity through a gas cloud. Whilethere is no analytic solution for this problem, we com-pare our findings with results of previous work. We con-sider a test case in which a single BH puncture is placedin a cloud with asymptotic sound speed a ∞ = 0 . /
3, with the puncture moving supersonically withspeed V ∞ = 0 .
25 relative to the gas cloud. We performour simulation in a frame in which the puncture is at restand the asymptotic fluid velocity is set to V ∞ . Figure 3shows a snapshot of the stationary-state density and ve-locity profile from our simulation. We measure ˙ M and FIG. 3: Snapshots of density contours for the BHL ac-cretion code test. The adiabatic index is Γ = 4 /
3, thesound speed is a ∞ = 0 .
1. Density contours are chosen at ρ = ρ , ∞ . j ( j = 1 , , .... compare with the results of [53]. In order to facilitatethis comparison, we define a canonical unit of rest-massaccretion flux to be˙ M can = 4 πλM ρ , ∞ ( V ∞ + a ∞ ) / , (55)We find ˙ M / ˙ M can = 2 .
6, which is slightly smaller thanthe value 3.0 reported in [53]. The most likely sources ofthe small discrepancy are the outer boundary condition(we use the stationary spherical Bondi solution with aLorentz boost, whereas [53] use constant asymptotic val-ues), and the fact that [53] impose an approximate innerboundary condition outside the horizon. Other differ-ences between these two codes are that we use a 3+1 codewith AMR, whereas [53] employ an axisymmetric code.Also, our outer boundary is placed at r max /M = 820for this test, whereas [53] place the outer boundary at r max /M = 140. VI. ACCRETION ONTO A BHBH BINARY
As discussed in Sec. II, the transonic radius for Bondiaccretion with a realistic asymptotic temperature of10 K is R a ∼ M . It is beyond the capability of cur-rent 3+1 GR simulations to evolve the binary from initialseparation d > R a all the way to merger, while resolvinga BH horizon of size ∼ M . The range of length scales is1 TABLE I: Parameters for BHBH simulationsrun run type case V ∞ V ∞ /a ∞ Γ T ∞ ( K ) a ∞ † R a /M †† d/M PA1 prototype Bondi 0 . . / . × .
148 22.7 40 , , , / . × PA3 5 / . × PB1 BHL 0 . . × .
148 15.6 40 , . × PB3 5/3 7 . × PC1 BHL 0 . . × .
148 2.74 40 , . × PC3 5/3 7 . × RA1 realistic Bondi 0 . ∗ × . × − . × / † R a = ( M/ / ( a ∞ + V ∞ ) / is the accretion radius for a single BH of mass M/ †† d = initial binary separation; simulations beginning at 10 M proceed to merger.TABLE II: Parameters for BHBH simulationsrun flow characteristics § emission characteristics §§ n max /n ∞ T max /T ∞ K max /K ∞ † ˙ M max / ˙ M a †† L maxff /L hν maxff †† L maxsyn /L hν maxsyn PA1 631 18.3 1.2 3.7PA2 1743 14.3 1.5 3.8PA3 164 30.1 1.6 2.6PB1 593 18.9 1.2 4.0PB2 1638 15.2 1.7 3.3 Not Physically RelevantPB3 153 30.8 1.9 2.5PC1 99.2 28.3 6.0 0.7PC2 158 32.2 14.0 0.7PC3 46.6 37.7 6.5 0.8RA1 1 . × . × × β − / (1 + z ) n / T − / β − / eVRA2 3 . × . × × β − / (1 + z ) n / T − / β − / eV † ˙ M a c = 4 . × λ / n T − / M erg s − is the accretion rate onto a single BH of mass M/ †† L = 10 n T − M erg s − § “max” label refers to the characteristic value at the moment of maximum luminosity. §§ “max” label refers to the maximum values of the luminosities during the merger and thecharacteristic frequencies at these times. too large and the total coalescence time too long for sucha task. We approach this issue by performing two typesof simulations. We perform simulations of binaries merg-ing in “realistic” gas clouds with asymptotic temperature T ∞ = 10 K , following only the last phase of the mergerin which the binary separation satisfies d ≪ R a . Ourfocus here is on identifying observable electromagneticsignals generated by the time-dependent shock heatingcaused by the binary motion. We also perform “proto-type” calculations with artificially high temperatures andsound speeds in order to study how the accretion flowchanges as the binary transitions between the followingregimes during the inspiral:1. “widely separated regime”, in which d > R a
2. “moderately separated regime”, in which d ≈ R a
3. “closely separated regime”, in which d < R a We set the asymptotic temperature T ∞ ∼ K in the“prototype” calculations to decrease the accretion radius R a to a value closer to d so that we can explore all threeregimes by combining the results of several numericalsimulations. All of our BHBH simulations are summa-rized in Table I. Important results are summarized inTable II. The initial BHs in the binary are all equal-mass and nonspinning. A. Scaling
For a given asymptotic gas temperature T ∞ or soundspeed a ∞ , our solution for binary Bondi flow can bescaled to arbitrary density n ∞ ≡ ρ , ∞ /m B (neglectingself-gravity of the gas) and black hole mass M . Hencea single simulation with an arbitrary n ∞ and M suf-2fices to determine the solution for any other n ∞ or M .Thus, for example, the accretion rate is proportional to n ∞ and M (see, e.g. Eq. (43)), while the 4-velocity u α as a function of coordinates x α /M are independent of n ∞ and M . If the asymptotic sound speed is sufficientlylow that a ∞ ≪
1, then the solution can also be scaledto arbitrary a ∞ or T ∞ . However, once a ∞ approachesthe speed of light and the transonic radius approachesthe horizon, scaling with a ∞ or T ∞ breaks down. Thisbehavior is already evident from the relativistic Bondisolution for accretion onto a single BH.The emergent electromagnetic luminosity and radia-tion spectrum also exhibit simple scaling. The formof the scaling relations depend on the temperature anddensity dependence of the adopted emissivities and willbe described later in Sec. VI C when we treat realisticasymptotic temperatures and sound speeds. B. Prototype Cases
1. Binary Bondi Accretion
The “prototype” calculations provide valuable qualita-tive insight into the evolution of the accretion flow as thebinary separation decreases, passing through the threeregimes identified above. For the case of “binary Bondiaccretion”, in which the gas at infinity is at rest relativeto the binary center of mass, we obtain a sequence of“snapshots” for different binary separations. Each snap-shot is generated by evolving the binary long enough toallow the gas to relax to a quasistationary flow, but notlong enough for the separation d to change significantly.By sandwiching together snapshots we can trace the evo-lution in accretion rate ˙ M and luminosity L as the bi-nary transitions from regime 1 to regime 2. The finaltransition from regime 2 to regime 3 is captured by asingle simulation which follows the binary through inspi-ral and merger. We perform these simulations for casesPA1 (Γ = 13 / / / M can be understood as follows. De-fine ˙ M a to be the accretion rate onto a single, isolatedBH of mass M/
2. Then it follows that the total accre-tion rate onto two infinitely separated BHs of mass M/ M = 2 ˙ M a . However, late in the inspi-ral, when the separation satisfies d ≪ R a , the binary canbe treated as a single gravitating object. From Eq. (43)we see the accretion rate is proportional to M , so weexpect that the accretion rate will approach ˙ M = 4 ˙ M a .During the final stage of the inspiral, mass-energy is ra-diated away in the form of gravitational waves. Thus weexpect that the final post-merger accretion rate will beapproximately given by ˙ M = 4 ˙ M a (1 − δM/M ) . In oursimulations, we find δM/M ≈ .
05, consistent with the
FIG. 4: Time evolution of ˙ M and δL for binary Bondi ac-cretion for the prototype case with Γ = 13 /
9. Time t/M ismeasured relative to the time at which the merger occurs.˙ M a and δL a are the accretion rate and luminosity enhance-ment over the background for a single isolated black hole withmass equal to the initial ADM mass of the binary. The left-hand box shows values from “snapshots” in regimes 1 and2. The right-hand box shows results for the final inspiraland merger from an initial separation of d = 10 M (regime 3).Solid lines denote numerical data, dashed lines denote extrap-olated data. Solid dots correspond to profile plots highlightedin Fig. 7, while open circles show expected values at t = ±∞ .The asymptotic sound speed is set at a ∞ = 0 . R a /M = 22 . values reported in [62]. We observe the expected behav-ior for both our PA1 (see Fig. 4) and PA2 (see Fig. 5)runs. For PA3, the accretion rate never reaches 4 ˙ M a be-fore the merger. We attribute this to the fact that forΓ = 5 /
3, gas is more efficiently heated, allowing P ∼ ρ v and a ∼ v , even in the absence of shocks. When shocksdo form, the flow is much more easily disrupted as thekinetic energy of the flow does not dominate the thermalenergy, contrary to cases PA1 and PA2. Thus, some mat-ter is swept away from the vicinity of the binary, causing3 FIG. 5: Same as Fig. 4, but for Γ = 4 / a lowering of the accretion rate (see Fig. 6). We note,however, that after the merger when the shocks have achance to dissipate, the accretion rate does settle to itsexpected value, taking mass loss into account.We also plot the luminosity enhancement due tobremsstrahlung and synchrotron emission in Figs. 4–6.Because the high-temperature homogeneous backgroundgas in our prototype simulations has an intrinsic, non-negligible emissivity, we subtract it from the total lumi-nosity measured. We define δL ≡ L − L bg , where L bg is the background luminosity which would be present inour computational domain for a homogeneous gas cloudof density ρ , ∞ and temperature T ∞ , with no BH present.We normalize δL by δL a , which we define as the lumi-nosity above background which would be present for asingle, isolated BH of mass M/
2. If we take the limitin which the binary separation d → ∞ , we expect that δL/δL a →
2. This limiting value is indicated by theopen circle at t/M = −∞ in Figs. 4–6. We also calculatethe expected value of δL/δL a for a single BH of mass M − δM and plot it for reference with an open circle at t/M = + ∞ . We see that for each Γ the luminosity en- FIG. 6: Same as Fig. 4, but for Γ = 5 / hancement increases by several orders of magnitude overthe course of the inspiral. While the numerical value ofthis variation is not physically meaningful due to the un-realistic temperatures used in these “prototype” calcula-tions, this behavior provides strong qualitative evidenceof a significant enhancement in luminosity that can beexpected to accompany such an inspiral. Mergers in re-alistic clouds, yielding realistic luminosities whose valuesare physically meaningful, will be treated in Sec. VI C.Figures 7 and 8 show snapshots of density and tem-perature contours for cases PA1 and PA3. We do notshow snapshots for the PA2 case because they look verysimilar to the PA1 case. We can see that in the earlyphases of the inspiral, the accretion flow resembles twoindependent spherical Bondi flows. As the separation de-creases and becomes comparable to the transonic radius,the orbital velocity of each BH becomes comparable tothe sound speed. Within R a shocks begin to form, whichgrow in strength until the merger. It is the heating fromthese shocks which contributes to the dramatic increasein the luminosity observed. Note that the final accre-tion flow near the BH is not spherically symmetric due4 FIG. 7: Snapshots of rest-mass density ρ and temperature T contours in the orbital plane for the binary Bondi “prototype”case with Γ = 13 /
9. First and second rows show density contours and velocity profiles, third and fourth rows show snapshotsof temperature T . Density contours are plotted at ρ = ρ , ∞ . j ( j = 1 , , ...., T = 10 . j K ( j = 1 , , ...., FIG. 8: Same as Fig. 7, but for Γ = 5 / FIG. 9: Snapshots of
K/K ∞ contours in the orbital plane for the binary Bondi “prototype” case with Γ = 13 /
9. Contours aredrawn for
K/K ∞ = 1 + 0 . j ( j = 1 , , ...., to the spin of the merged BH. The BH spin does notsignificantly change the final accretion rate predicted byEq. (43), as this quantity is determined by gas parame-ters at r ∼ R a ≫ M , where the effect of the BH spin isnegligible. This result conforms with the findings of [89].In order to highlight the role that shock heating playsduring phases 2 and 3 of the merger, we also presentcontours of K/K ∞ for our PA1 case (see Fig. 9). Here K ≡ P/ρ Γ0 and K ∞ is the value of K at infinity. Be-cause K/K ∞ = 1 everywhere for adiabatic flow in theabsence of shocks, this quantity serves as a useful tracerfor the amount of shock heating which is taking place.The quantity K = K ( s ), where s is the gas entropy, andis constant in the absence of shocks; shock heating yields K/K ∞ > K/K ∞ increases steeply near the shock front.For each snapshot, we compute the maximum value of K/K ∞ outside the horizon. We find that K max /K ∞ ini-tially increases as the separation decreases and shocksbecome stronger, as expected. This trend terminates inthe very late stage of the merger, when d ∼ M ∼ d ISCO [90]. This is likely due to the fact that kinetic energy dis-sipated as heat is confined to a small region and is beingquickly consumed by the BHs at this stage. After themerger, the gas relaxes to laminar spherical Bondi flowand
K/K ∞ returns to unity everywhere.
2. Binary Bondi-Hoyle-Lyttleton Accretion
In order to investigate additional electromagnetic sig-natures which may be present due to the motion of thebinary relative to the cloud, we have performed a seriesof “prototype” simulations of BHL accretion onto merg-ing binaries. We consider both subsonic cases ( V ∞ /a ∞ =0 .
7) and supersonic cases ( V ∞ /a ∞ = 2 . / / / d = 40 M ) to produce snapshots thatare quasistationary in the corotating frame of the binary,as well as close separation runs ( d = 10 M ) in which weevolve the binary from inspiral to merger.For our (asymptotically) subsonic cases, the departuresfrom the binary Bondi case are subtle. Figures 10 and 11show snapshots of density and velocity field for cases PB1and PB3. Snapshots for case PB2 are similar to PB1 andhence are now shown here. At wide separation d = 40 M ,we see an asymmetry in the accretion flow for cases PB1and PB2, as shocks develop around one BH as it movesagainst the flow of the gas, but not around the other asit moves in the same direction as the flow. This phe-nomenon continues up to the merger. At this separationthe orbital Keplerian velocity is V k ≈ .
08. Thus, these7
FIG. 10: Same as Fig. 7, but for subsonic BHL accretion with Γ = 13 / V ∞ = 0 .
1. The asymptotic velocity V ∞ is in the+ˆ x direction. FIG. 11: Same as Fig. 7, but for subsonic BHL accretion with Γ = 5 / V ∞ = 0 .
1. The asymptotic velocity V ∞ is in the+ˆ x direction. FIG. 12: Time evolution of ˙ M and δL for binary BHL inspi-rals with Γ = 13 / a ∞ = 0 . V ∞ = 0 .
1. The initialbinary separation is d = 10 M and the BHs evolve to merger.˙ M a and δL a are the accretion rate and luminosity enhance-ment over the background for a single isolated black hole withmass equal to the initial ADM mass of the binary. Dashedlines in ˙ M plot represent accretion rates onto individual BHs,the solid line represents the total accretion rate. Dots rep-resent the times highlighted in the last three snapshots inFig. 10. shocks are formed when one BH moves supersonicallyagainst the flow of the gas and ( V k + V ∞ ) /a > ∼
1. The rea-son that this behavior is not seen in case PB3 is that thegas is adiabatically heated more efficiently for Γ = 5 / /
3, ( V k + V ∞ ) /a < V ∞ = 0, but all shocks dissipate. InFigs. 12–14, we once again observe an increase in accre-tion rate and luminosity over the course of the merger.We have also plotted in these figures the individual ac-cretion rates onto each BH in order to demonstrate theeffect of the binary motion relative to the wind.For our supersonic cases (PC1-PC3), the departurefrom the binary Bondi case is more dramatic. Whenthe BHs are widely separated ( d > R a ), a bow shock FIG. 13: Same as Fig. 12, but for Γ = 4 / forms around each individual BH, as seen in Figs. 15 and16. Late in the inspiral, when d < R a , these shocksmerge and form a single bow shock surrounding the bi-nary, which persists after the merger. The final flow re-sembles the solution found by [53] for steady accretiononto a moving BH with V ∞ /a ∞ = 2 .
5, although in ourcase the remnant is spinning. During the inspiral, themodulation in ˙ M and δL is more pronounced than in thesubsonic case (see Figs. 17–19). C. Realistic Binary Bondi accretion
While the high temperature prototype runs give usvaluable insight into the general nature and differentphases of the accretion flows onto inspiralling BHBH bi-nary systems in gas clouds, we require simulations withmore realistic gas temperatures in order to identify ob-servational signatures. Accordingly, we have performedsimulations of the binary Bondi problem for a gas cloudwith asymptotic density n ∞ = 10 cm − and temperature T = 10 K . This choice is consistent with the proposed“cooling flow model of quasar fueling” describing the hotinterstellar gas found in galaxies [6, 46]. As mentionedpreviously, computational limitations demand that ourouter boundary be placed inside the transonic radius for0 FIG. 14: Same as Fig. 12, but for Γ = 5 / these simulations. As a result, we focus only on the finalphase of the inspiral and merger, in which d ≪ R a . Weuse both a Γ = Γ ∗ (case RA1) and Γ = 5 / / ∗ case sincethe Γ = 5 / / V ∞ = 0), andpostpone a study of binary BHL accretion in a realisticgas cloud for a later analysis. As in the prototype calcu-lations, we find evidence for a strong enhancement in theluminosity due to shock heating of the gas (see Fig. 20and Fig. 21). For these runs, we plot the luminosities incgs units.We note that for case RA1, strong shocks form, but areconfined to the immediate vicinity of the binary. This isbecause the inward gas flow onto the binary is stronglysupersonic, making it difficult for shocks to propagateoutward (see Fig. 22). For case RA2, on the other hand,the radial component of the 4-velocity u ∼ a everywhere,making it much easier for shocks to spread outward, asseen in Fig. 23. D. Scaling and Detectability
All of our quoted results for the accretion rate ˙ M are normalized by the value for a single black hole ofmass M/ M a c = 4 . × λ / n T − / M erg s − , (56)where we define λ / ≡ λ (Γ , a ∞ ) /λ (5 / , a ∞ ), n ≡ n ∞ / − , T ≡ T ∞ / K , and M ≡ M/ M ⊙ .Recall that in the Newtonian limit, when R a ≫ M , λ depends only on Γ.It is also possible to derive simple scaling relations forthe luminosities. In each of our simulations, most of theelectromagnetic luminosity is generated near the horizonsof the BHs, as the temperature and density both risesharply when approaching the horizon, and achieve theirmaximum values there. By examining Eqs. (B1)–(B5)(ignoring the logarithmic terms), and Eq. (B18), we seethat in the high temperature limit ( kT > m e c ), whichapplies near the horizon in all of our simulations, thebremsstrahlung and synchrotron emissivities depend onthe density and temperature according to, q ff ∝ n h T h (57) q syn ∝ n h T h β − (58)Here n h and T h refer to the density and temperature atthe horizon, respectively. Integrating Eqs. (57) and (58),we find L ff ≈ Z dV q ff ∝ n h T h M (59) L syn ≈ Z dV q syn ∝ n h T h β − M , (60)where we have assumed that the radiation is generatednear the horizon and hence the radiative volume scalesas M . To estimate n h , we use the fact that r h ≪ R a ,and so the fluid 4-velocity is approximately given by itsfree-fall value ˆ u ≈ (cid:18) Mr h (cid:19) / . (61)Substituting this into Eq. (34) and Eq. (43), we find that n h n ∞ ∼ (cid:18) R a r h (cid:19) / ∝ (cid:18) kT ∞ m B c (cid:19) − / (62)Here we have used a ≈ Γ P/ρ = 2Γ kT /m B c . Thisscaling should remain reasonably accurate, even in thepresence of strong shocks, as the shocks cause densityenhancements of < ∼ (Γ + 1) / (Γ −
1) (see Sec. 89 of [91]),which is of order unity.The temperature at the horizon T h , on the other hand,will be strongly affected by the presence of shocks. The1 FIG. 15: Same as Fig. 7, but for subsonic BHL accretion with Γ = 13 / V ∞ = 0 .
4. The asymptotic velocity V ∞ is in the+ˆ x direction. FIG. 16: Same as Fig. 7, but for subsonic BHL accretion with Γ = 5 / V ∞ = 0 .
4. The asymptotic velocity V ∞ is in the+ˆ x direction. FIG. 17: Same as Fig. 12, but for Γ = 13 / V ∞ = 0 . gas in any shocked region will be heated to kT ∼ m B v .Since v < ∼ c near the horizon, shock heating guaran-tees that kT h < ∼ m B c ∼ K independant of the gastemperature at infinity T ∞ . This result has importantconsequences for the scaling of the maximum luminosi-ties. Once the shock-heated gas is accreted following themerger, T h drops below this value for all Γ < / kT h < ∼ m B c forΓ = 5 /
3. This is because for this particular EOS, anappreciable fraction of the gravitational potential energyis converted into thermal energy (both scale as r − inside R a : kT ∼ M m B /r ) [50]. For Γ = Γ ∗ , it can be shownthat the temperature at the horizon for spherical Bondiflow is approximately given by [50] kT h ≈ (cid:18) (cid:19) / (cid:18) m e m B (cid:19) / m B c = 0 . m B c (63)Thus, we find that the temperature at the horizon, T h isindependent of gas parameters at infinity for both Γ =5 / ∗ . However, for Γ = const < / T h doesexhibit scaling with T ∞ in the absence of shocks. In this FIG. 18: Same as Fig. 12, but for Γ = 4 / V ∞ = 0 . case, P = Kρ Γ0 , so we find T h T ∞ = (cid:18) n h n ∞ (cid:19) Γ − ∝ (cid:18) kT ∞ m B c (cid:19) − − / (64)We can now apply these results to see how the lu-minosity scales in different phases of the inspiral. Dur-ing both the very early pre-merger phase of the inspiral,when d > R a , and during the post-merger phase after thefluid has settled to a quasiequilibrium state, there are noshocks present. We can therefore use the above resultsto see that in these phases, L ff ∝ (cid:26) n T − M , Γ = 5 / ∗ ,n T − (3Γ+1) / M , Γ = const < / , (65) L syn ∝ (cid:26) n T − β − M , Γ = 5 / ∗ ,n T − − / β − M , Γ = const < / . (66)Here β ≡ β/
10. During the late premerger phasewhen d ∼ M ≪ R a shocks will be present and T h will nolonger depend on T ∞ for any EOS. In this case, we find L ff ∝ n T − M , L syn ∝ n T − β − M . (67)Using the above scaling relations along with the resultsof our realistic temperature simulations, we find that the4 FIG. 19: Same as Fig. 12, but for Γ = 5 / V ∞ = 0 . peak luminosity shortly before the merger of an equalmass BHBH binary in a gas cloud for case RA1 (Γ = Γ ∗ )is given by L maxff ≈ × n T − M erg s − , (68) L maxsyn ≈ × n T − β − M erg s − . (69)Similarly, the peak luminosity for case RA2 (Γ = 5 /
3) isgiven by L maxff ≈ × n T − M erg s − , (70) L maxsyn ≈ × n T − β − M erg s − . (71)At the late post-merger phase, the fluid relaxes to astationary flow. The scaling relations (65) and (66) hold.Combining these scaling relations and our simulation re-sults, we find that during the post-merger phase L ff ≈ × n T − M erg s − , (72) L syn ≈ × n T − β − M erg s − (73)for case RA1, and L ff ≈ × n T − M erg s − , (74) L syn ≈ × n T − β − M erg s − (75)for case RA2. We note that in each case, the total lumi-nosity is dominated by the synchrotron emission. FIG. 20: Plots showing time evolution of ˙ M and ˙ L . Heretime is measured relative to the time at which the mergeroccurs. Asymptotic temperature is T = 10 K . Adiabaticindex given by Γ = Γ ∗ . n ≡ n ∞ / − , T ≡ T ∞ / K , M ≡ M/ M ⊙ . In each of our calculations, we have ignored the effectsof radiative cooling on the gas dynamics. We can esti-mate the error induced by this by comparing the rate ofthermal energy transport, ˙ E th , to the luminosity. Herewe define ˙ E th = ˙ M ǫ = ˙
M Pρ (Γ − . (76)Since the luminosity is dominated by emission nearthe horizon, we are primarily concerned about the regionnear the horizon. Using Eqs. (43), (62), (63), (72) and(73), we find that for case RA1, at late times after themerger when the flow has reached equilibrium, L syn + L ff ˙ E th ∼ . n T − / β − M . (77)Thus, we see that in these regimes, it is a good approxi-mation to neglect the effects of radiative cooling for ourcanonical model parameters. At the moment of maxi-mum luminosity, shortly before merger, we find that L syn + L ff ˙ E th ∼ n T − / β − M . (78)5 FIG. 21: Same as Fig. 20, but for Γ = 5 / Thus during the final stages of the merger, the validityof our assumption of adiabatic flow begins to break downfor canonical parameters. In future work, we will addressthis by including cooling terms in our gas evolution toaccount for energy losses due to radiation.We also note that in order for us to be able to neglectradiation pressure in the momentum equation, we requirethat the luminosity be small compared to the Eddingtonluminosity. We find that for case RA1, L maxsyn + L maxff L Edd ∼ . n T − β − M , (79)which suggests that radiation pressure may begin to playa role for parameters close to our canonical choices.In calculating the luminosity, we have assumed thatthe gas is optically thin. We can verify this assumptionby estimating the optical depth. For the gas parameterschosen for this study ( n ∞ = 10 cm − and T ∞ = 10 K),we find that the dominant source of opacity is electronscattering. We estimate that the optical depth for elec-tron scattering of synchrotron photons is τ es ≈ n h σ T R ∼ − n T − / M , (80)where n h ∼ n T − / is the density at the horizon, σ T = 0 . × − cm is the Thomson scattering cross- section, and R ∼ M cm is the characteristic size ofthe emission region. Thus, our assumption of an opticallythin gas is valid.To estimate the characteristic frequencies at whichthe emission occurs we again note that the maximumemission comes from near the horizons and compute thecharacteristic frequency produced in this region. Forbremsstrahlung emission, the characteristic observed fre-quency of the emission is given by hν ∼ kT h / (1 + z )for a source at redshift z . We measure the maximumtemperature near the horizon for our case RA1, and findthat at the moment of maximum luminosity in the latepre-merger phase, hν maxff ≈
150 MeV1 + z (RA1) , (81)After the merger, in the quasistationary phase, we findthat the characteristic frequency drops to hν ff ≈
10 MeV1 + z (RA1) . (82)Following the same procedure for case RA2, we find thatat the moment of maximum luminosity in the late pre-merger phase, hν maxff ≈
230 MeV1 + z (RA2) , (83)and in the post-merger phase, the frequency drops to, hν ff ≈
50 MeV1 + z (RA2) . (84)Thus, we see that the bremsstrahlung emission willbe predominantly in γ -rays, in agreement with [49].Given the bremsstrahlung luminosity calculated above,we estimate that the flux from this emission will be ∼ − erg cm − s − for a source at z = 1. Unfortu-nately, it is unlikely that this emission is strong enoughto be detectable.We can estimate the characteristic frequency of thesynchrotron emission by noting that Eq. (B10) is maxi-mized when x M ≈ .
09. For case RA1, this correspondsto an observed frequency hν maxsyn = 1 .
091 + z ehB πm e c (cid:18) kTm e c (cid:19) (85)= 801 + z n / T − / β − / eV (RA1) (86)during the late pre-merger phase at the moment of maxi-mum luminosity. In the post-merger phase, the frequencydrops to, hν syn = 0 . z n / T − / β − / eV (RA1) , (87)For case RA2, we find the characteristic synchrotron fre-quency to be hν maxsyn = 1001 + z n / T − / β − / eV (RA2) (88)6 FIG. 22: Snapshots of rest-mass density ρ and temperature T contours for the Γ = Γ ∗ case. First and second rowsshow snapshots of density contours and velocity profiles in the orbital plane. Third and fourth rows show snapshots of T . Density is plotted according to ρ = ρ , ∞ . j ( j = 1 , , ...., T = 10 . j K ( j = 1 , , ...., FIG. 23: Same as Fig. 22 but for Γ = 5 / hν syn = 0 .
751 + z n / T − / β − / eV (RA2) (89)during the post-merger phase.This corresponds to infrared and visible radiation. Fora binary at z = 1 with the luminosity calculated above,this source has an apparent magnitude of m = 24 andshould be observable by the proposed LSST instrument[92]. Our simulations follow the late stage of the inspiralin which the binary separation decreases from d = 10 M to merger. For a 10 M ⊙ binary, this corresponds to atimescale of ∆ t ∼ . V ∞ ≫ a ∞ , but can be derived in a similar fash-ion. For a realistic gas with T > ∼ K ( a ∞ ≈ VII. SUMMARY AND DISCUSSION
In this paper we have performed a set of fully gen-eral relativistic simulations of BHBH binary mergers ingaseous environments, using our relativistic hydrody-namics code with AMR capability. Our focus has beenon demonstrating a mechanism for an observable elec-tromagnetic signal which may accompany the gravita-tional wave signal from a black hole merger. We have re-stricted our attention to gas clouds which are asymptoti-cally uniform and either at rest (BHBH Bondi accretion),or moving (BHBH BHL accretion) with respect to thebinary center of mass. We have performed “prototype”high-temperature ( T ∞ ∼ K) simulations in orderto gain a qualitative understanding of the different flowregimes characterizing such systems. In these simula-tions, the accretion radius R a is placed within the compu-tational domain, which allows us to sample all the variousregimes, studying how the luminosity and accretion ratechange when the binary separation d goes from d > R a to d < R a . We have also performed “realistic” simula-tions with astrophysically plausible asymptotic temper-atures ( T ∞ = 10 K ) to identify observable electromag-netic signals. For each simulation, we have calculatedthe time-varying rest-mass accretion rate, as well as theelectromagnetic luminosity due to bremsstrahlung andsynchrotron emission. We also derive scaling relationsfor the luminosity for the realistic cases, enabling our re-sults to be used for different asymptotic gas parametersand BH masses.In each case, we find evidence for a time-varying elec-tromagnetic signature accompanying the BHBH binarymerger. For our “realistic” temperature simulations,we find that during the final inspiral from separation d = 10 M to merger, there is an enhancement in thetotal luminosity of ∼ − ∼ − n ∞ ∼
10 cm − and T ∞ ∼ K) someof the assumptions that go into our calculation may bebreaking down, like the assumptions of adiabaticity, sub-Eddington luminosities and an optically thin medium.These are complications we intend to address in futurework. We also hope to compare our findings to other re-cent investigations performed simultaneously with ours(e.g. [37]), although these treat a different scenario andadopt very different initial data, making a direct com-parison difficult. Finally, we hope to use the simulationsreported here as point of comparison with future simula-tions involving BHBH mergers in ambient gaseous disks.
Acknowledgments
We would like to thank Z. Etienne and C. Gammiefor useful discussions. We would also like to thankM. Ansorg for providing the
TWOPUNCTURES code forgenerating BHBH metric initial data. This paper wassupported in part by NSF Grants PHY02-05155 andPHY06-50377 as well as NASA Grants NNG04GK54Gand NNX07AG96G. Simulations were performed underTeraGrid Grant TG-MCA99S008.
APPENDIX A: DERIVATION OF ˙ M EXPRESSION
Consider a 3D hypersurface F given by f ( t, x, y, z ) = constant in a spacetime diagram (see Fig. 24). The inter-section of F with a t = constant time slice is a 2-surface S ( t ). Below we will identify S ( t ) to be the apparent hori-9 t + t ∆ S(t )t L L
S(t + t) ∆ F F
FIG. 24: Spacetime diagram depicting the hypersurfaces rel-evant to calculating the mass accretion rate. The 2D hyper-surface S ( t ) and S ( t + ∆ t ) (white regions) represent BHhorizons on neighboring time slices. They are enclosed byspacelike 3D hypersurfaces Σ t and Σ t +∆ t (shaded region)on these slices. zon of a black hole at time t . Let L be a 3D hypersurfacewhich is a worldtube enclosing F . Let us further defineΣ t to be the 3D region on the time slice t between thesurface S ( t ) and Ω( t ), the 2D cross section of L on thetime slice t . We can imagine a 4-volume V to be the re-gion bounded by the hypersurfaces Σ t , Σ t + δt , L, and F .Consider a 4-vector flux j µ with a vanishing divergence, ∇ µ j µ = 0. At a given time t , define the function q ( t )according to q ( t ) = Z Σ t j µ d Π µ (A1)From Gauss’s law,0 = Z V ∇ µ j µ d V = Z ∂ V j µ d Π µ = q ( t + δt ) − q ( t )+ Z F j µ d Π µ − Z L j µ d Π µ . (A2)To evaluate the integral on F , it is convenient tointroduce a coordinate system ( t, f, a, b ), where a = a ( t, x, y, z ) and b = b ( t, x, y, z ) are two other coordinates.We may write Z F j µ d Π µ = 13! Z j µ ǫ µνρσ dx ν dx ρ dx σ = 13! Z j f ǫ ftab dt ∧ da ∧ db = − Z p − g ′ j f dtdadb = − Z √− gj µ ∂ µ f J dtdadb , (A3)where g ′ is the determinant of the metric in the ( t, f, a, b )coordinate system, g is the determinant in the ( t, x, y, z )coordinate system, J is the Jacobian J = (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( t, f, a, b ) ∂ ( t, x, y, z ) (cid:12)(cid:12)(cid:12)(cid:12) − = (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( f, a, b ) ∂ ( x, y, z ) (cid:12)(cid:12)(cid:12)(cid:12) − , (A4) and j f = j µ ∂ µ f follows from the usual transformationformula for a vector field between the ( t, f, a, b ) and( t, x, y, z ) coordinate systems. Similarly, if L is given by l ( t, x, y, z ) = 0, then Z L j µ d Π µ = − Z √− gj µ ∂ µ lJ l dtdadb , (A5)where J l = | ∂ ( l, a, b ) /∂ ( x, y, z ) | − . Taking the limit δt →
0, we obtain dqdt = −F F + F L , (A6)where F F = − Z F √− gj µ ∂ µ f J dadb , (A7) F L = − Z L √− gj µ ∂ µ lJ l dadb . (A8)Consider a fluid accreting onto a BH. Let F be theBH horizon world tube. At any given time t , considerthe fluid in the region Σ t between the BH (apparent)horizon S ( t ) and a distant 2-surface Ω( t ) surrounding theBH. Let L be the 3D hypersurface formed by stacking Ωwith time. The continuity equation gives ∇ µ ( ρ u µ ) = 0.Setting j µ = ρ u µ , we have q ( t ) = Z Σ t ρ u µ d Π µ = Z Σ t ρ ∗ d x = M ( t ) (A9)is the rest mass bounded by the surface Ω and the hori-zon, where ρ ∗ = √− g ρ u . Hence we have dM dt = F L − F F . (A10)This equation states that the rate of change of the restmass is equal to the amount of rest mass flowing into Ωper unit time ( F L ) minus the amount of rest mass flowinginto the horizon per unit time ( F F ). Hence we define therest-mass accretion rate onto the BH according to˙ M ≡ F F = − Z F α √ γρ u µ ∂ µ f Jdθdφ , (A11)which is the expression given in Eq. (27). Here we haveused the identity √− g = α √ γ , and choose a and b to bethe spherical angular coordinates θ and φ with the originat the BH center.It is apparent from the definition that in general ˙ M depends on how the spacetime is sliced near the horizon.However, in some cases ˙ M may be time independent.Consider the cases where the fluid’s mass is negligiblecompared to the BH’s mass and there exists a timelikeKilling vector ξ = ∂/∂λ in the vicinity of a BH. ThisKilling vector could be the time Killing vector describinga stationary BH, or a helical Killing vector which approx-imates the BHBH spacetime in the inspiral phase. One0might choose to measure ˙ M using a coordinate system inwhich t = λ (at least locally). In this case,˙ M ≡ ˙ M λ = − Z F √− g λ ρ u f dθdφ , (A12)where g λ is the determinant of the spacetime metricin the ( λ, f, θ, φ ) coordinate system. Suppose the flowof the fluid also achieves a stationary state near thehorizon in which ∂ λ ( ρ u f ) = 0 everywhere on F , i.e. ρ u f = ρ u f ( θ, φ ) on F . Since ∂ λ g λ = 0, ˙ M λ is a con-stant independent of the coordinate time λ . On the otherhand, in a coordinate system in which λ is not the timecoordinate, we have˙ M = − Z F p − g ′ ρ u f dθdφ = − Z F ( ξ t ) − √− g λ ρ u f dθdφ , (A13)where g ′ = g λ / ( ξ t ) is the determinant of the spacetimemetric in the ( t, f, θ, φ ) coordinate system, and ξ t is thetime component of the Killing vector ξ . Note that both g λ and ρ u f are still time independent on F , but ˙ M istime dependent if ξ t is time dependent. The accretionrate ˙ M is time independent only if a gauge is chosenso that ∂ t ξ t = 0 on F . Furthermore, if ξ t is constanteverywhere on F , then ˙ M = ˙ M λ /ξ th , where ξ th is thevalue of ξ t on F . As an example, consider a stationaryaccretion flow onto a Kerr BH. In the boosted Kerr-Schildcoordinates, we have ξ t = γ b = 1 / p − v b everywhere inthe spacetime, where v b is the boost velocity. Hence wehave ˙ M = ˙ M λ /γ b in the boosted Kerr-Schild coordinates.We point out that in a numerical simulation, even ifa Killing vector exists, the adopted gauge (i.e. time slic-ing) may not correspond to the gauge in which ∂/∂t isthe Killing vector. However, in some situations there ex-ists a gauge in which a (quasi)stationary flow is expected.Such situations include the Bondi accretion onto a BHBHbinary in the inspiring phase and the Bondi accretiononto a single BH following the binary merger. In thesesituations, the mass accretion rate onto a distant, fixedsurface Ω (i.e. F L ) is time independent for any gaugechoices that give rise to a spacetime that is asymptoti-cally Minkowsky. A “well-behaved” gauge should give a F F (or the sum of two F F ’s in the BHBH case) equal to F L (when averaged over time); otherwise, Eq. (A10) im-plies that there will be an accumulation (if F F < F L ) ordepletion (if F F > F L ) of rest mass in the interior of Ω asa result of a pure gauge effect. In our numerical simula-tions, we do not see such a gauge effect. We compute F L on spherical surfaces of various radii. We find that afterthe flow reaches a (quasi)stationary state, F L is slowlychanging with time in the binary inspirling phase and isapproximately time independent after merger. The com-puted fluxes at various radii are also the same. Moreover,the sum of the computed fluxes at the BH horizons agreewith the value F L , indicating that our adopted puncturegauge conditions are well-behaved. APPENDIX B: EMISSIVITIES a. Bremsstrahlung emissivity
In order to estimate the electromagnetic emission dueto bremsstrahlung, we use the following expressions forelectron-ion, and electron-electron cooling rates given in[93] q ff = q ei + q ee (B1) q ei = n π α f r e c ) ( m e c ) F ei ( θ ) ergs cm − s − (B2) q ee = n ( α f r e c )( m e c ) F ee ( θ ) ergs cm − s − (B3)Here, α f = e / ¯ hc is the fine structure constant, r e = e /m e c is the classical electron radius, θ ≡ kT /m e c , n = ρ /m B is the baryon number density, and F ei ( θ ) = (cid:18) θπ (cid:19) / (1 + 1 . θ . ) θ < θ π [ln(1 . θ + 0 .
48) + 1 . θ > F ee ( θ ) = π / (44 − π ) θ / × (1 + 1 . θ + θ − . θ / ) θ < θ (ln 1 . θ + 1 . θ > . (B5) b. Synchrotron emissivity We use the estimates for synchrotron cooling ratesgiven by [94]: q ν, s = 4 πne ν √ cK (1 /θ ) I (cid:16) x M sin θ (cid:17) ergs cm − s − Hz − , (B6)where ν = eB πm e c = cyclotron frequency (B7) x M = 2 ν ν θ . (B8)From [95], we get the following approximation for I ( x M ). I ( x M ) = 2 . . x / M + 0 . x / M ! exp( − . x / M ) , (B9)and the angle averaged version, I ′ ( x M ) = 4 . x / M . x / M + 0 . x / M ! exp( − . x / M ) . (B10)We integrate over frequency to get the total coolingrate.1Let us denote, A = 4 . πne √ cK (1 / ( θ )) (B11) a = 23 ν θ (B12) a = 1 . . (B13)Thus, we see, Z ∞ q ν, s dν = Aa Z ∞ x / M . x / M + 0 . x / M ! × exp (cid:16) − a x / M (cid:17) dx M = 3 ACa (B14)Where C ≡ Γ(11 / a / + 0 . / a / + 0 . a = 2 . q s = Z ∞ q ν, s dν = 4 . C n e c √ K (1 /θ ) θ βm e c . (B16) Here we have used P = βP M ≡ β B π . (B17)As discussed in Sec. IV, we have chosen β = 10 for oursimulations based on simulations of magnetized accretionflows which have demonstrated that the magnetic fieldsdo not typically reach their full equipartition value [83].We further assume that the hydrodynamic flow is tur-bulent due to the magneto-rotational instability (MRI),causing the frozen-in field lines to be tangled, and ran-domly oriented [88]. This justifies our use of the angleaveraged function in Eq. (B10). Note that for x ≪ K ( x ) ≈ /x , so that for sufficiently large temperatureswe may approximate: q s ≈ Z ∞ q ν, s dν = 4 . C n e c √ θ βm e c . (B18) [1] D. Richstone, E. A. Ajhar, R. Bender, G. Bower,A. Dressler, S. M. Faber, A. V. Filippenko, K. Gebhardt,R. Green, L. C. Ho, et al., Nature (London) , A14+(1998), arXiv:astro-ph/9810378.[2] B. M. Peterson and A. Wandel, Astrophys. J. Lett. ,L13 (2000), arXiv:astro-ph/0007147.[3] L. Ferrarese and H. Ford, Space Science Reviews ,523 (2005), arXiv:astro-ph/0411247.[4] M. C. Begelman, R. D. Blandford, and M. J. Rees, Na-ture (London) , 307 (1980).[5] N. Roos, Astron. and Astrophys. , 218 (1981).[6] D. Merritt and M. Milosavljevi´c, Living Reviews in Rel-ativity , 8 (2005), arXiv:astro-ph/0410364.[7] M. Milosavljevi´c and E. S. Phinney, Astrophys. J. Lett. , L93 (2005), arXiv:astro-ph/0410343.[8] H. L. Maness, G. B. Taylor, R. T. Zavala, A. B. Peck,and L. K. Pollack, Astrophys. J. , 123 (2004),arXiv:astro-ph/0310663.[9] C. Rodriguez, G. B. Taylor, R. T. Zavala, A. B. Peck,L. K. Pollack, and R. W. Romani, Astrophys. J. ,49 (2006), arXiv:astro-ph/0604042.[10] C. Rodriguez, G. B. Taylor, R. T. Zavala, Y. M.Pihlstr¨om, and A. B. Peck, Astrophys. J. , 37 (2009),0902.4444.[11] H. J. Lehto and M. J. Valtonen, Astrophys. J. , 207(1996).[12] M. J. Valtonen, H. J. Lehto, A. Sillanp¨a¨a, K. Nilsson,S. Mikkola, R. Hudec, M. Basta, H. Ter¨asranta, S. Haque,and H. Rampadarath, Astrophys. J. , 36 (2006).[13] M. J. Valtonen, H. J. Lehto, K. Nilsson, J. Heidt, L. O.Takalo, A. Sillanp¨a¨a, C. Villforth, M. Kidger, G. Poyner, T. Pursimo, et al., Nature (London) , 851 (2008),0809.1280.[14] J. K. Adelman-McCarthy, M. A. Ag¨ueros, S. S. Allam,K. S. J. Anderson, S. F. Anderson, J. Annis, N. A. Bah-call, C. A. L. Bailer-Jones, I. K. Baldry, J. C. Barentine,et al., Astrophys. J. Supp. , 634 (2007), 0707.3380.[15] S. Komossa, H. Zhou, and H. Lu, Astrophys. J. Lett. , L81 (2008), 0804.4585.[16] T. A. Boroson and T. R. Lauer, Nature (London) ,53 (2009), 0901.3779.[17] B. Kocsis, Z. Haiman, and K. Menou, Astrophys. J. ,870 (2008), 0712.1144.[18] C. Deffayet and K. Menou, Astrophys. J. Lett. , L143(2007), 0709.0003.[19] D. E. Holz and S. A. Hughes, Astrophys. J. , 15(2005), arXiv:astro-ph/0504616.[20] B. Kocsis, Z. Frei, Z. Haiman, and K. Menou, Astrophys.J. , 27 (2006), arXiv:astro-ph/0505394.[21] P. Artymowicz and S. H. Lubow, Astrophys. J. , 651(1994).[22] R. G¨unther and W. Kley, Astron. and Astrophys. ,550 (2002), arXiv:astro-ph/0204175.[23] A. Escala, R. B. Larson, P. S. Coppi, and D. Mardones,Astrophys. J. , 152 (2005), arXiv:astro-ph/0406304.[24] A. I. MacFadyen and M. Milosavljevi´c, Astrophys. J. , 83 (2008), arXiv:astro-ph/0607467.[25] S. L. Shapiro, submitted to PRD.[26] L. R. Corrales, Z. Haiman, and A. MacFadyen, ArXive-prints (2009), 0910.0014.[27] S. M. O’Neill, M. C. Miller, T. Bogdanovi´c, C. S.Reynolds, and J. D. Schnittman, Astrophys. J. , 859 (2009), 0812.4874.[28] M. Megevand, M. Anderson, J. Frank, E. W.Hirschmann, L. Lehner, S. L. Liebling, P. M. Motl, andD. Neilsen, Phys. Rev. D , 024012 (2009), 0905.3390.[29] M. Campanelli, C. O. Lousto, Y. Zlochower, andD. Merritt, Physical Review Letters , 231102 (2007),arXiv:gr-qc/0702133.[30] T. Bogdanovi´c, C. S. Reynolds, and M. C. Miller, Astro-phys. J. Lett. , L147 (2007), arXiv:astro-ph/0703054.[31] M. Anderson, L. Lehner, M. Megevand, and D. Neilsen,ArXiv e-prints (2009), 0910.4969.[32] J. D. Schnittman and J. H. Krolik, Astrophys. J. ,835 (2008), 0802.3556.[33] G. A. Shields and E. W. Bonning, Astrophys. J. ,758 (2008), 0802.3873.[34] Z. Lippai, Z. Frei, and Z. Haiman, Astrophys. J. Lett. , L5 (2008), 0801.0739.[35] E. M. Rossi, G. Lodato, P. J. Armitage, J. E. Pringle,and A. R. King, Mon. Not. R. Astron. Soc. pp. 1726–+(2009), 0910.0002.[36] J. R. van Meter, J. H. Wise, M. C. Miller, C. S. Reynolds,J. M. Centrella, J. G. Baker, W. D. Boggs, B. J. Kelly,and S. T. McWilliams, ArXiv e-prints (2009), 0908.0023.[37] T. Bode, R. Haas, T. Bogdanovic, P. Laguna, andD. Shoemaker, ArXiv e-prints (2009), 0912.0087.[38] M. Shibata and T. Nakamura, Phys. Rev. D , 5428(1995).[39] T. W. Baumgarte and S. L. Shapiro, Physical Review D , 024007 (1999).[40] F. Pretorius, Classical and Quantum Gravity , 425(2005), arXiv:gr-qc/0407110.[41] F. Pretorius, Physical Review Letters , 121101 (2005),arXiv:gr-qc/0507014.[42] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo-chower, Physical Review Letters , 111101 (2006),arXiv:gr-qc/0511048.[43] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, andJ. van Meter, Physical Review Letters , 111102 (2006),arXiv:gr-qc/0511103.[44] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, Classicaland Quantum Gravity , 165003 (2008).[45] I. D. Novikov and K. S. Thorne, in Black Holes (LesAstres Occlus) (1973), pp. 343–450.[46] P. E. J. Nulsen and A. C. Fabian, Mon. Not. R. Astron.Soc. , 346 (2000), arXiv:astro-ph/9908282.[47] H. Bondi, Mon. Not. R. Astron. Soc. , 195 (1952).[48] F. C. Michel, Astrophys. Space Sci. , 153 (1972).[49] S. L. Shapiro, Astrophys. J. , 531 (1973).[50] S. L. Shapiro and S. A. Teukolsky, Black holes, whitedwarfs, and neutron stars: The physics of compact ob-jects (Wiley, New York, 1983).[51] F. Hoyle and R. A. Lyttleton, in
Proceedings of the Cam-bridge Philosophical Society (1939), vol. 35 of
Proceedingsof the Cambridge Philosophical Society , pp. 405–+.[52] H. Bondi and F. Hoyle, Mon. Not. R. Astron. Soc. ,273 (1944).[53] L. I. Petrich, S. L. Shapiro, R. F. Stark, and S. A. Teukol-sky, Astrophys. J. , 313 (1989).[54] R. Hunt, Mon. Not. R. Astron. Soc. , 141 (1971).[55] E. Shima, T. Matsuda, H. Takeda, and K. Sawada, Mon.Not. R. Astron. Soc. , 367 (1985).[56] J. A. Font and J. M. A. Ibanez, Astrophys. J. , 297(1998).[57] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys. Rev. D , 024028 (2005), arXiv:astro-ph/0503420.[58] J. R. van Meter, J. G. Baker, M. Koppitz, and D.-I. Choi,Phys. Rev. D , 124011 (2006), arXiv:gr-qc/0605030.[59] B. Fields and S. Sarkar, ArXiv Astrophysics e-prints(2006), arXiv:astro-ph/0601514.[60] M. Ansorg, B. Br¨ugmann, and W. Tichy, Phys. Rev. D , 064011 (2004), arXiv:gr-qc/0404056.[61] W. Tichy and B. Br¨ugmann, Phys. Rev. D , 024006(2004), arXiv:gr-qc/0307027.[62] W. Tichy and P. Marronetti, Phys. Rev. D , 081501(2008), 0807.2985.[63] J. A. Faber, T. W. Baumgarte, Z. B. Etienne, S. L.Shapiro, and K. Taniguchi, Phys. Rev. D , 104021(2007), 0708.2436.[64] .[65] E. Schnetter, S. H. Hawley, and I. Hawke, Classical andQuantum Gravity , 1465 (2004), arXiv:gr-qc/0310042.[66] J. Thornburg, Classical and Quantum Gravity , 743(2004), arXiv:gr-qc/0306056.[67] B. van Leer, Journal of Computational Physics , 276(1977).[68] A. Harten, P. D. Lax, and v. B. J., SIAM Rev. , 35(1983).[69] Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro,K. Taniguchi, and T. W. Baumgarte, Phys. Rev. D ,084002 (2008), 0712.2460.[70] R. W. Lindquist, Annals of Physics , 487 (1966).[71] K. S. Thorne, Mon. Not. R. Astron. Soc. , 439 (1981).[72] P. J. Schinder, Phys. Rev. D , 1673 (1988).[73] S. L. Shapiro, Phys. Rev. D , 1858 (1989).[74] P. J. Schinder and S. A. Bludman, Astrophys. J. ,350 (1989).[75] A. Mezzacappa and R. A. Matzner, Astrophys. J. ,853 (1989).[76] L. Rezzolla and J. C. Miller, Classical and QuantumGravity , 1815 (1994), arXiv:astro-ph/9406055.[77] L. Zampieri, J. C. Miller, and R. Turolla, Mon. Not. R.Astron. Soc. , 1183 (1996), arXiv:astro-ph/9607030.[78] S. L. Shapiro, Astrophys. J. , 308 (1996).[79] L. Rezzolla and J. C. Miller, Phys. Rev. D , 5411(1996).[80] S. Balberg, L. Zampieri, and S. L. Shapiro, Astrophys. J. , 860 (2000), arXiv:astro-ph/0004234.[81] M. Liebend¨orfer, O. E. B. Messer, A. Mezzacappa, S. W.Bruenn, C. Y. Cardall, and F.-K. Thielemann, Astro-phys. J. Supp. , 263 (2004), arXiv:astro-ph/0207036.[82] B. D. Farris, T. K. Li, Y. T. Liu, and S. L. Shapiro, Phys.Rev. D , 024023 (2008), 0802.3210.[83] J. C. McKinney and C. F. Gammie, Astrophys. J. ,977 (2004), arXiv:astro-ph/0404512.[84] Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys. Rev.D , 084017 (2007), 0706.2360.[85] Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro,and T. W. Baumgarte, Phys. Rev. D , 101503 (2007),0707.2083.[86] Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baum-garte, Phys. Rev. D , 044024 (2009), 0812.2245.[87] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.Rev. D , 061501 (2006), arXiv:gr-qc/0601091.[88] S. L. Shapiro, Astrophys. J. , 69 (1973).[89] S. L. Shapiro, Astrophys. J. , 343 (1974).[90] T. W. Baumgarte, Phys. Rev. D , 024018 (2000),arXiv:gr-qc/0004050. [91] L. D. Landau and E. M. Lifshitz, Fluid mechanics (Else-vier, Oxford, 1959).[92] .[93] R. Narayan and I. Yi, Astrophys. J. , 710 (1995),arXiv:astro-ph/9411059.[94] A. G. Pacholczyk,
Radio astrophysics. Nonthermal pro- cesses in galactic and extragalactic sources (Series ofBooks in Astronomy and Astrophysics, San Francisco:Freeman, 1970, 1970).[95] R. Mahadevan, R. Narayan, and I. Yi, Astrophys. J.465