© The Authors 2025

1. Introduction

Massive black hole binaries (MBHBs) are predicted to form in the aftermath of galaxy mergers (Begelman et al. 1980). These objects are expected to be among the loudest gravitational wave (GW) sources that will be identified by future space-borne interferometers (e.g. Amaro-Seoane et al. 2023) and pulsar timing array campaigns (Verbiest et al. 2016). Recently, evidence of a GW background, possibly produced by supermassive black hole binaries, has been detected by PTA experiments (see, e.g. Agazie et al. 2023; Antoniadis et al. 2023). Identifying these sources simultaneously through electromagnetic (EM) observations will increase our knowledge of black hole inspirals, accretion physics, cosmology, and general relativity (see, e.g. Baker et al. 2019; De Rosa et al. 2019; Amaro-Seoane et al. 2023).

While there are confirmed EM detections of dual active galactic nuclei (AGNs) at ∼kiloparsec-parsec separations (see De Rosa et al. 2019 for a review, Trindade Falcão et al. 2024 for a recent candidate with a separation of 100 pc and Rodriguez et al. 2006 for a candidate with a separation of 7.3 pc), strong evidence of gravitationally bound MBHBs at sub-parsec scales is still missing. For this reason, several indirect observational signatures have been proposed (see e.g. Nguyen & Bogdanović 2016; D’Orazio & Haiman 2017; Nguyen et al. 2019; Dotti et al. 2023; D’Orazio & Charisi 2023, for recent discussions).

Among these signatures, the photometric ones are of particular interest as they will be searched for in forthcoming large observational campaigns, such as the Vera Rubin Observatory’s 10-year Legacy Survey of Space and Time (LSST, see Ivezić et al. 2019) and the Roman Space Telescope’s High Latitude Time Domain Survey (HLTDS, see Haiman et al. 2023), with the appropriate depth and cadence. Examples of these signatures are the variability in the continuum emission caused by the interaction between the binary and a circumbinary disc (see MacFadyen & Milosavljević 2008) and by relativistic effects, such as Doppler boosting (see D’Orazio et al. 2015) and gravitational self-lensing (see Kelley et al. 2021).

The length and cadence of these new time domain surveys will allow us to reject sources that present false periodicities due to red noise observed in optical and UV light curves (e.g. Vaughan et al. 2016) for short period binaries, with P ≲ 1 month corresponding to separations of a ∼ 10⁻⁴ (M/10⁶ M_⊙)^1/3 pc, where M is the binary total mass. Even though a fraction of false candidates can be rejected through long-term monitoring, the effective presence of a binary needs to be tested in detail as single massive black holes (MBHs) may have a variable intrinsic luminosity mimicking modulation in the light curves expected from binaries (see e.g. the discussion in Sandrinelli et al. 2016). Our work proposes a new spectroscopic test to prove the binary nature of photometrically identified sub-milliparsec MBHB candidates. At such small separations, the hypothetical binaries will have mildly relativistic orbital speeds, approaching 0.1c for the more massive and compact systems, and are thought to be embedded in a common circumbinary broad-line region (BLR). If at least one of the two MBHs is accreting, the periodically modulated ionising flux impinging on the BLR is expected to affect the broad emission line (BEL) shape in unique ways, due to the beamed anisotropic irradiation. The BEL time evolution then encodes information about the presence of a binary.

In our work, we assume disc-like BLRs, and we model both single and double-peaked BELs by changing the inclination of the system. Non-asymmetric BELs can be modelled by including deviations from axisymmetry in the BLR emissivity profile (see Storchi-Bergmann et al. 2017). Under this assumption, we compute how the BEL responds to the Doppler-modulated ionising continuum, which varies to first order as O(v/c), in the case of an MBHB, as well as to an intrinsically isotropically variable continuum emitted by a single MBH that has the exact same photometric modulation as is observed in the MBHB case. The following discussion is tailored for the broad Hβ line, but the proposed test can be readapted to other optical-UV BELs by adjusting the ‘characteristic radius’ of the BLR. Our scenario has some similarities with the model proposed by Hutchings & McCall (1977) to explain the 3-hour periodic variation of the observed magnitudes and emission lines in the Nova VI500 Cyg. Unlike our model, in which the modulation is caused by the anisotropic emission due to the Doppler boost, the periodic emission is explained by the varying presentation aspect of a hot spot in a conventional cataclysmic variable model (see Robinson 1976). In their model, the nebula, which plays a role similar to that of the BLR in our work, responds to the varying impinging continuum and presents a multi-component, axisymmetric structure consisting of a disc and two polar blobs. Similarly to the binary signature that we are proposing (see Section 3), the observed spectroscopic features of the nova emission lines cannot be explained under the assumption of an isotropic source.

In Section 2, we describe the underlying assumptions of our model, the methods with which the BELs are constructed, and how the fluxes of interest are computed. In Section 3, we illustrate our main results. Finally, in Section 4, we draw our main conclusions and propose the next steps to generalise our work to more complex scenarios.

2. Methodology

In our work, we compare the response of a disc-like BLR (detailed in Section 2.2) to the flux emitted by either (i) two MBHs bound in a close binary, whose properties are detailed in Section 2.1, under the assumption that each MBH has a constant luminosity and a radiation pattern shaped by Doppler boosting, or (ii) a single MBH at rest at the centre of the BLR, with a time-dependent luminosity such that the intrinsic emission in the continuum would be indistinguishable from the binary case (see Section 2.3). The difference between the binary and single MBH arises because the binary’s emission is beamed and anisotropic, illuminating different parts of the BLR at different times, analogous to a ‘lighthouse’ (D’Orazio & Haiman 2017). To assess whether this allows us to distinguish binaries from single MBHs, we computed the BEL shape generated by the impinging flux in the two scenarios (see Section 2.3), we found the light curves from the red and blue components of the BEL, as is discussed in Section 2.4, and finally we searched for the delays of the red and blue light curves relative to the ionising continuum. The last step was done by cross-correlating the three light curves using PyCCF (Sun et al. 2018)¹, the Python version of the original algorithm discussed in Peterson et al. (1998). This algorithm allows one to search for time lags between unevenly sampled light curves through linear interpolation. Also, the code allows one to compute the uncertainties associated with the time lags through Monte Carlo iterations. The test we propose is similar to the one discussed in D’Orazio & Haiman (2017), in which the (photometric) response of the torus is modelled for MBHBs as well as for single MBHs with modulated continua. The most important difference between the two studies is that our analysis can distinguish the response of different parts of the reprocessing structure (the BLR, in our scenario), encoded in the spectral profile of the BELs, as is detailed below.

In what follows, we consider stochastic noise caused by intrinsic AGN variability as discussed in Sect. 3, and we implicitly assume the limit of infinite signal-to-noise ratio (S/N) in the observed spectra. This last assumption will be relaxed in future analyses.

2.1. Binaries and Doppler boosting

In the binary scenario, we model a binary orbiting within a common BLR. Even if from hydrodynamical simulations, binaries are expected to have a certain degree of eccentricity (e.g. Roedig et al. 2011; Muñoz et al. 2019; Zrake et al. 2021; D’Orazio & Duffell 2021) depending on their mass ratio (see Siwek et al. 2023), in our work, we start by assuming circular binaries. The effect of binary eccentricity will be explored in a follow-up study.

