Bayesian Synthesis of Astrometric Wobble and Total Light Curves in Close Binary Supermassive Black Holes

In the case of an unresolved binary system where the orbital motion of the components is significant, it is essential to consider the photocenter perturbation caused by this motion. This perturbation, ( δ _ph) approximated as the difference between the photocenter and the center of mass (see also Belokurov et al. 2020). Consequently, in CB-SMBH systems, the perturbation vector ${{\boldsymbol{\delta }}}_{\mathrm{ph}}=(\tilde{x}(t),\tilde{y}(t))$ (van Altena 2012) in CB-SMBH, can be approximated as follows:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\delta }}}_{\mathrm{ph}} & \propto & \displaystyle \frac{{M}_{1}{{\boldsymbol{\xi }}}_{1}+{M}_{2}{{\boldsymbol{\xi }}}_{2}}{{M}_{1}+{M}_{2}}-\displaystyle \frac{{F}_{1}{{\boldsymbol{\xi }}}_{1}+{F}_{2}{{\boldsymbol{\xi }}}_{2}}{{F}_{1}+{F}_{2}}\\ & = & (\mu -l(t))({{\boldsymbol{\xi }}}_{2}-{{\boldsymbol{\xi }}}_{1})\end{array}\end{eqnarray} \tag{ 1 }$$

where (M_i , F_i ), i = 1, 2 denote the masses and fluxes of the components in the binary. The mass and flux ratios are given by q = M₂/M₁ and f(t) = F₂/F₁, respectively. The reduced mass and flux are represented as $\mu =\tfrac{{M}_{2}}{{M}_{1}+{M}_{2}}$ and $l(t)=\tfrac{{F}_{2}}{{F}_{1}+{F}_{2}}$, respectively. The vectors ξ _i = (x_i , y_i ), i = 1, 2 represent the projected positions of the binary components on the sky. ⁴

In coordinates, the perturbation vector translates to:

