SEARCH FOR SUPERMASSIVE BLACK HOLE BINARIES IN THE SLOAN DIGITAL SKY SURVEY SPECTROSCOPIC SAMPLE

There is indirect evidence for the presence of SMBH pairs at ∼kpc scales in single galaxies detected via both direct imaging (Junkkarinen et al. 2001; Komossa et al. 2003; Ballo et al. 2004; Comerford et al. 2009b; Fabbiano et al. 2011) and double-peaked narrow emission lines (Comerford et al. 2009a; Wang et al. 2009; Smith et al. 2010; Liu et al. 2010a, 2010b; Shen et al. 2011a). A parsec-scale BH pair CSO 0402+379 was detected with the Very Long Baseline Array (Maness et al. 2004; Rodriguez et al. 2006), and OJ 287 is believed to be a sub-parsec scale SMBH binary due to its periodic luminosity variations with a period of ∼11.86 yr (Takalo 1994; Pursimo et al. 2000; Valtonen et al. 2008). Indirect evidence for sub-parsec-scale pairs includes the precession of radio jets (Begelman et al. 1980; Liu & Chen 2007), double–double radio galaxies (Schoenmakers et al. 2000; Liu et al. 2003), the X-shaped radio sources (Dennett-Thorpe et al. 2002; Merritt & Ekers 2002; Liu 2004), and "light deficits" in the cores of elliptical galaxies (Milosavljević et al. 2002; Ravindranath et al. 2002; Graham 2004; Merritt 2006; Kormendy & Bender 2009).

However, there are no confirmed detections of sub-parsec SMBH binaries to date (Eracleous et al. 2012), although they are theoretically long-lived. Gaskell (1996) interpreted 3C 390.3 with double-peaked broad emission lines as two SMBHs orbiting around each other, both associated with a broad-line region (BLR). However, long-term monitoring of this object disproved the binary BH interpretation (Eracleous et al. 1997; Shapovalova et al. 2001). In general, Doppler motions in an accretion disk are favored to explain the double-peaked broad emission line quasars (Eracleous & Halpern 2003; Strateva et al. 2003; Wu & Liu 2004; Gezari et al. 2007; Lewis et al. 2010).

As an alternative means to identify BH binaries via their orbital motion, there have been a number of searches for large velocity displacements between single-peaked broad Balmer (mostly Hβ) emission lines and the rest frame of the quasar host as traced (typically) by the narrow emission lines. These velocity shifts may signal orbital motion in a binary. Promising individual candidates include 3C 227 and Mrk 668 (Gaskell 1983, 1984), SDSS J092712.65+294344.0 (Dotti et al. 2009; Bogdanović et al. 2009; Vivek et al. 2009; Shields et al. 2009; Heckman et al. 2009; Decarli et al. 2009b), SDSS J153636.22+044127.0 (Boroson & Lauer 2009; Chornock et al. 2009, 2010; Gaskell 2009; Wrobel & Laor 2009, 2010; Decarli et al. 2009a), SDSS J093201.60+031858.7 (Barrows et al. 2011), 4C+22.25 (Decarli et al. 2010), and SDSS 0956+5128 (Steinhardt et al. 2012). For all of these initial extreme objects, alternate explanations to BH binaries have been proposed, including recoiling BHs (Komossa et al. 2008; Eracleous et al. 2012; Steinhardt et al. 2012) or a perturbed accretion disk around a single BH (Chornock et al. 2010; Gaskell 2010).

There have now been two systematic searches through the Sloan Digital Sky Survey (SDSS) for single-peaked broad emission lines with large velocity offsets from the quasar rest frame. Tsalmantza et al. (2011) found five new candidates with large velocity offsets in SDSS Data Release 7 (DR7; York et al. 2000) quasars at 0.1 < z < 1.5. Eracleous et al. (2012) found 88 candidates with a broad-line peak offset Δv > 1000 km s⁻¹ out of the 15,900 z < 0.7 quasars from SDSS DR7, corresponding to binaries with r  ∼  few  ×  0.1 pc. These surveys favor binary systems where only one of the BHs in the pair is actively accreting, and they are biased toward binary orbital phases where the line-of-sight (LOS) velocities are larger and thus LOS velocity drifts are smaller. It is theoretically plausible that only one BH in the pair should radiate since binary BHs will torque the circumbinary disk and open a cavity, analogous to type II planet migration (Lin & Papaloizou 1986). This cavity may at least partially screen ambient gas from the primary BH (Bogdanović et al. 2008; Cuadra et al. 2009; Hayasaki et al. 2007, 2008).