We started by exploring the binary parameter space to identify the binaries that were most likely to be diagnosed with our proposed method. In the upper panel of Figure 1, we show the orbital period, P, of the binary. Looking at the source for multiple periods is important as it can prune out contaminants from intrinsic red noise AGN variability. As can be seen, at separations smaller than ∼10⁻³ pc the period is shorter than a year for most binaries.

Fig. 1. Upper panel: Binary period for binaries with different total mass and separations. The horizontal dashed red line indicates where the period is P = 1 yr. Middle panel: Coalescence time of the binaries with different masses (see the colours in the upper panel) and mass ratios, q. The solid lines show binaries with q = 0.1, dashed lines show binaries with q = 0.25, and finally, dash-dotted lines show binaries with q = 1. The horizontal dashed red line indicates where the coalescence time due to GW emission is τcoal = 100 yr. Lower panel: Ratio between the characteristic radius of the BLR and the binary separation for different masses. In all three panels, the segments of the lines that meet the conditions of interest in this work (binaries with total masses of around 106 − 107 M⊙, an orbital period of P < 1 yr, and a coalescence time due to GW emission of τcoal > 100 yr) are drawn with a thicker line. In the middle panel, only the case of q = 0.1 has a thick segment to maintain readability.

In the middle panel of Figure 1, we show the coalescence time through GW emission, τ_coal, of the binary at a given initial separation for different masses and mass ratios. This quantity is defined as (Maggiore 2008)

$$\begin{array}{c}\hfill {\tau }_{\mathrm{coal}}=\frac{5}{256}\frac{c^5{a}^4}{G^3{M}^2\mu },\end{array}$$(1)

with μ and a being, respectively, the reduced mass and the orbital separation of the binary. From this plot, it becomes clear that at separations of a megaparsec or smaller, binaries with a high total mass will be rare in observations as they will coalesce quickly. Also, the majority of AGNs that LSST will observe will be at their faint detection threshold, powered by lower-mass (M ∼ 10⁶ M_⊙) MBHs (see Xin & Haiman 2021). For these reasons, we focus on lighter (10⁶ − 10⁷ M_⊙) binaries.

The emission from the two MBHs is assumed to be produced by a Novikov-Thorne disc (see Novikov & Thorne 1973; Page & Thorne 1974)², with the discs around the MBHs having an outer radius given by the Roche lobe radius. Considering the secondary MBH, this radius is well approximated by (see Eggleton 1983)

$$\begin{array}{c}\hfill R_{\mathrm{RL},2}\approx 0.49a\frac{q^{\frac{2}{3}}}{0.6q^{\frac{2}{3}}+ln(1+q^{\frac{1}{3}})},\end{array}$$(2)

where q ≡ M₂/M₁ ≤ 1 is the binary mass ratio. However, for very thin, close-to-Keplerian discs, the outer radius is smaller by a factor of 4−5 (see Runnoe et al. 2015).

The luminosity emitted from the two discs depends on the mass accretion rate of the two MBHs, which is related to the binary mass ratio (see Duffell et al. 2020). Expressing the accretion rate of the two MBHs in terms of their individual Eddington ratios, $f_{\mathrm{Edd}}\propto \eta \dot{M}/M$, one finds

$$\begin{array}{c}\hfill \frac{f_1}{f_2}=\frac{{\eta }_1}{{\eta }_2}q(0.1+0.9q),\end{array}$$(3)

where f₁ and f₂ are the Eddington ratios and η₁ and η₂ are the radiative efficiencies of the primary and secondary MBH. Due to the limited constraints on the distribution of MBH spins in the following, we assume, arbitrarily, the normalised angular momentum to be $\stackrel{\sim }{a}=0.8J/Mc$ for both MBHs and that the accretion disc is co-rotating with the MBH. Thus, the radiative efficiency is η ≈ 0.12.

Since the two MBHs are moving, their emission is subjected to Doppler boosting. The MBH rest-frame frequency, ν, of a photon emitted by its accretion disc is related to the frequency, ν′, seen by an observer as ν′/ν = D, where D is the Doppler coefficient. This quantity can be written in terms of the relative velocity between the emitting source and the observer as

$$\begin{array}{c}\hfill D=\frac{1}{\gamma (1-\beta cos\theta )},\end{array}$$(4)

where θ is the angle defined in Equation (7) below, β ≡ (υ_em − υ_el)/c is the relative velocity between the emitter and the BLR element normalised to the speed of light, and finally $\gamma \equiv 1/\sqrt{1-{\beta }^2}$ is the Lorentz factor. In β, the emitter’s velocity direction is not the one seen at time t but the direction of motion at the retarded time of emission defined in Equation (8) below. The Doppler boosted ionising flux is given by

$$\begin{array}{c}\hfill F_{\mathrm{ion}}^{\mathrm{boost}}=D^4{\displaystyle {\int }_{{\nu }_{\mathrm{ion}}/D}^{\infty }}{F}_{\nu /D}\mathrm{d}\nu ,\end{array}$$(5)

where hν_ion = 1 Ry and ν is the rest-frame frequency.

In Figure 2, we illustrate the emission pattern in the binary and single MBH scenarios. In the single MBH case, the emission is isotropic, while in the MBHB scenario, the orbital motion beams the emission along the direction of motion. Since the orbital motion of the binary is periodic, the variation in the ionising flux seen by a given observer (or fluid element of the BLR in our case) is expected to be periodic, as well (e.g. D’Orazio et al. 2015).

Fig. 2. Illustration of the emission patterns in two scenarios. In the binary scenario (left), the emission from the accretion disc is beamed along the direction of motion, while in the single MBH scenario (right) the emission is the same in all directions.

2.2. Broad-line region

In this work, the BLR that surrounds the binary has a disc-like geometry³. The BLR disc and the binary are assumed to be coplanar for two reasons: (i) because in gas-rich galaxy mergers both the gas flowing towards the centre of the galaxy remnant and the two MBHs binding in a binary are expected to be affected by the angular momentum of the merger, so that the binary is expected to form already within a coplanar gas reservoir (e.g. Capelo & Dotti 2017), and (ii) because even if a misaligned binary-BLR forms, gravitational torques are expected to realign the two structures on ∼1 − 10 Myr timescales (e.g. Miller & Krolik 2013). We nevertheless briefly comment on the impact of a relative MBHB-BLR misalignment at the end of Section 2.4. The whole MBHB-BLR system is inclined with an angle, i, relative to the line of sight as shown in Figure 3.

Fig. 3. Illustration of the BLR geometry. Here Rin and Rout are the inner and outer radii of the BLR, while RBLR is its characteristic radius. L is the angular momentum of the system, l.o.s. is the line of sight to an observer, and finally i is the inclination angle between L and the line of sight. M1 and M2 are the masses of the two MBHs, while a is the binary separation. The binary separation and the BLR characteristic radius are not drawn to scale.

The characteristic radius of the BLR is defined (see Kaspi et al. 2000; Bentz et al. 2009; Dotti et al. 2023) as

