Relativistic binary-disc dynamics and the timing of OJ-287’s flares Open Access

ABSTRACT

We revisit the precessing black hole binary model, a candidate to explain the bizarre quasi-periodic optical flares in OJ-287’s light curve, from first principles. We deviate from existing work in three significant ways: (1) Including crucial aspects of relativistic dynamics related to the accretion disc’s gravitational moments. (2) Adopting a model-agnostic prescription for the disc’s density and scale height. (3) Using Markov chain Monte Carlo methods to recover reliable system parameters and uncertainties. We showcase our model’s predictive power by timing the 2019 Eddington flare within 40 h of the observed epoch, exclusively using data available prior to it. Additionally, we obtain a novel direct measurement of OJ-287’s disc mass and quadrupole moment exclusively from the optical flare timings. Our improved methodology can uncover previously unstated correlations in the parameter posteriors and patterns in the flare timing uncertainties. According to the model, the 26th optical flare is expected to occur on the 21st of August 2023 ± 32 d, shifted by approximately a year with respect to previous expectations.

1 INTRODUCTION

The precessing black hole binary model (PBM)¹ originally presented in Lehto & Valtonen (1996) has arguably been the most successful in explaining and predicting several unique features of OJ-287’s luminosity curve (Wolf 1916; Browne 1971; Kinman & Conklin 1971; Craine & Warner 1973; Corso et al. 1984). With the first observations dating to the late 1880s, this peculiar object is now classified as a blazar, situated at a redshift of z = 0.306 (Sitko & Junkkarinen 1985; Carangelo et al. 2003), and is thus composed by a supermassive black hole surrounded by an accretion disc powering a relativistic jet (Antonucci 1993; Urry & Padovani 1995; Ghisellini et al. 1998; Dunlop et al. 2003). Crucially, OJ 287’s light curve features bright, doubly peaked optical flares occurring with an approximate periodicity of twice every 12 yr, as well as a slower modulation on a time-scale of ∼60 yr (Valtaoja et al. 2000; Fan et al. 2010; Tang, Zhang & Pang 2014). Inspired by the original work by Sillanpaa et al. (1988), the PBM proposes a scenario in which the periodicity is explained by the presence of a smaller secondary black hole, orbiting the primary on a highly relativistic, inclined and eccentric trajectory that can be matched to the available data (see also Karas & Vokrouhlicky 1994). The sharp optical flares are then associated to impacts between the secondary and the disc, producing the characteristic twice per 12 yr orbital period structure. The 6-yr modulation is instead associated to the relativistic periastron advance time-scale of the binary, thus fixing the system’s characteristic mass and size to ∼10¹⁰ M_⊙ and ∼0.1 pc, respectively.

From its original inception, the PBM has since gone through several iterations and improvements, including, for example, a more sophisticated description of the disc’s response to the impacts (Valtonen et al. 2006a,b), the use of more detailed post-Newtonian (PN) equations of motion (EoM; Valtonen et al. 2010, 2016; Kacskovics & Vasúth 2022), as well as incorporating input from numerical simulations (Sundelius et al. 1997) and ulterior electromagnetic data (Yanny, Jannuzi & Impey 1997; Dey et al. 2018; Titarchuk, Seifina & Shrader 2023). The strength of this model is exemplified in the confirmation by the Spitzer Space Telescope of the so called ‘Eddington flare’ (Laine et al. 2020), within a truly remarkable 4 h of the predicted epoch: the 31st of July 2019 ut ± 4.4 h (Dey et al. 2018). More recently the model has been challenged by the works of Komossa et al. (2023a,b), based on both the alleged lack of a flare in 2022 October and the use luminosity scaling relations. This specific date of the 2022 flare is mentioned in the pre-print arXiv document (Valtonen et al. 2022), while the different prediction of July is given in both previous (Valtonen 2007) and later (Valtonen et al. 2023b) published works, based on a more sophisticated treatment of the accretion disc’s scale height.

Notwithstanding the claims in Valtonen et al. (2023a), the recent debate and the inherent complexity of OJ-287’s system (along with some interesting discrepancies also discussed in e.g. Komossa et al. 2023a,b) warrant to take a second look at alternative models for OJ 287’s light curve. Many of the ones proposed in the literature still invoke the presence of a secondary black hole (Katz 1997; Villata et al. 1998; Valtaoja et al. 2000), but do not require such highly relativistic initial conditions. Other forgo the requirement of a binary entirely, explaining the variability via complex jet beaming and precession effects (Villforth et al. 2010; Britzen et al. 2018; Qian 2018; Butuzova & Pushkarev 2020). All of the alternative models above are physically plausible and able to qualitatively explain many features of the blazar’s emission, including variability in bands other than optical (see Valtaoja et al. 2000, in particular). However, it is important to note that, as of today, only the PBM has been able to repeatedly predict the timing of optical flares (occasionally with spectacular precision, see e.g. Gupta et al. 2017; Dey et al. 2018; Laine et al. 2020), and that some of the many additional components of OJ-287’s light curve (see e.g. Valtaoja et al. 1985; de Diego & Kidger 1990; Pihajoki, Valtonen & Ciprini 2013a; Pihajoki et al. 2013b) can still be associated to a highly precessing binary.

