Introduction

Blazars are a subclass of radio-loud active galactic nuclei (AGNs). They are further divided into two subclasses: flatspectrum radio quasars (FSRQs), with strong emission lines (e.g., R. D. Blandford & M. J. Rees 1978;G. Ghisellini et al. 1997), and BL Lacertae objects (BL Lacs), which have either no emission lines or very weak (equivalent width < 5 Å) emission lines (J. T. Stocke et al. 1991; M. J. M. Marcha et al. 1996). High brightness, high polarization, and extremely variable emission that is mostly nonthermal, spanning the whole electromagnetic (EM) spectrum, are the main characteristics of blazars. Typically, the emission is ascribed to the relativistic jet that is pointed near the line of sight (LOS) of the observer (C. M. Urry & P. Padovani 1995). Their multiwavelength (MW) spectral energy distribution (SED) is a double-humped structure. The low-energy hump, which is caused by synchrotron emission from nonthermal electrons in the jet, peaks somewhere in the infrared (IR) to soft-X-ray energy range, while the high-energy hump peaks in GeV to TeV γ-ray energies and is likely caused by inverse-Compton (IC) upscattering of synchrotron (synchrotron self-Compton) or external (external Compton) photons by the relativistic electrons responsible for producing the synchrotron emission (J. G. Kirk et al. 1998;H. Gaur et al. 2010).

Blazars are one of the best examples of persistent and highly variable but noncatastrophic sources in the era of MW transient astronomy. Studying the changes in the flux variability of blazars is a valuable way to uncover the physical processes behind the source's various states-whether low, high, or during outbursts. Simultaneous MW studies have been carried out in order to understand their emission mechanism spanning the whole EM spectrum (e.g., C. M. Urry et al. 1997; F. Aharonian et al. 2005Aharonian et al. , 2009; C. M. Raiteri et al. 2007Raiteri et al. , 2008Raiteri et al. , 2015;;S. Vercellone et al. 2009S. Vercellone et al. , 2010;;M. Villata et al. 2009;A. C. Gupta et al. 2017;A. Goyal et al. 2018;P. Kushwaha et al. 2018a;S. Komossa et al. 2020;MAGIC Collaboration et al. 2024, and references therein).

The BL Lac OJ 287 (α 2000.0 = 08 h 54 m 48. s 87, d = 2000.0 +  ¢  20 06 30. 64) is at redshift z = 0.306 (M. L. Sitko & V. T. Junkkarinen 1985). OJ 287 has been observed in optical bands since 1888 (M. J. Valtonen et al. 2024). A small fraction of the light curve was already available in 1982, when it was noticed that OJ 287 may exhibit a nearly periodic outburst about every 12 yr. The next outburst was expected in 1983 and it was indeed detected (A. Sillapää et al. 1988). The authors postulated a supermassive binary black hole (SMBBH) model to explain the 12 yr periodicity and predicted that the next outburst would take place in late 1994. A. Sillapää et al. (1988) also noted a possible shorter periodicity in the fades, the times of minimum light. Assuming that the difference in the periodicities arises from the procession of the major axis of the binary, A. Sillapää et al. (1988) calculated that the primary BH's mass was ∼5 × 10 9 M e , while the secondary's mass was estimated from the rapid variability over a 15.7 minutes timescale as ∼2 × 10 7 M e (E. Valtaoja et al. 1985). The anticipated outburst was observed in 1994, thanks to a global optical monitoring campaign of the source known as OJ-94 (A. Sillanpää et al. 1996). However, H. J. Lehto & M. J. Valtonen (1996) predicted that the outbursts should have a double-peaked structure and that the second peak should take place within a two-week interval in 1995 October. It was immediately verified by observations (A. Sillapää et al.1996).

B. Sundelius et al. (1997) calculated the binary model forward, to predict the next pair of outbursts in 2005 November and 2007 September. The increase of the two-flare interval is due to the orbit procession in the model and it improved the primary mass to ∼1.7 × 10 10 M e . Both flares were seen at expected times (M. Valtonen & A. Sillanpää 2011). Another set of flares, this times a triple set, was predicted for the years 2015, 2019, and 2022(B. Sundelius et al. 1997)). The model showed that the timing of the first flare was sensitive to the spin value of the primary. After it was observed, the spin value was calculated (M. J. Valtonen et al. 2016). The timing of the second flare was very precise (S. Laine et al. 2020). L. Dey et al. (2018) have developed a highly accurate SMBBH model that can forecast the times of the flares to within 4 hr. The last of the triple flares was not observable from the ground, since it was expected when OJ 287 was very close to the Sun. However, it was possible to infer the presence of the third flare from particular preflare activity (M. J. Valtonen et al. 2023). The BH binary model of L. Dey et al. (2018) yields the following values for OJ 287: primary BH mass = (18.35 ± 0.05) × 10 9 M e ; and secondary BH mass = (150 ± 10) × 10 6 M e .