$$\begin{eqnarray}&&{{\boldsymbol{\delta }}}_{\mathrm{ph}}\propto \left\{\begin{array}{l}\tilde{x}(t)=-(\mu -l(t)){s}_{x}(t)\\ \tilde{y}(t)=-(\mu -l(t)){s}_{y}(t)\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where, s (t) = (s_x (t) = x₂(t) − x₁(t), s_y (t) = y₂(t) − y₁(t)). A detailed explanation of this expression is provided in Appendix B (see Equation (B1)) and further context can be found in Wielen (1996).

The equations are rewritten in a more detailed form below (for further details, see Appendix B.1):

$$\begin{eqnarray}\begin{array}{rcl}\tilde{x}(t) & = & -(\mu -l(t)){s}_{x}(t)\\ & = & -(\mu -\,l(t))a\left[\left[\sin {\rm{\Omega }}\cos (\omega )\right.\right.\\ & & \left.+{ \mathcal I }\cos {\rm{\Omega }}\sin (\omega )\right]\left[\cos E(t)-e\right]\\ & & +\,\left[-\sin {\rm{\Omega }}\sin (\omega )\right.\\ & & \left.\left.+{ \mathcal I }\cos {\rm{\Omega }}\cos (\omega )\right]\right]\sqrt{1\,-{e}^{2}}\sin E(t),\\ \tilde{y}(t) & = & -(\mu -l(t)){s}_{y}(t)\\ & = & -(\mu -l(t))a\left[\left[\cos {\rm{\Omega }}\cos (\omega )\right.\right.\\ & & \left.-{ \mathcal I }\sin {\rm{\Omega }}\sin (\omega )\right]\left[\cos E(t)-e\right]\\ & & +\left[-\cos {\rm{\Omega }}\sin (\omega )\right.\\ & & \left.\left.-{ \mathcal I }\sin {\rm{\Omega }}\cos (\omega )\right]\right]\sqrt{1-{e}^{2}}\sin E(t)\end{array}\end{eqnarray} \tag{ 3 }$$

where a denotes the semimajor axis of the secondary, e is the eccentricity of the orbit, E(t) is the eccentric anomaly, and Ω, ω, and ${ \mathcal I }=\cos i$ are the Euler orbital angles as defined in Appendix B.1. The Euler angles Ω and ω, along with the orbital inclination i, play a crucial role in calculating the astrometric wobble as they influence the orientation of the binary orbit in space. Moreover, the quantity (μ − l(t))a, roughly corresponds to the semimajor axis of the photocenter orbit due to the Keplerian motion of the binary. These relations allow us to connect the physical parameters of the binary system (such as masses, fluxes, and orbital parameters) to the observable astrometric signatures, providing valuable insights into the dynamics of the system. In this formalism, the mass and light ratios (μ and l(t)) are important in shaping the astrometric wobble. An equivalence in the values of mass and luminosity ratios (μ ≈ l(t)) minimizes the astrometric wobble. In contrast, a pronounced disparity in these values (μ ≠ l(t)) amplifies the wobble.

The pioneering observations by the EHT of the shadow of a SMBH with mass (6.5 ± 0.7) × 10⁹ M_⊙ at the center of M87 (EHT Collaboration et al. 2019) have spurred interest in the search for binary companions within such massive systems. Safarzadeh et al. (2019) proposed that M87, situated in a region known for frequent galactic mergers (Hopkins et al. 2006), may host a binary companion. Prior to this, PSO J334.2028+01.4075, a radio-loud quasar at z = 2.060, was identified as a strong candidate for a CB-SMBH, supported by extensive photometric observations (Liu et al. 2015). Its periodic variability of 542 ± 15 days and the corresponding estimated black hole mass of $\mathrm{log}({M}_{\mathrm{tot}}/{M}_{\odot })=9.97\pm 0.50$ suggest an orbital separation of approximately $\sim {0.006}_{-0.003}^{+0.007}$ pc. This separation is indicative of a system that may be approaching the GW-driven orbital decay phase, located at the epoch of peak SMBH mergers. This underscores the capabilities of current surveys and foreshadows the potential of the upcoming Legacy Survey of Space and Time (LSST) in the discovery and characterization of CB-SMBHs.

Informed by these findings, we are encouraged to extend our exploration to binary CB-SMBHs with subannual orbital periods, within the upper mass range ∼10¹⁰ M_⊙ (Magorrian et al. 1998; Ghisellini et al. 2009; McConnell et al. 2012; Valtonen et al. 2012; Hlavacek-Larrondo et al. 2013; Walker et al. 2014; Ferré-Mateu et al. 2015; Ghisellini et al. 2015; Liu et al. 2015; Zuo et al. 2015; Saturni et al. 2016; Thomas et al. 2016; Yıldırım et al. 2016; Dullo et al. 2017; Ge et al. 2019; Mehrgan et al. 2019; Jeram et al. 2020; Onken et al. 2020; Mejía-Restrepo et al. 2022; Eilers et al. 2023; Nightingale et al. 2023).

The Figure 1 visualizes (see Section A) the estimated detectability CB-SMBHs with P → 1 yr, a range that is recognized as a "binary candidate desert" (see Charisi et al. 2016, and their Figure 12). We find that detectability of such systems is maximized for systems with greater total mass, lower redshifts, and orbits with higher inclinations, as these conditions yield stronger astrometric signals. Overlaid contour lines correspond to constant astrometric wobbles, establishing benchmarks relative to the detection capabilities of the astrometric instruments. The ngEHT is sensitive to wobbles ≥3 μas, while GRAVITY+ threshold lies ≥10 μas, which is also within the expected performance envelope of GAIA. It is important to note that astrometric wobbles ∼2 μas) are at the threshold of ngEHT detectability, indicating that objects producing such signals are on the cusp of observability. Given that the expected CB-SMBH orbital periods at this GW emission stage are months to years, these objects are accessible via multi-epoch observations with the ngEHT (D'Orazio & Loeb 2018). As for the GRAVITY+ candidates, CB-SMBH candidates with single, offset, moving broad emission lines (Dexter et al. 2020) are mostly observed at z ∼ 0.2, with apparent magnitudes of V ≲ 18 and K ≲ 15 (Runnoe et al. 2017; Guo et al. 2019). In our simulations, we sampled a binary configuration from this specific parameter space with orbital period of ≤1 yr, showcasing the capabilities of future synergies between astrometry and time-domain surveys, for identifying CB-SMBHs in regions once overlooked, thereby contributing to a more comprehensive understanding of binary supermassive black hole populations. However, we also examine one configuration with period 1.3 yr, to illustrate the dynamics of systems with periods marginally exceeding our primary interval, shedding light on their characteristics as they converge toward 1 yr from above.

We outline a two-tier Bayesian protocol (see also Kovačević et al. 2022b). The first tier involves a Simulator that generates observational multi-instrument data encompassing the astrometric wobble and optical light curves, predicated upon specified characteristics of the CB-SMBH. The second tier, the Bayesian Solver, processes this simulated data to extract "solutions" that reveal the decomposed flux of, for instance, the secondary component, along with the orbital parameters of the CB-SMBH.

2.2.1. First Tier: Simulator

2.2.2. Second Tier: Solver

As already mentioned, we use the Solver set of equations for the fitting procedure, but without the jitter term. Denoting the astrometric observations of photocenter motion by Θ^o , the flux observations by F^o _tot, then the posterior distribution of the orbital parameters of binary motion and flux of one of components factored out of total flux of binary (u = (a, e, P, i, ω, Ω, F_j ), j = 1, 2) via Gaussian process is:

$$\begin{eqnarray}&&\begin{array}{l}P(u| {{\rm{\Theta }}}^{o},{F}_{\mathrm{tot}}^{o})\propto P({{\rm{\Theta }}}^{o},{F}_{\mathrm{tot}}^{o}| u)\cdot P(\theta ,{\alpha }_{\mathrm{GP}})\\ =\,P({{\rm{\Theta }}}^{o}| \theta )\cdot P({F}_{j}| {F}_{1}+{F}_{2},{\alpha }_{\mathrm{GP}})\cdot P(\theta ,{\alpha }_{\mathrm{GP}}),j=1,2.\end{array}\end{eqnarray} \tag{ 4 }$$

The probability density functions contributing to the equation above could be obtained as follows. The function P(Θ^o ∣θ) is the log-likelihood of fitting the astrometry data:

$$\begin{eqnarray}&&{ \mathcal L }=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\left[\displaystyle \frac{{\left({x}_{i}-{x}_{\mathrm{model}}(\theta )\right)}^{2}}{{\sigma }_{{x}_{i}}^{2}}+\displaystyle \frac{{\left({y}_{i}-{y}_{\mathrm{model}}(\theta )\right)}^{2}}{{\sigma }_{{y}_{i}}^{2}}\right]\end{eqnarray} \tag{ 5 }$$

where N is the number of observed astrometric data points, x_i and y_i are the observed positions of the photocenter, x_model(θ) and y_model(θ) are the model-predicted positions of the photocenter using model parameters θ, and uncertainties of observations are ${\sigma }_{{x}_{i}}$ and ${\sigma }_{{y}_{i}}$. Next, a detailed derivation of P(F_i ∣F₁ + F₂, α_GP) is provided in Section C (see Equation (C27)). This probability density represents the log-likelihood of decomposing F_i from the total flux of a binary system using a Gaussian process with parameters α_GP. The prior probability density P(θ, α_GP) = P(θ)P(α_GP) is a simple product of priors on orbital and GP parameters (see Table 2). Then the posterior distribution $P(u| {{\rm{\Theta }}}^{o},{F}_{\mathrm{tot}}^{o})$, as given by Equation (5), is explored through Bayesian framework in Python PyMC package (Salvatier et al. 2016; Wiecki et al. 2022). This program samples the posterior probability of the model given the data, including prior probabilities for the parameters.

The priors for the orbital parameters are sampled from distributions as given in Table 2, with parameters based on inspection of both the total light curve and astrometric observation. For example, in the case of C1, the total light curve suggests a relatively stable and symmetric flux behavior, indicating a lower range for eccentricity (e) and a moderate range for orbital period (P). Astrometric observations reveal consistent positional changes indicative of a well-defined orbital plane, constraining parameters such as inclination (i). The priors for configuration C1 include a uniform distribution e = Uniform(0.1, 0.4) and a normal distribution for P centered around 2 yr. In addition, priors for i and Ω may be constrained to narrower ranges to reflect the stable orbital plane inferred from astrometric data. For the configuration C3, the total light curve exhibits significant flux variations and irregularities, suggesting a wider range for eccentricity (e) and a longer orbital period (P) to accommodate the observed dynamical behavior. Astrometric measurements reveal complex positional changes, so that a broader range (0, 2π) for ω and Ω priors is used. The priors for e span Uniform(0.0.01, 0.6) while P could follow a normal distribution centered around 7 yr. Additionally, priors for ω and Ω are Uniform(0, 2π) to capture the complexity observed in the astrometric data. The error priors are drawn from normal distributions that are centered at expected errors of artificial data and have a standard deviation of 5%.

For the GP model, we did not employ DRW for recovering light curves. The Bayesian framework allows the use of any type of GP kernels. While there is a plethora of GP kernels that can be applied, in this study, we specifically used the squared exponential kernel (SE), which has been tested by Kovačević et al. (2017):

$$\begin{eqnarray}&&k={\alpha }_{\mathrm{SE}}^{2}\exp \left(\displaystyle \frac{-{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{l}_{\mathrm{SE}}^{2}}\right)\end{eqnarray} \tag{ 6 }$$

where α_SE sets the amplitude of the Gaussian process and has the same units as the flux, and l_SE sets the length scale of the Gaussian process and has units of time t. The Bayesian framework provides a natural mechanism to determine the most probable values of α_SE and l through sampling priors of their ranges and optimization of Gaussian process models. For unfavorable configurations, more structured GP kernels could be used such as:

$$\begin{eqnarray*}&&k(X,X^{\prime} )=\left({\alpha }_{\mathrm{SE}}^{2}\exp \left(-\displaystyle \frac{{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{l}_{\mathrm{SE}}^{2}}\right)+\mathrm{diag}\left(D{\left(X,\theta \right)}^{(3+{\alpha }_{D})}\right)\right)\end{eqnarray*}$$

where covariance matrix of Doppler boost for binary parameters θ is given by $\mathrm{diag}\left(D{(X,\theta )}^{(3+{\alpha }_{D})}\right)$.

Our Bayesian synthesis method successfully estimated the orbital parameters of the binary quasar system across six configurations (Table 1), demonstrating the model's capability to accurately recover parameters (Table 3). The estimated values are generally in close agreement with the true values (Appendix Figures 7–9, with most parameters falling within the 1 HDI and confidence intervals ranging from 1 HDI to, in a few cases, 2–3 HDI for more challenging configuration of "q-accrete" C6 (Table 4). Configuration C1 benefits from an optimal match between observing cycles and the orbital period, supported by moderate eccentricity and mass ratio. This balance enables precise components' fluxes decomposition within 68% HDI, highlighting the effectiveness of the observational strategy in capturing the binary's stable dynamics over two complete cycles. C2 demonstrates the model's ability to accurately capture rapid dynamical changes within a shorter, 1 yr period with a single observing cycle. The higher eccentricity here necessitates dense observational coverage, which is achieved, resulting in flux decompositions aligning within 68% HDI. C3 faces challenges due to its extended 9 yr period and a lower L_oc = 2/3, which complicates the full capture of long-term dynamics, placing decomposed fluxes within 97.5% at the beginning of the observing cycle. This reflects the difficulty in tracking the binary's evolution with sparser data points. C4 and C5 contend with high eccentricity (e ∼ 0.7), impacting the flux decomposition's precision. The longer period and half-cycle observing rates introduce complexities in fully capturing the binary's dynamics, likely resulting in decomposed fluxes aligning within 97.5% sigma in the case of C4. However, due to the smaller period of C5 the cadence is enough granular to decompose fluxes within 68% HDI. Additionally, the relatively higher eccentricity and lower total mass in C5 compared to other configurations may contribute to more distinguishable flux variations that are easier to capture within the given observational sampling. C6, with a higher L_oc and a significant mass of "q-accrete" within a compact 1.3 yr period, faces intense dynamical variations. The alignment of decomposed fluxes within wider σ bands for C6 indicates the model's capacity to account for these rapid variations to a certain extent. However, the deviations from tighter sigma bands underscore the necessity for more nuanced observational strategies, such as increasing the sampling rate during critical phases or employing predictive modeling to enhance the temporal resolution of observations. Given the uncertainty surrounding the specific type of binary system at the outset of observation, the configurations examined here point out the importance of a flexible and adaptive observational approach, capable of accommodating the wide range of binary characteristics that might be encountered:

• 1.  
  Start with a baseline sampling rate, such as the 10 observations per year, and adjust based on early data analysis. If rapid flux variations or complex orbital dynamics are detected, increase the observation frequency accordingly.
• 2.  
  Multiwavelength observations can provide additional clues about the binary system's nature, aiding in the refinement of observational strategies to better target the detected features.
• 3.  
  Collaborations that allow for shared data and observations from different facilities and instruments can enhance the detection capabilities and flux decomposition.

An example of workflow could involve collecting LSST catalog data on variability, and providing this information to future astrometry facilities. These facilities can then use the catalog to optimize their observing strategy in real time based on modeling, ensuring efficient capture of dynamic binary behavior.

Here we predict and examine the astrometric wobbles of CB-SMBH across two distinct parameter spaces: the (q/(1 + q), l) plane and the total mass-period domain. Within each space, we incorporate variations in eccentricity and orbital phase. The purpose of the predictions is to assist and enhance multi-instrument observational strategies.

For the computations of the offset across a grid of binary parameters, a simplified equation Equation (B10) is employed, wherein the luminosity parameter l is assumed to be time invariant and both the luminosity and mass ratios are presumed to have a uniform distribution. This assumption allows us to focus exclusively on the dynamical parameters of the orbit, namely the semimajor axis a, eccentricity e, true anomaly ν, ⁶ inclination i, and orbital orientation ω. The formula for total photocenter displacement could have positive or negative sign depending on the values of μ and l. If μ > l, then the component with the larger mass contributes less to the total light of the system, and vice versa.

Given that the most significant effects are expected from objects within tens of Mpc (Dexter et al. 2020), we set the distance of the object to be roughly 70 Mpc. We then compute the offsets (Equation (B10)) for various binary parameter combinations in the (q/(1 + q), l) plane, under the assumption of a fixed inclination i = 0.

In Figure 5, the distinction between offsets for different values of the semimajor axis a is clear. In the top row, where a = 5 ld, the offsets are consistently lower than in the bottom row, with a = 10 ld, emphasizing the influence of semimajor axis on the offsets.

Zoom In Zoom Out Reset image size

The variation in eccentricity and true anomaly across the columns of Figure 5 illustrates the influence of these parameters on the resultant offsets. The uniformity in the first column (subplots a, d) for a circular orbit contrasts with the pronounced variation observed in the second (subplots b, e) and third columns (subplots c, f), especially at higher eccentricities. At a moderate eccentricity e = 0.5 and a true anomaly ν = π/2, the offset variation is more pronounced, especially when comparing it with the circular case in the first column. The slight differences between subplots (b) and (e) further highlight the enhanced offset sensitivity to the orbital phase at this eccentricity. With a high eccentricity, e = 0.9 and the orbital phase close to π, the offset values display a contrast in both subplots (c) and (f) with respect to other cases. This underscores the pronounced effects of having components near the apocenter of their elliptical orbit, producing more substantial displacements.

We now focus on the photocenter displacement, [ δ ph], as a function of the total mass and period, $[{\boldsymbol{\delta }}{ph}(\mathrm{log}[{M}_{\mathrm{tot}}/{M}_{\odot }],P)]$ (see Figure 6). This alternative phase plane is crucial due to its direct link to the binary's physical and orbital parameters. For this analysis, we reparametrize the binary's relative position, r, in Equation (B10) using the formula:

$$\begin{eqnarray}&&r=\displaystyle \frac{a(1-{e}^{2})}{1+e\cos \nu }={\left({P}^{2}{M}_{\mathrm{tot}}\right)}^{1/3}\displaystyle \frac{1-{e}^{2}}{1+e\cos \nu }.\end{eqnarray} \tag{ 7 }$$

This equation shows the dependencies on the total mass M_tot, orbital period P, eccentricity e, and the true anomaly ν. In an eccentric orbit, the radius r fluctuates, peaking at ν = π (the apocenter). Hence, the maximum displacements may occur when the binary system's orbit is highly eccentric (e = 0.9), and the orbital phase, defined by the true anomaly ν, is at π. This aligns with the apocenter, where the two bodies in the binary system are maximally separated.

Zoom In Zoom Out Reset image size

Based on Figure 6, we outline several trends dictated by dynamical parameters. In the top subplot, for circular systems encompassing total masses within 10⁷–10⁸ M_☉, there exists a balanced distribution of positive and negative offsets. However,this gets disrupted for masses in the range from 10⁸ to 10⁹ M_☉, wherein a pronounced tendency to positive offsets emerges. This suggests that as we progress to more massive systems, the distinction between the light-weighted center and mass-weighted center becomes more pronounced.

Increasing the eccentricity to 0.5 introduces deviations from circularity in the orbits. This scenario is set at ν = π/2. Here, the positive offset regions have enlarged compared to top subplot, especially for higher mass systems. This trend underscores the perturbative effects of higher eccentricities: as orbits deviate from perfect circles, gravitational interactions over the orbit become more varied, leading to pronounced photocenter displacements. The bottom subplot captures the systems at the most elliptical scenario, with the diametrically opposed components positions. The significant eccentricity profoundly impacts the offset distributions. The positive offset region has not only enlarged but also stretched toward higher masses. Notably, the boundaries defined by the yellow curve, representative of the barycenter parameter are most prominent here, illustrating how the offset distributions are inherently tied to this parameter, especially at such high eccentricities.

The plane of $(\mathrm{log}[{M}_{\mathrm{tot}}/{M}_{\odot }],P)$ reveals that the demarcation between positive and negative total offset values are delineated with the mass ratio μ = q/(1 + q). Specifically, if μ > l, then the component with the larger mass contributes less to the total emission of the system, and the total photocenter displacement [ δ ph] will be positive. This may suggest that the more massive component is less luminous (or active) or possibly obscured. On the other hand, if μ < l, then the component with the larger mass contributes more to the total light of the system, and the photocenter displacement [ δ ph] will be negative. This may suggest that the more massive component is more luminous (active) or less obscured. This relationship is particularly pronounced for total masses exceeding 10⁸ M_☉ and periods ≥1 yr, as larger masses and extended periods lead to larger orbital separations and, consequently, more significant photocenter displacements. Since our method provides the mass ratio and disentangles the light curve of the secondary, the total photocenter enables us to derive a preliminary indicator of the activity of the components in the CB-SMBH system.

Here we highlight further aspects that can modulate the interpretation of our results.

• 1.  
  Orbital Phase Dependancy: The photocenter offset exhibits a dependence on the orbital phase, with smaller variations at the pericenter compared to the apocenter. This could be seen from the ratio of photocenter offset values in pericenter Δ_p and apocenter Δ_ap inferred from Equation (B10):

  $$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Delta }}}_{{\rm{p}}}}{{{\rm{\Delta }}}_{\mathrm{ap}}}=\displaystyle \frac{1-e}{1+e}.\end{eqnarray} \tag{ 8 }$$

  However, we also see that, for CB-SMBH systems with mild eccentric orbits (0 < e < 0.3), the offset ratio tends to approach unity, implying that even at pericenter, the observed offset can be relatively nonneglibile. This phase-dependent behavior necessitates careful consideration when interpreting observed offsets, especially for systems with varying eccentricities (see related discussions in Kovačević et al. 2022b).
