Accretion disks around binary black holes of unequal mass:
GRMHD simulations near decoupling
Abstract
We report on simulations in general relativity of magnetized disks onto black hole binaries. We vary the binary mass ratio from 1:1 to 1:10 and evolve the systems when they orbit near the binarydisk decoupling radius. We compare (surface) density profiles, accretion rates (relative to a single, nonspinning black hole), variability, effective stress levels and luminosities as functions of the mass ratio. We treat the disks in two limiting regimes: rapid radiative cooling and no radiative cooling. The magnetic field lines clearly reveal jets emerging from both black hole horizons and merging into 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 and luminosities. Accretion streams near the horizons are among the densest structures; in fact, the 1:10 nocooling evolution results in a refilling of the cavity. The typical effective temperature in the bulk of the disk is yielding characteristic thermal frequencies . These systems are thus promising targets for many extragalactic optical surveys, such as LSST, WFIRST, and PanSTARRS.
pacs:
04.25.D, 04.25.dg, 47.75.+fI Introduction and Motivation
Supermassive black hole (SMBH) binaries can form in magnetized plasma following galaxy mergers, via barmode instability in rapidly rotating supermassive stars, or by other dynamical processes Begelman:1980vb (). After formation, a combination of dynamical friction and gasdriven migration is likely to catalyze the binary inspiral into the gravitational radiationdriven regime Ivanov:1998qk (); Haiman:2009te (); Rafikov:2012hd (); Begelman:1980vb (); Armitage:2002uu (); Armitage:2005xq (); Milosavljevic:2001vi (). The exact details of these processes, including the “lastparsec problem”, remain active areas of research (see e.g. Milosavljevic:2002ht (); Merritt:2004gc (); Khan:2013wbx () and references therein). As a result, event rates and population synthesis studies at this stage are highly uncertain Dotti:2011um ().
The exciting prospect of a simultaneous observation of both electromagnetic (EM) and gravitational waves (GWs) arising from accreting binary BHBHs makes these systems prime targets in the era of multimessenger astronomy. Such observations will enable us to determine the binary masses, BH spins, redshift and even determine the Hubble constant to better than 1% Schutz:1986gp (); Holz:2005df (); Nissanke13 ().
Gravitational waves (GWs) from SMBHs are expected to be detected by planned GW interferometers such as eLISA/NGO detectors AmaroSeoane:2012km (), sensitive to GW frequencies , and the currently operating Pulsar Timing Arrays Hobbs:2009yy (); Tanaka:2011af (); Sesana:2012ak (), sensitive to frequencies . As SMBHs are typically believed to have masses in the range , the GW strain at orbital separation  the value adopted in this work  is . Here is the total mass of the binary, and we normalized to a luminosity distance corresponding to redshift in a standard CDM Universe. The corresponding gravitationalwave frequency is Hz. The expected GW strain is above the eLISA strain sensitivity at these frequencies AmaroSeoane:2012km (), hence these systems will be detectable by eLISA. In particular, for , the value of at falls well within the eLISA sensitivity band even for larger redshifts, while for and large redshifts (), the value of at is marginally within the eLISA sensitivity band. However, the inspiral time from these separations is assuming equal masses. Hence, these systems will quickly enter the eLISA sweet spot, and EM precursor signals can trigger targeted GW searches with a convenient lead time of several days.
While awaiting the first detection of GWs, currently operating and future electromagnetic (EM) detectors such as LSST Abell:2009aa (), WFIRST Green:2012mj () and PanSTARRS Kaiser:2002zz () are promising instruments to identify accreting BHBH systems in the EM spectrum. Important steps have already been made toward realizing this goal.
Currently, we know of one spatially resolved SMBH binary candidate at an orbital separation pc: 0402+379 Rodriguez:2006th (). Other spatially unresolved, SMBH binary candidates have been found, including OJ287 Sillanpaa88 (); lehto96 (); Valtonen:2011ny (); Tanaka:2013cva () and SDSS J1536+0441 Boroson:2009va (); Chornock09 (). Binary AGN candidates have been singled out based on offsets in the broad line and narrow line regions, emission line profiles, and time variability Decarli:2013xwa (); Eracleous:2011ua (). Recently, verylong baseline interferometric observations interpreted ejection components from AGN cores as undulations caused by the precession of the accretion disks around a SMBH binary Roland:2013jda (). A simplified model was applied to two AGN sources; for 1823+568 their analysis yields pc, and a mass ratio in the range , while for 3C 279 pc, . However, given the lack of a robust circumbinary accretion disk theory these results are at best preliminary. Nevertheless, they still motivate a study of accretion flows onto supermassive BH binaries with different mass ratios, which we initiate in this work in general relativity (GR).
Other features identified as “smoking guns” for binary black holes include BH recoil/kicks Komossa:2012cy (), used to explain the large velocity offsets between emission lines in AGN spectra, as well as observed kinks in jets probably due to changes in BH spin (Xshaped radio sources), a past merger event, or precession effects Romero2000 (); Roos1993 (); Sudou:2003hv (). Modifications to the line profiles have also been proposed as promising characteristic features to distinguish binary BH AGN from classical, single BH AGN sources McKernan:2013cha ().
To assist and solidify all these detection efforts, it is crucial to identify and model possible electromagnetic “precursor” and “afterglow” signatures Milosavljevic:2004cg (); Kocsis:2007yu (); Haiman:2008zy (); Tanaka:2009iy (); Shapiro:2009uy (); Liu:2010mh (); Shapiro:2013qsa ().
Depending on the physical regime the properties characterizing the gas can differ considerably, and different accretion models are applicable. For example, if the gas has little angular momentum, the accretion flow resembles the binary analog of a BondiHoyleLyttleton solution bondi44 (); bondi52 (); 1989ApJ…336..313P () (see Farris:2009mt (); Zanotti:2011mb (); Giacomazzo:2012iv () for GR studies). If the gas has significant angular momentum, then the gas can become rotationally supported and form a disk.
For a BHBH embedded in a disk, one can identify several different regimes based on the time scales for migration of the binary. For a SMBH binary engulfed by a (thin) disk at large separation, the migration of the binary is initially governed by binarydisk angular momentum exchange mediated by (effective) viscosity Haiman:2009te (). At large enough separations, the reduced mass of the binary is less than the local interior disk mass ( where is the surface density of the disk). This leads to the socalled diskdominated type II migration occurring on the viscous time scale . As the migration proceeds, the reduced mass of the system becomes larger than the local disk mass and the migration enters the secondarydominated type II regime, which occurs on a longer time scale , where a constant of order if and otherwise. Ultimately, the binary enters a regime at smaller orbital separations where angular momentum losses due to GWs dominate, and the binary migrates on the GW time scale . In all regimes, the disk moves inwards on the viscous time scale. When or the disk can follow the inspiral of the binary and settle in a quasisteady state  this is called the predecoupling regime. When , the binary decouples from the disk, i.e. the inward migration of the binary outpaces the inward drift of the disk.
In this paper we focus on the phase near decoupling, while the postdecoupling regime will be the subject of a future paper. We note that unlike eLISA, Pulsar Timing Arrays are sensitive only to SMBH binaries well within the predecoupling regime.
Accretion onto a single BH has been studied in great detail for decades, and magnetohydrodynamic studies in GR have drastically improved our understanding of these flows (see lrr20131 () for a recent review). Many different disk models have been proposed in the literature. These models range from geometricallythin, optically thick disks Shakura73 (); Novikov73 () and slim disks Abramowicz1988 (), to geometrically thick, optically thin, radiatively inefficient accretion flows (RIAF) Ichimaru1977 (); Blandford:1998qn (); Igumenshchev:1997vi (); Esin:1997he (); Narayan:1996gp (). However, our understanding of accretion flows onto BHBHs remains poor and studies of these systems are still in their infancy.
The first analytic Newtonian model and smooth particle hydrodynamic simulation of a circumbinary accretion disk was given in Artymowicz:1994bw (). Since then, other Newtonian (semi)analytic studies Tanaka:2009iy (); Shapiro:2009uy (); Liu:2010mh (); Kocsis:2012ui (); Kocsis:2012cs () and hydrodynamic simulations in 2D Artymowicz:1996zz (); MacFadyen:2006jx (); D'Orazio:2012nz (); Farris:2013uqa (), and 3D Dotti:2006ef (); Cuadra:2008xn (); Rodig:2011jz (); Roedig:2012nc (); Hayasaki:2012wu () have followed. Newtonian magnetohydrodynamic (MHD) simulations were presented in Shi:2011us (), and PostNewtonian MHD simulations in Noble:2012xz (). Many of these earlier studies excluded the binary and most of the inner cavity from the computational domain, introducing an artificial inner boundary condition. The importance of treating the inner regions selfconsistently has been discussed in Roedig:2012nc (), and only full GR calculations can achieve this goal reliably.
A “GRhybrid” orbitaveraged model for thin disks, in which the viscous part is handled in GR and the tidal torques in Newtonian gravity was introduced in Shapiro:2013qsa (), and GR hydrodynamical simulations of accretion onto BHBH binaries  taking into account the dynamical spacetime  have been performed in Bode:2009mt (); Farris:2009mt (); Bogdanovic:2010he (); Bode:2011tq (). To date the only GR magnetohydrodynamic (GRMHD) simulations of disk accretion onto BHBHs that account both for the dynamical spacetime and the BH horizons were presented in Farris:2012ux ().
Using a different approach by assuming that the Bfield is anchored to a circumbinary disk outside the computational domain, Mosta:2009rr (); Neilsen:2010ax (); Palenzuela:2010xn (); Palenzuela:2010nf () modeled EM signatures by solving the GR forcefree electrodynamic equations.
Close to merger a single BH remnant is formed on a time scale much shorter than the dynamical time scale in the disk. The mass and angular momentum of the remnant BH is different from the total mass and angular momentum prior to merger due to GW emission, causing a quasiinstantaneous perturbation to the disk. This effect has been modeled using hydrodynamical and MHD simulations in Newtonian gravity and in GR Corrales:2009nv (); Rossi:2009nk (); Anderson:2009fa (); Megevand:2009yx (); Ponce:2011kv ().
A realistic and ideal 3D global model for a circumbinary disk around a SMBH binary near the decoupling radius requires: a) a fully relativistic treatment of gravitation in a dynamical spacetime, b) GRMHD for the plasma flow, c) realistic cooling processes, and d) radiative transfer in curved spacetime. Simulations incorporating these effects must also have high resolution and long integration times (several viscous time scales). However, including all these ingredients in one simulation would make these computations prohibitive, because the wallclock times required to integrate for even viscous time scales at the inner disk radius at high resolution are far too long with current supercomputer resources. Thus, some of these ingredients must be relaxed in order to obtain a qualitative understanding of these systems. For this reason, the models in this work feature a) and b). High resolution is only adopted for a few models. In addition, we model radiative cooling by a simple cooling function and consider the extreme opposite limits of “rapid” cooling and “no cooling” to bracket the possibilities.
In this paper we study the effects of the binary mass ratio on the disk near decoupling. We vary the BHBH binary mass ratio , considering 1:1, 1:2, 1:4, 1:8 and 1:10 mass ratios. The mass ratio regime is shown to be of high astrophysical relevance in Gergely:2007ny (); Sesana:2007sh (). Also, it has been argued that BHBHs forming in major galaxy mergers will typically have mass ratios in the range Sesana:2011zv (). These results motivate our choice of mass ratios. Note that Newtonian simulations studying the effects of mass ratio in the context of 2D thindisk models were performed in D'Orazio:2012nz (); Farris:2013uqa ().
The new aspects of our work are the inclusion of relativistic gravitation to resolve the crucial physics near the BH horizons, effective viscosity arising from the magnetorotational instability Balbus:1991ay (); Balbus:1998ja (), geometrically thick disks, three spatial dimensions (3D), a law equation of state (EOS), and effective cooling. We model the decoupling epoch by fixing the binary separation , evolving the spacetime by rotating our conformalthinsandwich initial data Pfeiffer:2002iy (); Cook:2004kt (); Caudill:2006hw (); BSBook () as we have done in Farris:2011vx (); Farris:2012ux (). The dependence of the flow variability, EM signatures, the magnetic field structure and the matter dynamics inside the lowdensity “hollow” on the mass ratio are investigated. We estimate thermodynamic properties of the gas and scale our results with the binary total mass and the luminosity in units of the Eddington luminosity where feasible. From this, emission characteristics including typical thermal frequencies and luminosities are given and relevant detectors are discussed.
In the absence of any observational constraints on the thermodynamic state of accreting BHBHs, the models we consider in this work adopt an adiabatic EOS governed by an adiabatic index , appropriate for radiation pressuredominated, optically thick disks. This choice is motivated in part both by theory and observations. First, accretion disk theory lrr20131 () predicts that the inner part of circumbinary disks around SMBH binaries are optically thick, radiation pressure dominated for a large set of possible disk parameters (see e.g. Novikov73 (); Rafikov:2012hd () and the next section). Furthermore, as discussed in Haiman:2009te (), steadystate disk models predict that radiation pressuredominated circumbinary disks will channel more material into the cavity for a given central mass. This finding makes radiation pressuredominated circumbinary disks promising sources for electromagnetic counterparts.
Second, AGN data from the Sloan Digital Sky Survey (SDSS) Heckman:2004zf (); Shen:2011ge (); Kelly:2012vz (), the AGN and Galaxy Evolution Survey (AGES) Kollmeier:2005cw (), XMMCOSMOS Lusso:2012yv () and others Schulze:2010wc () covering mostly type I AGNs and the local Universe to redshift of , reveal Eddington ratio distributions in the range with the tendency that higher values occur at higher (and possibly smaller central mass ). A comparison of accretion time scales with the age of the Universe suggests that earlier accretion episodes are closer to the Eddington limit, similar to the recently discovered quasar at with a central mass of Mortlock:2011va () accreting at . These surveys therefore motivate studies of disks accreting near the Eddington luminosity, for which radiation pressure is important. Moreover, radiation pressuredominated disks accreting near the Eddington limit are more likely to be detectable, even at large cosmological redshifts.
This paper is structured as follows: In Sec. II we present a qualitative discussion of the range of parameters and the associated physical regimes for which our simulations are appropriate. In Sec. III we describe our methods and techniques for evolving the spacetime, fluid, and magnetic fields, as well as our simple cooling prescription. We also present the different cases we consider and list our diagnostics for characterizing the accretion flow and EM signatures. In Sec. IV we show several tests we performed to motivate our numerical setup. In Sec. V we present the results from our simulations. We conclude in Sec. VI with a summary of the main results and a discussion of future work.
Here and throughout we adopt geometrized units, where , unless otherwise stated.
Ii Qualitative considerations
In this section we use the ShakuraSunyaev/NovikovThorne thindisk model Shakura73 (); Novikov73 () to make rough analytic estimates regarding the physical regime our disk models probe. While the ShakuraSunyaev/NovikovThorne model strictly applies to a viscous flow onto a single BH neglecting the binary tidal torques, we expect that our qualitative analysis will apply roughly to our circumbinary disk models, which have , where is the disk scale height. This expectation should be best realized outside the binary orbit because the ratio of the tidal to viscous torques decays quickly far away from the binary orbit (see e.g. Fig. 6 in Shapiro:2013qsa ()).
ii.1 Radiation pressure dominance
ii.1.1 ShakuraSunyaev/NovikovThorne model
The simulations reported here apply to any total (ADM) binary black hole mass and, since we neglect the selfgravity of the gas, to any restmass density , provided it has the same initial disk profile adopted here. The quantities that are fixed, in addition to the initial disk profile, are the adiabatic index appearing in the adopted ideal gas EOS, the ratio of the initial magnetictototal disk pressure, the initial magnetic field profile, the initial binary separation and the cooling law. In this section we use the steadystate thindisk solution for a ShakuraSunyaev/NovikovThorne disk about a single BH of mass to estimate the physical values of some of the gas dynamical quantities as functions of and luminosity . We show below that for a range of astrophysically relevant choices of these parameters the disk is thermal radiation pressuredominated, and this fact motivates our setting .
Neglecting the perturbative role of the secondary tidal torque, the steadystate accretion flow in a geometrically thin disk is uniquely specified by the ShakuraSunyaev viscosity parameter , the central mass and the accretion rate . The quantity is specified in turn by the disk luminosity and efficiency . The ShakuraSunyaev/NovikovThorne model Shakura73 (); Novikov73 () describes a disk that is radiation pressuredominated inside a radius . In this region the radiation to gas pressure ratio is lrr20131 ()
(1)  
where is the radiative efficiency, , and we chose the normalization for the parameter based on typical values found in our simulations. For a thindisk model the efficiency ranges from for a nonspinning BH to for a maximally spinning BH, so we scale to a value residing between these limits. The size of the inner region is determined by the condition , for which Eq. (1) yields
(2)  
The geometrically thick disks we evolve in this work extend radially out to . Hence, when scaling the accretion rate such that our models accrete near the Eddington limit, our models are fully immersed in this inner radiation pressuredominated region.
In this region the typical rest mass densities are lrr20131 ()
(4)  
As we will show later, these characteristic values for and are comparable to the values found in our simulations.
The ShakuraSunyaev/NovikovThorne model predicts that the effective optical depth to absorption is adopting the same normalizations as in the equation above. This wellknown inconsistency near the Eddington limit of the NovikovThorne model is removed with the generalization to the slim disk model Abramowicz1988 (). This model differs from a thin disk in that it allows for cooling to occur via advection, which dominates radiative cooling at high accretion rates (), thereby puffing up the disk. When scaling our models to accrete near the Eddington limit, they are closer to slimdisk models, which remain optically thick in this high luminosity limit Abramowicz1988 (); lrr20131 ().
As we discuss later, near the Eddington limit and , the effective optical depth satisfies in the bulk of our disk models, implying that the gas is optically thick to absorption and the photons eventually thermalize. Thus, these qualitative considerations motivate the adoption of an adiabatic index , appropriate for thermal radiation pressure dominance.
Note also that alternative disk models have been proposed, e.g. Rafikov:2012hd (). However, they largely share the prediction that the inner regions of the disk are radiation pressure dominated.
ii.1.2 Decoupling radius
We estimate the decoupling radius by equating and solve for the separation to find Farris:2012ux ()
(5) 
where we assumed that the inner disk edge radius settles to as typically found in our simulations, is the ShakuraSunyaev turbulent viscosity parameter, and (see also Armitage:2002uu ()). Notice that in contrast to geometrically thick disks, where the decoupling radius is a few tens of , geometrically thin disks have . The decoupling radius estimate (5) for the mass ratios considered in this work ranges from to . The normalizations in Eq. (5) are based on typical values obtained from our simulations. We choose the binary separation for all cases studied in this work, a value which is consistent with the crude estimate of Eq. (5). In the future we intend to start our evolutions at larger separations in order to dynamically determine the decoupling radius and evolve through it as in Shapiro:2013qsa ().
Iii Methods
The models we adopt here assume: a) circular binary orbits, neglecting the binary inspiral (justified for large separations; see Sec. II.1.2), b) non selfgravitating disks, which likely is a good assumption (see, e.g. Goodman:2002gv ()), c) ideal MHD, d) no radiative feedback, e) an effective cooling scheme that brackets no cooling and rapid cooling.
Some of these assumptions may not be obeyed strictly, e.g. the binary may become eccentric in the predecoupling regime Armitage:2005xq (); Dotti:2006ef (); Cuadra:2008xn (); Roedig:2011rn (); Rodig:2011jz (); Roedig:2012nc (), or radiative feedback may become important near Eddington accretion rates. However, our simulations still provide a qualitative understanding of the physics that will be useful in designing the next set of more realistic models of binary black holes immersed in circumbinary disks. In this section, we describe the initial data and computational methods we adopt to account for a)e).
iii.1 Initial data and AMR hierarchy
iii.1.1 Metric initial data
At large separation the binary inspiral time scale is much longer than the binary orbital period and the viscous time scale at the inner edge of the disk just beyond the binary orbit. Accordingly, the inspiral can be neglected over many orbital periods. To model this epoch in GR, we adopt quasiequilibrium conformalthinsandwich (CTS) solutions for the black hole binary Pfeiffer:2002iy (); Cook:2004kt (); Caudill:2006hw (); BSBook (). The spacetime initial data satisfying the CTS equations correspond to a circular orbit and possess a helical Killing vector. The CTS initial data have been generated using the spectral techniques described in Pfeiffer:2002wt () as implemented in the Spectral Einstein Code (SpEC) SpECwebsite () (see also Mroue:2013xna ()). We list the initial data parameters describing our spacetimes in Table 1.



