EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 709 (May 2026) A&A, 709 (2026) A246 Full HTML
Open Access
Issue
A&A
Volume 709, May 2026
Article Number A246
Number of page(s) 12
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202556979
Published online 25 May 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

In the hierarchical model of galaxy formation, smaller galaxies merge to form the larger galaxies (White & Rees 1978). As nearly all massive galaxies are expected to host a massive black hole (MBH) at their centre (Kormendy & Gebhardt 2001), massive black hole pairs (MBHPs) and, at even smaller separations, massive black hole binaries (MBHBs) are predicted to form in the aftermath of the merger (Begelman et al. 1980).

In our current view, two MBHs, still unbound from one another and residing at the centres of the merging galaxies, are expected to drift towards the centre of the newly formed galaxy because of dynamical friction. The dynamical friction efficiency at large MBH separations (0.001 − 10 kpc) is still poorly constrained (see Dotti et al. 2012), due to both the possible interaction between the MBHs and the galaxy substructures such as spiral arms, bars, and stellar or gaseous clumps (e.g. Fiacconi et al. 2013; del Valle et al. 2015; Souza Lima et al. 2020; Bortolas et al. 2020, 2022), and tidal effects. More specifically, tides can significantly increase the pairing timescale by decreasing the mass of the progenitor galaxy cores1, (e.g. Governato et al. 1994; Pfister et al. 2019; Varisco et al. 2024).

Past this stage, the evolution is expected to be driven by three-body scatterings with stars (Begelman et al. 1980; Mikkola & Valtonen 1992) and/or the interaction with the gas surrounding the two MBHs (Armitage & Natarajan 2002; Haiman et al. 2009). If the above-mentioned processes manage to efficiently remove orbital energy and angular momentum from the two MBHs, the resulting binary, i.e. the system where the two MBHs are gravitationally bound to each other, will coalesce due to gravitational wave emission (see Thorne & Braginskii 1976), thus becoming one of the most promising source populations to be detected by space-borne interferometers, for example the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2023), and pulsar timing arrays (see Verbiest et al. 2016). Recently, evidence for a gravitational wave background possibly generated by inspiralling MBHBs has been found in pulsar timing array data (see Agazie et al. 2023; Reardon et al. 2023; Xu et al. 2023; EPTA Collaboration and InPTA Collaboration 2024; Miles et al. 2025). Even though the background is likely generated by MBHBs, the uncertainties in stellar- and gas-driven evolution result in a significant scatter in the predicted gravitational wave signal (e.g. Varisco et al. 2021; Franchini et al. 2023, and references therein). Stringent EM constraints on MBHPs and MBHBs are needed for a comparison between theory and observation.

To date, there are only ∼100 confirmed observations of MBHPs resolved as dual AGNs (see De Rosa et al. 2019 for a review, Zhang et al. 2021 for a more recent collection of ∼10 kpc MBHPs, Hwang et al. 2020; Mannucci et al. 2022; Scialpi et al. 2024; Chen et al. 2023; Wu et al. 2024; Chen et al. 2025; Schwartzman et al. 2024 for recent searches of kiloparsec-scale MBHPs and Trindade Falcão et al. 2024 for a candidate with a ∼100 pc separation). A secure identification of spatially unresolved MBHPs and binaries is challenging. Current searches for unresolved signatures (see D’Orazio & Charisi 2023, for a review) have focused on close MBHBs where quasi-periodic modulations of the continuum (Valtonen et al. 2008; Ackermann et al. 2015; Graham et al. 2015; D’Orazio et al. 2015; Li et al. 2016; Charisi et al. 2016; Sandrinelli et al. 2016, 2018; Severgnini et al. 2018; Li et al. 2019; Liu et al. 2019; Hu et al. 2020; Chen et al. 2020; Zhu & Thrane 2020; Cocchiararo et al. 2024; Rigamonti et al. 2025a; Bertassi et al. 2025a) and of the broad emission line (BEL) profile (Shen & Loeb 2010; Tsalmantza et al. 2011; Eracleous et al. 2012; Decarli et al. 2013; Runnoe et al. 2017, 2025; Sottocorno et al. 2026; Rigamonti et al. 2025b; Bertassi et al. 2025b) are expected to occur on timescales of months to centuries.

At large enough separations, the search for photometric and spectroscopic features becomes either impossible (as the underlying assumptions of the models break) or unfeasible (as the required length of the observational campaigns to observe such modulations becomes too long). One exception to this is a test originally proposed by Gaskell (1988) and recently quantitatively detailed by Dotti et al. (2023), which assumes that the two MBHs are far enough from each other for their emission (both in the continuum and in the BELs) to be uncorrelated on observational timescales. In this case, different spectral regions of BELs react at two independent continua. This test identifies pairs and binaries with orbital periods of decades to centuries using reverberation-mapping campaigns lasting ≲1 yr (Dotti et al. 2023), but fails when the separation is so large (e.g. ∼0.5 pc for a binary of ∼106 M⊙) that the two independent broad-line regions have a too-small relative shift, due to the small orbital velocity of the binary (≲100 km/s, see Dotti et al. 2023).

In this work, we present a new photometric test for couples of active MBHs (either MBHBs or MBHPs) too widely separated to be identified through continuum periodicity (e.g., with separations ≳0.001 pc for a binary of ∼106 M⊙), but too close to be spatially resolved. For example, assuming an angular resolution of ≈ 0.65 arcsec in the r band for the Vera Rubin observatory (Ivezić et al. 2019), the maximum projected separation between the two MBHs for the test to succeed is ≈5 kpc at z ≈ 1 and ≈1.2 kpc at z ≈ 0.1. In the following, we refer to such systems as unresolved black hole duos (UBHDs).

The test leverages the intrinsic variability of AGN light curves (see Paolillo & Papadakis 2025, for a review) well modelled with a damped random walk (DRW, see Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010), a stochastic process with red-noise2, i.e. a frequency-dependent noise with power increasing towards the lower frequencies and becoming white noise, i.e. frequency independent at the lowest frequencies (see Section 2).

At the separations considered in this study, the short-timescale variability properties of the two MBHs were uncorrelated (as in the spectroscopic test proposed by Gaskell 1988 and Dotti et al. 2023). With simulations of UBHDs consisting of two DRW components, we tested whether these systems can be distinguished from the typical AGN behaviour described by a single DRW. For this, we employed a Bayesian model comparison.

By allowing the identification of unresolved UBHDs (associated with long residence timescales, which are expected to be abundant, see e.g. Haiman et al. 2009), the test has the potential to reveal large samples of UBHDs compared to the close MBHBs already searched for in time-domain surveys. These surveys are current, for example the Zwicky Transient Facility (ZTF; Bellm et al. 2018) and the Catalina Real-Time Transient Survey (CTRS; Djorgovski et al. 2011), and upcoming, for example the Vera Rubin Observatory’s ten-year Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) and the Roman Space Telescope’s High Latitude Time Domain Survey (HLTDS; Haiman et al. 2023).

In Section 2, we describe the properties of the DRW model. We also describe how mock light curves are generated, and discuss the methods with which the parameters are constrained and the model comparison is performed. In Section 3, we discuss the fraction of false positives expected, we quantify how well the injected parameters are retrieved, and identify the regions in the parameter space in which the test correctly identifies an UBHD. Our main conclusions are presented in Section 4. Appendix A discusses the properties of light-curve periodograms, and highlights the limitations that prompted the design of the analysis discussed in the main text. A less statistically solid test searching for UBHDs by fitting the power spectral density (PSD) obtained through the Lomb-Scargle periodogram (see VanderPlas 2018, and references therein) shape directly is discussed in Appendix B.

2. Methodology