In this work, we set out to revisit the PBM from first principles, deviating from existing methodology in three significant ways: First, we model the dynamical effects of the accretion disc’s gravitational moments and their PN cross-terms (see Section 2.2), a crucial element that was entirely missing in previous PBM iterations. Secondly, we adopt an agnostic description of OJ-287’s accretion disc based on density and scale height power laws, as opposed to assuming any specific model. Finally, we use Bayesian Markov chain Monte Carlo methods (MCMC, see Section 2.4) to recover reliable posterior distributions for the system’s parameters. As a consequence of our revised methodology, we obtain significantly diverging results from the established literature, which are discussed throughout Section 3.

2 METHODOLOGY

2.1 Basic setup and historical flare timings

In both the PBM and our work, the basic setup of OJ-287’s system consists in a primary black hole of mass M surrounded by a thin accretion disc, as illustrated in Fig. 1. A smaller secondary black hole of mass m orbits the primary on a highly eccentric, inclined orbit that intersects the disc twice every orbital period. The impact between the secondary and the disc deposits a large amount of energy into a localized region of gas (by means of relativistic Bondi accretion, Bondi 1952; Zanotti et al. 2011), which subsequently releases a flare of electromagnetic radiation. Given this basic picture, there are currently 11 detected optical flares in OJ 287’s light curve that can be clearly associated with impacts between the secondary and the disc (Dey et al. 2021; Valtonen et al. 2022, 2023b). Their epochs are listed together with the corresponding uncertainties in Table 1.

Over many refinements, the PBM has attempted to fit and predict the timings of the impact flares with higher and higher precision, eventually resulting in some extremely tight constraints for many parameters of the OJ-287 system, which we report in Table 2. Such tight bounds are made possible by the precision with which the historical flare timings have been observed (see Table 1), which occasionally amounts to a timing uncertainty below a few hours, or more typically ∼0.01 yr. In general, the observations vary in their precision as a consequence of the available instruments and the broadness of the flare’s luminosity curves, among other experimental constraints (Laine et al. 2020; Valtonen et al. 2021).

The goal of this work is to re-derive a precessing binary model for OJ-287’s system while also addressing some possible shortcomings of the PBM. The first ingredient is an extremely accurate description of the binary dynamics, collecting all contributions to the EoM that cause timing shifts comparable to the reported uncertainties (see Section 2.2). Then, we have to account for possible astrophysical mechanisms that can delay the emission of radiation after the impact occurs (see Section 2.3). Finally, we have to devise a strategy to efficiently explore a large parameter space of initial conditions, in order to fit the historical flare timings and recover parameters reliably (see Section 2.4). For the purposes of this work, we settle on a target accuracy of ∼0.01 yr regarding the dynamical evolution of the system, over a total integration time of approximately one century, close to the typical uncertainties in the data. Furthermore, we will show how the PBM is lacking of a proper description of the accretion disc’s gravitational potential, a modelling oversight that can induce timing shifts of order months to years, over the total integration time of approximately one century. A short list summarizing the differences between our model and the PBM cand be found in Table 3.

2.2 Binary, disc, and the EoM

2.2.1 PN binary dynamics

Taking the values reported in Dey et al. (2018) as a baseline, the black hole binary in OJ-287 can be characterized by a mass M ∼ 2 × 10¹⁰ M_⊙, a secondary to primary mass ratio q = m/M ∼ 0.01 and a semilatus rectum |$p \sim 20 \, r_{\rm S}$|⁠, where we defined the Schwarzschild radius of the primary r_S = 2GMc⁻². In order to precisely match the observed flare timings, we are required to model the dynamics of this highly relativistic binary with sufficient precision. As already noted in many PBM papers, a possibility is to use PN EoM. We can estimate the required PN order by the following simple considerations: corrections to Newtonian binary dynamics scale phenomenologically as powers of the dimensionless quantity r_S/p, i.e. the system’s Schwarzschild radius divided by the typical time-averaged orbital size (Einstein, Infeld & Hoffmann 1938; Blanchet 2014; Iorio & Zhang 2017; Maggiore 2018; Schäfer & Jaranowski 2018; Zwick et al. 2021). Over an integration time of ∼100 yr, it is therefore expected that PN corrections would produce a total accumulated shift in the flare timing of roughly 100/20ⁿ yr, where n is the order of the correction. Therefore, third-order PN corrections are required to achieve a timing accuracy of 0.01 yr. Beyond this simple estimate, we explicitly show that the higher order PN terms used in Dey et al. (2018) and Laine et al. (2020) are subdominant with respect to the modifications induced by the disc potential in Fig. 2. For the purposes of our work, we adopt the standard PN EoM for isolated binary black holes, which take a convenient form when expressed as an acceleration (see e.g. Blanchet 2014; Will & Maitra 2017; Bernard et al. 2018, for the explicit coefficients):

where v is the reduced mass’ velocity vector, n its unit vector while the coefficients |$\mathcal {A}$| and |$\mathcal {B}$| contain various PN modifications to the inverse square law. Of particular note are the 1.5 PN corrections required to model the primary’s spin (Barker & O’Connell 1975; Barker, O’Brien & O’Connell 1981), and the 2.5 PN corrections enforcing the the slow decay of orbital parameters caused by radiation reaction (Peters & Mathews 1963). The precession of the primary’s spin vector |$\vec{S}_1$| is mediated by the secondary’s mass (Barker et al. 1981; Kidder 1995; Faye, Blanchet & Buonanno 2006; Porto 2006). At lowest order in PN theory, it scales dimensionally as