• 2.  
  Parallax and Proper Motion: Our analysis did not account for the possibility of parallax and proper motion. While we assumed that most short-period CB-SMBH systems exhibit small and likely imperceptible parallax, still, we cannot exclude the existence of systems with larger parallax. Furthermore, proper motion of quasars may arise from various sources, such as differential variability between close sources (see Makarov & Secrest 2022), or due to the motion of a relativistic jet, as in the case of PKS 0119 + 11 (Lambert et al. 2021), which could affect proper motion measurements. These factors merit further investigation and consideration when interpreting photocenter offset measurements (see Lambert et al. 2021; Makarov & Secrest 2022). According to Makarov & Secrest (2022) many proper motion quasars may be unresolved dual systems displaying the variable imposed motion effect, while a smaller number may be coincidence alignments with foreground stars creating weak gravitational lensing.
• 3.  
  Low-mass Galaxies with Black Holes: Our study primarily focused on massive CB-SMBH systems; however, it is essential to recognize that low-mass galaxies hosting black holes with masses approaching ∼10⁵ M_☉ may require a separate examination. Several campaigns in recent years (see, e.g., Reines et al. 2013; Molina et al. 2021) have shown that 50%–80% of low-mass galaxies M_g ∼ 10⁹–10¹⁰ M_☉ may harbor in their core low-mass black holes with masses approaching ∼10⁵ M_☉ (see Xin & Haiman 2021 and references therein). These black holes follow scaling relations similar to SMBHs (see, e.g., Volonteri & Natarajan 2009) and could impact the study of binary candidates, given their unique characteristics and potential contribution to active galactic nuclei (see Kimbrell et al. 2021).
