© The Authors 2024

1. Introduction

Supermassive binary black hole (BBH) systems are natural products of the hierarchical galaxy mergers in the Λ cold dark matter (ΛCDM) cosmological frame (e.g., Begelman et al. 1980; Yu 2002) since almost every massive galaxy hosts a supermassive black hole (BH) at its center (e.g., Kormendy & Richstone 1995; Magorrian et al. 1998; Kormendy & Ho 2013). Recently, Pulsar Timing Array (PTA) observations reported evidence of the stochastic nanohertz gravitational-wave (GW) background with detected Hellings–Downs correlation (Hellings & Downs 1983) at the ∼2 − 4σ confidence level (e.g., Agazie et al. 2023a; Antoniadis et al. 2023; Reardon et al. 2023; Xu et al. 2023), which is compatible with the predicted GWs produced by the cosmic BBHs (e.g., Agazie et al. 2023b; Chen et al. 2023; Ellis et al. 2024). Individual subparsec BBHs with periods of about a year are also the main targets of PTA experiments (e.g., Burke-Spolaor et al. 2019) and are expected to be detected soon (e.g., Chen et al. 2023). Electromagnetic observational evidence of such BBHs have also been searched for over more than a few decades. A large number of candidates were selected according to various observational signatures, but definitive evidence for subparsec BBHs is still elusive (see a recent review by D’Orazio & Charisi 2023).

Current BBH candidates at subparsec separations are mainly found by indirect signatures, such as double-peaked and/or asymmetric broad emission lines (BELs) (e.g., Eracleous & Halpern 1994; Gaskell 1996; Tsalmantza et al. 2011; Eracleous et al. 2012; Popović 2012; Decarli et al. 2013; Ju et al. 2013; Shen et al. 2013; Liu et al. 2014; Guo et al. 2019), periodic optical/UV light curves (e.g., D’Orazio et al. 2015; Graham et al. 2015; Charisi et al. 2016; Liu et al. 2019; Chen et al. 2024), or other methods (e.g., Yan et al. 2015; Zheng et al. 2016). However, each type of signature could have alternative model interpretations other than the BBH model, which means that only taking one of the above signatures makes it hard to confirm the existence of a BBH system. For example, the Keplerian rotation of a disk can also produce double-peaked BELs as observed (e.g., Chen & Halpern 1989; Eracleous et al. 1997). The periodic light curves can also be produced by damped random walk of the quasi-stellar object (QSO) flux variation among a large QSO sample (e.g., Vaughan et al. 2016; Liu et al. 2019). Direct resolving and identifying of these BBH candidates may require a spatial resolution of several microarcsec or smaller, which is not achievable by current facilities.

Current numbers of BBH candidates selected from various methods are on the order of 10², which are mainly obtained by the SDSS QSO spectroscopic survey (Tsalmantza et al. 2011; Eracleous et al. 2012; Liu et al. 2014), the Catalina Real-time Transient Survey (CRTS), Panoramic Survey Telescope and Rapid Response System (Pan-STARRS), Palomar Transient Factory (PTF), and Zwicky Transient Facility (ZTF) photometric surveys (Graham et al. 2015; Charisi et al. 2016; Liu et al. 2019; Chen et al. 2024). Confirmation of these BBH candidates may require further observations and detailed theoretical modeling, in which observations include continuous photometric monitoring to verify the periodicity (e.g., Liu et al. 2018) and spectroscopic observations to reveal the response of BELs (e.g., Li et al. 2016), and detailed theoretical modeling is necessary to establish combined binary black hole and broad-line region (BBH+BLR) systems to study how the BELs vary with the co-rotating BBHs (e.g., Ji et al. 2021a,b; Songsheng & Wang 2023). The combination of photometric–spectroscopic monitoring and theoretical modeling may improve the interpretation to these BBH candidates at subparsec separations (e.g., D’Orazio et al. 2015; Song et al. 2020, 2021).

When applying theoretical models to interpret observed light curves or BEL profiles of QSOs, a simple assumption was made frequently that only the secondary BH is active and the two BHs are co-rotating on a circular orbit (e.g., D’Orazio et al. 2015; Song et al. 2020, 2021; Ji et al. 2021a,b; Songsheng & Wang 2023). However, both observations and simulations indicate that BBH systems may have high orbital eccentricities (e.g., Valtonen et al. 2021; Jiang et al. 2022; Lai & Muñoz 2023). Hydrodynamic simulations also suggest that the accretion mode of BBHs depends on several parameters, including the mass ratio, eccentricity, disk thickness, and kinematic viscosity of the accretion disk(s) (e.g., Miranda et al. 2017; Duffell et al. 2020; D’Orazio & Duffell 2021; Dittmann & Ryan 2022; Siwek et al. 2023a,b). The accretion rate of the secondary BH (with mass M₂) tends to be substantially higher than that of the primary BH (with mass M₁) when q_M = M₂/M₁ is substantially less than 1, and it becomes equal to that of the primary BH when q_M ∼ 1 (e.g., Farris et al. 2014; Kelley et al. 2019; Duffell et al. 2020; Dittmann & Ryan 2023).

The ongoing and planned photometric and spectroscopic surveys will find a large number of QSOs with periodic light curves and will obtain their multi-epoch spectra. These surveys include the Vera Rubin Observatory’s Legacy Survey of Space and Time (Rubin, Ivezić et al. 2019), WFST (Wang et al. 2023), ZTF (Bellm et al. 2019), DESI (DESI Collaboration 2016), the Nancy Grace Roman Space Telescope (Spergel et al. 2015), SDSS-V (Kollmeier et al. 2017), and the Sitian project (Liu et al. 2021). Combining the observed periodic light curves and the variation of the BEL profiles together could provide coherent evidence for the existence of BBHs in these QSOs. Therefore, it is of great importance to explore the dependence of the variation pattern of the periodic light curves and the BEL profiles on the orbital parameters of BBH systems.

In this paper we focus on investigating the variation pattern of periodic light curves and BELs of active BBH systems with different accretion modes and orbital eccentricities associated with a circumbinary BLR (cBLR). This paper is organized as follows. In Sect. 2 we describe the basic procedures for establishing the BBH model system. In Sect. 3 we show the main results of the light curve and BEL variations of hypothesized BBH systems with different accretion modes and eccentricities. In Sect. 4 we discuss the uncertainties and applications of the current model. The conclusions are given in Sect. 5.

2. Models of BBHs with circumbinary BLR systems

For the BBH candidates selected from periodic light curves, the typical values of the BBH total mass and orbital period are M_•• ∼ 10⁹ M_⊙ and T_orb ∼ 2 yr, respectively (e.g., Graham et al. 2015; Charisi et al. 2016). In this work we hence take BBH systems with this typical mass and period for case studies by focusing on how the continuum flux and BEL profiles are dependent on mass ratio (q_M), BH accretion mode (denoted by the UV photoionization luminosity ratio of the two components q_L), orbital eccentricities (e), arguments of periapsis (ω), true anomaly (ϕ), and inclination angles (i_obs). Here we define the luminosity ratio as $q_L=L_{\text{bol}}^2/L_{\text{bol}}^1$, with $L_{\text{bol}}^2$ and $L_{\text{bol}}^1$ representing the bolometric luminosities from the accretion disks around the secondary and primary BHs, respectively. Here the bolometric luminosity of the ith (primary/secondary) BH $L_{\mathrm{bol}}^i=ϵ{\dot{M}}_i{c}^2$, where ϵ is the mass-to-energy conversion efficiency, Ṁ_i is the mass accretion rate of the ith BH, and c is the speed of light (e.g., Yu & Tremaine 2002). Figure 1 shows the sketch diagram of the BBH+cBLR system. Here we simplify the system by setting the longitude of ascending node Ω = 0° to clarify variations of continuum light curves and BELs that are caused by e and ω.