$$\begin{array}{c}\hfill R_{\mathrm{BLR}}=c{\tau }_{\mathrm{LTT}}=11{\left(f_{\mathrm{Edd}}\frac{M}{{10}^6{M}_{\odot }}\right)}^{0.519}\mathrm{ld},\end{array}$$(6)

where τ_LTT is the light travel time across the BLR. Here we assume that τ_LTT is of the same order as the binary orbital period (P), ∼10 light-days in the example above. In such a case, we expect the response of the BEL profile to the orbit of the binary to be maximised. More specifically, if P ≪ τ_LTT the effect of the periodicity in the flux from the BLR would be averaged out over multiple periods, while if P ≫ τ_LTT, the length of an observational campaign aiming to detect a modulation in the BEL profiles will increase, and, for too large separations of the MBHB, the assumption of having a single BLR around the two MBHs will no longer be valid (see Dotti et al. 2023). In this work, the orbital period and the average light travel time are specifically assumed to be τ_LTT = 0.3 P.

The inner and outer radii are given by R_BLR, in = 0.5 R_BLR and R_BLR, out = 1.5 R_BLR. Once the binary masses and separation are chosen, R_BLR is fixed by the value of τ_LTT and the Eddington ratios for the two MBHs are set by combining Equations (3) and (6). The resulting values for the Eddington ratios of the individual BHs fall between 0.01 and ≈0.6, which is consistent with the observations (see Shankar et al. 2013).

In the lower panel of Figure 1, we show the ratio between the characteristic radius of the BLR and the binary separation. For lighter binaries (10⁶ − 10⁷ M_⊙), which are also the binaries of interest for our work as discussed in Section 2.1, the orbital separation is always much smaller than the characteristic radius of the BLR. For this reason, as well as for simplicity, we assume the binaries to be point-like sources at the centre of the BLR. Under these assumptions, the light travel time of a photon emitted by either of the two MBHs to reach an element of the BLR at a distance, r_el, will be well approximated by r_el/c and the effective emission time of a photon from the central source can be computed as in Equation (8) below.

Consider the Cartesian reference frame centred at the centre of the binary and also of the BLR shown in Figure 4. Assume that the inclination of the binary relative to an observer at an arbitrarily large distance, d₀, is the same as the inclination of the BLR disc, i, i.e. that the binary’s orbit is coplanar with its BLR disc, as is shown in Figure 3. The angle between the direction of motion, v_em, of an MBH and the direction towards a BLR element is given by

Fig. 4. Geometry of the system seen as face-on. The rotation of the angles relative to the Cartesian co-ordinate system is clockwise.

$$\begin{array}{c}\hfill {\theta }_{1,2}=\pi /2-{\theta }^{\prime }+{\phi }_{1,2},\end{array}$$(7)

where θ′ is the azimuthal angle of the BLR element with respect to a reference direction (the positive y axis in Fig. 4) and φ_1, 2 is the azimuthal position of the MBH (offset by an angle −π/2 from its orbital velocity vector).

The light emitted by the central pointlike source will reach the observer after a time of ∼d₀/c, where d₀ is the distance between the observer and the centre of the system. Since the BLR disc can be inclined relative to the line of sight, the light emitted at any time by the MBH’s accretion disc will reach the observer after being reprocessed by a BLR element with 2D co-ordinates (in the BLR plane), r′ and θ′, after a time delay, Δt, defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }t(r^{\prime },{\theta }^{\prime })=t-\frac{r^{\prime }}{c}(1+sinicos{\theta }^{\prime }),\end{array}$$(8)

where t ≡ d₀/c. The term r′sin i cos θ′/c comes from the effective distance of the element from the observer given by d(r′,θ′) = d₀ + r′ sin i cos θ′.

Our BLR model is simplified by the following assumptions: (i) In Equation (8), we are not taking into consideration the recombination time of the BLR, i.e. the time necessary for the line emission to respond to the new ionising continuum. The recombination time can be written as t_rec = (nα_r)⁻¹, where n is the number density and α_r is the recombination coefficient. Given the lower limit to the density n = 10¹² cm⁻³ given by the detection of strong and broad collisionally excited resonance lines in the UV (see Moloney & Shull 2014) and the typical BLR temperature (T ∼ 10⁴ K, resulting in a recombination coefficient, α_r ∼ 10⁻¹³ cm³ s⁻¹, e.g. Storey & Hummer 1995), the recombination time is expected to be t_rec ≲ 10 s, justifying our assumption;

(ii) We do not consider any effect due to structural changes in the BLR. Such changes would happen on the BLR dynamical time, t_dyn, of the order of the period of a circular orbit at the BLR radius. Within the assumptions of our study, the ratio between t_dyn and the binary period is t_dyn/P = (R_BLR/a_bin)^3/2, ranging from a few thousand to tens of thousands for the binary parameters considered here. Although a perturbation exerted by the MBHB on scales much smaller than R_BLR could in principle propagate to the BLR, the details of such propagation would be model-dependent and are outside the limited scope of this study.

Instead of solving for the time-dependent ionisation equilibrium, we assume a phenomenological description for the emissivity profile, ϵ(ξ, θ), which is able to reproduce the delays assuming the radius-luminosity relation discussed in Bentz et al. (2009). In our model, the BLR emissivity profile features a spiral pattern superimposed on an otherwise axisymmetric disc-like BLR. This prescription, originally proposed by Storchi-Bergmann et al. (2003) and recently modified by Sottocorno et al. (2025), is a proxy of a broad family of perturbations that make some portion of the disc brighter than the rest of the disc. The axisymmetric term is instead modelled as a Gaussian centred at a radius proportional to R_BLR following the radius-luminosity relation. The assumed emissivity profile is discussed in more detail in Appendix A.

2.3. BEL construction

The BEL profile observed at any given time was computed as the sum of the contributions from all the elements of the BLR. Dividing the BLR into N_r radial and N_θ′ azimuthal zones and sampling N_t times, N = N_t × N_r × N_θ′ computations (numerical integrations) of the ionising flux were needed to calculate the BEL shape evolution due to the Doppler-boost-driven anisotropic illumination of the BLR elements.

Since this operation is computationally expensive, we computed n ≪ N integrals for different values of the Doppler coefficient, D (see Eq. 4), and to fit the result using polynomials⁴. In Figure 5 we show an example of the results obtained for this fitting procedure for an MBH of 10⁶ M_⊙ at different Eddington ratios.

Fig. 5. Results for the fit of the relation between the boosted ionising flux, $F_{\mathrm{ion}}^{\mathrm{boost}}$, and the Doppler coefficient for an MBH of 106 M⊙ at Eddington ratios equal to fEdd = [0.25, 0.50, 0.75].

The algorithm proceeded as follows: (1) we evaluated D and the area (Σ) of the element, (2) the emissivity, ϵ(ξ, θ), of the BLR, and (3) the contribution of each element to the observed BEL at each time. This last step was sped up significantly by the fitting procedure for all the models considered. At the same time, step (2) benefited from our fitting too, as the chosen emissivity assumes that ξ_cent depends on the impinging ionising flux.

The observed flux for the considered BEL from the BLR at a given time, t, was computed as the sum of the contributions from each BLR element. We assumed that the contribution from each element is given by a Gaussian in the wavelength domain:

$$\begin{array}{c}\hfill F_{\mathrm{el}}(\lambda )\propto \frac{1}{2\pi \sqrt{{\sigma }_{\mathrm{line}}}}{e}^{-\frac{{(\lambda -{\lambda }_{\mathrm{obs}})}^2}{2{\sigma }_{\mathrm{line}}^2}}.\end{array}$$(9)

Here, the standard deviation of the Gaussian is given by σ_line = σλ_emit/c, where σ measures the local velocity dispersion due to small-scale turbulent motions within the BLR disc. In our work, this broadening parameter is given by σ = υ_Kep/10, where υ_Kep is the Keplerian velocity for that particular BLR element. Finally, λ_obs is the wavelength of the considered line in the rest frame properly shifted by the motion of the BLR element along the line of sight, λ_obs = λ_emit(1 + υ_los/c)⁵.

In this work, we do not consider general relativistic effects or higher-order special relativistic effects for the BLR elements, as their normalised distance is ξ ≫ 100 so that β ≪ 1. More specifically, the effects of the transverse redshift and the gravitational redshift of the elements of the BLR disc produce a fractional wavelength shift of the order of β², while the kinematic Doppler shift produces a fractional wavelength shift of the order of β ≫ β² so that the first two effects can be neglected. Also, we can neglect the effect of light bending for the light rays emitted from the BLR disc as the departure from a straight line is expected to be less than a degree (see Chen et al. 1989).

To take care of periodic modulations, we multiplied each Gaussian by the Doppler-boosted ionising flux, $F_{\mathrm{ion}}^{\mathrm{boost}}$, which strikes the element computed considering the correct retarded time, t′ (see Equation (8)), i.e. the current implementation we assume that each element responds instantaneously and linearly to the continuum flux reaching the element (i.e. we are implicitly assuming that the BLR clouds are optically thick and absorb all the incident ionising radiation). Assuming Σ to be the area of an element of the BLR and F_ion to be the observed ionising flux, the contribution from that particular element can be written as

$$\begin{array}{c}\hfill F_{\mathrm{el}}(\lambda )=F_{\mathrm{ion}}^{\mathrm{boost}}(t^{\prime },\xi ,{\theta }^{\prime })\frac{ϵ(\xi ,{\theta }^{\prime })\mathrm{\Sigma }}{2\pi \sqrt{{\sigma }_{\mathrm{line}}}}{e}^{-\frac{{(\lambda -{\lambda }_{\mathrm{obs}})}^2}{2{\sigma }_{\mathrm{line}}^2}},\end{array}$$(10)

where $F_{\mathrm{ion}}^{\mathrm{boost}}$ is the observed Doppler-boosted ionising flux computed by combining Equations (5), (7) and (8).

In the case of a single MBH at the centre of the BLR, we assume that the single MBH emits an intrinsic, periodically varying flux, F_ion, which mimics the ionising continuum light curve produced by the binary. The shape of the contribution from each element is the same as Equation (10) with the substitution of F_ion to $F_{\mathrm{ion}}^{\mathrm{boost}}$.

2.4. Light curve construction

The binary signature proposed in this work is based on the delays between the red and blue components of the BEL and the continuum light curve variations. One possible way to compute BEL-component light curves for this purpose is to divide the BEL in half, as is proposed in Dotti et al. (2023). Unfortunately, in the presence of a strong deviation from axisymmetry (such as a strong spiral perturbation) and for the short-period binaries considered in the current analysis, cutting the line in half might result in too large delays between the two BEL sides. As an example, consider the BLR geometry illustrated in Figure A.1, and assume that the region of the BLR with negative abscissa is receding from the observer, while the region with positive abscissa is approaching the observer. In this case, the spiral arm can introduce an arbitrary delay between the ionising continuum and the red light curve, as well as between the ionising continuum and the blue light curve. While the spiral arm perturbation does not affect the velocity field, the location where the bulk of emission is emitted in the receding and approaching sides of the BLR depends on the spiral parameters. To minimise the effect of the poorly constrained asymmetry of the BLR, we considered only the wings of the BEL – or, equivalently, we considered only the elements with a high line-of-sight velocity, highlighted in green in Figure A.1. For these elements, the light-crossing time difference is much smaller than half of a binary period, so that the time delays between the blue and red wings’ response will only be due to the phase shifts caused by the binary’s Doppler modulation. More specifically, the BEL was cut in the following way:

• Double-peaked lines: the line presents two maxima, one on the blue side (left of the rest-frame wavelength) and one on the red side (right of the rest-frame wavelength). We took the temporal mean of the BEL flux and defined two maxima, ${\overline{F}}_{\mathrm{blue}}$ and ${\overline{F}}_{\mathrm{red}}$, occurring at λ_max, blue and λ_max, red, respectively. The BEL was cut at the wavelengths λ_blue and λ_red, where the flux of the BEL is given by $\overline{F}({\lambda }_{\mathrm{blue}})={\overline{F}}_{\mathrm{blue}}/2$ and $\overline{F}({\lambda }_{\mathrm{red}})={\overline{F}}_{\mathrm{red}}/2$, with λ_blue < λ_max, blue and λ_red > λ_max, red. In Figure 6, we show an example of how a double-peaked line is cut.

  Fig. 6. Example of cuts for a double-peaked line. The black line represents the temporal mean of the BEL flux, while the vertical dashed red and blue lines represent the wavelengths at which the line is cut.

• Single-peaked line: the line presents only one maximum. Here, the wavelengths of the two cuts, λ_blue and λ_red, are found by finding the wavelengths where the mean flux is half its maximum value.

The two light curves (‘blue’ and ‘red’) were then computed by integrating the flux over the wavelengths below λ_blue and above λ_red, at each sampled time. Cutting the BEL as described above, the delays between the red and blue light curves are expected to be near zero in the case of a single MBH at the centre of the BLR, as in this scenario the light is emitted isotropically. In the case in which the central source is an MBHB, the emission is beamed along the direction of motion, and the delays between the red and blue light curves are expected to be nearly half the orbital period, independently of the properties of the spiral. This last expectation (confirmed in Section 3) depends on the assumption of a coplanar MBHB-BLR system. In the presence of a significant misalignment⁶ the beaming effect onto the BLR would be even clearer, as specific regions of the BLR would be exposed to Doppler-boosted flux only at two specific binary phases, when the MBH velocities point in the direction of the misaligned disc. The frequencies at which such an excess in the reprocessed luminosity would appear would, however, depend on the specific 3D orientation of the observer-binary-circumbinary disc. For circular orbit binaries (as considered in this analysis), there will always be two frequencies on opposite sides of the BEL rest frame centroid, but their exact location (i.e. whether they will be in the core of the lines or in their wings) would depend on the specific configuration.

3. Results