After identifying candidate BH binaries from the SDSS spectroscopy based on large velocity offsets, Eracleous et al. (and Decarli et al. 2013) have obtained a second epoch of spectroscopy to search for the expected radial velocity offsets from orbital motion. They identify 14 candidates with velocity shifts between the SDSS and second-epoch spectra due to orbital motion of the BHs, whose resulting accelerations lie between −120 and +120 km s⁻¹ yr⁻¹. Their methodology is very powerful, but comes with a few limitations. First, as we mentioned, the requirement of a large velocity offset favors binary systems at orbital phases with large LOS velocity offsets and thus small LOS accelerations. It is necessary to wait a long time for these systems to experience significant velocity changes to truly confirm their nature. Second, the method requires a robust determination of the systemic velocity of the system: at high redshift, where most of the quasars and most of the mergers are, narrow emission lines are weak to nonexistent, and so this technique is limited to z < 0.8. Third, and most important, the method selects extreme objects with highly nonstandard broad lines. While these may be the result of binary motion, they may also be the result of highly uncertain BLR physics. Fourth, much like double-peaked AGNs are seen to vary perhaps due to hot spots in the disk (e.g., Gezari et al. 2007), these sources show changes in their line shapes when the line luminosity changes, providing a significant source of noise in the determination of orbital motion (e.g., Decarli et al. 2013). For all of these reasons, we here pursue a complementary approach with no preselection on line shape.

We propose a complementary approach to that of Eracleous et al. (2012) and Decarli et al. (2013), to search for SMBH binaries using temporal changes in the velocity of broad emission lines (see also Shen et al. 2013). Specifically, we will search for radial velocity shifts in normal, single-peaked broad emission lines in all quasars with multi-epoch spectroscopy. In the past, such time-resolved spectroscopy was unavailable, leading to preselection of extreme velocity outliers. However, with modern multi-object spectrographs, time-domain surveys are becoming routine. In the SDSS main survey, >2000 quasars with Mg ii have been observed multiple times, and we focus on these in our pilot experiment here. We emphasize that while we begin with this modest sample, there are already >8000 quasars with multiple epochs of spectroscopy observed with the Baryon Oscillation Spectroscopic Survey (BOSS; Schlegel et al. 2009) to which we can apply the same method, and that upcoming surveys such as an extended BOSS survey and the Prime Focus Spectrograph (PFS) on Subaru (Takada et al. 2012) may provide many more.

Our approach is sensitive to scenarios where either each BH has its own BLR or only one of the BHs is accreting. In the first scenario, where each BH has its own BLR, the broad lines are double-peaked. In this case, the orbital motion of the binary would lead to shifts in velocity of the double peaks or at least distortion of the broad emission lines (e.g., from double-peaked to single-peaked lines when the LOS velocities of the binary BHs decrease from maximum to zero during one-fourth of their orbit; Shen & Loeb 2010). Our approach is sensitive to these types of line-shape changes in principle, while search methods looking for persistent double peaks or extremely high velocity offsets between broad and narrow lines may miss them.

If only one of the binary BHs is accreting, the single-peaked broad lines will trace the orbital motion of one of the binary BHs. Our method favors the phases of the binary orbit with large drifts in the projected radial velocity (i.e., acceleration) of the secondary BH, while previous studies are tuned to catch the maximum projected orbital velocity along the LOS. Since we do not require narrow lines to identify the velocity rest frame, we can in principle search AGNs at all redshifts by studying broad emission lines in the rest-frame optical and ultraviolet part of the AGN spectra (e.g., Mg ii, C iii], C iv, Lyα, etc.).

Of course, our approach has disadvantages and complications as well. Like Eracleous et al. (2012) and Decarli et al. (2013), we are sensitive to intrinsic variability in the broad emission lines. Also, because our targets are so luminous, we are heavily biased toward near-equal mass ratios, limiting the demographic reach of our survey. Finally, and most difficult to quantify, we rely on each BH having its own BLR. If the BLR surrounds both sources, we will no longer be sensitive to radial velocity shifts; given the major uncertainties in BLR physics, this is a major caveat.

