MOMO. VI. Multifrequency Radio Variability of the Blazar OJ 287 from 2015 to 2022, Absence of Predicted 2021 Precursor-flare Activity, and a New Binary Interpretation of the 2016/2017 Outburst

Based on our dedicated Swift monitoring program, MOMO, OJ 287 is one of the best-monitored blazars in the X-ray–UV–optical regime. Here, we report results from our accompanying, dense, multifrequency (1.4–44 GHz) radio monitoring of OJ 287 between 2015 and 2022 covering a broad range of activity states. Fermi γ-ray observations were added. We characterize the radio flux and spectral variability in detail, including discrete correlation function and other variability analyses, and discuss its connection with the multiwavelength emission. Deep fades of the radio and optical–UV fluxes are found to occur every 1–2 yr. Further, it is shown that a precursor flare of thermal bremsstrahlung predicted by one of the binary supermassive black hole (SMBH) models of OJ 287 was absent. We then focus on the nature of the extraordinary, nonthermal, 2016/2017 outburst that we initially discovered with Swift. We interpret it as the latest of the famous optical double-peaked outbursts of OJ 287, favoring binary scenarios that do not require a highly precessing secondary SMBH.


INTRODUCTION
Active galactic nuclei (AGN) are powered by accretion of matter onto supermassive black holes (SMBHs) at their centers (Lynden-Bell 1969).A fraction of AGN is radio-loud and harbors powerful, long-lived jets of relativistic particles that are launched in the immediate vicinity of their central SMBHs.The SMBH -accretion disk -jet interface represents one of the most extreme astrophysical environments where magnetic fields, high gas density, and (special and general) relativistic astrophysics all play a crucial role in shaping the multiwavelength electromagnetic emission and structural proper-ties of these systems (see Blandford et al. 2019, for a recent review).
In blazars, the jet is oriented close to the observer's line of sight, and its emission is enhanced by relativistic beaming.The spectral energy distribution (SED) of blazars is characterized by different emission processes at different wavelengths.SEDs show two broad emission humps (Abdo et al. 2010;Ghisellini et al. 2017): one at low energies peaking between the submillimeter and EUV, sometimes extending into the soft X-ray band, and explained as synchrotron radiation of a population of accelerating jet electrons, and a second maximum in the hard X-ray/γ-ray regime, usually explained as inverse Compton (IC) radiation from photons that scatter off the jet electrons.These photons are located either inside the jet (synchrotron-self-Compton radiation; SSC) or they are emitted by an external medium such as the broad-line region (BLR) or torus (external Comptonization; EC).Additionally, or alternatively, hadronic processes (ultra-relativistic protons) may contribute at high energies (e.g., Böttcher 2019).
OJ 287 is a key representative of the class of blazars.It is relatively nearby (redshift z = 0.306), very bright across the electromagnetic spectrum, highly variable, and highly polarized (e.g., Goddi et al. 2021;Komossa et al. 2021c;Valtonen et al. 2021;Abdollahi et al. 2022;Komossa et al. 2022a, and references therein).It exhibits exceptional optical flares that reached as bright as 12th magnitude (Kinman et al. 1971).OJ 287 has been detected in the γ-ray and very-high energy (VHE) regime (Abdo et al. 2009;O'Brien 2017).Due to a very dense and dedicated multiwavelength Swift (X-ray-UVoptical) monitoring campaign that we set up > 6 yrs ago in the course of the project MOMO (Multiwavelength Observations and Modelling of OJ 287; Komossa et al. 2021d), OJ 287 has one of the best-covered long-term Xray-UV-optical light curves and simultaneous SEDs of any blazar (Komossa et al. 2021c), along with the dense multifrequency radio monitoring presented here.
On the basis of the faintness of its optical emission lines from the broad-line region, only occasionally detected (Sitko & Junkkarinen 1985;Nilsson et al. 2010), OJ 287 is classified as a BL Lac object.OJ 287 caught particular attention following its bright and long-lasting optical outbursts in the 1970s and 80s.After strong evidence was presented that these outbursts repeat every ∼11-12 yr (Sillanpää et al. 1988) and are double-peaked with a peak separation of order 1 yr (Sillanpää et al. 1996), several binary black hole scenarios were suggested (see Komossa et al. 2021b, for a recent review), involving: (1) a mildly precessing secondary SMBH that is nearly in the accretion-disk plane and is tidally perturbing the primary SMBH's accretion disk near periastron (Sillanpää et al. 1988); or (2) a highly inclined, highly precessing secondary SMBH impacting the primary's accretion disk twice during each orbit (precessing binary (PB) model hereafter; Lehto & Valtonen 1996;Valtonen et al. 2021); or (3) a precessing accretion disk (Katz 1997); or (4) Doppler boosting of two jetted SMBHs orbiting each other (Villata et al. 1998); or (5) a secondary SMBH that impacts the primary's accretion disk only once during its orbit with subsequent tidal perturbations feeding the inner jet (Valtaoja et al. 2000).While the PB model required a very massive primary SMBH of order 10 10 M , Liu & Wu (2002) explored a variant of the Valtaoja et al. (2000) model that only required a primary's mass of order 10 8 M .
The pattern of variability of OJ 287 is not strictly periodic, but comes with deviations from strict periodicity of order 1 yr (Valtaoja et al. 2000).One of the largest deviations happened in [2005][2006] when the brightest outburst occurred in 2005 October-November (Villforth et al. 2010) while expected only in 2006.
The best explored binary scenario is the PB model, based on detailed orbital modelling of the system.Parameters of that model were adjusted each time a supposed 'main outburst' was observed at a different epoch than expected, then a new prediction was made, and so on.The model claims very high accuracy in the predicted timing of subsequent flares on the order of days -hours, for events that are separated by more than a decade.
Non-binary models producing semi-periodicity in AGN light curves in general have also been explored (e.g.Rieger 2004;Liska et al. 2018), some directly applied to OJ 287: Britzen et al. (2018) favored disk precession, leaving open the cause for precession (that could still be a binary), and Villforth et al. (2010) presented cautious comments on a binary interpretation and speculated about a model of magnetic breathing to explain the optical double-peaked outbursts of OJ 287.
Therefore, in addition to understanding disk-jet physics of this bright blazar, the project MOMO was carefully designed for rigorously distinguishing between different binary (and non-binary) SMBH scenarios of OJ 287 in general, and testing predictions of the PB model in particular, including distinctly different predictions made in the radio regime (Section 2).While most previous results focussed on the UV-optical, X-ray and γ-ray regime (e.g., Komossa et al. 2020Komossa et al. , 2021aKomossa et al. , 2022a), here we concentrate on the complete 2015-2022 multifrequency radio observations and their interpretation and implications.
This paper is organized as follows: In Section 2 we describe the key motivation and goals of the MOMO project with emphasis on the questions that will be addressed with the radio multifrequency observations presented here.Sections 3-6 provide a description of the data acquisition and analysis in the radio, optical, UV, X-ray, and γ-ray bands.In Section 7, the radio flux variability and spectral properties are established and discussed in detail, including DCF analyses, followed by a comparison with the multiwavelength light curves in Section 8. Results are discussed in Section 9 with particular focus on timing results, the extraordinary 2016/2017 outburst, and constraints on binary SMBH scenarios.A summary and conclusions are provided in Section 10.
Timescales and frequencies are given in the observer's frame when reporting measurement results, unless noted otherwise.We use a cosmology with H 0 =70 km s −1 Mpc −1 , Ω M =0.3 and Ω Λ =0.7 throughout this paper.At the distance of OJ 287, this corresponds to a scale of 4.5 kpc/arcsec (Wright 2006).

Overall introduction, goals, and previous results
In the course of the project MOMO (Komossa et al. 2017;Myserlis et al. 2018;Komossa et al. 2020Komossa et al. , 2021aKomossa et al. ,c,d, 2022a,b),b), we are carrying out dedicated, dense, long-term monitoring of OJ 287 at >16 frequencies from radio to X-rays.This is the densest monitoring project of OJ 287 so far carried out involving X-rays.MOMO provides spectral and timing information, and simultaneous SEDs, at all activity states of OJ 287.In particular, we are using the Effelsberg 100 m radio telescope and the Neil Gehrels Swift observatory (Swift hereafter; Gehrels et al. 2004) since 2015 December.Essentially all Swift observations of OJ 287 in recent years were obtained by us in the course of the MOMO program.The rich data sets and their interpretation are presented in a sequence of publications and so far include: (1) Our detection of two major nonthermal X-ray-UV-optical outbursts with Swift in 2016/17, and 2020 (Komossa et al. 2017(Komossa et al. , 2020)).(2) The detection of variable radio polarization in 2016 (Myserlis et al. 2018).(3) The rapid follow-up of the 2020 outburst with Swift, XMM-Newton and NuSTAR establishing the spectral components up to ∼ 70 keV, including an exceptionally strong soft X-ray excess of synchrotron origin and an unexpectedly steep spectrum in the NuSTAR band (Komossa et al. 2020).(4) XMM-Newton and Swift spectroscopy during the 2018 Event Horizon Telescope (EHT) campaign (Komossa et al. 2021a).( 5) A spectral analysis of all XMM-Newton observations during the last two decades (Komossa et al. 2021a), and a spectral and detailed timing analysis of all Swift observations during the last two decades (Komossa et al. 2021c), establishing OJ 287 as one of the most X-ray spectrally variable blazars known, and with strong implications for emission mechanisms.(6) Effelsberg observations between 2019 and 2021 in a period of high activity (Komossa et al. 2022a), and implying that the radio emission is dominated by the main jet.(7) First results on the absence of precursor flare activity (Komossa et al. 2022b), falsifying predictions of one of the binary SMBH model variants.
Data we obtain with Swift, the Effelsberg telescope, and at other facilities are analyzed quickly, typically within days, and the community is then rapidly alerted about exceptional flux states (outbursts, deep fades) and spectral states of OJ 287 in a sequence of Astronomer's Telegrams since 2015 (ATel #8411, #9629, #10043, #12086, #13658, #13702, #13785, #14052, #15145, and #15764).This allows independent data to be taken by the community at frequencies or with instruments not covered by our own program.For instance, the bright 2016/2017 X-ray outburst we detected with Swift (e.g., ATel #9629, #10043; Komossa et al. 2017Komossa et al. , 2020) ) was used to trigger a VERITAS very-high energy (VHE) observation (O'Brien 2017), and the bright 2020 X-ray outburst we detected with Swift (ATel #13658) triggered an AstroSAT observation (Prince et al. 2021;Singh et al. 2022) in addition to our own follow-ups with XMM-Newton and NuSTAR (Komossa et al. 2020).

MOMO-radio
Observations in the radio regime allow us to distinguish between different binary SMBH scenarios in general, and to test predictions of the well-developed PB model in particular, based on distinctly different predictions they make in the radio regime for the timing and properties of the first and second optical peak of the double-peaked outbursts and for any possible afterflares (Valtaoja et al. 2000).For instance, if major optical outbursts are thermal in nature, they will not be accompanied by radio (synchrotron) outbursts 1 .If both optical peaks are due to synchrotron emission, we expect two radio flares, with polarization evolution following synchrotron theory.These tests require long-term coverage over at least a significant fraction of one binary orbital period.The long-term goal of MOMO is to cover at least two orbital periods.
Further, independent of the binary hypothesis and of blazar physics directly traced by the Effelsberg observations themselves, we note that the flux and spectralindex measurements from the dense multifrequency radio coverage allow extrapolation to other radio frequencies of interest in deep radio imaging studies of OJ 287 that continue to be carried out with VLBI techniques (e.g., Lister et al. 2021;Weaver et al. 2022;Zhao et al. 2022).
Prior to MOMO, OJ 287 has been the target of other single-dish radio monitoring campaigns, for instance at  99-15, 19-16, 12-17, 13-18, 75-19, 65-20, 70-21, and 83-22).Note that some receivers have changed over the years, and at some epochs a larger number of frequencies was observed, partly depending on weather conditions.νc is the central frequency at which the flux density was extracted, ∆ν the band width, and HPBW the half power beam width.Since the employed central frequency of some bands changes with time, we assign a fixed representative frequency ν * (as listed in column 5) in the text and in labels of figures.For any analysis and plots, the exact frequencies are used.UMRAO (Hughes et al. 1998), Metsähovi (Hovatta et al. 2008), and Effelsberg (Fuhrmann et al. 2016).These did not include simultaneous UV and X-ray data.

