Authenticating the Presence of a Relativistic Massive Black Hole Binary in OJ 287 Using Its General Relativity Centenary Flare: Improved Orbital Parameters

Our approach to determine the parameters of the BBH central engine model for OJ 287 proceeds as follows. First, an approximate orbit of the secondary BH is calculated by numerically integrating the above-mentioned PN-accurate equations of motion (Equation (1)) while using some trial values of the independent parameters. This orbit produces a list of reference times at which the secondary crosses the y = 0 plane of the accretion disk (see Figure 3). However, these plane-crossing epochs are not the same as the observed outburst times. We need to take into consideration certain astrophysical processes that occur during the time interval between the BH impact and the observed optical outburst epoch. The effects of these processes are incorporated by adding a "time delay" to the plane-crossing times. These delays represent the time interval between the actual creation of a hot bubble of gas due to the disk impact and when it becomes optically thin and releases a strong burst of optical radiation. An additional temporal correction is required to model the fact that when the secondary black hole approaches the accretion disk, the disk as a whole is pulled toward the secondary. This ensures that the secondary BH impact occurs before it reaches the accretion-disk plane of the primary black hole, depicted by the y = 0 line in Figure 3. Therefore, we subtract a time interval, termed "time advance," from the plane-crossing time. This leads to a new list of corrected reference times.

Zoom In Zoom Out Reset image size

Let us digress briefly to introduce the way we model the time delay and the accretion disk of the primary BH. We use the accretion-disk model detailed by Lehto & Valtonen (1996), which is based on the α_g disk model of Sakimoto & Coroniti (1981), with scaling provided by Stella & Rosner (1984). In this model, the disk impacts are followed by thermal flares after a time delay t_d given by

$$\begin{eqnarray}&&{t}_{d}=d\,{{m}_{2}}^{1.24}\,{{v}_{\mathrm{rel}}}^{-4.23}\,{h}^{-0.29}\,{{\rm{\Sigma }}}^{0.91},\end{eqnarray} \tag{ 7 }$$

where the delay parameter d is a proportionality factor to be determined as part of the orbit solution. This also applies to the disk thickness h and the secondary BH mass m₂. The impact velocity of the secondary relative to the disk (v_rel) is known in the model for each impact and the fiducial values are those given by Lehto & Valtonen (1996). Furthermore, the scaling for the disk surface density Σ is

$$\begin{eqnarray}&&{\rm{\Sigma }}\approx {{\alpha }_{g}}^{-0.8}\,{\dot{m}}^{0.6},\end{eqnarray} \tag{ 8 }$$

where α_g is the viscosity coefficient and $\dot{m}$ is the mass accretion rate in Eddington units. Typical particle number density in the accretion disk in our model is $\sim {10}^{14}\,{\mathrm{cm}}^{-3}$ (see Table 2 of Lehto & Valtonen 1996 for detailed astrophysical properties of the disk). The value of α_g depends on the unbeamed total luminosity of OJ 287: α_g = 0.1, 0.3, 1.0 for the total luminosity of 2.5 × 10⁴⁵, 1.2 × 10⁴⁶, 5 × 10⁴⁶ erg cm⁻² s⁻¹, respectively. Since the observed (beamed) luminosity is ∼10⁴⁷ erg cm⁻² s⁻¹, the most likely α_g value is near the lowest quoted value, α_g ≈ 0.1, since the relativistic Doppler boosting factor is likely to be in excess of δ ≈ 20 (Worrall et al. 1982). Interestingly, the orbit determination provides a ${\alpha }_{g}-\dot{m}$ correlation as a side result while determining the orbit from impact flare timings. However, it is not possible to extract these two parameters individually, since the time delay is practically a function of $\dot{m}/{\alpha }_{g}$ and depends weakly on either parameter.

We now move to work on the above-mentioned corrected reference times. It is customary to normalize the list so that the 1983 outburst has the exact time of 1982.964. This is done by subtracting the difference between the 1983 corrected reference time, namely the "disk-crossing time plus time delay minus time advance," and the actual 1983 outburst time (1982.964), from all other reference times. Thereafter, we check the timing of a certain outburst, typically that of the 1973 outburst, by adjusting usually the initial orientation of the major axis of the binary. We pursue new trial solutions until the 1973 outburst time is within the observed time interval (1972.935 ± 0.012).

In the next step, the disk thickness parameter (h) is found by requiring that the 2005 outburst timing matches with the observations (2005.745 ± 0.015). This process is repeated until we determine all nine independent parameters of the BBH system. In other words, the procedure involves adjusting each parameter of the BBH central engine model so that some particular outburst happens within a certain time window. Further details of the solving procedure are described by Valtonen (2007). At each stage of the iterative procedure, we ensure that the previous conditions are still satisfied; if not, the procedures are repeated. When all the outburst timings, listed in Table 1, match within the listed uncertainties, we regard that solution as an acceptable solution.