As we noted in Section 2.1, in the binary scenario, the emission from the two MBHs is beamed along the direction of motion of the emitter, while in the single black hole scenario, this does not happen (see Figure 2). Consider the two innermost elements of the BLR along the line of nodes. If the radiation is beamed and one of the two MBHs emits much more than the other, one element (e.g. the element moving towards the observer) is expected to reprocess the maximum impinging ionising flux after a certain delay given by the BLR light travel time. The receding element is expected to reprocess the maximum ionising continuum half an orbital period after the approaching element. This means that, given the choice of BLR extent assumed here, the response of the red part of the BEL will have a time lag bigger than half a period relative to the ionising continuum light curve. In other words, the peak of the red response is nearer to the next observed peak of the continuum emission (i.e. the boosted contribution of the less luminous MBH), and thus the delay can be seen as negative. The lags of the blue and red light curves relative to the ionising continuum will have opposite signs: τ_blue > 0 and τ_red < 0. Cutting the line as explained in Section 2.4, we focus on the contribution from the inner part of the BLR near the line of nodes, and thus we expect to see the effect discussed in the example above by measuring the lags between pairs of light curves: the ionising continuum and the red light curve, as well as the ionising continuum and the blue light curve. Depending on the ratio between the BLR size and binary separation, there can be cases in which the delays have the same sign but their difference is nevertheless expected to be around half the binary period.

In the single black hole scenario, the light is not beamed. This means that both the approaching and receding elements will reprocess the impinging ionising flux nearly at the same time. Although the presence of a spiral arm in the BLR can introduce some differences, causing slight variations in the delays relative to the ionising continuum light curve, these lags will still be comparable and will have the same signs, and, in general, their difference should be negligible.

The lags between the pair of light curves can be combined in a parameter defined as

$$\begin{array}{c}\hfill \chi =\frac{|{\tau }_{\mathrm{red}}-{\tau }_{\mathrm{blue}}|}{|{\tau }_{\mathrm{red}}+{\tau }_{\mathrm{blue}}|}·\end{array}$$(11)

This quantity is sensitive to the relative sign between the two time lags, τ_red and τ_blue. If the signs are discordant, χ will be larger than 1, while if the signs are concordant, χ will be smaller than 1. In the binary scenario, χ is expected to spread around 1 and in most cases larger than 1, depending on the size of the BLR, while in the single MBH scenario, χ is expected to be around zero, regardless of the size of the BLR. This quantity discriminates between binaries and single MBHs with varying intrinsic luminosities:

$$\begin{array}{cc}\hfill \chi & ⪆1\Rightarrow \mathrm{binary}\hfill \\ \hfill \chi & \approx 0\Rightarrow \mathrm{single}\mathrm{MBH}.\hfill \end{array}$$(12)

In Figure 7 we show examples of light curves in the two cases. The upper panel of Figure 7 shows the binary scenario (with M₁ = 10⁷ M_⊙, q = 0.1 and f_Edd, 1/f_Edd, 2 ≈ 0.02) where, as expected, the blue light curve responds with a positive delay relative to the continuum light curve, in this case χ = 1.38. The lower panel of Fig. 7 shows the case of a single MBH with a varying luminosity; in this case, the red and blue light curves respond at nearly the same time, yielding χ = 0.004. In Figure 8 we show the distribution of χ in the two scenarios, keeping the same MBH masses as in Figure 7 and changing the BLR characteristic radius and the spiral parameters. We set R_BLR, so that the light travel time varies uniformly between 0.1 and 0.5 times the orbital period, while the spiral parameters vary uniformly in the ranges A = [0, 10], p = [ − π/2, π/2], δ = [0, 2π], and ϕ₀ = [0, 2π] (see Appendix A). When considering the binary scenarios of these realisations, we find χ < 0.5 in 0.6% of the realisations. In the single MBH scenario, instead χ ≈ 0, with χ > 0.05 in only 0.07% of the realisations.

Fig. 7. Light curves for the ionising continuum from the central source and for the red and blue BEL wings. The x axis represents time in units of the orbital period, while the y axis represents the normalised flux. Upper panel: Binary scenario with M1 = 107 M⊙ and M2 = 106 M⊙. The separation is a = 10−3.9 pc. Lower panel: Single MBH scenario with the same ionising continuum as in the binary case.

Fig. 8. Distribution of the statistic χ (Eq. 11) reflecting the difference in time delays in the blue and red BEL wings, in the single MBH versus the binary scenarios.

The test we propose works for unequal-mass binaries with a mass ratio of q ≲ 0.3 but becomes observationally challenging as the mass ratio increases. This is because, to first order, the amplitude of the variability induced by the Doppler boost decreases with an increasing mass ratio, as it results from the sum of sinusoids with similar amplitude but out of phase by half a period. If the two black holes emit similar spectra with small local spectral curvature (e.g. if the intrinsic emissions have power-law spectral energy distributions), then this results in a near-cancellation of the Doppler modulation. Concerning equal-mass binaries, the test does not work as discussed so far, but can be re-adapted if the BLR shows an axisymmetric emissivity profile. Specifically, a region closer to the bulk of the emission line needs to be considered to construct the red and blue light curves. A study of the effect of mass ratio is shown in Appendix B.

So far, the effect of intrinsic AGN variability has been neglected. Since we are considering MBH binaries, we expect that, in addition to the variability due to the Doppler boost, the continuum light curve, and consequently the red and blue fluxes from the BLR response, will all be affected by a stochastic red noise in the emission from both MBHs. The AGN variability is described well by a red noise, often modelled as a damped random walk (DRW, see Kelly et al. 2009). Since the DRW is an intrinsic property of the continuum source, its effect should be enhanced or diminished by the Doppler boost. The implementation of this effect would ideally enter the Doppler-boosted flux computation. We do not compute the Doppler-boosted flux through a double integration, as it would be too computationally expensive, but through a fitting procedure. Because of this, we cannot implement the DRW before the application of the Doppler boost. To study the effect of the DRW, we implement it more simply as a relative variation in the already Doppler-boosted flux. In Figure 9, we show examples of time series obtained by sampling exclusively from a DRW along with a comparison between the expected power spectrum and the periodogram computed from the sampled DRW time series. These were used to verify that the chosen sampling of the DRW light curve does not alter the power spectrum of the fluctuations.

Fig. 9. Left: Illustrative examples of DRW time series with different damping timescales and variation amplitudes. Here SFinf is the structure function of the process at infinitely large separations in time and is related to the root mean square brightness fluctuation amplitude, σCAR, by SF${}_{inf}=\sqrt{2}{\sigma }_{\mathrm{CAR}}$. Right: Expected power spectrum for the DRW superimposed on the periodogram obtained from the green time series of the left image. Note that in both panels, the y axis is reported in arbitrary units.

The parameters we assumed to model the DRW variability are: the mean value of the stochastic process, μ = 1, and the characteristic timescale and amplitude of the fluctuation, τ = [100, 500] d and σ_CAR = [0.21, 0.35] for light curves in Figure 9, τ = [30, 100] d and σ_CAR = [0.03, 0.5] for light curves in Figure 10, and finally τ = [100, 30] d and σ_CAR = [0.03, 0.01] for the MBHs whose light curves are shown in Figure 11 below. The light curves were sampled using the time_series.generate_damped_RW function from the AstroML library. We then multiplied the boosted flux from the two MBHs by these values to build the combined light curves, showing good agreement.

Fig. 10. Examples of light curves generated considering both the effects of the Doppler boost and of the variability induced by the damped random walk. In the lower panel, we show the case in which the S/N is 0.28, while in the upper panel, the S/N is 150, and the periodic variability induced by the Doppler boost can be clearly seen.