Our paper will proceed as follows. We describe the selection and basic properties of our sample in Section 2, and we introduce the cross-correlation method in Section 3. The results of our experiment are presented in Section 4, and we estimate the observability of SMBH binary systems with our detection strategy assuming various scenarios for their evolution in Section 5 and predict the observability with data from next-generation surveys in Section 6. We summarize and conclude in Section 7.

There are various sources of uncertainties in our measurements, including wavelength calibration, sky subtraction, continuum subtraction, and cross-correlation errors. The statistical error reported by the cross-correlation procedure is typically <20 km s⁻¹, with only ∼3% of objects having uncertainties greater than 100 km s⁻¹. We further create mock spectra using the SDSS variance arrays and perform input–output tests of our cross-correlation analysis. We recover a very narrow distribution of uncertainties in our output velocity shifts (<34 km s⁻¹ in all cases). Random errors alone cannot explain the width of the velocity distribution that we observe.

One source of systematic uncertainty comes from our choice of wavelength region for the cross-correlation. While the cross-correlation over the whole range [1700 Å, 4000 Å] helps us remove some false-positives, it also introduces further potential sources of uncertainty. For instance, absorption from the Ca H+K lines around the 4000 Å break, or errors in sky subtraction, may confuse our technique. As a sanity check, we did another set of cross-correlations for the seven candidates over [1700 Å, 3800 Å] to exclude Ca H+K absorption lines. Only three of the seven candidates (SDSS J002444.11+003221.4, SDSS J095656.42+535023.2, SDSS J004918.98+002609.4) survive this test. These are our three most reliable candidates. As a result, we consider seven as an upper limit in what follows.

There are many other contributions to the velocity uncertainties. The first effect is uncertainties in sky subtraction, as already discussed. Second, intrinsic changes in the quasar continuum may be a source of uncertainty. When we change the spectral region used for the continuum fit, the output velocities can change by tens of km s⁻¹, showing that our continuum fitting also contributes to the error budget. Third, there is irreducible scatter from temporal variations in the broad lines themselves (e.g., Sergeev et al. 2007; Gezari et al. 2007; Decarli et al. 2013). We do not yet know the full range of velocity shifts resulting from line-shape variability. Therefore, we cannot take the observed distribution of velocity offsets and infer anything about the binary population below our detection threshold, as in Badenes & Maoz (2012).

Let us start with the observability of a BH binary P_obs for a given total mass M_BH, mass ratio q, separation r, observing time interval Δt, and orientation on the sky. Note that the total mass M_BH is equal to (1 + q)/q times the BH mass shown in Figure 1(d), since we assume that we are observing the secondary BH. Although simulations show that very high eccentricity may arise through binary–star interactions (e.g., Sesana et al. 2008; Sesana 2010; Preto et al. 2011; Khan et al. 2012) and binary–disk interactions (e.g., Armitage & Natarajan 2005; Cuadra et al. 2009; Roedig et al. 2011), the eccentricity of the binary orbit is highly uncertain. For simplicity, we assume zero eccentricity throughout, so that the intrinsic velocity of the secondary BH is

$$\begin{equation} V_0 = \frac{1}{1+q} \sqrt{\frac{GM_{\rm BH}}{r}}. \end{equation} \tag{ 4 }$$

The observed velocity is V_obs = V₀sin ϕcos θ, where ϕ and θ are the azimuthal and polar angle respectively, of the observer on the sphere centered on the BH binary. θ = 0 when the binary orbit is edge-on relative to the observer, and ϕ = 0 when the secondary (accreting) BH lies between the observer and the center of mass of the binary. If this quasar is observed twice with a time interval of Δt, then the drift in velocity of the secondary BH is

$$\begin{eqnarray} |\Delta V_{{\rm obs}}| &=& V_0 \sin (\phi + \Omega \Delta t) \cos \theta - V_0 \sin \phi \cos \theta \nonumber \\ &\approx & V_0 |\cos \theta | |\cos \phi | |\sin (\Omega \Delta t)|, \end{eqnarray} \tag{ 5 }$$

where $\Omega = \sqrt{GM_{\rm BH}/r^3}$ is the angular velocity.