We performed 1000 trials for orbit solutions and at each time, as expected, with slightly different initial parameter values. It turns out that 285 cases converged to an acceptable solution, but the remaining trials were interrupted as the procedure exceeded the preset number of attempts while varying a parameter. The general experience with the code is that the convergence is not always found even if the trial is continued much longer. The average values of the parameters are listed in Table 2, as well as the 1σ scatter of these values as the uncertainty. The independent parameters of Table 2 are the two masses m₁ and m₂, the primary BH Kerr parameter χ₁, the apocenter eccentricity e₀, the angle of orientation of the semimajor axis of the orbit Θ₀ in 1856 (the starting year of the orbit calculation), the precession rate of the major axis per period ΔΦ, and an ambiguity parameter γ_obs that we employ to incorporate the effects of dominant-order hereditary contributions to GW emission in the BBH dynamics, as evident in Equation (5) (we use the subscript obs to distinguish the observationally determined γ from its GR-based estimate). Additionally, two independent parameters incorporate the effects of astrophysical processes that are associated with the accretion disk impact of the secondary BH. These are listed in Table 2 as d and h, where d is the delay parameter present in Equation (7), while the disk thickness parameter h is a scale factor with respect to the "standard" model of Lehto & Valtonen (1996). In other words, the average half thickness of the accretion disk is ∼3 × 10¹⁵ cm and we need to multiply the disk thickness, given in Table 2 of Lehto & Valtonen (1996), by disk thickness parameter h to get the actual thickness profile of the disk in our model. Note that these two parameters (d and h) are functions of the mass accretion rate $\dot{m}$ and the viscosity parameter α_g of the standard α_g accretion disk.

Table 2.  Independent and Dependent Parameters of the BBH System in OJ 287 According to our Orbit Solution

	Parameter	Value	Unit	Error
(1)	(2)	(3)	(4)	(5)
	m1	18348	${10}^{6}\,{M}_{\odot }$	±7.92
	m2	150.13	106 M⊙	±0.25
	χ1	0.381		±0.004
	h	0.900		±0.001
Independent	d	0.776		±0.004
	ΔΦ	38.62	deg	±0.01
	Θ0	55.42	deg	±0.17
	e0	0.657		±0.001
	γobs	1.304		±0.008
Derived	${P}_{\mathrm{orb}}^{2017}$	12.062	year	±0.007
	${\dot{P}}_{\mathrm{orb}}$	0.00099		±0.00006

Note. Note that γ_obs provides the observationally determined value for the γ parameter, invoked to incorporate heuristically the effect of dominant-order "tail" contributions to GW emission on Equation (5) for BBH dynamics.

Download table as:  ASCIITypeset image

Table 2 also lists certain derived parameters that characterize the BBH in OJ 287—namely, the present (redshifted) orbital period ${P}_{\mathrm{orb}}^{2017}$ and its rate of decrease due to the emission of GWs (${\dot{P}}_{\mathrm{orb}}$). We find the rate of orbital period shrinkage to be ∼10⁻³; in contrast, the measured ${\dot{P}}_{\mathrm{orb}}$ values for relativistic binary pulsar systems like PSR J1913+16 are ∼10⁻¹² (Wex 2014). This is roughly nine orders of magnitude smaller than in OJ 287, which demonstrates the strong-field relativistic nature of OJ 287. To probe the relevance of higher-order RR terms in Equation (1), we repeat the above detailed orbital fitting procedure while employing only the dominant 2.5PN order contributions to $\ddot{{\boldsymbol{x}}}$. This resulted in ${\dot{P}}_{\mathrm{orb}}=0.00106$, which indicates that the higher-order RR contributions reduce the quadrupolar-order GW flux by about 6.5%.

We demonstrate the predictive power of our BBH central engine model for OJ 287 with the help of Table 3, which lists all the epochs associated with past impacts as well as future impacts within the years 1886–2056 according to our model. The entries of Table 3 quantify many facets of our model. Column 1 provides the starting times of the outbursts (t_out) in a Julian year (J2000.0), while t_del indicates the time delay between the impact of the secondary BH with the accretion disk and the starting of the outburst. The listed t_del values differ from those of Lehto & Valtonen (1996) by the scale parameter d. The time advance (t_adv) arises from the bending of the accretion disk due to the presence of the secondary BH prior to the impact. The radial distance (R_imp) of the secondary and its orbital speed (v₀) at various impact flare epochs are also listed in Table 3. In a later section, we will explain the importance of the next predicted outburst.