• 4.  
  Large Binary Separations and Disk Alignment: The dynamics of CB-SMBH systems may change significantly at large binary separations or periods. Circumbinary gas accretion can dominate over gravitational wave emission, and the alignment of circumbinary disks with the orbital plane introduces complexities that deserve further investigation. Such systems may require a distinct treatment within the CB-SMBH framework (see Miller & Krolik 2013; Aly et al. 2015; D'Orazio & Di Stefano 2018). For example, any periodic emission from the binary minidisks or the circumbinary disk's inner edge will be obscured if the binary orbit is coplanar with a circumbinary disk.
• 5.  
  Variability of Accretion Disks: In cases where the source periodicity is primarily driven by the variability of the accretion disk rather than the Doppler boost, modeling the motion of the photocenter becomes challenging. The resulting curve may exhibit substantial noise and not form a closed ellipse, affecting our ability to model and interpret photocenter motion accurately (see also discussion in Kovačević et al. 2022b).
• 6.  
  Detection of Ultracompact Binaries: Extremely ultracompact binary systems with very short periods (see details in Xin & Haiman 2021) may present challenges for astrometric detection. The expected offsets for such systems may fall below the detection threshold for current astrometric instruments (e.g., GAIA threshold is 10 μas, see D'Orazio & Loeb 2019). Therefore, time-domain optical surveys like LSST and future gravitational wave observatories like LISA (Amaro-Seoane et al. 2017) may be necessary for detection and study of such systems.
• 7.  
  Newtonian and Relativistic Precession: Our analysis did not consider the effects of Newtonian precession or relativistic precession in compact binary systems. These effects can introduce additional complexities into the observed variability, and further investigation is required to assess their significance (see the example of OJ287 Laine et al. 2020).
• 8.  
  Contribution of Circumbinary Structures: The potential contribution of circumbinary structures, such as circumbinary disks or circumbinary broad-line regions, to the emission from CB-SMBH systems could impact the observed photocenter motion. The extent of this impact and its implications for binary framework studies remain subjects for future exploration (Dexter et al. 2020; Popović et al. 2021; Guo et al. 2022; Kovačević et al. 2022b; Simić et al. 2022).

The CB-SMBH is assumed to be a two-body system of SMBHs, such that their masses satisfy M₁ > M₂. The procedure for simulating astrometric data is discussed briefly below, and further information may be found in Kovačević et al. (2020b).

Naturally, dynamical parameters fully describe the SMBH position relative to the barycenter of system (see Table 5). The apparent relative orbit is the one of the secondary around the primary projected on the sky plane, and it can be determined from measurements of the relative position of the components obtained through astrometric imaging or interferometric observations. The projected spatial motion of the binary components is calculated within the reference frame centered on the primary component or barycenter and two axes in the plane tangent to the celestial sphere (Le Bouquin et al. 2013): the x-axis points north, while the y-axis points east. The z-axis runs parallel to the line of sight and points in the direction of rising radial velocities (positive radial velocities). Here we provide equation in the barycenter reference frame.

The transformations could be represented in the vectors P and Q (Thiele-Innes parameters) instead of utilizing cosine and sine terms of Euler rotations.

The vector of relative position s (t) = [s_x (t), s_y (t), s_z (t)] of an SMBH with respect to the barycenter is:

$$\begin{eqnarray}&&{\boldsymbol{s}}(t)={\boldsymbol{s}}({\bf{0}})+{\partial }_{t}{\boldsymbol{s}}(0)t+{\boldsymbol{P}}[\cos E(t)-e]+{\boldsymbol{Q}}\sqrt{1-{e}^{2}}\sin E(t),\end{eqnarray} \tag{ B1 }$$

where E(t) is the eccentric anomaly determined from the Kepler equation $E(t)-e\sin E(t)=2\pi n(t-{t}_{0}),n\,=\,{P}^{-1}$, and t₀ is set to 0 for simplicity. The vectors of barycenter position s (0) and velocity ∂_t s (0) are set to zero for simplicity. Thiele-Innes functions P and Q are defined as follows:

$$\begin{eqnarray}&&{\boldsymbol{P}}=a\left[{\boldsymbol{p}}\cos (\omega )+{\boldsymbol{q}}\sin (\omega )\right],\,\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&{\boldsymbol{Q}}=a\left[-{\boldsymbol{p}}\sin (\omega )\ +{\boldsymbol{q}}\cos (\omega )\right],\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&{\boldsymbol{p}}=(\sin {\rm{\Omega }},\cos {\rm{\Omega }},0),\,\end{eqnarray} \tag{ B4 }$$

$$\begin{eqnarray}&&{\boldsymbol{q}}=\left({ \mathcal I }\cos {\rm{\Omega }},-{ \mathcal I }\sin {\rm{\Omega }},\sin i\right),\,\end{eqnarray} \tag{ B5 }$$

where we assume that $a={M}_{1}{G}^{1/3}{(n\cdot {M}_{\mathrm{tot}})}^{1/3}$ is semimajor axis of secondary component, and $I=\cos i$. Therefore, the above equations describe relative orbit of secondary with respect to primary. The data matching the third coordinate of body position (s_z (t)) cannot be obtained.

The model for photocenter motion is now simply

$$\begin{eqnarray*}&&{\rm{\Theta }}(t)=[{{\rm{\Theta }}}_{x}(t),{{\rm{\Theta }}}_{y}(t)]={D}^{-1}(\mu -l(t))[{s}_{x}(t),{s}_{y}(t)],\end{eqnarray*}$$

where D is a distance, q = M₂/M₁, f(t) = F₂/F₁, $\mu \,=\tfrac{{M}_{2}}{{M}_{1}+{M}_{2}}=\tfrac{q}{1+q}$ and $l(t)=\tfrac{{F}_{2}}{{F}_{1}+{F}_{2}}=\tfrac{f(t)}{1\,+\,f(t)}$.

We assume that both masses and distance are known and that fluxes of both SMBH are obtained via simulation procedure given in Appendix B.2.2. We assume that the errors for each artificial observation are independent and identically distributed, resembling white noise at the level of 5% (see also Dexter et al. 2020). In order to avoid using the same model for both producing the observations and finding the solution, an additional jitter was added in the model when simulating the data.

B.1.1. Approximation of Photocenter Motion Using True Anomaly

In this section, we will reapproximate the motion of the photocenter using the true anomaly (ν). It represents the angle between the direction of periapsis (the point of closest approach to the central body) and the current position of the object. As such, it is a more intuitive parameter that can be visualized and understood easily. More importantly, the true anomaly proves particularly useful for predicting an object's future position on its orbit. This property is advantageous for the general exploration of the parameter grid space of CB-SMBH. The use of true anomaly allows for a more straightforward and less computationally intensive approximation.

Let the two components of binary be at the relative position given by $({{\boldsymbol{\xi }}}_{2}-{{\boldsymbol{\xi }}}_{1})=\tfrac{a(1-{e}^{2})}{1\,+\,e\cos \nu }{{ \mathcal P }}^{{\rm{T}}}(i,\omega ,{\rm{\Omega }})(\cos \nu ,\sin \nu ,0)$, where ^T marks transposition, ${ \mathcal P }$ is orbit orientation matrix of the binary, and $\tfrac{a(1-{e}^{2})}{1\,+\,e\cos \nu }(\cos \nu ,\sin \nu ,0)$ is polar coordinates in the orbital plane as given in Appendix A in Kovačević et al. (2020b). We then rewrite the components s_x , s_y given in Equation (B1) in terms of orbital phase ν:

$$\begin{eqnarray}&&{s}_{x}=r(\cos ({\rm{\Omega }})\cos (\nu +\omega )-\sin ({\rm{\Omega }})\sin (\nu +\omega )\cos i)\end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray}&&{s}_{y}=r(\sin ({\rm{\Omega }})\cos (\nu +\omega )+\cos ({\rm{\Omega }})\sin (\nu +\omega )\cos i)\end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray}&&r=\displaystyle \frac{a(1-{e}^{2})}{1+e\cos \nu }.\end{eqnarray} \tag{ B8 }$$