Fig. 1. Sketch of the configuration of a BBH+cBLR system, with the BBHs and BLR clouds rotating in the same direction. For the BBHs co-rotating in elliptical orbits, the longitude of ascending node is set to be Ω = 0° for simplicity, i.e., the sky plane (yellow) and the BBH orbital plane (gray) share the same x-axis in the barycentric coordinate system. The argument of periapsis is denoted as ω. The observer’s line of sight is set in the y − z plane with an inclination angle of iobs to the BBH orbital plane. The BLR has an opening angle of θo and its middle plane is the same as the BBH orbital plane.

The detailed model parameters are listed in Table 1. The mass ratios of the two BHs q_M are set to be either 0.1 or 0.5, which corresponds to a remnant by a minor or major merger, respectively. According to numerical simulations, the difference of mass accretion rates Ṁ₂ and Ṁ₁ can be simply represented as a function of q_M: Ṁ₂/Ṁ₁ = 1/(0.1 + 0.9)q_M) (Farris et al. 2014; Duffell et al. 2020). When q_M = 0.1, the value of Ṁ₂ is approximately five times higher than Ṁ₁. When q_M = 0.5, Ṁ₂ is approximately two times higher than Ṁ₁. In this work we simply assume two kinds of accretion modes for BBH systems. In the first, only the secondary BH is active (i.e., q_L = 1 : 0), which corresponds to those cases that the secondary BH dominates the flux variations in the Doppler boosting scenario. In the second mode both BHs are active with q_L = 1 : 1, which corresponds to those cases with a mass ratio of q_M ≲ 1 where the accretion rates of the two BHs are about the same. Considering that the real observations are more complex than simulations with a simplified parameter input, the analyses of the current case study includes the setup of (q_M, q_L) = (0.1, 1 : 0) and (0.5, 1 : 1), as suggested from simulations, and the cases of (q_M, q_L) = (0.1, 1 : 1) and (0.5, 1 : 0) for comparison.

We consider that the BELs (e.g., Hα and Hβ) are photoionized by UV photons emitted from the accretion disk. For the BBH+cBLR system, the BBHs and BEL profiles are correlated only if the accretion disk of the active BH could radiate UV photons. Here we simply assumed that the photoionization of BLR clouds by the UV photons emitted from the accretion disk of an active BH is effective only if the radius of the Roche lobe (R_crit) is larger than the typical radius of UV continuum emission (r_UV) (see Appendix A for a detailed analysis). Only with R_crit ≳ r_UV can we study the photoionization of BLR clouds by the central BHs. As shown in Fig. A.1, e ≲ 0.6 is a reasonable criterion for studying the Doppler boosting effect of BBH+cBLR systems with M_•• ∼ 10⁹ M_⊙ and T_orb ∼ 2 yr. In the case study we set e = 0 and 0.5 to compare the effects of circular and elliptical orbits of BBHs to the observed continuum and BEL flux variations.

We derived the semimajor axis a_BBH of BBH systems with given BBH total mass and orbital period as well as the orbital radii of the two BHs, the azimuthal and radial velocity components V_ϕ and V_r according to the elliptical motion of the two BHs (see Appendix B). According to the kinematic motion of the two BHs and the accretion mode settings, we can calculate the Doppler boosted continuum light curves and also the Doppler enhanced or weakened photoionization to BLR clouds. In the simplified BLR model (Appendix C), by assuming the BLR size of the BBH system R_BLR follows the same empirical relationship with optical luminosity as the single BH case, we find that R_BLR is roughly eight times larger than the separation of BBHs for those candidates selected by periodic light curves. This suggests that the BBH candidates have circumbinary BLRs. By setting a shifted Γ-distribution (e.g., Pancoast et al. 2014) of BLR geometry, we were able to calculate the response of BLR clouds to the continuum variation of the central ionizing source (Appendix D) by considering the following factors: 1) position variation of the two BH components, 2) time-dependent Doppler boosting or weakening effect of the ionizing flux to each BLR cloud, 3) the line emission from each BLR cloud with reverberation response to the central ionizing source, 4) the gravitational redshift of the photons emitted from each BLR cloud, and 5) the Doppler shifts of the photons emitting from each BLR cloud to the observer. By considering all these factors, we finally derived the coherent variation of BEL profiles with the continuum.

In the Doppler boosting scenario, the amplitudes of enhanced or weakened continuum and BEL fluxes depend on many parameters, such as those listed in Table 1. For an elliptical orbit (as shown in Fig. 1), by assuming the longitude of ascending node Ω = 0° for simplification, we sample the argument of periapsis ω ranging from 0° to 270° linearly with a step of 90° to show the dependence on the variation of the observed continuum light curves on it. We took six snapshots (i.e., ϕ = 0, π/3, 2π/3, π, 4π/3, and 5π/3) to monitor the phase evolution of BEL profiles modulated by the orbital motion of BBHs. The observer was set to have a viewing angle of either i_obs = 30° (close to face-on) or 60° (close to edge-on). The half opening angle of the BLR was set as 45° (flattened disk geometry; GRAVITY Collaboration 2018), with all the BLR clouds rotating in circular orbits. For the rotating direction of the two BHs and BLR clouds, all of them are assumed to be rotating counterclockwise.

3. Results

3.1. Variation of continuum light curves from BBHs in elliptical orbits

For a BBH system co-rotating in an elliptical orbit surrounded by a circumbinary BLR, the continuum light curve received by the observer depends not only on M_••, T_orb, q_M, and q_L, but also on e and ω, compared to the case of circular orbit.

3.1.1. q_L = 1 : 0

In the case of q_L = 1 : 0, Fig. 2 compares the continuum light curves obtained for cases with circular (e = 0, black dashed line) and elliptical orbits (e = 0.5 and ω = 0°–270°, colored from blue to red) with mass ratio q_M = 0.1 (left panel) and 0.5 (right panel). When only the secondary BH is active, its rotational velocity increases with decreasing q_M, which consequently leads to smaller variation amplitudes of the light curves for the system with q_M = 0.5 than those with q_M = 0.1.

Fig. 2. Comparison of continuum light curves emitted by BBH systems when only the secondary BH is active with mass ratio qM = 0.1 (left panel) and 0.5 (right panel) at iobs = 60°. The black dashed line shows the light curve from the system with circular orbit (e = 0), the other four lines colored from blue to red show those resulted from elliptical orbits (e = 0.5), with four different position angles of their major axes ω varying from 0° to 270°. The horizontal dotted line shows the light curve without the Doppler Boosting effect.

For a counterclockwise-rotating BBH+cBLR system, the evolution trend of the observed continuum light curve is mainly determined by the increase of decrease in the projected velocities to the observer. As shown in Fig. 1, for the cases of ω = 0° and 180°, the tangential lines at the pericenter are parallel to the observer’s plane (Y − Z plane for the current setup), which result in symmetric light curves at the first half and second half of each period, similar to that of the circular orbit but with different Doppler enhanced or weakened amplitudes and time durations. In addition to the counterclockwise rotation of the secondary BH, one can observe a continuum light curve with the strongest Doppler enhanced and the smallest Doppler weakened amplitudes at ω = 0°, and a light curve with the smallest Doppler enhanced and the largest Doppler weakened amplitudes at ω = 180°. For other cases (i.e., ω ≠ 0° or 180°) the observed continuum light curves are asymmetric. The shape of a continuum light curve is determined by the projected velocities of the secondary BH at the pericenter. When 0° < ω < 180°, the light curves show a variation pattern of a slow increase and rapid decrease (e.g., ω = 90° in Fig. 2). When 180° < ω < 360°, the light curves show a variation pattern of a rapid increase and slow decrease (e.g., ω = 270° in Fig. 2).