We are now in a position to explain why we were forced to introduce the γ parameter and obtain an estimate for it from the impact flare timings. Recall that the prediction and analysis of the GR centenary flare observations were pursued using fully 3PN-accurate orbital dynamics that incorporated the effects of dominant (Newtonian) order GW emission on BBH dynamics. Therefore, it is natural to extend the PN accuracy of our BBH dynamics by including the 3.5PN order contributions that are available in Iyer & Will (1995). However, it is not advisable to add only the 3.5PN order reactive contributions to BBH dynamics. This is because the fully 3.5PN order BBH equations of motion can provide extremely inaccurate secular orbital phase evolution during the inspiral regime of unequal mass compact binaries. The above statement is fully endorsed by the second column entries in Table 1 of Blanchet et al. (1995b), which showed that the accumulated number of GW cycles (${ \mathcal N }$) due to 1PN and 1.5PN tail contributions to GW emission tend to cancel each other for large mass-ratio binaries. A close inspection of the above table also reveals that the ${ \mathcal N }$ estimate, based on fully 1PN-accurate orbital frequency evolution ($\dot{\omega }$), substantially increases the expected number of GW cycles for the orbital revolutions of unequal mass BH-NS binaries. Recall that 1PN-accurate orbital frequency evolution corresponds to 3.5PN accurate orbital evolution in the PN description. These considerations are important for the BBH orbital evolution in OJ 287 as the fully 3.5PN accurate equations of motion can provide erroneous orbital phase evolution. An erroneous orbital phase evolution ensures that the observed impact flares' timings will not be consistent with the BBH central engine model. Indeed, we have verified that the use of fully 3.5PN accurate equations of motion leads to a loss of acceptable solutions, and the situation is not improved by adding the 4.5PN order contributions. Therefore, it is crucial for us to incorporate the effect of hereditary contributions to GW emission on our BBH dynamics that appear at 4PN order in Equation (1). This is also influenced by the earlier mentioned fact that ${ \mathcal N }$ estimates from 1PN contributions to $\dot{\omega }$ get essentially canceled by the 1.5PN "tail" contributions to $\dot{\omega }$ for unequal mass binaries. For this reason, it is rather important for us in our Equation (1) to incorporate terms that can lead to 1.5PN accurate tail contributions in $\dot{\omega }$.

Therefore, we introduce a heuristic way of incorporating the effect of dominant-order tail contributions to GW emission into BBH orbital dynamics. This is essentially due to the nonavailability of 4PN(tail) contributions in the form of Equation (3). The plan, as noted earlier, is to introduce an ambiguity parameter γ at the dominant-order RR terms such that the second line of Equation (1) can be replaced using Equation (5). Note that we can now introduce an additional parameter while describing the dynamics of the BBH in our model as we have, at our disposal, the tenth fixed point from the timing of the 2015 November impact flare. Indeed, the results we list in Tables 2 and 3 are obtained by such a prescription for the BBH dynamics.

Clearly, our procedure for evolving the BBH orbit requires further scrutiny. It is natural to ask if additional higher order RR contributions to $\ddot{{\boldsymbol{x}}}$ are required while evolving the BBH orbit in OJ 287. We gather from various numerical experiments, associated with the above detailed orbit-fitting procedure, that the velocity-dependent terms in Equation (1) are crucial for incorporating the effects of GW emission on the dynamics of the BBH system in OJ 287. A plot of appropriately scaled and dimensionless B coefficients that appear at 2.5PN, 3.5PN, 4PN, and 4.5PN orders in Equation (3) is displayed in Figure 4. The visible linear regression suggests that the further higher-order contributions to GW emission would not substantially influence the orbital evolution of the BBH in OJ 287. Note that the contribution appearing at 4PN order in Figure 4 is obtained by multiplying the scaled dimensionless 2.5PN order B coefficient by 0.304, which arises by subtracting unity from the our γ_obs = 1.304 estimate.

Zoom In Zoom Out Reset image size

In the next section, we show that the extracted γ_obs value is consistent with the general relativistic orbital phase evolution of the BBH in OJ 287 that explicitly incorporates the effects of dominant-order tail contributions to GW emission.

In a GW phasing approach, we are interested in an accurate evolution of the BBH orbital phase that takes into account the effects of the conservative and reactive terms present in first and second line of Equation (1), respectively. This approach accurately incorporates the effects of the first two lines of Equation (1) without explicitly invoking them. In our implementation, we employ a 3PN-accurate Keplerian-type parametric solution to model the conservative parts of the orbital dynamics, i.e., to model the first line of Equation (1) (Memmesheimer et al. 2004). This allows us to express the temporally evolving BBH orbital phase ϕ as

$$\begin{eqnarray}&&\phi -{\phi }_{0}=(1+k)\,l,\end{eqnarray} \tag{ 9 }$$