It is therefore suppressed by a factor of ∼100 with respect to other first-order PN terms. The time-scale on which the primary’s spin evolves is thus approximately 100 times longer than the relativistic perihelion precession time-scale, which already amounts to ∼120 yr for the system parameters. Over the available observation time, the primary’s spin vector may thus at most precess by a few degrees. We show explicitly in Fig. 2 that the contributions to the timings caused by such small changes in spin orientation are subdominant with respect to the effect of the disc’s potential, and we therefore neglect them. Note that Dey et al. (2018) take a more sophisticated approach, modelling the orientation and the evolution of the central BH’s spin vector with PN equations. However, also note that Dey et al. (2018) may be overestimating the magnitude of high-order PN terms, in particular radiation reaction. In the latter work, it is stated that the 2.5 PN terms cause a change in the orbital period of approximately ∼10⁻³ per orbit. Taking the notorious evolution equation from Peters & Mathews (1963):

where q is the mass ratio and f(e) is the eccentricity enhancement function. Evaluating the equations for a mass of 1.8 × 10¹⁰ M_⊙, a mass ratio of q ∼ 100, a semimajor axis of |$\sim 50 \, r_{\rm S}$| and an eccentricity of ∼0.7, we find a period loss per orbit of ∼5 × 10⁻⁵, which is also is also consistent with the numerical integrations shown in Fig. 2.² Radiation reaction terms of even higher order will only weakly affect the dynamics of the system, as their strength is similarly suppressed by the small-mass ratio. According to the simple estimate detailed above, the 3.5-PN radiation reaction term would only contribute to a 20-min shift over an integration time of 100 yr, which is actually comparable to the reported approximately hour shift in Laine et al. (2020). Schematically, our binary EoM can be summarized as follows:

To fully specify these equations in an appropriate frame of reference we have to provide three parameters, i.e. the primary’s mass M, the secondary’s mass m and the primary’s dimensionless spin parameter ξ.

2.2.2 Disc multipoles and cross-terms

Newtonian and PN point mass forces are not the only relevant contributions to this system’s dynamics. According to the estimates in Valtonen et al. (2019), derived self-consistently in the PBM, the accretion disc surrounding OJ-287’s primary has a typical scale height of ∼10¹⁵ cm and a typical number density of ∼10¹⁴ cm⁻³. Taking the values mentioned above, we obtain a reference gas mass enclosed within |$\sim 100 \, r_{\rm {S}}$| of order several 10⁸ M_⊙, close to |$1~{{\ \rm per\ cent}}$| of the system’s total mass. Clearly, the gravitational influence of such a massive disc is an important ingredient to faithfully capture the dynamics of the binary, as it can potentially produce timing shifts of order 1 yr over the expected integration time. For the purposes of this work, we adopt a simple parametrized disc model characterized by a density profile ρ and a scale height profile h. Crucially, we do not assume that the disc structure is determined by the α-viscosity prescription (Shakura & Sunyaev 1973), as is customary in the PBM. Rather, we allow for the two profiles to be independent power laws:

where l_s is an arbitrary length-scale (e.g. the innermost stable circular orbit at 3r_S) and R is the radial distance to the primary. We model the disc’s potential as a mean field gravitational monopole (DM) and a mean field gravitational quadrupole (DQ). Together, the two capture the fundamental influence of the disc on the trajectory of the binary, i.e. an additional radial force as well as an axisymmetric perturbation. Additionally, we make the assumption that the primary’s spin vector is closely aligned with the disc’s symmetry axis, as expected for a black hole that has grown through gas accretion (Natarajan & Pringle 1998; King et al. 2005; Volonteri et al. 2005; Barausse 2012)

Decomposing the disc’s potential allows us to write down explicit analytical formulae for the resulting accelerations, rather than having to perform numerical integrations of the disc’s mass distribution. Crucially, the latter would dramatically slow down likelihood function evaluations in our numerical pipeline (see Section 2.4). Another option would be to use analytical disc potentials. However, only few density-potential pairs exist (Toomre 1964; Binney & Tremaine 1987; Kuzmin & Malasidze 1987), defeating the purpose of being as general as possible. Thus, we adopt EoM of the following form:

where we aligned they system’s symmetry axis with the z-direction, and e_z is a unit vector. Note that, for similar reasons as originally presented in Lehto & Valtonen (1996), we can neglect subdominant frictional and accretion forces acting on the secondary while it is submerged within the disc. The disc’s enclosed mass M_d and mass quadrupole Q₂ are defined as

Both of these quantities only depend on three different combinations of the density and scale height profile parameters. For convenience, we thus define an enclosed mass profile and a dimensionless quadrupole moment J₂:

where:

and j_eff = j₁ + j₂. Then, the three parameters M_d, j_eff, and J₂ completely specify the gravitational mean field monopole and quadrupole of the disc. Note that in general, the disc’s monopole contribution is not degenerate with the primary’s mass, because its influence on an eccentric test particle varies with the orbital phase. Note also that, in the limit of thin discs, the dimensionless quadrupole moment J₂ reduces to recognizable formula:

leading to the classic result of J₂ = −1/2 for a flat, thin homogeneous disc.