Then we will calculate the total photocenter motion as:

$$\begin{eqnarray}&&[{{\boldsymbol{\delta }}}_{{ph}}]=(\mu -l(t))\sqrt{{s}_{x}^{2}+{s}_{y}^{2}}\end{eqnarray} \tag{ B9 }$$

$$\begin{eqnarray}&&\sim (\mu -l(t))r\sqrt{1-{\sin }^{2}(\nu +\omega ){{\sin }}^{2}i}.\end{eqnarray} \tag{ B10 }$$

At higher viewing angles, the detected binary radiation would be strongly modulated by relativistic effects such as gravitational lensing and Doppler boosting. Here, we consider the case when the emission coming from two minidisks surrounding each component in the system will appear Doppler boosted (Haiman et al. 2009; Hu et al. 2020). We schematically decompose the total observed flux ${\tilde{F}}_{\mathrm{tot}}$ of binary as (see, e.g., Hu et al. 2020; Gutiérrez et al. 2022)

$$\begin{eqnarray}&&{\tilde{F}}_{\mathrm{tot}}={\bar{F}}_{1}(1+{\rm{\Delta }}{F}_{{ \mathcal D }1})+{F}_{1}+{\bar{F}}_{2}(1+{\rm{\Delta }}{F}_{{ \mathcal D }2})+{F}_{2}\end{eqnarray} \tag{ B11 }$$

$$\begin{eqnarray}&&\begin{array}{l}=({\bar{F}}_{1}+{\bar{F}}_{2})\Space{0ex}{0.25ex}{0ex}(1+{\rm{\Delta }}{F}_{{ \mathcal D }1}\displaystyle \frac{{\bar{F}}_{1}}{{\bar{F}}_{1}+{\bar{F}}_{2}}\\ +{\rm{\Delta }}{F}_{{ \mathcal D }2}\displaystyle \frac{{\bar{F}}_{2}}{{\bar{F}}_{1}+{\bar{F}}_{2}}\Space{0ex}{0.25ex}{0ex})+{F}_{1}+{F}_{2}\end{array}\end{eqnarray} \tag{ B12 }$$

where the F_i , ${\bar{F}}_{i},{\rm{\Delta }}{F}_{{ \mathcal D }i},i=1,2$ account for stochastic variability, corresponding average values and Doppler boost term of each minidisk in the system, respectively.

The next two sections provide formulas for Doppler brightening and stochastic variability calculation for each minidisk.

B.2.1. Doppler Boost