where l is the mean anomaly defined to be l = 2 π (t − t₀)/P_b (Memmesheimer et al. 2004; Königsdörffer & Gopakumar 2006). Initial values of the orbital phase and its associated coordinate time are denoted by ϕ₀ and t₀, while P_b stands for the radial orbital period of a PN-accurate eccentric orbit. Furthermore, the dimensionless fractional periastron advance per orbit k ensures that we are dealing with a precessing eccentric orbit. A close comparison with Königsdörffer & Gopakumar (2006) reveals that we have neglected 3PN-accurate periodic contributions to the orbital phase in the above equation. This is justified as we are focusing our attention on the secular evolution of the BBH orbital phase. The fractional rate of periastron advance k is an explicit function of n and e_t where n = 2 π/P_b is the mean motion and e_t provides the eccentricity parameter that enters the PN-accurate Kepler equation of Memmesheimer et al. (2004). The 3PN-accurate expression for k is given in the Appendix. Note that it is by using the 3PN-accurate expression for k that we are incorporating the BBH orbital phase evolution at the 3PN-accurate conservative level into $\ddot{{\boldsymbol{x}}}$.

The next step requires us to model the influences of the reactive terms (second line of Equation (1)) on the conservative 3PN-accurate orbital phase evolution. This is done by providing differential equations that describe temporal evolutions for n and e_t due to emission of GWs, consistent with the "reactive" PN orders of Equation (1); details are provided by Königsdörffer & Gopakumar (2006). Following Blanchet et al. (1995a), it is convenient to split the fully 2PN-accurate differential equations for n and e_t into two parts, given by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{dn}}{{dl}}\right)}^{2\mathrm{PN}}={\left(\displaystyle \frac{{dn}}{{dl}}\right)}^{\mathrm{Inst}}+{\left(\displaystyle \frac{{dn}}{{dl}}\right)}^{\mathrm{Tail}},\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}^{2\mathrm{PN}}={\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}^{\mathrm{Inst}}+{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}^{\mathrm{Tail}},\end{eqnarray} \tag{ 10b }$$