EFFELSBERG DATA ACQUISITION AND REDUCTION
We started monitoring OJ 287 in the radio band with the Effelsberg 100 m telescope (Kraus et al. 2003) in December 2015 in the course of the MOMO program.
Here we discuss all the data obtained until June 2022.
OJ 287 is unobservable each year for about 3 months with ground-based optical telescopes, or space-based missions such as Swift, because of its close proximity to the Sun.Our radio observations have the advantage that OJ 287 can be observed over nearly the whole year.The Sun constraint of the Effelsberg telescope is only 2 degrees, so that sources at larger than 2 degrees separation from the Sun are observable.
The average cadence of our radio observations is 14 days (at 10.45 GHz).The cadence is higher at selected epochs, and our proposals contained a ToO mode to obtain denser coverage, should some exceptional variability be observed.We sometimes experience larger gaps, depending on times of telescope maintenance, availability of receivers, and weather conditions; the latter limiting the cadence for the observations at the highest frequencies (above ∼20 GHz).
Data from the first year of observations we already presented in Myserlis et al. (2018), and data between 2019 to 2021 May were discussed by Komossa et al. (2022a).These data are re-added here for completeness of discussion, and new analyses are performed on them.
Observations were carried out between 2.6 and 44 GHz, and occasionally at 1.4 GHz, switching between up to eight receivers.Some of the receivers changed in the course of the monitoring program (see Tab. 1 for details). 2Since the central frequency we employed varied in a few cases over the length of the observing program, we assigned a fixed representative frequency ν * throughout the text and in some figure labels (for instance, 2.595 and 2.64 GHz are referred to in the text as 2.6 GHz, etc).The exact frequency at any given date can be found in Tab. 1.
The cross-scan method (Kraus et al. 2003) was used to acquire the radio data.In the cross-scans, the telescope was moved in two perpendicular directions, azimuth and elevation, with the target source position at the center of the scans.For the low frequencies, 2-3 sub-scans were carried out per direction.For the high frequencies, it was 12-16 sub-scans per direction.
The data reduction and analysis was performed in a standard manner as described in, e.g., Kraus et al. (2003).A detailed description of the analysis procedures will be given by Kraus et al. (in preparation).In a first analysis step, the data of every single sub-scan were fit with a Gaussian profile.Bad sub-scans were identified and excluded from further analysis.Such bad data sets can be due to high pointing errors, radio frequency interference (RFI), or -in case of OJ 287 -disturbances by solar radiation around August 1st of each year due to the source's close proximity to the sun, for instance.
After correcting for small pointing errors of the telescope, the amplitudes of the individual sub-scans were averaged.In some cases, especially at the highest fre-quencies and in mediocre weather conditions, the subscans of one direction (azimuth/elevation) were averaged before the Gaussian profile was fit, in order to increase the signal to noise ratio.Next, corrections for the atmosphere's opacity were applied as well as for the gainelevation effect (the change of the sensitivity with elevation).Absolute flux calibration was then achieved by comparing the observed antenna temperatures with the expected flux densities of selected calibrators such as 3C 286 (Ott et al. 1994) 3 .The complete 2015-2022 Effelsberg radio light curves are shown in Fig. 1.

