###### Abstract

Successful observations of the seven predicted bremsstrahlung flares from the
unique bright blazar OJ 287 firmly point to the presence of a nanohertz
gravitational wave (GW) emitting supermassive black hole (SMBH) binary central
engine. We present arguments for the continued monitoring of the source in
several electromagnetic windows to firmly establish various details of the
SMBH binary central engine description for OJ 287. In this article, we explore
what more can be known about this system, particularly with regard to
accretion and outflows from its two accretion disks. We mainly concentrate on
the expected impact of the secondary black hole on the disk of the primary on
December 3, 2021 and the resulting electromagnetic signals in the following
years. We also predict the times of exceptional fades, and outline their
usefulness in the study of the host galaxy. A spectral survey has been carried
out, and spectral lines from the secondary were searched for but not found.
The jet of the secondary has been studied and proposals to discover it in
future VLBI observations are mentioned. In conclusion, the binary black hole
model explains a large number of observations of different kinds in OJ 287.
Carefully timed future observations will be able to provide further details of
its central engine. Such multi-wavelength and multidisciplinary efforts will
be required to pursue multi-messenger nanohertz GW astronomy with OJ 287 in
the coming decades.

###### keywords:

BL Lacertae objects: individual: OJ 287; quasars: supermassive black holes;
accretion, accretion discs; gravitational waves; galaxies: jets

1 1 0 https://doi.org/ Promise of persistent multi-messenger astronomy with
the blazar OJ 287 Promise of persistent multi-messenger astronomy with the
blazar OJ 287 Mauri J. Valtonen 1,2, Lankeswar Dey 3, A. Gopakumar3, Staszek
Zola4, S. Komossa5, Tapio Pursimo6, Jose L. Gomez7, Rene Hudec8,9, Helen
Jermak10 and Andrei V. Berdyugin2 Mauri Valtonen, Lankeswar Dey, A. Gopakumar,
Staszek Zola, S. Komossa, Tapio Pursimo, Jose L. Gomez, Rene Hudec, Helen
Jermak and Andrei Berdyugin Valtonen, M.; Dey, L.; Gopakumar, A.; et al.
Correspondence<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Supermassive black hole (SMBH) binary systems are expected to be common in the
Universe Begelman et al. (1980); Valtaoja et al. (1989); Mikkola & Valtonen
(1992); Valtonen (1996); Quinlan (1996); Milosavljevic & Merritt (2001);
Volonteri et al. (2003). Inspiral gravitational waves from SMBH binaries with
total mass $\gtrsim 10^{9}\,M_{\odot}$ and orbital periods between months and
years should lie in the nanohertz (nHz) - microhertz range (Sesana et al.,
2004; Burke-Spolaor et al., 2018, 2018). The routine detection of nHz
gravitational waves (GWs) from SMBH binary systems is the primary goal of
rapidly maturing Pulsar Timing Array (PTA) experiments Arzoumanian et al.
(2020); Goncharov et al. (2021); Chen et al. (2021). The eventual detection of
nHz GW sources will complement the vibrant field of GW astronomy, inaugurated
by the routine detections of hectohertz GWs from stellar-mass black hole (BH)
binaries by the LIGO-Virgo collaboration Abbott et al. (2021); The LIGO
Scientific Collaboration et al. (2021). This collaboration also heralded the
era of multi-messenger GW astronomy by the observations of hectohertz GWs from
a binary neutron star merger (GW170817) and its electromagnetic (EM)
counterparts (EM170817) Abbott et al. (2017, 2017). The eventual detection of
continuous nHz GWs from individual SMBHBs and observations of their
electromagnetic (EM) counterparts can make profound contributions to the
emerging field of multi-messenger GW astronomy Burke-Spolaor et al. (2018);
Xin et al. (2021). EM observations in recent years have suggested the
existence of many candidate binary SMBH systems with varying orbital periods
(e.g., OJ 287 with 12 yr period (Sillanpää et al., 1988); 3C 66B with 1 yr
period in the radio (Sudou et al., 2003, but also see Jenet et al. (2004);
Iguchi et al. (2010)); SDSS J1201+3003 identified from a tidal disruption
event optical lightcurve Liu et al. (2014); PG1302-102 identified in the
optical Catalina transient survey (Graham et al., 2015, see also Zhu & Thrane
(2020)); and see, e.g., Charisi et al. (2016) who found 33 further quasars
with statistically significant periodic variability). Some of these systems
and candidates are in the nHz GW regime which should be detectable in GWs by
the PTA consortium during the SKA era (Zhu et al., 2019; Dey et al., 2019;
Feng et al., 2020). However, none of these systems had as detailed modeling
and testable predictions as OJ 287, and no lightcurve has had the very long-
term coverage that was obtained for OJ 287.

Persistent multi-wavelength electromagnetic observations will be crucial to
pursue multi-messenger GW astronomy with OJ 287 in the coming decades. In what
follows, we provide a brief introduction to the binary black hole (BBH)
central engine description for OJ 287 and discuss the context and the need to
pursue a variety of electromagnetic observations. These efforts should allow
us to place tight observational constraints on the various aspects of the BBH
central engine description for OJ 287 that will be crucial for pursuing multi-
messenger GW astronomy with OJ 287. We now describe briefly our BBH
description for OJ 287 and ongoing multi-wavelength observational campaigns on
the source.

###  1.1 OJ 287 and its BH binary central engine description

The bright blazar OJ 287 (RA: 08:54:48.87, Dec: +20:06:30.6), situated at a
redshift of $z=0.306$, was first discovered as a radio source during the
course of the Ohio sky survey in the 1960s (Dickel et al., 1967). However, the
optical light curve of OJ 287 goes all the way back to the year 1888 as it had
been unintentionally photographed often in the past due to its proximity to
the ecliptic plane (Sillanpää et al., 1988; Hudec et al., 2013).
Interestingly, the optical light curve of OJ 287, extending over $130$ years,
exhibits unique quasi-periodic high-brightness flares (outbursts) with a
period of $\sim 12$ years (Sillanpää et al., 1988; Valtonen et al., 2006; Dey
et al., 2018). Further, the blazar shows a long-term variation in its apparent
magnitude with a period of $\sim 60$ years (Valtonen et al., 2006) (see Figure
1).

In 1982, the binary nature of OJ 287 was recognized by Aimo Sillanpää, while
constructing historical light curves for the quasars in the Tuorla - Metsähovi
variability survey which had begun two years earlier (Kidger, 2007). It was
obvious that the next major outburst was about to happen and the blazar
monitoring community was alerted. This first OJ 287 campaign was fruitful
(Smith et al., 1985; Sillanpää et al., 1985). A more detailed explanation of
the binary BH (BBH) scenario was subsequently proposed to explain the quasi-
periodic flares in OJ 287 by associating them with increased accretion flow to
the primary BH. This was due to the accretion disk perturbations of the
secondary BH during its pericenter passages in an orbit which lies in the disk
plane (Sillanpää et al., 1988). Moreover, a notable prediction of the model
was that OJ 287 should show an outburst again in the autumn of 1994. This was
verified by the second campaign called OJ-94 (Sillanpää et al., 1996a).

In 1994, Lehto & Valtonen (1996) recognized that the flares in OJ 287 were not
exactly periodic, and that the systematics of the past flares are better
understood if the flares are double-peaked and the two peaks are separated by
$\sim 1-2$ years. This led to the proposal of a new SMBHB central engine model
(LV96 model hereafter) for OJ 287 where a supermassive secondary BH is
orbiting the more massive primary BH in a relativistic eccentric orbit with a
redshifted orbital period of $\sim 12$ years and where the orbital plane is
inclined with respect to the accretion disk of the primary at a large angle.
In this model, the double-peaked flares occur due to the impacts of the
secondary BH on the accretion disk of the primary BH twice during each orbit.
The third campaign carried out by the OJ-94 group verified the flare on
October 1995. It came within the narrow two-week time window of the prediction
(Valtonen, 1996; Sillanpää et al., 1996b).

Each impact releases hot bubbles of gas from both sides of the disk which then
expand and cool down until they become optically thin (Lehto & Valtonen, 1996;
Ivanov et al., 1998). At this point, bremsstrahlung radiation from the entire
volume becomes observable causing the optical brightness of the blazar to
increase rapidly (Pihajoki, 2016). The resulting flare can last for a few
weeks and should reach its peak within the course of a few days.

The LV96 BBH central engine model of OJ 287 with subsequent fine-tuning of the
orbital and disk parameters has been very successful in explaining the past
observed optical flares, including many which have been discovered in
historical data since 1994. See Dey et al. (2018) for a collection of flare
light curves. Moreover, this model has accurately predicted the time of
arrival of future flares expected during the years 2005, 2007, 2015, and 2019
which triggered new observing campaigns (Sundelius et al., 1996, 1997;
Valtonen et al., 2006; Valtonen, 2007; Dey et al., 2018). The 2005 flare came
a week ahead of the orbital schedule calculated by Sundelius et al. (1996) and
corrected for the astrophysical time shift between the disk impact and the
flare by LV96 (Valtonen et al., 2006a; Valtonen & Sillanpää, 2011). The latest
prediction of the 2019 Eddington flare was successfully observed with the
Spitzer space telescope (Laine et al., 2020) which showed that the orbital
times are now known with the accuracy of $\pm 4$ hours.

The model also predicted that the excess emission during the impact flares is
unpolarised, and it has the bremsstrahlung spectrum at $\sim~{}3\times 10^{5}$
K. At the time (in the spring of 1995 when the model was constructed) it was
only known that during the rapid brightening of the 1983 flare the
polarisation decreased with rising flux: the first point at $\sim~{}16.5$ mJy
in the V-band had polarisation of $8.2~{}\%$ while the peak at $\sim~{}33$ mJy
had polarisation of $0.8~{}\%$ and the four points measured around the peak
had the average polarisation of $4~{}\%$ (Smith et al., 1985). This is exactly
what is expected if the excess emission is unpolarised. The same trend was
seen again in 2007 and 2015 (Valtonen et al., 2008, 2016). In 1984, 1994,
1995, 2005 and 2019 the polarisation measurements were not obtained for the
primary peak, partly due to seasonal and weather restrictions.

The spectral slope over a wide range of frequencies was measured during the
2005 flare and the expected spectral index of $\alpha_{\nu}\sim-0.2$ was
found. It is quite different from the usual spectral index of
$\alpha_{\nu}\sim-1.3$ over the same spectral range (Valtonen et al., 2012;
Kidger et al., 2018). The spectral index $\alpha_{\nu}\sim-0.2$ was confirmed
during the 2015 and 2019 flares. The complete change in the nature of emission
in the flare as compared with the usual out-of-flare emission excludes the
possibility that the flares are somehow related to Doppler boosting variations
in a turning jet (Villata et al., 1998; Rieger, 2004).

As an example, the 2015 flare was observed by the Swift observatory (Ciprini
et al., 2015). The first observation prior to the flare on JD 2457354 found
the flux 1.92 mJy in the UV channel M2. Almost simultaneously it was observed
in several optical bands, including the V-band, where the flux was recorded as
5.67 mJy (MacPherson et al., 2015). The highest M2 flux value in this set of
Swift observations was recorded on JD 2457360.65: it was 5.55 mJy, while the
nearly simultaneous V-flux (about 0.5 d earlier) was 10.0 mJy (Shappee et al.,
2015). Our R-band magnitude 13.36 on JD 2457360.65 is in agreement with this
V-flux. We take the first observations as the base level, and obtain the flare
contributions in M2 and V-bands of 3.63 mJy and 4.33 mJy, respectively. The
spectral index between these two wavebands becomes $\alpha_{\nu}\sim-0.2$, the
bremsstrahlung value for the temperature in question.

The observation of the bremsstrahlung optical flare during early September
2007 confirmed the highly-relativistic nature of the binary BH in OJ 287. This
is because the BBH model would have predicted the arrival of OJ 287’s impact
flare on Earth around early October 2007 if the effects of GW emission had
been ignored while describing the relativistic trajectory of the secondary BH
in the blazar (Valtonen, 2007; Valtonen et al., 2008). Thereafter, the
observation of the GR centenary flare during November-December 2015 was used
to constrain the dimensionless Kerr parameter of OJ 287’s primary BH Valtonen
et al. (2016); Dey et al. (2018). The updated BBH parameters of the OJ 287
system, determined using post-Newtonian (PN) equations of motion, are given
below: primary BH mass $m_{1}=18.35\times 10^{9}M_{\odot}$, secondary BH mass
$m_{2}=150\times 10^{6}M_{\odot}$, primary BH Kerr parameter $\chi_{1}=0.38$,
orbital eccentricity $e=0.657$, and orbital period (redshifted) $P=12.06$
years Dey et al. (2018). We now describe the various EM observational
campaigns that allowed us to extract additional astrophysical details of the
BBH central engine description for the blazar.

###  1.2 Multi-wavelength emission of OJ 287 and previous and ongoing
monitoring campaigns

