NASA Technical Reports Server (NTRS) 20140005676: Accretion Disks Around Binary Black Holes of Unequal Mass: GRMHD Simulations Near Decoupling

3.5 MB·English

by  NASA Technical Reports Server (NTRS)

Accretion disks around binary black holes of unequal mass:
GRMHD simulations near decoupling
Roman Gold1, Vasileios Paschalidis1, Zachariah B. Etienne1,2,3,4, Stuart L. Shapiro1,5
1Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801
2Physics Department & Joint Space-Science Institute,
University of Maryland, College Park, MD 20742
3Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771
4Department of Mathematics, West Virginia University, Morgantown, WV 26506
5Department of Astronomy & NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801
and
Harald P. Pfeiffer6,7
6Canadian Institute for Theoretical Astrophysics, University of Toronto, ON M5S 3H8, Canada
3 7Canadian Institute for Advanced Research, Toronto, ON M5G 1Z8, Canada
1
0 We report on simulations in general relativity of magnetized disks onto black hole binaries. We
2 varythebinarymassratiofrom1:1to1:10andevolvethesystemswhentheyorbitnearthebinary-
c disk decoupling radius. We compare (surface) density profiles, accretion rates (relative to a single,
e non-spinningblackhole),variability,effectiveα-stresslevelsandluminositiesasfunctionsofthemass
D ratio. We treat the disks in two limiting regimes: rapid radiative cooling and no radiative cooling.
Themagneticfieldlinesclearlyrevealjetsemergingfrombothblackholehorizonsandmerginginto