In the case of q_M = 0.5 (right panel of Fig. 2), the rotational velocities at each phase are systematically lower than that of q_M = 0.1 (left panel of Fig. 2), the corresponding Doppler enhanced or weakened amplitudes are hence also systematically smaller, but keeping the variation patterns the same as that of q_M = 0.1.

3.1.2. q_L = 1 : 1

When both the BHs are active, the amplitude of the continuum light curves depends on the luminosity ratio of the two BHs. The secondary BH has a higher rotational velocity than the primary BH, and hence contribute more to the Doppler enhanced or weakened amplitudes. As shown in the left column of Fig. 3, for q_L = 1 : 1 and q_M = 0.1, the emitted continuum light curve by the secondary BH is the same as that shown in Fig. 2, but with the assumed intrinsic flux decreased from 1 to 0.5. The amplitude of continuum light curve emitted by the primary BH is approximately ten times weaker than the secondary, and the shape of the light curve is opposite to that of the secondary BH due to the momentum conservation. With the mass ratio q_M increasing to 0.5 (right column of Fig. 3), the Doppler boosted amplitude of the secondary BH is systematically smaller than that of q_M = 0.1. On the other hand, the rotational velocity of the primary BH increases, which causes a higher amplitude of the light curve than the case of q_M = 0.1. The increasing mass ratio q_M indicates a decreasing (increasing) amplitude of light curves by the secondary (primary) BH, which means that at any time, the enhanced (weakened) flux modulated by the secondary BH always corresponds to the weakened (enhanced) flux by the primary BH. This means that the observed continuum light curves have smaller amplitudes than the case of q_L = 1 : 0. With increasing q_M, resolving the periodicity of these light curves hence requires a higher quality of photometric observations.

Fig. 3. Comparison of continuum light curves emitted by BBH systems when both the BHs are active, with qL = 1 : 1, iobs = 60°, qM = 0.1 (left column) and 0.5 (right column). The top panel shows the observed continuum light curve contributed by Doppler Boosting effect from the two co-rotating black holes. The Doppler boosted or weakened flux modulated by the secondary and the primary BHs are shown in the middle and bottom panels, respectively. The lines in different shapes and colors are the same as in Fig. 2.

With the interpretation on the behaviors of the light curve variations from BBH systems, we then investigate how the response of the circumbinary BLR varies for different cases.

3.2. The BEL profile variations for different BBH systems in circular orbits

To clarify how the elliptical orbits of BBHs can affect the response of BLR clouds, we first study the variation pattern of profile shapes and amplitudes with different mass and luminosity ratios of BBHs co-rotating in circular orbit. Figure 4 shows the profile variations in one period at six phases (ϕ = 0, π/3, 2π/3, π, 4π/3, and 5π/3) for the four models: (q_M, q_L) = (0.1, 1:0), (0.5, 1:0), (0.1, 1:1), and (0.5, 1:1). For each model, we show the flux variations relative to the mean profile obtained from the six phases in each corresponding bottom panel. In the case of BBHs co-rotating in the same direction with the BLR clouds, given an observer with i_obs > 0, the blueshifted clouds (moving toward the observer) contribute higher flux variations than the redshifted BLR clouds, due to the time-delay effect, and are characterized by the relative flux variations in each bottom panel of Fig. 4.

Fig. 4. BEL profile variations in a single orbital period for BBHs co-rotating in circular orbits with iobs = 60°. The two left columns show BBH systems when only the secondary BH is active, i.e., qL = 1 : 0, but with qM = 0.1 (left column) and 0.5 (middle-left column). The two right columns show profile variations when both the two BBHs are active with qL = 1 : 1, in the case of qM = 0.1 (middle-right column) and 0.5 (right column). In each column the top panel shows BEL profiles of the six phases (0 to 5π/3) in one period color-coded from blue to red (see legend), and the bottom panel presents the relative flux variation of BEL profiles compared to the mean profile obtained from that of the six phases.

When only the secondary BH is active (q_L = 1 : 0) the amplitude of the light curve for q_M = 0.5 is slightly smaller than that for q_M = 0.1 (Fig. 2). The BEL profile variations for q_M = 0.1 (left column of Fig. 4) and 0.5 (middle-left column of Fig. 4) are quite similar and the flux variation caused by the Doppler boosting effect for the model of q_M = 0.5 is also slightly smaller than that of q_M = 0.1.

The amplitudes of the BEL profile variations decrease significantly with increasing q_M for the case of q_L = 1 : 1. As shown in the two right columns of Fig. 4, the amplitude of flux variations for q_M = 0.5 is substantially smaller than that for q_M = 0.1. This dramatic decrease is not due to the Doppler boosting effect by the secondary BH, but is dominated by the ionizing photons emitted from the accretion disk of the primary BH. Although the profile variation caused by the secondary BH is slightly different between q_M = 0.1 and 0.5, the difference of rotating velocities and hence the Doppler boosting effect of their primary BHs for the two models increases dramatically from the cases with q_M = 0.1 to those with q_M = 0.5 (see the bottom panels of Fig. 3).

3.3. BEL profile variations for different BBH systems in elliptical orbits: dependence on ω and ϕ

When the two BHs are co-rotating in elliptical orbits, the two orbital parameters e and ω would make more complex profile variations compared to those in circular orbits. For a BBH+cBLR system with fixed M_•• and T_orb, the velocity difference between the apocenter and pericenter increases with increasing eccentricity. Given a viewing angle i_obs, the variation patterns of profile shapes and fluxes of observed BELs are mainly dominated by e and ω. Different values of ω (i.e., different projected velocities to the observer) and a different area of BLR clouds enhanced by the central BH, would make the observed BEL profiles have varied shapes at a certain phase. Instead, increasing e (i.e., rapidly increasing velocities at the pericenter; Fig. A.1) would increase the difference in profile shapes among different phases.

Figures 5 and 6 show the resulting BEL profiles from the Doppler boosting hypothesis of four BBH+BLR systems viewed at i_obs = 60° and 30°, respectively. In these two figures the columns from left to right show the cases with ω = 0°, 90°, 180°, and 270°. In each panel phase 0 means that the counterclockwise-rotated secondary BH is located at the pericenter, and phase π corresponds to that rotated to the apocenter.

Fig. 5. Comparison of continuum light curves emitted by BBH systems when only the secondary BH is active (top two rows) and both the two BH are active (bottom two rows) for e = 0.5 and iobs = 60°. In each row, the top panels show the profiles at six phases (0 to 5π/3) in one period, and the relative flux variations are shown in the corresponding bottom panels as that in Fig. 4. Columns from left to right show the cases with ω = 0° ,90° ,180°, and 270°, which reflects the profile variations caused by different eccentricity vectors of the elliptical orbits observed at a fixed viewing angle.

Fig. 6. Legends are the same as those for Fig. 5, except that iobs = 30°.