Our model for OJ-287’s accretion disc comes at a cost of one extra free parameter with respect to the the PBM, in which the disc is specified by a choice of the viscosity α and accretion rate |$\dot{M}$|⁠. While the latter choice certainly follows astrophysical expectations (at least as an effective model, see e.g. Gierliński & Done 2004; King, Pringle & Livio 2007; Kotko & Lasota 2012, often with values of α ∼ 0.1), relying on the α-disc model requires solving structural and thermal equations to obtain density and scale height profiles. In the context of the PBM, these equations have typically been separately solved beforehand, invoking a specific numerical model based on the work of Sakimoto & Coroniti (1981) and Stella & Rosner (1984). Considering how widely the behaviour of simulated accretion discs changes as a function of the chosen opacity tables or viscosity prescription (Abramowicz & Fragile 2013; Jiang, Davis & Stone 2016; Jiang & Blaes 2020; Chen et al. 2023), trading the computational and theoretical overhead of the PBM for a more general approach is of very clear benefit, despite the additional parameter. The downside of our choice is the impossibility to analytically model the rich dynamical interactions of the disc with the secondary, which include disc warping and scale height changes. The latter effects have been observed in hydrodynamical simulations of the system, and are required for the PBM to successfully match the flare timing in OJ-287’s light curve (Valtonen et al. 2023a,b). To our own surprise, we find that our model is still able to reproduce the correct timings optical flares without requiring a time advance prescription.

Finally, a consistent description of a highly relativistic PN system with monopole and quadrupole perturbations requires the inclusion of several additional cross-terms. The importance of PN cross-terms has been thoroughly discussed in Will (2014a, 2014b). In short, they are fundamental to assure that the system’s energy is properly conserved over relativistic perihelion advance time-scales, a factor that is crucial when one hopes to track the timing of consecutive flares. Thus, we must add both monopole and quadrupole PN cross-terms in our EoM, up to first PN order (see Will 2014a,b, for the explicit formulae).

2.2.3 EoM integration and plane intersections

To summarize, we adopt a set of EoM of the following schematic form:

Relative to the Newtonian contribution, we expect conservative PN accelerations to be smaller by a factor of roughly ∼20ⁿ on average over an orbit, where n is the PN order. The 2.5-PN contribution is further suppressed by a factor of ∼100 due to the mass ratio, making it de facto the least important term in the adopted EoM. In addition to point mass terms, the disc monopole and quadrupole both can contribute at |$\sim 1~{{\ \rm per\ cent}}$| relative level, given the baseline values in Valtonen et al. (2019). The PN cross-terms are suppressed by a further factor of ∼20, contributing at the |$\sim 0.1~{{\ \rm per\ cent}}$| level. Note that these rough expectations are confirmed a posteriori throughout our numerical integrations, an example of which is shown in Fig. 2. Crucially, both disc and cross-term contributions to the EoM had previously been neglected in the PBM, despite being orders of magnitude larger than even the 2.5-PN radiation reaction force.

Given the EoM and a set of initial conditions, our first goal is to efficiently determine the time at which the secondary black hole impacts the disc, or equivalently intersects the z = 0 plane. To achieve this, we integrate the EoM with the eighth-order Runge–Kutta implementation ‘DOP853’ of the scipy solve_ivp package (Virtanen et al. 2020) and extract the z component of the secondary’s position vector. To find the epochs of the impacts, we search for all local maxima of the function |1/z| with the scipy find_peaks function, further refining the impact times with the in-built interpolation feature ‘dense_output’ of solve_ivp. We perform several resolution studies, finding that a choice of the integrator parameter r_tol = 10⁻⁵ and a temporal interpolation resolution of 120 h assures that the numerical error in the impact timings remains below ∼1 h for all reasonable initial conditions. A single run of the integration and interpolation pipeline requires approximately 0.8 s of computation time, as opposed to several minutes when the same accuracy is requested without interpolation.

2.3 Time delays

2.3.1 Disc delay

As a result of the impact of the secondary with the disc, a hot bubble of gas forms, expands and releases a flare of thermal radiation after a certain time delay. For the purposes of this work, we settle on the original delay prescription used in Lehto & Valtonen (1996), derived from theoretical principles tailored to the radiation-dominated inner regions of accretion discs. As shown more recently, the physics upon which the prescription is based are also able to reproduce the luminosity, spectrum and duration of the optical impact flares (Pihajoki 2016; Dyba, Mach & Malec 2019), and may also be a plausible explanation for quasi-periodic emissions (Franchini et al. 2023). Note, however, that several different prescriptions for such delays have been experimented with throughout the early iterations of the PBM (Lehto & Valtonen 1996; Ivanov, Igumenshchev & Novikov 1998). A thorough and interesting discussion on the effects of varying the delay prescription can be found in Pietilä (1998). The disc time delay, τ_d, depends on both the properties of the disc at the impact site, as well as the properties of the impactor:

where the parameter δ describes the impactor’s orientation and velocity relative to a Keplerian disc (Pihajoki 2016). Similarly to Dey et al. (2018),³ we preserve the physical scaling laws of equation (21), but allow it to be modified by a constant proportionality pre-factor f_err, which is fitted as part of the orbital solution. Allowing for this freedom is important, since the disc’s response to impacts is the most complex hydro and thermodynamical ingredient of both the PBM and our model. As such, it is most likely not fully captured by a simple analytical formula. To summarize, our time delay prescription τ_d depends linearly on f_err, in addition to scaling as power laws with the general disc parameters M_d, j_eff and J₂. Contrast this with the PBM, where the delays are effectively a function of the viscosity and accretion rate parameters that determine a solution of a specific accretion disc model.