OJ 287 has been the target of intense multi-wavelength coverage and monitoring
in recent decades. Early optical observations first revealed very bright high-
states reaching up to $\sim$12 mag. Dedicated optical monitoring campaigns
(Pursimo et al., 2000; Villforth et al., 2010; Pursimo et al., 2021) and the
addition of photographic plate data dating back to the late 19th century
(Hudec et al., 2013) then revealed that bright optical maxima recur semi-
periodically. These observations led to first suggestions of the presence of a
binary SMBH in OJ 287 (Sillanpää et al., 1988) and its later confirmation in
subsequent campaigns, as mentioned above (Valtonen & Sillanpää, 2011; Valtonen
et al., 2017). An ongoing optical monitoring program of OJ 287 combines
optical data taken at multiple ground-based observatories and provides
coverage at sub-daily cadence except for epochs when OJ 287 is unobservable
due to its proximity to the Sun (Zola et al., 2016; Dey et al., 2018). Figure
1 shows the optical lightcurve of OJ 287 from late 19th century to the present
epoch. Based on the latest optical high-state discovered in this program in
late 2015 and interpreted as a binary impact flare (Valtonen et al., 2016), a
new broad-band radio – X-ray monitoring program of OJ 287 was initiated and is
described below. Most of the time, the optical emission of OJ 287 is dominated
by synchrotron emission from the jet, as known from SED modeling (Abdo et al.,
2009) and high levels of polarization (D’Arcangelo et al., 2009). However, it
shows a decrease in polarization, expected in the binary SMBH model, at epochs
where an additional thermal bremsstrahlung contribution from the disk impacts
contributes to the emission (e.g. Valtonen et al., 2008). Optical polarimetric
monitoring of OJ 287 was carried out during the epoch 2016-2017 following
changes in polarization and led to the interpretation of the 2016-17 optical
high-state as a binary after-flare (Valtonen et al., 2017).

Figure 1: Optical lightcurve of OJ 287 between 1888 and 2021. The figure
includes a large number of previously unpublished data points by two of the
authors (R.H., S.Z.). All point have been transformed to the V-band, including
data in Laine et al. (2020).

Based on the faintness of its optical Balmer lines from the broad-line region
(BLR; Sitko & Junkkarinen, 1985; Nilsson et al., 2010), OJ 287 is classified
as a BL Lac object. Given its spectral-energy distribution (SED), it is of LBL
type (low-frequency-peaked blazar; Padovani & Giommi, 1995). Optical emission
lines from the narrow-line region (NLR) like [OIII]$\lambda$5007 are only
faintly present (Nilsson et al., 2010).

As a blazar, OJ 287 is a bright radio emitter (Dickel et al., 1967) with a
structured, relativistic jet (Jorstad et al., 2005; Lee et al., 2020; Gomez et
al., 2021) that is highly polarized (Goddi et al., 2021). On average, the jet
is pointing close to our line of sight (Jorstad et al., 2005; Hodgson et al.,
2017), but it shows large swings in its position angle (Agudo et al., 2012;
Valtonen & Pihajoki, 2013) that can be understood in the context of the binary
SMBH model (Dey et al., 2021).

OJ 287 was first detected in the $\gamma$-ray regime with the CGRO/EGRET
satellite observatory (Shrader et al., 1996) and shows repeated flaring in the
Fermi $\gamma$-ray band (e.g., Abdo et al., 2009; Agudo et al., 2011; Hodgson
et al., 2017) along with more quiescent epochs not detected by Fermi. The
first very-high energy (VHE; $E>100$ GeV) detection was obtained in 2017 with
VERITAS (O’Brien, S. et al., 2017); an observation triggered by the
measurement of an exceptionally bright X-ray outburst (Figure 2) discovered
with Swift in the course of the dedicated long-term, multi-frequency
monitoring program MOMO (Multiwavelength Observations and Modelling of OJ 287;
(Komossa et al., 2017, 2020a)).

Figure 2: X-ray (0.3–10 keV) lightcurve of OJ 287 obtained in the course of
the MOMO program with Swift between 2015 December and 2020 March 1. Two bright
outbursts were discovered in 2016-17 and 2020. Adopted from Komossa et al.
(2021d).

OJ 287 is a bright and variable X-ray emitter observed with most major X-ray
observatories and first detected with the Einstein satellite (Madejski &
Schwartz, 1988). The X-ray spectrum of OJ 287, based on XMM-Newton
observations between 2005 and 2020 (energy range: 0.3–10 keV), is
characterized by a super-soft (low-energy) synchrotron component and a hard
(high-energy) inverse Compton (IC) component (Komossa et al., 2021a). The
synchrotron component dominates the X-ray outburst states (Komossa et al.,
2020a). High-energy emission up to 70 keV was detected with NuSTAR with an
unusually soft power-law spectrum of photon index $\Gamma_{\rm X}=2.4$
accompanying the 2020 outburst (Komossa et al., 2020a). OJ 287 shows an X-ray
jet extended by 20" (Marscher & Jorstad, 2011).

Starting in late 2015, OJ 287 has been the target of the dedicated multi-
wavelength monitoring program MOMO (Komossa et al., 2017, 2020a, 2021a, 2021c,
2021d) (see the review by Komossa et al. (2021b)). MOMO is aimed at
characterizing the blazar properties and facets of the binary SMBH in detail.
The program employs multiple frequencies in the radio band between 2-40 GHz
obtained with the Effelsberg observatory, UV–optical Swift data in 6 filters,
and the Swift X-ray band (0.3–10 keV). The monitoring cadence is as short as
1-4 days. Additional deeper observations with multiple ground and space-based
observatories are triggered in the case of exceptional outburst or low-states
discovered with Swift and/or in the radio band. The community is alerted in
form of Astronomer’s Telegrams about particular states of OJ 287 (ATel #8411,
#9629, #10043, #12086, #13658, #13702, #13785, and #14052). MOMO is the
densest multi-wavelength monitoring program of OJ 287 involving X-rays, with
$>$1000 data sets so far taken, providing spectra, broad-band SEDs, and
lightcurves, during all activity states of OJ 287 (Komossa et al., 2021b, d).

Some of the main results of the program MOMO are as follows: Two bright
outbursts were identified in the course of MOMO in 2016-17 and 2020 (Fig. 2);
the brightest so far recorded in X-rays (Komossa et al., 2017, 2020a, 2021d)).
These outbursts were non-thermal in nature, and they are consistent (Komossa
et al., 2020a) with after-flares predicted by the binary SMBH model, when disk
impacts of the secondary lead to the triggering of new jet activity of the
primary SMBH. OJ 287 has turned out to be one of the most spectrally variable
blazars in the X-ray regime (Figure 3) with power-law spectral indices as
steep as $\Gamma_{\rm X}=3$ (Komossa et al., 2021a, d). Supersoft synchrotron
emission dominates the outburst states. Another key features of the MOMO
lightcurve is a pronounced, symmetric UV–optical deep fade in 2017 November –
December (Figure 3) of still unknown nature (Komossa et al., 2020a, 2021d).
Structure function analyses of Swift data between 2015–2020 has revealed
characteristic break times of 4–39 days depending on activity state and
waveband (Komossa et al., 2021d). Discrete correlation analysis between UV,
optical and X-ray bands during the same epoch shows lags and leads in X-rays
w.r.t. the UV, attributed to an IC component, where external Comptonization is
favored over synchrotron-self-Compton emission (Komossa et al., 2021b, d).

Figure 3: Multi-wavelength X-ray–optical lightcurve of OJ 287 taken with Swift
since December 2015. The majority of data was obtained during the MOMO
program. Panels from top to bottom: Observed X-ray flux in the band (0.3–10)
keV in units of 10-11 erg/s/cm2, X-ray powerlaw photon index $\Gamma_{\rm x}$,
Swift UVOT UV and optical fluxes (corrected for Galactic extinction) in units
of 10-11 erg/s/cm2. Every year, observations are interrupted by $\sim$3 months
when OJ 287 is in close proximity to the Sun such that Swift cannot observe.
Selected epochs are marked in color: (1) The epoch near the 2015 binary impact
flare, (2) the bright 2016-17 outburst, (3) the 2017 optical-UV deep fade, and
(4) the 2020 April outburst. Epochs (2) and (4) were interpreted as possible
binary after-flares. Epochs like the deep fade in 2017 (feature 3) and the new
low-state in September 2020 can be used to obtain host galaxy imaging of OJ
287 when the blazar emission itself is least disturbing. Adopted from Komossa
et al. (2021d).

The MOMO program continues as OJ 287 nears its next predicted binary SMBH
impact flare. Results obtained so far including $>$1000 broad-band SEDs will
serve as a useful data base to identify new activity states of OJ 287 driven
by blazar variability in comparison to binary-driven activity, and allow
detailed modelling of the broad-band emission from the radio to the high-
energy regime at each epoch. We now take a closer look at the present
description of the primary BH accretion disk and propose a detailed scenario
to describe the outcome of secondary BH impact with the disk.

##  2 Accretion disk in OJ 287

Our binary model of OJ 287 is based on the interpretation of the huge optical
flares as impacts of a secondary black hole on accretion disk of the primary.
While the emission level of OJ 287 is between 1 and 2 mJy in the optical
V-band most of the time, the emission rises to 13 mJy (in the 2015 flare) or
similar values in the peaks that we interpret as impact flares, i.e they
represent an order of magnitude increase in the optical flux. The increase
occurs on the time scale of a day or so which is about a fraction of 0.0002 of
the 12 yr (9 yr redshift corrected) binary orbital period. This implies that
the emission region in the source responsible for these huge flux variations,
cannot be greater than one light day in size. Most of the work on the binary
black hole model has concentrated on the orbit of the secondary. However, it
is equally important that we also find out the basic properties of the
accretion disk from the way the disk reacts to the impacts, and how these one-
light-day emission regions are activated.

In the original binary black hole model of Lehto & Valtonen (1996) the disk
was taken as a standard disk with two parameters, viscosity $\alpha$ and
accretion rate $\dot{m}$. Fixed values were given for these parameters. When
more impacts were observed later, both in the future and in the past
historical light curve, it became possible to solve the values of these
parameters together with the orbit solution. Altogether it is a 9 parameter
problem which has a unique solution as soon as the exact times of 10 flares
are known. Mathematically it is equivalent to a ninth order equation which has
nine unique roots. The method of solution is iterative, similar to the Newton-
Raphson algorithm (Valtonen, 2007). From the detailed study of the 2015 flare,
the values of the two standard parameters were determined to be $\alpha=0.26$
and $\dot{m}=0.08$ in Eddington units.

In the model of Lehto & Valtonen (1996), a spherical bubble of hot plasma is
released from the disk by the impact. Ivanov et al. (1998) later demonstrated
that actually two bubbles are released, one on each side of the disk.
Therefore each impact is detectable via a radiating bubble on the front side
of the disk. The model further specified that the bubble expands isotropically
and therefore retains its spherical shape. The simulation of Ivanov et al.
(1998) casts doubt on whether this is true at an early stage of the bubble,
but this was not considered important, as the bubble becomes optically thin
only at a later time which results in the observed great increase in
brightness.

Valtonen et al. (2019) noted that the maximum brightness of the 2015 flare
should have been 26 mJy in the standard model. The fact that the peak value
was only half of this value is most easily explained if the bubble is
compressed along the line of sight. Thus it appears necessary that the
expansion of the bubble is not completely unconstrained, but that it is
influenced possibly by the general magnetic fields in the disk. In that case,
it would be interesting to know what happens to the flux of the bubble from
the very beginning if its expansion is constrained all the way. From Ivanov et
al. (1998) one could conclude that initially the expansion is more jet-like, a
flow perpendicular to the disk rather than a gentle release of bubbles. Ivanov
et al. (1998) called them fountains. In the following we consider what
difference this makes to the impact flare light curve.

##  3 Anisotropically expanding BH impact generated bubbles

According to Lehto and Valtonen (LV in the following, (Lehto & Valtonen,
1996)), the compressed plasma inside the accretion disk, following the impact,
bursts out of the disk with the speed $v_{perp}=(6/7)v_{sec}$ (LV Eq. 7) and
the time to come out of the disk is $\tau_{dyn}=h/v_{perp}$ (LV Eq. 8), where
$v_{sec}$ is the speed of the secondary black hole at the impact and $h$ is
the semi-height of the disk. The flux from this bubble follows
$S\sim\epsilon^{ff}V_{bubble}$ (LV Eq. 14), if the bubble is fully
transparent. The quantities $\epsilon^{ff}$ and $V_{bubble}$ are the free-free
emissivity and the volume of the bubble, respectively. However, we see only
the fraction $1/\tau$ of the bubble, where $\tau$ is its optical depth. Here
$\tau=d\epsilon^{ff}T^{-1}$, if $d$ is the geometrical depth and $T$ the
temperature of the bubble. Therefore the flux goes as

$S_{seen}\sim\epsilon^{ff}V_{bubble}/(d\epsilon^{ff}T^{-1})\,.$ (1)

We can substitute $T\sim V_{bubble}^{-1/3}$ in adiabatic expansion which means
that

$S_{seen}\sim V_{bubble}^{2/3}/d\,.$ (2)

We can put $V_{bubble}\sim r^{2}d$, where $r$ is the bubble radius ($d\sim r$
in spherical geometry). The adiabatic expansion speed
$v\sim~{}T^{1/2}\sim~{}V_{bubble}^{-1/6}$ and

$d\sim vt\sim V_{bubble}^{-1/6}\,t\,.$ (3)

