Accretion disks around binary black holes of unequal mass: GRMHD simulations near decoupling

Accretion disks around binary black holes of unequal mass:
GRMHD simulations near decoupling

Roman Gold, Vasileios Paschalidis, Zachariah B. Etienne, Stuart L. Shapiro Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801
Physics Department & Joint Space-Science Institute, University of Maryland, College Park, MD 20742
Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771
Department of Mathematics, West Virginia University, Morgantown, WV 26506
Department of Astronomy & NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801
Harald P. Pfeiffer
Canadian Institute for Theoretical Astrophysics, University of Toronto, ON M5S 3H8, Canada
Canadian Institute for Advanced Research, Toronto, ON M5G 1Z8, Canada

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 binary-disk decoupling radius. We compare (surface) density profiles, accretion rates (relative to a single, non-spinning 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 no-cooling 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.

04.25.D-, 04.25.dg, 47.75.+f

I Introduction and Motivation

Supermassive black hole (SMBH) binaries can form in magnetized plasma following galaxy mergers, via bar-mode instability in rapidly rotating supermassive stars, or by other dynamical processes Begelman:1980vb (). After formation, a combination of dynamical friction and gas-driven migration is likely to catalyze the binary inspiral into the gravitational radiation-driven regime Ivanov:1998qk (); Haiman:2009te (); Rafikov:2012hd (); Begelman:1980vb (); Armitage:2002uu (); Armitage:2005xq (); Milosavljevic:2001vi (). The exact details of these processes, including the “last-parsec 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 multi-messenger 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 gravitational-wave 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, very-long 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 (X-shaped 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 Bondi-Hoyle-Lyttleton 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 binary-disk 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 so-called disk-dominated 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 secondary-dominated 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 quasi-steady state - this is called the predecoupling regime. When , the binary decouples from the disk, i.e. the inward migration of the binary out-paces 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 lrr-2013-1 () for a recent review). Many different disk models have been proposed in the literature. These models range from geometrically-thin, 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 Post-Newtonian 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 self-consistently has been discussed in Roedig:2012nc (), and only full GR calculations can achieve this goal reliably.

A “GR-hybrid” orbit-averaged 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 B-field 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 force-free 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 quasi-instantaneous 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 wall-clock 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 thin-disk 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 conformal-thin-sandwich 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 low-density “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 pressure-dominated, optically thick disks. This choice is motivated in part both by theory and observations. First, accretion disk theory lrr-2013-1 () 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 (), steady-state disk models predict that radiation pressure-dominated circumbinary disks will channel more material into the cavity for a given central mass. This finding makes radiation pressure-dominated 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 (), XMM-COSMOS 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 pressure-dominated 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 Shakura-Sunyaev/Novikov-Thorne thin-disk model Shakura73 (); Novikov73 () to make rough analytic estimates regarding the physical regime our disk models probe. While the Shakura-Sunyaev/Novikov-Thorne 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 Shakura-Sunyaev/Novikov-Thorne model

The simulations reported here apply to any total (ADM) binary black hole mass and, since we neglect the self-gravity of the gas, to any rest-mass 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 magnetic-to-total disk pressure, the initial magnetic field profile, the initial binary separation and the cooling law. In this section we use the steady-state thin-disk solution for a Shakura-Sunyaev/Novikov-Thorne 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 pressure-dominated, and this fact motivates our setting .

Neglecting the perturbative role of the secondary tidal torque, the steady-state accretion flow in a geometrically thin disk is uniquely specified by the Shakura-Sunyaev viscosity parameter , the central mass and the accretion rate . The quantity is specified in turn by the disk luminosity and efficiency . The Shakura-Sunyaev/Novikov-Thorne model Shakura73 (); Novikov73 () describes a disk that is radiation pressure-dominated inside a radius . In this region the radiation to gas pressure ratio is lrr-2013-1 ()


where is the radiative efficiency, , and we chose the normalization for the parameter based on typical values found in our simulations. For a thin-disk model the efficiency ranges from for a non-spinning 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


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 pressure-dominated region.

In this region the typical rest mass densities are lrr-2013-1 ()


As we will show later, these characteristic values for and are comparable to the values found in our simulations.

The Shakura-Sunyaev/Novikov-Thorne model predicts that the effective optical depth to absorption is adopting the same normalizations as in the equation above. This well-known inconsistency near the Eddington limit of the Novikov-Thorne 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 slim-disk models, which remain optically thick in this high luminosity limit Abramowicz1988 (); lrr-2013-1 ().

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 ()


where we assumed that the inner disk edge radius settles to as typically found in our simulations, is the Shakura-Sunyaev 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 self-gravitating 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 quasi-equilibrium conformal-thin-sandwich (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
Table 1: CTS initial data parameters for the BHBH vacuum spacetime. Columns show mass ratio (), ADM mass () and angular momentum (), and irreducible masses (), and apparent horizon radii () for the two black holes. Diagnostics generating these quantities, but computed from vacuum, test simulations agree with these values to within one part in .

iii.1.2 Matter and B-field 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 non-spinning 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 B-field 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.

Figure 1: B-field lines (white curves) and volume rendering of rest-mass density normalized to its maximum value at . The two black dots at the center indicate the BH apparent horizons for the case.

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 adaptive-mesh-refinement (AMR) code Duez:2005sf (); Etienne:2010ui (); Etienne:2011re (), which adopts the Cactus/Carpet infrastructure Goodale2002a (); Carpet-cactusweb (); Carpet-carpetweb (). 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 flux-conservative formulation [see Eqs. (27)-(29) in Etienne:2010ui ()] employing a high-resolution-shock-capturing 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 zero-divergence 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 long-term GRMHD evolutions. The “algebraic” gauge condition used in the first GRMHD simulations adopting A-field 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 vice-versa, 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


where is the total pressure, the rest-mass density, and the specific internal energy. This EOS allows for shock heating. We choose appropriate for radiation pressure-dominated, 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 center-of-mass 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 slim-disk and ADAF models lrr-2013-1 (), is self-consistently 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 steady-state 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:

  1. Non-exponential cooling Shafee:2008mm (); Noble:2012xz ():


    where is an entropy parameter, and the initial target value. We call this emissivity “non-exponential”, because the effective cooling time scale for this scheme is not just , but depends on the internal energy (see Appendix A).

  2. Exponential cooling Paschalidis:2011ez (); Paschalidis:2012ff ()


    where is the internal energy calculated using , and . In this scheme the effective cooling time scale is .

Both emissivities dissipate all shock-induced 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 low-density regions (see e.g. Appendix B in Etienne:2008re ()). Thus, to stabilize the simulations with prescription I, one typically excludes cooling of the low-density regions or unbound matter Noble:2012xz (); Farris:2012ux (). As both the BH horizons and the low-density 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 no-cooling cases.

Case cooling?  levels(BH)  levels(bh) Grid hierarchy
1:1nc-hr 1:1 No 7 7 ,
1:1nc-mr 1:1 No 7 7 ,
1:1nc-lr 1:1 No 7 7 ,
1:2nc-lr 1:2 No 7 8 ,
1:4nc-lr 1:4 No 7 9 ,
1:8nc-lr 1:8 No 6 10 ,
1:10nc-lr 1:10 No 6 10 ,
0nc-hr 0 No 6 ,
0nc-mr 0 No 6 ,
0nc-lr 0 No 6 ,
1:1c-mr 1:1 Yes 7 7 ,
1:2c-lr 1:2 Yes 7 8 ,
1:4c-lr 1:4 Yes 7 9 ,
1:8c-lr 1:8 Yes 6 10 ,
1:10c-lr 1:10 Yes 6 10 ,
0c-lr 0 No 6 ,
Table 2: 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 M in all cases. From left to right the columns indicate the case label, mass ratio , whether cooling is included or not, the coarsest grid spacing , 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 , where is the level number such that corresponds to the coarsest level. A dash “–” indicates “no information available”.

iii.2.4 Evolution grids & models

We use a hierarchical, box-in-box adaptive mesh provided by the Cactus/Carpet infrastructure Carpet-cactusweb (); Carpet-carpetweb (). 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, non-spinning 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:1nc-hr means mass ratio , no-cooling, 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 36-node 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 orbit-averaged disk models. 4) Shakura-Sunyaev -stress parameter computed as where is the dominant orthonormal component to the Maxwell stress-energy 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 orbit-averaged 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 orbit-averaged 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 steady-state 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 (free-free 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:


where is the Stefan-Boltzmann 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 low-density “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 non-spinning 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 fastest-growing mode has to be resolved by gridpoints Shibata:2006hr (); Sano:2003bf (); Shiokawa:2012 (), and (III) the B-field must be sufficiently weak for . In other words the wavelength of the fastest growing mode should fit in the disk 111Condition (III) is satisfied everywhere shortly after our evolutions begin because magnetic winding converts poloidal magnetic field into toroidal, thereby decreasing . Our high-resolution 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 low-resolution grids such that condition (II) is satisfied at in the innermost parts of the disk.

Figure 2: Contours of the -quality factor in the equatorial plane, corresponding to the equal mass () medium resolution case [divide (multiply) by () for the low (high) resolution case] at . We resolve the fastest growing MRI mode by gridpoints over a large part of the disk (the blue ring stems from extremely small values of , when the vertical component of the B-field changes sign).

Regarding (III) we plot a meridional - slice of the rest-mass 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 disk-edge.

Figure 3: Rest-mass density contours (color coded) on a meridional slice, and (black solid line) at . The plot corresponds to equal mass () but is the same among all cases considered in this work. For the most part fits within disk.

iv.3 Cooling

We seek to compare two opposite limiting cases for each mass ratio: (I) No-cooling, 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 low-density cavity. We adopt in all our cooling simulations because it leads to rapid cooling, at least throughout the bulk of the disk.

Figure 4: for three runs with different at . Left panel: . Middle panel: . Right panel: No-cooling (). Cooling in the case cannot keep up with heating, causing to drift away from its initial value. Values much smaller than are required for rapid cooling. Steady state is not achieved in the cases.

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 rest-mass 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 rest-mass 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 rest-mass density in the disk. We choose the maximum rest-mass 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


and where is the radiative efficiency as computed via our simulations for our adopted cooling law in the single, nonspinning BH case. In the no-cooling 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 no-cooling 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).

Figure 5: Left panel: rest-mass accretion rate vs time for the no-cooling cases at low (lr), medium (mr), and high (hr) resolutions. Right panel: rest-mass accretion rate vs time for the no-cooling cases at low, medium, and high resolutions. The accretion rate in the () case is normalized by the mean accretion of the high resolution () run. The horizontal lines indicate the mean accretion rate at low, medium and high resolution.

v.1 Resolution study

Here we present the results of our resolution study. For the single BH, no-cooling and BHBH equal mass, no-cooling cases we use the three resolutions (see Table 2).

In the single BH-case 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 equal-mass 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 equal-mass 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 equal-mass 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 B-fields

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 equal-mass 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 time-dependent 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 Midplane-density

Figure 6 demonstrates the striking differences between no-cooling evolutions with and without magnetic fields at the same integration time. Magnetic-free 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.

Figure 6: Contours of rest-mass density normalized to the initial maximum rest-mass density in the equatorial plane at for two no-cooling cases. Left panel: . Right panel: . Notice the higher densities in spiral arms in the inner regions in the case.

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 MRI-driven 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 pile-up 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.

Figure 7: Surface density profiles for several cases: Initial profile (yellow, solid line), no-cooling (green circles) & cooling (green crosses), no-cooling (red, solid line) & cooling (black, dashed line) cases. Apart from the initial profile all other profiles are orbit-averaged over the last 10 orbits of evolution, beyond which the profile does not change appreciably. Notice the clear emergence of a pile-up near the inner disk edge in the cooling case, as well as the relatively small change relative to the initial data in the cases, and the similarity between and unmagnetized cases.
Figure 8: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the no-cooling case. Here .
Figure 9: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the cooling case. Here .
Figure 10: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the no-cooling case. Here .
Figure 11: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the cooling case. Here .
Figure 12: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the no-cooling case. Here .
Figure 13: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the cooling case. Here .
Figure 14: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the no-cooling case. Here .
Figure 15: Contours at select times of rest-mass density normalized to the initial maximum (log scale) in the equatorial plane. The plot corresponds to the cooling case. Here . The gas is denser everywhere compared to mass ratios closer to unity; compare to Figs. 8, 9, 10, 11, 12, 13.

v.2.4 Accretion rate

The ratio of the time-averaged accretion rate without magnetic fields to that with magnetic fields is (see also Sec. V.3). This result applies to both cooling and no-cooling 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 non-spinning BH case () to normalize and compare our results for the binary cases.

Figure 16: Contours of rest-mass density (linear color scale) normalized to the initial maximum in the equatorial plane at . Left panel: no-cooling. Right panel: no-cooling case.

v.3.1 Disk structure in the bulk and inside the cavity

We begin this section by discussing the qualitative evolution of the disk rest-mass 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. 815. Here we plot the rest-mass density contours in the equatorial plane. These streams are among the densest structures of the accretion flow, especially for the lowest non-zero 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. 815 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 no-cooling 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 rest-mass 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 back-sloshing 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 rest-mass 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 no-cooling 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 force-free electrodynamics approximation may be inadequate to globally describe the systems considered here.

For there is hardly a low-density 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 no-cooling 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. 815), the inner disk edge is closer than predicted by Newtonian thin-disk calculations Artymowicz:1994bw (); MacFadyen:2006jx (); Shi:2011us (). In the equal-mass 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 no-cooling 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.