An additional difference between our implementation and late iterations of the PBM are the latter’s use of ‘time advances’, which model the local warping of the disc towards the approaching secondary. In, for example, Dey et al. (2018), these time advances are pre-computed by means of numerical simulations and are reported to amount to weeks and even months. The secondary black hole can have a strong local influence on the trajectory of a fluid element only when its gravitational pull is comparable to the primary’s. For an impact at apoapsis, the gravitational interaction time between a fluid element and the ∼100 times lighter secondary will therefore be of the order of |$\sim (1/2) (50\, r_{\rm {S}}/\sqrt{100})/(c/\sqrt{50}) \sim 20$| days, found by equating gravitational forces and dividing by typical orbital speeds at |$\sim 50 \, r_{\rm {S}}$|⁠. By integrating the secondary’s gravitational acceleration for an interaction time, the total out-of-plane displacement of the fluid element would amount to only a few au, changing the impact timing by less than a day. Large time advances must therefore be caused by more complex tidal interactions between the binary and the disc, which have been typically modelled via hydrodynamical simulations in the PBM literature. In the latter framework, they are considered an essential part of a precise timing model. In this work, we have found that the timing of optical flares can be fit with good precision without requiring any modelling of disc warping or time advances. It is plausible to attribute this major difference to our implementation of the disc’s potential, which had been missing from previous iterations of the PBM. We do note that dynamical variations in the disc’s scale height caused by tidal interactions are certainly physical and may be an improvement to our model’s accuracy. However, we do not wish to rely on numerical simulations of specific accretion discs, nor do our estimates suggest that they are an important contribution to the timing of flares at the target precision of ∼0.01 yr. Most importantly, our analysis will also show that such advances are not required for our timing model to be predictive (see Section 3.3).

2.3.2 Relativistic delays

Several relativistic effects will influence the arrival times of radiation flares from massive systems at large distances (Karas & Vokrouhlicky 1994; Dai, Fuerst & Blandford 2010). The first and most obvious is cosmological redshift, which simply stretches the total period by a factor 1 + z = 1.306 (Lemaître 1931). In this work, we also implement the effect of Shapiro delays, τ_S, which produce impact radius-dependent shifts of order weeks for our baseline parameters (Shapiro 1964) and have also been neglected in the PBM:

where D is the distance corresponding to OJ-287’s redshift. Römer time-delays caused by the orientation of the disc could, in principle, contribute up to ∼50 d if the disc were edge-on, the impact occurred at apoapsis and the orbit’s major axis were also aligned with the line of sight. Considering, however, that the disc is expected to be only ∼4° off from a face-on configuration (Lehto & Valtonen 1996), the resulting shifts are suppressed by a small geometric factor, at most ∼sin (4°) ∼ 0.07. Thus, they are typically close to or smaller than our target accuracy of ∼0.01 yr. For the purposes of this work, we decide to neglect Römer delays in favour of having one less free parameter, which will allow us to fit our model to a subset of historical flare timings and ‘re-predict’ the 2019 July flare (see Section 3.3) without risking to overfit the data.

2.4 Fitting to the observed flare timings

2.4.1 Parameter space and initial conditions

To summarize, our updated flare timing model is determined by seven parameters, three of which characterize the binary and four of which characterize accretion disc and its response to impacts (see Table 2). Combined, these parameters represent a necessary compromise between the complexity required to properly capture the physics of OJ-287’s system and the requirement to reduce the total number of degrees of freedom. In addition, we must specify the secondary’s trajectory via several initial conditions. Since the timing of impacts is not affected by rotating the system along its symmetry axis, we can place the secondary in the mid-plane of the disc at the position (−x₀, 0, 0) at an arbitrary initial time t₀. We choose to define the secondary’s instantaneous velocity vector via a combination of its magnitude v₀, inclination Ω₀ and tilt ι₀ with respect to the z-axis (see Fig. 1 and Table 2). In total, a set of candidate flare timings thus depend on seven system parameters and four initial conditions. The historical data provides us with 11 measurements (see Table 1). Additionally, we assume that there were indeed 19 flares between 1912 and 2019, some of which went undetected. This consideration adds an additional constraint that essentially fixes the orbital period, excluding some extreme orbital solutions that would account for large gaps between flares. Our model is therefore overdetermined, and we have the additional opportunity to only use a subset of the data, containing, for example, 10 flares (see Section 3.3).

We must now devise a strategy to sample the 11-dimensional parameter space of our model and minimize the residuals between the candidate and the observed flare timings. In the PBM, the authors adopt what is essentially a manual minimization routine, in which individual model parameters are selected a priori to roughly reproduce the data and satisfy basic astrophysical expectations. Further parameters are then adjusted one by one iteratively until every single flare occurs within its observed timing window (Lehto & Valtonen 1996; Valtonen et al. 2006a; Dey et al. 2018; Valtonen et al. 2022). While this trial method can certainly produce valid solutions, it does not constitute a systematic search and gives no guarantee that the recovered parameters truly represent a global minimum. Most importantly, the resulting parameter uncertainties do not account for cross-correlations and degeneracies. They are therefore most likely vastly under-reported, even providing the validity of the underlying model.