There are many claims of detections of quasiperiodic oscillations (QPOs) from OJ 287 on a wide variety of timescales, from a few tens of minutes to decades and more, over multiple EM bands, aside from the well-established 12 yr and 55 yr periodicities in the optical band (M. J. Valtonen et al. 2006). N. Visvanathan & J. L. Elliot (1973) reported for the first time the detection of a ∼40 minutes optical QPO in OJ 287 using accurate optical photoelectric observations on 1972 March 18. Later, a few more optical QPOs were reported, with periods ranging from 23 to 40 minutes (A. Frohlich et al. 1974;L. Carrasco et al. 1985). In 1981 April, following observation of the source in the 37 GHz radio band, a ∼15.7 minutes QPO was reported (E. Valtaoja et al. 1985). Using recent advanced techniques, there are more claims of detections of QPOs in OJ 287 in different EM bands, on diverse timescales ranging from a few tens of days to months to years, over the different time spans of the data (e.g., P. Pihajoki et al. 2013;G. Bhatta et al. 2016;S. Britzen et al. 2018;P. Kushwaha et al. 2020, and references therein).

OJ 287 has been observed simultaneously in various flux, spectral, and polarization states, on several occasions with diverse timescales (e.g., H. Siejkowski & A. Wierzcholska 2017;A. Goyal et al. 2018;P. Kushwaha et al. 2018aP. Kushwaha et al. , 2018b;;P. Kushwaha 2020;S. Komossa et al. 2020S. Komossa et al. , 2021a;;R. Prince et al. 2021, and references therein). The source has shown a major γ-ray flare in a Fermi observation of 2009, which was studied to understand the high-energy emission mechanism during this episode (P. Kushwaha et al. 2013). The extensive X-ray flux and spectral variability of OJ 287 have been studied on several occasions, using various X-ray and MW space missions, and variabilities have been found on diverse timescales (e.g., E. Idesawa et al. 1997;N. Isobe et al. 2001;B. Kapanadze et al. 2018;P. Kushwaha et al. 2018b;M. Pal et al. 2020;S. Komossa et al. 2021aS. Komossa et al. , 2021b;;M. Mohorian et al. 2022;K. P. Singh et al. 2022;D. Zhou et al. 2024, and references therein).

When observing blazars at multiple epochs, simultaneous MW SEDs provide valuable information about the emission mechanisms of blazars and their various flux levels (e.g., R. M. Sambruna et al. 1996;E. Massaro et al. 2004E. Massaro et al. , 2006;;E. Nieppola et al. 2006;F. Massaro et al. 2008;B. Rani et al. 2011;J. Bhagwan et al. 2014; N. Sahakyan 2021; N. Sahakyan & P. Giommi 2022; N. Sahakyan et al. 2022, and references therein). Modeling the broadband SEDs of blazars is essential to understanding the extreme conditions within different emission regions. This approach helps us comprehend the dynamic phenomena shaping the observed behavior of blazars. In the ideal case, such studies require simultaneous data in multiple bands. In the present paper, by utilizing comprehensive data spanning radio, near-IR (NIR), optical, and ultraviolet (UV) bands for OJ 287, we construct multi-epoch flux-statespecific SEDs from nearly simultaneous observations, strictly maintaining temporal intervals of up to 10 days.

We describe the observations and data in Section 2 and the SED modeling in Section 3. The results are delivered in Section 4 and discussed in Section 5. We summarize our main results in Section 6. Throughout the paper, a flat ΛCDM cosmology with Ω Λ = 0.7, Ω m = 0.3, and H 0 = 70 km s -1 Mpc -1 is adopted.

Observations and Data

Multiband radio, NIR, optical, and UV data of the blazar OJ 287 are collected for the period of 1998-2023 from various public archives and observing facilities. The details of the data are provided in Table 1.

UVOT is one of the instruments on board the Swift observatory, capable of observing in six filters, namely V, B, U, w1, m2, and w2, covering optical to UV regions of the EM spectrum. We used all the observation IDs from 2005 to 2023 and analyzed them following the standard data reduction prescription, as mentioned in P. Kushwaha et al. (2021) and P. Kushwaha (2023).

Optical V-and R-band photometric observations of OJ 287 are obtained from the spectropolarimeters mounted on the 2.3 m Bok and 1.54 m Kuiper telescopes at Steward Observatory, University of Arizona, USA. OJ 287 data from 2008 October to 2018 June are taken from the public archive of the Steward Observatory. 14 The details of the instrument, observational program, observations, and data analysis procedures are provided in P. S. Smith et al. (2009).

Optical B-, V-, R-, and I-band photometric observations of OJ 287 were carried out from 2006 January to 2023 February at the Perkins telescope of the Perkins Telescope Observatory (Flagstaff, AZ, USA). The details of the instrument, observations, and data analysis methods are given in S. G. Jorstad et al. (2010).

