Differential Interferometric Signatures of Close Binaries of Supermassive Black Holes in Active Galactic Nuclei: II. Merged Broad Line Regions

Pairs of supermassive black holes (SMBHs) at different stages are natural results of galaxy mergers in the hierarchical framework of galaxy formation and evolution. However, identifications of close binaries of SMBHs (CB-SMBHs) with sub-parsec separations in observations are still elusive. Recently, unprecedented spatial resolutions achieved by GRAVITY/GRAVITY+ onboard The Very Large Telescope Interferometer through spectroastrometry (SA) provide new opportunities to resolve CB-SMBHs. Differential phase curves of CB-SMBHs with two independent broad-line regions (BLRs) are found to have distinguished characteristic structures from a single BLR \citep{songsheng2019}. Once the CB-SMBH evolves to the stage where BLRs merge to form a circumbinary BLR, it will hopefully be resolved by the pulsar timing array (PTA) in the near future as sources of nano-hertz gravitational waves. In this work, we use a parameterized model for circumbinary BLRs to calculate line profiles and differential phase curves for SA observations. We show that both profiles and phase curves exhibit asymmetries caused by the Doppler boosting effect of accretion disks around individual black holes, depending on the orbital parameters of the binary and geometries of the BLR. We also generate mock SA data using the model and then recover orbital parameters by fitting the mock data. Degeneracies between parameters contribute greatly to uncertainties of parameters but can be eased through joint analysis of multiple-epoch SA observations and reverberation mappings.