SUBMILLIMETER ARRAY (SMA) DATA ACQUISITION AND REDUCTION
Flux density data at 1.3 mm, and occasionally at 1.1 mm and 870 µm, were obtained at the Submillimeter Array (SMA) near the summit of Maunakea (Hawaii).OJ 287 is included in an ongoing monitoring program at the SMA to determine the fluxes of compact extragalactic radio sources that can be used as calibrators at mm wavelengths (Gurwell et al. 2007).Observations of available potential calibrators are from time to time observed for three to five minutes, with the measured source signal strength calibrated against known standards, typically solar system objects (Titan, Uranus, Neptune, or Callisto).Data from this program are updated regularly and are available at the SMA website 4 .Here, we show the data from 2015 October to 2022 June (Fig. 2).The majority of the data in the 1.3 mm band are taken at mean frequencies near 225 GHz, but the exact frequency can vary between 213 GHz and 240 GHz, depending on the science goal of the initial observation.Frequencies   is the exposure time of each single XRT or UVOT observation.Each UVOT band V:B:U:W1:M2:W2 is observed with a ratio of 1:1:1:2:3:4 of the total exposure time, respectively.
in the 1.1mm band of the data reported here range between 270-272 GHz, and in the 870µm band between 339-351 GHz.

SWIFT DATA ACQUISITION AND REDUCTION
The majority of Swift data we obtained in the course of the MOMO-UO (UV and optical) and MOMO-X (X-rays) program was presented in previous work since 2017, including a detailed discussion of the full 2005-2021 Swift light curve of OJ 287 by Komossa et al. (2021c).Here, we newly present the most recent Swift observations since 2022 mid-February (Tab.2).
The Swift UV-optical telescope (UVOT; Roming et al. 2005) and X-ray telescope (XRT; Burrows et al. 2005) data reduction is carried out in a standard manner and is described in great detail by Komossa et al. (2021c).It includes target source and background selection, instrumental corrections, and correction for Galactic extinction of the UVOT fluxes.In the X-ray band, absorptioncorrected fluxes are derived from single power-law spectral fits with fixed Galactic absorption.

