Determining Subparsec Supermassive Black Hole Binary Orbits with Infrared Interferometry

Central supermassive black holes in merging galaxies are thought to be efficiently driven to ≲10 pc separations by dynamical friction (Begelman et al. 1980). Their further evolution remains uncertain. Interactions with gas in a circumbinary accretion disk could either drive the binary closer together (Armitage & Natarajan 2002) or farther apart (e.g., Muñoz et al. 2019). Detections of subparsec supermassive black hole binaries (SMBHBs) would provide important input to galaxy formation models (Volonteri et al. 2003), estimates of the stochastic gravitational-wave background (e.g., Siemens et al. 2013), and the rate of individual merger events seen by the Laser Interferometer Space Antenna (LISA; e.g., Amaro-Seoane et al. 2012).

Growing numbers of dual active galactic nuclei (AGN) are seen on kiloparsec scales in interacting or postmerger galaxies (Comerford et al. 2009). The closest known supermassive black hole pair has a projected separation of ≃7 pc (Rodriguez et al. 2006), detected with radio very long baseline interferometry. Suggested evidence of subparsec binaries comes from AGN with double-peaked broad emission lines (Gaskell 1983), offset and moving broad emission lines (Eracleous et al. 2012), and periodically varying optical light curves (Graham et al. 2015).

Infrared interferometry with the Very Large Telescope Interferometer (VLTI) instrument GRAVITY (Gravity Collaboration et al. 2017) can now spatially resolve the broad emission line region (BLR) in the brightest AGN on sky by measuring its velocity-dependent photocenter offset from the hot dust continuum (Gravity Collaboration et al. 2018). For a system with double-peaked broad lines, an extension of this method could reveal the presence of an SMBHB (Songsheng et al. 2019). Several candidate double-peaked systems have been ruled out as binaries (Eracleous et al. 1997; Decarli et al. 2013), and both black holes are only expected to be actively accreting and retain their individual BLRs over a narrow region of parameter space (Bogdanović et al. 2008; Shen & Loeb 2010).

Monitoring campaigns have identified a number of candidates with single-peaked, offset, and moving emission lines (Runnoe et al. 2017; Guo et al. 2019). Here we consider the requirements for astrometrically confirming the presence of a binary in these systems. Over a relevant range of parameter space, relative astrometry between the BLR of an accreting secondary black hole and hot dust in the surrounding circumbinary disk could map out the binary orbit (Section 2). The observational requirements, given the current candidate systems, are well matched to the sensitivity of the planned upgrade of the GRAVITY instrument, GRAVITY+ (Section 3). A monitoring campaign over ≃5–10 yr could be sufficient to detect both proper motion and acceleration in these systems, constraining the system parameters and potentially providing robust detections of subparsec SMBHBs. Possible extensions of this study including the prospects of additional measurements and targets are discussed in Section 4.

We assume a binary system of total mass M = M₁ + M₂ and mass ratio of q = M₂/M₁ ≤ 1 in a circular orbit. The orbital period and semimajor axis on sky are then

$$\begin{eqnarray}&&P\simeq 95{\left(\displaystyle \frac{a}{0.1\mathrm{pc}}\right)}^{3/2}{\left(\displaystyle \frac{M}{{10}^{9}{M}_{\odot }}\right)}^{-1/2}\,\mathrm{yr},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\theta }_{a}\simeq 34\left(\displaystyle \frac{a}{0.1\,\mathrm{pc}}\right){\left(\displaystyle \frac{{D}_{A}}{600\mathrm{Mpc}}\right)}^{-1}\,\mu \mathrm{as},\end{eqnarray} \tag{ 2 }$$

