Toward the unambiguous identification of supermassive binary black holes through Bayesian inference

Supermassive binary black holes are thought to be common in the Universe as natural products of galaxy mergers (Begelman et al., 1980). They are likely to be the primary sources of nanohertz gravitational waves, which have been actively searched for using pulsar timing arrays in the past decade (see, e.g., Hobbs et al., 2010; Verbiest et al., 2016; Perera et al., 2019, and references therein). Current pulsar timing arrays are sensitive to individual binaries, with component black hole masses $\gtrsim 10^{9}M_{\odot}$ and binary orbital periods $\lesssim 10$ years (or orbital separations $\lesssim 0.01$ pc), up to distances of $\sim 200$ Mpc (Zhu et al., 2014; Babak et al., 2016; Aggarwal et al., 2019), though see Rosado et al. (2016) for the possibility of detecting high-redshift binaries. No detection has been made so far, although it has been suggested that the detection of a stochastic gravitational-wave background formed by the combined emission from binaries across the Universe is likely in a few years (e.g., Taylor et al., 2016; Cordes & McLaughlin, 2019).

The electromagnetic identification of sub-parsec supermassive binary black holes has proven to be challenging. Direct imaging of such close binaries is only possible for sources within $\sim 100$ Mpc via very-long-baseline radio interferometry observations. The closest binary with confident direct images, found in the radio galaxy 0402+379, has a separation of 7.3 pc (Rodriguez et al., 2006), corresponding to an orbital period of $\sim 10^{4}$ years (Bansal et al., 2017); see also Kharb et al. (2017) for a $0.35\,\mathrm{pc}$ binary candidate at $116\,\mathrm{Mpc}$ in the Seyfert galaxy NGC 7674. How binaries like these can reach milli-parsec orbital separations, where the emission of gravitational waves can efficiently drive the binary to merge, is referred to as the “final-parsec problem”, which remains an active area of theoretical investigations (e.g., Milosavljević & Merritt, 2001; Yu, 2002; Colpi, 2014; Khan et al., 2016; Ryu et al., 2018; Goicovic et al., 2018; Muñoz et al., 2020; Taylor et al., 2017; Chen et al., 2019). Finding a sub-parsec supermassive binary black hole would not only provide insights into the final-parsec problem, but also shed light on the expected stochastic gravitational-wave signal strength for pulsar timing arrays (Zhu et al., 2019; Goulding et al., 2019).

Periodically variable active galactic nuclei provide an interesting class of candidates for sub-parsec supermassive binary black holes. The most prominent object is OJ 287, a blazar which exhibits $12\,\mathrm{year}$ quasi-periodic outbursts in optical light curves that date back to the 19th century (Sillanpaa et al., 1988; Valtonen et al., 2008). This periodicity has been interpreted as the secondary black hole crossing the accretion disc of the primary in an eccentric orbit (Valtonen et al., 2016; Dey et al., 2019).¹¹1After this paper was submitted, Laine et al. (2020) reported Spitzer observations of the predicted Eddington flare from OJ 287, adding significant weight to the theory that it is, in fact, a supermassive binary black hole. In the past five years, with the growth of time-domain astronomy, there have been a large number of sub-parsec supermassive binary black hole candidates claimed by various groups.

Based on the Catalina Real-time Transient Survey (CRTS), Graham et al. (2015a) put forward 111 binary black hole candidates from an optical variability analysis of 243,500 quasars. Among them, PG1302$-$102 was the most significant candidate with a measured period of $1884\,\mathrm{days}$ (Graham et al., 2015b). This periodicity has been attributed to the relativistic Doppler boosting of emission from a mini-disc around the secondary black hole (D’Orazio et al., 2015). Charisi et al. (2016) reported 33 binary candidates from a periodicity search in a sample of 35,383 quasars in the photometric database of the Palomar Transient Factory. Liu et al. (2015) identified a periodic variability of $542\,\mathrm{days}$ in quasar PSO J334.2+01.4 using data from the Pan-STARRS1 Medium Deep Survey. However, such a detection was found to be insignificant in a subsequent analysis of extended data (Liu et al., 2016); see Liu et al. (2019) for a systematic search over 9000 quasars which resulted in one candidate in their extensive analysis.