FERMI γ-RAY DATA
In the γ-ray band, Fermi LAT (Large Area Telescope; Atwood et al. 2009) data are used.Unlike the Effelsberg and Swift observations that are proposed and analyzed by us, the γ-ray data were taken from the Fermi archive.They cover the energy regime 0.1-100 GeV.Publicly available Fermi LAT data of OJ 287 were retrieved from the Fermi LAT light-curve repository (Kocevski et al. 2021) 5 .We used weekly averages of the fluxes derived from the spectral model fit of a logarithmic parabolic power-law with fixed parameters (including a fixed photon index of 2.16) and free normalization, as described by Kocevski et al. (2021).This approach is preferred in our case over fits with free index (e.g.Hodgson et al. 2017;Kapanadze et al. 2018), because OJ 287 remains in a low γ-ray state most of the time during the epochs of interest..8Jy at 32 GHz.The radio outburst coincides with the bright and long lasting X-ray-UV-optical outburst of OJ 287, that was detected in the course of our MOMO Swift monitoring (Komossa et al. 2017(Komossa et al. , 2020)).During this time, we see the highestamplitude of radio variability in recent years.The flux density rises from a deep low-state in the radio emission in 2016 January by a factor of 3.2 (at 32 GHz) to the observed peak in February 2017.A similar rise is seen in the SMA band with an amplitude of a factor of 2.3 between January and the highest measured value in February.Due to the lower cadence with no coverage of the flare beyond 2017 February 16, the true peak of the flare cannot be measured at 225 GHz.
Several more radio flares are detected.The secondbrightest is observed in 2021 November to 2022 May (see also Komossa et al. 2022b).The 36 GHz flux density rises by a factor 2.2 from minimum to peak.As the flare evolves, the 36 and 44 GHz light curves show a fast decline and renewed rise to a second maximum in 2022 February.Sharp variability at that epoch is also seen in the SMA data.The radio flux starts to decline in 2022 May.
Every ∼2 years, deep radio minima are observed.The lowest states are found in late 2015, late 2017 (coincident with the optical-UV deep fade; Komossa et al. 2021c), late 2019, and mid 2021.
The most rapid flaring is seen at high frequencies.Flares are often smoothed into broader structures at low frequencies due to opacity effects; see, for instance, the rapid dipping and recovery at high frequencies during the 2022 flaring that is not observed at lower frequencies.In order to quantify the variability in each frequency band, a fractional variability amplitude analysis was carried out, which we discuss next.

Fractional variability amplitude
The fractional rms variability amplitude F var (Edelson et al. 2002;Vaughan et al. 2003) of the radio flux density was calculated separately for each radio frequency at annual time bins (Tab.3) and as a long-term average.
F var is given by: where S 2 is the variance of the light curve, σ 2 err is the mean square of the measurement errors, and x is the mean flux (Vaughan et al. 2003).
The error of F var was calculated according to Edelson et al. (2002), as: where N is the number of data points used in the computation of F var .Results are shown in Tab. 3 and Fig. 3.In the longterm average, F var increases from low to high frequency, between 0.17±0.01(at 2.6 GHz) and 0.28±0.03(at 43 GHz).Overall, F var is highest during the large outburst in 2016-2017.Results are further discussed below.

Intra-day variability (IDV)
Two epochs were covered with daily cadence at the Effelsberg telescope: One in December 2015, and one in April 2017 (when for technical reasons only sources at certain sky locations could be observed, and OJ 287 was covered repeatedly as it was outbursting at that time).Two more epochs were observed at high cadence with the SMA in 2022 January and March.These observations allow us to search for variability within ∼1 d.
Daily coverage at Effelsberg was obtained between 2017 April 6 and 16, with flux densities at 32 GHz ranging between 6.79 ± 0.04 Jy and 7.77 ± 0.05 Jy within 7 days.The most rapid change within a day was from 7.77 ± 0.05 Jy to 7.45 ± 0.08 Jy, on the order of 4% decline.
During the near-daily cadence from 2022 January 11 to 17, SMA radio flux densities at 225.5 GHz ranged between 3.63 ± 0.18 Jy and 4.03 ± 0.20 Jy with no evidence for short-term variability within the errors.No short-term variability within the errors is detected with SMA between 2022 March 2 and 8, either.
The lack of significant variability on the timescale of ∼ 1d beyond ∼4% is of interest for long-duration observations of OJ 287 on that timescale, for instance with RadioASTRON or EHT.

Radio SEDs
The radio SEDs of OJ 287 are highly variable.Turnover frequencies range between 10 and 36 GHz, or are outside the observed frequency band at >43 GHz at other epochs.The bright 2016/2017 outburst shows a strong spectral evolution.
The detailed shapes of all SEDs are displayed in Fig. 4 (with the exception of SED sequences between 2019 and 2022 that were already presented in Komossa et al. (2022a) and are not repeated here.Those cover the epochs MJD 58511 -58784, MJD 58799 -59028, and MJD 59034 -59426).The time period 2016 was already in our first paper on radio results obtained within the MOMO project (Myserlis et al. 2018) and is shown here again to visualize the spectral evolution of the outburst that continued into 2017.