If the observer has equal chances to be at any solid angle on this sphere, then ϕ is uniformly distributed in [0, 2π] and |cos θ| also has an equal chance to have any value between [0, 1] according to dΩ = cos θdθdϕ = d(sin θ)dϕ. Let us say that the smallest velocity drift we can measure is ΔV_lim (ΔV_lim = 340 km s⁻¹ for our whole sample and ΔV_lim = 280 km ⁻¹ for the high subsample). We now calculate the probability that |ΔV_obs| exceeds the observational limit.

For a given cos θ, the probability of |ΔV_obs| > ΔV_lim, i.e.,

$$\begin{equation} |\cos \phi | > \frac{\Delta V_{\rm lim} }{V_0 |\cos \theta | |\sin (\Omega \Delta t)|} \end{equation} \tag{ 6 }$$

in the range ϕ ∈ [0, 2π], is

$$\begin{equation} \frac{2}{\pi } \arccos \left(\frac{\Delta V_{\rm lim} }{V_0 |\cos \theta | |\sin (\Omega \Delta t)|} \right). \end{equation} \tag{ 7 }$$

If we weight this value with the probability distribution of |sin θ| in the range [0, 1], the total probability becomes

$$\begin{eqnarray} &&P_{{\rm obs}} (r, M_{\rm BH}, q, \Delta t) \nonumber \\ &=& \int _0^{|\sin \theta _{{\rm max}}|} \frac{2}{\pi } \arccos \left(\frac{\Delta V_{\rm lim} }{V_0 |\cos \theta | |\sin (\Omega \Delta t)|} \right) d(|\sin \theta |)\nonumber \\ &=& 1 - \frac{\Delta V_{\rm lim} }{V_0 |\sin (\Omega \Delta t)|}, \end{eqnarray} \tag{ 8 }$$

where θ_max satisfies ΔV_lim/(V₀|cos θ_max||sin (ΩΔt)|) = 1.

We calculate the observability above for each object in our sample and sum them, which gives us the total number of candidates. The dependence of the total expected number of candidates on the mass ratio and separation of the BH binary are shown in Figures 5 and 6, respectively. In Figure 5 we fix the separation to be 0.1 pc and vary the mass ratio q from 0.01 to 1. The observability is moderately sensitive to the mass ratio. The expected number of observable sources increases by a factor of 10 as q decreases by a factor of 100. This sensitivity is mainly due to the dependence of the total BH mass M_BH on q, as we assume that we observe the secondary BH. The number of observable candidates increases more steeply as the radius decreases (Figure 6), because the orbital velocity of the BHs is proportional to r^−1/2. No BH binaries are observable if the separation gets larger than ∼0.3 pc for q = 1 (or ∼1 pc if q = 0.1). We do additional calculations at fixed separation while truncating the observability at R_BLR for each object (dash-dotted and dotted lines in Figure 6). In Section 5.3 we will compare our number of candidates with the models in the scenario where all quasars live at a fixed separation forever.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We now calculate the second ingredient needed to compute P_obs using Equation (2)—the residence time t_res as a function of the binary separation r for a given set of binary characteristics.

From Table 1, our candidate SMBH binaries are at separations of 0.01–0.1 pc or closer, well inside the regime of "the final parsec." As mentioned in the Introduction, there are many mechanisms proposed to overcome the final parsec problem. One of the most promising is tidal interaction of the binary with a circumbinary disk, naturally leading to the orbital decay of the binary, which loses its angular momentum to the disk (e.g., Haiman et al. 2009; Cuadra et al. 2009; Rafikov 2013). It is generally expected that a massive binary clears a cavity in a disk around itself so that the mass inflow onto its components is considerably suppressed compared to the mass accretion rate $\dot{M}_\infty$ far out in the disk. Mass accumulation at the inner edge of the disk caused by the tidal barrier modifies the disk structure over time and accelerates binary inspiral.

In this work we compute t_res in the framework of this gas-assisted binary inspiral scenario. If we assume that all quasars are triggered by gas accretion at ∼0.1 pc and the inspiral of the central SMBH binaries is dominated by the interaction with a gas disk, then the probability of observing a pair at a given radius is just proportional to the residence timet_res ≡ |dln r/dt)|⁻¹. This time is obviously a function of the physical mechanism driving the orbital evolution of the binary. As a result, our observations in principle provide a way of constraining different inspiral mechanisms.