0nc-hr     444For ease of comparison, in the single BH case we normalize to even though it does not correspond to an orbital separation.
Table 3: Table summarizing our main results. Columns show case label, inner disk edge normalized to the binary separation , the mean accretion rate and its approximate standard deviation normalized to the no-cooling single BH mean accretion rate , main Fourier frequencies normalized by the binary orbital period in the Fourier analysis of the accretion rate, Shakura-Sunyaev -parameter333In some cases we quote a range of values because a single radially averaged value would overestimate ., the mean Poynting luminosity , and the mean cooling luminosity both normalized to the binary mean accretion rate. A dash “–”  indicates “no information available”.

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 no-cooling). 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


and where is the shear viscosity, with .

Figure 17: Mean accretion rate (normalized to the no-cooling single BH accretion rate), , and (normalized to twice the binary separation) as functions of for cooling and no-cooling cases (all at low resolution). In the absence of medium and high resolution runs for we place error bars in based on the resolution study. These error bars are chosen to be 50%, corresponding to the fractional difference between the high and low resolutions runs in the case (see Sec. V.1). The error bar in the case is 15% corresponding to the fractional difference between the high and low resolutions runs in the case.

No-cooling: 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 thin-disk 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 non-zero mass ratios, cooling yields a pile up of dense gas near the inner disk edge, which is absent in the no-cooling runs and is strongest for the equal-mass case. Apart from the pile-up 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 no-cooling, but still closer than twice the binary separation. The dependence on mass ratio is weaker than in the no-cooling cases, but still in some disagreement with the Newtonian thin-disk calculations.