Time delays and DCF
Opacity effects due to synchrotron self-absorption lead to delays in the rise times of flares (e.g., Lee et al. 2020).In order to quantify this effect across the long-term light curve at all frequencies between 2015 December and 2022 June we have computed the discrete correlation function (DCF).The DCF technique was designed to analyze unevenly sampled data sets (Edelson & Krolik 1988).We employ the DCF to search for time lags between the radio bands, measured against the 10.45 GHz data set.
We computed the DCFs as prescribed in Edelson & Krolik (1988), selecting the time step, τ , over which the DCFs were computed, to be twice the median time step across the entire light curve at each radio frequency.
To evaluate the significance level of measured lags we produced confidence contours for each DCF by simulating N = 10 3 artificial light curves in each band, following the prescription of Timmer & König (1995), assuming a power spectral density (PSD) of P (f ) ∝ f −α = f −1.97 .This value of α is based on the results of a struc- The artificial light curves were then used to compute artificial DCFs with the band-of-interest light curve, allowing for the computation of 90 th , 95 th , and 99 th percentiles based on the distribution of artificial DCFs at each time step.
To evaluate the error on the measured lags, we applied a least-squares minimization of a Gaussian function to the central peak, defined as the points with correlation strength corresponding to ≥ 50% of the peak value, within τ = ±400 days, which enabled measurement of the DCF centroid.To evaluate the 68 per cent confidence region of the lag (centroid) measurement, we employed the resampling technique outlined by Peterson (1998).Briefly, for each point in the two light curves used in a given DCF the flux density is resampled according to a Gaussian distribution with mean and standard deviation equal to the flux density and its error for that point.Then, a random subset of points from each light curve is selected and used to compute a new 'resampled' DCF, from which the centroid is measured in the aforementioned way.This procedure is conducted N = 10 3 times, allowing the 68 per cent confidence interval on the 'true' DCF centroid to be evaluated using the 16 th and 84 th percentiles of the 'resampled' DCF centroid distribution.Results are shown in Fig. 5 and 6, and Tab. 5.The large 2016/2017 X-ray-UV-optical outburst (Komossa et al. 2017, 2020) is well correlated with the radio emission that reaches its highest flux density at high frequencies during the whole MOMO monitoring period (Fig. 10).The X-ray-optical outburst shows  ,black crosses: SMA data between 213 and 351 GHz).The γ-ray flux (observed 0.1-100 GeV band; one-week averages), the absorption-corrected X-ray flux (observed 0.3-10 keV band), and the extinction-corrected UV flux at λ obs =1928 Å are given in units of 10 −11 erg s −1 cm −2 .The radio data are reported as flux densities in Jy.The gap in Swift observations between June to September each year is due to the close proximity of OJ 287 to the Sun, such that it becomes unobservable with Swift.Fermi and Effelsberg observations are not affected by this Sun constraint.Note that the single brightest γ-ray data point from 2015 December is off the scale (see Fig. 10 instead).some pronounced substructure with several peaks seen at all wavebands, and the optical-UV bands reached their highest peak flux earlier than the X-ray band.The brightest optical-UV peak was reached on 2016 October 20, whereas the brightest X-ray peak was only reached on 2017 February 2. Our monitoring cadence at that time was 1 day.The brightest radio peak was observed close in time to the brightest X-ray peak, on 2017 February 28 (MJD 57812 +4 −12 days, where the uncertainty reflects the cadence of the radio monitoring at that epoch).A second, earlier radio peak detected at all radio frequencies occurred on 2016 July 31 (MJD 57600 +34 −7 days), at an epoch when OJ 287 was in Swift Sun constraint, and no data in the optical-UV and Xrays were taken.
Fig. 8 compares the spectral indices α ν measured in selected radio bands with the optical-UV spectral index.The optical-UV spectral index, measured between V (5468 Å) and UVW2 (1928 Å), shows a 'bluer-whenbrightest' behavior, and is steepest during several lowstates.Optical and radio spectral indices are correlated.9. DISCUSSION

Radio variability
The radio light curves of OJ 287 show high-amplitude variability with multiple flares and several deep lowstates.On average, the data show an increase in F var from low-to high frequencies (Fig. 3), as expected in a plasma that becomes increasingly optically thin towards higher radio frequencies.An exception is seen around 19 GHz, where F var shows a dip.However, the bulk of the effect is caused by the lower number of 19 GHz observations at the beginning of the monitoring interval in 2016, when a strong rise into high-state occurred, and F var was highest.
Fractional variability amplitudes in the radio band are found to be lower than values obtained for the UVoptical and X-ray light curves (F var,V = 0.3 , F var,X = 0.4 respectively) at similar time intervals (Komossa et al. 2021c), reflecting decreasing opacity effects in the optical regime.The trend across frequency and the values of F var of OJ 287 are similar to those observed in other blazars (e.g., Schleicher et al. 2019).Values of F var are highest during the 2016/2017 outburst.This outburst will be further discussed below.
The majority of the radio and optical flares are closely correlated, with the optical leading the radio (see also Komossa et al. 2022a).This hints at a co-spatial origin of optical and radio emission.An exception is the latest 2021-2022 radio flare.Its bright state in late 2021 is not accompanied by major optical flaring.However, it is associated with two sharp γ-ray flares; the two bright-est in recent years (Fig. 11).One at the end of 2021 October and a second, broader one that starts in late December.This type of behaviour indicates a possible causal connection.It is conceivable that the first sharp γ-ray flare occurs when the jet is hitting the first recollimation shock.The second γ-ray flare is broader/slower, indicating a larger size of the emission region.It is possibly produced further downstream where the jet also becomes broader.We also note that the 2016/2017 outburst is accompanied in its rise phase by a γ-ray flare in 2016 August.