Our calculation of t_res closely follows that in Rafikov (2013), which provides a recipe for computing the fully coupled, time-dependent binary+disk evolution and should be referred to for details. We start binaries of all masses at a semi-major axis of 0.1 pc, in a circumbinary disk with viscosity characterized by α = 0.1. We assume that in a radiation-pressure-dominated regime (which is valid at small separations) viscosity in the disk is proportional to the total (rather than the gas) pressure. The value of the mass accretion rate in the disk far from the binary $\dot{M}_\infty$ is normalized by the Eddington rate $\dot{M}_{\rm Edd}$ computed for the full binary mass and radiative efficiency of ε = 0.1. Evolutionary calculations are always started with a constant $\dot{M}=\dot{M}_\infty$ disk extending all the way into the initial binary orbit at 0.1 pc. Such an initial setup naturally arises if the gas that has fallen into the center of the galaxy circularizes outside the binary orbit and then viscously spreads inward. For simplicity our calculations ignore the issue of gravitational stability of the circumbinary disk at large separations, which is poorly understood at the moment (Goodman 2003; Rafikov 2009).

Rafikov (2013) has shown that the torque acting on the binary and driving its inspiral is set by $\dot{M}_\infty$, the radius of influence r_infl out to which the binary torque has been propagated by the viscous stresses in the disk, and the degree of mass overflow across the orbit of the secondary, i.e., the fraction χ < 1 of $\dot{M}_\infty$ that penetrates into the cavity and gets accreted by the binary components. The dependence on the latter is rather weak and the orbital evolution is similar to the no-overflow case as long as χ ≲ 0.2, which we assume to be the case in this work for simplicity. Of course, in practice some overflow is needed to power the observed quasar activity, so χ is not strictly zero.

The time dependence of the radius of influence is self-consistently computed by following the viscous evolution of the circumbinary disk torqued by the binary, as described in Rafikov (2013). In calculating r_infl we keep track of the different physical regimes in the disk at r_infl. This is important since r_infl can vary by orders of magnitude with respect to the binary semi-major axis in the course of its inspiral, allowing the opacity behavior and relative role of radiation pressure at r_infl to change over time. This self-consistent calculation favorably distinguishes our present calculation (and that in Rafikov 2013) from the self-similar solution of Ivanov et al. (1999) or the SMBH inspiral calculations in Haiman et al. (2009).

One particular aspect of the t_res calculation that we wish to emphasize is that in many cases we find the system in so-called disk-dominated evolution, a regime when the "local disk mass"⁴M_d = Σr² at the inner edge of the disk (here Σ is the local surface density of the disk enhanced by the mass pileup near the cavity edge) is larger than the mass of the secondary M_s. This regime is often valid at the initial stages of evolution, especially for lower SMBH masses, because of the large starting semi-major axis of 0.1 pc that we adopt in this work. In this regime we follow conventional wisdom (Ward 1997; Haiman et al. 2009) and simply assume that the secondary passively follows viscous evolution of the disk, which is equivalent to setting t_res = t_ν, where t_ν = r²/ν is the viscous time at the inner edge of the circumbinary disk, computed using the local value of viscosity ν and accounting for the possibility of mass pileup at the inner disk edge. As the binary orbit shrinks, the system inevitably transitions into the secondary-dominated regime, when M_s > M_d. As soon as this happens, we let the disk respond viscously to the binary torque. On the other hand, Rafikov (2013) has shown that deep in the secondary-dominated regime (M_d ≪ M_s) the evolution (or residence) time is given by

$$\begin{eqnarray} &&t_{\rm res}=t_\nu \frac{1}{6\pi (1+q)}\frac{M_s}{M_d}. \end{eqnarray} \tag{ 9 }$$

Thus, t_res evaluated according to this formula becomes equal to t_ν only when M_d = M_s/[6π(1 + q)], which is significantly less than M_s. The precise details of the transition between the disk- and secondary-dominated regimes are not well understood at the moment, so in this work we assume for simplicity that t_res = t_ν for M_d > M_s/[6π(1 + q)] and switches to the behavior predicted by Equation (9) only for M_d < M_s/[6π(1 + q)]. This assumption makes the viscous evolution of the disk quite important and can potentially result in an overestimate of t_res at intermediate separations (∼0.01–0.1 pc) by a factor of several. We cannot be more accurate in this regime without better understanding the regime when M_d ∼ M_s.