For this work, we focused on UBHDs, i.e. MBHs separated enough to have independent emission variability, so that the PSD of the resulting time series is given by the sum of the power spectra of the two processes. Furthermore, we assumed that the only source of variability for each MBH is the intrinsic one produced by a DRW. We generated mock DRW light curves for each component MBH using celerite3 (see Foreman-Mackey et al. 2017) and created mock UBHD light curves by summing the realisations of the flux of two different MBHs (as detailed in Section 2.2). First, we considered evenly sampled time series, and then generalised the work to unevenly sampled time series, simulating the expected cadence of LSST. Next, we estimated the parameters and the evidence of a given model, either single MBH or UBHD, sampling the likelihood with a nested sampling algorithm (see Skilling 2006) implemented in RayNest4 (Veitch et al. 2024). Finally, we selected the best model by comparing the statistical evidence for each model.

2.1. Noise models

The behaviour of a quantity evolving under a DRW is described by the solution to the following stochastic differential equation (see Kelly et al. 2009):

d X ( t ) = − 1 τ X ( t ) d t + σ ̂ dt ϵ ( t ) + b d t with τ , σ ̂ , b > 0 . Mathematical equation: $$ \begin{aligned} dX(t) = -\frac{1}{\tau } X(t) \ dt + \hat{\sigma } \ \sqrt{dt} \ \epsilon (t) + \ b \ dt \ \ \mathrm {with} \ \tau , \ \hat{\sigma }, \ b \ >0. \end{aligned} $$(1)

Here X(t) is the process that is being observed as a function of time, τ is the damping timescale of the process (i.e. the time that is needed for the time series to become roughly uncorrelated), ϵ(t) is a white noise with zero mean and unit variance, and σ ̂ Mathematical equation: $ \hat{\sigma} $ describes the light curve variations at small timescales and is related to the variance of the process σ = σ ̂ τ / 2 Mathematical equation: $ \sigma=\hat{\sigma} \tau /2 $. Finally, b is related to the mean of the process given by bτ. This noise model has a PSD of the form:

P sing ( f ) = 2 σ ̂ 2 τ 2 1 + ( 2 π τ f ) 2 , Mathematical equation: $$ \begin{aligned} P_{\rm {sing}}(f) = \frac{2 \hat{\sigma }^2 \tau ^2}{1+(2 \pi \tau f)^2}, \end{aligned} $$(2)

while its covariance matrix can be written as:

C i , j sing = C i , j DRW = 1 2 σ ̂ 2 τ exp ( − Δ t i , j τ ) Mathematical equation: $$ \begin{aligned} C^\mathrm{{sing}}_{i,j}=C^\mathrm{{DRW}}_{i,j}=\frac{1}{2}\hat{\sigma }^2 \tau \exp \left(-\frac{\Delta t_{i,j}}{\tau }\right) \end{aligned} $$(3)

where Δti, j = |ti − tj| and i, j = 1…N, represent two arbitrary observations out of the total number of N pointings.

When considering an UBHD, assuming the light curve from one MBH to be independent of the other, the PSD can be written as the sum of the two PSDs:

P UBHD ( f ) = 2 σ ̂ 1 2 τ 1 2 1 + ( 2 π τ 1 f ) 2 + 2 σ ̂ 2 2 τ 2 2 1 + ( 2 π τ 2 f ) 2 Mathematical equation: $$ \begin{aligned} P_{\rm {UBHD}}(f) = \frac{2\hat{\sigma }_1^2 \tau _1^2}{1+(2 \pi \tau _1 f)^2} + \frac{2\hat{\sigma }_2^2 \tau _2^2}{1+(2 \pi \tau _2 f)^2} \end{aligned} $$(4)

where τ 1 , σ ̂ 1 Mathematical equation: $ \tau_1, \ \hat{\sigma}_1 $ and τ 2 , σ ̂ 2 Mathematical equation: $ \tau_2, \ \hat{\sigma}_2 $ are the damping timescales and short-term variability magnitudes of the two processes. The covariance function for this model can be written as the sum of two covariance functions:

C i , j UBHD = C i , j DRW , 1 + C i , j DRW , 2 = 1 2 σ ̂ 1 2 τ 1 exp ( − Δ t i , j τ 1 ) + 1 2 σ ̂ 2 2 τ 2 exp ( − Δ t i , j τ 2 ) . Mathematical equation: $$ \begin{aligned} C^\mathrm{{UBHD}}_{i,j }&=C^\mathrm{{DRW, 1}}_{i,j }+C^\mathrm{{DRW,2}}_{i,j}\nonumber \\&=\frac{1}{2}\hat{\sigma }_1^2 \tau _1 \exp \left(-\frac{\Delta t_{i,j}}{\tau _1}\right)+\frac{1}{2}\hat{\sigma }_2^2 \tau _2 \exp \left(-\frac{\Delta t_{i,j}}{\tau _2}\right). \end{aligned} $$(5)

This is consistent with the assumption of two independent, stationary, zero-mean processes being added, as the absence of cross-correlations between the two processes allows both the covariance functions and PSDs to be additive.

Several studies have been carried out to find a relation between the parameters of the stochastic equation and the physical parameters of the MBH (see Kelly et al. 2009; MacLeod et al. 2010; Burke et al. 2021; Arévalo et al. 2023, 2024; Su et al. 2024; Benati Gonçalves et al. 2025 for some examples), but there is currently no consensus regarding these relationships. The methodology followed in this work does not require any strong assumptions about the UBHD physical parameters. Since our objective is to explore the parameter space as broadly and agnostically as possible, we do not assume any relation between the parameters of the noise model and those of the MBHs.

2.2. Time series generation

Each light curve was generated by sampling the flux of the MBH using celerite. Once the kernel of the Gaussian process (GP) was chosen, the light curve could be sampled with an arbitrary sampling pattern and baseline. For this work, the injected parameters were drawn from log-uniform distributions for the damping timescales and uniform distributions for the variances of the two processes, as follows:

– Single MBH case:

10 d < τ < 500 d , 0 < σ ̂ 2 < 1 Mathematical equation: $$ \begin{aligned} 10 \, \mathrm{d}&< \tau < 500 \, \mathrm{d} , \\ 0&< \hat{\sigma }^2 < 1 \end{aligned} $$

– UBHD case:

10 d < τ 1 < 500 d , 10 d < τ 2 < τ 1 , 0 < σ ̂ 1 2 < 1 , 0 < σ ̂ 2 2 < 1 . Mathematical equation: $$ \begin{aligned} 10 \, \mathrm{d}&< \tau _1 < 500 \, \mathrm{d} , \\ 10 \, \mathrm{d}&< \tau _2 < \tau _1, \\ 0&< \hat{\sigma }^2_1 < 1, \\ 0&< \hat{\sigma }^2_2 < 1. \end{aligned} $$

For our work, we sampled the time series for 10 years to resemble the LSST nominal survey. As noted in Kozłowski (2017), to obtain good constraints on the flat part of the PSD, the process should be observed for at least 10 τ. With the assumed observational baseline, the flat part of the PSD can be observed particularly well for MBHs with a mass up to ∼108 M⊙ since they are expected to have a damping time of the order of 200 days (Burke et al. 2021). We initially sampled the light curves evenly in time every 3 days, and then generalised the test to more realistic light curves by omitting three to four months of observations each year to incorporate the inevitable gaps of ground-based observations. We assumed a relative flux uncertainty of 0.01. Such an uncertainty corresponds to a magnitude of ∼21 in the r band of a single LSST visit, based on the expected photometric precision Ivezić et al. (2019). In the case of an evenly sampled light curve, we assumed a time step of Δt = 3 day, roughly corresponding to the cadence of the Wide Fast Deep survey mode of LSST, which will cover about 95% of the survey time. Examples of light curves generated in such a way, along with their respective PSDs, are shown in Figure 1. Here it is possible to see how a change in the parameters changes the observed light curve and its power spectrum. More specifically, a change in the parameter σ ̂ Mathematical equation: $ \hat{\sigma} $ results in a uniform vertical shift of the PSD, with higher values of σ ̂ Mathematical equation: $ \hat{\sigma} $ leading to increased power at all frequencies, as shown in (b) panel of Figure 1. A change in the damping timescale τ causes a change in the break frequencies with a higher τ moving the flat part of the spectrum towards lower frequencies, as shown in panel (d) of Figure 1. The light curves for the UBHD case are constructed by summing two light curves generated as in the single MBH scenario with equal weight.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Examples of light curves generated using celerite varying the process parameters: damping timescale τ (panel a), and the variability amplitude σ ̂ Mathematical equation: $ \hat{\sigma} $ (panel c). In the right panels, the retrieved periodograms (computed as discussed in Appendix A) as well as the theoretical expected PSD (dashed lines) are shown. The time series ( F ̂ Mathematical equation: $ \hat{F} $) are centred at zero and have arbitrary units. The parameter σ and power (P) should be considered in terms of the same dimensional scale as F ̂ Mathematical equation: $ \hat{F} $, with P having units of the square of the units of F ̂ Mathematical equation: $ \hat{F} $.