2.4.2 MCMC sampling and numerical pipeline

MCMC sampling provides a framework in which both of the issues mentioned above are addressed naturally. For the purposes of this work, we use the MCMC software emcee (Foreman-Mackey et al. 2013) which supports efficient parallelization via the multiprocess package (McKerns et al. 2012). We adopt a standard likelihood function |$\mathcal {L}$| of the form:

where we attempt to minimize the least-squares difference between some trial flare timings and the historical data, while consistently accounting for the various data gaps. Additionally, we enforce the correct total number of flares by returning a vanishing likelihood when the appropriate condition is not met. To recover the parameter posteriors, we initially run 64 MCMC walkers on wide uninformed priors reported in Table 2. After ∼15 000 iterations, the walkers have thoroughly scouted parameter space and identified all likelihood maxima. In order to re-sample, the most interesting part of the posterior, we then re-initialize 32 walkers in a small neighbourhood of the global likelihood maximum, running further MCMC chains until the Gelman–Rubin criterion is met (Gelman & Rubin 1992; Vats & Knudson 2018) and the best-fitting values are converged within their standard deviations. Convergence typically requires approximately ∼30 000 iterations.

For the purposes of this work, we ran the numerical pipeline with a data set containing all 11 historical flare timings, forming the basis of our results in Sections 3.1 and 3.2. Additionally, we re-ran it on a data set which excludes the 2019 July flare, in order to confirm the predictive power of our model (see Section 3.3). In total, we evaluated just under 4000 000 trial orbits,⁴ for a total of ∼850 cpu h. Thanks to the magic of parallelization, the whole computation was performed on an 8-core Lenovo laptop over the course of several weeks.

3 RESULTS

3.1 Binary and disc parameter posteriors

3.1.1 General features and correlations

The results of our MCMC run, including all 11 historical flare timings, are a set of posterior distributions for the 7 parameters and 4 initial conditions of our model. They are visualized in Fig. 3 as two corner plots (Foreman-Mackey 2016). The best-fitting values and the uncertainties resulting from the marginalized posteriors are listed in Table 4 and compared to the most up-to-date results of the PBM. In general, the recovered binary and disc parameters qualitatively reproduce the broad expectations set by the framework of the PBM. However, they do differ in several important details. Here, we comment on the ones that shed most light on the physics of the system.

First, there is a clear bi-modality in the posterior for the initial inclination Ω₀. The two likelihood peaks are separated by around 3° around a mean value of ∼189|${_{.}^{\circ}}$|5. The common PBM assumption that the orbit be perfectly perpendicular to the disc plane is excluded, a constraint that arises from properly modelling the effects of a preferential symmetry plane, i.e. the disc’s gravitational quadrupole, on the flare timing. The bi-modality in Ω₀ is also responsible for widening the posteriors of many other parameters, including the primary mass M. This leads to factor 10 larger uncertainty with respect to the PBM result, despite the almost identical recovered best-fitting value of 1.834 × 10¹⁰ M_⊙. Additionally, we recover a significantly lower value for the primary’s spin with respect to the PBM. Its posterior distribution shows a clear correlation with the disc’s scale mass, which arises due to an interesting interaction with the disc’s monopole gravitational moment. The latter forces a precession-like effect at apoapsis, where the enclosed mass is large, rather than at periapsis where relativistic effects (such as spin induced frame dragging) are large. Furthermore, the correlation significantly broadens the range of likely spin values, leading to an uncertainty of |$\sim 50~{{\ \rm per\ cent}}$| of the best-fitting value rather than the remarkable (and perhaps overoptimistic) |$\sim 1~{{\ \rm per\ cent}}$| reported in the PBM. The best-fitting value for the secondary mass m is 0.82 × 10⁸ M_⊙ approximately half of the reported result in the PBM. This large discrepancy is related to the introduction of the disc’s quadrupole moment and there are hints of a correlation between the parameters m and J₂. Physically, this behaviour arises from small offsets in the system’s centre of mass being degenerate with spurious quadrupolar contributions to the potential.

3.1.2 Dynamical constraints on OJ-287’s accretion disc

Provided that our EoM are correct and that our delay prescriptions are accurate, our work provides a novel measurement of a quasar accretion disc’s mass profile and quadrupole moment that is not based on specific disc models, luminosity scaling relations or spectral data. We recover a disc mass scale M_d at |$100 \, r_{\rm S}$| of 1.2 × 10⁸ M_⊙, approximately five times lower than the best-fitting values in the PBM (which we extrapolate from Valtonen et al. 2019). Within the sampled radii of ∼10–50 r_S, the disc’s enclosed mass grows approximately linearly, scaling as |$r^{2-j_{\rm {eff}}}$|⁠, where the best-fitting value is given by j_eff = 0.92 ± 0.12. As mentioned previously, in the limit of thin discs the quadrupole moment J₂ and j_eff are related by a simple formula, given in equation (16). In our analysis, this turns out (a posteriori) to be an extremely good approximation, confirming that OJ-287’s accretion disc can be indeed modelled as thin, at least gravitationally. The individual parameters j₁ and j₂ describing the density and scale height profiles are therefore only very poorly constrained. This suggests that in effect our disc model only truly depends on two independent free parameters rather than three.