In Figure 7 we show maps of t_res computed according to the aforementioned prescription (Rafikov 2013) for several values of q and $\dot{M}_\infty /\dot{M}_{\rm Edd}$. Several features of these maps are worth discussing. At high masses (and small separations) color contours turn into straight lines, which reflects the fact that in this regime the binary inspiral is determined by GW emission. In this case t_res is simply given by

$$\begin{eqnarray} &&t_{\rm GW}(r_b)=\frac{5(1+q)^2}{8q}\frac{R_S}{c} \left(\frac{r_b}{R_S}\right)^4, \end{eqnarray} \tag{ 10 }$$

where R_S ≡ 2GM_BH/c² is the Schwarzschild radius of the BH with the combined mass M_BH. This timescale is very sensitive to M_BH and r, as can be seen in Figure 7. However, the quasars used in our study almost never probe this part of the parameter space because the size of the BLR R_BLR (shown by the dot-dashed line in t_res maps) exceeds the maximum separation at which GW emission dominates for all masses. Thus, our quasars are in the regime when the disk torque dominates binary inspiral.

Zoom In Zoom Out Reset image size

Figure 7 also reveals some obvious trends with q and $\dot{m}_{\rm Edd}$. Lower q results in faster inspiral since for a lower-mass secondary the binary contains less angular momentum but the disk torque is relatively insensitive to q (see Rafikov 2013). Higher $\dot{m}_{\rm Edd}$ at a fixed M_BH implies a circumbinary disk with higher surface density, which provides more torque, considerably accelerating the orbital evolution of the binary. Higher $\dot{m}_{\rm Edd}$ also extends the disk-dominated phase of the binary evolution. In particular, at low masses (and large separations) t_res exhibits rather weak dependence on M_BH, which is a result of the system evolving in a disk-dominated regime where the inspiral time is given by the local viscous timescale, which is only weakly dependent on M_BH.

Quite important for our work is the relatively weak dependence of t_res on separation in some cases, e.g., low $\dot{m}_{\rm Edd}$ and relatively high q, which results in rather long t_res for r ∼ 0.003–0.1 pc and increases the number of binaries in this range of r. This behavior is a result of two main factors. First, initially, when the system is in the disk-dominated regime, the inner parts of the disk around the binary are gas pressure dominated. Under these conditions the viscous timescale on which the binary orbit evolves in the disk-dominated case is a relatively weak function of r compared to the radiation-pressure-dominated case typical for high $\dot{m}_{\rm Edd}$ (Rafikov 2013), so that t_res is quite long at $r\mathrel {{\rlap{{{\sim }}}{{>}}}}0.01$ pc.

Second, even after evolution switches to the secondary-dominated regime, binary inspiral is slow and is still not very sensitive to the separation (e.g., see the behavior of t_res in Figure 8 of Rafikov (2013)). It is very important that here we self-consistently follow the time-dependent structure of the whole disk in the spirit of Ivanov et al. (1999). Previously, Haiman et al. (2009) found considerably faster binary evolution (lower t_res) for the same M_BH, r, etc., when employing the quasi-steady-state circumbinary disk solutions derived in Syer & Clarke (1995), which do not apply to real disks; see Ivanov et al. (1999) and Rafikov (2013) for details. The non-locality of the binary–disk coupling in the secondary-dominated regime (Rafikov 2013)—the fact that for a given value of $\dot{M}$ at the inner disk edge (equal to zero in the case of no overflow) the magnitude of the torque acting on the binary is set by the disk properties not at its inner edge, which is often radiation pressure dominated, but much farther out, at r_infl, where the disk is gas pressure dominated—is also very important for determining the residence time.