2.3. Parameter estimation and model selection

For this work, we combined nested sampling and GPs to estimate the posterior distributions of the parameters and to compare the two competing models (for works that use a similar approach, see Covino et al. 2020; Zhu & Thrane 2020; Covino et al. 2022; Witt et al. 2022). While GPs often employ a maximum likelihood approach to optimise kernel hyperparameters for computational convenience, we took a fully Bayesian approach by using nested sampling to marginalise over these parameters. More specifically, nested sampling transforms the multidimensional integral for the evidence into a one-dimensional integral over the prior volume and explores nested contours of increasing likelihood by iteratively replacing low-likelihood samples with higher-likelihood ones. Unlike Markov chain Monte Carlo (MCMC) algorithms, nested sampling allows us to capture the full posterior distribution and compute the model evidence, enabling a statistically robust comparison between models (Skilling 2006; Ashton et al. 2022; Buchner 2023). We note that, due to the parallelism between Gaussian process kernels and PSD, UBHDs can be searched for through the fitting of the periodograms of evenly spaced light curves. However, for uneven data samplings, a consistent Bayesian analysis is not possible (as detailed in Appendix A). An analysis directly performed on the PSDs is presented in Appendix B together with its shortcomings.

The fitting of the light curve through GPs is performed using the GP regressor celerite, using a kernel in the form

k single i , j = a e − c Δ t i , j . Mathematical equation: $$ \begin{aligned} k_{\mathrm{single} \ i, j} = a e^{-c \ \Delta t_{i,j}}. \end{aligned} $$(6)

Equation (6) becomes equal to the kernel in Equation (3) by setting a = 0.5τσ2 and c = 1/τ. In the context of celerite, a kernel such as the one in Equation (6) is known as RealTerm. The UBHD model can then be written as the sum of two RealTerm kernels:

k UBHD ( Δ t ) = a 1 e − c 1 Δ t + a 2 e − c 2 Δ t . Mathematical equation: $$ \begin{aligned} k_{\mathrm{UBHD} }(\Delta t) = a_{1} e^{-c_{1} \Delta t} + a_{2} e^{-c_{2} \Delta t}. \end{aligned} $$(7)

We note that we chose to use celerite due to its computational efficiency in evaluating the likelihood for kernels such as those presented in Equation (6). Specifically, the computational cost of celerite scales linearly with the number of data points, in contrast to the cubic scaling of standard GP methods.

The nested sampling is performed using RayNest. The priors are a log-uniform distribution in τ in between 10 and 500 d, and log-uniform in σ ̂ 2 Mathematical equation: $ \hat{\sigma}^2 $ in between 10−6 and 106. We emphasise that the prior distribution on the variance σ ̂ 2 Mathematical equation: $ \hat{\sigma}^2 $ is much wider than the distribution used to generate the light curves. This choice was made to ensure that the results and posteriors are not biased by an artificial truncation of the prior boundaries, which can affect the resulting posterior and evidence.

For model comparison, we computed the Bayes factor, defined as the ratio of the evidence of the two models:

B pair , sing = ∫ p ( D | M pair , θ pair ) p ( θ pair | M pair ) d θ pair ∫ p ( D | M sing , θ sing ) p ( θ sing | M sing ) d θ sing · Mathematical equation: $$ \begin{aligned} B_{\rm {pair,sing}}=\frac{\int p(D|M_{\rm {pair}},\boldsymbol{\theta }_{\rm {pair}} \ )\ p (\boldsymbol{\theta }_{\rm {pair}}|M_{\rm {pair}}) \ d\boldsymbol{\theta }_{\rm {pair}}}{\int p(D|M_{\rm {sing}},\boldsymbol{\theta }_{\rm {sing}} \ )\ p (\boldsymbol{\theta }_{\rm {sing}}|M_{\rm {sing}}) \ d\boldsymbol{\theta }_{\rm {sing}}}\cdot \end{aligned} $$(8)

Here p(D|M, θ, I) is the likelihood function and p(θ|M, I) is the prior; D is the data, i.e. the observed light curve (ti, Fi), with i, 1, 2, …N the number of data points and Fi being the flux observed at time ti; M is the considered model; and θ are the parameters for this particular model. In this work, we assume that the two models have the same prior probability. In the particular case of GPs, the likelihood of the data is given by (see Rasmussen & Williams 2006):

p ( D ∣ M , θ , I ) = 1 ( 2 π ) N 2 | Σ | 1 2 exp ( − 1 2 ( D − μ ) T Σ − 1 ( D − μ ) ) Mathematical equation: $$ \begin{aligned} p(D \mid M, \boldsymbol{\theta }, I) = \frac{1}{(2\pi )^{\frac{N}{2}} |\boldsymbol{\Sigma }|^{\frac{1}{2}}} \exp \left( -\frac{1}{2} (D - \boldsymbol{\mu })^{T} \boldsymbol{\Sigma }^{-1} (D - \boldsymbol{\mu }) \right) \end{aligned} $$(9)

where μ is the mean function of the GP, N is the number of observations, and Σ and |Σ| are the covariance matrix of the GP and its determinant. Here we consider the evidence for UBHD to be significant when Bpair, sing > 3, i.e. if it shows at least substantial evidence according to the Jeffreys scale (Jeffreys 1998).

3. Results

The aim of the test proposed in this work is to assess the presence of an UBHD directly from the analysis of the observed light curve. To achieve this, we analysed realistic time series using GPs and nested sampling.

In this section, we present the main findings of our study. First, we assessed the false positive rate of the proposed signature by analysing single MBH light curves with both models and computing the corresponding Bayes factor (Section 3.1) to infer the percentage of single MBH light curves wrongly identified as UBHDs. Next, we evaluated how accurately the nested sampling procedure recovers the injected parameters (Section 3.2), considering both the single MBH and UBHD scenarios, and using both evenly and unevenly sampled light curves. Finally, we explored which regions of the model parameter space allow a reliable identification of UBHDs. To do this, we randomly sampled parameters and evaluated the Bayes factor to determine the conditions under which the test is most effective (Section 3.3).

3.1. False positives

To quantify the robustness of the proposed signature, it is important to quantify how unique it is versus how often it can be confused with the typical DRW variability of a single MBH. For this, we assessed the fraction of false positives produced by the test by simulating 1000 evenly sampled and 1000 unevenly sampled single MBH light curves and calculating the Bayes factor for the UBHD versus the single MBH model, as described in Section 2.3. We then computed the fraction of light curves that our test identifies as UBHDs, i.e. the fraction of light curves for which the Bayes factor is greater than 3. We find the 0.2% and 0.59% of false positives in the evenly and unevenly sampled scenarios, respectively. The percentage of false positives remains smaller than 1%, increasing when the sampling becomes uneven, even if the length of the observational campaign remains the same. We emphasise that the overall false detection rate is very low (regardless of light curve quality).

3.2. Parameter recovery

Before searching for the best parameter space in which the analysis of the photometric light curve can show evidence of the presence of an UBHD, we tested how well our algorithm can recover the parameters of the injected light curve for each model. This comparison was not a simple consistency check between model and data realisations, as the processes we looked at are intrinsically stochastic, i.e. individual realisations would differ even for the same set of parameters. To assess the goodness of the parameter estimation, we computed the relative error between the injected θinj and recovered θrec parameters as:

δ θ = θ rec − θ inj θ inj × 100 % Mathematical equation: $$ \begin{aligned} \delta \theta =\frac{\theta _{\mathrm{rec} }- \theta _{\mathrm{inj} }}{{\theta _{\mathrm{inj} }}}\times 100 \% \end{aligned} $$(10)