Interestingly, the flare timing data selects a disc with profiles roughly following the expectations of the PBM α-disc, despite our model exclusively relying on dynamical information. This can bee seen as a strong validation of using such a disc in the first place, provided that the underlying model is correct. Indeed, the recovered disc in Valtonen et al. (2019) is shown to have an approximately constant scale height and a decreasing density, from which we can estimate values of j_eff ∼ 0.6 and J₂ ∼ 0.4. Given our setup, we cannot directly estimate an accretion rate nor a viscosity from our recovered parameters. However, assuming for simplicity the same viscous time-scales as in the PBM, the factor 5 reduction of the total mass budget more easily aligns with the luminosity and spectral constraints highlighted recently in Komossa et al. (2023a,b). This additionally showcases the advantage of not forcing the model to fit an α-disc, for which the aforementioned observational constraints are inconsistent with such a large primary mass. Finally, the disc delay coefficient takes a best-fitting value of ∼0.02, leading to delays of order days to months, also compatible with previous expectations from the PBM. Despite the many similarities, it is important to stress that our disc model is fundamentally different from the one used in the PBM, as it does not allow for warps and local scale height variations that cause timing advances. The fact that, with our model, the flare timings may be fitted without requiring such elements is a further indication of the importance of modelling the accretion disc’s gravitational potential. However, the assumption of a rigid disc is certainly a simplification, and is likely reflected in the best-fitting posterior parameters. Our disc is best understood as being an effective model, which suffices for the purposes of timing the optical flares. Determining whether this disc model can explain all other electromagnetic signatures of OJ-287 is significantly beyond the scope of this work, but an interesting avenue for further investigation.

For visualization purposes, a to-scale rendition of the best-fitting binary and disc components comprising of OJ-287’s system is shown in Fig. 3, along with the trajectory of the secondary over the last century.

3.2 Predictions and uncertainties for future flares

In stark contrast to the most recent PBM claims, our model predicts that the next flare, i.e. the 26th in the customary numbering system, should occur on the 21st of August 2023 ± 32 d. The expected date is thus shifted by by ∼13 months with respect to (Valtonen et al. 2023b). This large discrepancy is entirely caused by the addition of the gravitational influence of the disc in the EoM and the resulting differences in the recovered parameter posteriors. Indeed, comparing the best-fitting trajectory in Fig. 3 with, for example, fig. 3 in Dey et al. (2018), our model forecasts that the location of the 26th impact will also differ drastically from previous expectations, being much closer to the site of the 2015 impact. If detected, the upcoming flare is therefore likely to have similar characteristics to the well-studied 2015 Centenary flare when it comes to its luminosity, duration, and spectrum (Ciprini et al. 2015; Komossa et al. 2015; Shappee et al. 2015; Valtonen et al. 2016, 2019). Beyond the epoch itself, the uncertainty in our prediction is much larger than the typical timing uncertainties that have been reported in the PBM. While this is partially explained by the wider posteriors in our recovered parameters (see Table 4), further investigation reveals some interesting details.

In the top panel of Fig. 4, we construct the probability distribution function (PDF) of the 26th flare by drawing 500 posterior samples and collecting the timing results in a histogram. There is clear evidence that the multimodality present in the parameter posteriors (specifically the inclination Ω₀) is carried over in the timing PDF, which shows a primary peak accompanied by at least an additional local maximum shifted by approximately –30 d. In similar fashion, we compute the standard deviation in the timing of several past and future flares, and plot the results in the bottom panel of Fig. 4. We uncover a very clear alternating pattern, in which each consecutive flare is characterized by an unusually large or small uncertainty. Indeed, the uncertainty for the 27th flare on 2031 July 16 is only a few days, as is also shown in the top panel of Fig. 4. Furthermore, the pattern presents a modulation on a time-scale of roughly 60 yr, clearly associated to the system’s relativistic perihelion precession. The latter causes the system’s orientation to rotate fully in approximately 120 yr, forcing the trajectory of the secondary to extrude primarily above or below the accretion disc’s plane for roughly half as much time. The alternating pattern in the timing uncertainty is also associated to the obit’s orientation with respect to the disc. By inspecting the secondary’s trajectory, we observe that impacts for which the orbital nodes are aligned with the minor axis are characterized by either large (months) or small uncertainties (days). Conversely, impacts for which the nodes are aligned with the major axis typically present average uncertainties (weeks). Thus, the exact timing of the upcoming flare unfortunately happens to be the most uncertain of the last century, while the prospective 2031 July 16 flare is much more easily constrained, having an uncertainty of only ∼40 h. Note that both the timing PDF multimodality and the alternating pattern in the uncertainties are features that are clearly associated with the real physics of the system. They could only be uncovered with a sophisticated sampling of the posterior distributions, which in this work has been accomplished via MCMC methods.

We note that recent investigations of OJ-287’s light curve have revealed evidence for the secondary’s impact with the disc in early 2022, in the form of a blue flare (Valtonen et al. 2023a). If it is indeed produced by an impact, the observation of such a flare clearly invalidates both the results of our parameter recovery and our predictions for future flares. It does not, however, invalidate the requirement of any precessing binary model to properly include the gravitational potential of the accretion disc, which can produce timing shifts of order 1 yr. Furthermore, we have shown how a sophisticated sampling technique, for example, MCMC sampling, is required to truly understand the intricacies of the timings. We look forward to incorporating additional information, for example, from additional flares, and revisiting elements of our models (listed in Table 3) in future iterations of this work.