Now we have to specify the geometry of the bubble: For a tube of fixed radius
seen end-on, $V_{bubble}\sim dr^{2}$, $r=constant$, or $S_{seen}\sim
d^{-1/3}$, and since $d\sim d^{-1/6}t$, it follows that $d\sim~{}t^{6/7}$.
Thus

$S_{seen}\sim t^{-2/7}\,.$ (4)

Flux drops with time roughly as $S_{seen}\sim t^{-1/3}$.

For a tube of fixed depth seen end-on, $V_{bubble}\sim dr^{2}$, $d=constant$
or $S_{seen}\sim r^{4/3}$, and since $r\sim r^{-1/3}t$, it follows that $r\sim
t^{4/3}$ . Thus

$S_{seen}\sim t^{16/9}.$ (5)

Flux increases with time nearly as $S_{seen}\sim t^{2}$ (c.f. Figure 2 in
(Valtonen et al., 2019))

Figure 4: The optical light curve of OJ 287 during 2005 in V-band mJy units.
Plot shows the monthly averaged data points for the early part of 2005 while
the latter part, with faster variability, has 3 day averaged data points. The
dashed curve is the theoretical line based on present work and Valtonen et al.
(2019). Additionally, the theoretical curve provides the expected light-curve
of OJ 287 during 2021-2022 and we mark its epochs on the upper x-axis while
employing the time shift of the SMBH binary central engine description Dey et
al. (2018). This is motivated by the fact that the secondary BH impact
distances from the center should be roughly the same during 2005 and 2021. The
first big flare during 2022 is not expected to be visible from any Earth-based
facility while the secondary peak may become barely visible by early September
2022.

For standard spherical bubbles of radius $r$ the corresponding relation is
(Pihajoki, 2016)

$S_{seen}\sim V_{bubble}^{2/3}/r\sim r\sim t^{2/3}.$ (6)

Table 3 lists the disk impacts between 1971 and 2031. The first two columns
come from (Dey et al., 2018), the numbers from the third column onwards are
derived in the next section.

Times of disk impacts and bremsstrahlung flares as well as the radii of the
$H_{\alpha}$ bubbles, impact speed on the disk, the delay of the $H_{\alpha}$
flare and the time of the $H_{\alpha}$ flare Impact time time of flare
${R_{H{\alpha}}}$ (AU) ${v_{rel}/c}$ delay (yr) ${H{\alpha}}$ flare 1971.11
1971.13 2860 0.34 0.72 1971.83 1972.71 1972.93 3930 0.20 1.66 1974.37 1982.83
1982.96 3160 0.30 0.81 1983.74 1984.07 1984.12 3250 0.28 0.91 1984.98 1994.48
1994.59 3610 0.23 1.31 1995.79 1995.82 1995.84 2940 0.33 0.76 1996.58 2005.11
2005.74 4360 0.16 2.31 2007.42 2007.68 2007.69 2740 0.36 0.65 2008.33 2013.48
2015.87 5000 0.12 3.52 2016.80 2019.56 2019.57 2740 0.36 0.65 2021.21 2021.92
2022.55 4360 0.16 2.31 2024.23 2031.40 2031.41 2900 0.33 0.74 2032.14

In the spring of 2005 there was a major flare, 6.6 mJy at 2005.12, one of the
biggest deviations from the standard binary black hole model during years
1996-2010 (Sundelius et al., 1996, 1997; Valtonen et al., 2011). Figure 4
shows the observations (monthly averages) from 2005, prior to the main flare
in October of that year. The disk impact occurred at 2005.11 (February 9,
2005). With a time translation such that December 3, 2021 corresponds to
February 9, 2005, we generate a predicted light curve for the 2021/22 season.
For the values appropriate for 2005 and 2021 impacts, which take place at the
same radial distance from the center, $\tau_{dyn}\sim 20$ days. The flaring
related directly to the disk impact should start around December 23, 2021, and
the flares should be recognizable as bremsstrahlung flares in the optical - UV
region, with no counterpart in radio nor X-rays. The on-going and planned
observational campaigns should be able to provide excellent constraints on our
detailed modeling. So far, the period of early December 2021 was covered with
Swift in the ongoing MOMO project and with ground-based optical observations.
The optical–UV low-state detected in December 2–14 implies a small shift in
the onset of the flare (work in preparation).

We now move on to explore other observational implications of the expanding
impact bubble. This requires an extension of LV96 (Lehto & Valtonen, 1996)
where the brightness evolution of the bubbles was terminated at the epoch when
the bremsstrahlung flares emerge. However, at a later stage when the bubbles
have cooled enough, they may become visible again as hydrogen line flares. In
what follows, we study the timing of these proposed $H\alpha$ flares, and
connect them with one epoch when the hydrogen lines were seen to be
exceptionally bright, by a factor of ten above the normal.

##  4 Possible spectroscopic implications of the evolving BH impact generated
bubbles

After the bremsstrahlung flare the radiating bubble cools until it falls below
the jet continuum level and is not observed any more (Valtonen et al., 2019).
However, at a later time the bubble has cooled enough to become an $HII$
region. This may appear as a sudden increase of the strength of the Balmer
lines for a limited period of time. At the end of 1984, Sitko & Junkkarinen
(1985) observed such an increased line flux, approximately a factor of ten
higher than normal (Nilsson et al., 2010). To our knowledge, this was a unique
occasion when such an $H\alpha$ flare has been seen. Therefore it is tempting
to associate this epoch with the dying moments of one the bubbles that had
generated a bremsstrahlung flare at an earlier time.

In order to estimate when the bubbles left over from disk impact may reach the
$H\alpha$ stage, we calculate the required increase in the bubble radius until
its temperature has dropped to $10^{4}$ K. Using Table 3 of Dey et al. (2018)
and Table 1 of Valtonen et al. (2019), and assuming an adiabatically expanding
bubble, we obtain the radii of the bubbles listed in column 3 of Table 3. In
Valtonen et al. (2019) the expansion speed of the bubble is given as $\sim
v_{rel}$/4 where $v_{rel}$ is the speed at which the secondary impacts the
accretion disk. In fact, the expansion speed value $v_{rel}$/4.6 is used in
the following since it gives the correct timing for the 1984 $H\alpha$ flare,
1984.98 (Sitko & Junkkarinen, 1985). This value is well within the uncertainty
of the expansion speed, extracted from the numerical simulations by Ivanov et
al. (1998). These values are given in column 4 of Table 3, while column 5
gives the time of the expansion up to this stage. In column 6 we give the
estimated time of the beginning of the $H\alpha$ flare.

The epochs of observation in the line monitoring program of Nilsson et al.
(2010) and in the program by some of the present authors at the Nordic Optical
Telescope (T.P.) as well at the Rozhen Telescope of Bulgarian Academy of
Sciences (S.Z.) are given in Table 4.

The epochs of the line monitoring program. time telescope Strong lines?
secondary line 2005.26 VLT no (8059) 2005.89 VLT no (8092) 2006.25 VLT no
(8059) 2006.93 VLT no (8637) 2007.26 VLT no (9405) 2008.01 VLT no [7418]
2008.27 VLT no [7787] 2010.12 NOT no 6438 2010.81 NOT no 6494 2010.88 NOT no
6499 2011.01 NOT no 6506 2011.05 NOT no 6508 2011.11 NOT no 6512 2011.31 NOT
no 6521 2011.35 NOT no 6522 2011.75 NOT no 6538 2011.88 NOT no 6542 2012.14
NOT no 6549 2012.99 NOT no 6564 2013.18 NOT no 6566 2015.94 Rozhen no (8273)
2016.26 Rozhen no (8300) 2021.82 NOT no 6708

The expected flares in the time span and close to these observations were at
2007.42, 2008.33, 2016.80 and 2021.21. In the first two cases the
observational epochs were too early by 0.16 and 0.05 yrs, respectively, in
order to catch the beginning of the flares. The last two were too early by
0.54 yr and too late by 0.61 yr. Therefore, if the $H\alpha$ flares started as
expected, we would not have seen them, even if the flares lasted several
months.

In order to calculate the line flux, we proceed as follows. The emissivity of
the Balmer $H{\beta}$ line may be taken as

$\epsilon_{H{\beta}}=8.3\times 10^{-26}~{}n^{2}~{}erg~{}cm^{-3}~{}s^{-1}\,,$
(7)

where $n$ is the number density of Hydrogen atoms. We calculate the line flux
using the density for the 1984 line flare, $n=3.2\times 10^{8}$ and the bubble
radius 3250 AU from Table 3. The result is

$S_{H\beta}=2.2\times 10^{-15}~{}erg~{}s^{-1}.$ (8)

Sitko & Junkkarinen (1985) find $S_{H{\beta}}\sim 1.2\times
10^{-15}~{}erg~{}s^{-1}$, but state the value as uncertain. The more certain
$S_{H{\alpha}}=24\times 10^{-15}~{}erg~{}s^{-1}$, divided by some number of
order of 10 (Sitko & Junkkarinen, 1985) takes us close to the above estimate.
The radius and the number density are different in 2024 than in 1984, but the
different factors practically cancel each other. Therefore we expect the 1984
$H\alpha$ line flare to repeat itself in 2024.

##  5 Explaining observed gamma ray flares

Here we will distinguish the classical gamma-ray regime, as observed by the
Fermi satellite, with many gamma flares detected (typically in the $E\sim 100$
MeV regime); and the VHE emission ($E>100$ GeV) that this Section mostly
focuses on. At VHE, OJ 287 was only detected once, by VERITAS. The VERITAS
observation was triggered by the Swift detection of an exceptional X-ray high-
state in the course of the MOMO program (Komossa et al., 2017, 2020a).
Otherwise, OJ 287 is monitored in gamma-rays with Fermi, with many large and
small flares.

Classical gamma ray flares have been reported many times in OJ 287 (Agudo et
al., 2011; O’Brien, S. et al., 2017). For example, the 2009 classical flare
coincides with the last of the predicted tidal flares following the 2007
bremsstrahlung flare (Sundelius et al., 1997) while the 2017 VHE flare came
soon after the last of the predicted tidal flares following the 2015
bremsstrahlung flare (Pihajoki et al., 2013). We can predict the next VHE
gamma ray flare following the same principles.

In order to study the orbital phase of the 2017 VHE flare in our model, we
have calculated the response of the OJ 287 accretion disk to the 2013
secondary crossing, using a disk of 146880 particles. The particles leaving
the disk are counted, and the time ($time1$) and distance ($r$) from the
center are noted down. Assuming a constant transit speed to the jet (taken to
be at $r=3000$ AU), the jet reaction time ($time2$) is calculated. Figure 5
corresponds to the transit speed of $c/6$ in the observers frame which equals
to $0.22c$ in the OJ 287 frame. This transit speed value fits best with the
two observed light curve peaks at 2016.2 and 2016.8. The present calculation
differs from Pihajoki et al. (2013) in that the surface density of the disk
follows Lehto & Valtonen (1996) rather than a uniform density which was used
in Pihajoki et al. (2013). Also Pihajoki et al. (2013) did not take the
transit time into account, so that it would have to be done collectively for
all particles for a comparison with observations. Note that the second peak
becomes stronger when a realistic surface density model is used, in agreement
with observations.

Figure 5: Simulation of the response of the accretion disk to the crossing of
the secondary black hole through the disk before the 2015 bremsstrahlung
flare. Particles which escape the disk are counted, and they are placed on the
time axis assuming that the particle contributes to the light curve after
traveling the distance of $\sim 15,000$ AU with a constant speed of $0.22c$
where $c$ is the speed of light. The exact value of travel distance depends on
the radial distance from the center where the particle is released from the
disk. The particle count in each time box ($time2$) is turned into a
magnitude, assuming that the relation between the particle number and the
total brightness in each box is linear. The first flare in observations at the
end of 2015 is different (impact flare) and is not modeled by this process.

The tidal flares of the light curve of OJ 287 following the 2005 and 2022
bremsstrahlung flares were calculated by Pihajoki et al. (2013). No time
transformation was applied so that the plots in Figure 4 of Pihajoki et al.
(2013) are for $time1$, as defined above. Using the same transit speed of
particles to the jet as above we should add 1.0 yr to the time axis. For the
2005 impact the two main flares in the predicted light curve are at 2007.4 and
2008.6. Even though these are placed at the summer period of few or no
observations (OJ 287 is close the sun at the beginning of August), a strong
brightening of the source was observed in the following autumn (both 2007 and
2008) which agrees with the tidal flare explanation. The interpretation of the
light curve is made complicated by the fact that the bremsstrahlung flare of
2007.7 falls in the same time interval as the tidal flare (Valtonen &
Sillanpää, 2011).

As to the 2022 case, the predicted light curve has maxima at 2023.5 and 2024.2
(Pihajoki et al., 2013). Note that the nature of the simulation by Pihajoki et
al. (2013) is such that the tidal flares appear to get weaker with time. This
is an artifact of the simulation technique which does not replace the escaped
particles by new ones in the disk. Because of practically the same impact
distance in 2005 and 2022, the tidal responses should be rather similar. We
therefore expect that VHE gamma rays would again be detectable in the spring
of 2024. The prediction cannot be more precise until we see the corresponding
flare starting in the optical region. Further, we expect to see such flare
again in the X-ray regime, as in 2017 when the very bright X-ray flare
detected with Swift triggered the VHE observation. Given the ongoing
monitoring program with Swift with very dense coverage and $>$1000 SEDs
already obtained (e.g., Komossa et al., 2020a, 2021d), very detailed timing
and comparison with the 2017 flare will then be possible. This ongoing
project, or any X-ray monitoring with future X-ray missions, can then also be
used to trigger new VHE observations.