and the relative uncertainty range (σθ) defined as:

σ θ = q θ , 0.84 − q θ , 0.16 θ inj × 100 % Mathematical equation: $$ \begin{aligned} \sigma _{\theta }=\frac{q_{\theta ,0.84}-q_{\theta ,0.16}}{\theta _{\rm {inj}}} \times 100\% \end{aligned} $$(11)

where θinj is one of the parameters with which the light curve has been generated, θrec is the median value of the posterior distribution for that particular parameter, qθ, 0.16 and qθ, 0.84 are respectively the 0.16 and 0.84 quantiles of the posterior distribution.

In Figure 2, we show the distribution of the relative errors δθ for 1000 evenly sampled light curves for the single MBH case. The dashed vertical red line indicates the optimal scenario δθ = 0 while the dashed vertical blue line indicates the median value of δθ. We see that both parameters are recovered well as these distributions are scattered around δθ = 0 for both τ and σ ̂ Mathematical equation: $ \hat{\sigma} $. The values we found for the median relative error δθ and relative uncertainty range σθ for both the even and uneven scenarios are reported in Table 1.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Distribution of relative errors for the retrieved parameters in the single MBH evenly sampled case.

Table 1.

Parameter recovery for the single MBH model.

For the damping timescale, we see a tail at lower values, due to the fact that the length of the observational campaign is not long enough to efficiently sample the flat part of the PSD. In Figure 3, we do the same for the UBHD case. Here, the scatter from the true value is higher than that in the single DRW scenario, as can be seen from the values of the median relative error δθ and relative uncertainty range σθ reported in Table 2.

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Distribution of relative errors for the retrieved parameters in the UBHD evenly sampled case. Shown are the distributions for the cases in which the UBHD model is favoured (in green).

Table 2.

Parameter recovery for the UBHD model.

The green distributions in Figure 3 represent the distributions of δθ for the combination of parameters where the test correctly identifies an UBHD. The median values of the relative errors and the relative uncertainty range of such light curves are reported in Table 3. In this case, one can see that the scatter overall becomes smaller, but with τ1 still showing wide tails in the distributions, similarly to the single DRW case.

Table 3.

Parameter recovery for the correctly identified UBHD light curves.

In Table 4, we report the percentage of light curves for which all the parameters are retrieved with a relative error within the 10%, 20%, and 50% respectively for the single MBH, UBHD, and correctly identified (green distributions in Figure 3) cases.

Table 4.

Percentages of light curves for which all the parameters are retrieved with a relative error within the 10%, 20% and 50%.

3.3. True positives

Finally, we searched for the best region of the parameter space for which the test can efficiently identify UBHDs. More specifically, we sampled 1500 light curves (for both evenly and unevenly sampled light curves)varying the ratio of the two damping timescales qτ = τ2/τ1 with τ2 < τ1 and the ratio of the amplitudes of the two processes q σ = σ 2 / σ 1 = ( σ ̂ 2 τ 2 ) / ( σ ̂ 1 τ 1 ) Mathematical equation: $ q_{\sigma}=\sigma_2/\sigma_1=(\hat{\sigma}_2\sqrt{\tau_2})/(\hat{\sigma}_1\sqrt{\tau_1}) $ with σ2 < σ1. For each pair of qτ and qσ, we computed the Bayes factor in order to highlight where the test works the best, i.e. where the Bayes factor is greater than 3.

In Figure 4, we show the result in the case of unevenly sampled light curves observed for 10 years with a cadence of 3 d. The points are colour-coded by the value of log10B, and the set of parameters for which the test yields a Bayes factor greater than 3 is represented by stars. The result does not significantly change when considering evenly sampled light curves. More specifically, in the evenly sampled scenario, 4.65% of all the considered light curves are correctly identified as UBHDs, while in the unevenly sampled scenario, this fraction reduces to 3.77%.

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Map of Bayes factors for 1500 unevenly sampled realisations of UBHD light curves with a ten-year baseline. The colour scale represents the base-10 logarithm of the Bayes factor, while the axes show the ratio of the model parameters between the two light curves. The red square, diamond, triangle and circle correspond to the light curves used to compute the periodograms shown in panels a,b,c, and d, respectively, in Figures 6 and B.1.

Figure 4 shows clearly that the test we propose works only in a small region of the parameter space where qτ ≲ 0.2 and qσ > 0.2. In principle, setting a lower Bayes factor threshold could increase the fraction of correctly retrieved UBHD light curves and expand the region of the parameter space where the test works. For example, thresholds on B of 2, 1, and 0.5 would increase the fraction of correctly retrieved UBHD light curves from 3.77% (when the threshold is set to B = 3) to 5.3%, 11.3% and 73.9% in the case of unevenly sampled light curves and from 4.6% (for B = 3) to 5.4%, 12.5% and 87.6% in the case of evenly sampled light curves. The ratio of the damping timescales below which the test works moves from qτ = 0.21 for B = 3 to qτ = 0.31, qτ ∼ 0.65 and qτ ∼ 0.87 for B = 2, 1, and 0.5, respectively. At the same time, lowering the threshold would increase the false positive fraction. For the uneven sampling scenario, the false positive fraction would increase from 0.59% (for B = 3) to 1%, 6% and 82%, while in the even sampling scenario the false positives would increase from 0.59% (for B = 3) to 1.07%, 6%, and 82% for Bayes factor thresholds of 2, 1, and 0.5, respectively. Therefore, lower thresholds would make the statistical evidence in favour of the UBHD scenario less robust. For this reason, we chose a more conservative limit for the threshold.

In Figure 5, we show 500 new simulations in order to provide a zoomed-in image of this region of parameter space for UBHD light curves. In particular, UBHDs are identified when the two MBHs have very different damping timescales and similar variability amplitudes. This is not surprising if one considers the shape of the PSD. In Figure 6, we show the PSDs for four cases, one of which is identified (in panel a) as an UBHD. More specifically, the search identifies an UBHD if, after the first plateau of the PSD up to the break frequency related to the longer damping timescale, a second plateau starts appearing up to the second break frequency related to the shorter damping timescale.

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Same kind of map as reported in Figure 4 but for 580 light curves in the region of the parameter space where the UBHDs can be correctly identified.

Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

PSD of double DRW kernels for different parameter ratios. The Bayes ratio resulting from the GP fitting is reported above each panel. In all the panels, the light curves are sampled every 3 d for 10 yr. The power is reported in arbitrary units. Case (a): Similar variability amplitudes but considerably different damping timescales. The test can correctly identify the UBHD (shown as the red square in Figure 4). Case (b): Similar amplitudes and similar damping timescales (shown as the red diamond in Figure 4). Case (c): Considerably different amplitudes but similar damping timescales (shown as the red triangle in Figure 4). Case (d): Considerably different amplitudes and damping timescales (shown as the red circle in Figure 4).

Since the damping timescale positively correlates with the MBH luminosity and with the mass of the MBH, we can conclude that pairs where the two MBHs have very different luminosities or very different masses (see MacLeod et al. 2010), i.e. a low luminosity ratio or a low mass ratio, are expected to be more easily identified. The ability of this analysis to identify UBHDs large-scale photometric surveys (which provide unevenly sampled data) makes it relevant for currently available photometric surveys, most of which feature gaps and uneven pointings.

We repeated the same procedure, changing the length of the observational baseline. More specifically, we studied the cases of light curves observed for 6 and 3 years, keeping the uneven sampling discussed in Section 2.2. These baselines are relevant for current time-domain surveys (e.g. ZTF), and for the early years of LSST. As can be seen in Figure 7, as the length of the observation decreases, the region where the test is able to correctly retrieve the presence of an UBHD shrinks to a region of higher qσ. Specifically, by repeating simulations in the region with qτ < 0.2, analogous to those used to generate Figure 5, we find that the value of qσ above which 90% of the correctly retrieved light curves lie moves from qσ = 0.21 for the 10 yr baseline to qσ = 0.36 and qσ = 0.40 for baselines of 6 and 3 yr respectively.

Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Similar to the map in Figures 4 but with different observational baselines. In the upper panel, the light curves are observed for 6 years, while in the lower panel, the light curves are observed for 3 years.