0.98989  0.96865  0.50000  0.50000  0.42958  0.42958  
0.99097  0.85933  0.66667  0.33333  0.60192  0.27312  
0.99342  0.61603  0.80000  0.20000  0.75140  0.15832  
0.99589  0.37868  0.88889  0.11111  0.85640  0.08618  
0.99656  0.31652  0.90909  0.09091  0.88081  0.07022 
iii.1.2 Matter and Bfield initial data
For the fluid we use the same family of equilibrium disk models around single BHs as in Chakrabarti:1985 (); DeVilliers:2003gr (); Farris:2011vx (). We choose the initial inner disk edge radius and specific angular momentum around a nonspinning BH as in Farris:2012ux (). However, the disk model is not identical to Farris:2012ux () due to the different polytropic index ( here versus in Farris:2012ux ()).
We seed the initial disk with a small, purely poloidal Bfield using the same procedure as in Farris:2012ux (); Etienne:2011ea (); see Fig. 1. The field is dynamically unimportant initially: the initial maximum value for the ratio of magnetic to total pressure is 0.025. All cases we consider in this work are initialized with the same disk and magnetic field.
Although we will qualitatively discuss how the evolution of these matter initial data depends on the binary mass ratio, it is important to stress that the goal of our work is to assess how the final relaxed state of the disk depends on the binary mass ratio. The initial data, which governs the early evolution, has no physical significance.
iii.2 Evolution equations and techniques
iii.2.1 GRMHD evolution
We use the Illinois GRMHD adaptivemeshrefinement (AMR) code Duez:2005sf (); Etienne:2010ui (); Etienne:2011re (), which adopts the Cactus/Carpet infrastructure Goodale2002a (); Carpetcactusweb (); Carpetcarpetweb (). This code has been extensively tested and used in the past to study numerous systems involving compact objects and/or magnetic fields (see e.g. Paschalidis:2011ez (); Etienne:2011ea (); Etienne:2012te (); Paschalidis:2012ff (); Etienne:2013qia (); Paschalidis:2013jsa ()), including black hole binaries in gaseous media Farris:2009mt (); Farris:2011vx (); Farris:2012ux ().
The code solves the equations of ideal GRMHD in a fluxconservative formulation [see Eqs. (27)(29) in Etienne:2010ui ()] employing a highresolutionshockcapturing scheme (see Duez:2005sf (); Etienne:2010ui () for details), and including effective cooling source terms [see Eqs. (65) and (66) in Paschalidis:2011ez ()]. To enforce the zerodivergence constraint of the magnetic field, we solve the magnetic induction equation using a vector potential formulation [see Eq. (9) in Etienne:2011re ()]. As our EM gauge choice we use the generalized Lorenz gauge condition we developed in Farris:2012ux () and used in Etienne:2012te (); Paschalidis:2013jsa (). We choose the damping parameter . The advantage of this gauge condition is that it leads to damped, propagating EM gauge waves preventing spurious magnetic fields from arising near AMR boundaries even more effectively than the standard Lorenz gauge choice () Etienne:2011re (). The damping properties of the generalized Lorenz gauge are crucial for stable and accurate longterm GRMHD evolutions. The “algebraic” gauge condition used in the first GRMHD simulations adopting Afield evolution (see e.g. Etienne:2010ui (); Giacomazzo:2012iv (); Rezzolla:2011da ()) was shown in Etienne:2011re () to suffer from spurious conversion of EM gauge modes into physical modes and viceversa, due to interpolation at AMR boundaries. These spurious magnetic fields contaminate the evolution and the effect is exacerbated when matter crosses refinement boundaries.
To close the system of equations we use a law EOS
(6) 
where is the total pressure, the restmass density, and the specific internal energy. This EOS allows for shock heating. We choose appropriate for radiation pressuredominated, optically thick disks.
iii.2.2 Metric evolution
The metric evolution is treated under the approximation that the inspiral time scale due to GW emission is long compared to both the binary orbital period and the viscous time scale of the disk. Hence, we can neglect the inspiral for multiple binary orbits. The CTS initial data we adopt possess a helical Killing vector, which implies that the gravitational fields are stationary in a frame corotating with the binary. As a result, we can perform the metric evolution in the centerofmass frame of the binary by simply rotating the initial spacetime fields as was done in Farris:2009mt (); Farris:2012ux (). This technique simplifies our computations substantially. In addition, the rotating metric method facilitates our evolving dynamically to relaxed disk/magnetic field initial data for the inspiral.
To implement this method, we map the CTS solution from the spectral grid onto three grids corresponding to three partially overlapping spherical coordinate systems: One spherical coordinate system covers the entire evolution domain and is centered on the binary center of mass, and two smaller ones are centered on each BH. These new grids employ a logarithmic radial coordinate. We use the CTS solution stored on these spherical grids to interpolate the data onto the Cartesian evolution grids whenever we perform the rotation transformation.
We have checked that the mapping from the spectral grid to the spherical grids is implemented correctly by performing vacuum simulations that use the CTS solution stored in the spherical grids as initial data. More specifically, we have computed several diagnostic quantities which characterize the BHs and the global spacetime and compared them with the values known from the spectral CTS initial data (see Table 1). These agree to within 1 part in . Moreover, we have verified that a crude estimate for the orbital frequency of the binary (orbital trajectory traverses a full phase of ) as determined by a dynamical vacuum evolution agrees with the value given by the initial data () to within . We have computed the normalized norm of the Hamiltonian and momentum constraint violations as introduced in Eqs. (59) and (60) in Duez:2002bn (), with the modification that we split up the Laplacian operator into its individual components when computing normalized norms. We find the normalized norm of the Hamiltonian constraint to be dominant, with . We conclude that the CTS solutions are mapped correctly and accurately.
iii.2.3 Cooling
Without cooling, the binary tidal and the viscous torques act to gradually heat and puff up the disk. “Advective” cooling, which is crucial in slimdisk and ADAF models lrr20131 (), is selfconsistently accounted for in our simulations. However, adding radiative cooling may be necessary to achieve a steady state. Realistic radiative cooling based on actual physical mechanisms depends on complicated microphysics, which we do not model here, but intend to incorporate in future studies. To model steadystate solutions in this work, we introduce “radiative” cooling via an effective cooling leakage scheme. This scheme is, strictly speaking, valid in the optically thin regime. While this is a very crude approach to radiative cooling, treating both extremely rapid radiative cooling as well as no radiative cooling can help bracket the possibilities. Different formulae for the cooling emissivity have been proposed in the literature:

Nonexponential cooling Shafee:2008mm (); Noble:2012xz ():
(7) where is an entropy parameter, and the initial target value. We call this emissivity “nonexponential”, because the effective cooling time scale for this scheme is not just , but depends on the internal energy (see Appendix A).

Exponential cooling Paschalidis:2011ez (); Paschalidis:2012ff ()
(8) where is the internal energy calculated using , and . In this scheme the effective cooling time scale is .
Both emissivities dissipate all shockinduced thermal energy, driving the entropy of the gas to its initial value. We use prescription II, instead of I, because we have found that prescription I is prone to the development of a Courant instability, as the effective cooling time scale of this scheme depends on the amount of shock heating, which can be very strong in lowdensity regions (see e.g. Appendix B in Etienne:2008re ()). Thus, to stabilize the simulations with prescription I, one typically excludes cooling of the lowdensity regions or unbound matter Noble:2012xz (); Farris:2012ux (). As both the BH horizons and the lowdensity cavity is included 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 instability. In order to bracket the effect of cooling, this inner cavity needs to be cooled when cooling is enabled. In Appendix A, we present an analytic calculation illustrating the above considerations. We demonstrate that shock heating of matter in the cavity is important in Sec. IV.3.
To model “rapid” radiative cooling we set the cooling time scale equal to 10% of the local, Keplerian time scale , where is the cylindrical radial coordinate measured from the center of mass of the binary. In order to prevent the cooling time from becoming prohibitively small as we floor the cooling time at . Throughout this paper we refer to cases with as the cooling cases and as the nocooling cases.
Case  cooling?  levels(BH)  levels(bh)  Grid hierarchy  