Deep fades every 1-2 years
A remarkable sharp, symmetric, UV-optical deep fade was detected with Swift during 2017 October-December (Komossa et al. 2021c) that we now also cover in the radio band where it is seen as well at all frequencies.Similar, but less pronounced, UV-optical low-states are also detected in our light curves in late 2019, late 2020, and late 2021.Several of these have counterparts in the radio band in the form of deep dips of the radio emission at all bands, seen in late 2015, late 2017, and late 2019, whereas the radio deep fade in 2021 comes a few months earlier, in mid 2021.Another such radio dip is seen in the radio light curve of Lee et al. (2020) in late 2013.
With only a few well-documented occurrences in our data, a rigorous period search is not yet warranted.However, we note that the autocorrelation function (ACF) at 10.45 GHz (Fig. 9) does show evidence for periodicity at 410 days and at 650 days.It is also interesting to note that several previous studies of OJ 287 reported evidence for periodicities in the range of 1-1.7 yrs in the radio (Hughes et al. 1998;Hovatta et al. 2008;Britzen et al. 2018), and ∼400 days in the optical (Bhatta et al. 2016;Sandrinelli et al. 2016).Independent of the possibility of semi-periodicity, the question is raised what causes these deep fades.They may just represent the quiescent states inbetween epochs of flaring.However, their symmetric nature of fading and rising without remaining in quiescent states for significant amounts of time suggests otherwise.Deep fades could be caused by extreme absorption events, de-beaming due to de-collimation or wobbling of the jet, or they could be directly related to the electron launching process at accretion-disk scales in the first place.
The deep fades do not represent SAV (Symmetric Achromatic Variability; Vedantham et al. 2017) events that were speculated to be caused by milli-lensing from intermediate-mass black holes, since they are not achromatic but show a systematic change in radio spectral indices (Fig. 8).Some deep fades show a U-shaped structure that is particularly similar to the SAVs (e.g., the 2017 UV-optical deep fade or the 2021 radio deep fade), but since they go deeper than the surrounding quiescent level and are not achromatic, the observations of OJ 287 also show, that there are mechanisms other than possible milli-lensing to produce these characteristic structures in blazar lightcurves.
Independent of the interpretation of these events, the prediction of their future occurrences is of interest for deep imaging studies of the host galaxy of OJ 287 (as done shortly after the 2017 optical deep fade by Nilsson et al. 2020) with JWST, the Hubble Space Telescope, or ESO's Enhanced Resolution Imaging Spectrograph (ERIS), and for high-resolution optical spectroscopy in search of host absorption features for direct SMBH mass measurements from stellar velocity dispersion σ * .At epochs of deep fades the blazar glare itself is least interfering with the measurements.9.3.Absence of predicted 2021 December precursor flare activity It was speculated that one of the optical flares in the 2005 light curve of OJ 287 was driven by binary SMBH activity in the form of a 'precursor flare' preceding the main outburst.Based on the PB model it was predicted that the precursor flare would repeat in 2021, on December 23, with an optical thermal bremsstrahlung spectrum of index α ν = −0.2(and without any X-ray or radio counterpart; Valtonen et al. (2021), their Section 3).However, OJ 287 was found in a deep optical-UV low-state throughout December 2021, and any meaningful optical-UV flaring activity at that time period could be excluded (Komossa et al. 2022b).Further, neither the rise in emission after the deep low-state, nor any other UV-optical activity until June 2022, showed the predicted thermal bremsstrahlung spectrum.Observed spectral indices α ν,opt−UV range between -1.2 and -1.54 (Fig. 8) whereas the predicted thermal bremsstrahlung spectrum would have had α ν ∼ −0.2.We therefore conclude that the bremsstrahlung flare predicted by the PB model did not happen, neither in 2021 December nor at any other epoch.
We now turn our attention to the remarkable 2016/2017 outburst and its interpretation in the context of previous variants of binary models that do not involve a strong orbital precession, in contrast to the PB model.9.4.The 2016/2017 outburst, its interpretation as the latest double-peaked outburst of OJ 287, and implications for binary SMBH models 9.4.1.MWL properties of the 2016/2017 outburst A high-amplitude outburst occurred in 2016/2017, detected in the course of our Swift monitoring in all bands from the X-rays to the optical with a systematic X-ray 'softer-when-brighter' variability pattern with steep Xray spectral indices Γ x 3 at peak brightness (Komossa et al. 2017(Komossa et al. , 2020)).While Komossa et al. (2017) first speculated about a thermal nature of this outburst, detailed follow-ups then established the nonthermal, Synchrotron, nature of this event, based on the following arguments: First, with Swift, we detected X-ray flux doubling timescales as short as 4 days, shorter than the light-crossing time at the last stable orbit of the accretion disk around a primary SMBH if its mass was M BH = 1.8 × 10 10 M as required in the PB model, ruling out a primary's disk origin (Komossa et al. 2021a) 6 .Second, optical-UV DCF results (Komossa et al. 2021c) are consistent with synchrotron theory, but the lags are too small for accretion-disk reverberation (Kammoun et al. 2021) of a SMBH with a mass as low as ∼ 10 8 M .Third, the Swift X-ray spectra are well explained by a soft synchrotron emission component and show the same softer-when-brighter variability pattern also seen during other outbursts (Komossa et al. 2020(Komossa et al. , 2022a;;Idesawa et al. 1997).Fourth, the early phase of the outburst was already associated with an increase in radio emission (Myserlis et al. 2018;Kapanadze et al. 2018;Lee et al. 2020), the levels of radio polarization were high (Myserlis et al. 2018;Goddi et al. 2021), and VHE emission near the X-ray peak was detected by VERITAS (O'Brien 2017) [only the Fermi γ-ray regime lacked significant flaring (Kapanadze et al. 2018;Komossa et al. 2022a)].
Here, we provide the coverage of the whole outburst at multiple radio frequencies between 2.6 and 225 GHz, strongly re-confirming the nonthermal nature of the outburst.The radio emission is closely correlated with the MWL emission, and the brightest radio and X-ray peak coincide closely and are both reached in 2017 February (Fig. 10).The radio peak flux in 2017 is among the highest so far recorded from OJ 287.It is comparable to the brightest state in the 1972-1996 light curve reported by Valtaoja et al. (2000).
In the next Section, we will argue that the extraordinary 2016/2017 outburst matches the properties of the previous famous semi-periodic double-peaked outbursts of OJ 287 of the 1970s -2000s and that it is the latest such event in the sequence.First assuming strict periodicity of the previous double peaks with ∆t = 11.86 yr and using the timing of the well-observed 1994 burst (Valtaoja et al. 2000), the latest one should have occurred in 2018.However, no burst was observed at that epoch, but in 2016/2017 instead.Therefore, next, we drop the requirement of strict periodicity as also done in most previous works and for the additional reasons given below, and we then discuss implications for binary SMBH models of OJ 287.