In addition, OJ 287 observations can reveal interesting astrophysical
information when it enters low luminosity state and this is what we focus on
in the next section.

##  6 Optical multicolor photometry at minimum light

From time to time, the brightness of OJ 287 falls to a very low light level.
It is hard to predict when these times will happen again, even though it is
thought that the low light level may be associated with the secondary moving
in front of the primary (Takalo et al., 1990). The absolute minimum light
observed so far took place in the spring of 1989. In Figure 6 we re-plot the
observations (Takalo et al., 1990), and fit smooth curves through them. When
the brightness of the AGN falls towards the brightness of the host galaxy, the
color of OJ 287 is expected to change gradually from the typical jet color to
the color of the host galaxy as the fade proceeds. Takalo et al. (1990)
conclude that at the brightness level V = 17.4, the host galaxy contributes at
least half of the brightness of the source. Using the usual conversion V $-$ K
= 3.2, K $-$ correction and evolution correction of 0.6 and the distance
modulus m - M = 41, this implies the absolute K-band magnitude of
$M_{K}=-26.8\pm 0.3$ (Valtonen et al., 2021).

Figure 6: The observations of the 1989 fade of OJ 287 with the Nordic Optical
Telescope and the Jacobus Kapteyn Telescope at La Palma, with $2\sigma$
errorbars. Here we plot the magnitudes in three filters B, V and R: the I band
data were too few to show a significant variation. The lines are fifth order
polynomial fits through the data points. We note that the dip on May 18 was
deepest in the B band and less deep going to longer wavelengths. There is a
possible time delay in the sense that the dip starts from longer wavelengths
and happens later at shorter wavelengths. At 37 GHz the minimum took place on
May 3, 1989

Using surface photometry in the K-band with the Nordic Optical Telescope
Nilsson et al. (2020) find $K=-14.65\pm 0.6$, leading to $M_{K}=-25.85\pm 0.6$
with the K-correction -0.5, which is somewhat, but not significantly, fainter
than the value measured by the color method. However, this may not be the
total brightness of the host. It is well known that galaxies of this
brightness category have extensive haloes which contribute as much as one
magnitude to the host. For example, the absolute K-magnitude of NGC 4889 is
$M_{K}=-25.4$ (Sandage, 1972) (using m $-$ M = 34.7 and V $-$ K = 3.2)
measured out to roughly the same 30 kpc linear distance as the surface
photometry of OJ 287 in Nilsson et al. (2020) while its total magnitude out to
100 kpc distance is about $M_{K}=-26.8$ (Graham & Scott, 2013). Typically the
outer envelope contributes about 0.8 magnitudes to the total brightness (De
Vaucouleurs & De Vaucouleurs, 1970; Läsker et al., 2014; Huang et al., 2013).
Actually, if taken to the redshift of $z=0.3$, NGC 4889 would probably appear
very much like the OJ 287 host. Even the estimated mass of its central black
hole is practically the same as the primary BH mass in OJ 287 (Graham & Scott,
2013).

We also note that the color photometry method is more likely to capture the
faint halo light than the surface photometry method since the latter is
seriously influenced by the high background light level generated by the
fringes of the strong point source in the center (Yanny et al., 1997). Also
the signal-to-noise ratio is better at longer wavelengths than in the
traditional optical range which may be the reason why it is difficult to get a
host magnitude measurement in the shorter wavelength bands below the infrared
K-band (Nilsson et al., 2020). For example, the Hubble Space Telescope image
gives the range from $M_{K}=-24.7$ to $M_{K}=-26.3$, converted from the I-band
to the K-band in the usual way (Yanny et al., 1997). The magnitude of the
i-band image taken at Gran Telescopio Canarias corresponds to $M_{K}=-25.2\pm
0.3$ (Nilsson et al., 2020).

It is interesting to see where this magnitude by multicolor method places OJ
287 in the black hole mass - host galaxy mass diagram. Considering that OJ 287
is at a much higher redshift than galaxies normally displayed in this diagram
(Saglia et al., 2016) and that the correlation is redshift dependent
(Portinari et al., 2012), we find that OJ 287 follows almost exactly the well
established correlation (Figure 7).

Figure 7: The galaxy sample of Saglia et al. (2016) with three additional
galaxies in the black hole mass vs. galaxy stellar mass diagram (stars), while
the position of OJ 287 is given as a solid square. The dashed line represents
the correlation found by Portinari et al. (2012) for galaxies in the upper -
right hand part of the graph around redshift 0.5.

When do we expect to be able to repeat the host galaxy brightness measurement
by the photometric method? It depends on the reason behind the fade. At the
time in 1989 the secondary was near apocentre of its orbit, and the major axis
of the orbit was close to perpendicular to the disk. This configuration
minimizes the tidal perturbation of the disk by secondary. The next time a
similar configuration arises is around 2032 (see Figure 8).

Takalo et al. (1990) suggests that during the fade the secondary moves between
us and the primary, and this causes a temporary misalignment of the jet,
leading to decreased Doppler boosting that shows up as a fade. This is a
sensible proposition as we may find as follows: Suppose the secondary passes
within 100 AU of the jet. Then the jet is gravitationally deflected by $\sim
4^{\circ}$. Using the parameters from Hovatta et al. (2009), and adding
$4^{\circ}$ to the viewing angle, we find that the jet becomes fainter by
about 3 magnitudes. That is, if we have the normal V-band magnitude $\sim 15$,
in the fade it may drop to $\sim 18$ which is actually what happened during
the 1989 fade. We cannot perform a more detailed calculation since we do not
know the inclination of the orbit of the secondary. This proposition could be
studied further, as the jet swinging back and forth should appear also in
radio maps, with the proper delay for each frequency, as well as in optical
polarisation with a small delay. A comparison of optical and radio light
curves at the fades should also establish the relative positions of the radio
and optical cores in the jet. The 1989 fade suggests that they lie close to
each other.

We may ask when this configuration repeats itself. Figure 8 and its caption
lists such times. Since the optical emission region is in the jet a little way
from the center, the times given by the model are about half-a-year earlier
than the response in the optical emission. The 1998 and 2008 fades were
discussed in Pietilä et al. (1999) and Valtonen et al. (2011) while the
2020.05 alignment, with the same response time as in 1989, would place the
fade in the summer period of 2020 when we had no observations (Figure 8).

Figure 8: Temporal evolution of the angular separation between the primary jet
in OJ 287 and the line joining the primary and secondary BHs. The separation
is calculated assuming that the jet precession follows the disk model of Dey
et al. (2021). Minimum separation angle occurs at 1980.73, 1988.95, 1997.61,
2008.45, 2020.05, and 2031.85.

Obviously a fade would be a good time to repeat the surface photometry study
of the host, too. In addition, epochs of very faint intrinsic emission from
this blazar (jet) can be used to trigger host galaxy imaging because the
central point source and its point spread function (PSF) is least affecting
the host measurements at such epochs. For instance, in November 2017 and
September 2020 such optical–UV low-states were identified in Swift monitoring
(Fig. 3) and new ones will be detected in the course of the ongoing MOMO
program as well as in ground-based optical monitoring campaigns.

##  7 The radio jet

OJ 287 has a prominent radio jet which is known for its large swings in the
position angle in the sky (Agudo et al., 2012). These swings are well
understood if they result from the variations in the orientation of the inner
accretion disk, and if these variations are translated into the direction of
the jet (Valtonen & Pihajoki, 2013; Dey et al., 2021). Actually, these large
swings were predicted by Valtonen et al. (2011) before they were observed. The
high resolution data at 86 GHz is compared with a theory line in Figure 9 (see
Dey et al. (2021) for details). A 3.4 yr delay is found between the change of
the disk axis and the direction of the jet. This is because the point at which
the jet position angle is measured is typically 200 microarcsec from the radio
core, and it takes some time to transit the information of the disk
orientation change to this distance in the jet. The actual jet speed in OJ 287
is discussed e.g. in Valtonen et al. (2012a). In this case we find the
relativistic $\Gamma~{}\sim~{}7$. This is somewhat lower than the values
measured by other methods. For example, Savolainen et al. (2010) give
$\Gamma~{}\sim~{}9.3$, Jorstad et al. (2005) find $\Gamma~{}\sim~{}16.5$, and
Hovatta et al. (2009) find $\Gamma~{}\sim~{}15.4$. Anyway, our
$\Gamma~{}\sim~{}7$ is close enough, and gives further support for the model.
If the disk model was not accurate, or if the principles of jet-disk
connection were not valid, there is no reason why the value coming out of this
model would be anywhere near the values obtained via quite different methods.

Figure 9: The position angle of the radio jet in 86 GHz observations (points)
compared with the model line. The best-fit parameters are taken from Table 1
of Dey et al. (2021) and we have used the disk model to calculate the
theoretical model line. The viewing angle varies with time but is typically
similar to the value found in Hovatta et al. (2009). The delay time between
the change of the jet direction at its origin and the appearance of the change
in 86 GHz observations is 3.4 yr.

The region of the jet which we observe depends on the frequency because the
resolution of the global VLBI maps is proportional to the observing frequency.
Unfortunately this technique cannot be extended as far as the optical
frequencies. However, some information about the jet may be obtained by
monitoring of optical polarization. Valtonen et al. (2012) demonstrate that
the optical polarization angle is a good proxy for the position angle of the
radio jet. A big position angle jump occurred in 1995 in the optical
polarization while a similar jump took place in 2009 in 15 GHz radio maps.
Figure 10 shows how the optical polarization is modeled by our swinging jet
model. The time delay is one quarter of the 86 GHz delay, so the distance of
the region emitting in optical should be closer to the center by the same
amount.

Figure 10: The annually averaged optical position angles of OJ 287 (points)
together with the same jet model as in the radio jet observations. Note that
the position angle has large fluctuations after the time of large optical
flares, 1983/84 and 1994/95. The model parameters are the same as for the
radio flares, except that time delay is 0.92 years, i.e. the optical emission
region is closer to the center than the observed radio jet region at 86 GHz.
This region cannot be resolved in optical observations.

In the high resolution map from April/May 2014 (Gomez et al., 2021), the 86
GHz component C2 at 90 microarcsec from the center has the position angle
$-8\pm 45^{\circ}$. The latter component more or less coincides with the point
in the jet where our model places the optical emission region (Figure 11). The
innermost component C1b at 40 microarcsec lies at the position angle $13\pm
10^{\circ}$, interestingly close to the theoretically calculated position
angle the secondary jet at the position angle $14.6\pm 4^{\circ}$ of Dey et
al. (2021) (see Figure 11). Below we explore other possibilities of searching
for observational evidence of the secondary black hole.

Figure 11: The distance and position angle PA of radio components in the
Radioastron VLBI map by Gomez et al. (2021) with error bars. The 22 GHz data
are used when available. The component J1 comes from 15 GHz observations.
Superposed are the theoretical positions of the corresponding components at
two epochs (stars and triangles) in the spin model (Dey et al., 2021; Valtonen
& Pihajoki, 2013). The position of the expected secondary jet is also marked
(cross). Its position angle is independent of epoch.

##  8 On the presence of secondary BH accretion disk

The mass ratio of the two black holes in the OJ 287 binary is 122. Thus for
the same $\dot{m}$ and $\alpha$, the primary disk should be 122 times
brighter. However, if the mass accretion rate is 122 times greater in the
secondary than in the primary, both disks could be equally bright. This would
put the brightness of the secondary over the Eddington limit, which may be
possible at least in some accretion disk models (Heinzeller & Duschl, 2007).
It is well known that tidal perturbations increase the accretion flow
considerably, see e.g. Sillanpää et al. (1988).

We may use Sundelius et al. (1997) and Pihajoki et al. (2013) to estimate the
accretion rate to the secondary. In these simulations one particle represents
$\sim 100$ solar mass of gas and the gas is accreted at the average rate of
$\sim$ 5 particles per year, i.e. the mass accretion rate is about $5\times
10^{2}$ solar mass per year. This is approximately the mass contained within
the Bondi-Hoyle accretion radius for the secondary, and much of it is likely
to remain bound to the secondary when it leaves the disk. At what rate this
matter is accreted to the black hole itself is another problem: it may well be
that only some fraction of the gaseous envelope around the secondary black
hole finds its way to the center before it suffers another collision with the
disk. On the other hand, the steady state accretion rate of the primary is
much lower, $\sim 0.6$ solar mass per year (Valtonen et al., 2019). Thus the
higher accretion rate to the secondary may compensate for its lower mass, in
which case the secondary may not be all that much fainter than the primary.

However, the secondary is truncated by tides which may reduce its optical
emission below the brightness of the primary, since the optical emission is
associated primarily with larger radii in standard disk models. Also the
accretion by the secondary is not steady, and its brightness may vary
considerably depending on the epoch of observation. Not seeing it at one epoch
does not mean that it could not be seen at all.