4. Discussion and conclusions

In this paper, we studied the possibility of detecting the presence of an UBHD using AGN light curves. In particular, we examined whether we can distinguish the presence of such a system from the typical variability of an AGN with a single MBH (typically described with a DRW model). For this signature, we focused on the optical variability, where we will soon have large samples of quasar light curves.

We generated both single MBH mock light curves and UBHD mock light curves by summing two single MBH time series with different cadences (either evenly sampled light curves with a three-day cadence or including observational gaps of a few months every year) and different lengths in terms of observation time (10, 6, or 3 yr). We then modelled the simulated time series employing GPs, using either a single DRW kernel or the sum of two uncorrelated DRW kernels. Finally, we checked whether the presence of a double DRW could be effectively observed through Bayesian model comparison. In particular, we used nested sampling to compute the evidence of the two models (single and double DRW) and calculated the Bayes factor, which allowed us to assess which model describes the simulated data better.

We studied how many false positives, i.e. light curves generated as single MBH, are better described by the double DRW model according to the Bayes factor (Section 3.1). To do so, we generated single MBH light curves with different cadences, and for each of the time series, we computed the Bayes ratio. We considered as a false positive a time series that presented a Bayes ratio higher than 3, i.e. the Bayes ratio indicates substantial evidence in favour of the double DRW model following the Jeffreys scale. More specifically, 1000 evenly sampled light curves with a 10 yr baseline show a false positive rate of 0.2%, while the same number of unevenly sampled light curves with the same baseline length result in 0.59%, demonstrating that the proposed signature has a low confusion rate with the typical AGN variability.

We assessed whether nested sampling successfully recovered the injected parameters, both for the single MBH scenario and for UBHDs (Section 3.2). To do so, we computed the relative per cent error δθ (i.e. the normalised difference between the injected and recovered parameter, see Equation (10)) and the relative per cent uncertainty range (i.e. the normalised uncertainty range between by the injected parameter, see Equation (11)). We showed that 51.16% of evenly sampled single MBH light curves, 50.90% of unevenly sampled single MBH light curves, 7.13% of evenly sampled UBHD light curves and 5.35% of unevenly sampled UBHD light curves with a 10 yr baseline show recovered parameters within 20% of the injected parameters.

We studied which region of the parameter space allows us to detect an UBHD (Section 3.3). In particular, we generated light curves with different ratios of the two damping timescales, and variances of the two DRW and constructed maps that show which combinations of these two ratios provide substantial evidence for the double DRW model. We found that our test is able to identify UBHDs that show very different damping timescales, qτ < 0.2, and similar variability amplitudes qσ > 0.2 for both evenly and unevenly sampled light curves observed for 10 years. For the correctly identified UBHDs, the fraction of light curves with δθ < 20% grows to 14.29% (8.31%) for the evenly (unevenly) sampled data. The low ratio in the damping timescales implies that only UBHDs with very different luminosities (and masses) will be correctly classified. Furthermore, we changed the length of the observation, maintaining an uneven sampling, and found that shortening the length of the observational baseline increases the qσ needed to correctly identify an UBHD.

In the Appendices, we explored another possible method to detect the presence of an unresolved UBHD by studying the shape of the periodogram applicable only for evenly sampled data, which is not readily available. In Appendix A, we showed how the periodogram can be computed and discussed the statistical properties of the classical periodogram. In Appendix B, we showed that by using the Whittle likelihood (see Whittle 1951) to describe the distribution of the retrieved periodogram with respect to the true PSD, one can assess whether a single or double DRW better describes the observed periodogram shape. This method works only in the case of evenly sampled time series. In the case of unevenly sampled time series, this approach cannot be used as the statistical properties are different from the classical properties and depend on the particular noise considered.

Once candidates are identified, confirming the UBHD nature of objects with intermediate values of Bpair, sing and discerning whether the two MBHs are bound in a binary or not will require follow-up observations. High-resolution imaging would have the chance to directly identify MBHPs (e.g. Hwang et al. 2020; Mannucci et al. 2022; Scialpi et al. 2024; Chen et al. 2023), while optical spectroscopy allows us to check for the presence of two sets of narrow emission lines (NELs), supporting the UBHD scenario (e.g. Comerford et al. 2009). If a single set of NELs is observed, the identification of BELs that are asymmetric or offset with respect to the host rest frame in optical/UV spectra would be an indication of the presence of a close binary (e.g. Shen & Loeb 2010; Tsalmantza et al. 2011; Eracleous et al. 2012). In this case, reverberation mapping studies could be used to further support binary interpretation, searching for uncorrelated variability of different parts of single BELs (Gaskell 1988; Dotti et al. 2023).

It should be noted that our analysis was able to identify systems composed of two MBHs whose hosts are not in interaction, for example, chance superpositions of two AGNs at different redshifts or within a galaxy cluster (e.g. Dotti & Ruszkowski 2010). This particular scenario can be tested by looking at the frequency shift between two sets of NELs and by looking at the local density of galaxies in the UBHD–candidate neighbourhood (e.g. Decarli et al. 2014). Another scenario that could be similar to the one explored in this work is the case of a strongly lensed quasar with angularly unresolved multiple images. However, such a scenario cannot produce false UBHDs as the light curves of the different images generated by gravitational lensing will have the same damping timescale, generated by the same quasar.

The only differences between the lensed light curves would be the amplitude, due to the different magnification of the images, and a time delay. To test this scenario, we simulated 2000 unresolved lensed quasar light curves by simulating realisations of a single DRW and duplicating them. The two duplicated light curves are then multiplied by different magnifications computed assuming an un-truncated Singular Isothermal Sphere (SIS) as the matter distribution for the lens, and sampling from a uniform distribution between 0 and 1 for the ratio between the source-lens angular separation and the lens Einstein radius. A delay sampled from a log-uniform distribution between a minute and two weeks is then introduced between the two light curves. The Bayes factor for all the simulated light curves is smaller than B = 1, thus favouring the single DRW model with respect to the UBHD model, as expected.

Finally, in this work, we examined our proposed signature in single-band AGN light curves. However, the upcoming LSST will deliver multi-band light curves, which may be advantageous to combine. Therefore, we plan to extend the analysis presented in this paper to light curves obtained in different wavelength bands (as explored in Covino et al. 2022, for the identification of periodicities in AGN light curves).

All of the above discussion and conclusions rest on the assumption that AGN intrinsic variability is accurately described by the DRW model. However, deviations from the DRW have been reported (e.g. Zu et al. 2013; Kasliwal et al. 2015; Wang & Shi 2019), and alternative models have been proposed, such as the damped harmonic oscillator (DHO) model (Yu et al. 2022). If AGN light curves follow a DHO process rather than a DRW, our results may be affected. In the over-damped DHO case (see Yu et al. 2022), the covariance function reduces to a sum of two exponential kernels, resembling our UBHD model. Nevertheless, the two models predict different parameter correlations, which could, in principle, be used to distinguish between them. We plan to investigate the impact of including the DHO model on our results in future work.

Acknowledgments

LB acknowledges ISCRA for awarding this project access to the LEONARDO supercomputer, owned by the EuroHPC Joint Undertaking, hosted by CINECA (Italy). LB wishes to thank the “Summer School for Astrostatistics in Crete” for providing training on the statistical methods adopted in this work. MC acknowledges support by the European Union (ERC, MMMonsters, 101117624). 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. MD acknowledge funding from MIUR under the grant PRIN 2017-MB8AEZ, financial support from ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU, and support by the Italian Ministry for Research and University (MUR) under Grant ‘Progetto Dipartimenti di Eccellenza 2023-2027’ (BiCoQ). We acknowledge a financial contribution from the Bando Ricerca Fondamentale INAF 2022 Large Grant, Dual and binary supermassive black holes in the multi-messenger era: from galaxy mergers to gravitational waves. FR acknowledges support from ELSA. ELSA Euclid Legacy Science Advanced analysis tools” (Grant Agreement no. 101135203) is funded by the European Union. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or Innovate UK. Neither the European Union nor the granting authority can be held responsible for them. UK participation is funded through the UK Horizon guarantee scheme under Innovate UK grant 10093177.