In Fig. 5, when the BBH+BLR system is viewed close to edge-on (i_obs = 60°), double-peaked profiles appear for different ω and phases. As analyzed in Sects. 3.1 and 3.2, the system with (q_M, q_L) = (0.1, 1:0) has the highest light curve variation amplitude with respect to the other three systems and the greatest BEL profile variations. The amplitudes of BEL profile variations in one period for elliptic BBH systems are different from those circular ones (Fig. 4). In the case of e = 0.5, the BBHs at phases from ϕ = π/2 to 3π/2 that cross the pericenter have shorter rotating timescale and higher velocities than those from 3π/2 to π/2. For the case when all the BLR clouds are rotating in circular orbits, given the fixed inclination angle of an observer, the BEL profile at a certain phase would not vary. However, for BBHs in elliptical orbits the varying ω can cause profile variation at a certain phase.

With the argument of periapsis ω varying from 0° to 360° (left to right columns of Fig. 5), BEL profiles have large variations at a given phase as the positions of each of the two BH components relative to the observer at that phase are different for cases with different ω. Co-adding with the time-delay effect, the envelope of BEL profile variations in each panel are hence different for those cases with different ω. For example, at ϕ = 2π/3 (green line) the double-peaked profile has similar fluxes at the blue and red peaks at ω = 0°. When ω changes to 90°, the blue peak has higher flux than the red peak. At ω = 180°, the total flux of the profile increases, and the red and blue peak fluxes become similar again. While for ω = 270°, the flux of the red peak is significantly higher than that of the blue peak. Therefore, in the elliptic BBH system, two periodic signatures appear. The first is a periodic profile variation in one period (e.g., ϕ = 0 to 5π/3 shown in each panel), which is dominated by rotation of the central BHs. The second signature is the periodic profile variation in a certain phase (e.g., ω = 0° to ∼270° at ϕ = 2π/3), which is caused by ω. In the case when the BBH systems have no orbital precession or the timescale of the precession is hugely longer than the rotating timescale, it is only possible to see the periodic BEL profile variation in one period.

When observing the BBH+cBLR system at a viewing angle close to face-on (e.g., i_obs = 30° in Fig. 6), the BEL profiles only present Gaussian or asymmetric shapes, and the amplitudes of flux variation decrease significantly for all four model systems compared to that observed for i_obs = 60°. However, the relative flux variation trends are still similar to that shown in Fig. 5.

3.4. Varied amplitudes of BELs for different continuum light curves and orbital eccentricities of BBHs

With the interpretation on the evolution trend of BEL profiles for different ω, we explore how the amplitude of BEL profiles (i.e., $A_{\mathrm{BEL}}=1.25(log{L}_{\mathrm{BEL}}^{\max }-log{L}_{\mathrm{BEL}}^{\min })$) vary with the intrinsic continuum emission from the central ionization sources in different amplitudes (i.e., $A_{\mathrm{conti}}=1.25(log{L}_{\mathrm{conti}}^{\max }-log{L}_{\mathrm{conti}}^{\min }))$, mass ratio (q_M), and orbital eccentricity (e).

When only the secondary BH is active, as shown in the two left panels of Fig. 7, at fixed e, higher q_M means lower velocity of the secondary BH, and hence lower A_BEL and A_conti. For the case of M_•• = 10⁹ M_⊙, when the eccentricity e increases from 0 to 0.6, the rotation velocity at the pericenter for the secondary BH increases, which corresponds to both increased A_conti and A_BEL. For A_BEL, its increasing trend slows down with increasing e. The variation of A_BEL is caused by two competing factors: (1) the increasing rotation velocity caused by increasing e at phases ϕ = 3π/2 to π/2 (at the pericenter ϕ = 0) could contribute stronger Doppler enhancement to the ionization of BLR clouds; (2) the decreasing time duration at phases ϕ = 3π/2 to π/2 would decrease the fraction of BLR clouds with stronger Doppler enhancement.

Fig. 7. Variation of BEL amplitudes (ABEL) with different continuum amplitudes (Aconti) and orbital eccentricities (e) of BBHs with M•• = 109 M⊙ at ω = 0°. The two left columns show ABEL as a function of Aconti (left panel) and e (middle-left panel) for the case when only the secondary BH is active (qL = 1 : 0), and the two right columns show ABEL vs. Aconti (middle-right panel) and ABEL vs. e (right panel) for the case when both the BHs are active with qL = 1 : 1. In each panel, the points connected by lines of the same color show results of a fixed mass ratio, and the four different mass ratios qM = 0.1, 0.3, 0.5, and 0.7 are colored in blue, green, violet, and red, respectively. Each group of four points with different qM but the same e are linked by gray lines.

For the case when both the BHs are active with q_L = 1 : 1 (two right panels of Fig. 7), the response of BLR clouds to the central ionization source is affected by both the secondary and primary BHs. The ionization behavior of the secondary BH is the same as that of q_L = 1 : 0, while the behavior of the primary BHs is opposite to that of the secondary BH. The differences between flux amplitudes and variation trends shown in the two left and two right panels are caused by the activity of the primary BH. The intensity ratio of the two beams depends on the luminosity ratio q_L and on the mass ratio q_M. The higher luminosity fraction taken by the primary BH or higher mass ratio of the two BHs indicate enhanced ionization of BLR clouds by the primary BH, which appear in the Doppler weakened regions that are modulated by the secondary BH. Therefore, the total fluxes of BELs in the case of q_L = 1 : 1 increase and A_BEL decreases compared to the case of q_L = 1 : 0. This explains why A_BEL of q_L = 1 : 1 is systematically smaller than that of q_L = 1 : 0.

As shown in the middle-right panel of Fig. 7, the variation of A_conti with e is different for increasing q_M. At q_M = 0.1, A_conti first decreases and then increases with increasing e. At q_M ≳ 0.3, A_conti decreases monotonically with e at higher q_M. The value of A_conti is determined by the differential amplitude between the primary and secondary BHs (Fig. 3) as a function of e since the velocities of the two BHs increase with increasing e (see $V_{\mathrm{peri}}^{\sec }$ vs. e in the right panel of Fig. A.1). The varying trends of A_conti versus e at different q_M reflect the competitive contribution of continuum emission by the two BHs with the Doppler boosting factor D⁵ (Eq. (D.4)).

Although the activity of the primary BH decreases both the observed A_conti and A_BEL (the left and middle-right panels of Fig. 7), the BEL profiles still have an increasing A_BEL with increasing e (the right panel of Fig. 7) since the response of BLRs to the central ionizing sources is still dominated by the secondary BH.

For q_L varying from 1 : 0 to 0 : 1, the increasing luminosity fraction of the primary BH corresponds to decreasing A_conti and A_BEL when the secondary BH dominates the Doppler boosting effect. Then the variation trends and phases of the continuum light curve and BEL profiles become inverse when the primary BH dominates the Doppler Boosting effect.

When changing the viewing angle i_obs from 0° (face-on) to 90° (edge-on), a similar variation pattern of A_conti and A_BEL with e appear at different i_obs. With increasing e, the value of A_conti and A_BEL increase. On the other hand, the values of A_conti and A_BEL also increase with increasing i_obs for both q_L = 1 : 0 and 1 : 1 cases (Fig. 8), which means that the periodic variations of continuum and BEL fluxes are more detectable with larger inclination angles.

Fig. 8. Variation of continuum (Aconti) and BEL amplitudes (ABEL) with different inclination angles for BBHs with M•• = 109 M⊙, ω = 0°, and e = 0.0–0.6. The two left panels show the case when only the secondary BH is active (qL = 1 : 0), and the two right panels show the case when both the BHs are active with qL = 1 : 1. In each panel, the points colored from dark blue to red represent the orbital eccentricity of BBHs e = 0.0–0.6, respectively.