At present there are several ways by which the secondary could show up. In
spectral observations there could be extra lines that match the redshift of
the secondary at that particular instant in time. Because of the high orbital
speed, these lines would be found quite far apart from the primary redshift in
the spectrum. For example in November 2021 the secondary should have had
redshift 0.38, in contrast to 0.307 in the primary. However, the relative
redshift changes constantly, and moreover, the secondary spends a great deal
of time behind the accretion disk and cannot be seen. After the December 3
crossing to the far side of the disk we should not see it for another ten
years, and even after that only briefly. This is because precession now lines
up the secondary orbit such that the pericentre is towards us, and the
secondary spends only about $1/10$ of its orbital time in the pericenter part
of the orbit.

Table 4, column 4, gives the expected wavelengths (in angstroms) of the
secondary $H\beta$ or $H\alpha$ lines (the latter in square brackets). During
the periods when the secondary should be behind the disk, we give the line for
the opposite viewing direction, since we cannot be absolutely sure that we
have identified the viewing direction correctly. Our interpretation in this
respects relies heavily on Takalo et al. (1990). These values are in curved
brackets, and refer to $H\alpha$. In Figure 12, we show the theoretically
calculated redshift of the secondary BH.

Figure 12: The redshift of the secondary black hole as a function of time. Two
inclination values are shown: $90^{\circ}$ (edge-on) and $50^{\circ}$. The
model does not specify the inclination. The shaded regions are the time
intervals when the secondary is in front of the accretion disk and thus
visible to us. The identification of the "front side" is based on taking the
secondary to be in front of us during the 1989 deep fade (Takalo et al., 1990)

Among the VLT data, there is only one epoch where the expected secondary line
would be within the observed wavelength range, and no line is seen. However,
in this case we would need to have had the viewing direction wrong. Thus there
is nothing in the VLT spectra that excludes the possibility of secondary lines
in our current model. In the NOT and Rozhen spectra we are in the right
wavelength region for $H\beta$, and the NOT spectra are taken at the time when
we should see the secondary lines (see Figure 13). We do not see any of them.
In fact, only a single primary line is detected weakly. Thus this non-
detection does not give strong observational limits.

Figure 13: Optical spectra of OJ 287 taken at the Nordic Optical Telescope.
The continuum has been fitted with a low-order polynomial and subtracted and
the spectra have been shifted in vertical direction for clarity. The positions
of $\lambda$4861 Hβ and $\lambda$50006 [OIII] at the redshift of OJ 287 are
indicated by vertical lines. The red curve below shows the positions of
prominent atmospheric emission lines. Apparent deep absorption structures in
the spectrum are telluric lines. The epochs are (from top to bottom) :
12-02-2010, 23-10-2010, 18-11-2010, 04-01-2011, 17-01-2011, 01-10-2011,
11-02-2011, 23-04-2011, 07-05-2011, 16-11-2011, 19-02-2012, 27-12-2012, and
06-03-2013.

The second way to find signs of the secondary is to look for the jet of the
secondary. It may arise as soon as a new accretion disk has formed after a
disk impact. Since the new accretion disk is formed from the material pulled
out of the primary accretion disk, its axis should share the orientation of
the axis of the mean disk of the primary. This is in contrast to the primary
jet, which wobbles considerably around the axis of the mean disk under the
influence of the companion (Dey et al., 2021). A search for radio emission in
a constant direction at the position angle $\sim 15^{\circ}$ may reveal an
indication of such a jet. The length of the jet should scale with the mass of
the black hole, so that the secondary jet may be a factor $\sim 122$ shorter
than the primary jet (Ghosh & Abramowicz, 1997; Chen et al., 2015). In the
primary jet the outermost component J1 at 15 GHz is found at 3400 microarcsec
from the center; thus the secondary jet may be expected in the scale of $\sim
30$ microarcsec. This is the scale which has just been reached at VLBI (Gomez
et al., 2021) and will be studied in coming years.

Also, some of the flares may arise in the secondary jet rather than in the
primary jet. They could be recognized for their short variability time scale.
However, since the secondary jet is expected to represent only a fraction of
the flux of the primary jet, the fast variations would appear at a low
amplitude, consistent with Valtaoja et al. (1985) who found $1-10\%$ peak-to-
peak amplitude variations during the 1983/84 flaring season, with a 15.7 min
period. Pihajoki et al. (2013) find that variations down to 3.5 d may be
associated with the primary while the variability of the secondary could be
faster by at least a factor of 122, depending on the secondary spin. For a
maximally spinning secondary the 15.7 min period may arise (Pihajoki et al.,
2013). It is the period of the innermost stable orbit about the secondary
black hole in the jet reference frame.

Since the secondary disk is likely to overflow its Roche lobe, we may also
find interesting X-ray emission from the ring in the primary disk where the
overflow stream terminates, akin to X-ray binaries. It is also worth searching
for the direct accretion disk signature of the primary in the SED including a
strong disk contribution in X-rays persistently detected in AGN with SMBHs of
that mass.

##  9 Are there more OJ 287-like systems?

The remaining lifetime of the OJ 287 black hole binary before a merger is
$10^{4}$ yr. Comparing this with the merger time scale of $10^{7}$ yr, there
is only a chance of $\sim 10^{-3}$ to catch the system at this late stage of
evolution. Is this a problem? It depends on the parent population of OJ
287-like systems. As an order of magnitude estimate, it has been deduced that
there are something like $10^{4}$ binary black hole systems in the sky in the
OJ 287 brightness category (Valtonen et al., 2010). If this is the parent
population, some ten of them would be found in the final stages of merger
(Krolik et al., 2019).

How could we find other OJ 287-like systems? We define an OJ 287-like binary
black hole system such that it shows repeated bremsstrahlung flares due to the
impacts on the accretion disk. The main requirement is that the mass ratio of
the two black holes is rather large, so that the smaller black hole does not
disrupt the accretion disk entirely in the process (Valtonen et al., 2012).
This is not very restrictive as large galaxies act as attractors of mergers
with smaller galaxies, and together with the black hole mass - galaxy mass
correlation this produces typically unequal pairs of black holes.

OJ 287 is a blazar with a jet pointing nearly toward us which increases its
brightness many-fold, by about 3 magnitudes by our estimate above. The jets of
other similar systems would likely point in some other direction. Thus it is
unlikely that other OJ 287-like systems are found in lists of brightest radio
sources, such as the Ohio Survey which gives the initial to the name of OJ
287.

However, if such a source is discovered through the repeated bremsstrahlung
flares, it will take time before enough flares are seen to solve the orbit. In
the case of OJ 287 we are lucky that the source lies in the ecliptic zone
which has been frequently photographed over 130 yrs. No such luck can be
expected for the other OJ 287-like systems. It is also rather likely that
these other systems contain smaller binary components by mass than OJ 287,
with weaker bremsstrahlung flares. While other candidate binary black hole
systems have been identified in recent years (review by Komossa & Zensus,
2016), none of them was explained by involving optical bremsstrahlung flares,
and none of them so far has the dense polarimetry monitoring that OJ 287 has
seen. Thus we may find more OJ 287-like systems, but not before full-sky
monitoring programs have been operational for some time. We now briefly list
GW aspects of our description for OJ 287.

##  10 Pulsar Timing Array response to nHz GWs from OJ 287’s SMBH binary
engine

The above listed observational campaigns in many electromagnetic observational
windows clearly point to the presence of a nHz GW emitting SMBH binary in OJ
287. Therefore, it is important to develop prescriptions to model the expected
pulsar timing array (PTA) response to nHz inspiral GWs from systems like OJ
287 (Burke-Spolaor et al., 2018). The fact that the orbital eccentricity and
primary BH spin can substantially influence the BH binary orbital evolution,
as demonstrated in (Dey et al., 2018; Klein et al., 2018), prompted us to
employ an improved version of GW phasing approach to model temporally evolving
GW polarization states from such systems Damour et al. (2004); Susobhanan et
al. (2020); Dey et al. (2021). This is required as the PTA relevant GW-induced
pre-fit pulsar timing residual, namely PTA response to inspiral GWs, induced
by an isolated GW source at a given epoch $t$ may be written as (Anholm et
al., 2009)

$\displaystyle
R(t)=\int_{0}^{t}\left[h(t^{\prime}-\tau_{P})-h(t^{\prime})\right]\
dt^{\prime}\,,$ (9)

where $\tau_{P}$ is the geometric time delay in the solar system barycenter
(SSB) frame of the underlying millisecond pulsar while $t$ and $t^{\prime}$
are the SSB relevant coordinate time variables. Further, $h(t)$ provides the
effective GW strain induced by an isolated source at an epoch $t$ and depends
on the intrinsic parameters of the GW source like the total mass, spin values,
orbital eccentricity and angles that specify its sky location, as detailed in
Susobhanan et al. (2020). In other words, PTA response of a pulsar to inspiral
GWs from OJ 287 like system at a time $t$ depends not only on the strength of
the inspiral GWs at the SSB at $t$ but also on the GW amplitude at location of
the pulsar at time $t-\tau_{P}$. This is another reason why we need to employ
PN approximation to track the orbital evolution and the resulting GW
polarization states while modeling PTA response to inspiral GWs from SMBH
binaries (Susobhanan et al., 2020).

In an ongoing effort, we are modeling temporally evolving quadrupolar order
$h(t)$ from comparable mass BH binaries inspiraling along PN-accurate
eccentric orbits under the influence of dominant order spin-orbit interactions
due to the primary BH spin. This is being done by adapting and extending the
GW phasing efforts of Damour et al. (2004); Susobhanan et al. (2020). It may
be recalled that this approach usually employs Keplerian-type parametric
solution to model the conservative dynamics of SMBH binaries in eccentric
orbits. Thereafter, the effects of GW emission are incorporated by adapting
the method of variation of constants, as detailed in (Damour et al., 2004). In
practice, we invoke the arguments of energy and angular momentum balance that
lead to differential equations for orbital frequency and eccentricity with the
help of orbital averaged PN-accurate far-zone energy and angular momentum
fluxes Dey et al. (2021).

Figure 14: Timing residuals induced by inspiral GWs from OJ 287’s SMBH binary
on a typical PTA pulsar PSR J0751+1807. We let the pulsar distance to be 1.4
Kpc Perera et al. (2019) while SMBH binary parameters as listed in Section 1.1
with its luminosity distance to be 1.6 Gpc. Timing of milli-second pulsars
with SKA are expected to achieve similar timing residuals.

In Figure 14, we plot the expected PTA response of PSR J0751+1807 to inpsiral
nHz GWs associated with our SMBH binary central engine description for OJ 287
as detailed in (Dey et al., 2018). The plot is based on Equation 9 where we
have employed quadrupolar order $h(t)$ expression for spinning BH binaries in
eccentric orbits from Königsdörffer & Gopakumar (2005) while the orbital
parameters are as detailed in Dey et al. (2018). The sharp features depend on
the orbital eccentricity and pulsar distance and provide unique signatures of
the expected PTA response to nHz GWs from such SMBH binaries. Unfortunately,
the expected timing residuals are substantially below what we can achieve with
the current timing campaigns on PSR J0751+1807. However, it should be possible
to find and time milli-second pulsars with such exquisite timing precision
during the SKA era Feng et al. (2020). Therefore, it should be possible to
monitor persistent inspiral GWs from SMBH binaries like OJ 287 during the SKA
era. This development should naturally allow us to pursue persistent multi-
messenger GW astronomy with crucial inputs from the ongoing, planned and
proposed multi-wavelength electromagnetic observations of OJ 287. This is
because such observational campaigns should allow us to constrain various
astrophysical parameters of the binary BH central engine description for the
unique blazar OJ 287.

##  11 Discussion

This article emphasizes the importance of continuing monitoring of OJ 287 at
all wavelengths of the electromagnetic spectrum as well as in gravitational
waves. Even though the binary orbit has now been solved (Dey et al., 2018) and
confirmed (Laine et al., 2020), there are many aspects of the system which
still need verification and from which we can learn interesting astrophysics.

(1) The confirmation of the 2022 optical flare: This flare does not appear in
any of the simple schemes that have been proposed over the years as
alternatives to the LV96 binary black hole model. Difficult as it is to
observe from the ground, we should just be able to observe the tail end of
this flare. Polarisation observations play an important role here since the
tail end of the flare should very highly polarised (Smith et al., 1985;
Villforth et al., 2010; Valtonen et al., 2016).

(2) The detection of the signals from the December 3, 2021 disk impact: This
is energetically a huge event which has at previous times been studied only by
occasional photometric measurements. This time we expect a full coverage from
optical polarisation to X-rays which are the expected key indicators of the
impact. So far, in ongoing MOMO/Swift and ground-based monitoring
observations, OJ 287 was found in a low-state in early December 2021 (Komossa
et al., and Zola et al., in prep.). Also the timing of the impact depends on
the detailed model which will be improved after the 2021/22 observing season.

(3) The $H\alpha$ flares: This has been caught only once, accidentally by
Sitko & Junkkarinen (1985). Later spectroscopic surveys, some of which are
reported in this paper for the first time, have not been equally lucky. Now we
give a rather precise prediction for the $H\alpha$ flare which follows the
December 3, 2021, disk impact.