Optical B-, V-, and R-band data as well as NIR J-and K-band data for OJ 287 are taken from the public archive of the Small and Moderate Aperture Research Telescope System (SMARTS) from 2008 February to 2017 April. 15 SMARTS consists of 0.9 m, 1.0 m, 1.3 m, and 1.5 m telescopes at the Cerro Tololo Inter-American Observatory in Chile. These telescopes observed the blazars at both NIR and optical wavelengths that Fermi-LAT monitors. The SMARTS telescopes, detectors, observations, and data analysis details are provided in E. Bonning et al. (2012) and M. M. Buxton et al. (2012).

The J, H, and K s NIR-band observations of OJ 287 from 1995 October to 2021 November were carried out with the 2.12 m telescope that is equipped with an NIR camera named the Cananea Near-Infrared Camera of the Guillermo Haro Astrophysical Observatory (OAGH), located in Cananea, Sonora, Mexico. The details of the instrument, observations, and data analysis procedures are provided in, e.g., J. A. Cardelli et al. (1989), L. Carrasco et al. (2017), and A. C. Gupta et al. (2022), while the photometric data have already been published in A. C. Gupta et al. (2022).

The University of Michigan Radio Astronomy Observatory (UMRAO) flux density data of OJ 287 at 4.8, 8.0, and 14.5 GHz from 2007 November to 2012 June are obtained from the Michigan 26 m equatorially mounted, prime-focus paraboloid, as part of the University of Michigan extragalactic variable-source monitoring program (H. D. Aller et al. 1985). The radio data of OJ 287 at 15 GHz are taken from the blazar monitoring program of the 40 m telescope of the Owens Valley Radio Observatory (OVRO) for the period from 2008 January to 2023 August. The details of this observational program, observations, and data analysis procedures are provided in J. L. Richards et al. (2011).

Using the 14 m radio telescope at Aalto University Metsähovi Radio Observatory in Finland, observations of OJ 287 at 37.0 GHz were conducted. H. Teräesranta et al. (1998) provided a thorough explanation of the Metsähovi data reduction and analysis process.

The Very Long Baseline Array (VLBA)-Boston University (BU) BLAZAR monitoring effort involves about monthly VLBA observations of a sample of AGNs identified as γ-ray sources at 43 GHz and 86 GHz. The observations and data analysis of OJ 287 at 43 GHz and 86 GHz are presented in detail (S. G. Jorstad et al. 2017; Z. R. Weaver et al. 2022, and references therein). OJ 287 is a very compact core-dominated source at radio wavelengths, especially at high radio frequencies, such as 43 and 86 GHz. As described in S. G. Jorstad et al. (2017), for each epoch, we calculated the total flux density in the images of several sources in the sample that are known to have very weak emission outside the angular size range of the VLBA images (0235+164, 0420-014, 0716+714, OJ 287, and 1156+295). These values were compared with the total flux densities obtained by interpolating in time the measurements of these sources by monitoring programs carried out at the Very Large Arrayfoot_3 and the Effelsberg telescope at 43 GHz

SED Modeling

The observed SED covering the UV to radio bands was modeled with a parabola in the logarithms of the variables (hereafter, a log-parabola or LP, for short). The simplified model assumes that radiation comes from a single region in the jet, filled with chaotic magnetic fields and electrons, moving relativistically at a small angle to the observer's LOS. Note that for blazars, the location of the radio core varies significantly with frequency, particularly across the range of 4.8-86 GHz, which we will use in this work. However, this variation is generally smaller for BL Lacs. In the specific case of OJ 287, A. B. Pushkarev et al. (2012) estimate that the 15 GHz core is located within 4.1 pc of the BH, with the positional difference between the 15 and 8 GHz cores being less than 0.05 mas. Consequently, all cores at frequencies higher than 15 GHz should lie within 4 pc of the BH. Given this, if the emission region spans about 4 pc, it can reasonably be treated as a single region for modeling purposes. As a result, the observed radiation experiences Doppler boosting, described by the Doppler factor

, where β is the velocity of the source divided by the light velocity, Γ is the Lorentz factor, and θ is the angle between the LOS of the observer and the direction of motion of the source.

An LP distribution is not only a simple mathematical tool for spectral modeling, but also relates to the physics of the electron acceleration processes. Both the statistical and stochastic acceleration mechanisms can reproduce the electron energy distribution as an LP law, resulting in an LP SED approximately (E. Massaro et al. 2004;A. Tramacere et al. 2007A. Tramacere et al. , 2011;;F. Massaro et al. 2008; L. Chen 2014, and references therein). The LP function for SED modeling has three spectral parameters and can be defined as