In this paper, as a case study, we focus on the continuum and BEL profile variations of the BBH+cBLR system by assuming Ω = 0° (i.e., cosΩ = 1, for simplification). For those observed periodic QSOs with arbitrary Ω ≠ 0°, if all other parameters are fixed, increasing Ω means decreased amplitudes of both A_conti and A_BEL. However, their coherent variation behaviors remain the same as those analyzed above.

4. Discussion

We investigated the coherent variations of light curves and BEL profiles for co-rotating BBH systems with different orbital eccentricities, by assuming either that only the secondary is active or that both the BHs are active with equal luminosity. However, the observed cases should be more complex. In the following we discuss the complexity in the BBH+cBLR modeling and the observational strategy on how to identify BBH systems.

4.1. Different BLR geometries

To explain the observed BEL profiles that vary from single-peaked to asymmetric or double-peaked shapes, the inclination angle could be changed from being viewed face-on to edge-on, and the BLR may also be changed from elliptical (e.g., Shen & Loeb 2010) to disk-like (e.g., Eracleous & Halpern 1994) geometries. For the observed line asymmetry, the disk-like BLRs with clouds in elliptical orbits are preferred for modeling (e.g., Kovačević et al. 2020). Some complex BLR structures have also been investigated, for example BLRs in thin and thick disks with inflow or outflow substructures (e.g., Wang et al. 2018) or with two flattened disks perpendicular to each other (e.g., Ji et al. 2021b). If considering co-planar BBH+cBLR systems with the two BHs co-rotating in circular orbit, the BLRs varying from spherical to thin disk indicate an increasingly sensitive response of BLR clouds to the central source (e.g., Ji et al. 2021a).

4.2. Orbital decay of BBHs

The orbital evolution of two BHs at different semimajor axes is controlled by different decay mechanisms, for example from disk-dominated viscous evolution at large separations to the secondary-dominated viscous evolution stage (the distance of BBHs at the pericenter D_peri ≲ 10⁴R_Sch for q_M ∼ 0.1), then the gravitational wave-dominated (GW-dominated) evolution stage (D_peri ≲ 500R_Sch), and finally the gas-disk decoupled stage at D_peri ≲ 100R_Sch, where R_Sch is the Schwarzschild radius corresponding to the total mass of BBHs (e.g., Haiman et al. 2009).

For BBHs dominated by the GW-driven orbital decay mechanism, the GW decay timescale is

$$\begin{array}{c}\hfill t_{\mathrm{GW}}=(8.2\times {10}^3\mathrm{yr}){\left[\frac{q}{{(1+q)}^2}\right]}^{-1}{\left(\frac{M_{\bullet \bullet }}{{10}^9{M}_{\odot }}\right)}^{-5/3}{\left(\frac{T_{\mathrm{orb}}}{2\mathrm{yr}}\right)}^{8/3}F(e),\end{array}$$(1)

where F(e) is the enhancement factor caused by orbital eccentricity e:

$$\begin{array}{c}\hfill F(e)=\frac{1+\frac{73}{24}{e}^2+\frac{37}{96}{e}^4}{{(1-e^2)}^{7/2}}.\end{array}$$(2)

From Eq. (1) we can derive that at M_•• = 10⁹ and T_orb = 2 yr the GW-driven orbital decay timescales are ∼1.2 × 10⁵ and ∼4.5 × 10⁴ yr for q_M = 0.1 and 0.5 at e = 0.5, respectively. The orbital decay does not affect the periodicity of BBHs in this work significantly.

4.3. An optimized way of identifying BBH candidates

Currently, there are two efficient ways of searching for BBH candidates: by selecting periodic light curves (e.g., D’Orazio et al. 2015; Graham et al. 2015; Charisi et al. 2016; Liu et al. 2019; Chen et al. 2024) and based on asymmetric or double-peaked BEL profiles (e.g., Eracleous & Halpern 1994; Gaskell 1996; Tsalmantza et al. 2011; Eracleous et al. 2012; Popović 2012; Decarli et al. 2013; Ju et al. 2013; Shen et al. 2013; Liu et al. 2014; Guo et al. 2019). However, both the dynamical modeling of BBHs (e.g., D’Orazio et al. 2015; Jiang et al. 2022) and BEL profile modeling (e.g., Shen & Loeb 2010; Nguyen et al. 2019, 2020) are difficult to distinguish from the single-BH case, due to a set of model degeneracies (see reviews by Wang & Li 2020; D’Orazio & Charisi 2023). Limited by photometric and spectroscopic observations, current efforts on the selection and identification of BBH candidates have mainly focused on the continuum light curves or BEL profiles separately. In Ji et al. (2021a,b), we analyze the detailed response of BLR clouds and BEL profiles to the central BBHs that are co-rotating in circular orbits by focusing on the cases where the BBH orbital plane is co-planar (Ji et al. 2021a) or misaligned (Ji et al. 2021b) with the middle plane of the circumbinary BLR. In this paper we investigate the coherent variation of BBH continuum light curves and BEL profiles for BBHs orbiting at different eccentricities with simplified orbital orientation, which indicates that the joint analyses of periodic light curves and multi-epoch observed BELs could lead to successful identification of BBH candidates.

As analyzed in this work, the amplitudes of periodic continuum light curves and BELs are correlated with each other (see Figs. 7 and 8). Because of the momentum conservation, the active primary BH could cause non-sinosoidal continuum light curves and decreased A_conti (Figs. 2 and 3). The BEL profiles at different phases in one period not only have significant differences in shapes and fluxes (Figs. 4–6) caused by orbital eccentricity and argument of periapsis of BBHs at certain (q_M, q_L), but also have a decrease in flux variation amplitudes with increasing mass ratio and luminosity contribution by the primary BH (Figs. 7 and 8). When taking into account the inclination angle, the amplitudes of continuum light curves and BEL profiles both increase with increasing inclination angles. This means that modeling the continuum light curves and multi-epoch observed BEL profiles together would improve the probability of identifying BBH candidates, and the BBH scenario can be confirmed if a clear coherent variation of the BELs with the continuum can be found.

On the other hand, the ongoing (e.g., WFST, Wang et al. 2023) and future photometric surveys (e.g., the Rubin and Sitian project, Ivezić et al. 2019; Liu et al. 2021) would find more BBH+cBLR systems occupied at much wider parameter space compared with those currently obtained by CRTS, Pan-STARRS, PTF, and ZTF projects (e.g., Graham et al. 2015; Charisi et al. 2016; Liu et al. 2019; Chen et al. 2024). There will be more spectroscopic data obtained from DESI (DESI Collaboration 2016), the Nancy Grace Roman Space Telescope (Spergel et al. 2015), and SDSS-V (Kollmeier et al. 2017) surveys. This makes the simultaneous modeling of periodic light curves and BEL profiles being possible for many BBH candidates.

5. Conclusion

In this paper we investigated the response of the line emission from BLR clouds to the continuum emission from central BBH systems co-rotating in elliptical orbits by focusing on the evolution trends of both observed light curves and BEL profiles. For QSOs with the central BBHs surrounded by a circumbinary BLR (BBH+cBLR) system, we investigated two cases of BH activity, one where only the secondary BH is active (q_L = 1 : 0), and the second where two BHs are active with equal luminosity (q_L = 1 : 1), with simplified orbital orientation (i.e., Ω = 0°).