At a given photon frequency ν_p , assuming that the minidisk emission is ${F}_{i}\propto {\nu }_{p}^{\alpha }$ the Doppler effect is (D'Orazio & Di Stefano 2018)

$$\begin{eqnarray}&&{\rm{\Delta }}{F}_{{ \mathcal D }i}=(3-{\alpha }_{{\nu }_{p}})({v}_{i}/c)\cos \phi \sin I,i=1,2\end{eqnarray} \tag{ B13 }$$

where ${\alpha }_{{\nu }_{p}}$ I, ϕ, and v_i are the power-law spectrum index, orbital inclination, phase of binary, and three-dimensional velocity of component, respectively. The projection of the velocity vector onto the line of sight is given as

$$\begin{eqnarray}&&{v}_{i}={K}_{i}(\cos (\omega +\nu )+e\cos (\omega ),i=1,2\end{eqnarray} \tag{ B14 }$$

$$\begin{eqnarray}&&{K}_{1}=\displaystyle \frac{2\pi }{P}\displaystyle \frac{q}{1+q}\displaystyle \frac{a\sin I}{\sqrt{1-{e}^{2}}}\end{eqnarray} \tag{ B15 }$$

$$\begin{eqnarray}&&{K}_{2}=\displaystyle \frac{2\pi }{P}\displaystyle \frac{1}{1+q}\displaystyle \frac{a\sin I}{\sqrt{1-{e}^{2}}}\end{eqnarray} \tag{ B16 }$$

where ω is the argument of pericenter and ν is the true anomaly. We solve for the radial velocity for both binary components, a prescription given in Kovačević et al. (2020a, 2022b) for binary mass ratio q < 1. Both components are assumed to be Doppler-modulated with corresponding line-of-sight velocities and that emission from both black holes share the same spectral index with a value ${\alpha }_{{\nu }_{{\rm{p}}}}\sim -0.44$ of the composite quasar spectrum (Vanden Berk et al. 2001).

B.2.2. Stochastic Variability

It is predicted that the quasars high-quality light curves obtained from the LSST and other large time-domain facilities will be not only shifted in time but also distorted at different wavelengths (Chan et al. 2020; Kovačević et al. 2022a). We adopted this approach to model the stochastic emission F_i of each minidisk in the context of the lamp-post model, thin-disk geometry, and damped random walk (DRW) for the driving function, where the size of the transfer function expands with wavelength.

For this purpose, we will introduce auxiliary scaling factors. To approximate the mass accretion rate ${\dot{M}}_{\mathrm{tot}}$ of the binary system in terms of the total mass of the binary M_tot, we can use the Eddington accretion rate:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{tot}}=4\pi {m}_{{\rm{p}}}{{GM}}_{\mathrm{tot}}/\epsilon c{\sigma }_{{\rm{T}}}\propto \rho {M}_{\mathrm{tot}}\end{eqnarray} \tag{ B17 }$$

where m_p is the proton mass, σ_T is the Thompson scattering cross-section,  is a typical radiative efficiency of quasar, and ρ is a constant. For the symmetrical binary systems q = 1, it is expected that ${\dot{M}}_{1}={\dot{M}}_{2}$, and as the q decreases, the ratio $\eta ={\dot{M}}_{2}/{\dot{M}}_{1}$ increases (Raimundo & Fabian 2009; Shankar et al. 2009; Duffell et al. 2020). Therefore, if the implied accretion rate of the binary is ${\dot{M}}_{\mathrm{tot}}={\dot{M}}_{1}+{\dot{M}}_{2}$, then the accretion rate is split between the two black holes with ${\dot{M}}_{i}={f}_{i}{M}_{\mathrm{tot}},i=1,2$ where f₁ = (1 + η)⁻¹ ρ and ${f}_{2}={(1+{\eta }^{-1})}^{-1}\times \rho$. We assume that η = (1 + 0.9q)⁻¹ as given in Duffell et al. (2020), if not otherwise stated. Each minidisk in a binary is set up as a nonrelativistic thin disk emitting blackbody radiation (Shakura & Sunyaev 1973; Thorne 1974), with the central source treated as a "lamp post" placed near the black hole (Cackett et al. 2007; Starkey et al. 2017). Then, the temperature profiles T_i , i = 1, 2 of each minidisk are described as (Cackett et al. 2007; Starkey et al. 2017; Chan et al. 2020; Kovačević et al. 2022a; Davelaar & Haiman 2022)

$$\begin{eqnarray}&&{T}_{i}\propto {\left({\left(\displaystyle \frac{r}{{R}_{0}}\right)}^{-3}\left(1-\sqrt{\displaystyle \frac{{r}_{0}}{r}}\right)\right)}^{1/4}\end{eqnarray} \tag{ B18 }$$

where the scale radius is

$$\begin{eqnarray}&&{R}_{0}={\left({\left(\displaystyle \frac{k\lambda }{{hc}}\right)}^{4}\displaystyle \frac{3{{GM}}_{i}{\dot{M}}_{i}}{8\pi \sigma }\right)}^{1/3}={\left({\left(\displaystyle \frac{k\lambda }{{hc}}\right)}^{4}\displaystyle \frac{3{{GM}}_{\mathrm{tot}}^{2}{f}_{i}\tilde{q}}{8\pi \sigma }\right)}^{1/3},\end{eqnarray} \tag{ B19 }$$

and λ is the rest-frame wavelength, σ is the Stefan-Boltzmann constant, k is Boltzmann constant, G is a gravitational constant, M_tot is total mass of the binary, r₀ = 6GM_i /c² is the inner radius of minidisk surrounding mass M_i , and

$$\begin{eqnarray}\tilde{q}=\left\{\begin{array}{cc}\displaystyle \frac{1}{1+q} & \mathrm{for}{M}_{1}\\ \displaystyle \frac{q}{1+q} & \mathrm{for}{M}_{2}\end{array}\right..\end{eqnarray} \tag{ B20 }$$

Using Equations (B18)–(B20) and following prescription in Chan et al. (2020), Kovačević et al. (2022a), the transfer function for each disk at given λ is:

$$\begin{eqnarray}&&\begin{array}{l}{\psi }_{i}(\tau ,\lambda )\propto {\displaystyle \int }_{{r}_{0}}^{{r}_{{out}}}{\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{\zeta {e}^{\zeta }}{{\left({e}^{\zeta }-1\right)}^{2}}{rdrd}\theta \cos i\times \\ \delta (\tau -(r/c)(1-\sin i\cos \theta ))\end{array}\end{eqnarray} \tag{ B21 }$$

$$\begin{eqnarray}&&\zeta =\displaystyle \frac{{hc}}{k\lambda {T}_{i}}\end{eqnarray} \tag{ B22 }$$

where the outer disk radius is approximated as ${r}_{\mathrm{out}}\,\sim 0.5{r}_{0}{({M}_{i})}^{2/3}$ (Jiménez-Vicente et al. 2014), θ is azimuth of reprocessing site, and i is inclination of the disk.

Finally, assuming that each minidisk has an own driving-source emission ${{ \mathcal F }}_{i},i=1,2$, the total minidisk flux F_i (λ, t), i = 1, 2 is obtained by convolution Tie & Kochanek (2018):

$$\begin{eqnarray}&&{F}_{i}(\lambda ,t)\sim \displaystyle \sum _{\tau ={\tau }_{0}}^{{\tau }_{\max }}{\psi }_{i}(\tau ,\lambda ){{ \mathcal F }}_{i}(t-\tau ){\rm{\Delta }}\tau .\end{eqnarray} \tag{ B23 }$$

Here, ${{ \mathcal F }}_{i}(t-\tau )$ is the fractional luminosity variability "lagged" by the light travel time τ = r/c from the disk center. We represent this driving light curve by the Damped Random Walk (DRW) model with the characteristic amplitude σ_DRW and timescale $\tilde{\tau }$ (Kelly et al. 2009; Kozłowski 2017). The DRW is constructed from the procedure given in Kovačević et al. (2021).

Here we provide an auxiliary formula for the multivariate Gaussian conditional distribution x _A ∣ x _B , which is used for likelihood of additive decomposition of Gaussian process (see Mardia et al. 1979; Rasmussen 2004; Duvenaud et al. 2013).

Let x follow a multivariate normal distribution:

$$\begin{eqnarray}&&x\sim { \mathcal N }(\mu ,{\rm{\Sigma }}).\end{eqnarray} \tag{ C1 }$$

Then, the conditional distribution of any subset vector x_A , given the complement vector x_B , is also a multivariate normal distribution:

$$\begin{eqnarray}&&{x}_{A}| {x}_{B}\sim { \mathcal N }({\mu }_{A| B},{{\rm{\Sigma }}}_{A| B})\end{eqnarray} \tag{ C2 }$$

where the conditional mean and covariance are:

$$\begin{eqnarray}\begin{array}{rcl}{\mu }_{A| B} & = & {\mu }_{A}+{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}({x}_{B}-{\mu }_{B})\\ {{\rm{\Sigma }}}_{A| B} & = & {{\rm{\Sigma }}}_{{AA}}-{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}\end{array}\end{eqnarray} \tag{ C3 }$$

with blockwise mean and covariance defined as:

$$\begin{eqnarray}\begin{array}{rcl}\mu & = & \left[\begin{array}{l}{\mu }_{A}\\ {\mu }_{B}\end{array}\right]\\ {\rm{\Sigma }} & = & \left[\begin{array}{ll}{{\rm{\Sigma }}}_{{AA}} & {{\rm{\Sigma }}}_{{AB}}\\ {{\rm{\Sigma }}}_{{BA}} & {{\rm{\Sigma }}}_{{BB}}\end{array}\right].\end{array}\end{eqnarray} \tag{ C4 }$$

To prove the concept of C3, we will assume that C4 and following equation holds:

$$\begin{eqnarray}&&x=\left[\begin{array}{l}{x}_{A}\\ {x}_{B}\end{array}\right]\end{eqnarray} \tag{ C5 }$$

where x_A is an n_A × 1 vector, x_B is an n_B × 1 vector, and x is an n_A + n_B = n × 1 vector. Then, we recall that by construction, the joint distribution of x_A and x_B is:

$$\begin{eqnarray}&&{x}_{A},{x}_{B}\sim { \mathcal N }(\mu ,{\rm{\Sigma }}).\end{eqnarray} \tag{ C6 }$$

Moreover, the marginal distribution of x_B follows from a multivariate normal distribution C1 and C4 as

$$\begin{eqnarray}&&{x}_{B}\sim { \mathcal N }({\mu }_{B},{{\rm{\Sigma }}}_{{BB}}).\end{eqnarray} \tag{ C7 }$$

According to the Bayes law of conditional probability, it holds that

$$\begin{eqnarray}&&p({x}_{A}| {x}_{B})=\displaystyle \frac{p({x}_{A},{x}_{B})}{p({x}_{B})}.\end{eqnarray} \tag{ C8 }$$

Applying C6 and C7 to C8, we have:

$$\begin{eqnarray}&&p({x}_{A}| {x}_{B})=\displaystyle \frac{{ \mathcal N }(x;\mu ,{\rm{\Sigma }})}{{ \mathcal N }({x}_{B};{\mu }_{B},{{\rm{\Sigma }}}_{{BB}})}.\end{eqnarray} \tag{ C9 }$$

Using the probability density function of the multivariate normal distribution, we get:

$$\begin{eqnarray}\begin{array}{rcl}p({x}_{A}| {x}_{B}) & = & \displaystyle \frac{\tfrac{1}{\sqrt{{\left(2\pi \right)}^{n}| {\rm{\Sigma }}| }}\cdot \exp \left[-\tfrac{1}{2}{\left(x-\mu \right)}^{{\rm{T}}}{{\rm{\Sigma }}}^{-1}(x-\mu )\right]}{\tfrac{1}{\sqrt{{\left(2\pi \right)}^{{n}_{B}}| {{\rm{\Sigma }}}_{{BB}}| }}\cdot \exp \left[-\tfrac{1}{2}{\left({x}_{B}-{\mu }_{B}\right)}^{{\rm{T}}}{{\rm{\Sigma }}}_{{BB}}^{-1}({x}_{B}-{\mu }_{B})\right]}\\ & = & \displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{n-{n}_{B}}}}\cdot \sqrt{\displaystyle \frac{| {{\rm{\Sigma }}}_{{BB}}| }{| {\rm{\Sigma }}| }}\cdot \exp \left[-\displaystyle \frac{1}{2}{\left(x-\mu \right)}^{{\rm{T}}}{{\rm{\Sigma }}}^{-1}(x-\mu )\right.\\ & & +\,\left.\displaystyle \frac{1}{2}{\left({x}_{B}-{\mu }_{B}\right)}^{{\rm{T}}}{{\rm{\Sigma }}}_{{BB}}^{-1}({x}_{B}-{\mu }_{B})\right].\end{array}\end{eqnarray} \tag{ C10 }$$

Writing the inverse of Σ as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}^{-1}=\left[\begin{array}{l}{{\rm{\Sigma }}}^{{AA}}{{\rm{\Sigma }}}^{{AB}}\\ {{\rm{\Sigma }}}^{{BA}}{{\rm{\Sigma }}}^{{BB}}\end{array}\right]\end{eqnarray} \tag{ C11 }$$