log log log p log p p, 1 2 where b measures the curvature around the SED peak, νp is the peak frequency, and n n f log p p is the peak flux (B. Rani et al. 2011; L. Chen 2014; A. C. Gupta et al. 2016; J. H. Yang et al. 2022).

The statistical acceleration mechanism framework requires either an energy-dependent acceleration probability (p a ) or variations in the fractional acceleration gain (ò). Studies by E. Massaro et al. (2004) and F. Massaro et al. (2008) demonstrate that an LP spectrum can be obtained when the probability of particle acceleration is energy-dependent. This scenario naturally occurs when particles are confined by a magnetic field whose efficiency decreases as the gyration radii of the particles increase (B. Rani et al. 2011). Additionally, in cases where there are fluctuations in the energy gain parameter ò, an LP spectrum can also form under specific conditions if ò is treated as a random variable centered around a systematic value (A. Tramacere et al. 2011).

Moreover, an LP spectrum can result from the stochastic acceleration mechanism, described by the Fokker-Planck equation with an included momentum diffusion term (A. Tramacere et al. 2007Tramacere et al. , 2011)). In this framework, an LP distribution of electron energy can be derived from a "quasimonoenergetic" injection (N. S. Kardashev 1962).

By maintaining temporal intervals of up to 10 days, we successfully constructed 106 SEDs spanning from UV to radio bands. The choice of a 10 days interval is primarily motivated by the need to balance the quantity of SEDs and the simultaneity of the MW data comprising these SEDs. This choice allows for a 10% to 18% increase in the number of SEDs compared to intervals of 4-8 days. However, extending the interval beyond 10 days yields less than a 5% increase in SEDs, while compromising the simultaneity across different data filters. Additionally, 10 days correspond to the typical observational window for OJ 287 during a month, especially around the new moon.

These SEDs cover the MJD range from 54850 (2009-01-19) to 59227 (2021-01-13). The durations of the MW data in each band, used for constructing the SEDs in this work, are listed in Columns (4) and (5) of Table 1. The numbers of data points of each band within the applied durations are also shown in Column (6), from which we select the data points for constructing the SEDs. Each SED includes at least one data point in the following seven series of bands:

(i) UV bands (w2, m2, w1, and u);

(v) partial optical plus partial NIR bands (I, J, and H); (vi) additional NIR bands (K and Ks); and (vii) radio bands (86 GHz, 43 GHz, and 37 GHz).

For each band within each series of bands, if multiple measurements are available from one observatory or different observatories, the final flux for that band is calculated as the median of these measurements. Among the 106 SEDs, 72 SEDs includes at least one data point in additional radio bands (15.0 GHz, 14.5 GHz, 8.0 GHz, and 4.8 GHz). Galactic extinction correction was performed for the data in the NIR to UV bands (J. A. Cardelli et al. 1989;D. J. Schlegel et al. 1998), and redshift correction was subsequently performed for the constructed SEDs.

We fit the SEDs using the LP model with the maximum likelihood method, which minimizes the negative log-likelihood. This is implemented using the optimize.minimize function (P. Virtanen et al. 2020). The negative log-likelihood function is

where y means the observed n n f log values, y model means the predicted y value obtained from the model shown in Equation (1), and log_f is a parameter representing an additional scatter beyond the measurement error y err . The total uncertainty, σ 2 , is calculated as the square of the measurement uncertainty, y err 2 , plus an additional term that scales with the model value and the exponential of log_f. With log_f > 0, the model uncertainty allows for more flexibility to account for additional scatter not captured by y err . In practice, incorporating log_f during the fitting helps balance between underfitting, by providing too little flexibility, and overfitting in the model.

The significance of our SED model fitting is first evaluated using the reduced χ 2 , the sum of the squared normalized residuals divided by the degrees of freedom. In our case, the reduced χ 2 is greater than 1, indicating larger residuals than expected from the uncertainties. For further investigation, we used the probability based on log_f as a model-adjusted flexibility measure to account for additional scatter. The log_f values, ranging from 0 to 0.013, with a median of 0.005, result in a small increase in the total uncertainties σ 2 of up to 2.6%, with a median of 1.0%. These low-log_f values confirm that the larger reduced χ 2 values are likely caused by slightly underestimated measurement uncertainties. Overall, the fitting remains generally significant.

A Monte Carlo approach is applied to estimate the uncertainties of the fitting parameters. For each SED, 50 random mock SEDs are generated, by introducing Gaussian noise into the original SED. At each frequency n log , in a given mock SED, the noise term is randomly drawn from a normal distribution, with the observed n n f log error as the standard deviation. We then fit each mock SED using the same fitting strategy. The 1σ dispersion of the measurements relative to the original values is taken as the corresponding uncertainty. Together with visual check, we define SEDs with b > 0.02 as well fit by the LP function. All the 106 SEDs can be well fit by the LP model, as shown in Figures 1 and 2.