Figure 18: Disk surface density -profile for different mass ratios as a function of cylindrical radius (normalized to the binary separation ). Left panel: no-cooling cases. Right panel: cooling cases. Except for the initial profile all other curves correspond to the relaxed state averaged over the last 10 binary orbital periods. In all these cases the relaxed profile deviates strongly from the initial one. In the sequence the density reaches further toward small radii with decreasing . In the sequence, one can see the pile-up near the inner disk edge, which is strongest for large .

v.3.4 Effective -stress

Figure 19 shows the averaged effective Shakura-Sunyaev parameter profiles for all mass ratios, both for cooling and no-cooling 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 no-cooling 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 no-cooling cases is . The higher stress for cooling cases results from the additional gas compression and associated amplification of magnetic fields when cooling is allowed.

Figure 19: Shakura-Sunyaev profiles for all mass ratios. Left panel: no-cooling cases. Right panel: cooling cases. Profiles have been averaged over the last 10 binary orbital periods. Cooling cases tend to yield larger values. In all cases the stress increases inside the cavity.
Figure 20: Left panels: Late evolution of mass accretion rate onto the binary (solid, black lines), primary (red, dashed lines), and secondary (solid-dotted, yellow lines). Right panels: Fourier transform of the accretion rates shown on the left. From top to bottom we show the cooling and no-cooling, cooling and no-cooling, cooling and no-cooling, and no-cooling cases.

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 (time-averaged) single, non-spinning BH accretion rate.

For a given maximum rest-mass 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 tidal-torque 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 rest-mass density and total BHBH mass, while the average accretion rate in the equal-mass 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 anti-phased, 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 Bondi-Hoyle-Lyttleton accretion: The secondary has a smaller mass and moves faster on its orbit reducing its effective capture cross section as suggested by Bondi-Hoyle-Lyttleton 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 no-cooling case a dense persistent structure co-orbits 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 no-cooling 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 self-consistent 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.

Figure 21: Contribution to from the inner cavity and outer disk for the cooling case. The notation means the cooling luminosity integral , where is the volume within a coordinate sphere of radius . Note that about of the total cooling luminosity arises within .

v.3.8 Opacities

We estimate the Thompson scattering () and free-free absorption () optical depths in all cases. We do not find a strong dependence of , and on the mass ratio. Our crude analysis shows