References

  1. Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 813, L41 [NASA ADS] [CrossRef] [Google Scholar]
  2. Agazie, G., Anumarlapudi, A., Archibald, A. M., et al. 2023, ApJ, 951, L8 [NASA ADS] [CrossRef] [Google Scholar]
  3. Amaro-Seoane, P., Andrews, J., Sedda, M. A., et al. 2023, Liv. Rev. Relat., 26, 2 [NASA ADS] [CrossRef] [Google Scholar]
  4. Arévalo, P., Lira, P., Sánchez-Sáez, P., et al. 2023, MNRAS, 526, 6078 [CrossRef] [Google Scholar]
  5. Arévalo, P., Churazov, E., Lira, P., et al. 2024, A&A, 684, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Armitage, P. J., & Natarajan, P. 2002, ApJ, 567, L9 [NASA ADS] [CrossRef] [Google Scholar]
  7. Ashton, G., Bernstein, N., Buchner, J., et al. 2022, Nat. Rev. Methods Primers, 2, 39 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barret, D., & Vaughan, S. 2012, ApJ, 746, 131 [NASA ADS] [CrossRef] [Google Scholar]
  9. Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature, 287, 307 [Google Scholar]
  10. Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2018, PASP, 131, 018002 [Google Scholar]
  11. Benati Gonçalves, H., Panda, S., Storchi Bergmann, T., Cackett, E. M., & Eracleous, M. 2025, ApJ, 988, 27 [Google Scholar]
  12. Bertassi, L., Charisi, M., Buscicchio, R., et al. 2025a, A&A, submitted [arXiv:2512.13688] [Google Scholar]
  13. Bertassi, L., Sottocorno, E., Rigamonti, F., et al. 2025b, A&A, 702, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Bloomfield, P. 2004, Fourier Analysis of Time Series: An Introduction (Wiley) [Google Scholar]
  15. Bortolas, E., Capelo, P. R., Zana, T., et al. 2020, MNRAS, 498, 3601 [Google Scholar]
  16. Bortolas, E., Bonetti, M., Dotti, M., et al. 2022, MNRAS, 512, 3365 [NASA ADS] [CrossRef] [Google Scholar]
  17. Buchner, J. 2023, Stat. Surv., 17, 169 [NASA ADS] [CrossRef] [Google Scholar]
  18. Burke, C. J., Shen, Y., Blaes, O., et al. 2021, Science, 373, 789 [NASA ADS] [CrossRef] [Google Scholar]
  19. Callegari, S., Mayer, L., Kazantzidis, S., et al. 2009, ApJ, 696, L89 [NASA ADS] [CrossRef] [Google Scholar]
  20. Charisi, M., Bartos, I., Haiman, Z., et al. 2016, MNRAS, 463, 2145 [Google Scholar]
  21. Chen, Y.-C., Liu, X., Liao, W.-T., et al. 2020, MNRAS, 499, 2245 [Google Scholar]
  22. Chen, Y.-C., Liu, X., Lazio, J., et al. 2023, ApJ, 958, 29 [NASA ADS] [CrossRef] [Google Scholar]
  23. Chen, Y.-C., Gross, A. C., Liu, X., et al. 2025, ApJ, 988, 126 [Google Scholar]
  24. Cocchiararo, F., Franchini, A., Lupi, A., & Sesana, A. 2024, A&A, 691, A250 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Comerford, J. M., Gerke, B. F., Newman, J. A., et al. 2009, ApJ, 698, 956 [Google Scholar]
  26. Covino, S., Landoni, M., Sandrinelli, A., & Treves, A. 2020, ApJ, 895, 122 [NASA ADS] [CrossRef] [Google Scholar]
  27. Covino, S., Tobar, F., & Treves, A. 2022, MNRAS, 513, 2841 [Google Scholar]
  28. De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, New Astron. Rev., 86, 101525 [Google Scholar]
  29. Decarli, R., Dotti, M., Fumagalli, M., et al. 2013, MNRAS, 433, 1492 [Google Scholar]
  30. Decarli, R., Dotti, M., Mazzucchelli, C., Montuori, C., & Volonteri, M. 2014, MNRAS, 445, 1558 [NASA ADS] [CrossRef] [Google Scholar]
  31. del Valle, L., Escala, A., Maureira-Fredes, C., et al. 2015, ApJ, 811, 59 [NASA ADS] [CrossRef] [Google Scholar]
  32. Djorgovski, S. G., Drake, A. J., Mahabal, A. A., et al. 2011, ArXiv e-prints [arXiv:1102.5004] [Google Scholar]
  33. D’Orazio, D. J., & Charisi, M. 2023, ArXiv e-prints [arXiv:2310.16896] [Google Scholar]
  34. D’Orazio, D. J., Haiman, Z., & Schiminovich, D. 2015, Nature, 525, 351 [Google Scholar]
  35. Dotti, M., & Ruszkowski, M. 2010, ApJ, 713, L37 [Google Scholar]
  36. Dotti, M., Sesana, A., & Decarli, R. 2012, Adv. Astron., 2012, 940568 [CrossRef] [Google Scholar]
  37. Dotti, M., Rigamonti, F., Rinaldi, S., et al. 2023, A&A, 680, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. EPTA Collaboration and InPTA Collaboration, (Antoniadis, J., et al.) 2024, A&A, 690, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Eracleous, M., Boroson, T. A., Halpern, J. P., & Liu, J. 2012, ApJS, 201, 23 [Google Scholar]
  40. Fiacconi, D., Mayer, L., Roškar, R., & Colpi, M. 2013, ApJ, 777, L14 [NASA ADS] [CrossRef] [Google Scholar]
  41. Foreman-Mackey, D., Agol, E., Angus, R., & Ambikasaran, S. 2017, AJ, 154, 220 [NASA ADS] [CrossRef] [Google Scholar]
  42. Franchini, A., Lupi, A., Sesana, A., & Haiman, Z. 2023, MNRAS, 522, 1569 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gaskell, C. M. 1988, in Active Galactic Nuclei, eds. H. R. Miller, & P. J. Wiita, 307, 61 [Google Scholar]
  44. Governato, F., Colpi, M., & Maraschi, L. 1994, MNRAS, 271, 317 [Google Scholar]
  45. Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74 [Google Scholar]
  46. Haiman, Z., Kocsis, B., & Menou, K. 2009, ApJ, 700, 1952 [CrossRef] [Google Scholar]
  47. Haiman, Z., Xin, C., Bogdanović, T., et al. 2023, ArXiv e-prints [arXiv:2306.14990] [Google Scholar]
  48. Hu, B. X., D’Orazio, D. J., Haiman, Z., et al. 2020, MNRAS, 495, 4061 [Google Scholar]
  49. Hwang, H.-C., Shen, Y., Zakamska, N., & Liu, X. 2020, ApJ, 888, 73 [NASA ADS] [CrossRef] [Google Scholar]
  50. Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111 [Google Scholar]
  51. Jeffreys, H. 1998, The Theory of Probability (Oxford: OUP), Oxford Classic Texts Phys. Sci. [Google Scholar]
  52. Jenkins, G., & Watts, D. 1968, Spectral Analysis and Its Applications (Holden-Day) [Google Scholar]
  53. Kasliwal, V. P., Vogeley, M. S., & Richards, G. T. 2015, MNRAS, 451, 4328 [CrossRef] [Google Scholar]
  54. Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895 [Google Scholar]
  55. Kormendy, J., & Gebhardt, K. 2001, AIP Conf. Ser., 586, 363 [Google Scholar]
  56. Kozłowski, S. 2017, A&A, 597, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Kozłowski, S., Kochanek, C. S., Udalski, A., et al. 2010, ApJ, 708, 927 [CrossRef] [Google Scholar]
  58. Li, Y.-R., Wang, J.-M., Ho, L. C., et al. 2016, ApJ, 822, 4 [NASA ADS] [CrossRef] [Google Scholar]
  59. Li, Y.-R., Wang, J.-M., Zhang, Z.-X., et al. 2019, ApJS, 241, 33 [Google Scholar]
  60. Liu, T., Gezari, S., Ayers, M., et al. 2019, ApJ, 884, 36 [Google Scholar]
  61. Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
  62. MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, ApJ, 721, 1014 [Google Scholar]
  63. Mannucci, F., Pancino, E., Belfiore, F., et al. 2022, Nat. Astron., 6, 1185 [NASA ADS] [CrossRef] [Google Scholar]
  64. Mikkola, S., & Valtonen, M. J. 1992, MNRAS, 259, 115 [NASA ADS] [CrossRef] [Google Scholar]
  65. Miles, M. T., Shannon, R. M., Reardon, D. J., et al. 2025, MNRAS, 536, 1489 [Google Scholar]
  66. Paolillo, M., & Papadakis, I. 2025, Nuovo Cimento Riv. Ser., 48, 537 [Google Scholar]
  67. Papadakis, I. E., & McHardy, I. M. 1995, MNRAS, 273, 923 [NASA ADS] [Google Scholar]
  68. Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101 [Google Scholar]
  69. Priestley, M. 1981, Spectral Analysis and Time Series (Academic Press) [Google Scholar]
  70. Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning [Google Scholar]
  71. Reardon, D. J., Zic, A., Shannon, R. M., et al. 2023, ApJ, 951, L6 [NASA ADS] [CrossRef] [Google Scholar]
  72. Rigamonti, F., Bertassi, L., Buscicchio, R., et al. 2025a, A&A, 702, A242 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Rigamonti, F., Severgnini, P., Sottocorno, E., et al. 2025b, A&A, 693, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Runnoe, J. C., Eracleous, M., Pennell, A., et al. 2017, MNRAS, 468, 1683 [Google Scholar]
  75. Runnoe, J. C., Eracleous, M., Bogdanović, T., Halpern, J. P., & Sigurðsson, S. 2025, ApJ, 984, 17 [Google Scholar]
  76. Sandrinelli, A., Covino, S., Dotti, M., & Treves, A. 2016, AJ, 151, 54 [NASA ADS] [CrossRef] [Google Scholar]
  77. Sandrinelli, A., Covino, S., Treves, A., et al. 2018, A&A, 615, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
  79. Schwartzman, E., Clarke, T. E., Nyland, K., et al. 2024, ApJ, 961, 233 [NASA ADS] [Google Scholar]
  80. Scialpi, M., Mannucci, F., Marconcini, C., et al. 2024, A&A, 690, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Severgnini, P., Cicone, C., Della Ceca, R., et al. 2018, MNRAS, 479, 3804 [Google Scholar]
  82. Shen, Y., & Loeb, A. 2010, ApJ, 725, 249 [NASA ADS] [CrossRef] [Google Scholar]
  83. Skilling, J. 2006, Bayesian Anal., 1, 833 [Google Scholar]
  84. Sottocorno, E., Ogborn, M., Bertassi, L., et al. 2026, A&A, 708, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Souza Lima, R., Mayer, L., Capelo, P. R., Bortolas, E., & Quinn, T. R. 2020, ApJ, 899, 126 [Google Scholar]
  86. Su, Z.-B., Cai, Z.-Y., Sun, M., et al. 2024, ApJ, 969, 78 [Google Scholar]
  87. Thorne, K. S., & Braginskii, V. B. 1976, ApJ, 204, L1 [NASA ADS] [CrossRef] [Google Scholar]
  88. Trindade Falcão, A., Turner, T. J., Kraemer, S. B., et al. 2024, ApJ, 972, 185 [Google Scholar]
  89. Tsalmantza, P., Decarli, R., Dotti, M., & Hogg, D. W. 2011, ApJ, 738, 20 [Google Scholar]
  90. Uttley, P., McHardy, I. M., & Papadakis, I. E. 2002, MNRAS, 332, 231 [Google Scholar]
  91. Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Nature, 452, 851 [Google Scholar]
  92. van der Klis, M. 1989, in Fourier Techniques in X-Ray Timing, eds. H. Ögelman, & E. P. J. van den Heuvel (Dordrecht, Netherlands: Springer), 27 [Google Scholar]
  93. VanderPlas, J. T. 2018, ApJS, 236, 16 [Google Scholar]
  94. Varisco, L., Bortolas, E., Dotti, M., & Sesana, A. 2021, MNRAS, 508, 1533 [NASA ADS] [CrossRef] [Google Scholar]
  95. Varisco, L., Dotti, M., Bonetti, M., Bortolas, E., & Lupi, A. 2024, A&A, 689, A279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Vaughan, S., Edelson, R., Warwick, R. S., & Uttley, P. 2003, MNRAS, 345, 1271 [Google Scholar]
  97. Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145 [Google Scholar]
  98. Veitch, J., Del Pozzo, W., Lyttle, A., et al. 2024, https://doi.org/10.5281/zenodo.12801702 [Google Scholar]
  99. Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267 [Google Scholar]
  100. Vio, R., Andreani, P., & Biggs, A. 2010, A&A, 519, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Wang, H., & Shi, Y. 2019, Ap&SS, 364, 27 [Google Scholar]
  102. White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341 [Google Scholar]
  103. Whittle, P. 1951, Hypothesis Testing in Time Series Analysis (Almqvist& Wiksells boktr.) [Google Scholar]
  104. Witt, C. A., Charisi, M., Taylor, S. R., & Burke-Spolaor, S. 2022, ApJ, 936, 89 [NASA ADS] [CrossRef] [Google Scholar]
  105. Wu, Q., Scialpi, M., Liao, S., Mannucci, F., & Qi, Z. 2024, A&A, 692, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Xu, H., Chen, S., Guo, Y., et al. 2023, Res. Astron. Astrophys., 23, 075024 [CrossRef] [Google Scholar]
  107. Yu, W., Richards, G. T., Vogeley, M. S., Moreno, J., & Graham, M. J. 2022, ApJ, 936, 132 [NASA ADS] [CrossRef] [Google Scholar]
  108. Zhang, Y.-W., Huang, Y., Bai, J.-M., et al. 2021, AJ, 162, 276 [NASA ADS] [CrossRef] [Google Scholar]
  109. Zhu, X.-J., & Thrane, E. 2020, ApJ, 900, 117 [NASA ADS] [CrossRef] [Google Scholar]
  110. Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, ApJ, 765, 106 [NASA ADS] [CrossRef] [Google Scholar]