Between MJD 55558 and 55621, although the K-band data points of six SEDs deviate most significantly from the model fits, our analysis shows that their inclusion does not significantly impact the overall results, as comparisons of fits with and without these data points reveal little difference in the derived parameters. The relative bump in the K s band is generally attributed to thermal emission from dust at a wide range of temperatures (B. Wilkes 2004), contribution of the host galaxy of the blazar, and IR contribution of the torus (P. Giommi et al. 2024).

With the number of data points in radio bands of each SED shown in the upper-left corner of each panel in Figures 1 and 2, we investigate whether the obtained b values are influenced by the number of data points in radio bands. We find that only 12 epochs have a single data point in radio. When all the 106 SEDs are ranked in order of decreasing b, none of the 30 highest epochs have only a single radio data point, but seven of the lowest epochs do. To further check the influence of the number of radio data points, as 105 SEDs except for the first SED include data points at 37 GHz, we reduce the number of radio data points at all epochs to one, i.e., the data point at 37 GHz is selected, if available, otherwise the closest data point in frequency at 43 GHz is chosen. This simplification results in updated b values ranging from 0.034 to 0.210, with a median of 0.104 ± 0.043, as shown by the green histogram in Figure 3. For comparison, the original b values ranged from 0.038 to 0.234, with a median of 0.117 ± 0.045, as shown by the gray histogram. The difference between the median values is -0.013, smaller than the standard deviation around the median. The Kolmogorov-Smirnov test yielded a test statistic of 0.18 and a p-value of 0.06, indicating that at the 5% significance level, there is no statistically significant difference between the two distributions.

Results

Flare and Quiescent States

As shown in the upper panel of Figure 4, the target OJ 287 shows optical variability in the V band across the MJD range of 54000-60000. For almost three months each year, OJ 287 is not visibly accessible to the optical telescopes used to collect the data for this study. By analyzing the V-band flux distribution of OJ 287, we identified a distinct log-normal profile, as shown in Figure 5. We first determined the optimal number of Gaussian components using the Bayesian Information Criterion, which indicated that a single Gaussian component was most appropriate. We then fit a Gaussian Mixture Model using this optimal number of components and extracted the mean and standard deviation (σ) from the fitted profile. Using these parameters, we established a flux limit based on the mean plus half the σ of the distribution, resulting in a value of 10 -10.49 erg cm -2 s -1 (shown as the right edge of the green region in Figure 5 and also as the horizontal green line in the upper panel of Figure 4). The cumulative distribution function at the flux limit is 0.69, indicating the probability that a randomly selected sample will have a value less than or equal to the flux limit.

We defined a "flare" segment in the V-band light curve as any observation period containing more than three consecutive data points with flux exceeding a specified limit, with segments not meeting this criterion designated as "quiescent." Testing variations from one to six consecutive data points revealed that the number of flare segments fluctuated only slightly, by two to four segments, without affecting the number of SEDs within the flare segments or the duration proportion of the flare segments. This indicates that the choice of consecutive data points does not influence the subsequent analysis of the SEDs in flare versus quiescent segments. Our choice of three consecutive data points strikes a balance, by minimizing misclassification from isolated outliers, ensuring genuine flare detection, and maintaining enough segments for meaningful analysis, making it an optimal threshold. The flare segments in the V band are shaded in green in Figure 4; the start and stop dates of the individual flare segments are listed in Table 2.

This categorization results in 19 flare segments, with durations ranging from 1 to 312 days. As shown in Figure 4, there are no constructed SEDs available for 10 of the 19 identified flare segments, including the first and longest flare segment, with a duration of 312 days. This absence is partly due to the scarcity of data points in the UV to radio bands.

The flare segments summarized in Table 2 differ from those discussed in the introduction, which are used for orbit determination. Only two of these segments-SEDs with central MJDs of 57362 and 57367 (2015 December 6 and 11)coincide with the time range of the predicted and confirmed flare in 2015, as shown in the fifth row and third to fourth columns of Figure 2 and as listed as the thirteenth row in Table 2. Both flares are exceptionally bright and exhibit rapid variability, requiring higher temporal resolution for studying their spectral changes (M. J. Valtonen et al. 2016). Although they could have been excluded from the adopted flare segments, their inclusion as two single epochs does not affect the results of this study.

The total duration of the flare segments spans 954 days, representing 15% of the 6475 days time span of the V-band light curve analyzed in this work. Among the 106 SEDs constructed from nearly simultaneous multiband photometric data, 30 SEDs occur during flare segments (the modeled SEDs in green in Figures 1 and 2), while 76 SEDs are in quiescent states (the modeled SEDs in gray).

SEDs at Different States