For BBHs co-rotating in elliptical instead of circular orbits, the Doppler boosting effect caused by their orbital eccentricity and argument of periapsis could affect both the shape of continuum light curve, for example sharply (slowly) increasing but slowly (sharply) decreasing trends in a period, and two kinds of periodicities in BEL profile variations in the Doppler boosting hypothesis. On the other hand, BBHs in different active behaviors can result in quite different amplitudes of light curves and BEL profiles. At q_L = 1 : 0 the variation amplitudes of continuum light curves and BELs both increase with increasing eccentricity, but decrease with increasing mass ratio of BBHs. Instead, for q_L = 1 : 1, the activity of the primary BH can cause systematically decreased amplitudes of continuum light curves and BELs, due to its enhanced ionization to cBLR clouds in the Doppler weakened area that are modulated by the secondary BH. Currently, we only consider BHs that accrete with stable accretion rate; in the future, we will also explore how BBH systems vary with different fluctuations in accretion rates, such as the damped random walk variation.

The confirmation of BBH candidates becomes more complex when considering the elliptically co-rotated orbits and accretion activities of BBHs rather than those cases with simply assumed circular orbits or single BH activity. With the theoretical understanding on BBH systems at varied mass ratios, luminosity ratios, orbital eccentricities, and BH accretion activities, the incoming huge amount of photometric observations (e.g., by Rubin, Sitian, and WFST) and multi-spectroscopic observations (e.g., by DESI, Roman space telescope, and SDSS-V) would search for BBH candidates with more complete parameter spaces than mentioned in this work, based on which we could constrain the merger rate of massive black hole binaries.

Acknowledgments

We would like to thank the anonymous referee for the suggestions that helped us to improve this paper. This work is supported by the National Key program for Science and Technology Research and Development (grant Nos. 2020YFC2201400 and 2022YFC2205201), the Beijing Municipal Natural Science Foundation (No. 1242032), the National Natural Science Foundation of China (NSFC) (grant Nos. 11903046, 12273050, 11991052, 11988101, and 11933004), and the Strategic Priority Program of the Chinese Academy of Sciences (Grant No. XDB23040100). JG acknowledges support from the Youth Innovation Promotion Association of the Chinese Academy of Sciences (No. 2022056) and the science research grants from the China Manned Space Project, and JFL acknowledges support from the New Cornerstone Science Foundation through the New Cornerstone Investigator Program and the XPLORER PRIZE.

References

Appendix A: Stability of the accretion disks for BBH systems

For the active BH in the BBH system, we simply assume that the BH accretion disk is steady inside the Roche lobe. When studying the photoionization process of UV photons that radiated from the accretion disk to BLR clouds, the minimum radius of the accretion disk should make sure the quantity of radiated UV photons is enough for ionizing BLR clouds and result in matchable continuum and BEL fluxes, as observed from those single BH cases; in other words, the radius of Roche lobe R_crit is equal to or larger than the half light radius of the accretion disk at UV band (r_UV). Only with R_crit ≳ r_UV could the BEL profiles be observed and the dynamics of the BBH system encoded. For BBHs systems (Eggleton 1983), the radius of the Roche lobe for the BH could be appoximated within 1% of its value with the formula

$$\begin{array}{c}\hfill R_{\mathrm{crit}}=\frac{0.49q_M^{2/3}}{0.6q_M^{2/3}+ln(1+q_M^{1/3})}.\end{array}$$(A.1)

At the pericenter, the minimum radius of the Roche lobe for the secondary BH is $R_{\mathrm{crit}}^{\min }=R_{\mathrm{crit}}(1-e)$. If the radius of the accretion disk emitting UV photons is smaller than the radius of the Roche lobe, the BBH+cBLR system would be reasonably established for reflecting the BBH features by reverberation mapping procedures. For a BH with M_• = 10^8.6 M_⊙, the typical radius of the UV continuum emission at 1350Å in the accretion disk is $r_{\mathrm{UV}}^{1350\mathrm{\AA }}=21.4R_{\mathrm{Sch}}$ for the QSO SDSS J1339+1310 (Shalyapin et al. 2021). If estimating r_UV at 2500 Å with the empirical correlation derived from 14 QSOs (Morgan et al. 2018), the BH masses ranging from 10⁸ to 10⁹ M_⊙ correspond to $r_{\mathrm{UV}}^{2500\mathrm{\AA }}\sim 20$ to 30 R_Sch. Given a fixed BH mass and accretion rate, $r_{\mathrm{UV}}^{1350\mathrm{\AA }}$ is smaller than $r_{\mathrm{UV}}^{2500\mathrm{\AA }}$. Figure A.1 shows R_crit as a function of e, from which we can see that when studying a BBH system with M_•• = 10⁹ M_⊙ and T_orb = 2 yr, focusing on those BBHs with orbital eccentricities e ≲ 0.6 would be reasonable (left and middle panels) if using $r_{\mathrm{UV}}^{1350\mathrm{\AA }}$ of SDSS J1339+1310 as a criteria, in which the rotational velocity of the secondary BH ranging from ∼10⁴ km/s to ∼3 × 10⁴ km/s (right panel) could result in a significant Doppler boosting effect to the photoionization of BLR clouds and to the flux variation received by observers.

Fig. A.1. Radii of Roche lobes of BBHs and rotating velocity of the secondary BH at the orbital pericenter as a function of orbital eccentricity of the two BHs. The left and middle panels show the radii of Roche lobes ($R_{\mathrm{crit}}^{\mathrm{pri}}$ and $R_{\mathrm{crit}}^{\sec }$) in units of the Schwarzschild radii for the primary ($R_{\mathrm{Sch}}^{\mathrm{pri}}$) and secondary ($R_{\mathrm{Sch}}^{\sec }$) BHs, respectively. The typical radius of the UV emission (rUV) in the accretion disk for a M• = 108.6 M⊙ is shown as a horizontal dashed line for comparison. The right panel presents the rotating velocity of the secondary BH at the pericenter ($V_{\mathrm{peri}}^{\sec }$), which is the maximum velocity in each period. In each panel are shown the results of M•• = 108 (blue), 108.5 (green), and 109 M⊙ (red) with Torb = 2 yr and qM = 0.5. Considering that the UV emission might disappear when Rcrit ≲ rUV, the current model configuration only focuses on those systems with Rcrit ≳ rUV, i.e., e ≲ 0.6 for M•• = 109 M⊙ as a case study.

Appendix B: Dynamics of the BBH system

Assuming a BBH system with an orbital period of T_orb = 2 yr, which is roughly the median intrinsic period of BBH candidates in the sample of Charisi et al. (2016), the semimajor axis of the two black holes with total mass M_•• = 10⁹ M_⊙ is

$$\begin{array}{c}\hfill a_{\mathrm{BBH}}=7.7\times {10}^{-3}{\left(\frac{M_{\bullet \bullet }}{{10}^9{M}_{\odot }}\right)}^{1/3}{\left(\frac{T_{\mathrm{orb}}}{2\mathrm{yr}}\right)}^{2/3}\text{pc},\end{array}$$(B.1)

In the case of elliptical orbits (e.g., Valtonen et al. 2008; Charisi et al. 2022), the velocity of the secondary BH depends on its location at different position angle. Given the secondary BH rotating in elliptical orbit with eccentricity e and mass ratio q_M in the X-Y plane (Figure 1), the radius of the secondary BH is

$$\begin{array}{c}\hfill r_{\sec }=\frac{a_{\mathrm{BBH}}(1-e^2)}{(1+q_M)(1+ecos\varphi )},\end{array}$$(B.2)

and the velocity components V_ϕ and V_r of the secondary BH are

$$\begin{array}{c}\hfill V_{\varphi }=r_{\sec }\frac{d\varphi }{\mathit{dt}}=\sqrt{\frac{G{M}_{\bullet \bullet }}{a_{\mathrm{BBH}}(1-e^2)}}\frac{1+ecos\varphi }{1+q_M},\end{array}$$(B.3)