Following the periodicity report of PG1302$-$102, Vaughan et al. (2016) cautioned that the stochastic variability of normal quasars (i.e., those that do not host binary black holes) can resemble a periodic feature in light curves that span only a few periods and highlighted the importance of careful evaluation of false alarm rate. Through Bayesian model comparison, Vaughan et al. (2016) found a red-noise model is significantly favoured over a sinusoidal model for PG1302$-$102 based on the CRTS data. More recently, Liu et al. (2018) revisited the periodicity of PG1302$-$102 using additional data from the All-Sky Automated Survey for Supernovae (ASAS-SN). They employed a maximum likelihood method to search for a periodic signal in the presence of quasar red noise and found that the inclusion of ASAS-SN data reduced the periodicity significance. Therefore, Liu et al. (2018) concluded that the binary black hole model was disfavoured for PG1302$-$102. Kovačević et al. (2019) proposed a model that posits a cold spot in the accretion disc of the primary black hole. Such a model can produce a perturbed sinusoidal feature and thus explain the apparent decrease in periodicity significance found by Liu et al. (2018).

Here, we propose a fully Bayesian framework for the identification of supermassive binary black hole candidates in time-domain electromagnetic surveys. It is capable of dealing with generic signal forms in the presence of red noise. Our work improves on previous studies in several ways. First, our method allows the inference of noise properties and signal parameters simultaneously and thus accounts for potential covariance between a periodic signal and quasar red noise. Second, it is robust to offsets in data collected with different surveys and possible over/under-estimation of measurement uncertainties. Third, we adopt a general form of quasar red noise and search for deviation from the commonly assumed damped random walk (DRW) model (Kelly et al., 2009). We also compare the sinusoidal signal hypothesis against a quasi-periodic oscillation (QPO) model based on behavior observed in many stellar-mass black hole X-ray binaries (see, e.g., Motta, 2016, for a recent review).

The remainder of this paper is organized as follows. In Section 2, we describe the Bayesian inference framework and introduce our signal and noise models. In Section 3, we apply the method to PG1302$-$102 with data from CRTS, ASAS-SN, and the Lincoln Near-Earth Asteroid Research (LINEAR) survey, and present our analysis results. In Section 4 we discuss various aspects of a periodicity search through simulations. Last, we provide concluding remarks and outline directions for future work in Section 5.

In this section, we describe the framework of Bayesian inference and model selection for the analysis of time-series data. We focus on the case of searching for periodicity in quasar light curves—quasar brightness measurements as a function of time. Assuming stationary Gaussian noise, the likelihood function for a quasar light curve is

where $\mathbf{d}$ is the time-series light curve data with length $N$, $\vec{\theta}_{n}$ and $\vec{\theta}_{s}$ include the noise and signal parameters, respectively. In this work we use measurements of optical magnitudes (i.e., the logarithmic light curve). The data $\mathbf{d}$ is modeled as

where $\mathbf{n}$ is the noise vector, which contains measurement uncertainties and additional intrinsic stochastic quasar variability; $\mathbf{m}$ is a constant vector with identical entries of $m$, and $m$ accounts for the mean magnitude and any constant offset (e.g., constant level of contamination due to host galaxy light); the signal vector $\mathbf{s}$ is given by

The signal parameters are $\vec{\theta}_{s}=\{A,\phi,f_{0}\}$, where the signal frequency $f_{0}$ is related to the period by $f_{0}=1/T_{0}$. In Equation (2), $\mathbf{C}_{ij}=\langle\mathbf{n}_{i}\mathbf{n}_{j}\rangle$ is the noise covariance matrix, which contains two components, $\mathbf{C}=\mathbf{C}^{w}+\mathbf{C}^{r}$, where $\mathbf{C}^{w}$ is a diagonal matrix that represents white noise, and $\mathbf{C}^{r}$ accounts for the stochastic quasar variability, usually termed “red noise.”

The white noise matrix takes the following form

where $\sigma_{i}$ is the measurement uncertainty for the $i$th observation, $\nu$ is a scale factor used to quantify the over/under-estimation of measurement uncertainties²²2This is equivalent to the EFAC (Error FACtor) parameter used in pulsar timing data analysis. An additional EQUAD (Error added in QUADrature) parameter could be introduced to account for potential systematic errors that are not included in measurement uncertainties., $\delta_{ij}$ is the Kronecker delta function.

In the DRW model (Kelly et al., 2009), the covariance function is an exponential function. The red-noise covariance matrix is