We compared the SEDs in the flare and quiescent segments based on three key parameters: the peak intensity ( n n f log p p), the SED curvature (b), and the peak frequency ( n log p). As shown in Figure 6, the median n n f log p p for flare segments is 0.37 ± 0.22 dex higher than for quiescent segments. The median curvature b is slightly larger in flare segments (0.14) compared to quiescent segments (0.11). However, this difference is negligible when considering the uncertainty of the median value (∼0.04). Similarly, n log p values remain consistent, with 14.02 ± 0.70 for flare segments and 13.95 ± 0.79 for quiescent segments, respectively. Here the uncertainties of the median values are derived from the 1σ dispersion of the distribution of the corresponding parameter.

Color Variability

As variations in the optical flux of blazars are accompanied by spectral changes, studying the color index-magnitude (CM) relation can help in understanding the origin of the variability in blazars. Earlier studies have found significant bluer-whenbrighter (BWB)/redder-when-brighter (RWB) and achromatic trends on diverse timescales on the CM diagram (e.g., M. F. Gu et al. 2006;H. Gaur et al. 2012;A. Agarwal et al. 2016A. Agarwal et al. , 2019A. Agarwal et al. , 2021, and references therein), and references therein).

Due to the potential for non-negligible magnitude fluctuations when switching filters during nonsimultaneous observations, making accurate color measurements difficult, it is necessary to obtain very dense and precise simultaneous multiband observations to detect weak CM relationships. Based on the V-band magnitudes and B -V color indices of the 106 SEDs, as shown in Figure 7, we find there is a weak BWB relation, with a Spearman correlation coefficient r ∼ 0.28 at a confidence level above 99.6%, for the SEDs. Further excluding the three outlier points with significant B -V difference and constraining 0 < B -V < 0.6, we achieved r ∼ 0.26 with a confidence level exceeding 99.2%. This is also confirmed by a weak anticorrelation between the SED peak frequency n log p and the V-band magnitude (Figure 8), showing the peak frequency being higher at brighter magnitude, i.e., r ∼ -0.19 at a confidence level above 94.6%.

Optical emission from blazars typically consists of contributions from both the relativistic jet and the accretion disk, with the jet often being dominant. When the synchrotron radiation from the relativistically boosted jet outshines the emission from the disk, the BWB trend can be attributed to either the acceleration of relativistic particles or the injection of fresh electrons with an even harder energy distribution (J. G. Kirk et al. 1998;A. Mastichiadis & J. G. Kirk 2002;M. Fiorucci et al. 2004;A. C. Gupta et al. 2016). For the RWB trend, the contribution of the accretion disk to the total emission could be significant. In general, BWB and RWB trends were found in BL Lacs and FSRQs, respectively (e.g., H. Gaur et al. 2012;A. Agarwal et al. 2019A. Agarwal et al. , 2021, and references therein), but sometimes the opposite trend is also noticed (e.g., H. Gaur et al. 2012, and references therein).

In Figure 7, we find a stronger BWB trend during flares (green symbols) compared to quiescent states (gray symbols), indicating the dominance of the jet over the accretion disk in the flare segments. In the flare and quiescent segments, the correlation coefficient is r ∼ 0.40 at a confidence level over 97.3% and r ∼ 0.20 at a confidence level above 91.1%, respectively. This pattern is corroborated by anticorrelations between n log peak and the V-band magnitude (Figure 8), with r ∼ -0.53 and a confidence level above 99.7% in flares, compared to r ∼ -0.30 and a confidence level over 99.2% in quiescent segments.

The Doppler factor variations are also usually attributed to achromatic behavior, and this interpretation is most likely supported by the geometric scenario (e.g., M. Villata et al. 2002). I. Liodakis et al. (2021) estimated the Doppler factor versus frequency in log-log space for 61 blazars, including OJ 287. They used data from five radio bands from 4.8 to 37 GHz and found there was a linear relation with slope 0.22 - + 0.29 0.29 , intercept 1.07 - + 0.35 0.32 , and a Pearson correlation coefficient of 0.58. This linear relation may be extended from the NIR to UV bands to estimate the Doppler factor in these EM bands.

Discussion

LP SEDs and Statistical Particle Acceleration

The study of LP SEDs in blazars has uncovered significant correlations. For a sample of 60 blazars, radio-to-X-ray SEDs were well fitted by the LP model, where the peak frequency was found to be anticorrelated with the bolometric luminosity (R. M. Sambruna et al. 1996). In contrast, for a sample of 300 BL Lacs, SEDs from radio to X-ray also fit the LP model, showing an anticorrelation between the peak frequency and flux at radio (5 GHz) and optical (5500 Å), but no such anticorrelation was observed with the X-ray flux (E. Nieppola et al. 2006).

Additional studies have also explored the connection between peak frequency and curvature. By fitting the SED from radio to optical with an LP model for a sample of 18 blazars, R. Landau et al. (1986) found an anticorrelation between the peak frequency and curvature for the 15 blazars that can be well fit by an LP model. Similar results were obtained in recent works (B. Rani et al. 2011;L. Chen 2014;J. H. Yang et al. 2022).