1

The efficiency of tidal stripping depends on the possible deepening of the potential well of the core, as a consequence of merger driven gas inflows (e.g. Callegari et al. 2009).

2

Notably, this red noise presents a major limitation in the photometric periodicity searches for close MBHBs, leading to false positives (Vaughan et al. 2016).

3

The documentation about the package can be found at the following link: https://celerite.readthedocs.io/en/stable/python/install/

4

The documentation about this package can be found at this link: https://github.com/wdpozzo/raynest

5

We emphasise again that the ratio obtained from fitting the Lomb-Scargle periodogram in unevenly sampled data should not be considered as a statistically robust indicator of the presence of an UBHD.

Appendix A: Periodogram computation

The periodogram is a useful tool to estimate the true PSD of a process (see Priestley 1981; Bloomfield 2004). In its classical definition, the periodogram can be defined as the square of the discrete Fourier transform (see Scargle 1982). Panels (b) and (d) of Figure 1 show examples of the periodograms computed using the classical definition.

The sampled periodogram is not a smooth function of the frequencies but is scattered around the true PSD (see Vaughan et al. 2003). More specifically, the periodogram computed at a given frequency is distributed as a χ2 distribution with two degrees of freedom, i.e. it is exponentially distributed as (see Scargle 1982; van der Klis et al. 1989):

p ( I j / S j ) = 1 S j e − I j S j Mathematical equation: $$ \begin{aligned} p(I_j/S_j) = \frac{1}{S_j}e^{-\frac{I_j}{S_j}} \end{aligned} $$(A.1)