$$\begin{array}{c}\hfill V_r=\frac{d{r}_{\sec }}{\mathit{dt}}=\sqrt{\frac{G{M}_{\bullet \bullet }}{a_{\mathrm{BBH}}(1-e^2)}}\frac{esin\varphi }{1+q_M},\end{array}$$(B.4)

where G is the acceleration of gravity and ϕ is the true anomaly of the secondary BH. The position and velocity of the primary BH can then be derived by momentum conservation.

Appendix C: The simplified BLR model

For a BBH system with a cBLR (BBH+cBLR), we simply assume its BLR size has the same relation with the luminosity of the central ionization source as that of single-BH systems (Kaspi et al. 2000, 2007; McLure & Jarvis 2002; Peterson et al. 2004), although it could be different. The BLR size of the BBH system can then be derived by following the empirical relationship between the BLR size and optical luminosity (e.g., Lu et al. 2016). For convenience of model construction, we convert the optical luminosity to the Eddington ratio (λ_Edd) and M_•• by assuming λ_Edd = 0.1 for the BBHs:

$$\begin{array}{c}\hfill R_{\mathrm{BLR}}=0.127\times {({\lambda }_{\mathrm{Edd}}/0.1)}^{0.54}{(M_{\bullet \bullet }/{10}^9{M}_{\odot })}^{0.54}\mathrm{pc}.\end{array}$$(C.1)

The typical BLR size is approximately eight times larger than the separation of the BBHs, significantly far away from the central two BHs, and can be taken as a cBLR affected little by the central BBH dynamics. Therefore, in this work, we consider the configuration of BBHs surrounded by a cBLR, for which we take a shifted Γ-distribution (see details in Pancoast et al. 2014) as the radial emissivity distribution of BLR clouds:

$$\begin{array}{c}\hfill R_{\mathrm{ga}}=R_{\mathrm{S}}+F{R}_{\mathrm{BLR}}+g(1-F){\beta }^2{R}_{\mathrm{BLR}}.\end{array}$$(C.2)

Here R_S is the Schwarzschild radius, R_BLR is the mean BLR radius, F = R_min/R_BLR is the fractional inner BLR radius (R_min), β is the shape parameter, and g = p(x|1/β², 1) is drawn randomly from a Gamma distribution.

As to the angular distribution of the BLR clouds, for a BLR with an half opening angle θ_o (Figure 1), we assume

$$\begin{array}{c}\hfill \theta =arccos[cos{\theta }_{\mathrm{o}}+(1-cos{\theta }_{\mathrm{o}})U],\end{array}$$(C.3)

where θ represents the angle between the orbital angular momentum of a cloud and the normal to the BLR middle plane, and U is a random number ranging from 0 to 1 and is set to describe the randomly distributed clouds, but constrained with θ ≤ θ_o (see Pancoast et al. 2014, for details).

Appendix D: Response of BLR clouds to the central source

For the response of BLR clouds to the central BBH systems co-rotating in elliptical orbits, the equations are similar to that with BBHs co-rotating in circular orbits as introduced in our previous work (Ji et al. 2021a), except for differential rotating time in each orbital phase and the variation of projected major axis of the elliptical orbit.

For the continuum emission dominated by the active secondary black hole, the time in the observer’s frame (t_obs) is determined by the intrinsic time (t_in) in the BBH mass center rest frame (hereafter BBH frame) and the cosmological time dilation:

$$\begin{array}{cc}\hfill t_{\mathrm{obs}}-t_{\mathrm{obs},0}& =[(t_{\mathrm{in}}-t_{\mathrm{in},0})+\frac{\overrightarrow{\mathit{OS}}·\widehat{\mathbf{l}}}{c}](1+z).\hfill \end{array}$$(D.1)

Here t_obs, 0 is the starting time for the observation; t_in and t_in, 0 are the times in the BBH frame corresponding to t_obs and t_obs, 0, respectively; $\overrightarrow{\mathit{OS}}$ is the position vector of the secondary black hole at the time t_in; and $\widehat{\mathit{l}}$ is the unit vector along the line of sight. For convenience, we set t_obs, 0 = 0 and t_in, 0 = 0 when the secondary BH is located at the apocenter of the elliptical orbit.

For photons emitted from a cloud C in the BLR located at (x₀, y₀, z₀), the observation time t_obs is related to the intrinsic time $t_{\mathrm{in}}^{\prime }$ in the BBH frame when the secondary BH is located at S′:

$$\begin{array}{cc}& t_{\mathrm{obs}}-t_{\mathrm{obs},0}=[(t_{\mathrm{in}}^{\prime }-t_{\mathrm{in},0})+\frac{\overrightarrow{\mathit{OC}}·\widehat{\mathbf{l}}+|\overrightarrow{{\mathit{S}}^{\prime }\mathit{C}}|}{c}](1+z).\hfill \end{array}$$(D.2)

Here $\overrightarrow{\mathbf{OC}}$ is the position vector the cloud C and $\overrightarrow{{\mathit{S}}^{\prime }\mathit{C}}$ is the distance between the secondary BH and the cloud C. At a specific moment of observation t_obs, we can derive the emission time difference between the continuum observed by the observer and that received by the BLR cloud as

$$\begin{array}{c}\hfill \tau =t_{\mathrm{in}}-t_{\mathrm{in}}^{\prime }=(\overrightarrow{\mathbf{OC}}·\widehat{\mathbf{l}}+|\overrightarrow{{\mathit{S}}^{\prime }\mathit{C}}|-\overrightarrow{\mathbf{OS}}·\widehat{\mathbf{l}})/c.\end{array}$$(D.3)

The above equations are the same for the BBH systems co-rotating both in circular and elliptical orbits.

For the modeling of the BBH+cBLR system under the Doppler boosting scenario, we assume the accretion rates are constant (e.g., D’Orazio et al. 2015; Duffell et al. 2020), whenever only the secondary BH is active or when both BHs are active. Five factors are taken into account to derive the finally observed BEL profiles (see details in Ji et al. 2021a,b):

• the position variation of the two BH components;

• the time-dependent Doppler boosting effect on the ionizing flux emitted from the primary or secondary BH accretion disks and recieved by the BLR clouds;

• the reverberation process of the BLR clouds to the central source(s);

• the gravitational redshift of the emission from BLR clouds (Tremaine et al. 2014);

• the Doppler blue- and redshifts of photons from the BLR clouds to the observer.

We assume that the intrinsic continuum flux emitted from the disk around the primary or secondary BH can be described by a power law $F_{\nu ,\mathrm{e}}(\nu ;t)\propto {\nu }^{\alpha }$, in which we assume the spectral index α = −2 (e.g., D’Orazio et al. 2015; Song et al. 2020, 2021; Ji et al. 2021a,b) of the flux emitted from each BLR cloud C_i can be derived. We can then obtain the BLR radiation (see Ji et al. 2021a, for more details)

$$\begin{array}{cc}\hfill L(v,t_{\mathrm{obs}})& \propto {{\displaystyle \sum _{i=1}^{N_{\mathrm{tot}}}}{F}_{\nu ,\mathrm{e}}(\nu ;t_{\mathrm{in}}^{\prime })D_{{\mathrm{C}}_{i}2}^{3-\alpha }\frac{|{\mathit{r}}_{{\mathrm{C}}_i}{|}^2}{|{\mathit{r}}_{{\mathrm{C}}_i}-{\mathit{r}}_{\mathrm{BH}}{|}^2}|}_{v=v_{{\mathrm{C}}_i}^{\mathrm{tot}}},\hfill \end{array}$$(D.4)