Re-interpretation of the 2016/2017 outburst
The fact that the precursor flare predicted by the PB model was not observed, and the fact that the most pronounced outburst in recent years was the one in 2016/2017 (not expected in the PB model), leads us to re-consider alternative binary SMBH models of OJ 287.We ask if the timing of the 2016/2017 outburst fits any of the previous models, that involve less extreme orbital precession and that lead to the expectation of more regularly-spaced outbursts than the PB model, with a known variance on the order of 1-2 years as directly ob-served previously during the major double-peaked outbursts.Most authors reported time intervals of ∆t=11-12 yrs between the outbursts (e.g., Sillanpää et al. 1988;Valtaoja et al. 2000) directly from observations, with additional uncertainties in ∆t of at least several months based on (1) gaps in the cadence of observations during the epochs of dedicated monitoring since 1970, (2) additional gaps of ∼3 months each year when OJ 287 was unobservable due to its solar proximity, (3) the increasing availability of photographic plate archives that cover selected time intervals since the 1880s (so that the timing of single peaks could shift back or forth with the availability of new data sets), ( 4) each authors' preference for the reliability of one or another time period of photographic plate data, and/or (5) different methods of single peak determinations based on the same data sets over time.
If we take the timing of the first peak in the 1970s when the dedicated monitoring started, t peak =1971.13(Valtonen et al. 2021), and the peak of the first 2016 V-band flare at 2016 February 29 (Fig. 8; upper panel), then this results in an average ∆t of 11.3 yr for 4 cycles.
How do the properties of the 2016/2017 outburst compare with those of previous double-peaked outbursts of OJ 287 ?First, the 2016/2017 outburst shows a second, bright optical peak on 2016 October 20 (Fig. 10), 8 months after the first peak, matching the wellestablished characteristic historical double-peak structure of past flares, with the second peak appearing ≈1 yr after the first peak (e.g., Valtaoja et al. 2000).Further, the 2016/2017 outburst agrees with the observation of Valtaoja et al. (2000) in that the second peak is associated with radio flaring.However, we also find evidence that the first peak has a radio counterpart at high frequencies.It is possible that it was missed in earlier observations because it was self-absorbed or due to gaps in the coverage.The whole 2016/2017 outburst is associated with radio emission.Another similarity of the 2016/2017 outburst with previous double-peaked outbursts is the fact that no major X-ray flare is associated with the first peak (Idesawa et al. 1997).Apart from the 2016/2017 outburst, no other such bright and longlived outburst has been observed since the previous one in late 2005-2007.Our results are also consistent with other evidence based on polarimetry iin 2016 that the outburst originated in the central region, at a projected size of ≤0.67 pc Myserlis et al. (2018).
In the 2010s, for the first time we achieved a very dense light curve coverage beyond the optical band, including several UV bands and X-rays up to 10 keV (in addition to multiple optical bands, and multiple radio frequencies).Therefore, predictions of a recurrence of a similar outburst in a time frame of 10-12 years from now will be testable in many wavebands.Given the complex gas physics and magnetic fields involved, the best way to identify binaries will continue to be a phenomenological one based on semi-periodicities [or direct spatial resolution of the components, if both of them are radio emitting (Rodriguez et al. 2006;Komossa & Zensus 2016)].Possible emission mechanisms within the OJ 287 binary scenarios to explain the 2016/2017 outburst will be discussed further in future work.
In summary, we followed a phenomenological approach in identifying the latest double-peaked outburst, since high-precision modelling is still challenging given the huge numbers of free parameters and the complicated physics of this problem of a gas-rich binary system with disk impacts, and where strict periodicity is not expected even if the secondary is not highly precessing.Additional motivation for our approach came from the fact the high-precision prediction of the PB model of characteristic flaring activity in 2021 was not confirmed by our observations.The identification of the 2016-2017 outburst as the latest optical double-peaked outburst of OJ 287 is consistent with alternative binary models, matching several aspects of the models of, e.g., Valtaoja et al. (2000) and Liu & Wu (2002).A revision of those models would still be needed to account for the fact, that the latest bursts came early, and for the fact that the whole double-peaked multi-months burst structure is associated with radio emission.Such a revised model is beyond the scope of this publication, and will be discussed in future work.
Our new prediction is, that we should observe a new double-peaked outburst with properties matching the characteristic X-ray-UV-optical and radio properties of the 2016-2017 outburst, and with the first peak in 2026-2028.In contrast, the PB model predicted the latest outburst to happen in October 2022 (Valtonen et al. 2022).At that epoch we found OJ 287 in a low-state instead (Komossa et al. 2023).
Finally, it is interesting to note that the brightness of past double-peaked outbursts has systematically declined over the last few decades (Valtaoja et al. 2000;Dey et al. 2018;Komossa et al. 2021c), no longer reaching the exceptional magnitudes of the 1970s and 1980s, and therefore making it more difficult to distinguish between blazar flaring in general, and binary-driven flaring in particular, based on peak flux alone.However, including pre-1970 data in the long-term light curve, a possible period of ∼60 years in the brightness of double peaks has been identified (Valtonen et al. 2006).If it continues, we then expect that future double peaks will become systematically brighter again.

SUMMARY AND CONCLUSIONS
We have presented densely covered multifrequency radio light curves of OJ 287 between 2015 and 2022, along with optical, UV, X-ray and γ-ray observations.All data except in the γ-ray band were obtained by us in the course of the project MOMO.OJ 287 displayed a broad range of extraordinary activity including multiple flare events, deep fades, and strong spectral variability.Our main results can be summarized as follows: • We have characterized in detail the radio flux and spectral variability of OJ287 between 2015 and 2022, including turn-over frequencies, spectral indices, fractional variability amplitudes and DCFs.
In particular, we densely cover the large nonthermal 2016/2017 MWL outburst that is accompanied by strong radio flaring.
• We find that deep low-states repeat every 1-2 years in the optical and radio regime.Independent of their nature, predicting the next such events is important, as optical deep fades will facilitate host galaxy imaging spectroscopy for independent SMBH mass measurements, e.g., via stellar velocity dispersion.
• The two brightest γ-ray flares in recent years coincide with the sharp rise and re-rise of the bright 2021-2022 radio flare, suggesting a causal connection.
• The light curve evolution and spectroscopic information have been used to test predictions of binary SMBH models.Precursor flare activity, predicted by the precessing binary model to occur in 2021 December, is absent.Neither the flare, nor the thermal bremsstrahlung spectrum were observed; neither in 2021 December nor any other time until 2022 June.
• We interpret the big 2016/2017 outburst as the latest of the characteristic semi-periodic doublepeaked outbursts that OJ 287 is famous for.This favors binary SMBH models of OJ 287 that involve a non-(or mildly)precessing binary.A scenario involving a highly precessing binary (PB model) with a particularly high primary SMBH mass of order 10 10 M is then no longer required.
This interpretation leads to the prediction of the next double-peaked outburst in the period 2026-2028.Instead, the PB model predicted this outburst to happen in 2022 October.However, no outburst was observed at that time (Komossa et al. 2023).The MOMO project continues with the ultimate goal of covering one to two decades of dense multiwavelength observations of OJ 287.
It is our great pleasure to thank the Swift team for carrying out the observations of OJ 287 that we proposed, and for very useful discussions regarding the observational set-up and related questions throughout the years.This work is partly based on data obtained with the 100-m telescope of the Max-Planck-Institut für Radioastronomie at Effelsberg.The Submillimeter Array near the summit of Maunakea is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica.The authors recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community.We are most fortunate to have the opportunity to conduct observations from this mountain.This research has made use of the XRT Data Analysis Software (XRT-DAS) developed under the responsibility of the ASI Science Data Center (SSDC), Italy.This work has made use of Fermi-LAT data supplied by Kocevski et al. (2021) at https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/.This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Figure 1 .
Figure 1.Multifrequency radio light curves of OJ 287 between 2015 December and 2022 June obtained with the Effelsberg telescope in the course of the MOMO program.Note that some receivers and/or frequencies have slightly changed in the course of the monitoring.See Tab. 1 for details.
7. RADIO FLUX AND SPECTRAL VARIABILITY 7.1.Radio light curves The 2015-2022 radio light curves are shown in Fig. 1.Several bright radio flares are present.The brightest occurred in 2016-2017 with peak in February 2017 at a flux density of 10