2 one common jet at large distances. The magnetic fields give rise to much stronger shock heating
than the pure hydrodynamic flows, completely alter the disk structure, and boost accretion rates
] andluminosities. Accretionstreamsnearthehorizonsareamongthedenseststructures;infact,the
E
1:10no-coolingevolutionresultsinarefillingofthecavity. Thetypicaleffectivetemperatureinthe
H bulk of the disk is ∼ 105(M/108M(cid:2))−1/4(L/Ledd)1/4K yielding characteristic thermal frequencies
h. ∼1015(M/108M(cid:2))−1/4(L/Ledd)1/4(1+z)−1Hz. Thesesystemsarethuspromisingtargetsformany
p extragalactic optical surveys, such as LSST, WFIRST, and PanSTARRS.
-
o PACSnumbers: 04.25.D-,04.25.dg,47.75.+f
r
t
s
a I. INTRODUCTION AND MOTIVATION h ∼ 10−16(d/10M)−1(M/108M(cid:2))(D/16Gpc)−1. Here
[ M is the total mass of the binary, and we normal-
1 Supermassive black hole (SMBH) binaries can form ized to a luminosity distance D corresponding to red-
v shift z (cid:3) 2 in a standard ΛCDM Universe. The
in magnetized plasma following galaxy mergers, via bar-
0 corresponding gravitational-wave frequency is f ∼
mode instability in rapidly rotating supermassive stars, (cid:2) (cid:3) GW
60 or by other dynamical processes [1]. After formation, a 10−4 M/108M(cid:2) −1(d/10M)−3/2[(1+z)/3]−1Hz. The
0 combination of dynamical friction and gas-driven migra- expected GW strain is above the eLISA strain sensitiv-
. tionislikelytocatalyzethebinaryinspiralintothegrav- ity at these frequencies [15], hence these systems will be
12 itationalradiation-drivenregime[1–7]. Theexactdetails detectable by eLISA. In particular, for M ∼ 106M(cid:2),
3 of these processes, including the “last-parsec problem”, the value of fGW at d = 10M falls well within the
1 remain active areas of research (see e.g. [8–10] and refer- eLISA sensitivity band even for larger redshifts, while
: ences therein). As a result, event rates and population for M ∼ 108M(cid:2) and large redshifts (z ∼ 10), the value
iv synthesis studies at this stage are highly uncertain [11]. of fGW at d=10M is marginally within the eLISA sen-
X The exciting prospect of a simultaneous observation sitivity band. However, the inspiral time from these sep-
r of both electromagnetic (EM) and gravitational waves arations is tGW ∼ 20(M/108M(cid:2)) days assuming equal
a (GWs)arisingfromaccretingbinaryBHBHsmakesthese masses. Hence, these systems will quickly enter the
systems prime targets in the era of multi-messenger as- eLISA sweet spot, and EM precursor signals can trig-
tronomy. Such observations will enable us to determine ger targeted GW searches with a convenient lead time of
thebinarymasses,BHspins,redshiftandevendetermine several days.
the Hubble constant to better than 1% [12–14]. While awaiting the first detection of GWs, currently
Gravitationalwaves(GWs)fromSMBHsareexpected operating and future electromagnetic (EM) detectors
to be detected by planned GW interferometers such as such as LSST [19], WFIRST [20] and PanSTARRS [21]
eLISA/NGO detectors [15], sensitive to GW frequencies are promising instruments to identify accreting BHBH
10−5−1Hz, and the currently operating Pulsar Timing systems in the EM spectrum. Important steps have al-
Arrays [16–18], sensitive to frequencies 10−9 −10−6Hz. ready been made toward realizing this goal.
As SMBHs are typically believed to have masses in the Currently, we know of one spatially resolved SMBH
range 106 − 109M(cid:2), the GW strain at orbital separa- binary candidate at an orbital separation d ∼ 7pc:
tion d = 10M - the value adopted in this work - is 0402+379 [22]. Other spatially unresolved, SMBH bi-
2
nary candidates have been found, including OJ287 [23– t < t < t the disk can follow the inspiral of the
vis GW sd
26] and SDSS J1536+0441 [27, 28]. Binary AGN can- binary and settle in a quasi-steady state - this is called
didates have been singled out based on offsets in the the predecoupling regime. When t < t , the binary
GW vis
broad line and narrow line regions, emission line pro- decouples from the disk, i.e. the inward migration of the
files, and time variability [29, 30]. Recently, very-long binary out-paces the inward drift of the disk.
baselineinterferometricobservationsinterpretedejection In this paper we focus on the phase near decoupling,
components from AGN cores as undulations caused by while the postdecoupling regime will be the subject of a
the precession of the accretion disks around a SMBH bi- futurepaper. WenotethatunlikeeLISA,PulsarTiming
nary [31]. A simplified model was applied to two AGN Arrays are sensitive only to SMBH binaries well within
sources; for 1823+568 their analysis yields d ∼ 0.42pc, the predecoupling regime.
and a mass ratio q in the range 0.095 < q < 0.25, while Accretion onto a single BH has been studied in great
for 3C 279 d ∼ 2.7pc, q ∼ 0.36. However, given the detail for decades, and magnetohydrodynamic studies in
lackofarobustcircumbinaryaccretiondisktheorythese GRhavedrasticallyimprovedourunderstandingofthese
results are at best preliminary. Nevertheless, they still flows (see [50] for a recent review). Many different disk
motivate a study of accretion flows onto supermassive modelshavebeenproposedintheliterature. Thesemod-
BH binaries with different mass ratios, which we initiate els range from geometrically-thin, optically thick disks
in this work in general relativity (GR). [51, 52] and slim disks [53], to geometrically thick, opti-
Otherfeaturesidentifiedas“smokingguns”forbinary cally thin, radiatively inefficient accretion flows (RIAF)
black holes include BH recoil/kicks [32], used to explain [54–58]. However, our understanding of accretion flows
the large velocity offsets between emission lines in AGN onto BHBHs remains poor and studies of these systems
spectra, as well as observed kinks in jets probably due are still in their infancy.
to changes in BH spin (X-shaped radio sources), a past The first analytic Newtonian model and smooth par-
mergerevent,orprecessioneffects[33–35]. Modifications ticle hydrodynamic simulation of a circumbinary accre-
to the line profiles have also been proposed as promis- tion disk was given in [59]. Since then, other Newto-
ingcharacteristicfeaturestodistinguishbinaryBHAGN nian(semi-)analyticstudies[40–42,60,61]andhydrody-
from classical, single BH AGN sources [36]. namicsimulationsin2D[62–65],and3D[66–70]havefol-
To assist and solidify all these detection efforts, it is lowed. Newtonian magnetohydrodynamic (MHD) simu-
crucial to identify and model possible electromagnetic lationswerepresentedin[71],andPost-NewtonianMHD
“precursor” and “afterglow” signatures [37–43]. simulations in [72]. Many of these earlier studies ex-
Depending on the physical regime the properties char- cluded the binary and most of the inner cavity from
acterizing the gas can differ considerably, and different the computational domain, introducing an artificial in-
accretion models are applicable. For example, if the gas ner boundary condition. The importance of treating the
has little angular momentum, the accretion flow resem- inner regions self-consistently has been discussed in [69],
bles the binary analog of a Bondi-Hoyle-Lyttleton solu- and only full GR calculations can achieve this goal reli-
tion [44–46] (see [47–49] for GR studies). If the gas has ably.
significant angular momentum, then the gas can become A“GR-hybrid”orbit-averagedmodelforthindisks,in
rotationally supported and form a disk. which the viscous part is handled in GR and the tidal
ForaBHBHembeddedinadisk, onecanidentifysev- torquesinNewtoniangravitywasintroducedin[43],and
eral different regimes based on the time scales for mi- GRhydrodynamicalsimulationsofaccretionontoBHBH
gration of the binary. For a SMBH binary engulfed by binaries - taking into account the dynamical spacetime -
a (thin) disk at large separation, the migration of the havebeenperformedin[47,73–75]. TodatetheonlyGR
binary is initially governed by binary-disk angular mo- magnetohydrodynamic(GRMHD)simulationsofdiskac-
mentum exchange mediated by (effective) viscosity [3]. cretionontoBHBHsthataccountbothforthedynamical
At large enough separations, the reduced mass μ of the spacetime and the BH horizons were presented in [76].
binary is less than the local interior disk mass (4πd2Σ UsingadifferentapproachbyassumingthattheB-field
where Σ is the surface density of the disk). This leads is anchored to a circumbinary disk outside the computa-
to the so-called disk-dominated type II migration occur- tionaldomain,[77–80]modeledEMsignaturesbysolving
ring on the viscous time scale t . As the migration the GR force-free electrodynamic equations.
vis
proceeds,thereducedmassofthesystembecomeslarger Close to merger a single BH remnant is formed on a
than the local disk mass and the migration enters the time scale much shorter than the dynamical time scale
secondary-dominated type II regime, which occurs on a in the disk. The mass and angular momentum of the
longer time scale t ≡(4πd2Σ/μ)−k×t ≥t , where remnant BH is different from the total mass and angular
sd vis vis
k a constant of order ∼ 0.4 if 4πd2Σ/μ < 1 and 0 oth- momentumpriortomergerduetoGWemission,causing
erwise. Ultimately, the binary enters a regime at smaller a quasi-instantaneous perturbation to the disk. This ef-
orbital separations where angular momentum losses due fect has been modeled using hydrodynamical and MHD
to GWs dominate, and the binary migrates on the GW simulations in Newtonian gravity and in GR [81–85].
time scale t . In all regimes, the disk moves inwards A realistic and ideal 3D global model for a circumbi-
GW
on the viscous time scale. When t ≤ t < t or nary disk around a SMBH binary near the decoupling
vis sd GW
3
radius requires: a) a fully relativistic treatment of grav- section). Furthermore, as discussed in [3], steady-state
itation in a dynamical spacetime, b) GRMHD for the disk models predict that radiation pressure-dominated
plasma flow, c) realistic cooling processes, and d) radia- circumbinary disks will channel more material into the
tive transfer in curved spacetime. Simulations incorpo- cavity for a given central mass. This finding makes radi-
rating these effects must also have high resolution and ation pressure-dominated circumbinary disks promising
longintegrationtimes(severalviscoustimescales). How- sources for electromagnetic counterparts.
ever, including all these ingredients in one simulation
Second, AGN data from the Sloan Digital Sky Survey
would make these computations prohibitive, because the
(SDSS) [96–98], the AGN and Galaxy Evolution Survey
wall-clock times required to integrate for even ∼ 5 vis-
(AGES)[99],XMM-COSMOS[100]andothers[101]cov-
coustimescalesattheinnerdiskradiusathighresolution
ering mostly type I AGNs and the local Universe to red-
are far too long with current supercomputer resources. shiftofz (cid:2)5,revealEddingtonratiodistributionsinthe
Thus, some of these ingredients must be relaxed in order range 0.01(cid:2)L/L (cid:2)1 with the tendency that higher
Edd
to obtain a qualitative understanding of these systems. L/L values occur at higher z (and possibly smaller
Edd
For this reason, the models in this work feature a) and central mass M). A comparison of accretion time scales
b). High resolution is only adopted for a few models. In
with the age of the Universe suggests that earlier ac-
addition, we model radiative cooling by a simple cool-
cretion episodes are closer to the Eddington limit, sim-
ing function and consider the extreme opposite limits of ilar to the recently discovered quasar at z ∼ 7.1 with
“rapid” cooling and “no cooling” to bracket the possibil- a central mass of M ∼ 2 × 109M(cid:2) [102] accreting at
ities. L/L ∼ 1.2±0.5. These surveys therefore motivate
Edd
In this paper we study the effects of the binary mass studiesofdisksaccretingneartheEddingtonluminosity,
ratio on the disk near decoupling. We vary the BHBH for which radiation pressure is important. Moreover, ra-
binary mass ratio q ≡ M1/M2 ≤ 1, considering 1:1, 1:2, diation pressure-dominated disks accreting near the Ed-
1:4, 1:8 and 1:10 mass ratios. The mass ratio regime dington limit are more likely to be detectable, even at
0.1(cid:2)q (cid:2)1isshowntobeofhighastrophysicalrelevance large cosmological redshifts.
in[86,87]. Also,ithasbeenarguedthatBHBHsforming
This paper is structured as follows: In Sec. II we
in major galaxy mergers will typically have mass ratios
present a qualitative discussion of the range of param-
in the range (0.01,1) [88]. These results motivate our
eters and the associated physical regimes for which our
choice of mass ratios. Note that Newtonian simulations
simulations are appropriate. In Sec. III we describe our
studying the effects of mass ratio in the context of 2D
methodsandtechniquesforevolvingthespacetime,fluid,
thin-disk models were performed in [64, 65].
and magnetic fields, as well as our simple cooling pre-
The new aspects of our work are the inclusion of rela-
scription. Wealsopresentthedifferentcasesweconsider
tivisticgravitationtoresolvethecrucialphysicsnearthe
and list our diagnostics for characterizing the accretion
BH horizons, effective viscosity arising from the magne-
flowandEMsignatures. InSec.IVweshowseveraltests
torotationalinstability[89,90],geometricallythickdisks,
weperformedtomotivateournumericalsetup. InSec.V
three spatial dimensions (3D), a Γ-law equation of state
wepresenttheresultsfromoursimulations. Weconclude
(EOS), and effective cooling. We model the decoupling
inSec.VIwithasummaryofthemainresultsandadis-
epoch by fixing the binary separation d, evolving the
cussion of future work.
spacetime by rotating our conformal-thin-sandwich ini-
Here and throughout we adopt geometrized units,
tial data [91–94] as we have done in [76, 95]. The depen-
where G=c=1, unless otherwise stated.
denceoftheflowvariability,EMsignatures,themagnetic
field structure and the matter dynamics inside the low-
density“hollow”onthemassratioqareinvestigated. We
estimate thermodynamic properties of the gas and scale
our results with the binary total mass and the luminos-
II. QUALITATIVE CONSIDERATIONS
ity in units of the Eddington luminosity where feasible.
Fromthis,emissioncharacteristicsincludingtypicalther-
mal frequencies and luminosities are given and relevant In this section we use the Shakura-Sunyaev/Novikov-
detectors are discussed. Thorne thin-disk model [51, 52] to make rough analytic
In the absence of any observational constraints on the estimates regarding the physical regime our disk mod-
thermodynamicstateofaccretingBHBHs,themodelswe els probe. While the Shakura-Sunyaev/Novikov-Thorne
consider in this work adopt an adiabatic EOS governed model strictly applies to a viscous flow onto a single BH
by an adiabatic index Γ=4/3, appropriate for radiation neglecting the binary tidal torques, we expect that our
pressure-dominated, optically thick disks. This choice qualitative analysis will apply roughly to our circumbi-
is motivated in part both by theory and observations. nary disk models, which have H/R = O(0.1), where H
First, accretion disk theory [50] predicts that the inner is the disk scale height. This expectation should be best
partofcircumbinarydisksaroundSMBHbinariesareop- realized outside the binary orbit because the ratio of the
ticallythick,radiationpressuredominatedforalargeset tidaltoviscoustorquesdecaysquicklyfarawayfromthe
of possible disk parameters (see e.g. [4, 52] and the next binary orbit (see e.g. Fig. 6 in [43]).
4
A. Radiation pressure dominance In this region the typical rest mass densities are [50]
(cid:4) (cid:5) (cid:4) (cid:5)
r 3/2 α −1
1. Shakura-Sunyaev/Novikov-Thorne model ρ0 ∼5.5·10−12 20M 0.1
(cid:6) (cid:7)(cid:8) (cid:9)−2
M M˙ g
The simulations reported here apply to any total (3)
(ADM) binary black hole mass M and, since we neglect 108M(cid:2) 2.3M(cid:2)/yr cm3
tphreovsiedlef-dgritavhiatsytohfetshaemgeaisn,ittioaladniyskrpesrto-fimleasasdodpetnesdityheρr0e,. ∼5.5·10−12(cid:4) r (cid:5)3/2(cid:4) α (cid:5)−1
20M 0.1
The quantities that are fixed, in addition to the initial (cid:6) M (cid:7)−1(cid:6) L (cid:7)−2(cid:4) ε (cid:5)2 g
disk profile, are the adiabatic index Γ appearing in the . (4)
adopted ideal gas EOS, the ratio of the initial magnetic- 108M(cid:2) Ledd 0.1 cm3
to-total disk pressure, the initial magnetic field profile, As we will show later, these characteristic values for
the initial binary separation and the cooling law. In P /P and ρ are comparable to the values found in
rad gas 0
thissectionweusethesteady-statethin-disksolutionfor our simulations.
a Shakura-Sunyaev/Novikov-Thorne disk about a single TheShakura-Sunyaev/Novikov-Thornemodelpredicts
BH of mass M to estimate the physical values of some thattheeffectiveopticaldepthtoabsorptionisτ∗ (cid:3)0.02
of the gas dynamical quantities as functions of M and adopting the same normalizations as in the equation
luminosity L. We show below that for a range of astro- above. This well-known inconsistency near the Edding-
physically relevant choices of these parameters the disk ton limit of the Novikov-Thorne model is removed with
is thermal radiation pressure-dominated, and this fact thegeneralizationtotheslimdiskmodel[53]. Thismodel
motivates our setting Γ=4/3. differs from a thin disk in that it allows for cooling to
Neglectingtheperturbativeroleofthesecondarytidal occur via advection, which dominates radiative cooling
torque, the steady-state accretion flow in a geometri- at high accretion rates (L (cid:3) 0.3L ), thereby puffing
Edd
cally thin disk is uniquely specified by the Shakura- up the disk. When scaling our models to accrete near
Sunyaev viscosity parameter α, the central mass M the Eddington limit, they are closer to slim-disk models,
and the accretion rate M˙ . The quantity M˙ is spec- whichremainopticallythickinthishighluminositylimit
ified in turn by the disk luminosity L and efficiency [50, 53].
ε. TheShakura-Sunyaev/Novikov-Thornemodel[51,52] As we discuss later, near the Eddington limit and
describes a disk that is radiation pressure-dominated in- M = 108M(cid:2), the effective optical depth satisfies τ(cid:2) (cid:3) 1
side a radius rinner. In this region the radiation to gas in the bulk of our disk models, implying that the gas
pressure ratio is [50] is optically thick to absorption and the photons eventu-
(cid:4) (cid:5) (cid:4) (cid:5) ally thermalize. Thus, these qualitative considerations
Prad ∼5.4×105 r −21/8 α 1/4 motivate the adoption of an adiabatic index Γ = 4/3,
P 20M 0.1 appropriate for thermal radiation pressure dominance.
gas
(cid:6) M (cid:7)1/4(cid:6) L (cid:7)2(cid:4) ε (cid:5)2 Note also that alternative disk models have been pro-
, (1) posed,e.g.[4]. However,theylargelysharetheprediction
108M(cid:2) Ledd 0.1 that the inner regions of the disk are radiation pressure
dominated.
whereε≡L/M˙ c2istheradiativeefficiency,L (cid:3)1.3×
Edd
1046(M/108M(cid:2))erg/s, and we chose the normalization
for the α parameter based on typical values found in our 2. Decoupling radius
simulations. For a thin-disk model the efficiency ranges
from0.057foranon-spinningBHto0.42foramaximally We estimate the decoupling radius ad by equating
spinningBH,sowescaletoavalueresidingbetweenthese t =t and solve for the separation to find [76]
GW vis
limits. The size of the inner region is determined by the
(cid:6) (cid:7) (cid:6) (cid:7) (cid:6) (cid:7)
condition Prad/Pgas =1, for which Eq. (1) yields ad ≈13.3 α −2/5 H/R −4/5 η˜ 2/5 (5)
(cid:4) α (cid:5)2/21(cid:6) M (cid:7)2/21 M 0.1 0.3 1
rinner/M ∼3(cid:6)000L 0(cid:7).116/21(cid:4) ε1(cid:5)08−M16(cid:2)/21 αwtlehsiesrteothwreienSa≈hsasku1u.m5rdead-aSstuhntayytpaeitcvhaeltlyuinrfnboeuurlneddnitsiknvioesdcuogrseistirymadupilauartsaiomsneets--,
. (2)
L 0.1 ter, and η˜ ≡ 4q/(1+q)2 (see also [5]). Notice that in
edd
contrast to geometrically thick disks, where the decou-
The geometrically thick disks we evolve in this work pling radius is a few tens of M, geometrically thin disks
extendradiallyouttorout ∼100M−200M. Hence,when have ad/M (cid:3) 100. The decoupling radius estimate (5)
scaling the accretion rate such that our models accrete for the mass ratios considered in this work ranges from
near the Eddington limit, our models are fully immersed ad(q =0.1)≈8.5M to ad(q =1)≈13.3M. The normal-
in this inner radiation pressure-dominated region. izations in Eq. (5) are based on typical values obtained
5
from our simulations. We choose the binary separation
d∼10M forallcasesstudiedinthiswork, avaluewhich TABLE I. CTS initial data parameters for the BHBH vac-
uum spacetime. Columns show mass ratio (q), ADM mass
is consistent with the crude estimate of Eq. (5). In the
(MADM) and angular momentum (JADM), and irreducible
future we intend to start our evolutions at larger separa- masses (Mi ), and apparent horizon radii (ri ) for the two
irr hor
tions in order to dynamically determine the decoupling
blackholes. Diagnosticsgeneratingthesequantities,butcom-
radius and evolve through it as in [43]. puted from vacuum, test simulations agree with these values
to within one part in 104.
q MADM JADM Mi1rr Mi2rr rh1or rh2or
III. METHODS
1:1 0.98989 0.96865 0.50000 0.50000 0.42958 0.42958
1:2 0.99097 0.85933 0.66667 0.33333 0.60192 0.27312
The models we adopt here assume: a) circular binary
1:4 0.99342 0.61603 0.80000 0.20000 0.75140 0.15832
orbits, neglecting the binary inspiral (justified for large
1:8 0.99589 0.37868 0.88889 0.11111 0.85640 0.08618
separations;seeSec.IIA2),b)nonself-gravitatingdisks,
whichlikelyisagoodassumption(see,e.g. [103]),c)ideal 1:10 0.99656 0.31652 0.90909 0.09091 0.88081 0.07022
MHD, d) no radiative feedback, e) an effective cooling
scheme that brackets no cooling and rapid cooling.
not identical to [76] due to the different polytropic index
Someoftheseassumptionsmaynotbeobeyedstrictly,
(Γ=4/3 here versus Γ=5/3 in [76]).
e.g. the binary may become eccentric in the predecou-
Weseedtheinitialdiskwithasmall,purelypoloidalB-
pling regime [6, 66–69, 104], or radiative feedback may
field using the same procedure as in [76, 110]; see Fig. 1.
becomeimportantnearEddingtonaccretionrates. How-
Thefieldisdynamicallyunimportantinitially: theinitial
ever, our simulations still provide a qualitative under-
standing of the physics that will be useful in designing maximum value for the ratio of magnetic PM to total
pressure P is 0.025. All cases we consider in this work
thenextsetofmorerealisticmodelsofbinaryblackholes
are initialized with the same disk and magnetic field.
immersed in circumbinary disks. In this section, we de-
scribe the initial data and computational methods we
adopt to account for a)-e).
A. Initial data and AMR hierarchy
1. Metric initial data
At large separation the binary inspiral time scale is
much longer than the binary orbital period and the vis-
cous time scale at the inner edge of the disk just be-
yond the binary orbit. Accordingly, the inspiral can
be neglected over many orbital periods. To model this
epochinGR,weadoptquasi-equilibriumconformal-thin-
sandwich (CTS) solutions for the black hole binary [91–
FIG. 1. B-field lines (white curves) and volume rendering of
94]. The spacetime initial data satisfying the CTS equa-
rest-mass density normalized to its maximum value at t=0.
tions correspond to a circular orbit and possess a helical The two black dots at the center indicate the BH apparent
Killingvector. TheCTSinitialdatahavebeengenerated horizons for the q=1 case.
usingthespectraltechniquesdescribedin[105]asimple-
mented in the Spectral Einstein Code (SpEC) [106] (see Although we will qualitatively discuss how the evolu-
also[107]). Welisttheinitialdataparametersdescribing tion of these matter initial data depends on the binary
our spacetimes in Table I. mass ratio, it is important to stress that the goal of our
workistoassesshowthefinalrelaxedstateofthediskde-
pends on the binary mass ratio. The initial data, which
governs the early evolution, has no physical significance.
2. Matter and B-field initial data
B. Evolution equations and techniques
For the fluid we use the same family of equilibrium
disk models around single BHs as in [95, 108, 109]. We 1. GRMHD evolution
choosetheinitialinnerdiskedgeradiusrin,0 =18M and
specific angular momentum l(rin,0) = 5.15M2 around a WeusetheIllinoisGRMHDadaptive-mesh-refinement
non-spinning BH as in [76]. However, the disk model is (AMR)code[111–113], whichadoptstheCactus/Carpet
6
infrastructure [114–116]. This code has been extensively One spherical coordinate system covers the entire evolu-
tested and used in the past to study numerous sys- tiondomainandiscenteredonthebinarycenterofmass,
tems involving compact objects and/or magnetic fields andtwosmalleronesarecenteredoneachBH.Thesenew
(see e.g. [110, 117–121]), including black hole binaries in gridsemployalogarithmicradialcoordinate. Weusethe
gaseous media [47, 76, 95]. CTS solution stored on these spherical grids to interpo-
The code solves the equations of ideal GRMHD in a latethedataontotheCartesianevolutiongridswhenever
flux-conservativeformulation[seeEqs. (27)-(29)in[112]] we perform the rotation transformation.
employingahigh-resolution-shock-capturingscheme(see We have checked that the mapping from the spectral
[111, 112] for details), and including effective cooling grid to the spherical grids is implemented correctly by
source terms [see Eqs. (65) and (66) in [117]]. To en- performing vacuum simulations that use the CTS solu-
forcethezero-divergenceconstraintofthemagneticfield, tion stored in the spherical grids as initial data. More
we solve the magnetic induction equation using a vec- specifically, we have computed several diagnostic quan-
tor potential formulation [see Eq. (9) in [113]]. As our tities which characterize the BHs and the global space-
EM gauge choice we use the generalized Lorenz gauge time and compared them with the values known from
condition we developed in [76] and used in [118, 121]. the spectral CTS initial data (see Table I). These agree
We choose the damping parameter ξ = 8/M. The to within 1 part in 104. Moreover, we have verified that
advantage of this gauge condition is that it leads to a crude estimate for the orbital frequency of the binary
damped, propagating EM gauge waves preventing spu- (orbital trajectory traverses a full phase of 2π) as deter-
rious magnetic fields from arising near AMR boundaries mined by a dynamical vacuum evolution agrees with the
even more effectively than the standard Lorenz gauge value given by the initial data (MΩ = 0.028) to within
choice (ξ =0) [113]. The damping properties of the gen- ∼ 10%. We have computed the normalized L2,N norm
eralized Lorenz gauge are crucial for stable and accurate of the Hamiltonian and momentum constraint violations
long-term GRMHD evolutions. The “algebraic” gauge as introduced in Eqs. (59) and (60) in [123], with the
condition used in the first GRMHD simulations adopt- modificationthatwesplituptheLaplacianoperatorinto
ing A-field evolution (see e.g. [49, 112, 122]) was shown its individual components when computing normalized
in [113] to suffer from spurious conversion of EM gauge norms. WefindthenormalizednormoftheHamiltonian
modes into physical modes and vice-versa, due to inter- constrainttobedominant,withL2,N(H)∼2%. Wecon-
polation at AMR boundaries. These spurious magnetic clude that the CTS solutions are mapped correctly and
fields contaminate the evolution and the effect is exacer- accurately.
bated when matter crosses refinement boundaries.
To close the system of equations we use a Γ-law EOS
3. Cooling
P =ρ (cid:9)(Γ−1), (6)
0
where P is the total pressure, ρ the rest-mass density, Without cooling, the binary tidal and the viscous
0
and (cid:9) the specific internal energy. This EOS allows for torques act to gradually heat and puff up the disk. “Ad-
shock heating. We choose Γ = 4/3 appropriate for radi- vective” cooling, which is crucial in slim-disk and ADAF
models[50],isself-consistentlyaccountedforinoursimu-
ation pressure-dominated, optically thick disks.
lations. However,addingradiativecoolingmaybeneces-
sarytoachieveasteadystate. Realisticradiativecooling
based on actual physical mechanisms depends on com-
2. Metric evolution
plicated microphysics, which we do not model here, but
intendtoincorporateinfuturestudies. Tomodelsteady-
The metric evolution is treated under the approxima-
state solutions in this work, we introduce “radiative”
tion that the inspiral time scale due to GW emission is
cooling via an effective cooling leakage scheme. This
long compared to both the binary orbital period and the
scheme is, strictly speaking, valid in the optically thin
viscous time scale of the disk. Hence, we can neglect
regime. While this is a very crude approach to radiative
the inspiral for multiple binary orbits. The CTS initial
cooling, treating both extremely rapid radiative cooling
data we adopt possess a helical Killing vector, which im-
as well as no radiative cooling can help bracket the pos-
pliesthatthegravitationalfieldsarestationaryinaframe
sibilities. Different formulae for the cooling emissivity Λ
corotating with the binary. As a result, we can perform
have been proposed in the literature:
the metric evolution in the center-of-mass frame of the
binary by simply rotating the initial spacetime fields as I. Non-exponential cooling [72, 124]:
was done in [47, 76]. This technique simplifies our com-
putations substantially. In addition, the rotating metric Λ= ρ0(cid:9) (cid:2)ΔS +(cid:10)(cid:10)ΔS(cid:10)(cid:10)(cid:3), (7)
method facilitates our evolving dynamically to relaxed τ S S
cool 0 0
disk/magnetic field initial data for the inspiral.
To implement this method, we map the CTS solution where S ≡ K = ρPΓ is an entropy parameter, and
from the spectral grid onto three grids corresponding to S the initial targe0t value. We call this emissiv-
0
three partially overlapping spherical coordinate systems: ity “non-exponential”, because the effective cooling
7
time scale for this scheme is not just τ , but de- The grid spacing is also motivated by both MRI resolu-
cool
pends on the internal energy (cid:9) (see Appendix A). tion requirements and the results gleaned from test runs
involving hydrodynamic disk evolutions around a single,
non-spinning BH, which we report on in Sec. IV.
II. Exponential cooling [117, 119]
In Table II we also list the distinguishing characteris-
Λ= ρ0(cid:9)th = ρ0(cid:9)0 (cid:2)ΔS +(cid:10)(cid:10)ΔS(cid:10)(cid:10)(cid:3), (8) tlaicbselosfatrheecdhiffoseernenttocadseessigwneatceonthsiedemrainsstrhaitsiow,owrkh.etThheer
τ 2τ S S
cool cool 0 0 cooling is applied or not, and the resolution. For ex-
where (cid:9) is the internal energy calculated using S , ample, the label 1:1nc-hr means mass ratio q = 1, no-
0 0
and(cid:9)th =(cid:9)−(cid:9)0. Inthisschemetheeffectivecooling cooling, and high resolution.
time scale is τ . We point out that the disparity in length scales (hori-
cool
zon vs. disk size) and time scales (the Courant condi-
Bothemissivitiesdissipateallshock-inducedthermalen- tion vs. viscous time scale) intrinsic to the circumbinary
ergy, driving the entropy of the gas to its initial value. BHsdiskproblemintroducesalargecomputationalcost.
We use prescription II, instead of I, because we have Mostofoursimulationswereruncontinuously(excluding
found that prescription I is prone to the development of queue waiting times) for ∼ 2 months on Blue Waters,
a Courant instability, as the effective cooling time scale Kraken,Lonestar,aswellastheIllinoisRelativitygroup
of this scheme depends on the amount of shock heating, 36-nodeBeowulfcluster. TheCPUhoursuseddepended
which can be very strong in low-density regions (see e.g. on the computer cluster, but we estimate that the simu-
Appendix B in [125]). Thus, to stabilize the simulations lations presented here required ∼2×106 CPU hours.
with prescription I, one typically excludes cooling of the
low-density regions or unbound matter [72, 76]. As both
the BH horizons and the low-density cavity is included
C. Diagnostics
in our computational domain (unlike earlier studies), we
find that the strong shock heating inside the cavity in
conjunction with emissivity I leads to a numerical insta- Wedistinguishtwotypesofdiagnostics. Thefirsttype
bility. Inordertobrackettheeffectofcooling, thisinner consists of probes of the MHD flow, including the den-
cavity needs to be cooled when cooling is enabled. In sityandvelocityprofiles,accretionrates,luminosityesti-
AppendixA,wepresentananalyticcalculationillustrat- mates, magnetic field profiles, the establishment of MRI
ingtheaboveconsiderations. Wedemonstratethatshock turbulenceandjets,etc. Thesecondtypeconcernsprop-
heatingofmatterinthecavityisimportantinSec. IVC. erties of the plasma such as local temperatures, optical
To model “rapid” radiative cooling we set the cooling depths, characteristic frequencies of emitted radiation,
timescaleequalto10%ofthelocal,Kepleriantimescale etc. The first type are straightforward to calculate from
τ /M =0.1τ /M =0.1·2π(r/M)3/2, where r is the oursimulationdata, asthey dependon theoverallMHD
cool Kep
cylindricalradialcoordinatemeasuredfromthecenterof behavior of the disk, and once we have chosen an EOS
mass of the binary. In order to prevent the cooling time andcoolingprescription, areindependentofdetailedmi-
from becoming prohibitively small as r →0 we floor the crophysics. We are quite careful and confident in re-
cooling time at τ ≥ 10M. Throughout this paper we porting these diagnostic quantities. The second type de-
cool
refer to cases with Λ(cid:9)=0 as the cooling cases and Λ=0 pends crucially on the specific physical values we assign
as the no-cooling cases. to our nondimensional input parameters (e.g. BH mass
and disk density) and to the microphysics that is not
incorporated in our calculation (e.g. realistic radiative
cooling and radiative transport). Nevertheless, we can
make crude estimates for the latter quantities, and do so
4. Evolution grids & models below. Althoughconsiderablecautionmustbeappliedto
theseestimates,theymayserveasusefulguidestosubse-
quent, more detailed investigations and to astronomical
We use a hierarchical, box-in-box adaptive mesh pro-
instrumentsthatmaybeabletoobservethescenarioswe
videdbytheCactus/Carpetinfrastructure[115,116]. We
are simulating.
constructed two sets of nested boxes, with one set cen-
tered on each BH, on which we discretize the GRMHD The first type of diagnostics includes: 1) Accretion
evolution equations. The coarsest level has an outer rate M˙ as defined in [47]. We compute the total accre-
boundary at r = 240M. Due to a range of resolution tionrateontothebinary,andalsotheaccretionrateonto
requirements related to the different sizes of the BHs, each individual BH. 2) Fourier analysis of the accretion
different models use different number of refinement lev- rate FFT(M˙ ), targeted to identify possible quasiperi-
els,whichinturnyieldsdifferentfinestgridspacings(see odic signatures in the accretion flow. 3) Surface den-
Table II). We set up the locations of our AMR bound- sity profile Σ(r) as defined in [95]. This diagnostic is
aries such that the computational grids resolve both the also useful to compare with studies of 1D orbit-averaged
BHs and the inner disk region for the given resources. disk models. 4) Shakura-Sunyaev α-stress parameter
8
TABLE II. List of grid parameters for all models. Equatorial symmetry is imposed in all cases. The computational mesh
consists of two sets of nested AMR grids, one centered on each BH, with the outer boundary at 240M in all cases. From left
torightthecolumnsindicatethecaselabel,massratio q,whethercoolingisincludedornot,thecoarsestgridspacingΔxmax,
number of AMR levels around the primary (BH) and the secondary (bh), and the half length of each AMR box centered on
each BH. The grid spacing of all other levels is Δxmax/2n−1, n=1,2,...,nmax, where n is the level number such that n=1
corresponds to the coarsest level. A dash “–” indicates “no information available”.
Case q cooling? Δxmax levels(BH) levels(bh) Grid hierarchy
1:1nc-hr 1:1 No 4.8M 7 7 240M/2n−1, n=2,...5, 240M/2n, n=6,7
1:1nc-mr 1:1 No 6.0M 7 7 240M/2n−1, n=2,...5, 240M/2n, n=6,7
1:1nc-lr 1:1 No 8.0M 7 7 240M/2n−1, n=2,...5, 240M/2n, n=6,7
1:2nc-lr 1:2 No 8.0M 7 8 240M/2n−1, n=2,...5, 240M/2n, n=6,7,8
1:4nc-lr 1:4 No 8.0M 7 9 240M/2n−1, n=2,...5, 240M/2n, n=6,...,9
1:8nc-lr 1:8 No 8.0M 6 10 240M/2n−1, n=2,...5, 240M/2n, n=6,...,10
1:10nc-lr 1:10 No 8.0M 6 10 240M/2n−1, n=2,...5, 240M/2n, n=6,...,10
0nc-hr 0 No 4.8M 6 – 240M/2n−1, n=2,...5, 240M/2n, n=6
0nc-mr 0 No 6.0M 6 – 240M/2n−1, n=2,...5, 240M/2n, n=6
0nc-lr 0 No 8.0M 6 – 240M/2n−1, n=2,...5, 240M/2n, n=6
1:1c-mr 1:1 Yes 6.0M 7 7 240M/2n−1, n=2,...5, 240M/2n, n=6,7
1:2c-lr 1:2 Yes 8.0M 7 8 240M/2n−1, n=2,...5, 240M/2n, n=6,7,8
1:4c-lr 1:4 Yes 8.0M 7 9 240M/2n−1, n=2,...5, 240M/2n, n=6,...,9
1:8c-lr 1:8 Yes 8.0M 6 10 240M/2n−1, n=2,...5, 240M/2n, n=6,...,10
1:10c-lr 1:10 Yes 8.0M 6 10 240M/2n−1, n=2,...5, 240M/2n, n=6,...,10
0c-lr 0 No 8.0M 6 – 240M/2n−1, n=2,...5, 240M/2n, n=6
computed as α ∼ αEM ≡ (cid:10)(cid:4)T(cid:4)rˆPEφˆM(cid:5)(cid:5)(cid:11)t where TrˆEφˆM is the radiative luminosity Lb =LEM+Lcool.
dominantorthonormalcomponenttotheMaxwellstress- The second type of diagnostics √includes: 10) Opti-
cal depth to true absorption τ∗ = 3τ τ (Eddington
energy tensor evaluated using the tetrad defined in [126] es ff
approximation), where we assume pure hydrogen and
(the brackets denote an orbit averaged quantity). More
specifically, we report an azimuthally- and z- averaged where τes (τff) is the optical depth to electron scatter-
α=α(r)profile, whichcanbeusedin1Dorbit-averaged ing (free-free absorption), calculated using the Thomp-
disk models. 5) Disk scale height H = Σ/ρ (z = 0). son scattering opacity κes = 0.4cm2/g (Rosseland mean
6) Inner disk edge radius r : In all Σ(r)-profil0es we ob- opacity κ¯ff = 6.45 · 1023ρ0T−3.5cm2/g) as τes = κesΣ
serve that inside the cavitiyn Σ(r) declines rapidly with (τff =κffΣ) [127]. τ∗ >1 implies the matter is optically
decreasing r and becomes convex. We fit a fifth or- thicktoabsorption. 11)Localtemperatureofthematter,
der polynomial to the orbit-averaged Σ(r) in the con- calculated by solving (cid:9) = aT4/ρ0 +2kBTρ0/mp, where
vex region at small r, and define r as the radius where mp is the proton mass and kB the Boltzmann constant.
in
the curvature of the Σ(r) fitting function is maximized, 12) In the cases with cooling, the effective disk temper-
[Σ(cid:6)(cid:6)(r)/(1 + Σ(cid:6)(r)2)3/2](cid:6) = 0 where (cid:6) ≡ d/dr. 7) EM ature (in cooling cases), estimated by assuming that all
Poynting luminosity (L ) as defined in Eq. (1) of cooling luminosity is emitted as black body radiation:
EM
[121]. 8) En(cid:11)ergy loss rate carried off by the outflowing
gas Lgas = sT0,r(gas)dS. This surface integral must be Teff =[Lcool/4πσ(ro2ut−ri2n)]1/4, (9)
performed in the asymptotically flat regime. Given that
we do not perform the integration at an infinite radius, whereσ istheStefan-Boltzmannconstantandr isthe
as a crude approximation to L we include in the in- out
gas disk outer radius. 13) Characteristic observed thermal
twehgircahtioant loanrglye rmaadtiiteEr t=ha−tuis0 >unb1.ouWnde,cio.em.,pumtaet7te)rafnodr rthadeiPatlaionnckfrceoqnusetnacnitesanνdbbz=thekcBoTsemffo/lho(g1ic+alzr)e,dwshhifetr.eThhiiss
8) at several radii including 90M,140M,210M. 9) For is calculated only when Λ (cid:9)= 0. 14) Cyclotron emission.
cases where our cooling scheme i(cid:11)s enabled, we compute While we find the bulk of the disk to be optically thick
the cooling luminosity Lcool = sT0,r(rad)(cid:12)dS, wh√ich we near Eddington accretion rates, the low-density “cavity”
estimate via the volume integral Lcool = Λu0α γd3x. is optically thin. From these regions cyclotron lines may
Thevolumeintegrationisexactlyequaltothesurfacein- be detectable. 15) In cases where Λ(cid:9)=0 we compute the
tegration at steady-state and in spacetimes possessing a characteristic cyclotron frequencies ν =eB/mc(1+z),
cy
timelikeKillingvectorandwhenweignoreanyradiation where e is the electron charge, m the electron mass and
captured by the BHs. We also compute the bolometric B the magnitude of the magnetic field.
9
IV. TESTS AND RESOLUTION
REQUIREMENTS
In this section we describe tests we performed that
motivate the choices for the grid resolution and cooling
time scale.
A. Hydrodynamical evolutions with B =0
To set a lower limit on the necessary resolution to
perform our GRMHD simulations, we found the mini-
mum resolution required so that our code maintains sta-
ble equilibrium of an unmagnetized disk around a single
non-spinningBHforseveralthousandsofM ofevolution.
The equilibrium disk solution we use is described at the FIG. 2. Contours of the λMRI-quality factor Q=λMRI/dx in
theequatorialplane,correspondingtotheequalmass(q=1)
beginning of Sec. IIIA2. Our study indicates that for
medium resolution case [divide (multiply) by 1.33 (1.25) for
thelowresolution(Δx =8.0M)listedinTableIIthe
max the low (high) resolution case] at t = 0M. We resolve the
surface density profile of the initial disk is maintained to fastest growing MRI mode by (cid:3) 10 gridpoints over a large
within 2% throughout the disk for at least t∼5000M. part of the disk (the blue ring stems from extremely small
values of λMRI, when the vertical component of the B-field
changes sign).
B. MRI resolution requirements
HerewechecktheconditionsforMRItooperateinour
diskmodels. Forthistobethecasethreeconditionshave
to be satisfied: (I) A magnetic field configuration must
be present that violates the stability condition for MRI
dΩ/dr ≥ 0, (II) The wavelength of the fastest-growing
mode λ has to be resolved by (cid:3) 10 gridpoints [128–
MRI
130], and (III) the B-field must be sufficiently weak for
λ (cid:2)2H. Inotherwordsthewavelengthofthefastest
MRI
growing mode should fit in the disk [131].
Regarding (I) our initial disk model is unstable to the
MRI because of the outwardly decreasing angular veloc-
ity and the presence of an initially small poloidal mag-
netic field. Regarding (II) we plot the quality factor
Q ≡ λ /dx where dx is the local grid spacing (which
MRI
jumpsbyafactoroftwoatAMRboundaries);seeFig.2.
FIG.3. Rest-massdensitycontours(colorcoded)onamerid-
One can see that we satisfy the criterion Q (cid:3) 10 over a ionalslice,andλMRI/2(blacksolidline)att=0M. Theplot
ratherlargeportionofthediskinitiallyexceptforthere- corresponds to equal mass (q=1) but is the same among all
gionneartheradiuswherethepoloidalfieldchangessign cases considered in this work. For the most part λMRI/2 fits
and becomes very small. We chose our low-resolution within disk.
grids such that condition (II) is satisfied at t = 0 in the
innermost parts of the disk.
Regarding(III)weplotameridional(x,z)-sliceofthe C. Cooling
rest-mass density overlayed by a line plot showing λ
MRI
as a function of x in Fig. 3. The plot shows that for We seek to compare two opposite limiting cases for
the most part λ fits within the disk. At the inner each mass ratio: (I) No-cooling, τ (cid:13) τ for which
MRI cool Kep
disk edge the MRI is likely to be suppressed initially, Λ=0. Here τ is the local Keplerian time scale which
Kep
butastheevolutionproceedsmagneticwindingconverts is comparable to the (shock) heating time scale; (II) ex-
poloidalfieldsintotoroidalonesloweringλ andeven- tremelyrapidcooling, τ (cid:14)τ whichwemodelwith
MRI cool Kep
tually triggering the MRI near the inner disk-edge. the effective emissivity Λ= ρ0(cid:5)th.
τ
cool
10
To determine the value for τ at which cooling be- and where ε ≡ L /M˙ = 0.08 is the radiative efficiency
cool b
comes rapid (at least in the bulk of the disk), we per- as computed via our simulations for our adopted cooling
formed the q = 1 cooling case using different cooling lawinthesingle,nonspinningBHcase. Intheno-cooling
time scales. We concluded that rapid cooling requires cases the only radiation luminosity estimate we have is
a cooling time scale significantly shorter than the lo- the Poynting luminosity which is expected to be a small
cal Keplerian time scale. In Fig. 4, we plot K/K , fraction of the total radiative luminosity. Hence, in the
0
where the entropy parameter S ≡ K = P/(ρΓ) and no-cooling cases we do not scale our results with the Ed-
0
S = K = K(t = 0) for a run with τ = 0.1τ dingtonratio. Instead,wechooseafiducialaccretionrate
0 0 cool Kep
(leftpanel),arunwithτ =τ (middlepanel)anda similar to the one given in Eq. (10).
cool Kep
run without cooling τ =∞ (right panel). It becomes
cool
apparentthatwhenτ =τ ,K isnotdrivenbackto
cool Kep
itsinitialvalue. Physically,thismeansthatnotallshock A. Resolution study
generated entropy is radiated away, hence τ = τ
cool Kep
does not correspond to rapid cooling and steady state
Here we present the results of our resolution study.
is not achieved. For τ = 0.1τ we find K/K ∼ 1
cool Kep 0 ForthesingleBH,no-coolingandBHBHequalmass,no-
in the bulk of the disk. The values depart from unity
cooling cases we use the three resolutions (see Table II).
only inside the cavity where low density gas exists and
In the single BH-case the average accretion rate varies
can be shock heated to very high K/K , demonstrating
0 littlewithresolution(seeleftpanelinFig.5). Themaxi-
thatshockheatingisextremelystronginthelow-density
mum fractional difference of the mean accretion rate be-
cavity. We adopt τ =0.1τ in all our cooling simu-
cool Kep tweendifferentresolutionsis15%. Otherquantitiesshow
lationsbecauseitleadstorapidcooling,atleastthrough-
a similar behavior. These results indicate that the res-
out the bulk of the disk.
olutions used in this case are high enough for the main
MRI effects to be captured and the results to be qualita-
tivelyindependentofresolution. However,theresolution
V. RESULTS
is not sufficiently high to perform a formal convergence
test.
In this section we present the results of our numerical
For the equal-mass case we observe a different behav-
simulations. First, in Sec. VA we show results from our
ior. The mean accretion rates appear to converge (see
resolution study. In Sec. VB, we directly compare the
right panel in Fig. 5), but the low resolution run under-
q =1 binary case with B (cid:9)=0 to the B =0 case. Lastly,
estimatestheaccretionratebyalmostafactorof2,while
we report on the variation of our diagnostics with mass
the medium and high resolution runs are in good agree-
ratio for all magnetized cases in Sec. VC.
ment. For the latter resolutions mean accretion rates
Our simulations have two parameters that scale out agree to within 9%. We conclude that for the q =1 case
of the problem: the total mass of the BHBH binary M
our adopted medium and high resolutions are sufficient
and the rest-mass density of the disk. Alternatively, we
for drawing qualitative conclusions and that higher res-
can exchange one of these parameters with another pa-
olutions are necessary for accurate quantitative results
rameter that depends on these two free parameters. So,
that reside in the convergent regime.
instead of the rest-mass density, in the results we quote The results of the resolution study in the q = 0 and
we choose theEddington ratioLb/LEdd, where Lb is the q =1casesdifferbecauseofthedistributionofmatterin
bolometric EM luminosity described in Sec. IIIC. The both cases and our grid setup. In the q =0 case there is
relation between these parameters is determined as fol-
more matter close to the BH, where very high resolution
lows: the accretion rate through the horizon must scale refinement levels reside, whereas in the q = 1 case the
like M˙ ∝ ρ0,refM2, where ρ0,ref is a reference rest-mass bulk of the matter remains outside the inner edge of the
density in the disk. We choose the maximum rest-mass
disk,wherethegridresolutionisnotashigh. Asweshow
density at t=0 as the reference density, and our simula- later, this is not the case in our q < 1 models. There,
tions determine the proportionality constant. For exam-
more matter resides closer to the BHs, and hence closer
ple, in the single, nonspinning BH case with cooling we
to the high resolution levels.
find
Based on our resolution study we conclude that the
(cid:6) (cid:7)(cid:6) (cid:7)
L M −1 low resolution used in the equal-mass case is not suffi-
ρref =9×10−12 LEbdd 108M(cid:2) g cm−3, ciently high to yield reliable results, but for the other
mass ratios we consider in this work, our adopted low
wherewehavereplacedtheaccretionrateviathefollow- resolution suffices for a qualitative discussion. Thus, in
ing expression theequal-masscasesresultswillbereportedmainlyfrom
(cid:6) (cid:7) our medium resolution runs.
L
M˙ = L ε−1 = b L ε−1 We stress here that these simulations, which include
b L Edd
(cid:6) (cid:7)E(cid:6)dd (cid:7) (10) all relativistic effects and resolve the BHs, are computa-
L M tionallyverydemanding(seeSec.III)andincreasingthe
≈ 2.75 LEbdd 108M(cid:2) M(cid:2) yr−1, resolution even further will incur a very large computa-