where $\hat{\sigma}^{2}$ is the intrinsic variance between observations on short timescales ($\sim 1$ day), and $\tau_{0}$ is the damping timescale, and $\tau_{ij}\equiv|t_{i}-t_{j}|$. The noise power spectral density of the DRW model is

To account for the possibility that the quasar red noise deviates from the DRW model (e.g., Zu et al., 2013; Guo et al., 2017), we also consider the following form of noise covariance matrix:

This is also called the stretched exponential function³³3In the context of describing relaxation in disordered systems, it is called the Kohlrausch-Williams-Watts function.. Note that the special case of $\gamma=1$ corresponds to the DRW model, $\gamma=0$ means white noise, and $\gamma=2$ reduces to the Gaussian function.

We further extend our red-noise covariance matrix to account for the QPO phenomenon

where $T_{q}$ is the period of the QPO. In the case of $\gamma=1$, the power spectral density of the above covariance function becomes

where $f_{q}=1/T_{q}$ is the QPO frequency. Note that in the limit of $T_{q}\rightarrow\infty$, the above equation reduces to Equation (6). In practice, this model is indistinguishable from the DRW model once $T_{q}$ is longer than the data span. The noise parameters are $\vec{\theta}_{n}=\{\nu,\hat{\sigma},\tau_{0},\gamma,T_{q}\}$. We illustrate the red-noise models and discuss their phenomenology in detail in Appendix A.

We use Bayesian model selection to quantify the statistical significance of the presence of periodic signals in quasar light curves. We start with Bayes’ theorem, which states that

Here $P(\vec{\theta}|\mathbf{d},\mathcal{H})$ is the posterior probability distribution function of parameters $\vec{\theta}$ given data $\mathbf{d}$ and hypothesis $\mathcal{H}$; $\mathcal{L}(\mathbf{d}|\vec{\theta},\mathcal{H})$ is the likelihood function given in Equation (2), which describes the probability of observing data given the hypothesis $\mathcal{H}$ and parameters $\vec{\theta}$. Meanwhile, $P(\vec{\theta}|\mathcal{H})$ is the prior distribution of parameters $\vec{\theta}$ while ${\cal Z}(\mathbf{d}|\mathcal{H})$ is the Bayesian evidence for hypothesis $\mathcal{H}$

Given the observational data, we wish to compare two hypotheses: $\mathcal{H}_{n}$ = the data are consistent with only noise (i.e., $\mathbf{s}=0$), and $\mathcal{H}_{s}$ = there is a periodic signal present in the data. The ratio of posterior probability, called the odds ratio, between hypotheses $\mathcal{H}_{s}$ and $\mathcal{H}_{n}$ is:

where $P(\mathcal{H}_{n})$ and $P(\mathcal{H}_{s})$ are the prior probability for hypotheses $\mathcal{H}_{n}$ and $\mathcal{H}_{s}$, respectively. Assuming equal prior probability for both hypotheses, Bayesian model selection is usually performed by computing the Bayes factor:

In this work, we are concerned with the support for a sinusoidal signal quantified by $\mathcal{B}^{s}_{n}$. We also wish to compute the support for the presence of QPO by comparing two models involving Equations (8) and (7), and for the deviation from the DRW model, by comparing two models involving Equations (7) and (5). Following Kass & Raftery (1995), the interpretation of Bayes factors is as follows. In a natural logrithmic scale, $0<\ln\mathcal{B}<1$ indicates that the support (for the hypothesis in the numerator) “worth no more than a bare mention,” $1<\ln\mathcal{B}<3$ implies positive support, $3<\ln\mathcal{B}<5$ implies strong support and $\ln\mathcal{B}>5$ implies very strong support. These thresholds are somewhat arbitrary; $\ln\mathcal{B}=8$ is often adopted as the detection threshold in gravitational-wave astronomy (e.g., Thrane & Talbot, 2019). Throughout this paper, the log evidence or log Bayes factor refer to the natural logarithm.