where Sj and Ij are respectively the true PSD and the value of the periodogram computed at the frequency fj = j/NΔt with j = 1...N/2. The shape of this distribution is due to the fact that the real and imaginary parts of the discrete Fourier transform are normally distributed for a stochastic process (see Jenkins & Watts 1968; Priestley 1981). As an example, in Figure A.1, we show an example of a white noise light curve, its periodogram with the true PSD and the distribution of the periodogram around the true PSD, along with the theoretical χ2 distribution. The classical definition works particularly well for evenly sampled time series, but it is not without flaws. The periodogram is not a consistent estimator as the scatter does not decrease as the number of data points in the light curve increases (Jenkins & Watts 1968). Also, periodograms measured from finite data tend to be biased by windowing effects (see Vaughan et al. 2003, and references therein).

Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Example of a distribution for a white noise periodogram. Upper panel: noisy time series for a white noise with unity variance. Middle panel: comparison between the theoretical PSD (red dashed line) and the retrieved periodogram. Lower panel: distribution of the individual frequency bin estimates of the retrieved periodogram around the true PSD compared with the expected exponential distribution shown in Equation A.1 (red dashed lines).

When the sampling is uneven, one can use the Lomb-Scargle periodogram (see Lomb 1976; Scargle 1982). The Lomb-Scargle periodogram is a generalisation of the classical periodogram, but its statistical properties are different. In particular, the observed periodogram is distorted from the underlying PSD by the sampling pattern of the underlying light curve. More specifically, in the context of red noise, effects such as aliasing and red noise leakage become important (see Papadakis & McHardy 1995; Uttley et al. 2002). Furthermore, the distribution of periodogram peaks is generally unknown, and correlations between the power at different frequency bins arise unless the light curve is evenly sampled and the periodogram is computed at the Fourier frequencies or another orthonormal frequency grid (see Vio et al. 2010; Covino et al. 2022).

Appendix B: Fitting directly the PSD

The best parameters of the stochastic process can be estimated directly from the periodogram. Also in this case, one can use nested sampling to do both evidence and best parameter estimation. In particular, the likelihood for the sampled periodogram points is given by the Whittle likelihood (see Barret & Vaughan 2012, for an application to X-ray power spectra). Given the probability (A.1), the Whittle log-likelihood can be written as:

log L = ∑ j = 1 N − 1 I j S j + log S j . Mathematical equation: $$ \begin{aligned} \log L= \sum _{j = 1}^{N-1} \frac{I_j}{S_j}+ \log S_j \ . \end{aligned} $$(B.1)

This approach can be used only in light curves where the time is evenly sampled. In the case of uneven sampling, the computation of the periodogram through the Lomb-Scargle procedure introduces correlations between different frequencies so that the Whittle likelihood does not approximate the true likelihood of the periodogram anymore. In particular, the specific powers at each frequency bin are not always independent. In principle, the presence of an UBHD can be searched for by fitting a Lomb-Scargle periodogram from an unevenly sampled light curve, but the statistical issues mentioned above will prevent a robust statistical estimate of the evidence of the UBHD hypothesis. In the following, we present for completeness a few examples of evenly sampled mock UBHD light curves fitted with two red-noise PSDs, indicating the Bayes ratios obtained under the two assumptions5. In Figure B.1, we show the result of the fitting procedure using the PSD and Whittle likelihood for a time series sampled every 3 d for 10 yr with the same parameter ratios used in Figure 6. The dashed green line represents the true PSD, while the solid red and blue lines show the best fit from the UBHD and single MBH models, respectively. The red and blue shaded areas represent the 68% and 95% confidence intervals of the UBHD and single MBH models, respectively. Also in this case, we reported the Bayes ratio resulting from the nested sampling.

Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Results of the direct fit of the PSD for evenly sampled light curves with a 3 d cadence for 10 yr. The grey line indicates the light curve periodogram. The red and blue solid lines show the median model for the double and single DRW, respectively. The red and blue shaded areas show confidence intervals (16th − 84th and 2.5th − 97.5th ) for the UBHD and single DRW, respectively.

Another approach to assess the presence of an UBHD without a fully Bayesian analysis, which we do not explore here, is to generate a large number of simulated curves assuming a large grid of values for the parameters of the assumed noise model. By defining a suitable test-statistic, such as the χ2, one can find the set of parameters that best describe the observed light curve, and given a goodness-of-fit statistic, one can compare the results for different noise models. We note, though, that while such a procedure may reduce the computational cost of the fitting procedure, it typically provides only point estimates of the parameters and approximations to their uncertainties, and it does not fully account for correlations or degeneracies between parameters, in contrast to the fully Bayesian approach presented in this work.

All Tables

Table 1.

Parameter recovery for the single MBH model.

In the text
Table 2.

Parameter recovery for the UBHD model.

In the text
Table 3.

Parameter recovery for the correctly identified UBHD light curves.

In the text
Table 4.

Percentages of light curves for which all the parameters are retrieved with a relative error within the 10%, 20% and 50%.

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Examples of light curves generated using celerite varying the process parameters: damping timescale τ (panel a), and the variability amplitude σ ̂ Mathematical equation: $ \hat{\sigma} $ (panel c). In the right panels, the retrieved periodograms (computed as discussed in Appendix A) as well as the theoretical expected PSD (dashed lines) are shown. The time series ( F ̂ Mathematical equation: $ \hat{F} $) are centred at zero and have arbitrary units. The parameter σ and power (P) should be considered in terms of the same dimensional scale as F ̂ Mathematical equation: $ \hat{F} $, with P having units of the square of the units of F ̂ Mathematical equation: $ \hat{F} $.
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Distribution of relative errors for the retrieved parameters in the single MBH evenly sampled case.
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

Distribution of relative errors for the retrieved parameters in the UBHD evenly sampled case. Shown are the distributions for the cases in which the UBHD model is favoured (in green).
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Map of Bayes factors for 1500 unevenly sampled realisations of UBHD light curves with a ten-year baseline. The colour scale represents the base-10 logarithm of the Bayes factor, while the axes show the ratio of the model parameters between the two light curves. The red square, diamond, triangle and circle correspond to the light curves used to compute the periodograms shown in panels a,b,c, and d, respectively, in Figures 6 and B.1.
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Same kind of map as reported in Figure 4 but for 580 light curves in the region of the parameter space where the UBHDs can be correctly identified.
In the text
Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

PSD of double DRW kernels for different parameter ratios. The Bayes ratio resulting from the GP fitting is reported above each panel. In all the panels, the light curves are sampled every 3 d for 10 yr. The power is reported in arbitrary units. Case (a): Similar variability amplitudes but considerably different damping timescales. The test can correctly identify the UBHD (shown as the red square in Figure 4). Case (b): Similar amplitudes and similar damping timescales (shown as the red diamond in Figure 4). Case (c): Considerably different amplitudes but similar damping timescales (shown as the red triangle in Figure 4). Case (d): Considerably different amplitudes and damping timescales (shown as the red circle in Figure 4).
In the text
Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

Similar to the map in Figures 4 but with different observational baselines. In the upper panel, the light curves are observed for 6 years, while in the lower panel, the light curves are observed for 3 years.
In the text
Thumbnail: Fig. A.1. Refer to the following caption and surrounding text. Fig. A.1.

Example of a distribution for a white noise periodogram. Upper panel: noisy time series for a white noise with unity variance. Middle panel: comparison between the theoretical PSD (red dashed line) and the retrieved periodogram. Lower panel: distribution of the individual frequency bin estimates of the retrieved periodogram around the true PSD compared with the expected exponential distribution shown in Equation A.1 (red dashed lines).
In the text
Thumbnail: Fig. B.1. Refer to the following caption and surrounding text. Fig. B.1.

Results of the direct fit of the PSD for evenly sampled light curves with a 3 d cadence for 10 yr. The grey line indicates the light curve periodogram. The red and blue solid lines show the median model for the double and single DRW, respectively. The red and blue shaded areas show confidence intervals (16th − 84th and 2.5th − 97.5th ) for the UBHD and single DRW, respectively.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.