Figure 3 .
Figure 3. Average fractional variability amplitude between 2.6 and 43 GHz during 2015 December to 2022 June.(Note the scarcity of data at 19 and 24 GHz during the outburst epoch 2016-2017 that partially explains the lower average value of Fvar at these two frequencies).

Figure 4 .
Figure 4. Radio SEDs of OJ 287 between 2015 and 2022 taken at the Effelsberg telescope.The upper two panels cover the bright 2016/2017 outburst between MJD 57816 and MJD 58062 [left panel: rise phase until the main peak in 2017 February (black; decline phase overplotted in green for comparison), right panel: decline phase until the minimum in 2017 November (black; rise phase overplotted in green)].The bottom left panel covers the time interval 2017 November to 2019 January (exact dates: MJD 58076-58510), and the bottom right panel covers the most recent bright radio flare during the time interval 2021 May to 2022 June (exact dates: MJD 59359-59750).See Komossa et al. (2022a) for the time intervals MJD 58511 -58784, MJD 58799 -59028, and MJD 59034 -59426, that are not repeated here.ture function (SF) analysis (following Gallo et al. 2018) of the 10.45 GHz light curve, that yielded an index of β = 0.97 for timescales of 240 days, and subsequently taking α = β + 1.The artificial light curves were then used to compute artificial DCFs with the band-of-interest light curve, allowing for the computation of 90 th , 95 th , and 99 th percentiles based on the distribution of artificial DCFs at each time step.To evaluate the error on the measured lags, we applied a least-squares minimization of a Gaussian function to the central peak, defined as the points with correlation strength corresponding to ≥ 50% of the peak value, within τ = ±400 days, which enabled measurement of the DCF centroid.To evaluate the 68 per cent

Figure 5 .
Figure 5.The S10.45GHz −Sν DCF (black solid line) was computed for each flux density between ν=2.6 and 36 GHz as labelled in each panel.Filled regions indicate the ±90 th (light blue), ±95 th (medium blue), and ±99 th (dark blue) percentiles from the N = 10 3 light-curve simulations.Horizontal dotted lines indicate DCF = 0, and vertical dotted lines mark τ = 0 days.Negative time-axis values indicate S10.45GHz leading the other band, positive values indicate lagging.The vertical red line indicates the measured time lag, with filled regions corresponding to the 68 per cent confidence interval of the lag.In a next step, we have re-run all analyses on the two brightest outbursts of 2016/2017 and 2021/2022 only, in the time interval MJD 57370-58077 and MJD 59426-59751, respectively.Overall, the two events show similar time lags.Longer opacity delays at low frequencies are found for the 2016/2017 outburst with δt=-83 +11 −12 days (2.6 GHz) and -38±9 days (4.85 GHz).

8.
MWL LIGHT CURVES MWL light curves [Fermi γ-rays, Swift X-rays, Swift UV-W2 (representative for the 3 optical and 3 UV bands that are always strongly correlated, even though the spectral index does change)] are shown in Figs. 7, 10, 11 (see Komossa et al. (2021c) for the full Swift long-term light curve including all optical and UV bands between 2005 and 2021 March).

Figure 6 .
Figure 6.DCF lag versus radio frequency.The violinshaped regions represent the distribution of the centroid measurements at each frequency.The reference band is 10.45 GHz.

Figure 7 .
Figure7.MWL flux light curve of OJ 287 between 2015 December and 2022 June.From top to bottom: Fermi γ-rays, Swift X-rays, Swift UV (at 1928 Å), and Effelsberg and SMA radio observations at selected frequencies (green circles: 10.45 GHz, red circles: 32-36.25 GHz, black crosses: SMA data between 213 and 351 GHz).The γ-ray flux (observed 0.1-100 GeV band; one-week averages), the absorption-corrected X-ray flux (observed 0.3-10 keV band), and the extinction-corrected UV flux at λ obs =1928 Å are given in units of 10 −11 erg s −1 cm −2 .The radio data are reported as flux densities in Jy.The gap in Swift observations between June to September each year is due to the close proximity of OJ 287 to the Sun, such that it becomes unobservable with Swift.Fermi and Effelsberg observations are not affected by this Sun constraint.Note that the single brightest γ-ray data point from 2015 December is off the scale (see Fig.10 instead).

Figure 8 .
Figure8.Spectral indices αν in selected optical and radio bands (lower panels; Sν ∝ ν αν ), and the Swift V flux at 5468 Å in units of 10 −11 erg s −1 cm −2 for reference (upper panel).The optical-UV spectral index is measured between V (5468 Å) and UVW2 (1928 Å).The double-peaked 2016/2017 outburst is marked in red.The first maximum is reached in 2016 February, the second and brighter one in 2016 October.
A. ACF The autorcorrelation function (ACF) for the reference band of 10.45 GHz is shown in Fig. 9. B. MWL LIGHT CURVES MWL light curves are shown here in higher resolution, zooming on selected time intervals around the main outbursts (Fig. 10 and 11).
C. RADIO DATA Effelsberg and SMA flux density measurements are reported in Tab.6 and 7.The full tables will be available in the online material.

Table 1 .
Receivers used in our Effelsberg 100 m telescope observations of OJ 287 in the course of the MOMO program since 2015 (Effelsberg program IDs

Table 2 .
Log of Swift observations since 2022 February 19 with

Table 3 .
Fractional variability amplitude Fvar of the radio flux densities of OJ 287.Note that some frequencies were not or barely covered during selected years.When none or too few data points were available, no entry in the table is given; when the value is based on just few data points, the table entry is followed by ':'.

Table 4 .
Observed range of spectral indices.

Table 5 .
Time lags with respect to the flux density at 10.45 GHz.Negative values mean S10.45 is leading.
radio frequency in GHz, (2) modified Julian date, (3) flux density in Jansky and (4) its error.The full table is available in the online material.