where we used n dt = dl to obtain (dn/dl, de_t/dl) expressions from their (dn/dt, de_t/dt) counterparts of Königsdörffer & Gopakumar (2006). We call those contributions that depend only on the state of the binary at the usual retarded instant ${T}_{{\rm{r}}}$ as the "instantaneous" contributions and refer those terms that are a priori sensitive to the BBH dynamics at all previous instants to T_r as the "tail" (or hereditary) terms. Moreover, these instantaneous contributions appear at the "absolute" 2.5PN, 3.5PN, and 4.5PN orders like the reactive contributions in Equation (1) for $\ddot{{\boldsymbol{x}}}$, while the "tail" contribution enters at the "absolute" 4PN order. Therefore, the differential equations for the evolution of n and e_t may also be symbolically written as

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{\left(\displaystyle \frac{dn}{dl}\right)}_{\,2{\rm{P}}{\rm{N}}}^{{\rm{e}}{\rm{x}}{\rm{a}}{\rm{c}}{\rm{t}}} & = & {\left(\displaystyle \frac{dn}{dl}\right)}_{2.5{\rm{P}}{\rm{N}}}+{\left(\displaystyle \frac{dn}{dl}\right)}_{3.5{\rm{P}}{\rm{N}}}\\ & & +{\left(\displaystyle \frac{dn}{dl}\right)}_{4.5{\rm{P}}{\rm{N}}}+{\left(\displaystyle \frac{dn}{dl}\right)}_{4{\rm{P}}{\rm{N}}},\end{array}\end{array}\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}\begin{array}{rcl}{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{2\mathrm{PN}}^{\mathrm{exact}} & = & {\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{2.5\mathrm{PN}}+{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{3.5\mathrm{PN}}\\ & & +{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{4.5\mathrm{PN}}\,+{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{4\mathrm{PN}}.\end{array}\end{eqnarray} \tag{ 11b }$$

The explicit expressions for these 2PN-accurate contributions are listed in the Appendix.

A few comments are in order. Note that these orbital-averaged equations are derived by computing time derivatives of 2PN-accurate n and e_t expressions in the modified harmonic gauge, available in Memmesheimer et al. (2004). We apply the heuristic arguments that the time derivatives of the binary orbital energy and angular momentum should be balanced by their far-zone GW energy and angular momentum fluxes. Close inspection of Königsdörffer & Gopakumar (2006) and the above equations reveals that there exists no strict one-to-one correspondence between different PN orders present in the second line of Equation (1) and the above equations for n and e_t. In other words, 3.5PN-order contributions to the differential equations arise from both 2.5PN and 3.5PN terms in Equation (1), and this will be relevant for us while estimating the value γ_GR. Furthermore, it is the use of certain rational functions of e_t that allowed us to write closed-form expressions for the tail contributions to dn/dl and de_t/dl, as explained by Tanay et al. (2016).

We can now obtain secular orbital phase evolution of the BBH in OJ 287 due to the action of fully 2PN-accurate reactive dynamics on the fully 3PN-accurate conservative dynamics. It should be noted that this is what the first and second-line contributions in Equation (1) aim to provide while implementing our BBH central engine model, provided we ignore periodic contributions to its solution. We obtain the desired orbital phase evolution by numerically imposing on $\phi -{\phi }_{0}\,=(1+k)l$, given by Equation (9), the secular evolutions in n(l) and e_t(l) due to Equations 11(a) and (b). The resulting phase evolution is treated as ϕ_exact (as we use exact 2PN-accurate equations for dn/dl and de_t/dl) in our effort to obtain γ_GR associated with the BBH in OJ 287. Let us emphasize that the use of fully 2PN-accurate differential equations for n(l) and e_t(l) ensures that the orbital phase evolution incorporates the effects of the dominant-order hereditary contributions to the GW emission.

To obtain γ_GR, we repeat the above procedure while employing slightly different differential equations for n and e_t that do not contain the tail contributions. This is expected as we are trying to model the BBH orbital phase evolution defined by the first line of Equation (1), while the second-line contributions are replaced by Equation (5) to model the reactive contributions to $\ddot{{\boldsymbol{x}}}$. In other words, the instantaneous contributions to dn/dl and de_t/dl are modified to include the effect of the γ factor, present in Equation (5). The resulting Newtonian (quadrupolar or 2.5PN) contributions to dn/dl and de_t/dl are multiplied by the γ_GR factor. However, the 3.5PN-order instantaneous contributions now consist of two parts. The first part contains γ_GR as a common factor and arises from the 2.5PN terms in $\ddot{{\boldsymbol{x}}}$. The second part is independent of γ_GR. The resulting differential equations for n and e_t can be written as

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{{dn}}{{dl}}\right)}_{\ 2\mathrm{PN}}^{\mathrm{approx}}={\gamma }_{\mathrm{GR}}{\left(\displaystyle \frac{{dn}}{{dl}}\right)}_{2.5\mathrm{PN}}\\ +\,{\gamma }_{\mathrm{GR}}{\left(\displaystyle \frac{{dn}}{{dl}}\right)}_{3.5\mathrm{PN}}^{(2.5\mathrm{PN})}+{\left(\displaystyle \frac{{dn}}{{dl}}\right)}_{3.5\mathrm{PN}}+{\left(\displaystyle \frac{{dn}}{{dl}}\right)}_{4.5\mathrm{PN}},\end{array}\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{2\mathrm{PN}}^{\mathrm{approx}}={\gamma }_{\mathrm{GR}}{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{2.5\mathrm{PN}}\\ +\,{\gamma }_{\mathrm{GR}}{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{3.5\mathrm{PN}}^{(2.5\mathrm{PN})}+{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{3.5\mathrm{PN}}+{\left(\displaystyle \frac{{{de}}_{t}}{{dl}}\right)}_{4.5\mathrm{PN}},\end{array}\end{eqnarray} \tag{ 12b }$$

where the 2.5PN, 3.5PN, and 4.5PN terms stand for the relative Newtonian, 1PN, and 2PN contributions due to GW emission, and the influence of the tail contribution needs to be captured by certain values of γ for the BBH in OJ 287. The appearance of γ_GR at 2.5PN and 3.5PN orders in Equation (12) is due to the introduction of γ_GR at the 2.5PN order in Equation (1). This is because the time derivatives of PN-accurate orbital energy and angular momentum using Equation (5) are required to obtain Equation (12) as detailed in Königsdörffer & Gopakumar (2006). We do not list explicitly these contributions, and as expected these contributions add to the instantaneous contributions in dn/dl and de_t/dl, given by Equations 11(a) and (b), when we let γ = 1. However, it is straightforward to modify Equations (28)–(33) of Königsdörffer & Gopakumar (2006) to obtain γ versions of Equations (32) of Königsdörffer & Gopakumar (2006), and they provide the 2.5PN and 3.5PN contributions to our Equation (12). As to the 4.5PN contributions to the above-listed approximate equations for dn/dt and de_t/dt, we employed 2PN contributions in Equations (B2) and (B3) of Königsdörffer & Gopakumar (2006).

We are now in a position to repeat the numerical procedure that gave us ϕ_exact while employing Equations 12(a) and (b) for dn/dl and de_t/dl. The resulting accumulated orbital phase evolution is termed ϕ_approx (as we are using an approximate formula for dn/dl and de_t/dl that involves γ_GR). To obtain γ_GR, we demand that ${\rm{\Delta }}\phi ={\phi }_{\mathrm{exact}}-{\phi }_{\mathrm{approx}}$ should be smaller than ${({\rm{\Delta }}\phi )}_{\mathrm{tol}}$: an observationally extracted tolerable value for the accumulated orbital phase for the BBH in OJ 287. This is justified as uncertainties in the observed outburst timings constrain the accuracy with which we can track the orbital phase evolution of the BBH in OJ 287. An estimate for (Δϕ)_tol is obtained by noting that at best the uncertainty in the observed outburst timing is roughly 0.0005 year. In the BBH model, the associated minimum ϕ uncertainty occurs at the apocenter as the secondary moves slowly there. Invoking a Keplerian orbit and the expression for dϕ/dt at the apocenter, we write

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{d\phi }{{dt}} & = & \displaystyle \frac{2\pi }{{P}_{\mathrm{orb}}}\displaystyle \frac{\sqrt{1-{e}^{2}}}{{(1+e)}^{2}}\\ & = & \displaystyle \frac{2\pi \times 0.263}{{P}_{\mathrm{orb}}}.\end{array}\end{eqnarray*}$$

This leads to an observationally relevant (Δϕ)_tol estimate as

$$\begin{eqnarray}&&{({\rm{\Delta }}\phi )}_{\mathrm{tolerance}}=\displaystyle \frac{0.0005\times 2\pi \times 0.263}{9.385}=0.000088.\end{eqnarray} \tag{ 13 }$$

We have checked that the inclusion of the PN corrections to dϕ/dt does not substantially vary the above estimate for (Δϕ)_tol.

To estimate γ_GR, we equate the earlier estimated difference in the accumulated BBH orbital phase Δϕ to (Δϕ)_tol while considering roughly 12 orbital cycles of the BBH in OJ 287. This number of orbital cycles roughly corresponds to the span of observational data on OJ 287 (∼130 years) that we used to determine γ_obs from impact flare timings. In practice, we let (Δϕ)_tol  vary between −10⁻⁴ and +10⁻⁴ while varying γ_GR during (Δϕ) estimations. We infer that (Δϕ)_tol  forces the γ_GR estimates to be 1.2917 ± 0.0045. This is a very encouraging result as our GR-based estimate for γ (γ_GR) is fairly close to our earlier estimate γ_obs, listed in Table 2, that was purely based on the observed impact flare timings while employing Equation (5) to model the reactive contributions to $\ddot{{\boldsymbol{x}}}$.

We have verified that the inclusion of the spin–orbit contributions to the PN-accurate orbital phase evolution does not affect the above estimate. In Figure 5 we plot Δϕ(l) as a function of l for different γ_GR values to show its secular variations. These plots show that the differences between our two (exact and approximate) estimates for the accumulated BBH orbital phase remain within the (Δϕ)_tol  values while varying γ_GR between 1.287 and 1.296. Note that the quantity "${\gamma }_{\mathrm{GR}}-1$" provides the present comparative strength of the dominant-order tail contributions to BBH dynamics in comparison with the quadrupolar order gravitational RR terms that appear at 2.5PN order. The expected GW emission-induced hardening of the BBH orbit ensures that the effects of higher order contributions such as the tail terms can be more prominent during later epochs, which implies that the value of γ will be epoch dependent. The present detailed analysis opens up the possibility of evolving the BBH in OJ 287 with the help of Equations (1) and (5) while using the GR-based γ estimate obtained above. In the next section, we explore the observational benefits of such an approach after taking a careful look at the light curves of OJ 287 in the vicinity of impact flare epochs.

Zoom In Zoom Out Reset image size

This subsection mostly deals with the extended historical data sets on OJ 287 near impact flare epochs, predicted from our BBH central engine model.

The availability of optical data sets on OJ 287 extends to the late 19th century because of the fact that OJ 287 lies close to the ecliptic, and consequently was often unintentionally photographed in the past. New historical data points have been found by searching photographic plate archives for images containing OJ 287. The main data set used in this work was compiled by Milan Basta and one of the authors (H.L.,http://altamira.asu.cas.cz/iblwg/data/oj287) in 2006, using the previous compilations by other authors (A.S., L.O.T., T.P.). For the last 12 years our main database has been compiled by K.N. and S.Z. One of us (P.P.) made a light-curve compilation in 2012 including much of the data up to that point in time (Pihajoki et al. 2013; Valtonen & Pihajoki 2013). For the historical part, there has been a huge increase of data from the photographic plates of the Harvard plate collection, studied by R.H. He was able to evaluate numerous (almost 600) additional HCO plates not used before, partly because the object brightness was close to the plate magnitude limit, and he provided 364 additional measurements and 209 upper limits. Some earlier historical data were published by Hudec et al. (2013) (HCO data) and (2001) (Sonneberg Observatory data).

In Figure 1 we display the summary of the present state of available observational optical data points on OJ 287. The figure includes unpublished photographic data measurements obtained by R.H. from various astronomical photographic archives. The collection of the photometric data for the OJ 287/2015 campaign is described by Valtonen et al. (2016). S.Z. has collected the data and harmonized them in a uniform system.

As noted earlier, we are mostly interested in observational data sets around a number of impact flare epochs, predicted in our model. Let us emphasize that many of these data sets are not employed while determining the parameters of our BBH central engine model for OJ 287. This influenced us to find possible signatures of impact flares in the historical data sets on OJ 287. For this purpose, we compare sections of observed data points on OJ 287 around the predicted impact flare epochs with a template light curve and search for the presence of possible patterns. The light curve from late 2015 to early 2017 is used as the standard template to compare with less-complete data sets from earlier major flare epochs that were created by secondary BH impacts according to our model. This is what we pursue in Figures 6–10. For smooth comparisons, we shift the 2015 outburst light curve, shown as line plots, backward in time by certain amounts such that the start times of outbursts coincide with the epochs listed in Table 3. In these figures, observational data sets are usually marked by points. For many epochs we only have upper limits for the brightness of OJ 287; these data points are marked by the tips of red thin arrows pointing downward.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Influenced by Valtonen et al. (2011b), we introduce a certain "speed factor" (s.f.) parameter that indicates how fast the earlier outburst has proceeded in comparison with the 2015 outburst. Termed the f-parameter by Valtonen et al. (2011b), it depends on the velocity of the secondary BH and the distance of the impact site from the primary. It turned out that the rate at which the outburst takes place is about three times faster for impacts near the pericenter than for impacts at a larger distance from the primary (Lehto & Valtonen 1996). This is an important aspect, quantified by the dynamical timescale t_dyn of Table 3 in Lehto & Valtonen (1996), while making detailed comparisons of our template flare light curve with actual observational data points. We have also shifted the template light curve vertically by a certain magnitude (mentioned in the plot legends at the top-right corner of each plot) for matching the base levels of two outbursts. The variations in base levels arise because of long-term variations in the optical light curve as is evident from Figure 1.

We begin by displaying in Figures 6 and 7 our comparisons of the 2015 template light curve with the well-documented outbursts of Table 1. To quantify these comparisons, we compute Pearson correlation coefficients between the 2015 and the earlier outburst data sets, implemented using the Correlation routine of  Mathematica (Wolfram Research 2018). These coefficients are computed for restricted data sets that span two months and their values are listed below the panels. The bottom panels in Figure 7 require further explanation. We do not compute the correlation coefficient for the 1947 outburst, as we do not have sufficient data points to calculate it (this is despite the fact that the crucial sudden rise in brightness is well covered by observations during 1947). In panel (d) of Figure 7, we display the expected light curve for the 2019 July outburst, and this is obtained by moving the 2007 light-curve data points forward in time. We invoked the data set of the 2007 flare, as the predicted 2019 outburst is expected to be a periastron flare like the 2007 one in our BBH central engine model. An additional point worth noting is the low correlation coefficient of 0.65 between the 2007 and 2015 outburst data sets; this is mainly due to the smaller rise in brightness of the first peak during the 2007 outburst as visible in panel (b) of Figure 7.

In Figures 8 and 9, we compare a few less-well-covered outburst data sets with the 2015 light curve. For the 1906 and 1945 outbursts shown, respectively, in panels (b) and (d) of Figure 8, the upper limits provide good constraints on outburst timings. However, we required a shifted 2007 impact flare light curve to obtain a visual match with very few data points of the 1898 outburst, as shown in panel (a) of Figure 8. We employed the 2007 light curve as our standard outburst light curve; the 2015 data set turned out to be not adequately long, as we compare data sets spanning roughly two years in panel (a) of Figure 8 (a more detailed study of the 1898 outburst is available in Hudec et al. 2013).

We now show comparisons for the 1959, 1964, and 1971 outburst data sets in panels (a), (b), and (c) of Figure 9. Visual inspection reveals that impact flares most likely occurred at these predicted epochs. Unfortunately, we do not have archival data points for the primary peaks of these impact outbursts. However, available data points from neighborhood epochs are consistent with our 2015 light curve. It is reasonable to infer from these figures that the available data sets are consistent with 17 outburst epochs, and this is seven more than what is necessary to describe the BBH orbit. Finally, in Figure 10, we pursue comparisons of the remaining 6 outbursts that are listed in Table 3. Unfortunately, we have fairly poor observational data associated with the relevant epochs. Visual inspections suggest that the data points are not inconsistent with the model light curve. However, the available data points are too far from the relevant primary peaks to make any meaningful timing estimates.

A closer look at Figures 6 and 7 reveals that the outbursts with sufficient observed data points give high correlation coefficients with our template light curve. This clearly requires us to invoke the speed factor (s.f.). The time evolution of the bubble of gas that is released as a result of the impact of the secondary on the accretion disk depends on the local disk conditions, the thickness of the disk, and the internal sound speed of the released bubble of gas; these quantities are encapsulated in the dynamical timescale or the speed factor (s.f.). Once the speed factor is included, the listed high correlations suggest the possibility of predicting the general shape of the optical light curves associated with the future impact outbursts.

It should be noted that the major axis of the OJ 287 BBH eccentric orbit happens to lie in the disk plane during the disk crossing associated with the 2015 outburst (Valtonen et al. 2016). This indicates that the BBH orbit is symmetric while going backward and forward from its 2013 configuration. Therefore, the 2019 impact is expected to occur in a manner similar to the 2007 impact with the same speed and at the same impact angle. This should result in an optical light curve that has a high resemblance to the well-documented 2007 impact flare light curve. These are our main arguments behind panel (d) of Figure 7, where we are essentially predicting the shape of the light curve for the 2019 impact outburst. Extending these arguments, we may state that the 2022 impact event should resemble the 2005 impact outburst, while the 2031 impact light curve should resemble the one associated with the 1995 event in OJ 287, and so on. It applies also to the other future impact events, listed in Table 3. These considerations open up interesting possibilities during the next outburst, and this is tackled in the next subsection.

This subsection revisits an idea that was explored in detail by some of us earlier (Valtonen et al. 2011b): testing a formulation of the BH no-hair theorem that requires that the dimensionless quadrupole moment of a BH (q₂) should fulfill the relation q₂ = −q χ² where q = 1 in GR (Thorne 1980). For obtaining observationally relevant constraints on q, we invoke our GR-based estimate for γ, namely γ = 1.2917, and treat the above q as a free parameter while numerically determining the BBH orbit. Recall that the q parameter enters the BBH dynamics via the classical spin–orbit coupling term ${\ddot{{\boldsymbol{x}}}}_{Q}$ in the equations of motion.

It turns out that the 2019 outburst time is highly correlated with the q parameter. We plot the correlation of this "no-hair" parameter q against the time of the rapid rise of the flare (t₂₀₁₉) in the left panel of Figure 11, where the slanted lines show the expected correlation and the 1σ deviation from it. For this numerical experiment, we employed the BBH parameters listed in Table 2. Additionally, we extracted a numerical fit that provides the accuracy with which we can estimate q in terms of t₂₀₁₉:

$$\begin{eqnarray}&&q=1.0-311.857({t}_{2019}-2019.569).\end{eqnarray} \tag{ 14 }$$

This formula clearly shows that an accurate t₂₀₁₉ measurement should allow us to determine the q parameter below the 10% level. Indeed, this provides a substantial improvement over the present q estimate that confines it only in the 0.5–1.5 range. Currently, there exist no other observational constraints on q even at this rather relaxed accuracy level.

Zoom In Zoom Out Reset image size

We now turn our attention to the feasibility of observing the next predicted outburst. In our BBH central engine model, the next outburst is expected to peak on 2019 July 31. However, observing OJ 287 during the expected 2019 outburst window is practically impossible from the ground: the angular distance in the sky between the Sun and OJ 287 is only ∼6° at the beginning of this event, and it goes down to 4° by the time of peak brightness. Unfortunately, objects at small angular distances from the Sun are difficult to observe even from a satellite in Earth's orbit, owing to the high background caused by intense sunlight.

In any case, it will be useful to have good light curves for the 2018–2019 and 2019–2020 observing seasons. We provide the expected light curve of OJ 287 that spans a few weeks around 2019 July 31, in panel (b) of Figure 11. Especially important would be to determine the epoch of rapid flux rise associated with the primary peak of the 2019 impact flare. However, even if the primary peak is missed, the secondary peak can still provide some information on the flare timing. This was the case in observations of the 1994 outburst. We infer from panel (d) of Figure 9 that the primary peak of OJ 287 was missed owing to the closeness of OJ 287 to the Sun during the 1994 outburst. However, the visual correlation of the secondary flares with our 2015 template light curve allows us to state that the starting time of the outburst was around 1994.59, and this is consistent with our BBH model. Similarly, it will be useful to organize an observational campaign to observe the secondary and subsequent peaks using ground-based facilities around the 2019 impact flare epoch.

The flare of 2019 July 31 would be best observed from a satellite observatory that is far from Earth, such as the STEREO-A solar telescope or the Spitzer Space Telescope. As an observational target, OJ 287 is relatively easy, as it brightens to ∼13 mag in the optical R band. The only problem is to find a space telescope that is able to do optical photometry during the rapid rise of flux, on 2019 July 29–31!