The main signature that will identify a binary candidate at separations smaller than 10⁻³ pc is the periodic continuum variability. Red noise introduces two problems. First, it might produce fake periodicities. A distinction between a real periodicity and a fake one can be found by observing the object for a long enough time (see Xin & Haiman 2024). The second problem is associated with the S/N. If the variability caused by the binary nature of the source is weaker than the variability induced by red noise, it will be difficult to identify the object as a binary candidate using its continuum emission. This implies that binaries with separations small enough to have a Doppler boost variability stronger than the noise will have a higher chance of being identified.

In Figure 10 we show two light curves in the case of S/N < 1 and S/N > 1, where we consider as a signal the average over five orbital periods of the flux obtained by our model ($\overline{f}(t)$), while the noise was computed as the variance of the noisy light curve,

$$\begin{array}{c}\hfill s={\displaystyle \sqrt{\frac{\sum _{i=1}^N{(f(t)-\overline{f}(t))}^2}{N}}},\end{array}$$(13)

where f(t) is the observed noisy flux and N is the number of samples.

Once the period of variation is determined, this information can be used to apply the test proposed in this work to determine whether the source of the light curve is a binary or a single MBH. Specifically, phase-folding the light curve over multiple cycles enhances the periodic features, reduces noise, and allows for a more accurate estimation of the delays between pairs of light curves and their associated uncertainties. Also, by phase-folding the light curve, it is possible to determine whether the periodicity is physical or induced by noise.

In Figure 11 we show the results of such a procedure for a binary with q = 0.1, M = 10⁷ M_⊙, a = 10^−4.2 pc observed for 5, 10, and 20 periods, with the S/Ns being, respectively, 181, 222, and 242 (a binary with these parameters would merge in about 100 yr; thus, it would be rarely observed). As can be seen, after phase-folding the light curves become smoother, and uncertainties can be associated with them, allowing for the determination of errors in the measurement of the delays and of the χ parameter. From the plots, it is clear that the red and blue fluxes show a smoother behaviour compared to the ionising light curve. This is caused by the implementation we are using. As is discussed in Section 2.4, the red and blue light curves were computed from an emission line built as a sum of the contributions from each element of the BLR. In Section 2.2, we showed that each BLR element responds to the ionising flux emitted at a time that depends on the position of the element. This means that the red and blue fluxes already smooth the noisy continuum, explaining the behaviour shown in the light curves. In this implementation, we did not consider the effect of the S/N due to the measurement noise of such a small part of the spectrum. Taking into account this effect, the uncertainties associated with the red and blue light curves are expected to increase.

Fig. 11. Mean ionising continuum, red and blue light curves, and associated uncertainties resulting from the phase-folding procedure observing the source for 5, 10, and 20 periods, respectively (left to right panels), for a binary with q = 0.1, M = 107 M⊙, a = 10−4.2 pc.

As was mentioned in Section 2, PyCCF can be used to estimate the uncertainties on the time lags between different light curves. Thus, one can use the light curves shown in Figure 11 to compute the time lags and the uncertainties associated with them, and the χ parameter to assess whether the source is a binary. In Figure 12 we show the normalised distribution of χ, computed through Monte Carlo simulations, for a noisy light curve observed for 20 periods. Here, the points of the light curves are resampled in a way similar to the one described in Peterson et al. (1998), for the light curves shown in the lower plot of Figure 11, showing that the median value of χ = 1.67 correctly identifies the source as a binary.

Fig. 12. Normalised χ distribution found through Monte Carlo resampling of fluxes within their observed uncertainties to estimate the median value and uncertainty for a binary observed for 20 periods with q = 0.1, M = 107 M⊙, a = 10−4.2 pc.

4. Conclusions and prospects

The goal of this work was to develop an observational signature for MBHBs at orbital separations of the order of ≤10⁻³ pc where a binary is thought to be surrounded by a common BLR. We studied the effect of the Doppler boost caused by the orbital motion of the binary on the BEL evolution. We assumed that the two MBHs are on a circular orbit and that the common BLR has a disc-like structure with a non-axisymmetric emissivity (see Appendix A), where the asymmetry in the BLR is modelled as a spiral arm superimposed on an axisymmetric Gaussian emissivity profile (see Sottocorno et al. 2025). The spiral arm can be used to model asymmetric and double-peaked lines (Storchi-Bergmann et al. 2003), and could, in principle, imprint a different time delay in the different regions of the BEL.

Since the BELs respond to the impinging luminosity emitted by the central engine, information about the nature of the source is expected to be encoded in the shape of the BELs and can be found in reverberation mapping campaigns (with a ∼daily cadence). In this work, we compared the shape evolution of the Hβ line in two scenarios: in the first case, the central source is a binary, while in the other, the central source is a single MBH with a varying intrinsic luminosity that mimics the continuum emitted in the binary scenario. We focus on the wings of the BEL, i.e. the regions of the line that are mostly affected by the BLR elements with a high line-of-sight velocity. Doing so reduces the effect of the spiral arm on the emissivity. Integrating the total flux in the selected regions of the BEL, we constructed red and blue light curves and computed the time lags between either of these two time series, relative to the ionising continuum light curve coming directly from the central source.

The main result of this work is that the presence of a binary is encoded in these time lags. Defining the parameter χ = |τ_red − τ_blue|/|τ_red + τ_blue|, one finds that if χ ≳ 1, the central object is a MBHB, while if χ ≈ 0, the central source is a single MBH.

After discussing the effectiveness of the test, we introduced stochastic AGN variability, modelled as a damped random walk, to the synthetic light curves. The results show that if the period of variability can be estimated, the light curve can be phase-folded to enhance the periodic features, smooth out the noise, and determine whether the observed periodicity is physical or whether it results from long-term correlations induced by red noise. We show that if the variability is caused by physical mechanisms, the proposed test can correctly identify the source as a binary.

In Appendix B, we show how the results change by varying the binary mass ratio. The main conclusion is that the signature we propose works efficiently for unequal mass binaries with a mass ratio of less than q ∼ 0.3. As q approaches unity, the use of this signature becomes observationally challenging as the sum of the contributions to the continuum from the two MBHs would make the first-order effect of the Doppler boost smaller and smaller by comparison. In the limiting case of an equal mass binary, the continuum is expected to be constant to first order. Thus, the only variability that is shown in Appendix B is due to second-order effects that are too small to be observed in the presence of noise. Although realistically such effects would be difficult to observe, we showed that the test would still work for unequal mass binaries, while, in the case of an equal mass binary the test can be adapted for an axisymmetric BLR.

The signature proposed in this work must be tested in more general scenarios. First, the Doppler boost is not the only source of variability expected for MBHBs: hydrodynamical simulations show that the accretion rate from a circumbinary disc changes periodically and effects of general relativity, such as gravitational lensing, should be implemented. Second, MBHBs are expected to have some eccentricity, while we considered only circular binaries in this work. Different geometries and emissivity profiles can be used to model the BLR, while more realistic accretion disc models can be implemented, and the effect of the S/N of such a small part of the spectrum should be considered. We plan to test the proposed signature in such scenarios in future work.