In computing observability, we consider two basic scenarios for the BH binary: no temporal evolution, and radial migration via a gaseous accretion disk. From each, combined with our upper limit of seven candidates, we can infer something about the BH binary population in luminous z ≈ 1–2 QSOs. In the first scenario, all SMBH binaries live at fixed separation with no time evolution as we discussed in Section 5.1. The expected numbers of observable sources from our high-S/N subsample versus binary separation is shown in Figure 6. Comparison of this theoretical distribution with our observational upper limit (N ⩽ 7) rules out that all the ∼10⁹ M_☉ BHs that we observe at z ≈ 1–2 exist in binaries in the separation range of ∼[0.03, 0.2] pc. Binaries closer than 0.03 pc are not observable due to the truncation or destruction of the BLR, while binaries with r > 0.2 pc are not observable given the Δt distribution of our sample. We note that even this basic result is of cosmological interest. For instance, in the merger-tree models of Volonteri et al. (2003), many major mergers are expected in the massive halos that host our massive BHs between 1 < z < 2. We can rule out that the lifetime of the resulting tight binaries is comparable to the Hubble time in all cases (Section 7).

In the second scenario, the binaries evolve through a gas disk. With the residence time $t_{\rm res}(r, M_{\rm BH}, \dot{m}_{\rm Edd})$ calculated in the previous section (Section 5.2, Figure 7), the probability $P_{{\rm obs}}(M_{\rm BH}, q, \Delta t, \dot{m}_{\rm Edd})$ is evaluated according to Equation (2) for each object in our sample using its observed M_BH and Δt and assuming that all of the objects have the same mass ratio q. The sum of P_obs over the objects is then the total expected number of observable sources (N_obs). As we discussed in Section 2, the virial masses of the secondary BHs shown in Figure 1(d) imply an $\dot{m}_{\rm Edd}$ distribution peaked at ∼0.25. For the first calculation of N_obs, we use t_res in the case that $\dot{m}_{\rm Edd} \sim 0.1$ and the virial BH masses. Then, to bracket the known large uncertainties in M_BH measurements (e.g., Vestergaard & Peterson 2006; Shen et al. 2008), we also do the calculations assuming that $\dot{m}_{\rm Edd}=1$ and adopting M_BH accordingly. In both of the calculations, we truncate the binary migration at R_BLR for each object according to Equation (2), using the virial and Eddington-limited BH masses, respectively. We show the range of expected N_obs in our S/N pixel⁻¹ >10 subsample for $\dot{m}_{\rm Edd}=0.1, 1$ and q = 1, 0.1 in Table 3.

Table 3. Expected Number of Observable Sources in the S/N pixel⁻¹ > 10 Subsample

	$\dot{m}_{\rm Edd}=0.1$	$\dot{m}_{\rm Edd}=1.0$
q = 1.0	44	0.34
q = 0.1	282	81

Notes. The four sets of numbers show the expected number of observable sources in our S/N pixel⁻¹ > 10 subsample, adopting four sets of ($q, \dot{m}_{\rm Edd}$) in the scenario of gas-assisted sub-parsec evolution of the SMBH binaries. We take the virial BH masses for the $\dot{m}_{\rm Edd}=0.1$ cases and Eddington-limited BH masses for the $\dot{m}_{\rm Edd}=1$ cases. For comparison, in the no-evolution case with r = 0.1 pc, we expect 40–98 sources for q = 1 and 498–655 sources for the unrealistic q = 0.1 case.

Download table as:  ASCIITypeset image

If we adopt the most likely parameters $\dot{m}_{\rm Edd}=0.1$ and q = 1, the final expected number of SMBH binaries in the high-S/N subsample is N_obs = 44. Comparison of this expectation with our seven candidates implies that <16% of quasars host SMBH binaries with r < 0.1 pc. A more extreme mass ratio of q = 0.1, which seems unlikely for the bulk of our sample, increases the expected N_obs by a factor of ∼6.4.

Note, however, that our results are very sensitive to the extrapolated R_BLR–L_bol scaling relation, which is highly uncertain. Since our method is only sensitive to binaries that are less than ∼0.1 pc apart (see Figure 6), the final expected number of observable candidates is actually determined by the detectability within the range [R_BLR, 0.1 pc]. Candidates in this range are the most massive BHs with the corresponding largest velocity drifts, so their BLR size is approaching ∼0.1 pc (e.g., R_BLR ∼ 0.07 pc for 10⁹ M_☉ assuming $\dot{m}_{\rm Edd}=0.1$). A slightly larger R_BLR makes the valid range [R_BLR, 0.1 pc] much narrower. We will test the sensitivity to R_BLR quantitatively in Section 6.