PG1302$-$102 is a nearby bright quasar with a median V-band magnitude of $15.0$ at a redshift of $z=0.2784$. We apply our method to V-band magnitude measurements collected with CRTS (Drake et al., 2009), ASAS-SN (Shappee et al., 2014; Jayasinghe et al., 2019), and LINEAR (Sesar et al., 2011). The CRTS data are publicly available⁴⁴4http://nesssi.cacr.caltech.edu/DataRelease/ as part of the Catalina Surveys Data Release 2, including 290 phototmetric measurements taken between 6 May 2005 and 30 May 2013. The ASAS-SN data are downloaded from its photometry database⁵⁵5https://asas-sn.osu.edu/photometry as part of their Data Release 9, including 232 measurements acquired between 23 November 2013 and 27 November 2018. The LINEAR data contain 626 measurements made between 23 January 2003 and 12 June 2012. The median measurement uncertainties are 0.06, 0.05 and 0.01 mag for CRTS, ASAS-SN and LINEAR data, respectively. The span of the combined data set is $15.84\,\mathrm{years}$. To reduce the computational cost, we average LINEAR data with an interval of $1\,\mathrm{day}$, which reduces the number of data points to $138$. The bin size is chosen such that any stochastic time-correlated variations over a timescale greater than $1\,\mathrm{day}$ are preserved. The new uncertainty for binned data is computed as the standard deviation of original measurements inside the binning window. The median uncertainty for averaged LINEAR data is $0.015\,\mathrm{mag}$.

Figure 1 shows the data used in our analysis: LINEAR in pink, CRTS in black and ASAS-SN in blue. The error bars indicate the 1-$\sigma$ measurement uncertainties. We also show the 68% credible region of the sinusoidal signal (grey shaded band) with a period of $5.66_{-0.10}^{+0.17}\,\mathrm{years}$ inferred from the data, and a red-noise realization (thin light blue line) that contains a QPO with a period of $5.56\,\mathrm{years}$ and might fit the data well.

Before we present details of our analysis results, we show in Figure 2 the Lomb-Scargle periodogram of PG1302$-$102 data (solid blue line), along with the averaged periodogram taken from $10^{3}$ noise realizations for three models: red noise (dotted red line), red noise plus a QPO (dash-dotted green line), and red noise plus a sinusoid (dashed cyan line). The vertical black dash line indicates the data span of $15.84\,\mathrm{years}$ for PG1302$-$102. We simulate data using the same sampling in time and error bars as real data to obtain the average periodogram for three models. The following median posterior probability model parameters as inferred from real data are adopted: 1) $\ln\hat{\sigma}^{2}=-2.63\,\text{mag}^{2}\,\text{yr}^{-1}$, $\ln(\tau_{0}/\textrm{yr})=2.05$, $\gamma=0.74$, for the red-noise hypothesis; 2) $T_{0}=5.56\,\mathrm{yr}$, $\ln\hat{\sigma}^{2}=-3.47\,\text{mag}^{2}\,\text{yr}^{-1}$, $\ln(\tau_{0}/\textrm{yr})=2.94$, $\gamma=0.52$, for the “red noise + QPO” hypothesis; 3) $A=0.13\,\mathrm{mag}$, $T_{0}=5.66\,\mathrm{yr}$, $\phi=1.64\,\mathrm{rad}$, $\ln\hat{\sigma}^{2}=-2.87\,\text{mag}^{2}\,\text{yr}^{-1}$, $\ln(\tau_{0}/\textrm{yr})=0.95$, $\gamma=0.57$, for the “red noise + sinusoid” hypothesis. The purpose of Figure 2 is to provide some insights into how different models might fit the data, instead of a proper accounting of the goodness of these models. The periodogram of real data exhibits a strong peak at $5.6\,\mathrm{years}$, with secondary peaks at around 1 year. As one can see, the major peak can be reproduced with both a QPO and sinusoidal model but not the pure red-noise model. The secondary peaks at around 1 year are probably artefacts of the irregular sampling of the data, since simulated data that contain no signal at that period also result in apparent peaks there.

Sub-parsec supermassive binary black holes are crucial in our understanding of galaxy evolution. They are also the primary sources of nanohertz gravitational waves highly anticipated by international pulsar timing array campaigns. While nearly impossible to resolve through direct imaging, such close binaries are expected to produce periodic variations in light curves of active galactic nuclei. New candidates of periodicity are reported on a nearly monthly basis.

In this work, we propose a fully Bayesian method for the identification of periodicity in astronomical time-series that exhibit red noise. We apply this method to one of the most promising periodicity candidates, PG1302$-$102, using data from CRTS, ASAS-SN and LINEAR surveys. Our main findings are:

1. 1.

   There is very strong support for the presence of either a sinusoidal signal ($\ln\mathcal{B}=12.7$) or a QPO ($\ln\mathcal{B}=14.5$) at a period of $5.6\,\mathrm{yr}$. The data slightly favours the QPO hypothesis with an odds ratio of 6.

2. 2.

   The inclusion of ASAS-SN data reduces the log Bayes factor that supports the presence of a sinusoidal signal by 1.6. This, combined with the fact that the posterior distribution of sinusoidal period is unstable when we combine new data with CRTS, provides further evidence against the sinusoidal hypothesis.

3. 3.

   There is also conclusive evidence for deviation from the DRW red-noise model, with $\ln\mathcal{B}=11.4$. The noise power spectral density is shallower than a power law with an index of 2, and the red-noise damping timescale is long $\gtrsim 10\,\mathrm{yr}$.

We perform simulations of sinusoidal signals and red noise and demonstrate the following:

1. 1.

   The growth of the periodicity significance with additional data is not guaranteed because of the stochastic fluctuations of red noise. For data sets like LINEAR, CRTS and ASAS-SN, the log Bayes factor is expected to grow once $\ln\mathcal{B}\gtrsim 8$ is achieved with initial data.

2. 2.

   Our method is capable of distinguishing a sinusoidal signal from a QPO. However, this is only possible for strong signals or a large number of wave cycles. A strong prior constraint on the red noise would also help.

3. 3.

   The use of data binning can reduce our ability to detect periodicity or deviation from the DRW model, because the binning throws away high-frequency information that helps estimate the spectral slope of the red noise (see Fig. 8).

We show that the data of PG1302$-$102 favours the presence of a quasi-periodicity at around $5.6\,\mathrm{years}$ against the noise hypothesis with a Bayes factor of $2\times 10^{6}$. However, given that PG1302$-$102 was selected as the most periodic candidate out of $2\times 10^{5}$ quasars, it cannot be ruled out that such a quasi-periodicity is caused by pure red noise. Assuming that it is indeed a QPO and the central black hole is $10^{8.5}M_{\odot}$ (Graham et al., 2015b), we infer a period corresponding to a Keplerian orbit at 344 Schwarzschild radii ($\approx 0.01\,\mathrm{pc}$). We note that the QPO found in PG1302$-$102 roughly corresponds to the low-frequency QPOs seen in stellar-mass black hole X-ray binaries, although the QPO frequency of PG1302$-$102 lies one order of magnitude below the linear scaling relation between black hole mass and QPO frequency fitted to several candidates in Smith et al. (2018). The physical origins of QPOs are poorly understood and therefore such a scaling relation might not exist if there are different mechanisms driving QPOs in small and big black holes. Nevertheless, the QPO hypothesis of PG1302$-$102 can be tested with continued observations, and if the quasi-periodicity signature is confirmed, it may have significant implications for accretion physics of supermassive black holes.

Finally, our Bayesian framework can be adopted to establish unambiguous binary black hole detections with the following extensions:

1. 1.

   Apply this method to various physical models for the periodicity, such as periodic mass accretion rate of the binary (e.g., Farris et al., 2014) or jet precession of a single or binary supermassive black holes (e.g., Abraham & Carrara, 1998; Kudryavtseva et al., 2011; Britzen et al., 2018), in addition to relativistic Doppler boosting;

2. 2.

   Use more sophisticated signal models that account for a cold-spot-induced perturbation in the accretion disc for PG1302$-$102 (Kovačević et al., 2019) or post-Newtonian orbital evolution for OJ 287 (Dey et al., 2018);

3. 3.

   Combine multi-wavelength observations (e.g., Xin et al., 2020), and other signatures associated with a binary black hole, for example, Doppler velocity offsets in broad emission line profiles (Eracleous et al., 2012; Bon et al., 2012; Li et al., 2016), and flux deficits in the spectral energy distribution (Gültekin & Miller, 2012; Guo et al., 2020). See Zheng et al. (2016) for a binary black hole candidate for which different types of signatures are analyzed;

4. 4.

   Use astrophysically-motivated population priors. For example, the period distribution of binary black hole population is expected to be dominated by long periods, and the distribution of red-noise parameters can be obtained by applying our method to a large number of active galactic nuclei using hierarchical inference.