(4) The secondary black hole: The secondary black hole has never been directly
detected. One of the problems is the short separation between the two black
holes in the sky, typically 10 microarcseconds, depending on its orbital phase
in an eccentric orbit. It will be several years before this resolution is
achieved in radio VLBI. However, an attempt was made to detect the secondary
by its spectral lines, but so far we only report a null result. It is a
difficult task due to the high background level of the jet component. New
attempts should be made during periods of low light level. Here we have made a
prediction when such occasions may arise.

(5) Radio jet of the secondary: Even though the immediate surroundings of the
secondary are hard to resolve, the current resolution in radio VLBI should be
just enough to see the jet. A strong radio component was recently reported
(Gomez et al., 2021) at the expected position of the secondary jet (Dey et
al., 2021). It remains to be seen if this component behaves as expected. One
problem is that the current direction of the primary jet is not very far from
the suspected secondary jet. However, if the primary jet behaves according to
our models (Dey et al., 2021), the primary jet will move away from the
secondary jet, and the identification becomes more clear. Thus the continuing
radio VLBI monitoring of OJ 287 is essential.

(6) Radio jet of the primary: The radio jet of OJ 287 has been mapped in radio
VLBI since 1980’s. The origin of the jet and its connection to the central
black hole and its accretion disk are still open problems. In the LV96 binary
black hole model it is possible to calculate the direction of the spin axis of
the primary as well as the normal to the innermost disk. They follow each
other more or less (Dey et al., 2021), with the disk normal showing a little
more short time scale wobble than the spin axis. If we make the reasonable
assumption that the initial jet direction is determined either by the spin
axis (spin model) or the disk normal (disk model), we obtain a unique model
for the jet which wobbles predictably. When projected onto the sky, we obtain
position angles at any point in time. And since the resolution at higher
frequencies is better than at lower frequencies in the global VLBI, we get
snapshots of the jet at different distances from the center. After solving for
the unknown viewing angles, we get a full picture of the bulk flow in the jet.

(7) Optical polarisation angle: It is usually assumed that the optical
emission arises in the jet closer to the black hole than the radio emission.
We have applied the same model for the optical polarisation that works for the
radio position angle, and find that the region of optical emission is about
four times closer to the central black hole than the inner 86 GHz emission
knot. It is interesting also that the optical polarisation gives the direction
of the innermost jet. This is not surprising since the jet has been found to
have helical magnetic field geometry (Gomez et al., 2021). When compressed
along the jet axis in shock fronts, the field lies practically perpendicular
to the jet axis and the polarisation is lined up with the jet. It is important
to continue regular polarisation monitoring programs to learn more about the
inner jet.

(8) Optical fades: The times of low light level in the jet are useful also for
the detection of the host galaxy. The absolute magnitude of the host has been
measured by surface photometry out to $\sim 30$ kpc from the galaxy center,
beyond which it is difficult to measure any light above the background
generated by the bright point source, the AGN itself. Its brightness measured
out to this distance is among the highest ever detected, and therefore it is
likely that the total brightness of OJ 287 host is also among the highest. The
brightness may also be measured by multicolor photometry at low light, because
the color of the host is expected to be very different from the color of the
jet. Thus the combined color of the source becomes more and more like galaxy
colors when the source fades. This method suggests a brightness increase by
one magnitude if the outer envelope of the host galaxy is included.

(9) High-energy emission and broad-band SEDs: Very dense, simultaneous
X-ray–UV–optical monitoring has allowed us to identify the latest binary
after-flares in 2017 and 2020, based on the brightest X-ray flares so far
discovered from OJ 287 and other arguments (Komossa et al., 2020a). Ongoing
monitoring will reveal new outstanding X-ray–UV outbursts. Spectroscopy, SED
modelling, cross-correlation analyses and structure function analyses (Komossa
et al., 2021b, d) will allow us to probe in detail the binary and outburst
physics across the electromagnetic spectrum, and will allow to search for an
accretion disk contribution of the secondary at select epochs. Ongoing high-
energy monitoring could also serve as trigger criterion for new VHE
observations.

(10) Very high energy gamma rays: VHE gamma rays have been detected in OJ287
only once. This occurred after the two big tidal flares which followed the
2015 impact flare. Even though we do not know the mechanism which produced
this event, it may useful to monitor OJ 287 again after the tidal flares which
follow the 2022 impact flare. More frequent monitoring would be also useful to
see if the two phenomena are really connected.

(11) Fast variability: The shortest time scale of periodic or quasi-periodic
variability that may be associated with the primary jet is 3.5 days. This has
been seen in many occasions. The time scale of variability should be $\sim
122$ times shorter in the secondary jet, due to the smaller mass of the
secondary by this factor, or it could be even as low as $\sim 15.7$ minutes if
the secondary spins much faster than the primary. Thus the study of fast
variability in radio, optical and X-rays may reveal details of the secondary
jet which are otherwise hard to recover.

12) Pulsar Timing Array: It will be important to develop accurate and
efficient algorithms to model the expected PTA signals from SMBH binaries like
OJ 287. Further, data analysis approaches that can efficiently look for such
nHz GW sources in the PTA data sets of SKA era will need to be developed
Burke-Spolaor et al. (2018). Additionally, it will be desirable to search for
milli-second pulsars close to the Sky location of OJ 287 as they are expected
to be more sensitive to inspiral GWs from OJ 287. Investigations should be
pursued to figure out the possibility of tracking the trajectory of secondary
BH during the ngEHT era.

13) It has been argued that OJ 287 may produce a few muon neutrinos with
energies more than 100 TeV in instruments like IceCube-Gen2 during ten years
of Fermi flaring epochs, influenced by the modeling of the 2017 flare of TXS
0506+056 (Oikonomou et al., 2019; IceCube Collaboration et al., 2018). The
neutrino studies will be an important part of multi-messenger astronomy with
OJ 287 in the coming decades.

M.J.V. and L.D. have been responsible for the main design of this article as
well as for theoretical developments. L.D and A.G. were in charge of the
gravitational wave part, S.Z. of optical photometry, S.K. of the MOMO project
and X-ray parts, T.P. and S.Z. of optical spectroscopy, R.H. of historical
magnitude measurements, H.J. and A.V.B. of optical polarimetry and J.L.G. of
radio mapping.

L.D. and A.G. acknowledge the support of the Department of Atomic Energy,
Government of India, under project identification # RTI 4002. A.G. is grateful
for the financial support and hospitality of the Pauli Center for Theoretical
Studies and the University of Zurich. S.Z. acknowledge the grant NCN
2018/29/B/ST9/01793. R.H. acknowledges support by GACR grant 13-33324S and by
the MSMT Mobility project 8J18AT036.

###### Acknowledgements.

We thank Kari Nilsson for a valuable contribution in data reduction and
producing one of the figures. We also appreciate the help that Vilppu Piirola
has provided and continues to provide in optical polarisation observations of
OJ 287. We also thank Gene Byrd and Shirin Haque for reading the manuscript
and for valuable inputs. The authors declare no conflict of interest.
References