To date, we have not applied the test to real data, as we are not aware of any reverberation mapping campaign on an MBHB candidate selected through a periodically modulated light curve. Multi-epoch spectroscopic campaigns studying AGNs with strong evidence of periodic light curves in current and future time-domain surveys will be the natural test bed for our newly proposed test.

Acknowledgments

MD acknowledge funding from MIUR under the grant PRIN 2017-MB8AEZ, and financial support from ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU. This study is supported by the Italian Ministry for Research and University (MUR) under Grant ‘Progetto Dipartimenti di Eccellenza 2023-2027’ (BiCoQ). FR acknowledges the support from the Next Generation EU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures, Project IR0000012 – CTA+ – Cherenkov Telescope Array Plus. DJD acknowledges support from the Danish Independent Research Fund through Sapere Aude Starting Grant No. 121587 ZH acknowledges support from NSF grant AST-2006176 and NASA grants 80NSSC22K0822 and 80NSSC24K0440.

References

Appendix A: BLR emissivity profile

The BLR emission is modelled using the modification of the prescription of Storchi-Bergmann et al. (2003) proposed by Sottocorno et al. (2025) (see also Rigamonti et al. 2025 for an application of the model to real data). In this model, the BLR emissivity features a spiral pattern superimposed on an otherwise axisymmetric disc-like BLR. More specifically, the assumed emissivity profile is given by

$$\begin{array}{cc}\hfill ϵ(\xi ,\theta )=& \frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{cent}}}\frac{1}{\xi }exp[-\frac{{(\xi -{\xi }_{\mathrm{cent}})}^2}{2{\sigma }_{\mathrm{cent}}^2}]\{1+\frac{A}{2}{exp}^{[-\frac{4ln2}{{\delta }^2}{(\varphi -{\psi }_0)}^2]}\hfill \\ \hfill & +\frac{A}{2}exp[-\frac{4ln2}{{\delta }^2}{(2\pi -\varphi +{\psi }_0)}^2]\},\hfill \end{array}$$(A.1)

where ξ is the radial distance from either the single MBH or the MBHB centre of mass in units of gravitational radii⁷, while ϕ is the azimuthal co-ordinate in the BLR plane. Here, the spiral perturbation is intended as a proxy of a broad family of perturbations that make some portion of the disc brighter than the rest of the disc. Spirals are however particularly appealing, as their development in discs are observed in a variety of systems, from spiral galaxies to discs around cataclysmic variables (e.g. Patterson et al. 1993; Steeghs et al. 1997), and their presence allows for the description of asymmetric broad lines as seen in observations (e.g. Storchi-Bergmann et al. 2017; Ward et al. 2024, 2025, and references therein). The axisymmetric term of the profile is given by a Gaussian, centred at the normalised radius ξ_cent assumed to be equal to the characteristic radius of the BLR (ξ_cent = R_BLR), divided by the normalised distance. Using this profile, we automatically obtain a characteristic radius of the BLR consistent with the empirically inferred luminosity-radius relation (Bentz et al. 2009), once the relation between ξ_cent and ionising source luminosity is specified (see Sottocorno et al. 2025).

The spiral arm is parametrised by its azimuthal width (or FWHM, δ) and an azimuthal reference co-ordinate ψ₀(ξ). Here ψ₀ = ϕ₀ + log(ξ/ξ_sp)/tan(p), whereϕ₀ is the azimuthal position of the spiral at its innermost normalised radius ξ_sp and p is the spiral pitch angle. In this work, we assume ξ_sp to coincide with the innermost radius of the BLR. Finally, the parameter A represents the brightness contrast between the spiral arm and the underlying, axisymmetric disc, so that in total the spiral shape of the BLR is specified by five parameters. In Figure A.1, we show an example of the emissivity profile with a spiral pattern for the set of parameters reported in Table A.1.

Fig. A.1. Example of a BLR with a spiral arm. The emissivity is reported in units of its minimum value (on a linear colour gradient) and is set to 0 in regions where the BLR is not present. The grid elements responsible for the considered red and blue parts of the BEL are highlighted in green (see Section 2.4).

Appendix B: Effect of the mass ratio

In our model, the binary mass ratio determines the accretion rate of the two MBHs. The Eddington ratios at which the two MBHs are emitting are fixed by equation 3. As the mass ratio approaches unity, the emission from the two MBHs is therefore more and more similar. Hence, as q increases, the emission from the MBH with the lower accretion rate will become more important. Since the emission from the two MBHs will vary sinusoidally with the same frequency, the orbital one, but with a phase difference of π, the observed light curves will have smaller and smaller amplitudes. If the amplitude becomes small enough, second-order effects will become visible. In the context of Doppler boosting for a binary, the second-order effect is a sinusoid with a periodicity equal to half the orbital period. Thus, when summing the contributions from the two MBHs, both the continuum and BEL light curves will not necessarily look like sinusoids with a period equal to the orbital period of the binary anymore. Since the second-order effects are small, identifying such systems as binary candidates by the observation continuum light curves will be difficult in real scenarios, as the instrumental noise, as well as the intrinsic AGN variability, may hide the effect of the Doppler boost. Figure B.1 shows the light curves for binaries (left panels) with a separation a = 10^−4.2 pc with the primary MBH mass of 10⁶ M_⊙ and mass ratios q = [0.25, 0.55, 0.75] and compares them with the BLR response to a single MBH showing the same continuum light curve (right panels). The y-axis range changes to show smaller and smaller relative variations for increasing values of q.

Fig. B.1. Light curves from BELs for binaries with a = 10−4.2 pc, M1 = 106 M⊙ and mass ratio q = [0.25, 0.55, 0.75].

When the light curves resemble the example shown in the middle and lower rows of Figure B.1, the presence of two very similar peaks in the light curves, introduced by second-order effects, could prevent an efficient time lag computation with PyCCF. In particular, in the presence of noise, as the mass ratio approaches unity, a higher S/N would be required to correctly compute the time delays, making the proposed test more observationally challenging. A critical scenario is found when an equal-mass binary is considered. In this case, to first order, the continuum light curve will not be variable, and the only variability that will be visible in the absence of noise is the one induced by the second-order effects. Considering this variability, the binary and single black hole scenarios are almost indistinguishable.

More detailed analyses of the time evolution of the whole BEL (and of a varying continuum) might still identify the presence of a binary. In the following, we propose a specific test for the limiting q = 1 case in the case of axisymmetric BLRs (i.e. in the absence of a significant spiral pattern in the BLR emissivity).

The test can, however, be adapted for equal mass binaries in the specific case of symmetric BLRs. Consider, for example, two opposite elements of the BLR at the same distance R_el from the central source misaligned with respect to the line of nodes, as shown in the lower panel of figure B.2. In the binary scenario, when one of the two MBHs is illuminating the farthest approaching element of the BLR, the other black hole will be illuminating the nearest receding element, inducing different delays between the red and blue light curves relative to the ionising continuum. In the single MBH scenario, where the light is not beamed, no delay between the red and blue light curves is expected. An example of the resulting light curves obtained selecting the BEL wavelength range highlighted in the upper panel of figure B.2 is shown in figure B.3.