and applying C4–C10, we get:

$$\begin{eqnarray}&&\begin{array}{l}p({x}_{A}| {x}_{B})=\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{n-{n}_{B}}}}\cdot \sqrt{\displaystyle \frac{| {{\rm{\Sigma }}}_{{BB}}| }{| {\rm{\Sigma }}| }}\cdot \\ \exp \left[-\displaystyle \frac{1}{2}{\left(\left[\begin{array}{l}{x}_{A}\\ {x}_{B}\end{array}\right]-\left[\begin{array}{l}{\mu }_{A}\\ {\mu }_{B}\end{array}\right]\right)}^{{\rm{T}}}\left[\begin{array}{ll}{{\rm{\Sigma }}}^{{AA}} & {{\rm{\Sigma }}}^{{AB}}\\ {{\rm{\Sigma }}}^{{BA}} & {{\rm{\Sigma }}}^{{BB}}\end{array}\right]\left(\left[\begin{array}{l}{x}_{A}\\ {x}_{B}\end{array}\right]-\left[\begin{array}{l}{\mu }_{A}\\ {\mu }_{B}\end{array}\right]\right)\right]\\ \exp \left[+\displaystyle \frac{1}{2}\,{\left({x}_{B}-{\mu }_{B}\right)}^{{\rm{T}}}\,{{\rm{\Sigma }}}_{{BB}}^{-1}\,({x}_{B}-{\mu }_{B})\right].\end{array}\end{eqnarray} \tag{ C12 }$$

Rearranging C12, we have

$$\begin{eqnarray}\begin{array}{rcl}p({x}_{A}| {x}_{B}) & = & \displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{n-{n}_{B}}}}\cdot \sqrt{\displaystyle \frac{| {{\rm{\Sigma }}}_{{BB}}| }{| {\rm{\Sigma }}| }}\cdot \\ & & \exp \left[-\displaystyle \frac{1}{2}\left({\left({x}_{A}-{\mu }_{A}\right)}^{{\rm{T}}}{{\rm{\Sigma }}}^{{AA}}({x}_{A}-{\mu }_{A})\right.\right.+\\ & & 2{\left({x}_{A}-{\mu }_{A}\right)}^{{\rm{T}}}{{\rm{\Sigma }}}^{{AB}}({x}_{B}-{\mu }_{B})+{\left({x}_{B}-{\mu }_{B}\right)}^{{\rm{T}}}{{\rm{\Sigma }}}^{{BB}}({x}_{B}-{\mu }_{B})\\ & & +\displaystyle \frac{1}{2}\left.\left.{\left({x}_{B}-{\mu }_{B}\right)}^{{\rm{T}}}{{\rm{\Sigma }}}_{{BB}}^{-1}({x}_{B}-{\mu }_{B})\right]\right].\end{array}\end{eqnarray} \tag{ C13 }$$

Based on the inverse of a block matrix, the inverse of Σ in C11 is:

$$\begin{eqnarray}\begin{array}{l}{\left[\begin{array}{ll}{{\rm{\Sigma }}}_{{AA}} & {{\rm{\Sigma }}}_{{AB}}\\ {{\rm{\Sigma }}}_{{BA}} & {{\rm{\Sigma }}}_{{BB}}\end{array}\right]}^{-1}=\left[\begin{array}{ll}X & Y\\ Z & W\end{array}\right],\end{array}\end{eqnarray} \tag{ C14 }$$

where block matrices are:

$$\begin{eqnarray*}\begin{array}{rcl}X & = & {\left({{\rm{\Sigma }}}_{{AA}}-{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}\right)}^{-1},\\ Y & = & -{\left({{\rm{\Sigma }}}_{{AA}}-{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}\right)}^{-1}{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1},\\ Z & = & -{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}{\left({{\rm{\Sigma }}}_{{AA}}-{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}\right)}^{-1},\\ W & = & {{\rm{\Sigma }}}_{{BB}}^{-1}+{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}{\left({{\rm{\Sigma }}}_{{AA}}-{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}\right)}^{-1}{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}.\end{array}\end{eqnarray*}$$

We also recall that the determinant of a block matrix is:

$$\begin{eqnarray}\left|\begin{array}{ll}{{\rm{\Sigma }}}_{{AA}} & {{\rm{\Sigma }}}_{{AB}}\\ {{\rm{\Sigma }}}_{{BA}} & {{\rm{\Sigma }}}_{{BB}}\end{array}\right|=| {{\rm{\Sigma }}}_{{BB}}| | {{\rm{\Sigma }}}_{{AA}}-{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}| .\end{eqnarray} \tag{ C15 }$$

Plugging C14 and C15 into C13 and doing some algebra operations, we have:

$$\begin{eqnarray}\begin{array}{rcl}p({x}_{A}| {x}_{B}) & = & \displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{n-{n}_{B}}}}\cdot \sqrt{\displaystyle \frac{| {{\rm{\Sigma }}}_{{BB}}| }{| {\rm{\Sigma }}| }}\cdot \\ & & \exp \left[-\displaystyle \frac{1}{2}\left[{\left({x}_{A}-{\mu }_{A}\right)}^{{\rm{T}}}{M}^{-1}({x}_{A}-{\mu }_{A})\right.\right.\\ & & -2{\left({x}_{A}-{\mu }_{A}\right)}^{{\rm{T}}}{M}^{-1}{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}({x}_{B}-{\mu }_{B})\\ & & \left.\left.+{\left({x}_{B}-{\mu }_{B}\right)}^{{\rm{T}}}N({x}_{B}-{\mu }_{B})\right]\right],\end{array}\end{eqnarray} \tag{ C16 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}M & = & {{\rm{\Sigma }}}_{{AA}}-{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}},\\ N & = & {{\rm{\Sigma }}}_{{BB}}^{-1}+{{\rm{\Sigma }}}_{{BB}}^{-1}{{\rm{\Sigma }}}_{{BA}}{M}^{-1}{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}.\end{array}\end{eqnarray*}$$

With ${{\rm{\Sigma }}}_{{BA}}={{\rm{\Sigma }}}_{{AB}}^{{\rm{T}}}$, because Σ is a covariance matrix, and n − n_B = n_A , we finally arrive at

$$\begin{eqnarray}\begin{array}{rcl}p({x}_{A}| {x}_{B}) & = & \displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{{n}_{A}}| M| }}\\ & & \exp \left[-\displaystyle \frac{1}{2}\left[{\left({x}_{A}-{\mu }_{A}\right)}^{{\rm{T}}}{M}^{-1}({x}_{A}-{\mu }_{A})\right.\right.\\ & & -2\left.\left.{\left({x}_{A}-{\mu }_{A}\right)}^{{\rm{T}}}{M}^{-1}{{\rm{\Sigma }}}_{{AB}}{{\rm{\Sigma }}}_{{BB}}^{-1}({x}_{B}-{\mu }_{B})\right]\Space{0ex}{3ex}{0ex}\right].\end{array}\end{eqnarray} \tag{ C17 }$$