where D_A is the source angular diameter distance. We further assume that the SMBHB is surrounded by a circumbinary gas disk, which is centered on the system center of mass and truncated at a radius ≃2a (Artymowicz & Lubow 1994). Accretion proceeds through a central, low-density cavity via thin streams, forming "mini-disks" around the two black holes (e.g., Cuadra et al. 2009; Noble et al. 2012; D'Orazio et al. 2013; Bowen et al. 2018).

2.1. Relevant Parameter Regime

Mapping out the binary orbit requires an astrometric measurement of a light source centered on one of the black holes. With near-infrared observations, the most promising candidate is a broad emission line from ionized gas bound to one of the black holes. For concreteness, we assume that this is the secondary black hole M₂. Many simulations find a much higher accretion rate onto the secondary (e.g., Cuadra et al. 2009; D'Orazio et al. 2013; Muñoz et al. 2019; Duffell et al. 2020). This is also the assumption made by recent radial-velocity studies (e.g., Runnoe et al. 2017), allowing for a direct comparison. With infrared interferometry, we also need a reference source. Here we consider the method in current use, where the broad emission line is phase-referenced to the continuum emission radiated by the surrounding hot dust.

Two requirements to make this measurement are that (1) the BLR is bound to the secondary black hole (R_BLR < R_L, where R_L is the Roche-lobe radius as approximated by Eggleton 1983), while (2) hot dust is not (R_sub > R_L, where R_sub is the sublimation radius). We estimate R_BLR ≃ 0.07L_2,46 pc and R_sub ≃ 0.4L_2,46 pc using scaling relations with luminosity as measured separately for the BLR (e.g., Bentz et al. 2013) and near-infrared continuum (Suganuma et al. 2006; Kishimoto et al. 2011; Gravity Collaboration et al. 2020a). The luminosity of the secondary black hole is L₂ = L_Edd (qM), where L_Edd is the Eddington luminosity and  = 0.1 is the assumed Eddington ratio of the secondary. Any viable candidates identified by the radial-velocity method would by definition have a BLR bound to the secondary. Even large graphite grains, often assumed responsible for the NIR continuum (e.g., Kishimoto et al. 2007), should be sublimated within the Roche lobe of the secondary for binary orbital periods of ≲10³ yr (below the gray lines in Figure 3).

The major uncertainty in this scenario is where the near-infrared continuum emission originates, and how it evolves over the course of a binary orbit. We consider two scenarios (Figure 1). (i) If the continuum emission is stationary, e.g., tracing the inner edge of the circumbinary disk, then relative astrometry of the BLR measures the secondary's orbit. (ii) Empirically, the near-infrared emission size scales with that expected for the sublimation radius. It seems possible that the continuum emission could instead preferentially originate in the regions of the circumbinary disk closest to the secondary, where the irradiating flux is strongest and dust temperatures highest. In that case, both the line and continuum emission could track the binary orbit, although we have not tested this using radiative transfer calculations including dust heating, anisotropic emission, or obscuration along the line of sight.

Zoom In Zoom Out Reset image size

We have developed a simple geometric model for the second "evolving continuum" scenario. Hot dust is assumed to form outside the binary and at the sublimation radius of the secondary. The possible emission locations are then along a circle of radius R_sub centered on the position of the secondary. When the sublimation radius intersects the circumbinary disk, we assume that hot dust emission is produced with equal intensity everywhere along the circle where it intersects the circumbinary disk. When the sublimation radius is smaller than the distance from the secondary to the edge of the circumbinary disk, we assume that some small region (e.g., in an accretion stream) at a distance of ≃R_sub will form and radiate hot dust instead. The astrometric shift is then the offset between the secondary black hole and the continuum photocenter.

The expression is derived in the Appendix and the result is shown in Figure 2. At very small R_sub/a the offset is small because hot dust forms close to the secondary. Once R_sub becomes large enough to heat dust all along the circumbinary disk, the continuum photocenter is at the position of the secondary black hole, and the astrometric shift vanishes. For a range of 0.5 ≲ R_sub/a ≲ 2.5, the relative offset is similar in magnitude to the true orbital offset. We plot this parameter space constraint as the dark red lines in Figure 3. It is more restrictive than simply requiring that hot dust cannot be bound to the secondary. In particular, for q ≪ 1 the available parameter space shrinks until a minimum q_min ≃ 6 × 10⁻³ where no solutions are possible. Still, the geometric model suggests that relative astrometry might trace the binary orbit over much of the relevant parameter space, even if the near-infrared continuum is tracking the motion of the secondary.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.2. Supermassive Black Hole Binary Astrometry

We next consider the radial velocity and astrometric position of the secondary black hole on sky. Following Eracleous et al. (2012), we write the radial velocity as

$$\begin{eqnarray}&&{u}_{2}(t)={\left(\displaystyle \frac{2\pi G\tilde{m}}{P}\right)}^{1/3}\sin i\sin \phi (t).\end{eqnarray} \tag{ 3 }$$

where P is the orbital period, i the inclination, and ϕ(t) = 2πt/P + ϕ₀ where t is the current time and ϕ₀ is the orbital phase. For a position angle (PA) on sky measured E of N, the astrometric positions are

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{x}}(t) & = & {\left(\displaystyle \frac{G\tilde{m}{P}^{2}}{4{\pi }^{2}}\right)}^{1/3}\\ & & \times \ \left[\begin{array}{c}-\cos \,\mathrm{PAs}\,\mathrm{in}\phi (t)-\cos i\,\sin \,\mathrm{PA}\,\cos \phi (t)\\ \sin \,\mathrm{PA}\,\sin \phi (t)-\cos i\,\cos \,\mathrm{PA}\,\cos \phi (t)\end{array}\right].\end{array}\end{eqnarray} \tag{ 4 }$$

With only observations of the secondary's motion, the measurable combination of masses is $\tilde{m}=M/{(1+q)}^{3}$, resulting in a factor of 8 range in allowed total mass M. Assuming the hot dust emission is centered on the binary center of mass, a single measurement of the offset (x, y) provides a lower limit to the semimajor axis a on sky. The astrometric offset should be large when the radial-velocity offset is near maximum, as selected by Eracleous et al. (2012). A proper motion measurement can be compared with the radial-velocity offset, and a second derivative of either quantity measures the orbital period P. Combining positions and proper motions with radial-velocity measurements provides enough information to constrain an orbit.

If the hot dust emission is stationary (e.g., uniform or asymmetric around the circumbinary disk), Δx(t) = x(t) + x₀ would be the measured quantity, with x₀ a potentially constant offset of the dust emission. If instead the hot dust emission follows the motion of the secondary as in the geometric model above, then Δx(t) ≃ −(0.5–1.5) (1 + q) x(t). An unknown prefactor would produce additional scatter by a factor of ≃10 in the inferred value of M, but with weak dependence on q.

Currently known candidate SMBHBs with single, offset, moving broad emission lines are generally found at z ≃ 0.2, with apparent magnitudes of V ≲ 18 and K ≲ 15 (Runnoe et al. 2017; Guo et al. 2019). For a semimajor axis of a ≃ 0.1 pc, the size on sky θ_a ≃ 30 μas, while the BLR size is a factor of several smaller. These properties are well matched to the expected sensitivity of the planned upgrade to the GRAVITY instrument, GRAVITY+. Through a combination of ongoing and near-future upgrades including new grisms, improved VLTI vibration control, new adaptive optics systems, and laser guide stars the goal is to reach limiting magnitudes K ≲ 14–15 with comparable astrometric accuracy as is currently possible for K ≲ 10–11.⁶

3.1. Differential Phase Astrometry

The astrometric offset of an emission line of strength 1 + f relative to the normalized hot dust continuum is measured by the differential phase Δϕ = ϕ(λ) − ϕ_c,

$$\begin{eqnarray}&&{\rm{\Delta }}\phi =-2\pi \displaystyle \frac{f}{1+f}(u{\rm{\Delta }}x+v{\rm{\Delta }}y)\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&| {\rm{\Delta }}\phi | \simeq 0\buildrel{\circ}\over{.} 3\left(\displaystyle \frac{{f}_{\mathrm{line}}}{0.1}\right)\left(\displaystyle \frac{a}{0.1\,\mathrm{pc}}\right){\left(\displaystyle \frac{{D}_{A}(z)}{500\mathrm{Mpc}}\right)}^{-1}{\left(1+q\right)}^{-1},\end{eqnarray} \tag{ 6 }$$

with Δx(t) and Δy(t) the astrometric offsets discussed above. The line strength f is normalized to the continuum flux, f_line = f/(1 + f), and D_A(z) is the angular diameter distance. The differential phase signal of a wavelength-independent (x, y) offset has the shape of the emission line itself, with an amplitude depending on the (u,v) coordinates of each baseline.

3.2. A Case Study with SDSS J1402+2631

From the parent radial-velocity samples of Runnoe et al. (2017) and Guo et al. (2019), we have listed properties of some SMBHB candidate targets visible from the VLTI (Dec < 30°) with K < 15 and 0.09 < z < 0.25 in Table 1. For those redshifts, the Pa α line is redshifted into the GRAVITY K band. All 10 targets have predicted phase signatures of ≳03 for a 0.1 pc binary orbit. As such they form a promising first set of candidates for GRAVITY+ astrometry.

Table 1.  Some Candidate SMBHB GRAVITY+ Targets

SDSS ID	z	K	V	θa,1 (μas)	Δϕ (deg)	amin (pc)	amax (pc)	References
SDSS J001224.02-102226.2	0.2287	13.7	17.1	27.0	0.28	0.10	0.43	1
SDSS J015530.01-085704.0	0.1648	12.7	16.8	35.0	0.37	0.08	0.33	1
SDSS J091928.69+143202.6	0.2072	14.5	17.6	29.1	0.31	0.07	0.30	1
SDSS J093844.45+005715.7	0.1707	13.8	17.2	34.0	0.36	0.07	0.29	1
SDSS J111230.90+181311.4	0.1952	14.5	18.4	30.5	0.32	0.04	0.19	2
SDSS J115158.90+122128.9	0.1697	14.5	17.9	34.2	0.36	0.05	0.21	1
SDSS J125142.28+240435.3	0.1887	14.0	17.6	31.4	0.33	0.06	0.26	1
SDSS J140251.19+263117.5	0.1877	12.5	16.9	31.5	0.33	0.09	0.37	1
SDSS J153705.95+005522.8	0.1365	13.5	17.3	40.9	0.43	0.05	0.22	2
SDSS J155654.47+253233.5	0.1645	13.9	18.0	35.0	0.37	0.04	0.19	1

Note. Targets are selected as those with K < 15, decl. > 30°, and 0.09 < z < 0.25 from the offset radial-velocity SMBHB candidates identified by (1) Runnoe et al. (2017) and (2) Guo et al. (2019). The estimated astrometric size θ_a,1 is scaled to a semimajor axis of 0.1 pc using angular diameter distances from the target redshifts. The phase signal is calculated according to Equation (5) assuming a Pa α line strength of f_line = 0.1 and q = 0.1. The allowed range of semimajor axis for astrometric measurements is inferred from the optical luminosity as described in the text.

Download table as:  ASCIITypeset image

We have further used the observed optical luminosity, radial-velocity offset, and minimum periods for the sample to constrain the parameter space where astrometric monitoring might be feasible. Following Section 2, we calculate the allowed range of semimajor axis from R_BLR = R_L(q, a_min) and a_max = 2R_dust. As shown in Table 1, we are sensitive to binary semimajor axes of ≃0.05–0.4 pc. This range depends on the mass ratio q, in the sense that a_min increases with decreasing q. The full range is feasible for nearly equal mass binaries with q ≲ 1. We can impose further constraints to estimate allowed total binary mass ranges. We require a total binary mass that (1) results in P_min < P < 10³ yr, where P_min is the minimum period obtained from fitting the measured radial-velocity curves (Runnoe et al. 2017), (2) can match the observed radial-velocity offset u₂ (Equation (3)), and (3) results in an Eddington ratio of 10⁻³ < L₂/L_Edd < 3 for the secondary. All of those constraints are satisfied for total masses of M ∼ 10^7–10M_⊙.

As one example, we consider the object SDSS J1402+2631. We have measured the Pa α emission line profile of this quasar (Figure 4) using the TripleSpec instrument at the Apache Point Observatory 3.5 m telescope. Observations were taken in 2020 June with the 11 slit in a standard nodding ABBA sequence of 8 × 120s exposures. The seeing was 1''. The data were reduced using a modified version of the Spextool package (Cushing et al. 2004), and an A0V star was used for telluric correction (Vacca et al. 2003). We detect broad emission lines of Pa α, β, γ, δ, at a redshift of z ≃ 0.188 in the JHK band spectra. The continuum flux corresponds to K = 12.8, similar to the K = 12.5 measured by 2MASS. Figure 4 shows a decomposition of the Pa α emission line into Gaussian broad and narrow components, where the broad-line component has a velocity width σ ≃ 3300 km s⁻¹ and peak relative line strength of f ≃ 0.12. The line width is consistent with the reported range of Hβ FWHM (Runnoe et al. 2015).

Zoom In Zoom Out Reset image size

We use the broad-line component model to simulate GRAVITY+ data, adopting a phase error of 01 per baseline as achieved in observations of bright (K ∼ 10–11) AGN to date with GRAVITY (Gravity Collaboration et al. 2018; GRAVITY Collaboration et al. 2020b). We take VLTI (u, v) coordinates of this northern target from Aspro (Bourgès et al. 2013). The top right panel of Figure 4 compares the measured line profile and simulated differential phase signals for fiducial parameters of $\tilde{m}={10}^{9}\,{M}_{\mathrm{sun}}$, a = 0.1 pc, P = 100 yr, using the model described in Equations (3) and (4) and assuming a stationary continuum photocenter. The differential phase is averaged over the three longest baselines. Fitting Equation (5) for the offset (x, y) results in errors of ≃2 × 4 μas. The measured offsets and errors are shown compared to the underlying model in the bottom panel of Figure 4. Both proper motion and acceleration would be detected from astrometric monitoring, resulting in confirmation of the target as an SMBHB and allowing estimates of $\tilde{m}$ and P, in combination with radial-velocity measurements. The parameters ϕ₀ and PA are difficult to constrain, likely due to the low inclination angle.

3.3. Mass and Period Estimates from a Parameter Survey

We next perform a mock parameter survey to see how well binary mass and period information might be recovered. We consider periods of 3 × 10^2–4 yr and $\tilde{m}=3\times {10}^{7-9}{M}_{\odot }$. We generate 10 epochs of simulated radial-velocity data taken over a 25 yr time baseline (since current candidates have ≲15 yr time baselines) with errors of 100 km s⁻¹ intended to mimic the "jitter" noise that dominates the error budget in many current candidates (Runnoe et al. 2017). We generate 10 epochs of astrometric data over 8 yr, adopting errors of 4 μas in both the x and y (R.A. and decl.) coordinates.

For each combination of P and $\tilde{m}$, we generate N = 300 realizations of mock data, varying the random error realization as well as the parameters of i, ϕ₀, and PA. The inclination is constrained to be i < 75°, while ϕ₀ and PA are varied over their full ranges. We use a least-squares method to identify the best-fitting parameters in each case. The initial guess for least squares is fixed to fiducial values of $\tilde{m}={10}^{8.5}\,{M}_{\odot }$ and P = 300 yr. The median parameter bias and scatter over the N = 300 simulations for each parameter combination are shown in Figure 5, excluding the ≃2% of simulations where the minimization method fails. We recover the input parameters with errors of ≲0.6 dex for periods of P ≲ 10³ yr. For longer periods, second derivatives are usually not detected in radial velocity or astrometry. For P ≳ 300 yr and $\tilde{m}\gtrsim {10}^{8}{M}_{\odot }$, the recovered parameters show bias, in that they systematically find shorter periods and smaller $\tilde{m}$ than the input values.

Zoom In Zoom Out Reset image size

Current subparsec SMBHB candidates with single, moving, broad emission lines have K ∼ 12–15 and sizes on sky of θ_a ≃ 30 μas (a/0.1 pc). It may be possible to trace the binary orbit in these systems with the upgraded near-infrared interferometry instrument GRAVITY+ at the VLTI. A monitoring campaign over 5–10 yr could reveal proper motions and accelerations, resulting in robust detections of the progenitors of merging supermassive black holes and constraining their system parameters.

As an example, using the Pa α profile of one current candidate and current GRAVITY phase noise we find astrometric errors of ≲4 μas. We simulate a combined radial velocity and astrometric campaign, resulting in robust detections of binaries with ≲0.5 dex measurements of $\tilde{m}$ and P for systems with P ≲ 10³ yr where accelerations can be measured.

With radial-velocity data alone, generally P can still be well constrained, since radial-velocity changes (accelerations) can usually be measured over our assumed 25 yr of monitoring. Constraining $\tilde{m}$ requires astrometry. We also note that the complicating issues of line profile changes and jitter noise would not impact the astrometric offset measurement. The differential phase signal is proportional to the ratio of line to total flux, even for a variable line profile.

In principle, combining astrometric and radial-velocity data we can fit for the angular diameter distance D_A without using the redshift. The result would then provide a cosmological constraint. As expected, fitting directly for the distance results in a strong correlation between $\tilde{m}$ and D_A, while P remains well measured. In our tests, precise measurements of both $\tilde{m}$ and D_A require astrometric errors of ≲1 μas and/or astrometric campaigns of ≳25 yr. This may be feasible for short period systems, and/or if even higher astrometric precision becomes possible.

The differential phase measurement is referenced to the continuum photocenter position. The continuum near-infrared emission is due to hot dust, whose origin and time evolution in the SMBHB scenario is unclear. We have considered two extreme cases. In one case, the continuum is stationary with a photocenter at the center of mass of the binary. In this case, relative astrometry directly measures the orbital position of the secondary black hole. We have used this model to generate synthetic data above.

We have considered a simple geometric model of the second case, assuming that the hot dust emission originates in the circumbinary disk at the distance of the sublimation radius away from the secondary. In this case, the hot dust photocenter tracks the orbital motion of the binary. Remarkably, over a large portion of the relevant parameter space (Figure 3) the relative offset in this model is opposite in sign and comparable in amplitude to the orbital motion of the secondary (Figure 2). In the evolving continuum model, it is possible in principle to measure both $\tilde{m}$ and q, e.g., the two black hole masses M₁ and M₂. This seems to require lower measurement errors and longer campaigns than we have assumed.

A time-variable central luminosity will produce fluctuations of the hot dust photocenter due to differential light travel time delays (reverberation; Shen 2012; D'Orazio & Haiman 2017). For relatively small fluctuations, the maximum amplitude of this effect has comparable contributions from changes in the hot dust emission radius and intensity (Δx/a ≲ 10% each for ΔL/L ≃ 20% at i = 30°). We evaluate the possible impact of uncertainties in the hot dust structure and its time variability using experiments with fake data. We consider models with (i) a constant hot dust offset (e.g., due to asymmetry), (ii) a fluctuating hot dust offset due to luminosity variations of ΔL/L ≃ 20% using a measured R-band light curve of 3C 273 (Fan et al. 2014), and (iii) an evolving offset tracking the orbit according to the geometric model described above. In each case, we run 30 trials of fitting the static dust orbital model (with no continuum photocenter offset; Equation (4)) to the generated data and errors. Data are generated with $\tilde{m}={10}^{8.5}{M}_{\odot }$ and P = 100 yr, and errors of 100 km s⁻¹ in radial velocity and 4 μas in astrometry. As in Section 3.3, we identify the best-fitting parameters using a least-squares method. Distributions of the identified best-fitting P and $\tilde{m}$ are shown in Figure 6. For the constant and fluctuating offset cases, the mass parameter is overestimated. Depending on the choice of parameters, we have also found underestimates. For the evolving offset case, the mass parameter is well recovered, while the orbital period is overestimated. These biases are introduced by the use of an incorrect hot dust emission model. In all cases, proper motions and accelerations can still be detected.

Zoom In Zoom Out Reset image size

In the evolving dust scenario, the hot dust emission region size is smaller and concentrated on one side of the circumbinary disk. The amplitude of the reverberation offset will be smaller as a result. However, the light travel time delay will cause the offset vector Δx between the BLR and hot dust to point slightly away from the center of mass. In principle, current GRAVITY observations could detect both the fluctuating sublimation radius size and reverberation effect using differential amplitude and phase data (e.g., Gravity Collaboration et al. 2020a) from different epochs where the continuum luminosity varies.

The same interferometry measurements proposed here could help distinguish scenarios for the hot dust continuum emission and its time variability in SMBHB candidates. The evolving continuum model would generically predict a smaller size than the stationary dust model for R_sub ≲ 2a. Candidates in that regime should show larger (smaller) dust sizes than expected from the radius–luminosity relation according to the stationary (evolving) dust emission models. The evolving continuum model might also show time-variable, asymmetric structure. Further constraints on both hot dust and BLR evolution would be possible if more distant, narrow emission line components of Pa α or Si [VI] are present, since they could be used as independent, static phase references.

We have focused on targets with offset, moving broad emission lines and assumed that the broad emission line originates from atomic gas centered on the secondary black hole. In the model of Nguyen et al. (2020), the larger BLR size around the primary could result in substantial contributions from its own line flux. If most of the atomic line emission is from around the primary, the astrometric signals considered here will be suppressed by a factor of q/(1 + q), and interferometry measurements would be most sensitive to large mass ratios of q ≳ 1/3. Our simulations have also used circular binary orbits. The same measurements are in principle possible if the binary is driven to high eccentricity. Additional time variability of the accretion luminosity and circumbinary disk size and shape could result in larger fluctuations of the hot dust photocenter location in this case.

Photometric candidates showing sinusoidal optical variations (e.g., PG 1302−102; Graham et al. 2015) should also be sufficiently bright to detect with GRAVITY+. For the very short periods ≲30 yr accessible with photometric data to date, a single complex BLR structure might surround both black holes (e.g., Shen & Loeb 2010). The astrometric signature in that case is unclear. Songsheng et al. (2019) calculated velocity-dependent photocenter signatures of a binary system with two active black holes, each with its own BLR. They further assumed a continuum photocenter at the center of mass, and identical Eddington ratios for both black holes. Relaxing either of those assumptions (e.g., D'Orazio & Loeb 2019) would result in an additional velocity-independent astrometric offset like that discussed here. Kovacevic et al. (2020) presented a first exploration of the combined effects for a somewhat different parameter regime than explored here. Both an overall offset of the hot dust and BLR photocenters, and velocity-resolved kinematics of the BLR have been detected recently in IRAS 09149−6206 (GRAVITY Collaboration et al. 2020b), providing independent measurements of the photocenter offset and BLR size.

We identified 10 possible candidates, which show evolution in radial velocity consistent with binary motion in ≥3–5 epochs over 5–15 yr (Runnoe et al. 2017; Guo et al. 2019). Large spectroscopic surveys will likely add additional candidates in the next several years. For example, the SDSS-V black hole Mapper program plans to take between 3 and 13 spectra of each of 25,000 quasars (Kollmeier et al. 2017). Additional candidates in the southern sky would be particularly promising for GRAVITY+ observations, since deep integrations of ≃4–8h may be required to achieve the astrometric accuracy needed to confirm candidate systems as SMBHBs and map out their orbits.

J.D. thanks C. Gray, A. Kowalski, K. Davis, D. Dewitt, M. Jacobs, and J. Crowley for their help in obtaining the APO data used here. We thank T. Bogdanovic, D. D'Orazio, S. Hönig, Y. Shen, J. Runnoe, and the anonymous referee for helpful comments that improved this manuscript. J.D. was supported in part by NSF grant AST-1909711 and an Alfred P. Sloan Research Fellowship.

Consider two circles, one describing the inner edge of the circumbinary disk of radius 2a centered on the center of mass, and one with radius R_sub centered on the secondary black hole (the "sublimation ring"), offset (without loss of generality) in the −x direction by a distance a/(1 + q) (see the right panel of Figure 1). When the two circles intersect, we calculate the offset between the line and continuum emission as the centroid of the arc of the sublimation ring that intersects the circumbinary disk. The centroid of the arc is

$$\begin{eqnarray}&&{\rm{\Delta }}x=-{R}_{\mathrm{sub}}\displaystyle \frac{\sin \alpha }{\alpha },\end{eqnarray} \tag{ A1 }$$

where the offset Δx = x_BLR − x_dust is negative, and α is the half-angle of the arc,

$$\begin{eqnarray}&&\sin \alpha =\displaystyle \frac{a(1+q)}{2{R}_{\mathrm{sub}}}\xi (q,{R}_{\mathrm{sub}}/a),\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&\xi (q,r)=\sqrt{({q}_{o}^{2}-{r}^{2})({r}^{2}-{q}_{i}^{2})},\end{eqnarray} \tag{ A3 }$$

where q_i = (1 + 2q)/(1 + q) and q_o = (3 + 2q)/(1 + q) bound the range of solutions where the two circles intersect. To calculate the correct half-angle, we need to switch solutions at a transition point ${r}_{t}=\sqrt{{q}_{i}{q}_{o}}$ given by ∂[ξ(q, r)/2r]/∂r = 0, where r = R_sub/a. A piecewise expression for the offset is then

$$\begin{eqnarray}\begin{array}{l}{\rm{\Delta }}x\ =\ \left\{\begin{array}{ll}-{ar} & r\lt {q}_{i}\\ -a(1+q)\displaystyle \frac{\xi (q,r)/2}{\arcsin \left[(1+q)\xi (q,r)/2r\right]} & {q}_{i}\leqslant r\leqslant \sqrt{{q}_{i}{q}_{o}}\\ -a(1+q)\displaystyle \frac{\xi (q,r)/2}{\pi -\arcsin \left[(1+q)\xi (q,r)/2r\right]} & \sqrt{{q}_{i}{q}_{o}}\lt r\leqslant {q}_{o}\\ 0 & r\gt {q}_{o}.\end{array}\right.\end{array}\end{eqnarray} \tag{ A4 }$$

The expression can be written more compactly using arcsin z = zR_C(1 − z², 1) for −1 ≤ z ≤ 1, where R_C(x, y) is the Carlson (1979) circular function:

$$\begin{eqnarray}{\rm{\Delta }}x=\left\{\begin{array}{ll}-{ar} & r\lt {q}_{i}\\ -\displaystyle \frac{{ar}}{{R}_{C}(1-{z}^{2},1)} & {q}_{i}\leqslant r\leqslant \sqrt{{q}_{i}{q}_{o}}\\ -\displaystyle \frac{{arz}}{\pi -{{zR}}_{C}(1-{z}^{2},1)} & \sqrt{{q}_{i}{q}_{o}}\lt r\leqslant {q}_{o}\\ 0 & r\gt {q}_{o},\end{array}\right.\end{eqnarray} \tag{ A5 }$$

and z = (1 + q) ξ(q, r)/2r.

When r = q_i, z = 0 and R_C(1, 1) = 1 so that Δx = −ar. When r² = q_iq₀, z = 1 and R_C(0, 1) = π/2 and the solutions again match on smoothly. When r = q₀, z = 0 and Δx = 0. We have verified the expressions for the arc centroid through comparison with a direct numerical calculation using discretized circles.