Considering the application of the test to real observations, we suggest that the test should be performed by cutting the BEL in different regions, as the binary mass ratio q is not known a priori, and the choice of the region where the line should be cut is not straightforward. Considering different BEL regions, the light curves built considering the tails are expected to be more sensitive to unequal mass binaries, while the light curves built using a cut similar to the readapted test described above are expected to be more sensitive to equal mass binaries.

Another point that is worth stressing is that testing the presence of a binary by cutting the line as discussed in section 2.4 works particularly well for the BLR size we assumed (see section 2.2). When the inner and outer radii of the BLR are changed, the fraction of the maximum flux at which we find the limiting wavelength of the cut can be different from ${\overline{F}}_{\max (\mathrm{red}/\mathrm{blue})}/2$. For example, consider the case in which the outer radius of the BLR is bigger than 1.5 R_BLR. In this scenario, we expect more BLR elements to contribute to the BEL flux near the rest-frame wavelength. Cutting the line as in 2.4, the resulting light curves are affected by elements with a lower line-of-sight velocity relative to the elements of interest. Considering unequal mass binaries, performing the test by cutting the BEL at different wavelengths can be beneficial in finding the regions that are mostly affected by elements with a high line-of-sight velocity. In particular, when the outer radius of the BLR is bigger than the one we assumed, the test can be performed multiple times cutting the line at wavelengths in which the BEL flux is a smaller and smaller fraction of ${\overline{F}}_{\max (\mathrm{red}/\mathrm{blue})}$ up until the considered region is not dominated by noise.

Fig. B.2. Upper panel: Example of a cut in a different region of the BEL that is used to highlight the presence of a binary in the equal-mass case. Lower panel: qualitative example of the elements to be considered to highlight the presence of a binary in the case of an equal-mass binary with no spiral pattern. The red elements are the ones that are used in the method proposed for unequal-mass binaries, while the green ones are the ones that will encode the binary signature for the equal-mass binary.

Fig. B.3. Example of how a cut in a different region of the BEL (those highlighted in figure B.2) can highlight the presence of a binary with the red and blue light curves separated by a half-period lag.

All Tables

All Figures

Fig. 1. Upper panel: Binary period for binaries with different total mass and separations. The horizontal dashed red line indicates where the period is P = 1 yr. Middle panel: Coalescence time of the binaries with different masses (see the colours in the upper panel) and mass ratios, q. The solid lines show binaries with q = 0.1, dashed lines show binaries with q = 0.25, and finally, dash-dotted lines show binaries with q = 1. The horizontal dashed red line indicates where the coalescence time due to GW emission is τcoal = 100 yr. Lower panel: Ratio between the characteristic radius of the BLR and the binary separation for different masses. In all three panels, the segments of the lines that meet the conditions of interest in this work (binaries with total masses of around 106 − 107 M⊙, an orbital period of P < 1 yr, and a coalescence time due to GW emission of τcoal > 100 yr) are drawn with a thicker line. In the middle panel, only the case of q = 0.1 has a thick segment to maintain readability.

In the text

	Fig. 2. Illustration of the emission patterns in two scenarios. In the binary scenario (left), the emission from the accretion disc is beamed along the direction of motion, while in the single MBH scenario (right) the emission is the same in all directions.
In the text

Fig. 3. Illustration of the BLR geometry. Here Rin and Rout are the inner and outer radii of the BLR, while RBLR is its characteristic radius. L is the angular momentum of the system, l.o.s. is the line of sight to an observer, and finally i is the inclination angle between L and the line of sight. M1 and M2 are the masses of the two MBHs, while a is the binary separation. The binary separation and the BLR characteristic radius are not drawn to scale.

In the text

	Fig. 4. Geometry of the system seen as face-on. The rotation of the angles relative to the Cartesian co-ordinate system is clockwise.
In the text

	Fig. 5. Results for the fit of the relation between the boosted ionising flux, $F_{\mathrm{ion}}^{\mathrm{boost}}$, and the Doppler coefficient for an MBH of 106 M⊙ at Eddington ratios equal to fEdd = [0.25, 0.50, 0.75].
In the text

	Fig. 6. Example of cuts for a double-peaked line. The black line represents the temporal mean of the BEL flux, while the vertical dashed red and blue lines represent the wavelengths at which the line is cut.
In the text

Fig. 7. Light curves for the ionising continuum from the central source and for the red and blue BEL wings. The x axis represents time in units of the orbital period, while the y axis represents the normalised flux. Upper panel: Binary scenario with M1 = 107 M⊙ and M2 = 106 M⊙. The separation is a = 10−3.9 pc. Lower panel: Single MBH scenario with the same ionising continuum as in the binary case.

In the text

	Fig. 8. Distribution of the statistic χ (Eq. 11) reflecting the difference in time delays in the blue and red BEL wings, in the single MBH versus the binary scenarios.
In the text

Fig. 9. Left: Illustrative examples of DRW time series with different damping timescales and variation amplitudes. Here SFinf is the structure function of the process at infinitely large separations in time and is related to the root mean square brightness fluctuation amplitude, σCAR, by SF${}_{inf}=\sqrt{2}{\sigma }_{\mathrm{CAR}}$. Right: Expected power spectrum for the DRW superimposed on the periodogram obtained from the green time series of the left image. Note that in both panels, the y axis is reported in arbitrary units.

In the text

	Fig. 10. Examples of light curves generated considering both the effects of the Doppler boost and of the variability induced by the damped random walk. In the lower panel, we show the case in which the S/N is 0.28, while in the upper panel, the S/N is 150, and the periodic variability induced by the Doppler boost can be clearly seen.
In the text

	Fig. 11. Mean ionising continuum, red and blue light curves, and associated uncertainties resulting from the phase-folding procedure observing the source for 5, 10, and 20 periods, respectively (left to right panels), for a binary with q = 0.1, M = 107 M⊙, a = 10−4.2 pc.
In the text

	Fig. 12. Normalised χ distribution found through Monte Carlo resampling of fluxes within their observed uncertainties to estimate the median value and uncertainty for a binary observed for 20 periods with q = 0.1, M = 107 M⊙, a = 10−4.2 pc.
In the text

	Fig. A.1. Example of a BLR with a spiral arm. The emissivity is reported in units of its minimum value (on a linear colour gradient) and is set to 0 in regions where the BLR is not present. The grid elements responsible for the considered red and blue parts of the BEL are highlighted in green (see Section 2.4).
In the text

	Fig. B.1. Light curves from BELs for binaries with a = 10−4.2 pc, M1 = 106 M⊙ and mass ratio q = [0.25, 0.55, 0.75].
In the text

Fig. B.2. Upper panel: Example of a cut in a different region of the BEL that is used to highlight the presence of a binary in the equal-mass case. Lower panel: qualitative example of the elements to be considered to highlight the presence of a binary in the case of an equal-mass binary with no spiral pattern. The red elements are the ones that are used in the method proposed for unequal-mass binaries, while the green ones are the ones that will encode the binary signature for the equal-mass binary.

In the text

	Fig. B.3. Example of how a cut in a different region of the BEL (those highlighted in figure B.2) can highlight the presence of a binary with the red and blue light curves separated by a half-period lag.
In the text