1:1nchr  1:1  No  7  7  ,  
1:1ncmr  1:1  No  7  7  ,  
1:1nclr  1:1  No  7  7  ,  
1:2nclr  1:2  No  7  8  ,  
1:4nclr  1:4  No  7  9  ,  
1:8nclr  1:8  No  6  10  ,  
1:10nclr  1:10  No  6  10  ,  
0nchr  0  No  6  –  ,  
0ncmr  0  No  6  –  ,  
0nclr  0  No  6  –  ,  
1:1cmr  1:1  Yes  7  7  ,  
1:2clr  1:2  Yes  7  8  ,  
1:4clr  1:4  Yes  7  9  ,  
1:8clr  1:8  Yes  6  10  ,  
1:10clr  1:10  Yes  6  10  ,  
0clr  0  No  6  –  , 
iii.2.4 Evolution grids & models
We use a hierarchical, boxinbox adaptive mesh provided by the Cactus/Carpet infrastructure Carpetcactusweb (); Carpetcarpetweb (). We constructed two sets of nested boxes, with one set centered on each BH, on which we discretize the GRMHD evolution equations. The coarsest level has an outer boundary at . Due to a range of resolution requirements related to the different sizes of the BHs, different models use different number of refinement levels, which in turn yields different finest grid spacings (see Table 2). We set up the locations of our AMR boundaries such that the computational grids resolve both the BHs and the inner disk region for the given resources. The grid spacing is also motivated by both MRI resolution requirements and the results gleaned from test runs involving hydrodynamic disk evolutions around a single, nonspinning BH, which we report on in Sec. IV.
In Table 2 we also list the distinguishing characteristics of the different cases we consider in this work. The labels are chosen to designate the mass ratio, whether cooling is applied or not, and the resolution. For example, the label 1:1nchr means mass ratio , nocooling, and high resolution.
We point out that the disparity in length scales (horizon vs. disk size) and time scales (the Courant condition vs. viscous time scale) intrinsic to the circumbinary BHs disk problem introduces a large computational cost. Most of our simulations were run continuously (excluding queue waiting times) for months on Blue Waters, Kraken, Lonestar, as well as the Illinois Relativity group 36node Beowulf cluster. The CPU hours used depended on the computer cluster, but we estimate that the simulations presented here required CPU hours.
iii.3 Diagnostics
We distinguish two types of diagnostics. The first type consists of probes of the MHD flow, including the density and velocity profiles, accretion rates, luminosity estimates, magnetic field profiles, the establishment of MRI turbulence and jets, etc. The second type concerns properties of the plasma such as local temperatures, optical depths, characteristic frequencies of emitted radiation, etc. The first type are straightforward to calculate from our simulation data, as they depend on the overall MHD behavior of the disk, and once we have chosen an EOS and cooling prescription, are independent of detailed microphysics. We are quite careful and confident in reporting these diagnostic quantities. The second type depends crucially on the specific physical values we assign 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 below. Although considerable caution must be applied to these estimates, they may serve as useful guides to subsequent, more detailed investigations and to astronomical instruments that may be able to observe the scenarios we are simulating.
The first type of diagnostics includes: 1) Accretion rate as defined in Farris:2009mt (). We compute the total accretion rate onto the binary, and also the accretion rate onto each individual BH. 2) Fourier analysis of the accretion rate , targeted to identify possible quasiperiodic signatures in the accretion flow. 3) Surface density profile as defined in Farris:2011vx (). This diagnostic is also useful to compare with studies of 1D orbitaveraged disk models. 4) ShakuraSunyaev stress parameter computed as where is the dominant orthonormal component to the Maxwell stressenergy tensor evaluated using the tetrad defined in Penna:2010hu () (the brackets denote an orbit averaged quantity). More specifically, we report an azimuthally and  averaged profile, which can be used in 1D orbitaveraged disk models. 5) Disk scale height . 6) Inner disk edge radius : In all profiles we observe that inside the cavity declines rapidly with decreasing and becomes convex. We fit a fifth order polynomial to the orbitaveraged in the convex region at small , and define as the radius where the curvature of the fitting function is maximized, where . 7) EM Poynting luminosity () as defined in Eq. (1) of Paschalidis:2013jsa (). 8) Energy loss rate carried off by the outflowing gas . This surface integral must be performed in the asymptotically flat regime. Given that we do not perform the integration at an infinite radius, as a crude approximation to we include in the integration only matter that is unbound, i.e., matter for which at large radii . We compute 7) and 8) at several radii including . 9) For cases where our cooling scheme is enabled, we compute the cooling luminosity , which we estimate via the volume integral . The volume integration is exactly equal to the surface integration at steadystate and in spacetimes possessing a timelike Killing vector and when we ignore any radiation captured by the BHs. We also compute the bolometric radiative luminosity .
The second type of diagnostics includes: 10) Optical depth to true absorption (Eddington approximation), where we assume pure hydrogen and where () is the optical depth to electron scattering (freefree absorption), calculated using the Thompson scattering opacity (Rosseland mean opacity ) as () Shapiro:1983book (). implies the matter is optically thick to absorption. 11) Local temperature of the matter, calculated by solving , where is the proton mass and the Boltzmann constant. 12) In the cases with cooling, the effective disk temperature (in cooling cases), estimated by assuming that all cooling luminosity is emitted as black body radiation:
(9) 
where is the StefanBoltzmann constant and is the disk outer radius. 13) Characteristic observed thermal radiation frequencies , where h is the Planck constant and the cosmological redshift. This is calculated only when . 14) Cyclotron emission. While we find the bulk of the disk to be optically thick near Eddington accretion rates, the lowdensity “cavity” is optically thin. From these regions cyclotron lines may be detectable. 15) In cases where we compute the characteristic cyclotron frequencies , where is the electron charge, the electron mass and the magnitude of the magnetic field.
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.
iv.1 Hydrodynamical evolutions with
To set a lower limit on the necessary resolution to perform our GRMHD simulations, we found the minimum resolution required so that our code maintains stable equilibrium of an unmagnetized disk around a single nonspinning BH for several thousands of of evolution. The equilibrium disk solution we use is described at the beginning of Sec. III.1.2. Our study indicates that for the low resolution () listed in Table 2 the surface density profile of the initial disk is maintained to within throughout the disk for at least .
iv.2 MRI resolution requirements
Here we check the conditions for MRI to operate in our disk models. For this to be the case three conditions have to be satisfied: (I) A magnetic field configuration must be present that violates the stability condition for MRI , (II) The wavelength of the fastestgrowing mode has to be resolved by gridpoints Shibata:2006hr (); Sano:2003bf (); Shiokawa:2012 (), and (III) the Bfield must be sufficiently weak for . In other words the wavelength of the fastest growing mode should fit in the disk ^{1}^{1}1Condition (III) is satisfied everywhere shortly after our evolutions begin because magnetic winding converts poloidal magnetic field into toroidal, thereby decreasing . Our highresolution grids satisfy condition (II) in the innermost parts of the disk even after the evolution begins..
Regarding (I) our initial disk model is unstable to the MRI because of the outwardly decreasing angular velocity and the presence of an initially small poloidal magnetic field. Regarding (II) we plot the quality factor where is the local grid spacing (which jumps by a factor of two at AMR boundaries); see Fig. 2. One can see that we satisfy the criterion over a rather large portion of the disk initially except for the region near the radius where the poloidal field changes sign and becomes very small. We chose our lowresolution grids such that condition (II) is satisfied at in the innermost parts of the disk.
Regarding (III) we plot a meridional  slice of the restmass density overlayed by a line plot showing as a function of in Fig. 3. The plot shows that for the most part fits within the disk. At the inner disk edge the MRI is likely to be suppressed initially, but as the evolution proceeds magnetic winding converts poloidal fields into toroidal ones lowering and eventually triggering the MRI near the inner diskedge.
iv.3 Cooling
We seek to compare two opposite limiting cases for each mass ratio: (I) Nocooling, for which . Here is the local Keplerian time scale which is comparable to the (shock) heating time scale; (II) extremely rapid cooling, which we model with the effective emissivity .
To determine the value for at which cooling becomes rapid (at least in the bulk of the disk), we performed the cooling case using different cooling time scales. We concluded that rapid cooling requires a cooling time scale significantly shorter than the local Keplerian time scale. In Fig. 4, we plot , where the entropy parameter and for a run with (left panel), a run with (middle panel) and a run without cooling (right panel). It becomes apparent that when , is not driven back to its initial value. Physically, this means that not all shock generated entropy is radiated away, hence does not correspond to rapid cooling and steady state is not achieved. For we find in the bulk of the disk. The values depart from unity only inside the cavity where low density gas exists and can be shock heated to very high , demonstrating that shock heating is extremely strong in the lowdensity cavity. We adopt in all our cooling simulations because it leads to rapid cooling, at least throughout the bulk of the disk.
V Results
In this section we present the results of our numerical simulations. First, in Sec. V.1 we show results from our resolution study. In Sec. V.2, we directly compare the binary case with to the case. Lastly, we report on the variation of our diagnostics with mass ratio for all magnetized cases in Sec. V.3.
Our simulations have two parameters that scale out of the problem: the total mass of the BHBH binary and the restmass density of the disk. Alternatively, we can exchange one of these parameters with another parameter that depends on these two free parameters. So, instead of the restmass density, in the results we quote we choose the Eddington ratio , where is the bolometric EM luminosity described in Sec. III.3. The relation between these parameters is determined as follows: the accretion rate through the horizon must scale like , where is a reference restmass density in the disk. We choose the maximum restmass density at as the reference density, and our simulations determine the proportionality constant. For example, in the single, nonspinning BH case with cooling we find
where we have replaced the accretion rate via the following expression
(10) 
and where is the radiative efficiency as computed via our simulations for our adopted cooling law in the single, nonspinning BH case. In the nocooling cases the only radiation luminosity estimate we have is the Poynting luminosity which is expected to be a small fraction of the total radiative luminosity. Hence, in the nocooling cases we do not scale our results with the Eddington ratio. Instead, we choose a fiducial accretion rate similar to the one given in Eq. (10).
v.1 Resolution study
Here we present the results of our resolution study. For the single BH, nocooling and BHBH equal mass, nocooling cases we use the three resolutions (see Table 2).
In the single BHcase the average accretion rate varies little with resolution (see left panel in Fig. 5). The maximum fractional difference of the mean accretion rate between different resolutions is . Other quantities show a similar behavior. These results indicate that the resolutions used in this case are high enough for the main MRI effects to be captured and the results to be qualitatively independent of resolution. However, the resolution is not sufficiently high to perform a formal convergence test.
For the equalmass case we observe a different behavior. The mean accretion rates appear to converge (see right panel in Fig. 5), but the low resolution run underestimates the accretion rate by almost a factor of 2, while the medium and high resolution runs are in good agreement. For the latter resolutions mean accretion rates agree to within . We conclude that for the case our adopted medium and high resolutions are sufficient for drawing qualitative conclusions and that higher resolutions are necessary for accurate quantitative results that reside in the convergent regime.
The results of the resolution study in the and cases differ because of the distribution of matter in both cases and our grid setup. In the case there is more matter close to the BH, where very high resolution refinement levels reside, whereas in the case the bulk of the matter remains outside the inner edge of the disk, where the grid resolution is not as high. As we show later, this is not the case in our models. There, more matter resides closer to the BHs, and hence closer to the high resolution levels.
Based on our resolution study we conclude that the low resolution used in the equalmass case is not sufficiently high to yield reliable results, but for the other mass ratios we consider in this work, our adopted low resolution suffices for a qualitative discussion. Thus, in the equalmass cases results will be reported mainly from our medium resolution runs.
We stress here that these simulations, which include all relativistic effects and resolve the BHs, are computationally very demanding (see Sec. III) and increasing the resolution even further will incur a very large computational cost. With increasing computer power and larger computer allocations we plan to improve our results in the near future.
v.2 Significance of Bfields
Previous hydrodynamic simulations and (semi)analytic models of circumbinary accretion disks using the simplified disk model (e.g. Artymowicz:1994bw (); MacFadyen:2006jx ()) showed that the main feature in the equalmass binary case is that the density inside and near the binary’s orbit remains substantially lower than in the single BH case (see, e.g., Fig. 6). Such an inner cavity or “hollow” can have important consequences for the emergent radiation, such as line emission due to small optical depth and small bolometric luminosity from the hollow. Any such difference between single BH vs. circumbinary accretion disks can provide a path to distinguishing a binary AGN versus a classical, standard (single BH) AGN Tanaka:2013oju (). The explanation for the existence of a hollow is that the binary tidal torques for , are strong enough to push most matter away from the binary orbit Haiman:2009te (). The effect is most prominent in geometrically thin disks, which arise when radiative cooling is highly efficient.
However, even in the absence of viscosity or magnetic fields, the timedependent tidal field strips off matter from the inner disk edge, giving rise to an accretion pattern consisting of two streams which penetrate the inner cavity and extend to the horizons of the BHs Haiman:2009te (); Farris:2009mt (); Bode:2009mt (); Bode:2011tq (); Bogdanovic:2010he (); Farris:2011vx ().
Furthermore, recent MHD simulations Shi:2011us (); Noble:2012xz (); Farris:2012ux () universally revealed that the reduction of density inside the cavity in the binary case is not as substantial as previously thought. Such simulations explore regimes in which the disk is geometrically thick, which partially accounts for the difference.
Here we present a comparison between magnetized and unmagnetized circumbinary accretion disks onto an equal mass BHBH, while all other physical and numerical parameters remain identical, to illustrate the importance of magnetic fields in filling the hollow with dense material.
v.2.1 Midplanedensity
Figure 6 demonstrates the striking differences between nocooling evolutions with and without magnetic fields at the same integration time. Magneticfree hydrodynamic evolutions severely underestimate both the density in the inner regions and the overdensity due to the spiral arms. This indicates that the amount of matter stripped off by tidal torques is small compared to the amount of matter flowing into the hollow due to MHD turbulence.
v.2.2 profiles,
Next we compare the surface density profiles of magnetized vs unmagnetized disks. We find different profiles between the two cases as shown in Fig. 7. The model remains relatively close to the initial data, apart from a slow, mild expansion due to tidal heating and shocks. When magnetic fields are present, the final disk profile is completely different from the initial data even though the binary torques are identical. This implies that the MRIdriven viscous torques have a much larger impact on the global disk structure than the binary tidal torques, except perhaps near the inner disk edge.
v.2.3 Sensitivity to cooling
We find a fundamental difference between and evolutions regarding their sensitivity to cooling. In the case both and evolutions lead to essentially the same profile (see Fig. 7). This is in stark contrast to the (magnetized) cases shown in Fig. 7, for which cooling has a strong impact, leading to a matter pileup near the inner disk edge. As our particular choice of serves to dissipate heat from shocks, we conclude that magnetic fields lead to far stronger shock heating in the disk than the binary tidal torques.
v.2.4 Accretion rate
The ratio of the timeaveraged accretion rate without magnetic fields to that with magnetic fields is (see also Sec. V.3). This result applies to both cooling and nocooling cases.
In summary, evolutions underestimate the material inside the cavity and accretion rates by orders of magnitude. Hence, incorporating magnetic fields is paramount for a proper treatment of circumbinary accretion disks.
v.3 Trend with mass ratio,
In this section we discuss the dependence of our multiple diagnostics on the binary mass ratio for our cases. We use results from the single nonspinning BH case () to normalize and compare our results for the binary cases.
v.3.1 Disk structure in the bulk and inside the cavity
We begin this section by discussing the qualitative evolution of the disk restmass density . For , the early evolution is similar for the different mass ratios, but departs strongly in the subsequent evolution depending on . The onset of accretion occurs through two spiral streams, which remain attached to the horizons throughout the evolution, as shown in Figs. 8–15. Here we plot the restmass density contours in the equatorial plane. These streams are among the densest structures of the accretion flow, especially for the lowest nonzero mass ratio case (see Figs. 14 and 15). Spiral density waves are launched near the inner edge of the disk, which propagate and dissipate into the outer disk. This feature can also be seen in Figs. 8–15 for all mass ratios.
Late in the evolution we find that when , the density of the matter inside the “cavity” in cases is larger than that in cases (see Figs. 8 and 9). Hence, the amount of matter in the hollow is smaller when we allow rapid cooling. This dependence on cooling arises because of the larger disk thickness in the nocooling cases, which leads to a reduction in tidal torques Lin:1993 (); Armitage:2002uu () near the binary orbit, allowing matter to overflow more easily.
We also find that the smaller the mass ratio the more matter pours into the cavity. This is anticipated from the Newtonian expression for the binary tidal torque Lin:1993 (); Armitage:2002uu (), which decreases with decreasing . As expected, the largest contrast arises between the and cases, which becomes obvious by simply inspecting the restmass density contours in the equatorial plane as shown in Fig. 16. One can see the main difference: The presence of a region of lowered density with two accretion streams near the BHs in the binary case – the “hollow” (or “cavity”) – and the absence of these features in the case (left panel).
The relaxed disk structure in the predecoupling regime is not axisymmetric. The backsloshing of material towards the inner disk edge occurs mainly in two opposing directions and leads to a gradual overdense feature in the disk, which has been referred to as a “lump” Noble:2012xz (); Shi:2011us (), and its presence has been linked with a growth in disk eccentricity (see also MacFadyen:2006jx ()). The nonaxisymmetric feature is stronger for models with cooling. Hence, a more realistic calculation with radiative transfer is necessary to assess the strength of nonaxisymmetric structure in a circumbinary disk.
As expected, the restmass density contrast between the two accretion streams becomes larger with smaller mass ratios. This effect is easily seen when comparing the rest mass density contours in the equatorial plane for , and cases (see Figs. 12 and 15).
In the and nocooling cases we observe time variations in the density of the streams relative to each other: for about half an orbit one stream is stronger than the other.
We find that the supply of material channeled onto the BHs is sufficient to keep the BHs immersed in a persistent gaseous environment with . This means that the forcefree electrodynamics approximation may be inadequate to globally describe the systems considered here.
For there is hardly a lowdensity hollow (see Fig. 14). This is also revealed by the inner disk edge being close to or inside the orbit of the secondary, especially in the nocooling case. In the limit no hollow appears, as expected. However, we observe a region of lowered surface density near and inside the ISCO of the primary BH.
v.3.2 Inner disk edge
In Newtonian 2D studies of geometrically thin disks and large binary separations, the location of the inner disk edge is found to be roughly twice the binary separation, independent of ; see, e.g., Table I in Artymowicz:1994bw (). For the geometrically thick disks and binary separations we are considering, we find the inner disk edge in the relaxed state to be dependent on and whether cooling is enabled.
In all cases (see the snapshots in Figs. 8–15), the inner disk edge is closer than predicted by Newtonian thindisk calculations Artymowicz:1994bw (); MacFadyen:2006jx (); Shi:2011us (). In the equalmass cooling case the inner edge remains closer to the initial one (see Fig. 7). The trend is such that the smaller the mass ratio the closer the disk edge is to the binary orbit. In the nocooling case the inner disk edge effectively coincides with the orbit of the secondary. We report the value for found in each case in Table 3, and plot vs in Fig. 17.
Case  