The sensitivity of our results to $\dot{m}_{\rm Edd}$ has two origins. First, the BH masses we use scale with $\dot{m}_{\rm Edd}$ given the observed luminosity. Second, the gas-assisted evolution of BH binaries depends on $\dot{m}_{\rm Edd}$: in the lower $\dot{m}_{\rm Edd}$ case the circumbinary disk is less massive, meaning (1) more significant dominance of the secondary BH over the local disk (larger M_s/M_d) and (2) less torque acting on the binary (which is proportional to $\dot{M}_\infty \propto \dot{m}_{\rm Edd}$). These factors increase t_res for R_BLR < r < 0.1 pc in our scenario of the gas-assisted evolution, boosting the expected N_obs for lower $\dot{m}_{\rm Edd}$.

Our results are also very sensitive to the details of the gas-assisted evolution model in the intermediate separation range [R_BLR, 0.1 pc]. Evolutionary models where BH binaries spend a longer time in this radial range predict a larger number of observable sources N_obs. In our calculation of t_res following Rafikov (2013) in Section 5.2, this separation range falls in the secondary-dominated regime, in which matter piles up outside the binary orbit and binary evolution is slow and is described by Equation (9). Faster orbital evolution in this separation range, as suggested by Haiman et al. (2009) based on the circumbinary disk solution by Syer & Clarke (1995), would result in a smaller number of detectable SMBH binaries. Thus, it should in principle be possible in the future to discriminate between the different evolutionary scenarios.

We restate that all the results in this study are based on the assumption of circular SMBH binary orbits for simplicity. The effects of eccentric orbits on SMBH binary detection are not certain due to the uncertainty of BLR physics. For an eccentric binary orbit (1) the individual BLRs around the SMBHs would either be truncated earlier during inspiral and thus more compact, making the velocity shifts easier to detect, or be tidally destroyed earlier, making a detection more challenging; and (2) larger velocity drifts characteristic for eccentric orbits would in principle allow us to sample a broader mass range of SMBHs relative to the circular orbits.

Another important mechanism that may solve the final parsec problem is that stars in a triaxial galaxy nucleus could continuously replenish the loss cone and drive subsequent inspiral evolution (Merritt & Poon 2004; Merritt & Vasiliev 2011). Merritt & Poon (2004) found that stars on chaotic (centrophilic) orbits in a triaxial galactic nucleus could refill the loss cone and feed the central BHs at quite a high rate and thus drive binary evolution at sub-parsec scales. Merritt & Vasiliev (2011) found that there exists a family of stellar orbits in triaxial galaxies called pyramid orbits. The angular momentum of stars can reach zero at the "corners" of the pyramid, at which point they are easily captured by the central BHs. The lifetime of the binary then becomes a function of the structure of the merged remnant (Khan et al. 2012) but may take billions of years. In this scenario, quasars could be ignited at any radius, since the lifetime of the quasar is likely much shorter than the inspiral timescale of the BH binaries. However, without more assumptions about the dynamical state of the galaxy centers to determine t_res, we cannot easily predict the detectability of SMBH binaries in this scenario.

In the most pessimistic evolutionary scenario, the SMBH binaries stop spiraling inward due to depletion of stars and stall there forever. Quinlan (1996) did scattering experiments and found that the orbital decay of the binary stalls when most of the stars have been ejected, and this stalling separation is 0.01–1 pc depending on the mass and mass ratio of the BH binary. The stalling separation we extract from Quinlan (1996, see their Figure 7) is

$$\begin{equation} \log _{10}\left(\frac{R_{\rm stall}}{\rm pc}\right) = 0.082 \log _{10}\left(\frac{M_{\rm BH}}{M_\odot }\right) + 0.448 \log _{10} q - 0.474. \end{equation} \tag{ 11 }$$

For major mergers (q ∼ 1), the stalling radius is around ∼1 pc. If there is no subsequent driving mechanism and all binary BHs stall at ∼1 pc, our experiment would not be sensitive to any binary BHs according to Figure 6. With time baselines of only a few years, as probed here, we cannot hope to see radial velocity changes in parsec-scale binaries. However, as we will see in the next section, next-generation surveys could significantly increase the observability even in this scenario due to longer observing time intervals.