The above equation is the probability density function of a conditional multivariate normal distribution:

$$\begin{eqnarray}&&p({x}_{A}| {x}_{B})={ \mathcal N }({x}_{A};{\mu }_{A| B},{{\rm{\Sigma }}}_{A| B})\end{eqnarray} \tag{ C18 }$$

with the mean μ_A∣B and covariance Σ_A∣B given by C3.

The CB-SMBH observed light curve is F_tot(t) = F₁(t) + F₂(t) and a priori unknown and independent emissions from components are given as F_i (t), i = 1, 2 at observed time instances t = e₁,..e_n . We assume that the light curves have the following Gaussian Process representation (Rasmussen 2004):

$$\begin{eqnarray}&&{F}_{1}(t)\sim G({\mu }_{1}(t),{k}_{1}(t,t))\end{eqnarray} \tag{ C19 }$$

$$\begin{eqnarray}&&{F}_{2}(t)\sim G({\mu }_{2}(t),{k}_{2}(t,t)\end{eqnarray} \tag{ C20 }$$

$$\begin{eqnarray}&&\begin{array}{l}{F}_{\mathrm{tot}}(t)={F}_{1}(t)+{F}_{2}(t)\sim G({\mu }_{1}(t)+{\mu }_{2}(t),{k}_{1}(t,t)\\ \,+{k}_{2}(t,t))\end{array}\end{eqnarray} \tag{ C21 }$$

where G is a Gaussian distribution with mean function ${\mathbb{E}}({F}_{i}(t))\,={\mu }_{i}(t),i\,=\,1,2$ and covariance functions ${Cov}({F}_{i}(t),{F}_{i}(t^{\prime} )\,={k}_{i}(t,t^{\prime} )$ (also called the kernels). Any number of components in multiple systems can be summed this way. In what follows, we distinguish between flux values at training locations t = e₁,..e_n and the flux values at some set of query locations ${t}^{\star }={e}_{1}^{\star },..{e}_{n}^{\star }$

The joint prior distribution over two unknown component light curves drawn independently from GP priors and their sum is given by Duvenaud et al. (2013):

$$\begin{eqnarray}&&\begin{array}{l}\left[\begin{array}{l}{F}_{1}(t)\\ {F}_{1}({t}^{\star })\\ {F}_{2}(t)\\ {F}_{2}({t}^{\star })\\ {F}_{1}(t)+{F}_{2}(t)\\ {F}_{1}({t}^{\star })+{F}_{2}({t}^{\star })\end{array}\right]\sim \\ \end{array}\begin{array}{l}G\left(\begin{array}{l}\left[\begin{array}{l}{\mu }_{1}\\ {\mu }_{1}^{\star }\\ {\mu }_{2}\\ {\mu }_{2}^{\star }\\ {\mu }_{1}+{\mu }_{2}\\ {\mu }_{1}^{\star }+{\mu }_{2}^{\star }\end{array}\right],\left[\begin{array}{lccccc}{\kappa }_{1} & {\kappa }_{1}^{\star } & 0 & 0 & {\kappa }_{1} & {\kappa }_{1}^{\star }\\ {\kappa }_{1}^{\star T} & {\kappa }_{1}^{\star \star } & 0 & 0 & {\kappa }_{1}^{\star } & {\kappa }_{1}^{\star \star }\\ 0 & 0 & {\kappa }_{2} & {\kappa }_{2}^{\star } & {\kappa }_{2} & {\kappa }_{2}^{\star }\\ 0 & 0 & {\kappa }_{2}^{\star T} & {\kappa }_{2}^{\star \star } & {\kappa }_{2}^{\star } & {\kappa }_{2}^{\star \star }\\ {\kappa }_{1} & {\kappa }_{1}^{\star T} & {\kappa }_{2} & {\kappa }_{2}^{\star T} & {\kappa }_{1}+{\kappa }_{2} & {\kappa }_{1}^{\star }+{\kappa }_{2}^{\star }\\ {\kappa }_{1}^{\star T} & {\kappa }_{1}^{\star \star } & {\kappa }_{2}^{\star T} & {\kappa }_{2}^{\star \star } & {\kappa }_{1}^{\star T}+{\kappa }_{2}^{\star T} & {\kappa }_{1}^{\star \star }+{\kappa }_{2}^{\star \star }\end{array}\right]\end{array}\right)\end{array}\end{eqnarray} \tag{ C22 }$$

where

$$\begin{eqnarray}&&{\kappa }_{i}={k}_{i}(t,t)\end{eqnarray} \tag{ C23 }$$

$$\begin{eqnarray}&&{\kappa }_{i}^{\star }={k}_{i}(t,{t}^{\star })\end{eqnarray} \tag{ C24 }$$

$$\begin{eqnarray}&&{\kappa }_{i}^{\star \star }={k}_{i}({t}^{\star },{t}^{\star }).\end{eqnarray} \tag{ C25 }$$

Using the formula for Gaussian conditional (see Equation (C17)–C18), the conditional distribution of unknown light curves (F_i , i = 1, 2) conditioned on their sum is given as:

$$\begin{eqnarray}&&\begin{array}{l}p({F}_{i}({t}^{\star })| {F}_{1}(t)+{F}_{2}(t))\\ \qquad =\,G\left({\mu }_{i}^{\star }+{\kappa }_{i}^{\star T}{\left({\kappa }_{1}+{\kappa }_{2}\right)}^{-1}({F}_{1}(t)+{F}_{2}(t)-{\mu }_{1}-{\mu }_{2}),\right.\\ \,\left.\qquad {\kappa }_{i}^{\star \star }-{\kappa }_{i}^{\star T}{\left({\kappa }_{1}+{\kappa }_{2}\right)}^{-1}){\kappa }_{i}^{\star }\right).\end{array}\end{eqnarray} \tag{ C26 }$$

From the above equation, the log-likelihood of the conditional distribution of ${F}_{2}^{\star }$ given F_tot at the query time instances t^⋆ is easily obtained:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal L }}_{\mathrm{GP}} & = & -\displaystyle \frac{1}{2}\left[\displaystyle \sum _{i=1}^{n}\mathrm{log}(2\pi )+\displaystyle \sum _{i=1}^{n}\mathrm{log}({\sigma }_{{F}_{2}| {F}_{\mathrm{tot}}(t)}^{2}({e}_{i}^{\star }))\right.\\ & & \left.+\displaystyle \sum _{i=1}^{n}\displaystyle \frac{{\left({F}_{\mathrm{tot}}({e}_{i})-{\mu }_{{F}_{2}| {F}_{\mathrm{tot}}}({e}_{i}^{\star })\right)}^{2}}{{\sigma }_{{F}_{2}| {F}_{\mathrm{tot}}(t)}^{2}({e}_{i}^{\star })}\right]\end{array}\end{eqnarray} \tag{ C27 }$$

where the number of observed points is n and the variance and mean of F₂(t^⋆) given F_tot(t) are

$$\begin{eqnarray*}&&{\sigma }_{{F}_{2}| {F}_{\mathrm{tot}}(t)}^{2}({e}_{i}^{\star })={\kappa }_{2}^{\star \star }-{\kappa }_{2}^{\star T}{\left({\kappa }_{1}+{\kappa }_{2}\right)}^{-1}{\kappa }_{2}^{\star }\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\mu }_{{F}_{2}| {F}_{\mathrm{tot}}}({e}_{i}^{\star })={\mu }_{2}^{\star }+{\kappa }_{2}^{\star T}{\left({\kappa }_{1}+{\kappa }_{2}\right)}^{-1}({F}_{\mathrm{tot}}({e}_{i})-{\mu }_{1}-{\mu }_{2}),\end{eqnarray*}$$

respectively. This log-likelihood of the component of the signal is integrating (marginalizing) over the possible configurations of the other components.