It is also the case that, at the time of revising this manuscript (2023 September 22), Earth-based tracking of OJ-287’s light curve has recommenced after the yearly period (early July to early September) in which the Sun’s position affects observations. An up-to-date visualization of the system’s light curve is curated online.⁵ It can be seen that observations via the KRK instrument of the Jagellonian University recommence in the early days of September. They show hints of a preceding luminosity maxima, somewhere around mid-August, that may have occurred during the observational gap. However, with the given data, it is clearly impossible to determine whether an optical flare may or may not have actually taken place in the period predicted by our work, at least by naive analyisis. This is unfortunate, since the timing of the 26th flare sets a clear distinction between our work and the predictions of the PBM. Nevertheless, it is interesting in its own right that both models instead give very similar results for the previous optical flare, despite the very different dynamical and methodological ingredients.

3.3 Mock prediction of the 2019 Eddington flare

As discussed thoroughly in Section 2.4, our model only requires 10 historical flare timings to be completely determined. Therefore, we are able run our numerical pipeline on a subset of the historical timing data that does not contain the 2019 July flare, in order to test whether our model would have been able to predict the latter with comparable accuracy to Dey et al. (2018) and Laine et al. (2020). Once again, the result of our numerical runs is a set of parameter posteriors, differing only slightly from what is shown in Fig. 3. While uncertainties are generally increased due to more pronounced correlations, the new posteriors still preserve the patterns in the timing uncertainties discussed in Section 3.2. Thus, conditioning our model to this mock data set leads to a relatively well-constrained ‘mock prediction’ for the 2019 July flare. Fig. 5 shows the PDF of the 25th impact flare, computed by collecting the timing results from 500 draws of the new posterior distributions. According to the PDF, our model would have been able to predict the Eddington flare within 40 h of the actual confirmed detection (Laine et al. 2020). The timing PDF is peaked around 29th of July ± 33 h and overlaps with the actual detection epoch.

4 CONCLUSION

In this paper, we have constructed an alternative precessing binary model to explain the quasi-periodic optical flares in OJ-287’s light curve. In contrast to the model first proposed by Lehto & Valtonen (1996) (and then iterated upon over many years), our work consistently accounts for the gravitational influence of the accretion disc by decomposing it into a monopole and a quadrupole contribution. Furthermore, we have shown how sophisticated sampling techniques are required to uncover realistic correlations and uncertainties in both parameter posteriors and flare timings. The results of our work are discussed throughout Section 3, and compared with the established literature. However, it is important to highlight the predictions of our model which may differ from previous expectations:

• Our model can reproduce the timing of the 2019 July flare within 40 h of it’s detection, when conditioned with timing data only available prior to it.

• In stark contrast with established literature, our model predicts that the 26th optical flare should have occurred on the 21st of August 2023 ± 32 d. Unfortunately, the majority of this range lies in the period when OJ-287 is unobservable to Earth-based telescopes.

• The location of the 26th impact is extremely close to the previous impact responsible for the 2015 Centenary flare. The flare is therefore likely to have similar luminosity, duration and spectral properties.

The crucial insight, however, is that despite the very different ingredients, both the PBM and our model have (or would have) been able to successfully predict the timing of the great 2019 Eddington flare. The inclusion of the accretion disc’s potential makes the binary model predictive without requiring the complex numerical modelling of disc and its response to the perturber as required in the PBM. However, from the pattern revealed in the timing uncertainties (see Section 3.2 and Fig. 4), we have shown that the epoch of the 25th flare is incredibly robust to changes in the system’s parameters. By analogy, it must also be robust to changes in the chosen underlying model. This clearly indicates how the simpler framework for a precessing black hole binary adopted in this work is able to explain and re-predict past flare timings in OJ-287’s light curve. Additionally, it validates the prospect of using flare timing data to directly measure the gravitational potential of quasar accretion discs (see Section 3.1), revealing their structure and composition in an unprecedented way. What is perhaps more difficult to interpret is the difference in the predicted epoch for the 26th flare. Within the assumptions of our model, including the dynamical influence of the accretion disc causes a timing shift of almost a full year with respect to previous expectation. Interestingly, the shift occurs despite an almost identical recovered primary black hole mass. Clearly, only a thorough investigation of the differences between the PBM and our own model can help shed light on the truly essential elements of the precessing binary framework. Another possibility is that the characteristics of OJ-287’s light curve are conclusively shown not to be associated with a binary system. Then, it would be time to revise and refine alternative frameworks, such that they may become predictive rather than only descriptive. To be able to conclusively answer these questions, we would like to strongly encourage the observational community to keep tracking this fascinating system over the months of August and September and the upcoming years.

ACKNOWLEDGEMENTS

The authors acknowledge support from the Swiss National Science Foundation under Grant 200020_192092. LZ acknowledges Pedro R. Capelo, Mudit Garg, and Andrea Derdzinski for several helpful discussions. LZ acknowledges the institution and participants of GALSTAR UZH.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the authors.

Footnotes

References