Mainly, there are two different scenarios explaining the correlation between the peak frequency and curvature. The first scenario is within the framework of statistical acceleration. For the case of the energy-dependent acceleration probability (p a ), E. Massaro et al. (2004) showed that if p a is inversely related to the particle's energy, the resulting SED naturally adopts an LP form, where the curvature b can be inversely correlated with the peak energy or frequency ( n log p), described by / µ b 1 / n 5 2 log p. In contrast, considering the fluctuations in the fractional energy gain (ò), A. Tramacere et al. (2011) demonstrated that treating ò as a random variable around a systematic energy gain also leads to an inverse relationship between b and n log p, following the relation / µ b 1 / n 10 3 log p.

The second scenario is within the framework of the stochastic acceleration mechanism, which can predict an anticorrelation between b and n log p , described by the relation Tramacere et al. 2007Tramacere et al. , 2011)). Using a sample of 10 low-to-intermediate synchrotronpeaked blazars, B. Rani et al. (2011) found an anticorrelation between b and ν p and suggested that the LP SED shape is likely characterized by a full statistical acceleration mechanism acting on the emitting electrons, while using a large sample of 48 blazars, L. Chen (2014) found that the slope of the correlation between 1/b and ν p as 2.04 ± 0.03 is consistent with the prediction of the stochastic acceleration scenario ( ∼2). This is further confirmed by M. S. Anjum et al. (2020), who found that BL Lacs show a strong signature of stochastic acceleration compared to FSRQs.

In Figure 9, we observe a strong anticorrelation between b and n log p, with r = -0.95 at a confidence level above 99.9%, as shown by the green and gray points. This tendency is also evident in the lower two panels of Figure 4, where both 1/b and n log p vary with the MJD values in similar patterns. By performing a linear fit to 1/b and n log p, we derive the relation

6.20 0.08 log p 77.82 1.03) (represented by the solid line). If we exclude the seven data points in the upper-right corner with significant 1/b differences and constrain 1/b < 22, the slope decreases to 5.79 ± 0.06 (represented by the dashed line). This revised slope aligns more closely with the predicted value of 10/3 from the statistical acceleration mechanism, which accounts for fluctuations in the fractional acceleration gain ò, though a notable discrepancy remains, suggesting the influence of additional factors.

The observed slope is significantly steeper than both the previously reported value of 2.04 by L. Chen (2014) and the theoretical prediction of 10/3 (A. Tramacere et al. 2011). Observationally, L. Chen (2014) derived the slope using fewer simultaneous data spanning a broad wavelength range (radio to γ rays) from 48 blazars, including both BL Lacs and FSRQs, while our study focuses on the single blazar OJ 287, with nearly simultaneous data but limited to a narrower range (radio to UV bands). Moreover, the steeper slope in our results may stem from uncertainties in estimating the SED peak frequency and curvature, due to sparse data between the radio and UV bands. A much closer estimate of the observed slope to the theoretical one may be achieved with a large number of blazar SEDs with much denser data coverage in frequency and time in the future.

Nevertheless, the discrepancy with the theoretical predictions may reflect additional physical factors beyond the standard statistical acceleration mechanism, which primarily considers electron acceleration processes in the jet. For example, deviations from idealized conditions or radiative contributions, such as thermal emission from the accretion disk in the optical-UV bands, especially during quiescent states, could contribute to the observed steeper slope. These results highlight the importance of incorporating additional complexities and/or exploring alternative explanations, rather than strictly adhering to the standard statistical acceleration mechanism.

The Cause of SED Changes

Under the frame of the statistical acceleration mechanism, there are several possible reasons that can explain the changes in the low-energy-peak SEDs along with time. If the changes in the SEDs are primarily caused by a gradual change in the electron energy density distribution, due to the synchrotron and IC losses, with no other injections during the period, one would There is a BWB relation, with r ∼ 0.28, at a confidence level above 99.6%. Further excluding the three outlier points with significant B -V difference, by constraining 0 < B -V < 0.6 (the data points between the two dashed horizontal lines), we achieved r ∼ 0.26, with a confidence level exceeding 99.2%. The BWB trend is stronger during flare segments (green symbols) compared to quiescent ones (gray symbols), with r ∼ 0.40 at a confidence level over 97.3% vs. r ∼ 0.20 at a confidence level over 91.1%. n log p vs. V-band magnitude, where r is -0.19 at a confidence level above 94.6%. The trend is stronger during flare segments (green symbols) compared to quiescent ones (gray symbols), with r ∼ -0.53 at a confidence level over 99.7% vs. r ∼ -0.30 at a confidence level over 99.2%. expect a positive relation between the peak-intensity and the peak-frequency changes (B. Rani et al. 2011). Adopting the epoch where OJ 287 is faintest in the V band among all the considered SEDs as the reference epoch, we calculate the differences in the peak intensity ( n n D f log p p) and the peak frequency ( n D log p) relative to the reference epoch. In Figure 10, we find that there is a significant anticorrelation between n n D f log p p and n D log p, i.e., the Spearman correlation coefficient r ∼ -0.38 at a confidence level above 99.9%. The discrepancy from the prediction suggests that the evolution of the electron energy or electron injection may not be the primary driver of the SED changes. Note that the anticorrelation shows a hint of OJ 287 following the blazar sequence -the anticorrelation between n log p and n n f log p p for a blazar sample (G. Fossati et al. 1998), related to the physical conditions in the jet (G. Ghisellini et al. 1998).