where N_tot is the total number of BLR clouds and $D_{{\text{C}}_{i}2}^{3-\alpha }$ represents the Doppler boosted or weakened ionizing flux $F_{\nu ,\mathrm{e}}(\nu ;t_{\mathrm{in}}^{\prime })$ emitted at $t_{\mathrm{in}}^{\prime }$ and received by the BLR cloud at t_obs, due to the relativistically orbital motion of the primary or secondary BH. The term $|{\mathit{r}}_{{\text{C}}_i}{|}^2/|{\mathit{r}}_{{\text{C}}_i}-{\mathit{r}}_{\text{BH}}(t_{\text{in}}^{\prime }){|}^2$ reflects the continuum flux variation due to the moving primary or secondary BH; $v_{{\mathrm{C}}_i}^{\mathrm{tot}}\simeq v_{{\mathrm{C}}_i}^{\mathrm{D}}+v_{{\mathrm{C}}_i}^{\mathrm{g}}$ represents the summation of the Doppler redshift or blueshift $v_{{\mathrm{C}}_i}^{\mathrm{D}}=\{{\left[{\gamma }_{{\mathrm{C}}_i\mathrm{obs}}(1-{\mathit{\beta }}_{{\mathrm{C}}_i\mathrm{obs}}·\widehat{\mathit{l}})\right]}^{-1}-1\}c$ and the gravitational redshift of photons emitted from an individual BLR cloud received by the distant observer $v_{{\mathrm{C}}_i}^{\mathrm{g}}=G{M}_{\bullet \bullet }/(r_{{\mathrm{C}}_{\mathrm{i}}}c)$ (Tremaine et al. 2014); β_Ciobs = v_Ciobs/c is the observed velocity of the cloud C_i in units of the light speed, ${\gamma }_{{\mathrm{C}}_i\mathrm{obs}}=1/\sqrt{1-|{\mathit{\beta }}_{{\mathrm{C}}_i\mathrm{obs}}{|}^2}$; and r_Ci is the distance vector of the BLR cloud to the mass center of the BBH system. In each model, the line emissivity from each cloud is assumed to be proportional to the flux received by the cloud.

All Tables

All Figures

Fig. 1. Sketch of the configuration of a BBH+cBLR system, with the BBHs and BLR clouds rotating in the same direction. For the BBHs co-rotating in elliptical orbits, the longitude of ascending node is set to be Ω = 0° for simplicity, i.e., the sky plane (yellow) and the BBH orbital plane (gray) share the same x-axis in the barycentric coordinate system. The argument of periapsis is denoted as ω. The observer’s line of sight is set in the y − z plane with an inclination angle of iobs to the BBH orbital plane. The BLR has an opening angle of θo and its middle plane is the same as the BBH orbital plane.

In the text

Fig. 2. Comparison of continuum light curves emitted by BBH systems when only the secondary BH is active with mass ratio qM = 0.1 (left panel) and 0.5 (right panel) at iobs = 60°. The black dashed line shows the light curve from the system with circular orbit (e = 0), the other four lines colored from blue to red show those resulted from elliptical orbits (e = 0.5), with four different position angles of their major axes ω varying from 0° to 270°. The horizontal dotted line shows the light curve without the Doppler Boosting effect.

In the text

Fig. 3. Comparison of continuum light curves emitted by BBH systems when both the BHs are active, with qL = 1 : 1, iobs = 60°, qM = 0.1 (left column) and 0.5 (right column). The top panel shows the observed continuum light curve contributed by Doppler Boosting effect from the two co-rotating black holes. The Doppler boosted or weakened flux modulated by the secondary and the primary BHs are shown in the middle and bottom panels, respectively. The lines in different shapes and colors are the same as in Fig. 2.

In the text

Fig. 4. BEL profile variations in a single orbital period for BBHs co-rotating in circular orbits with iobs = 60°. The two left columns show BBH systems when only the secondary BH is active, i.e., qL = 1 : 0, but with qM = 0.1 (left column) and 0.5 (middle-left column). The two right columns show profile variations when both the two BBHs are active with qL = 1 : 1, in the case of qM = 0.1 (middle-right column) and 0.5 (right column). In each column the top panel shows BEL profiles of the six phases (0 to 5π/3) in one period color-coded from blue to red (see legend), and the bottom panel presents the relative flux variation of BEL profiles compared to the mean profile obtained from that of the six phases.

In the text

Fig. 5. Comparison of continuum light curves emitted by BBH systems when only the secondary BH is active (top two rows) and both the two BH are active (bottom two rows) for e = 0.5 and iobs = 60°. In each row, the top panels show the profiles at six phases (0 to 5π/3) in one period, and the relative flux variations are shown in the corresponding bottom panels as that in Fig. 4. Columns from left to right show the cases with ω = 0° ,90° ,180°, and 270°, which reflects the profile variations caused by different eccentricity vectors of the elliptical orbits observed at a fixed viewing angle.

In the text

	Fig. 6. Legends are the same as those for Fig. 5, except that iobs = 30°.
In the text

Fig. 7. Variation of BEL amplitudes (ABEL) with different continuum amplitudes (Aconti) and orbital eccentricities (e) of BBHs with M•• = 109 M⊙ at ω = 0°. The two left columns show ABEL as a function of Aconti (left panel) and e (middle-left panel) for the case when only the secondary BH is active (qL = 1 : 0), and the two right columns show ABEL vs. Aconti (middle-right panel) and ABEL vs. e (right panel) for the case when both the BHs are active with qL = 1 : 1. In each panel, the points connected by lines of the same color show results of a fixed mass ratio, and the four different mass ratios qM = 0.1, 0.3, 0.5, and 0.7 are colored in blue, green, violet, and red, respectively. Each group of four points with different qM but the same e are linked by gray lines.

In the text

Fig. 8. Variation of continuum (Aconti) and BEL amplitudes (ABEL) with different inclination angles for BBHs with M•• = 109 M⊙, ω = 0°, and e = 0.0–0.6. The two left panels show the case when only the secondary BH is active (qL = 1 : 0), and the two right panels show the case when both the BHs are active with qL = 1 : 1. In each panel, the points colored from dark blue to red represent the orbital eccentricity of BBHs e = 0.0–0.6, respectively.

In the text

Fig. A.1. Radii of Roche lobes of BBHs and rotating velocity of the secondary BH at the orbital pericenter as a function of orbital eccentricity of the two BHs. The left and middle panels show the radii of Roche lobes ($R_{\mathrm{crit}}^{\mathrm{pri}}$ and $R_{\mathrm{crit}}^{\sec }$) in units of the Schwarzschild radii for the primary ($R_{\mathrm{Sch}}^{\mathrm{pri}}$) and secondary ($R_{\mathrm{Sch}}^{\sec }$) BHs, respectively. The typical radius of the UV emission (rUV) in the accretion disk for a M• = 108.6 M⊙ is shown as a horizontal dashed line for comparison. The right panel presents the rotating velocity of the secondary BH at the pericenter ($V_{\mathrm{peri}}^{\sec }$), which is the maximum velocity in each period. In each panel are shown the results of M•• = 108 (blue), 108.5 (green), and 109 M⊙ (red) with Torb = 2 yr and qM = 0.5. Considering that the UV emission might disappear when Rcrit ≲ rUV, the current model configuration only focuses on those systems with Rcrit ≳ rUV, i.e., e ≲ 0.6 for M•• = 109 M⊙ as a case study.

In the text