INTRODUCTION
The existence of supermassive black holes (SMBHs) at the center of nearly every galaxy has been confirmed by both star/gas dynamics of local galaxies (Genzel et al. 2010;Kormendy & Ho 2013) and reverberation mapping (RM) of distant active galaxies (e.g. Kaspi et al. 2000;Du et al. 2016). Galaxies grow through frequent mergers according to the hierarchical model of galaxy formation and evolution (White & Rees 1978;Lacey & Cole 1993;Patton et al. 2002;Lin et al. 2004;Conselice 2014), leading to the formation of binary SMBHs. Dual AGNs with separations ranging from several parsecs to several kiloparsecs have been identified through direct imaging (e.g. Rodriguez et al. 2006;Comerford et al. 2009). However, close binaries of SMBHs (CB-SMBHs) with sub-parsec separations have not been identified conclusively. The evolution mechanism of SMBH binaries after its separation decreases below a few parsecs is still debated (e.g. Milosavljević & Merritt 2003;Escala et al. 2005;Cuadra et al. 2009), and the existence and distribution of CB-SMBHs in the universe remain uncertain. Identification of CB-SMBHs would be a litmus test for the theory of galaxy mergers, resolving the long-standing "final parsec problem" completely. Additionally, CB-SMBHs are dominating sources of nano-hertz gravitational waves (GWs) background in our universe, which can be detected by pulsar timing arrays (PTAs) (Sazhin 1978;Jenet et al. 2006;Antoniadis et al. 2022). Sufficiently strong gravitational waves produced by massive CB-SMBHs can even be resolved individually (Sesana et al. 2009). If the GW source can also be identified by electromagnetic characteristics in advance, we can not only improve the sensitivity of the PTA through target search (Arzoumanian et al. 2020;Songsheng et al. 2021), but also restrict orbital parameters of the binary more precisely to test the physics of nano-hertz GWs through joint analysis of multimessenger signals.
Several electromagnetic characteristics have been proposed to search CB-SMBHs indirectly, among which the most widely used are periodic variations of the flux. In CB-SMBHs, periodic accretion of the gas in circumbinary disks or the Doppler boosting effect of the mini-disks around individual black holes (D'Orazio et al. 2015) can modulate the radiation flux periodically. Hundreds of candidates of CB-SMBHs with periodic variations have been selected from several photometric surveys, such as Catalina Real-time Transient Survey (Graham et al. 2015), Palomar Transient Factory (Charisi et al. 2016), Pan-STARRS1 Medium Deep Survey (Liu et al. 2019) and Zwicky Transient Facility (Chen et al. 2022). However, stochastic variabilities of quasars driven by thermal fluctuations (Kelly et al. 2009(Kelly et al. , 2011 or other random processes can mimic periodic signals in a finite time span (Vaughan et al. 2016;Goyal et al. 2018). Generally, the light curve should span at least three times the claimed period to avoid false periodicities . Other approaches mainly rely on the double-peaked or asymmetric profiles of broad emission lines (Shen & Loeb 2010;Popović 2012;Nguyen & Bogdanović 2016), which reflect the complex dynamics of the gas in the gravitational fields of the binary black holes. However, the broad-line region (BLR) around a single black hole can also be complicated enough to produce double-peaked or asymmetric profiles (Eracleous & Halpern 1994;Wang et al. 2017), and more information is needed to break the degeneracy. For example, long-term spectroscopic monitoring of NGC 4151 (Bon et al. 2012), NGC 5548 (Li et al. 2016), and Ark 120  show that their line profiles vary periodically, making them promising candidates for CB-SMBHs. Particularly, binary candidates with relatively extreme mass ratios and at separations between 0.1 and 1 pc can be selected from large populations of AGNs by detecting the time-dependent velocity offsets through multi-epoch spectroscopy (Kelley 2021). RM is also an efficient way to probe the geometry and velocity field of BLR through 2D transfer functions, which reflects the response lags of clouds with different line-of-sight (LOS) velocities to continuum variation (Blandford & McKee 1982). 2D transfer functions of BLRs in binary black holes are distinguished from those around single black holes since the broad emission line will respond to the continuum variation asymmetrically or out-of-sync at the red and blue wing Songsheng et al. 2020;Kovačević et al. 2020b;Ji et al. 2021).
The difficulty to identify CB-SMBHs lies in their extremely small angular sizes, which can not be resolved through direct imaging. Recently, the GRAVITY onboard The Very Large Telescope Interferometer (VLTI) has achieved an astrometric accuracy of ∼ 10 µas (Gravity Collaboration et al. 2017) and spatially resolved BLRs of several AGNs using spectroastrometry (SA) (Gravity Collaboration et al. 2018, 2020, 2021. By measuring the variation of the interferometric phase with wavelength across the broad emission line, SA detects the angular displacements of photocenters of BLR clouds with different LOS velocities, which further reflects the geometry and dynamics of the BLR. The differential phase curve of a single BLR around a SMBH usually shows a symmetric "S-like" shape with a peak and valley. When it is applied to CB-SMBHs with separated BLRs, the binary BLRs can be partially resolved through phase plateaus or extra peaks between the major peak and valley in interferometric phase curves (Songsheng et al. 2019;Kovačević et al. 2020a). For SMBHs with masses of 10 8 − 10 9 M , the size of their BLRs are usually larger than 0.1 pc (Bentz et al. 2013). If BLRs of the two black holes are still detached, frequencies of GWs emitted by them are too low to be detected using PTA. For multimessenger observation of CB-SMBHs, we must involve those with merged BLRs. On the one hand, CB-SMBH in a circumbinary disk can be diagnosed through their complex emissionline profiles generated by two mini-disks that are gravitationally bound to the individual black holes and the circumbinary disk (Nguyen & Bogdanović 2016;Nguyen et al. 2019Nguyen et al. , 2020. On the other hand, the orbital velocities of these CB-SMBHs usually reach a few percent of the speed of light, and the relativistic Doppler boosting effect of the mini-disks around individual black holes will be observably (Charisi et al. 2022). As a result, the circumbinary BLR will be illuminated anisotropically, and the distribution of emissivities of clouds will also be asymmetric. The RM of circumbinary BLR has been studied thoroughly in Ji et al. (2021), showing that the broad emission lines' blue and red wings are enhanced and weakened periodically but with a phase difference. Due to the coupling between timescales of light propagation and orbital motion, the RM cannot be fully described by a transfer function, making the identification not straightforward. Fortunately, the SA probes the dynamics of the BLR by taking a snapshot and so is free of the complications in the time domain, providing a different perspective to recognize CB-SMBHs.
In this paper, we established a parametric model for circumbinary BLRs and studied the dependence of differential phase curves on model parameters. We also fit the model to the simulated interferometric data and recover the orbital parameters of the CB-SMBHs. The paper is scheduled as the following. We present the fundamental formulations of the circumbinary BLR model in §2. The differential phase curves under different model parameters and the recovery of model parameters using mock data are provided in §3. Discussions on the issues of binary black holes are given in §4. We conclude in the last section. The geometry and kinematics of the BLR around a single SMBH have been thoroughly studied through RM in the last three decades. The velocity-resolved delay map constructed by maximum entropy method (MEM; Horne 1994;Xiao et al. 2018) or pixon-based method  can be used to identify the structure of the BLR, and demonstrate that flattened disk with Keplerian rotation is common in Seyfert galaxies Du et al. 2016). Multiple RMs of a few AGNs (such as in NGC 5548, 3C 390.3, and NGC 7469;Lu et al. 2016;Wandel et al. 1999) also show evidence for virialized BLRs. If the separation between the binary black hole is much larger than the sizes of both BLRs, the BLR of each black hole can be described by a detached flattened disk and the velocity of the cloud will be the superposition of the individual virial motion and orbital motion of the binary system (Shen & Loeb 2010;Popović 2012;Wang et al. 2018;Songsheng et al. 2019Songsheng et al. , 2020. As the binary of black holes evolves closer, two BLRs will merge to form a common circumbinary BLR and the dynamics of clouds will be much more complicated in the gravitational field of the binary. The geometry and kinematics of the BLR could be obtained by integrating the equation of motion of the circular restricted three-body problem until a quasi-equilibrium is established (Shen & Loeb 2010). However, the numeric integration is time-consuming and orbits of clouds may be chaotic and even become unbound in the long run.
The long-term stability of test particles in the gravitational field of a binary system has been examined both analytically (Szebehely 1980) and numerically (Holman & Wiegert 1999) for various mass ratios and binary eccentricities (e). For the binary with e = 0, test particles in circumbinary orbits will remain stable after 10 4 binary periods if the orbital radius is larger than 2.3 times the binary semimajor axis. So we assume the circumbinary BLR as a flattened disk with an inner radius larger than 2.3A, where A is the binary semimajor axis. The rotational velocity of the cloud is simply given by V = GM • /R, where G is the gravitational constant, R is the distance of the cloud to the center of mass of the binary, and M • is the total mass of the binary. To check the self-consistency of our simplified BLR model, we let the system evolves dynamically under Newtonian gravity (Murray & Dermott 1999). Initially, 2000 clouds distribute uniformly from R = 2A to R = 5A. After integrating the system dynamically for 1000 binary periods using Runge-Kutta-Fehlberg method, we obtain the result shown in Fig. 1. In Fig.  1(a), all clouds with initial distances larger than 2.3A are still bound after evolution but move a little bit outward collectively, while some clouds within 2.3A either become unbound or accreted by the black hole. In Fig. 1(b) and (c), rotational velocities of survival clouds still follow Keplerian law tightly and their radial velocities are small compared to the orbital velocity of the binary. All clouds are limited in the orbital plane of the binary and their velocities also have no component along the direction perpendicular to the orbital plane in the test, but the model is also valid as long as the opening angle of the BLR is small. We neglect the interaction between clouds, but our model is generally consistent with the model of a thin circumbinary gaseous disk governed by gas pressure and viscosity (Artymowicz & Lubow 1994), except that the mass flow through the gap is neglected here (Artymowicz & Lubow 1996).
Finally, we assume the radial emissivity distribution of clouds in BLR is generated by a shifted Γ-distribution (Pancoast et al. 2014), where R is the distance of the cloud to the center of mass of the binary, R BLR is the mean radius, F = (R in − 2.3A)/(R BLR − 2.3A) is the modified fraction of the inner to the mean radius, R in is the inner radius, β is the shape parameter, and Γ 0 = p(x|β −2 , 1) is a random number drawn from a Γ-distribution with α = β −2 and x 0 = 1. For a typical local massive CB-SMBH, we assume its angular size distance, total mass and separation are D A = 200 Mpc, M • = 10 9 M and A = 0.02 pc respectively. On the one hand, the angular separation of the binary is which can be resolved by GRAVITY/GRAVITY+ through SA. On the other hand, the frequency and strain amplitude of GWs emitted by the binary is and where z is the cosmological redshift of the binary, q ≡ M 1 /M 2 is the mass ratio between the primary (M 1 ) and secondary (M 2 ) black hole. So GWs emitted by the binary are expected to be detected by PTA in the era of SKA (Moore et al. 2015).  At the time tij, the accretion disk around the secondary black hole at rBH,2 emits an ultraviolet photon, whose trace is represented by the purple arrow . The photon then reaches the BLR cloud at rC,j and excites the hydrogen atom. The recombination line emission propagating along the line of sight (the green arrow) finally reaches the distant observer.

Doppler boosting effect
For our typical CB-SMBH, the ratio between the orbital velocity and speed of light is given by As the orbital velocity of the binary reaches a few percent of the light speed, the Doppler boosting will make the accretion disk a non-isotropic source of ionization radiation and the surface brightness of the BLR will also become asymmetric. We quantify this effect in a similar framework as that in Ji et al. (2021).
If the accretion disk of the i-th black hole emits an ultraviolet photon at time t ij , it will reach the j-th cloud at r C,j with velocity v C,j at time t ij + |r C,j − r BH,i |/c and excite the hydrogen atom there. The frequency of the photon will be boosted by where is the projection of the j-th cloud's relative velocity along the direction of its relative position to the i-th black hole and is the corresponding Lorentz factor. According to the conservation of I ν /ν 3 for the propagation of light in a vacuum, where ν is the frequency of the photon and I ν is the specific intensity at ν, the ionizing flux received by the j-th cloud is then given by where F ADi (ν) ∝ ν αν,i is the intrinsic SED of the accretion disk of the i-th black hole, and the term |r C,j | 2 /|r C,j − r BH,i | 2 accounts of the variation of the opening angle spanned by the j-th cloud due to the motion of the i-th black hole. Note that we assume the intrinsic variation of F ν,ADi is negligible compared to the variation of the Doppler boosting factor. The photon of the recombination line reaches the observer at a time where n obs is the unit vector pointing from the AGN to the distant observer, and T 0 is an overall delay for a photon propagating from the center of mass of the binary to the observer. The Doppler boosting factor due to the relative motion of the cloud relative to the observer is given by where β j = −v C,j · n obs /c and γ j = (1 − |v C,j | 2 /c 2 ) −1/2 . So the frequency of the photon received by the observer is ν c D j (1 + z) −1 , where ν c is the frequency of the emission line in the rest frame. The corresponding flux is Eq. (13) describes relative emissivities of clouds at different positions, which mainly depends on the Doppler factor D ij . In the direction where the accretion disk moves toward the cloud, the relative emissivity will boost significantly.

Spectroastrometry
Suppose the optical interferometer takes the observation at time T obs . For each black hole and cloud, we should solve the equation T ij (t ij ) = T obs for t ij . The equation can be further reorganized as The last term of the Eq. (14) is a small quantity and so the equation can be solved iteratively as Once all φ ij are obtained, the wavelength-dependent photocenter of the BLR can be figured out as where λ j = (1 + z)λ c D −1 j is the wavelength of the photon emitted by the j-th cloud, λ c is the wavelength of the emission line at rest frame; θ j = [r C,j − (r C,j · n obs )n obs ]/D A is the angular displacement of the cloud in the tangent plane of the sky.
The continuum emission in the K-band, which is covered by GRAVITY, mainly originates from the hot dust near the sublimation radius (Kishimoto et al. 2009), which is systematically larger than the size of the BLR by about a factor of four to five from infrared RM and interferometry (Koshida et al. 2014;Gravity Collaboration et al. 2020). Since the size of the dust torus is much larger than the separation between the binary, the photocenter of the continuum c will keep constant as the binary rotates. c is non-dynamic and independent of wavelength, only causing vertical shifts of the interferometric phase curve. To highlight the characteristics of the differential phase curve of the circumbinary BLR, we assume the location of the continuum photocenter to be the mass center of the binary, i.e. c = 0. The differential phase curve will be where f = F /(F + F c ) is the ratio of the emission line flux F to the total flux F + F c in each wavelength channel and B is the sky-projected baseline of the interferometer (Gravity Collaboration et al. 2018). All independent parameters used in the circumbinary BLR model are shown in Table 1. For convenience, we define a dimensionless BLR size r BLR ≡ R BLR /A in Table 1. Once the angular size distance of the AGN is fixed, other model parameters, such as binary separation, binary period, and total mass, can be calculated directly.

Dependence on model parameters
Firstly, we use the fiducial values in Table 1 as parameters of the merged BLR model to calculate the profile and differential phase curve. As shown in Fig. 3 (a), both the profile and the phase curve become asymmetric due to the Doppler boosting effect. Distributions of relative intensities and line-of-sight velocities of BLR clouds projected onto the celestial sphere are presented in Fig. 3 (b) and (c) respectively. The blue-shifted clouds are boosted by about 10 % while the red-shifted clouds are diminished, so the blue-shifted peak in the profile is higher than the red-shifted peak. The velocity gradient field of clouds is still symmetric, but the asymmetric f (λ) makes the blue-shifted peak in the phase curve high and narrow and the red-shifted valley shallow and wide. Now we increase the mean radius of the BLR from 4a to 8a and 16a, and the result is shown in Fig. 4. As the BLR gets larger, both the profile and phase curve become less asymmetric. When the size of the BLR is much larger than the separation between the binary, velocities of BLR clouds will be much smaller than that of the binary, and the equation 8 and 14 can be approximated To further simplify the equation, we consider only the clouds in the equatorial plane and suppose the inclination of the LOS is small. For the secondary black hole, we have  where ϕ C,j is the azimuthal angle of the cloud in the equatorial plane. For clouds having the same ϕ C,j + β a |r C,j |/A, their Doppler boosting factors will be the same, so a spiral structure will develop in the emissivity distribution of clouds, as shown in the second row of Fig. 4. We can also infer that the larger the rotational velocity of the binary is, the tighter the spiral arm is wound. As the spiral arm appears in the emissivity distribution, clouds with similar LOS velocities will have a wide range of emissivities, smearing the asymmetry in the profile and phase curve.
Next, we adjust the rotational velocity of the binary from 0.05 to 0.03 and 0.10, and the effect is presented in Fig. 5. A larger rotational velocity will lead to a more significant boosting effect, and so a more asymmetric profile and phase curve. However, we should note that if the rotational velocity is large enough to generate a spiral structure in the surface brightness distribution, the asymmetry will weaken again.
The effect of the orbital phase of the binary on the profile and phase curve is shown in Fig. 6. When the direction of the LOS velocity gradient is perpendicular to that of the emissivity gradient, as illustrated in Fig. 6(a), the profile and phase curve will be symmetric. As the binary rotates, the emissivity gradient will also change its direction accordingly, while the velocity gradient is always parallel to the long axis of the BLR's projection. The profile and phase curve will reach the maximum asymmetry when the gradient of emissivity is parallel with that of LOS velocity. So the profile and phase curve vary over the same period as the orbital period.
When the inclination angle increased, the projected velocity of clouds will get larger and so broaden the line profile and phase curve, as shown in Fig. 7. The doppler boosting effect mainly depends on the projected velocity of the binary "seen" by the BLR clouds, and so the asymmetry of the profile and phase curve is insensitive to the observer's inclination.
Finally, we consider the effect of the opening angle of the BLR. As the BLR gets thicker, the emissivities and LOS velocities of clouds will vary across the vertical direction of the BLR. Subsequently, the inhomogeneity of emissivities and velocities will be smoothed out when projected onto the celestial sphere, leading to a more symmetric line profile and phase curve. Meanwhile, double peaks in the line profile will also disappear for large opening angles, which is similar to the case of a single black hole. We should emphasize that our model only applies to thin BLRs. When clouds deviate significantly from the orbital plane of the binary, the stabilities of their orbits need to be further clarified, which may lead to a more complex BLR model.

Mock data analysis
The dependence of line profiles and differential phase curves on model parameters of circumbinary BLR indicates the possibility to infer the orbital parameters of the CB-SMBH through SA. To see that quantitatively, we use the fiducial values in Table  1 to generate mock data for SA observations 1 . The details of the mock data generation can be found in the Appendix. Given model parameters {Θ}, the likelihood function of data set D is where N λ and N b are numbers of wavelength channels and interferometer baselines respectively, F data,j and F model,j (Θ) are input and predicted value of the flux at the j-th wavelength channel, φ data,ij and φ model,ij (Θ) are input and predicted value of the differential phase at the j-th wavelength channel of the i-th baseline, and σ i and σ ij are corresponding uncertainties of the line profile and differential phase curve. The prior distribution of the model parameters P (Θ) are listed in Table 1. In light of Bayes's theorem, the posterior probability distribution for model parameters is which can be sampled numerically using diffusive nested sampling (DNest) algorithm (Brewer et al. 2011). The result is shown in Fig. 9. For most model parameters, the true values are within or deviate somewhat from the 1σ interval of the posterior probability distribution. The relative size of BLR r BLR and the angular size of the binary separation ξ a are highly correlated. Only their product, the angular size of the BLR, can be determined quite accurately from the amplitude of the differential phase curve. In our simulation, r BLR is overestimated by by a factor of 2 while ξ a is underestimated by about 50 %. The inclination of LOS  i 0 and opening angle of BLR θ opn are also correlated since both decreasing i 0 and increasing θ opn can make the system more symmetric and weaken the double-peak structure in the profile. The rotational velocity of the binary β a can be measured fairly accurately from the degree of asymmetry of the profile and phase curves. Its correlation with the inclination angle is weak since clouds in the BLR always view the accretion disks edge-on, independent of the LOS of the observer. However, the width of the profile is highly dependent on the inclination. The severe uncertainty of the inclination makes it arduous to determine the relative size of the BLR accurately.

Parameter space and uncertainty of orbital parameters
In this work, we mainly focus on the interferometric signatures of circumbinary BLRs around CB-SMBHs. Binary separations are required to be large enough for nonuniformity in BLR to be resolved spatially, but also small enough such that Doppler boosting caused by orbital velocities is detectable, which will be of 0.01 − 0.1pc generally. To detect the GWs from the binary through PTA in the near future, we further require massive enough (M • > 10 8.5 M ) binaries with non-extreme mass ratios. Circumbinary BLRs around CB-SMBHs may also be detected through line profiles (Nguyen & Bogdanović 2016) and RM (Ji et al. 2021). In both cases, we only require the orbital velocities to be large enough to be resolved spectroscopically or lead to a significant Dopper boosting effect, and so may explore a larger parameter space. For binary with separation 0.1pc, the BLR of each individual black hole will be detached. If the mass ratio is not extreme and each of the binary BLRs has similar radiation flux, we can identify them through RM ) and SA (Songsheng et al. 2019). While for binaries with extreme mass ratio, they can be selected from large populations of AGNs spectroscopically through time-dependent velocity offsets in line profiles (Kelley 2021).
The orbital parameters of the binary, such as binary separation and total mass, can be recovered from the posterior distribution of model parameters. The result is shown in Fig. 10. Input binary separation, orbital period, and total mass are all within the 2σ confidence interval of the posterior distribution, but they are underestimated by about 50% systematically.
The major origin of the bias is the degeneracy between the relative BLR size and binary separation. With SA observation, only their product, the absolute size of the BLR, can be determined accurately. Extra information or data set is needed to break their  degeneracy. RM provides a supplementary perspective to probe the structure of the BLR along the LOS. It has been demonstrated that joint analysis of RM and SA data can improve the measurement accuracy of black hole mass in the case of a single black hole . We would apply the joint analysis to check whether the uncertainty of the binary separation can be improved in a separate paper.
The minor origin of the bias is the large uncertainty of the inclination angle, which is always hard to measure accurately using line profiles. The uncertainty of the inclination also contributes to the scatter of virial factor in RM (Collin et al. 2006), hindering our accuracy in measuring masses of black holes through single-epoch spectrometry. For AGNs with large-scale radio jets, the inclination can be constrained by superluminal motion (e.g. Jorstad et al. 2017). However, the orbital motion of the binary may make jets from the two black holes wind with each other and form X-type lobes, leading to a more complicated scenario (e.g. Merritt & Ekers 2002;Cheung 2007). Another way to measure the inclination angle is through spectropolarimetry. Observations show that there are equatorial scattering regions in Type I AGNs (Smith et al. 2005). Unpolarized photons from the BLR are scattered in the direction of the LOS by electrons in the scattering regions and gain certain degrees of polarization. As a result, the width of the polarized line profile is determined by the LOS velocities viewed by equatorial scattering regions and independent of the inclination of the observer's LOS (Baldi et al. 2016;. By comparing the width of total and polarized line profiles, the inclination angle can be constrained properly. Another way to increase the measurement accuracy of orbital parameters is to use the periodicity of the binary. By monitoring the continuum variation modulated by the Doppler boosting effect, the orbital period can be determined accurately. Fixing the period when fitting the BLR model to the SA data, the uncertainties of orbital parameters will decrease. A more efficient method to constrain the orbital period is to conduct two SA observations several years apart and detect the difference between orbital phases. Moreover, if the binary is close or massive enough, the GW signals generated by them could stand above the stochastic background and be resolved individually by PTA (Sesana et al. 2009). In this circumstance, the data of timing residuals and SA can be analyzed jointly to achieve multimessenger observations. When the mean size of BLR is only a few times the separation between the binary, the region enhanced by the Doppler boosting is like a "hot spot" in the BLR. A similar feature has also been reported in the RM observation of Arp 151 and can be explained by a BLR with warped-disk geometry (Bentz et al. 2010). The warp structure can expose or shield gas according to its position, causing the required azimuthal structure-enhancing responses. It is hard to distinguish the circumbinary model from the warpeddisk model by direct imaging of the BLR, let alone by a single epoch SA observation. Besides detecting the periodicity of the light curve through the long-term monitoring campaign, it is possible to tell the two models apart using RM. In the circumbinary BLR model, the source of ionizing radiation is in orbital motion rather than at rest in the center. The light travel time from the inner accretion disk to the BLR varies with azimuthal angle and time, displaying special features in the responses of broad emission lines to the continuum (Ji et al. 2021).

Alternative model
If the size of BLR exceeds ∼ 10 times the separation between the binary, spiral structures develop in the distribution of emissivities of BLR clouds. Similar structures can also be generated by density waves of ionized gas in BLRs if the selfgravitating effects are important (Wang et al. 2022). In this case, not only the distribution of emissivities but also velocity gradients are modified by the "spiral arms", which is slightly different from the effect of Doppler boosting. It is possible to distinguish them if the kinematics of the BLR gas is well constrained using SA or RM observations with high fidelity.

Further improvement of the model
In our formulation of the Doppler boosting effect, the intrinsic radiation flux variation of the accretion disk is neglected for simplicity. Actually, the radiation from the i-th accretion disk F ADi in Eq. (13) depends on the departure time t ij of the photon ionizing the j-th clouds. For a rotating BLR, the difference of the departure time t ij for blueshifted and redshifted clouds is about A/c. Assuming the intrinsic variation can be described by a damped random walk, the characteristic amplitude of the variation on a timescale of A/c should be σ A/c, where σ 2 represents the variance in the light curve on short timescales (Kelly et al. 2009). If the mass of the black hole is about 10 9 M , σ 2 will be of the order 10 −4 mag 2 day −1 . Taking A = 10 −2 pc, we have ∆M int ∼ 5 × 10 −2 mag(∆F int /F ∼ 5 %), which is much smaller than the anisotropy caused by the Doppler boosting effect. However, the scatter of σ 2 is large among SMBHs with similar masses. If the effect of the intrinsic variation is comparable to that of the Doppler boosting, we should include it in the model, and fit the SA data and continuum light curve simultaneously.
Secondly, the geometry of the circumbinary BLR is simply a circular disk in our model. Nevertheless, the radiation of the accretion disk is highly anisotropic, so the ionization parameter varies with azimuthal angles. In the direction where the radiation is enhanced, the distances of line-emitting clouds will be large, and vice versa. As a result, the width of the line profile is also asymmetric in the red wing and blue wing, and the asymmetry of the differential phase curve is also enhanced. To address the effect, the geometry of the BLR can be parameterized by an egg shape (e.g. Narushin et al. 2021) rather than a circular.
Finally, we suppose that the dynamics of the BLR are stable and can be approximated by a Keplerian rotating disk. Small deviations from this can be treated by adding random velocities to clouds, parameterized by "turbulent" velocities σ turb . If the timescale for the BLR to reach the equilibrium state is longer than that of the orbital decay, we must include the evolution of the binary orbit in the dynamical simulation of cloud orbits to find a consistent velocity distribution for BLR clouds. Furthermore, clouds in the BLR may be continuously replenished by torus (Wang et al. 2017) or accretion disk (Goad et al. 2012), and so there could be clouds in dynamically unstable regions. A more physical model is needed to connect the torus, BLR, and accretion disk and give a complete description of the gas environment around the SMBHs. In this work, we establish a parameterized model for circumbinary BLRs, where the spatial and velocity distribution of clouds can be calculated conveniently and reliably. Due to the Doppler boosting effect of accretion disks around individual black holes and finite travel time of light, structures like "hot spot" or "spiral arm" will appear in the distribution of emissivities, leading to asymmetric line profiles and differential phase curves. Shapes of emission lines and phase curves vary with orbital parameters of the binary and geometries of the BLR. As a result, orbital parameters can be recovered from mock SA data by fitting the model to the data, demonstrating its ability to identify CB-SMBHs and constrain their orbits in the future. Uncertainties of parameters are mostly contributed by correlations between parameters. Through joint analysis of multiple SA observations, RMs and PTAs, correlations can be broken to realize robust measurements of CB-SMBHs.
We are grateful to the members of the IHEP AGN group for enlightening discussions. JMW thanks the support of the National Key R&D Program of China through grant -2016YFA0400701, by NSFC through grants NSFC-11991050, -11991054, -11833008, -11690024, and by grant No. QYZDJ-SSW-SLH007 and No.XDB23010400. The line profile and differential phase curves for mock data analysis are calculated using the fiducial values of model parameters listed in Table 1. The angular distance of the AGN is 200 Mpc, and the corresponding redshift is 0.05. So the emission line used for SA observation in the K band is Brγ. We normalize the profile of the Brγ line so that its flux relative to the continuum is about 0.1 at the peak. The profile is then convolved using a Gaussian with FWHM of 4 µm to account for the effect of instrument broadening, which corresponds to the middle spectral resolution mode of GRAVITY. The profile is sampled in 40 wavelength bins between 2.22 µm and 2.32 µm. At each bin, the uncertainty of flux measurement is 5 × 10 −3 relative to the continuum flux, which is also the typical value of line profiles obtained by GRAVITY. The result is shown in Fig. 11.
To get differential phase curves, we need to project the displacement of the photocenter onto baselines in the u − v plane, as indicated by Eq. 17. For simplicity, we assume there are 20 baselines with length of 100 m, and their azimuthal angles in the u − v plane vary from −90°to 81°uniformly. Wavelength bins of differential phase curves are the same as those of the line profile. At each bin, the uncertainty of interferometric phase measurement is 0.1°. The result is shown in Fig. 12.
Once the mock data is available, we can sample the posterior probability distribution of model parameters in the Bayesian framework by fitting the model to the mock data using DNest. We select 100 groups of parameters from the posterior sample randomly and calculate corresponding line profiles and differential phase curves, as shown by gray lines in Fig. 11 and 12. We also select the parameters with minimum chi-square from the posterior sample, and the best-fit curves predicted by them are shown as red lines in Fig. 11 and 12.  Figure 12. Mock differential phase curves and fitting results in each baseline. The blue dots with error bars are data points with 1σ uncertainties. The thick red line corresponds to the best-fit curve, while the thin gray lines are fittings using model parameters randomly drawn from the posterior sample.

B. SAMPLING THE POSTERIOR DISTRIBUTIONS OF MODEL PARAMETERS
The degeneracies between model parameters are the potential cause of discrepancies between the input and reconstructed values of model parameters when sampling the posterior distributions. As shown in Fig. 9, the i 0 − θ opn degeneracy and the r BLR − ξ a degeneracy are both significant. Firstly, we fix the inclination to the input value and then sample the probability distribution of the remaining parameters. The result is shown in Fig. 13. Just as we expect, the reconstructed value of the opening angle is exactly the same as the input one. However, although the reconstructed values of the BLR size and binary separation are closer to the input ones than those in Fig. 9, discrepancies between the input and reconstructed values are still larger than the 1σ uncertainties of the probability distributions. So the degeneracy between the inclination and opening angle is not the primary cause of the biases in sampling the posterior distributions.
Next, we fix the binary separation to the input value and conduct the sampling. The result is shown in Fig. 14. The input value of the BLR size coincides with the center of the posterior distribution, and the differences between reconstructed and input values for the inclination and the opening angle are also within the 1σ uncertainties of the probability distributions. We can conclude that the biases are mainly due to the degeneracy between the binary separation and the BLR size.
We further fix the period of the binary to the input value and rerun the program. The result is shown in Fig. 15. The i 0 − θ opn degeneracy and the r BLR − ξ a degeneracy still exist and biases are larger than the 1σ uncertainties of the posterior distributions. But the uncertainties of posterior distributions of the BLR size and binary separation are improved by a factor of 3 and so the estimation becomes more accurate. We now change the relative size of the BLR from the fiducial value of r BLR = 4 to 8, simulate the new line profile and differential phase curve, and rerun the sampling process. The result is shown in Fig. 16. The BLR size and binary separation are over-estimated and under-estimated by about 1σ respectively. But the uncertainties of their posterior distributions are smaller than those in Fig. 9. Meanwhile, the reconstructed values of the inclination and the opening angle also match the input values within 1σ uncertainties. So for larger BLR sizes, biases caused by the degeneracy between the binary separation and the BLR size are partially relieved, maybe because the amplitude of the phase curve gets larger.
Finally, we apply another sampling method, dynamical nested sampling method (Higson et al. 2019;Speagle 2020), to construct the posterior distributions of model parameters for comparison. The result is shown in Fig. 17. Similarly, the reconstructed values of most parameters differ from input values by more than one sigma. However, the posterior distributions of the BLR size, the inclination, and the binary separation show prominent multiple peaks, and the major peaks match the input values quite well. The result demonstrates the fitting of the model to the data has degenerate solutions, which cause biases when inferring the model parameters.