Moreover, the changes of other parameters, including the Doppler-boosting factor and the magnetic field, may also cause the changes in the SED along with time. By studying the correlation between the change in the peak intensity and the change in both the Doppler-boosting factor and the magnetic field for a sample of 10 blazars, B. Rani et al. (2011) conclude that the change in the Doppler factor is a strong driver of the SED changes, whereas the changes in the magnetic field strength may influence only BL Lacs but not all blazars.

It is found that the Doppler factor is substantially higher in the flaring states of blazars, which may cause the strong increase in the Compton dominance, as the external photon density in the comoving frame of the jet depends on the Doppler-boosting factor (N. Sahakyan 2021). Either a fresh injection or reacceleration of the faster-moving emitting region during propagation could cause it to emit near the central source. Geometrical effects, such as when the jet zones have various orientations, such as in the cases of the jets in a jet model (D. Giannios et al. 2009) or twisted inhomogeneous jet model (C. M. Raiteri et al. 2017), could also be the cause of the Doppler-boosting-factor rise. Therefore, the Doppler-boosting factor is increased because during the flares, photons may be emitted in a zone observed at smaller angles than the entire jet.

The effects of magnetic field topology on the SEDs in blazars demonstrate that, in the case of a purely oblique field, the synchrotron component is annulated if the magnetic field is aligned along the LOS (in the plasma frame; M. Joshi et al. 2020). However, the impact of an oblique field is diminished and the same effect is not noticed in the presence of a disordered component (M. Joshi et al. 2020).

In the case of the BL Lac OJ 287, we find that the median n n f log p p for flare segments is 0.37 ± 0.22 dex higher than that for quiescent segments, while the median n log p and b remain consistent within their uncertainties, suggesting that flaring might be mainly caused by the Doppler factor and the ambient magnetic field. This also explains the observed stronger BWB trend in the flare segments compared to the quiescent ones. Therefore, we argue for the possible contribution from variations in the Doppler factor and the magnetic field strength to the observed SED changes. However, we are not able to quantify their contributions based solely on the SED changes.

Conclusion

Using nearly simultaneous radio to NIR to UV data with temporal intervals up to 10 days, we have conducted SED studies of the blazar OJ 287. We have constructed 106 SEDs covering from MJD 54850 (2009 January 19) to 59227 (2021 January 13) and modeled them in the n log -n n f log diagram using an LP synchrotron model. The main results are summarized as follows:

1. All the constructed 106 SEDs can be well fit by the LP model, with curvature b > 0.02. The b values ranged from 0.038 to 0.234, with a median of 0.117 ± 0.045. 2. We classified the observational periods into flare and quiescent segments based on whether the flux values at the V band exceed or fall below 10 -10.49 erg cm -2 s -1 , the mean plus half the standard deviation of the V-band flux distribution. We found that the median flux at the peak frequency of the SEDs during flare segments was 0.37 ± 0.22 dex higher than during quiescent segments, while no significant differences were observed in the median values of the curvature parameter b or the peak frequency n log p. 3. There is a significant relation between the V-band magnitude and B -V color index for the 106 SEDs, confirming a BWB relation. A stronger BWB trend is found in the flare segments compared to the quiescent ones, as further supported by the anticorrelation between the SED peak frequency and the V-band magnitude. 4. We found a significant anticorrelation between the SED curvature b and the peak frequency n log p of the synchrotron component. The slope of the correlation between 1/b and n log p, measured as 5.79, aligns more closely with the prediction of the statistical acceleration scenario than with the stochastic acceleration scenario, though a notable discrepancy persists. This discrepancy indicates that additional factors-such as deviations from idealized conditions or radiative contributions, such as thermal emission from the accretion disk in the optical-UV range during quiescent states-may play a role in producing the observed steeper slope. 5. Within the framework of the statistical acceleration mechanism, we considered potential factors influencing the observed SED changes in blazars. No positive correlation was found between the changes in peak intensity and peak frequency, suggesting that the change in the electron energy distribution is unlikely to be the primary driver. Other factors, such as changes in Doppler-boosting factor and/or magnetic fields, may contribute to the observed SED changes.