1:1ncmr  –  
1:2nclr  –  
1:4nclr  –  
1:8nclr  –  
1:10nclr  –  –  
0nchr  ^{4}^{4}4For ease of comparison, in the single BH case we normalize to even though it does not correspond to an orbital separation.  –  –  
1:1cmr  
1:2clr  
1:4clr  
1:8clr  
1:10clr  –  
0clr  – 
v.3.3 Surface density
In Fig. 18 we show the surface density () profiles of the relaxed disks, averaged over the last 10 binary orbital periods (for all mass ratios, cooling and nocooling). For all cases the relaxed state deviates strongly from the initial profile, highlighting the importance of evolving the system for at least a viscous time scale at the radii of interest, where
(11) 
and where is the shear viscosity, with .
Nocooling: The evolutions for and are similar in terms of their profiles. The other cases () yield similar which extend further in than for and . The surface density diagnostic clearly demonstrates that the Newtonian thindisk feature that the cavity edge is at twice the binary separation and independent of the mass ratio does not hold in this class of runs.
Cooling: For nonzero mass ratios, cooling yields a pile up of dense gas near the inner disk edge, which is absent in the nocooling runs and is strongest for the equalmass case. Apart from the pileup at small radii, all evolutions have a rather similar profile at larger radii. In the cooling cases the cavity edge is farther out than for nocooling, but still closer than twice the binary separation. The dependence on mass ratio is weaker than in the nocooling cases, but still in some disagreement with the Newtonian thindisk calculations.
v.3.4 Effective stress
Figure 19 shows the averaged effective ShakuraSunyaev parameter profiles for all mass ratios, both for cooling and nocooling models in the relaxed state. The average is taken over the last 10 binary orbits. In Table 3 we also quote characteristic values obtained by additionally averaging over radii from the inner disk edge out to the location of the density maximum. We plot vs in Fig. 17.
In all cases we observe larger values for in the cooling cases than for the nocooling cases, and there is always a steep increase in stress at smaller radii. For the cooling case we find and for all other cooling cases. A typical value for all nocooling cases is . The higher stress for cooling cases results from the additional gas compression and associated amplification of magnetic fields when cooling is allowed.
v.3.5 Accretion rates
We compute the accretion rates through the individual BHs as well as the total accretion rate onto the binary, and normalize these by the (timeaveraged) single, nonspinning BH accretion rate.
For a given maximum restmass density and total BH mass, we find the highest accretion rate in the single BH case. This is consistent with the expectation that the absence of a tidaltorque barrier will allow matter to flow more easily toward the BH(s).
By contrast, the tidal torque is maximized for (all else being equal), so the expectation is that the accretion rate will be minimum for . In agreement with this expectation we find the accretion rate in the cases to be smaller than all other mass ratios we consider here.
In Table 3 we list the average accretion rate vs mass ratio, and plot the results in Fig. 17. The general trend is that lower mass ratios have higher accretion rates. In the case the average accretion rate is about that of the single BH case with the same initial maximum restmass density and total BHBH mass, while the average accretion rate in the equalmass case is roughly of that in the case.
For and both black holes accrete at comparable rates whether cooling is applied or not (see four upper rows, left panels Fig. 20). However, we observe that often the accretion rates on the individual BHs are antiphased, i.e., accretion occurs for half an orbit primarily on one BH and then for the second half of the orbit on the other BH. This behavior in the relaxed state is due to an “alternating” pattern in which denser material primarily plunges first through one stream and then through the other.
In Fig. 20 we also plot the accretion rates for and , with and without cooling (three lower rows, left panels). It is apparent that for the accretion rate onto the primary is comparable to that onto the secondary when . However, when the dominant contribution to the total accretion rate comes from the primary. In the case we observe that the dominant contribution to the total accretion rate comes from the primary whether cooling is applied or not, and the same holds true for the case.
This result in the case can be qualitatively understood through a rough analogy to BondiHoyleLyttleton accretion: The secondary has a smaller mass and moves faster on its orbit reducing its effective capture cross section as suggested by BondiHoyleLyttleton accretion (see e.g. Farris:2009mt ()). Also the surface area of the secondary is roughly a factor smaller than that of the primary. Note however, that the secondary plays a role in stripping matter off the inner disk edge effectively as it orbits closest to (or even through) the disk, so the accretion rate onto the secondary is not generally expected to be 100 times smaller than the accretion rate onto the primary.
In particular, in the nocooling case a dense persistent structure coorbits with the secondary (see Fig. 14). The density in this structure exceeds the density of matter near the primary by more than a factor of two.
v.3.6 Variability
We now report results from the Fourier analysis of the accretion rate. These can be seen in the right panels of Fig. 20. A summary of the primary Fourier modes for the different cases is presented in Table 3.
In the case a Fourier analysis of (accretion rate onto the primary) and (accretion rate onto the secondary) reveals a characteristic frequency near , in agreement with Farris:2012ux (). The analysis of (the total binary accretion rate) gives a dominant Fourier mode with a frequency twice as high. We observe a peak at the binary period only for and , and only for the cooling cases (see Fig. 20). These results indicate that if the variability in the accretion rate is directly translated into a variability of EM signatures, inferring the binary frequency from EM observations may not be straightforward. We observe variability for other mass ratios as well. The frequencies for the equal mass case also appear in other cases, in addition to other weaker contributions, but a clear trend is not evident. The most prominent and clean periodic signature occurs in the nocooling case (see Fig. 20 and discussion in D'Orazio:2012nz ()). Other strong Fourier modes are observed in cooling and the cases. For and no significant periodicities are observed.
In the case the variability is dominated by variations in the accretion flow onto the primary. The Fourier analysis yields rather irregular accretion, i.e. not very pronounced frequencies. The secondary accretes at several pronounced frequencies, but the amplitude of the variations is much smaller.
In D'Orazio:2012nz (); Farris:2013uqa () the dependence on of the variability was studied for geometrically thin (“locally isothermal” disks) and proposed as a key feature to observationally distinguish accreting BHBHs from standard, single BH AGNs. In our Fourier analysis the individual peaks are less significant than in D'Orazio:2012nz (); Farris:2013uqa () and the Fourier spectrum yields a more complex structure. This discrepancy is likely due to a combination of additional effects including differences in the viscosity prescription (i.e. MHD turbulence vs. viscosity), cooling prescriptions, thin vs. thick disks, 2D vs 3D, and the EOS. These differences result in geometrically thin vs thick disks, which are seen to have different variability.
v.3.7 Luminosities
We compute the cooling () and Poynting () luminosities, as well as the energy loss rate due to outflowing matter () for all cases. We typically find is comparable to and quite smaller than . We highlight representative cases and summarize all values for the different characterizing the relaxed state in Table 3.
The large variability seen in the accretion rate is only partially reflected in the cooling luminosities and not reflected in the Poynting luminosities. However, these conclusions need to be confirmed with a selfconsistent treatment involving radiative transfer.
In all cooling cases we find that the cooling luminosity is significantly larger than the Poynting luminosity and the energy loss rate due to outflowing matter by almost a factor of 10 (see Table 3).
We estimate and compare the contributions to the cooling luminosity from various regions in the disk: the outer disk, the inner edge, and the cavity (see Fig. 21). We find that although the outer disk gives the largest contribution, the inner edge and cavity interior are a substantial portion () of the total cooling luminosity. Therefore, the activity in the cavity cannot be ignored, as has been done in earlier studies.
v.3.8 Opacities
We estimate the Thompson scattering () and freefree absorption () optical depths in all cases. We do not find a strong dependence of , and on the mass ratio. Our crude analysis shows