## References

  * Begelman et al. (1980) Begelman, M.C., Blandford, R.D. & Rees, M.J. Massive black hole binaries in active galactic nuclei. Nature 1980, 287, 307
  * Valtaoja et al. (1989) Valtaoja, L., Valtonen, M.J. & Byrd, G.G. Binary Pairs of Supermassive Black Holes: Formation in Merging Galaxies. ApJ 1989, 343, 47
  * Mikkola & Valtonen (1992) Mikkola, S. & Valtonen, M.J. Evolution of binaries in the field of light particles and the problem of two black holes. MNRAS 1992, 259, 115
  * Valtonen (1996) Valtonen, M.J. Triple black hole systems formed in mergers of galaxies. MNRAS 1996, 278, 186
  * Quinlan (1996) Quinlan, G.D. The dynamical evolution of massive black hole binaries I. Hardening in a fixed stellar background. NewAstr 1996, 1, 35
  * Milosavljevic & Merritt (2001) Milosavljevic, M. & Merritt, D. Formation of Galactic Nuclei. ApJ 2001, 563, 34
  * Volonteri et al. (2003) Volonteri, M., Haardt, F. & Madau, P. The Assembly and Merging History of Supermassive Black Holes in Hierarchical Models of Galaxy Formation. ApJ 2003, 582, 559
  * Sesana et al. (2004) Sesana, A., Haardt, F., Madau, P.& Volonteri, M. Low-Frequency Gravitational Radiation from Coalescing Massive Black Hole Binaries in Hierarchical Cosmologies. ApJ 2004, 611, 623
  * Burke-Spolaor et al. (2018) Burke-Spolaor, S., Blecha, L., Bogdanovic, T., Comerford, J.M., Lazio, J., Liu, X., Maccarone, T.J., Pesce, D., Shen, Y.& Taylor, G. Supermassive Black Hole Pairs and Binaries. ASPC 2018, 517, 677
  * Burke-Spolaor et al. (2018) Burke-Spolaor, S., Taylor, S.R., Charisi, M., Dolch, T., Hazboun, J.S., Holgado, A.M., Kelley, L.Z., Lazio, T.J.W., Madison, D.R., McMann, N. et al. The astrophysics of nanohertz gravitational waves. A&ARv 2019, 27, 5
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P.T., Blumer, H., Becsy, B., Brazier, A., Brook, P.R., Burke-Spolaor, S., Chatterjee, S., Chen, S., Cordes, J.M. et al. The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background. ApJ 2020, 905, L34
  * Goncharov et al. (2021) Goncharov, B., Shannon, R.M., Reardon, D.J., Hobbs, G., Zic, A., Bailes, M., Curylo, M., Dai, S., Kerr, M., Lower, M. E. et al. On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array. ApJ 2021, 917, L19
  * Chen et al. (2021) Chen, S., Caballero, R. N., Guo, Y. J., Chalumeau, A., Liu, K., Shaifullah, G., Lee, K.J., Babak, S., Desvignes, G., Parthasarathy, A. et al. Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search. MNRAS 2021, 508, 4970
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Abraham, S., Acernese, F., Ackley, K., Adams, A., Adams, C., Adhikari, R.X., Adya, V.B., Affelt, C. et al. GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run. Physical Review X 2021, 11, 021053
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run. 2021, arXiv:2111.03606
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., Ackley, K., Adams, C., Adams, T., Addesso, P., Adhikari, R.X., Adya, V.B. et al. GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral. Phys. Rev. Lett. 2017, 119, 161101
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., Ackley, K., Adams, C., Adams, T., Addesso, P., Adhikari, R.X., Adya, V.B. et al. Multi-messenger Observations of a Binary Neutron Star Merger. ApJ 2017, 848, L12
  * Xin et al. (2021) Xin, C., Mingarelli, C. M. F., & Hazboun, J. S. Multimessenger Pulsar Timing Array Constraints on Supermassive Black Hole Binaries Traced by Periodic Light Curves. ApJ 2021, 915, 97
  * Sillanpää et al. (1988) Sillanpää, A., Haarala, S., Valtonen, M.J., Sundelius, B. & Byrd, G.G. OJ 287: Binary Pair of Supermassive Black Holes. ApJ 1988, 325, 628
  * Sudou et al. (2003) Sudou, H., Iguchi, S., Murata, Y.& Taniguchi, Y. Orbital Motion in the Radio Galaxy 3C 66B: Evidence for a Supermassive Black Hole Binary. Science 2003, 300, 1263
  * Jenet et al. (2004) Jenet, F. A., Lommen, A., Larson, S. L. & Wen, L. Constraining the Properties of Supermassive Black Hole Systems Using Pulsar Timing: Application to 3C 66B. ApJ 2004, 606, 799
  * Iguchi et al. (2010) Iguchi, S., Okuda, T., & Sudou, H. A Very Close Binary Black Hole in a Giant Elliptical Galaxy 3C 66B and its Black Hole Merger. ApJ 2010, 724, L166
  * Liu et al. (2014) Liu, F. K., Li, S., & Komossa, S. A Milliparsec Supermassive Black Hole Binary Candidate in the Galaxy SDSS J120136.02+300305.5. ApJ 2014, 786, 103
  * Graham et al. (2015) Graham, M. J., Djorgovski, S. G., Stern, D., Glikman, E., Drake, A.J., Mahabal, A.A., Donalek, C., Larson, S. & Christensen, E. A possible close supermassive black-hole binary in a quasar with optical periodicity. Nature 2015, 518, 74
  * Zhu & Thrane (2020) Zhu, X.-J. & Thrane, E. Toward the Unambiguous Identification of Supermassive Binary Black Holes through Bayesian Inference. ApJ 2020, 900, 117
  * Charisi et al. (2016) Charisi, M., Bartos, I., Haiman, Z., Price-Whelan, A.M.,Graham, M.J., Bellm, E.C., Laher, R.R. & Marka, S. A population of short-period variable quasars from PTF as supermassive black hole binary candidates. MNRAS 2016, 463, 2145
  * Zhu et al. (2019) Zhu, X.-J., Cui, W., & Thrane, E. The minimum and maximum gravitational-wave background from supermassive binary black holes. MNRAS 2019, 482, 2588
  * Dey et al. (2019) Dey, L., Gopakumar, A., Valtonen, M., Zola, S., Susobhanan, A., Hudec, R., Pihajoki, P., Pursimo, T., Berdyugin, A., Piirola, V. et al. The Unique Blazar OJ 287 and Its Massive Binary Black Hole Central Engine. Universe 2019, 5, 108
  * Feng et al. (2020) Feng, Y., Li, D., Zheng, Z. & Tsai, C.-W. Supermassive binary black hole evolution can be traced by a small SKA pulsar timing array. Phys.Rev. D 2020, 102, 023014
  * Dickel et al. (1967) Dickel J.R., Yang, K.S., McVittie, G.C. & Swenson, G.W., Jr. A survey of the sky at 610.5 MHz. II. The region between declinations +15 and +22 degrees. AJ 1967, 72, 757
  * Hudec et al. (2013) Hudec R., Basta, M., Pihajoki, P. & Valtonen, M. The historical 1900 and 1913 outbursts of the binary blazar candidate OJ287. A&A 2013, 559, 20
  * Valtonen et al. (2006) Valtonen, M.J., Lehto, H.J., Sillanpää, A., Nilsson, K., Mikkola, S., Hudec, R., Basta, M., Teräsranta, H., Haque, S. & Rampadarath Predicting the Next Outbursts of OJ 287 in 2006-2010. ApJ 2006, 646, 36-48
  * Dey et al. (2018) Dey, L., Valtonen, M.J., Gopakumar, A., Zola, S., Hudec, R., Pihajoki, P., Ciprini, S., Matsumoto, K., Sadakane, K., Kidger, M. et al. Authenticating the Presence of a Relativistic Massive Black Hole Binary in OJ 287 Using Its General Relativity Centenary Flare: Improved Orbital Parameters. ApJ 2018, 866, 11
  * Kidger (2007) Kidger, M. Cosmological Enigmas: Pulsars, Quasars, and Other Deep-Space Questions. The Johns Hopkins University Press 2007
  * Smith et al. (1985) Smith, P.S., Balonek, T.J., Heckert, P.A., Elston, R. & Schmidt, G. D. UBVRI field comparison stars for selected active quasars and BL Lacertae objects. AJ 1985, 90, 1184-1187
  * Sillanpää et al. (1985) Sillanpää, A., Teerikorpi, P., Haarala, S., Korhonen, T., Efimov, I. S., & Shakhovskoi, N. M. Similar structures in the outbursts of OJ 287 in 1972 and 1983. A&A 1985, 147, 67-70
  * Sillanpää et al. (1996a) Sillanpää, A., Takalo, L.O., Pursimo, T., Lehto, H.J., Nilsson, K., Teerikorpi, P., Heinämäki, P., Kidger, M., de Diego, J.A., Gonzalez-Perez, J.N. et al. Confirmation of the 12-year optical outburst cycle in blazar OJ 287. A&A 1996, 305, L17
  * Lehto & Valtonen (1996) Lehto, H.J., & Valtonen, M.J. OJ 287 Outburst Structure and a Binary Black Hole Model. ApJ 1996, 460, 207
  * Valtonen (1996) Valtonen, M.J. The OJ 287 binary model and the expected outburst in November 1995. Workshop on Two Years of Intensive Monitoring of OJ 287 and 3C 66A 1996, Proc. meeting at Oxford, England, 11-14 September, 1995. Ed. Leo O. Takalo, U.Turku, 64
  * Sillanpää et al. (1996b) Sillanpää, A., Takalo, L.O., Pursimo, T., Nilsson, K., Heinämäki, P., Katajainen, S., Pietilä, H., Hanski, M., Rekola, R., Kidger, M. et al. Double-peak structure in the cyclic optical outbursts of blazar OJ 287. A&A 1996, 315, L13-L16
  * Ivanov et al. (1998) Ivanov, P.B., Igumenshchev, I.V. & Novikov, I.D. Hydrodynamics of Black Hole-Accretion Disk Collision. ApJ 1998, 507, 131-144
  * Pihajoki (2016) Pihajoki, P. Black hole accretion disc impacts. MNRAS 2016, 457, 1145-1161
  * Sundelius et al. (1996) Sundelius, B., Wahde, M., Lehto, H.J. & Valtonen, M.J. Long-time Brightness Variations of 0J287 in the Binary Black Hole Model. ASP-CS 1996, 110, 99
  * Sundelius et al. (1997) Sundelius, B., Wahde, M., Lehto, H.J. & Valtonen, M.J. A Numerical Simulation of the Brightness Variations of OJ 287. ApJ 1997, 484, 180-185
  * Valtonen (2007) Valtonen, M.J. New Orbit Solutions for the Precessing Binary Black Hole Model of OJ 287. ApJ 2007, 659, 1074-1081
  * Valtonen et al. (2006a) Valtonen, M.J., Nilsson, K., Sillanpää, A., Takalo, L.O., Lehto, H.J., Keel, W.C., Haque, S., Cornwall, D. & Mattingly, A. The 2005 November Outburst in OJ 287 and the Binary Black Hole Model. ApJ 2006, 643, L9-L12
  * Valtonen & Sillanpää (2011) Valtonen, M.J. & Sillanpää 2005-2010 Multiwavelength Campaing of OJ287. AcPol 2011, 51, 76
  * Laine et al. (2020) Laine, S., Dey, L., Valtonen, M., Gopakumar, A., Zola, S., Komossa, S., Kidger, M., Pihajoki, P., Gomez, J.L., Caton, D. et al. Spitzer Observations of the Predicted Eddington Flare from Blazar OJ 287. ApJL 2020, 894, L1
  * Valtonen et al. (2008) Valtonen, M.J., Lehto, H.J., Nilsson, K., Heidt,J., Takalo, L.O., Sillanpää, A., Villforth, C., Kidger, M., Poyner, G., Pursimo, T. et al. A massive binary black-hole system in OJ287 and a test of general relativity. Nature 2008, 452, 851-853
  * Valtonen et al. (2016) Valtonen, M. J., Zola, S., Ciprini, S., Gopakumar, A., Matsumoto, K., Sadakane, K., Kidger, M., Gazeas, K., Nilsson, K., Berdyugin, A. et al. (OJ287-15/16 Collaboration). Primary Black Hole Spin in OJ 287 as Determined by the General Relativity Centenary Flare. ApJ 2016, 819, L37
  * Valtonen et al. (2012) Valtonen, M. J., Ciprini, S. & Lehto, H. J. On the masses of OJ287 black holes. MNRAS 2012, 427, 77-83
  * Kidger et al. (2018) Kidger, M., Zola, S., Valtonen, M., Lähteenmäki, A., Järvelä, E., Tornikoski, M., Tammi, J., Liakos, A. & Poyner, G. Far-infrared photometry of OJ 287 with the Herschel Space Observatory. A&A 2018, 610, A74
  * Villata et al. (1998) Villata, M., Raiteri, C.M, Sillanpää, A. & Takalo, L.O. A beaming model for the OJ 287 periodic optical outbursts. MNRAS 1998, 293, L13-L16
  * Rieger (2004) Rieger, F.M. On the Geometrical Origin of Periodicity in Blazar-type Sources. ApJ 2004, 615, L5-L8
  * Ciprini et al. (2015) Ciprini, S., Perri, M., Verreccha, F. & Valtonen, M. Fermi-LAT detection of hard spectrum and enhanced gamma-ray emission from the BL Lac object PKS 1717+177. ATel 2015, 8401, 1
  * MacPherson et al. (2015) MacPherson, E., Ister, J.C., Urry, M., Coppi, P., Bailyn, C. & Dincer, T. SMARTS Enhanced Optical & Infrared Activity in Blazar OJ 287. ATel 2015, 8392, 1
  * Shappee et al. (2015) Shappee, B.J., Stanek, K.Z., Holoien, T.W.-S., Brown, J.S., Kochanek, C.S., Godoy-Rivera, D., Basu, U., Prieto, J.L., Bersier, D., Dong, S. et al. Strong Optical Flare from Blazar OJ 287 Detected by ASAS-SN. ATel 2015, 8372, 1
  * Pursimo et al. (2000) Pursimo T., Takalo, L.O., Sillanpää, A., Kidger, M., Lehto, H.J., Heidt, J., Charles, P.A., Aller, H., Aller, M., Beckmann, V. et al. Intensive monitoring of OJ 287. A&AS 2000, 146, 141-155
  * Villforth et al. (2010) Villforth C., Nilsson, K., Heidt, J., Takalo, L.O., Pursimo, T., Berdyugin, A., Lindfors, E., Pasanen, M., Winiarski, M., Drozdz, M. et al. Variability and stability in blazar jets on time-scales of years: optical polarization monitoring of OJ 287 in 2005-2009. MNRAS 2010, 402, 2087-2111
  * Pursimo et al. (2021) Pursimo T., Takalo, L.O., Sillanpää, A., Kidger, M., Lehto, H.J., Heidt, J., Charles, P.A., Aller, H., Aller, M., Beckmann, V. et al. VizieR Online Data Catalog: Intensive monitoring of OJ 287 (Pursimo+, 2000). VizieR On-line Data Catalog: J/A+AS/ 2000, 146, 141
  * Valtonen et al. (2017) Valtonen, M. J., Zola, S., Jermak, H., Ciprini, S., Hudec, R., Dey, L., Gopakumar, A., Reichart, D., Caton, D., Gazeas, K. et al. Polarization and Spectral Energy Distribution in OJ 287 during the 2016/17 Outbursts. Galaxies 2017, 5, 83
  * Zola et al. (2016) Zola, S., Valtonen, M. J., Bhatta, G., Goyal, A., Debski, B., Baran, A., Krzesinski, J., Siwak, M., Ciprini, S., Gopakumar, A. et al. A Search for QPOs in the Blazar OJ287: Preliminary Results from the 2015/2016 Observing Campaign. Galaxies 2016, 4, 41
  * Abdo et al. (2009) Abdo A.A., Ackermann, M., Ajello, M., Atwood, W.B., Axelsson, M., Baldini, L., Ballet, J., Barbiellini, G., Bastieri, D., Baughman, B.M. et al. Bright Active Galactic Nuclei Source List from the First Three Months of the Fermi Large Area Telescope All-Sky Survey. ApJ 2009, 700, 597-622
  * D’Arcangelo et al. (2009) D’Arcangelo F.D., Marscher A.P., Jorstad S., Smith, P.S., Larionov, V.M., Hagen-Thorn, V.A., Williams, G.G., Gear, W.K., Clemens, D.P., Sarcia, D. et al. Synchronous Optical and Radio Polarization Variability in the Blazar OJ287. ApJ 2009, 697, 985-995
  * Sitko & Junkkarinen (1985) Sitko, M.L., & Junkkarinen, V. T. Continuum and line fluxes of OJ 287 at minimum light. PASP 1985, 97, 1158-1162
  * Nilsson et al. (2010) Nilsson, K., Takalo, L.O., Lehto & Sillanpää, A. H-alpha monitoring of OJ 287 in 2005-08. A&A 2010, 516, A60
  * Padovani & Giommi (1995) Padovani P.& Giommi P. The Connection between X-Ray– and Radio-selected BL Lacertae Objects. ApJ 1995, 444, 567
  * Jorstad et al. (2005) Jorstad, S.G., Marscher, A.P., Lister, M.L., Stirling, A.M., Cawthorne, T.V., Gear, W.K., Gomez, J.L., Stevens, J.A., Smith, P.S., Forster, J.R & Robson, E.I. Polarimetric Observations of 15 Active Galactic Nuclei at High Frequencies: Jet Kinematics from Bimonthly Monitoring with the Very Long Baseline Array. AJ 2005, 130, 1418-1465
  * Lee et al. (2020) Lee, J.W., Lee, S.S., Algaba, J.-C., Hodgson, J., Kim, J.-Y., Park, J., Kino, M., Kim, D.-W., Kang, S., Yoo, S. et al. Interferometric Monitoring of Gamma-Ray Bright AGNs: OJ 287. ApJ 2020, 902, 104
  * Gomez et al. (2021) Gomez, J.L., Traianou, E., Krichbaum, T.P., Lobanov, A., Fuentes, A., Lico, R., Zhao, G.-Y., Bruni, G., Kovalev, Y.Y., Lähteenmäki, A. et al. Probing the innermost regions of AGN jets and their magnetic fields with RadioAstron. V. Space and ground millimeter-VLBI imaging of OJ 287. arXiv:2111.11200 2021
  * Goddi et al. (2021) Goddi C., Marti-Vidal, I., Messias, H., Bower, G.C., Broderick, A.E., Dexter, J., Marrone, D.P., Moscibrodzka, M., Nagai, H., Algaba, J.C. et al. Polarimetric Properties of Event Horizon Telescope Targets from ALMA. ApJL 2021, 910, L14
  * Hodgson et al. (2017) Hodgson J. A., Krichbaum, T.P., Marscher, A.P., Jorstad, S.G., Rani, B., Marti-Vidal, I., Bach, U., Sanchez, S., Bremer, M., Lindqvist, M. et al. Location of $\gamma$-ray emission and magnetic field strengths in OJ 287. A&A 2017, 597, A80
  * Agudo et al. (2012) Agudo, I., Marscher, A.P., Jorstad, S.G., Gomez, J.L., Perucho, M., Piner, P.G., Rioja, M. & Dodson, R. Erratic Jet Wobbling in the BL Lacertae Object OJ287 Revealed by Sixteen Years of 7 mm VLBA Observations. ApJ 2011, 747, 63
  * Valtonen & Pihajoki (2013) Valtonen, M. J. & Pihajoki, P. A helical jet model for OJ287. A&A 2013, 557, A28
  * Dey et al. (2021) Dey, L., Valtonen, M.J., Gopakumar, A., Lico, R., Gomez, J., Susobhanan, A., Komossa, S. & Pihajoki, P. Explaining temporal variations in the jet PA of the blazar OJ 287 using its BBH central engine model. MNRAS 2021, 503, 4400-4412
  * Shrader et al. (1996) Shrader C.R., Hartman R.C., Webb J.R., Probable detection of high-energy gamma-ray emission from OJ 287 during a major optical flare. A&AS 1996, 120, 599-602
  * Agudo et al. (2011) Agudo, I., Jorstad, S.G., Marscher, A.P., Larionov, V.M., Gomez, J.L., Lähteenmäki, A., Gurwell, M., Smith, P.S., Wiesemeyer, H., Thum, C. et al. Location of $\gamma$-ray Flare Emission in the Jet of the BL Lacertae Object OJ287 More than 14 pc from the Central Engine. ApJ 2011, 726, L13
  * O’Brien, S. et al. (2017) O’Brien, S. et al. for the VERITAS Collaboration VERITAS detection of VHE emission from the optically bright quasar OJ 287. arXiv 2017, 1708.02160
  * Komossa et al. (2017) Komossa S., Grupe D., Schartel N., Gallo L., Gomez J.L., Kollatschny W., Kriss G., Leighly K., Longinotti A.L., Parker, M., Santos-Lleo M., Wilkins D. & Zetzl M. The Extremes of AGN Variability. IAUS: New Frontiers in Black Hole Astrophysics 2017, 324, 168-171
  * Komossa et al. (2020a) Komossa S., Grupe D., Parker M.L., Valtonen M. J., Gómez J.L., Gopakumar A. & Dey L. The 2020 April-June super-outburst of OJ 287 and its long-term multiwavelength light curve with Swift: binary supermassive black hole and jet activity. MNRAS 2020, 498, L35
  * Komossa et al. (2021d) Komossa S., Grupe D., Gallo L.C., Gonzalez A., Yao S., Hollett A.R., Parker M.L. & Ciprini S. MOMO IV: The complete Swift X-ray and UV/optical light curve and characteristic variability of the blazar OJ 287 during the last two decades. ApJ 2021, 923, 51
  * Madejski & Schwartz (1988) Madejski G.M.& Schwartz, D.A. Studies of BL Lacertae Objects with the Einstein Observatory: The Soft X-Ray Spectra of OJ 287 and PKS 0735+178. ApJ 1988, 330, 776
  * Komossa et al. (2021a) Komossa S., Grupe D., Parker M.L., Gómez J.L., Valtonen M.J., Nowak M.A., Jorstad S.G., Haggard D., Chandra, S., Ciprini, S. et al. 2021a, MNRAS, 504, 5575
  * Marscher & Jorstad (2011) Marscher A.P. & Jorstad S.G. 2011, ApJ, 729, 26
  * Komossa et al. (2021c) Komossa S., Ciprini S., Dey L., Gallo L.C., Gomez J.L., Gonzalez A., Grupe D., Kraus A., Laine, S., Parker, M.L. et al. Supermassive Binary Black Holes and the Case of OJ 287. Publ. Astron. Obs. Belgrade 2021, 100, 29-42
  * Komossa et al. (2021b) Komossa S., Grupe D., Kraus A., Gallo L.C., Gonzalez A.G., Parker M.L., Valtonen M.J., Hollett A.R., Bach, U., Gomez, J.L. et al. Project MOMO: Multiwavelength Observations and Modeling of OJ 287. Universe 2021, 7, 261
  * Valtonen et al. (2011) Valtonen, M.J., Lehto, H.J., Takalo, L.O. & Sillanpää, A. Testing the 1995 Binary Black Hole Model of OJ287. ApJ 2011, 729, 33
  * Valtonen et al. (2019) Valtonen, M. J., Zola, S., Pihajoki, P., Enestam, S., Lehto, H.J., Dey, L., Gopakumar, A., Drozdz, M., Ogloza, W., Zejmo, M. et al. Accretion Disk Parameters Determined from the Great 2015 Flare of OJ 287. ApJ 2019, 882, 88
  * Pihajoki et al. (2013) Pihajoki, P., Valtonen, M. & Ciprini, S. Short time-scale periodicity in OJ 287. MNRAS 2013, 434, 3122-3129
  * Takalo et al. (1990) Takalo, L. O., Kidger, M., de Diego, J. A., Sillanpää, A., Piirola, V., Teräsranta, H. A sudden fade of OJ 287. A&AS 1990, 83, 459
  * Valtonen et al. (2021) Valtonen, M. J., Zola, S., Ciprini, S., Kidger, T.,Pursimo, T., Gopakumar, A., Matsumoto, K., Sadakane, K., Caton, D.B., Nilsson, K. et al. Host galaxy magnitude of OJ287 from its colours at minimum light 2021, submitted to MNRAS
  * Saglia et al. (2016) Saglia, R. P., Opitsch, M., Erwin, P., Thomas, J., Beifiori, A., Fabricius, M., Mazzalay, X., Nowak, N., Rusli, S. P. & Bender, R., The SINFONI Black Hole Survey: The Black Hole Fundamental Plane Revisited and the Paths of (Co)evolution of Supermassive Black Holes and Bulges. ApJ 2016, 818, 47
  * Nilsson et al. (2020) Nilsson, K., Kotilainen, J., Valtonen, M., Gomez, J.L., Castro-Tirado, A.J., Gopakumar, A., Jeong, S., Kidger, M., Komossa, S., Mathur, S. et al. The Host Galaxy of OJ 287 Revealed by Optical and Near-infrared Imaging. ApJ 2020, 904, 102
  * Sandage (1972) Sandage, A. Absolute Magnitudes of E and so Galaxies in the Virgo and Coma Clusters as a Function of U - B Color. ApJ 1972, 176, 21
  * Graham & Scott (2013) Graham, A. & Scott, N., The M BH-L spheroid Relation at High and Low Masses, the Quadratic Growth of Black Holes, and Intermediate-mass Black Hole Candidates. ApJ 2013, 764, 151
  * De Vaucouleurs & De Vaucouleurs (1970) De Vaucouleurs, G., De Vaucouleurs, A., Photometry of Intergalactic Matter in the Coma Cluster. ApLett 1970, 5, 219
  * Läsker et al. (2014) Läsker, R., Ferrarese, L. & van de Ven, G. Supermassive Black Holes and Their Host Galaxies. II. The Correlation with Near-infrared Luminosity Revisited. ApJ 2014, 780, 70
  * Huang et al. (2013) Huang, S., Ho, L. C., Peng, C. Y., Li, Z.-Y. & Barth, A. J. The Carnegie-Irvine Galaxy Survey. III. The Three-component Structure of Nearby Elliptical Galaxies. ApJ 2013, 766, 47
  * Yanny et al. (1997) Yanny, B., Jannuzi, B.T. & Impey, C. Hubble Space Telescope Imaging of the BL Lacertae Object OJ 287. ApJ 1997, 484, L113-L116
  * Portinari et al. (2012) Portinari, L., Kotilainen, J., Falomo, R. & Decarli, R. On the cosmological evolution of the black hole-host galaxy relation in quasars. MNRAS 2012, 420, 732-744
  * Pietilä et al. (1999) Pietilä, H., Takalo, L.O., Tosti, G., Benitez, E., Chiattelli, B., Corradi, R.L.M., Cox, G., de Diego, J.A., de Francesco, G., Dultzin-Hacyan, D. et al. OJ 287 and the predicted fade of 1998. A&A 1999, 345, 760-768
  * Valtonen et al. (2012a) Valtonen, M. J., Gopakumar, A., Mikkola, S., Wiik, K. & Lehto, H. J. Black hole binary OJ287 as a testing platform for general relativity. arXiv 2012, 1208.4524
  * Savolainen et al. (2010) Savolainen, T., Homan, D.C., Hovatta, T., Kadler, M., Kovalev, Y.Y., Lister, M.L., Ros, E. & Zensus, J.A. Relativistic beaming and gamma-ray brightness of blazars. A&A 2010, 512, A24
  * Hovatta et al. (2009) Hovatta, T., Valtaoja, E., Tornikoski, M. & Lähteenmäki, A. Doppler factors, Lorentz factors and viewing angles for quasars, BL Lacertae objects and radio galaxies. A&A 2009, 494, 527-537
  * Valtonen et al. (2012) Valtonen, M. J. & Wiik, K. Optical polarization angle and VLBI jet direction in the binary black hole model of OJ287. MNRAS 2012, 421, 1861-1867
  * Heinzeller & Duschl (2007) Heinzeller, D. & Duschl, W.J. On the Eddington limit in accretion discs. MNRAS 2007, 374, 1146-1154
  * Ghosh & Abramowicz (1997) Ghosh, P. & Abramowicz, M.A. Electromagnetic extraction of rotational energy from disc-fed black holes: the strength of the Blandford-Znajek process. MNRAS 1997, 292, 887-895
  * Chen et al. (2015) Chen, Y.-Y., Zhang, X., Xiong, D. & Yu, X. Black Hole Mass, Jet Power, and Accretion in AGNs. AJ 2015, 150, 8
  * Valtaoja et al. (1985) Valtaoja, E., Lehto, H., Teerikorpi, P., Korhonen, T., Valtonen, M., Teräsranta, H., Salonen, E., Urpo, S., Tiuri, M., Piirola, V. & Saslaw, W.C. A 15.7-min periodicity in OJ287. Nature 1985, 314, 148-149
  * Pihajoki et al. (2013) Pihajoki, P., Valtonen, M., Zola, S., Liakos, A., Drozdz, M., Winiarski, M., Ogloza, W., Koziel-Wierzbowska, D., Provencal, J., Nilsson, K. et al. Precursor Flares in OJ 287. ApJ 2013, 764, 5
  * Valtonen et al. (2010) Valtonen, M. J., Mikkola, S., Merritt, D., Gopakumar, A., Lehto, H.J., Hyvönen, T., Rampadarath, H., Saunders, R., Basta, M. & Hudec, R. Measuring the Spin of the Primary Black Hole in OJ287. ApJ 2010, 709, 725
  * Krolik et al. (2019) Krolik, J. H., Volonteri, M., Dubois, Y. & Devriendt, J. Population Estimates for Electromagnetically Distinguishable Supermassive Binary Black Holes. ApJ 2019, 879, 110
  * Komossa & Zensus (2016) Komossa S., Zensus J.A. Radio-loud narrow-line Seyfert 1 galaxies with high-velocity outflows. IAUS:Star Clusters and Black Holes in Galaxies across Cosmic Time 2016, 312, 63-65
  * Klein et al. (2018) Klein, A., Boetzel, Y., Gopakumar, A., Jetzer, P. & de Vittori, L. Fourier domain gravitational waveforms for precessing eccentric binaries. Phys. Rev. D 2018, 98, 104043
  * Damour et al. (2004) Damour, T., Gopakumar, A., & Iyer, B. R. Phasing of gravitational waves from inspiralling eccentric binaries. Phys. Rev. D 2004, 70, 064028
  * Susobhanan et al. (2020) Susobhanan, A., Gopakumar, A., Hobbs, G. & Taylor, S.R. Pulsar timing array signals induced by black hole binaries in relativistic eccentric orbits. Phys. Rev. D 2020, 101, 043022
  * Dey et al. (2021) Dey, L., Susobhanan, A., Gopakumar, A., Valtonen, M., et al. Pulsar timing array signals due to inspiralling spinning massive black hole binaries in relativistic eccentric orbits. in preparation
  * Anholm et al. (2009) Anholm, M., Ballmer, S., Creighton, J. D. E., Price, L.R. & Siemens, X. Optimal strategies for gravitational wave stochastic background searches in pulsar timing data. Phys. Rev. D 2009, 79, 084030
  * Perera et al. (2019) Perera, B. B. P., DeCesar, M. E., Demorest, P. B., Kerr, M., Lentati, L., Nice, D.C., Oslowski, S., Ransom, S.M., Keith, M.J., Arzoumanian, Z. et al. The International Pulsar Timing Array: second data release. MNRAS 2019, 490, 4666-4687
  * Königsdörffer & Gopakumar (2005) Königsdörffer, C. & Gopakumar, A. Post-Newtonian accurate parametric solution to the dynamics of spinning compact binaries in eccentric orbits: The leading order spin-orbit interaction. Phys. Rev. D 2005, 71, 024039.
  * Oikonomou et al. (2019) Oikonomou F., Murase K., Padovani P., Resconi E.& Mészáros P. High-energy neutrino flux from individual blazar flares. MNRAS 2019, 489, 4347
  * IceCube Collaboration et al. (2018) IceCube Collaboration, Aartsen M. G., Ackermann M., Adams J., Aguilar J. A., Ahlers M., Ahrens M., Al Samarai, I., Altmann, D., Andeen, K., Anderson, T. et al. Multimessenger observations of a flaring blazar coincident with high-energy neutrino IceCube-170922A. Science 2018, 361, eaat1378