Revealing the Variation Mechanism of ON 231 via the Two-component Shock-in-jet Model

The variation mechanism of blazars is a long-standing unresolved problem. In this work, we present a scenario to explain diverse variation phenomena for ON 231, where the jet emissions are composed of the flaring and the less variable components (most probably from the post-flaring blobs), and the variation is dominated by shock-in-jet instead of the Doppler effect. We perform correlation analysis for the multiwavelength light curves and find no significant correlations. For the optical band, ON 231 exhibits a harder when brighter (HWB) trend, and the trend seems to shift at different periods. Correspondingly, the correlation between the degree of polarization and flux exhibits a V-shaped behavior, and a similar translation relation during different periods is also found. These phenomena could be understood via the superposition of the flaring component and slowly varying background component. We also find that the slopes of the HWB trend become smaller at higher flux levels, which indicates the energy-dependent acceleration processes of the radiative particles. For the X-ray band, we discover a trend transition from HWB to softer when brighter (SWB) to HWB. We consider that the X-ray emission is composed of both the synchrotron tail and the synchrotron self-Compton components, which could be described by two log-parabolic functions. By varying the peak frequency, we reproduce the observed trend transition in a quantitative manner. For the γ-ray band, we find the SWB trend, which could be explained naturally if a very-high-energy γ-ray background component exists. Our study elucidates the variation mechanism of intermediate synchrotron-peaked BL Lac objects.


INTRODUCTION
Blazars are a subclass of active galactic nuclei (AGN) with relativistic jets pointing to our line of sight (Urry & Padovani 1995).They display violent flux, polarization, and spectral variability at all wavelengths.According to the optical emission line features, they can be divided into two subclasses, including the flat-spectrum radio quasars (FSRQs) and BL Lacaerte objects (BL Lacs).The spectral energy distributions (SEDs) of blazars exhibit two bumps.It is generally accepted that the low energy bump is dominated by synchrotron radiation, while the high energy bump is produced by the inverse Compton scattering process (Sikora et al. 2009;Böttcher et al. 2013;Romero et al. 2017).The radiation could be from multiple background components such as the accretion disk, dusty torus, and broad line region (Ghis-ellini et al. 2019), and the variation could be affected by many factors including geometric effects, magnetic reconnection, and shock.Thus, the variation mechanism of blazars at various timescales is complex and under intensive debate.A comprehensive study of various aspects of the variation, including flux, spectrum, polarization degree (PD), and polarization angle (PA) can help us constrain the variation mechanism and understand the essential nature of blazars.
ON 231 (also known as W Comae, W Com, 1219+285) is a violently variable BL Lac object in Coma Berenices.It was discovered as a variable star by Wolf (1916) and identified as the counterpart of a radio source by Browne (1971).An accurate redshift (z = 0.102) was determined by Nieppola et al. (2006).ON 231 was classified as an intermediate synchrotron-peaked BL Lac (IBL) object, since its low energy bump of SED peak at optical band (Nieppola et al. 2006;Abdo et al. 2010a).The very long baseline interferometry (VLBI) observations of ON 231 revealed its structure with a core and a multi-component jet elongated in Position Angle ∼ 110 • (Weistrop et al. 1985;Gabuzda et al. 1992Gabuzda et al. , 1994;;Gabuzda & Cawthorne 1996;Massaro et al. 2001).The historical optical light curves of ON 231 exhibited variability at various timescales ranging from a few hours to several years (Xie et al. 1992;Smith & Nair 1995;Gaur et al. 2012).Since the beginning of the last century, the source has shown several optical outbursts in 1940, 1953, and 1968(Tosti et al. 1998)).After that, the source started brightening slowly and reached its most luminous state in 1988 with R = 12.2 mag (Massaro et al. 1999).After the exceptional optical outburst occurred in 1998, a large bend in the jet at about 10 mas from the core was observed by the European VLBI Network plus MERLIN at 1.6 GHz and 5 GHz (Massaro et al. 2001).ON 231 was discovered as an extragalactic very-highenergy γ-ray emitter by VERITAS in 2008 March, and it entered into a second outburst state in June of the same year with the γ-ray flux about three times brighter than the previous one detected in 2008 March (Acciari et al. 2008(Acciari et al. , 2009)).Sorcia et al. (2014) reported a large rotation of PA from 78 • to 315 • (∆θ ∼ 237 • ) that coincided in time with the second γ-ray outburst.Combined with the Stokes parameters, they inferred that two optically thin synchrotron components contributed to the polarized flux.Rajput et al. (2022) studied the correlation between optical flux and polarization during ten cycles observed by the Steward Observatory.The source exhibited significant positive correlations only in two cycles.At short-term timescales, they also found two significant correlations (one is positive and the other is negative).Additionally, many studies show that ON 231 exhibits evident bluer when brighter (hereafter BWB) behaviors, just as most BL Lacs (Meng et al. 2018;Cheng et al. 2013;Tosti et al. 1998;Zhang et al. 2008).Zhang et al. (2023) found the source became bluer and then gradually stabilized when it brightened, which was briefly named the bluer-stable-when-brighter (BSWB).They inferred that the optical emission was mainly composed of two stable-color components.Generally speaking, ON 231 exhibits complex variation behaviors at various timescales, and there is still a lack of a convincing theoretical framework to elucidate the variations.
In this work, we make a comprehensive study on various variation phenomena for this target.We also propose a theoretical framework, which could consistently explain these phenomena and reveal the variation mechanism in an analytical manner.This paper is organized as follows.In Section 2, the multiwavelength data are collected and reduced.In Section 3, the multiwavelength light curves from γ-ray to radio are presented in Section 3.1.The results of the local cross-correlation analysis are presented in Section 3.2.In Section 3.3, we discuss various variation phenomena, including the γ-ray photon index (hereafter PI), X-ray PI, optical spectral index, color index (hereafter CI), and broadband SED.A theoretical framework is presented to describe these behaviors and reveal the variation mechanism.In Section 3.4, we study the correlation between PD with flux and the rotation behavior of PA on qu-plane.The above behaviors at short-term timescales are also investigated in Section 3.5.Finally, we review the complex variation behavior at various timescales and present the conclusion in Section 4.

DATA COLLECTION
In this work, we collect multiwavelength data from public data archives, including the γ-ray data from the Fermi LAT Light Curve Repository (LCR; Abdollahi et al. 2023), the X-ray data from Swift (Gehrels et al. 2004), the optical data from the All-Sky Automated Survey for Supernovae (ASAS; Shappee et al. 2014), the optical photometry, spectroscopy, and polarization data from Steward Observatory (SO; Smith et al. 2009), and the radio 15 GHz data from the Owen Valley Radio Observatory (OVRO; Richards et al. 2011).
Fermi LAT data: We collect nearly 14 years (from 2008 August 8 to 2022 June 8) γ-ray data of ON 231 from the LCR (Abdollahi et al. 2023) 1 .The LCR is a database of multi-cadence flux calibrated light curves for over 1500 sources deemed variable in the 10-year Fermi LAT point source (4FGL-DR2) catalog (Ballet et al. 2020).The photons are selected from a 12 • radius energy-independent region of interest (ROI) centered on the target.The photon energy range is from 0.1 to 100 GeV.The data analysis is performed with the standard Fermi LAT science tools version v11r5p3 using a maximum likelihood analysis (Abdo et al. 2009).The instrument response functions (IRFs), Galactic diffuse background, and isotropic γ-ray background files are P8R2 SOURCE V6, gll iem v07.fits, and iso P8R3 SOURCE V3 v1.txt, respectively.We choose one week as the time bin and use the test statistics (TS) to filter the data.Finally, we obtain 723 data points with TS >10.We also obtain the γ-ray PI fitted by the power-law model.The formula is given in dN/dE = N 0 (E/E 0 ) −Γ , here, the PI is defined as −Γ.From top to bottom, panels exhibit the light curves of γ-ray, X-ray, optical V -band, optical R-band, PD, PA, and radio 15GHz, respectively.Unfortunately, we are unable to show radio data obtained from the OVRO directly in the publication due to some data usage policy.C1∼ C10 represent the periods divided by the SO sampling interval.Different background colors also represent different periods.The purple zone is named period (A), the blue zone is named period (B), the green zone is named period (C), the red zone is named period (D), and the orange zone is named period (E).

X-ray data:
We collect the X-ray data from the XRT (Burrows et al. 2005) mounted on the Swift satellite (Gehrels et al. 2004), which operates in the energy range of 0.3 − 10 keV.All the observations are performed in Photon Counting (PC) mode.We analyze the data with the xrtpipeline task, utilizing the latest CALDB and response files provided in HEASOFT version 6.32.The source spectra are extracted from a circular area with a radius of 50 ′′ , positioned at the center of the source, whereas, the background spectra are extracted from an annulus region with an inner radius of 70 ′′ and an outer radius of 130 ′′ .We perform the X-ray spectral analysis within XSPEC (Arnaud 1996) and employ an absorbed simple power law model with a Galactic neutral hydrogen column density N H = 1.95 × 10 20 cm −2 (HI4PI Collaboration et al. 2016;Kalberla et al. 2005;Dickey & Lockman 1990).Finally, we obtain the X-ray PI and flux in different energy bands, including 0.3 − 10 keV, 0.3 − 1.5 keV, and 2.5 − 10 keV.Here, we define the soft band as 0.3−1.5 keV, and the hard band as 2.5−10 keV.The hardness ratio (HR) is calculated by HR = H/S, where H and S are the flux in the defined hard and soft bands, respectively.
Optical data: The photometry and polarization data are mainly taken from the SO (Smith et al. 2009) 2 which is an important part of the Ground-based Observational Support of the Fermi Gamma-ray Space Telescope.The duration of V -band, R-band, and polarization data observed by the 2.3 m Bok Telescope and the 1.54 m Kuiper Telescope is from 2008 October to 2018 July.All these data obtained by the SPOL CCD Imaging or Spectropolarimeter have been calibrated.We also retrieve the V -band and g-band datas from the ASAS (Shappee et al. 2014) 3 .For the above magnitude data, we convert them into flux data ulteriorly.
Additionally, we download the optical spectroscopy data from the SO and correct the Galactic interstellar extinction and reddening of all spectra by the procedure given in Cardelli et al. (1989).We chose a line free window at 5100 Å to extract the density flux F λ and translate it to F ν [mJy].To obtain the optical spectral index α o , we represent the data in the log ν -log νF ν plane and manually select six frequencies from the time integrated spectrum.These six frequencies are in the line (either emission or absorption) free regime, and the fluxes at them are of local minimum.We take the slope of the linear fitting result as α o for each spectrum.We also abandon the air absorption window at 6800 -7500 2 http://james.as.arizona.edu/∼ psmith/Fermi/ 3 https://asas-sn.osu.edu/Å for analysis.The spectral behavior for this target is exhibited in the plot of α o versus log F ν in Figure 4.
Radio 15 GHz data: The radio 15 GHz data from 2008 January 28th to 2020 January 25th are collected from the OVRO 40 m monitoring program (Richards  et al. 2011) 4 .The data point sampling of this light curve is relatively continuous.There are 600 data points in total with at least one data point per week on average.However, we couldn't display them directly in Figure 1 due to the data usage policy reasons.The data may be available on request to the OVRO 40 m collaboration.

Character of Light Curves
The multiwavelength light curves of ON 231 from MJD 53550 to 59800 are plotted in Figure 1.The light curve of γ-ray is shown on the top panel of Figure 1.Combined with previous research, it can be seen that the γ-ray flux gradually drops down to a low state (F γ <10 −7 Phs cm −2 s −1 ) after a series of very high-energy flares during 2008 (Acciari et al. 2008(Acciari et al. , 2009)).At the same time, the sparse X-ray data exhibits similar behavior, i.e., the source gradually becomes quiet after the outburst occurred in MJD 54625.On the third and fourth panels of Figure 1, the optical observation from the SO shows that the flux of R-band (F R ) varies from 1.4 mJy to 8.6 mJy while the flux of V -band (F V ) varies from 1.3 mJy to 8.3 mJy.It has a minimum brightness with F R = 1.42 mJy (F V = 1.41 mJy) around February 26th, 2009, and a maximum brightness with F R = 8.56 mJy (F V = 8.27 mJy) around June 5, 2013.The maximum flux changes are ∆F V ∼ 6.86 mJy and ∆F R ∼ 7.14 mJy, which occurs on a timescale of ∼ 1560 days (∼ 4.3 yr).According to the data from the ASAS, we notice that there is a remarkable flare with the flux of g-band up to 12.23 mJy on July 5th, 2018 which is brighter than the last outburst on February 26th, 2009.Throughout the photometric history, ON 231 has exhibited several notable outbursts in 1940, 1953, 1968, 1975, 1987, 1998, 2008(Pollock et al. 1979;;Xie et al. 1989;Massaro et al. 1992;Tosti et al. 1998;Acciari et al. 2008).From this historical record, it is likely that there is an optical periodicity of about ten years.Additionally, it's worth noting that there are many short-term flares superimposed on long-term trends.Given this, we use the trigonometric function and impulsive function to fit the long-term variation in the optical band and radio band, respectively.We plot the fitted trigonometric function curve with red solid lines on the third panel of Figure 1 and the fitted impulsive function curve with blue solid lines on the last panel of Figure 1 (not shown due to the OVRO data policy).
To study the variation behavior at various timescales, we select five periods from the time series based on the optical sampling interval and activity states.In Figure 1, period (A) from MJD 54650 to 55000 is marked with purple background, during which there is a giant flare at high luminosity state while period (C) from MJD 56275 to 56428 is marked with green background, during which there is a weak flare at low luminosity state.Period (D) from MJD 56600 to 56900 marked with red background can represent a stable flux state with a mini flare.Period (B) from MJD 55500 to 55800 marked with blue background can represent a decaying part of flare while period (E) from MJD 58050 to 58330 marked with orange background can represent a rising part of flare.Additionally, we divide the optical light curves into ten cycles from C1 to C10 based on the SO sampling interval.In summary, each period represents a different state.We will analyze the variation at various timescales in the following.

Cross-Correlation Analysis
To our best knowledge, for two separate uneven sampled time series, there are quite a few analysis methods such as the interpolated cross-correlation function (ICCF; Gaskell & Peterson 1987), the discrete Fourier transform (DFT; Scargle 1989), and the z-transformed discrete correlation function (ZDCF; Alexander 1997).Most researchers mainly use the discrete correlation function (DCF; Edelson & Krolik 1988) and the localized cross-correlation function (LCCF; Welsh 1999) to calculate the correlation coefficient between the multiwavelength light curves.Comparing the DCF and LCCF, previous studies suggest that the LCCF is more efficient than the DCF in picking up the physical signal (Max-Moerbeck et al. 2014a).Besides, considering that the LCCF exhibits a constraint on the [-1,1] range of the correlation coefficient, we finally select the LCCF to perform the correlation analysis between light curves at different bands.
The Monte Carlo (MC) procedure is applied to estimate the significance of the time lag (Shao et al. 2019).For two light curves, we will first choose the better sampling one and simulate 10000 artificial light curves by the method given in Timmer & Koenig (1995) (TK95).For ON 231, the spectral slope of the power spectral density (PSD) is computed by the program based on the actual light curves, among which β γ = 0.75, β optical = 0.21, β radio = 0.66.During the simulated procedure, we take the time interval of artificial light curves as one day.Afterward, we calculate the LCCFs between the artificial light curves and observed light curves.We obtain the 1σ, 2σ, and 3σ confidence levels corresponding to the 68.27%, 95.45%, and 99.73% chance probabilities to Notes.The negative lags indicate that the former leads the latter.τp is the lag corresponding to the highest peak of LCCF while τc is the centroid lag of the LCCF defined as τc = i τiCi/ i Ci where Ci is the correlation coefficient satisfying Ci > 0.8 LCCF (τp) (Shao et al. 2019).τ is the average of τp and τc.
evaluate the significance.The correlation peak is significant if it exceeds the 3σ level.
Based on the above methods, we perform correlation analysis on the multiwavelength light curves and the results are shown in Figure 2. To avoid spurious signals caused by the profile of the light curves, we acquire the detrend light curves of optical and radio bands.Then, we perform the cross-correlation analysis based on the detrend light curves.We find that the LCCF values usually increase and some peaks with physical significance will appear after subtracting the long-term trends.In Figure 2, a peak at about 400 days appears on panel (c), and a peak at about 250 days appears on panel (e), even though the correlation peaks remain below the 3σ level.Moreover, a negative correlation exists between PD and optical V -band at near-zero time lag.The correlation between X-ray and other bands is too scattered to be shown in Figure 2. As a whole, we fail to identify any significant (3σ) correlations between different bands.This indicates that the variation in multiwavelength may be caused by the superposition of emission from multi-components.
Based on the 2σ correlation results, we take the flux randomization and the random subset selection processes into account (Peterson et al. 1998;Larsson 2012), and calculate two kinds of time lag, namely τ p and τ c .τ p is the lag corresponding to the highest peak of LCCF while τ c is the centroid lag of the LCCF defined as (Shao et al. 2019).The time lags are listed in Table 1.In general, the γ-ray band leads the V -band while the V -band leads the radio band.The physics model to explain the time lags between the light curves of different bands is that the photons from different bands are emitted to the outside in different regions of jet (Kudryavtseva et al. 2011;Max-Moerbeck et al. 2014b).The distance between different emission regions is calculated according to the formula: D = β app c∆T /(1 + z) sin ξ, where β app is the apparent velocity in the observer frame, c is the light speed, ∆T is the time lag, z is the redshift, and ξ is the viewing angle between jet axis and observing line of sight.Therefore, results indicate that the optical emitting region may be upstream of the radio emitting region and downstream of the γ-ray emitting region.

Spectral Variability
Based on the above cross-correlation analysis in Section 3.2, we can notice that the correlation between different time series is not significant.Thus, the onezone emission model is not natural to explain the complex variability of the target.We consider the twocomponent model and try to explain the behavior of the γ-ray PI, X-ray PI, optical spectral index, CI and broadband SED in the framework of this model.Here, the model is to say that the observed flux contains two components, where F j is a flaring component and F b is a less variable background component (most probably from the postflaring blobs).For most BL Lacs, the shape of bumps in the log νF ν vs log ν plot is characterized by a rather smooth curvature extending through several frequency decades, which can be well described by a log-parabolic function.Massaro et al. (2004) provides a physical explanation in terms of statistical particle acceleration.In the jet-comoving frame, the log-parabolic model to describe the SED is given in where ν ′ p is the intrinsic peak frequency corresponding to the maximum ν ′ p F ν ′ p of the SED and b is the curvature index.Considering the Doppler enhancement, the observed flux can be described as F ν = δ m F ′ ν ′ , and the peak frequency in the observing frame is given as ν p = δν ′ p , where the δ is the Doppler factor.m takes the value 3. The flux contributed by jet in the observing frame is given in Jiang et al. ( 2024) here, ν can be taken as the different observing frequencies, including optical, γ-ray, and X-ray.

γ-ray Phonton Index and Optical Spectral Index
On the left panel of Figure 4, the γ-ray PI versus log F γ is plotted.It is obvious that PI decreases as F γ increases, showing a softer when brighter (hereafter SWB)  trend.The one-zone emission model is not suitable for producing the SWB trend, therefore, a two-component model could be a viable alternative.Similar to the redder when brighter (hereafter RWB) trend of the optical variation for FSRQs, it can be understood that the increase of the flaring component attenuates the observed spectral index if there is a stable hard background component.The origin of this background component could be studied by very-high-energy γ-ray (MAGIC Collaboration et al. 2020).Due to the fact that LAT instrument response is biased toward hard photon indexes at low fluxes, there is also a possibility that the SWB behavior is caused by this.
Additionally, we pair the optical spectral index α o with log F ν and plot them on the right panel of Figure 4.In contrast to the SWB of the γ-ray, the optical variation exhibits a harder when brighter (hereafter HWB) trend.These variation behaviors are studied by a linear fitting based on the least-squares method and the results are shown in Table 2.For period (E), it is divided into two sub-periods (E1) and (E2) when studying the behavior of the optical spectral index because it shows different behaviors in the two sub-periods.Specifically, there is a SWB trend at low luminosity (E1) and a HWB trend at high luminosity (E2).
In addition, we can figure out whether the dominant factor of the variation is extrinsic (the Doppler factor δ) or intrinsic (the intrinsic peak frequency ν ′ p ) by the slope in the α o versus log F ν plane.We consider that the optical emission is dominated by the synchrotron process in the jet and the flux can be described by Equation 3. By definition, the optical spectral index of the jet component α J ν can be derived as here, α J ν is a function of the curvature parameter b, the Doppler factor δ, and the intrinsic peak frequency ν ′ p .To further simplify, we take an empirical relation p , where n is usually taken as 0 (Jiang et al. 2024).The function of log F J ν relying on α J ν can be derived as here, C is a constant independent of ν p .In the case of the variation dominated by the Doppler factor δ, the slope is derived as K δ = (m + 1)/2b + log ν/δν ′ p −1 .The slope is derived as Since ON 231 is an intermediate synchrotron peaked blazar, the synchrotron peak frequency (δν ′ p ) is near the observed frequency ν = 10 14.78 Hz (5100 Å).We can ignore the part log(ν/δν ′ p ) ∼ 0 and make the slope formula simpler.With m = 3 and n = 0, we have K δ ≈ b/2 and K ν ′ p ≈ 2b.Based on the above scenario, if the jet component dominates the variation, the whole process will exhibit a HWB trend whether the variation is modulated by δ or ν ′ p .The only difference is that the process dominated by ν ′ p will show a steeper HWB trend compared to that dominated by δ.
The fitting slopes of α o versus log F ν for each period shown in Table 2 all exceed 1 except period (E1).We also calculate the overall fitting slope for the entire period with K obs = 0.57±0.08.In Figure 3, we exhibit the variation of the theoretical slope of α o versus log F ν with log ν p when the variation is modulated by ν ′ p and δ.We also consider different scenarios when b is taken as 0.2, 0.4, and 0.6 because the curvature parameter b usually takes different values in previous studies (Kalita et al. 2019;Anjum et al. 2020).Based on the cases in Figure 3, we can see that the theoretical slope K ν p ′ vary from 0.30 to 8.82, and K δ vary from 0.10 to 0.38.It is evident that the theoretical slope will match the observed slope better if the variation is modulated by ν ′ p rather than δ.It indicates that this variation mechanism tends to be explained by the shock-in-jet model, i.e., the relativistic electrons are accelerated by the shock in the jet, causing the intrinsic synchrotron peak frequency ν ′ p to shift to higher energy.On the right panel of Figure 4, the fitted lines on different periods seem to shift horizontally relative to each other.This can be explained by that the optical emission is not only from the flaring    2) against the average flux on the corresponding period are plotted.Different colors represent different periods, consistent with those in Figure 4.In this panel, period (E) is eliminated due to its V-shaped behavior.
component in the jet but also includes a slightly varying background.As shown in the third panel of Figure 1, short-term flares are superimposed on the long-term variable component.Such a component plays the role of the constant background when we study short-term optical spectral behaviors.Similarly, this phenomenon that the slope at long-term timescale is lower than that on each period can be the stacking effect of the event on each period.Thus, the long-term behavior may be not fundamental to reveal the variation mechanism of this target.For the V-shaped behavior in period (E), it can be interpreted that the background component dominates the emission from the jet and causes a trend of SWB before MJD 58150, and then the flaring component becomes dominant and shows the HWB trend.In Figure 5, we exhibit the slopes of α o versus log F ν (see Table 2) against the average flux on the corresponding period.Among these, period (E) is eliminated due to its V-shaped behavior.We notice that the slope decreases as the source becomes brighter.This also implies that the curvature parameter is not a fixed value in different activity states.The specific explanation of this phenomenon will be discussed in Section 3.3.2 in relation to the behavior of CI.

Color Index
Besides the spectroscopy data, we also consider the photometric data to investigate and confirm the optical variation behaviors, since the photometric data has better sampling.Many works have investigated the color variation behaviors of ON 231, which exhibit a BWB trend (Zhang et al. 2008;Cheng et al. 2013;Zhang et al. 2023).We pair the magnitudes in V -band and R-band with the time bin less than 1 day and calculate the CI V − R. According to optical data sampling of the SO, the optical light curve can be divided into ten cycles, seeing on the fourth panel of Figure 1.We select six cycles (C1, C4, C6, C7, C8, C10) for presenting the CI V − R versus V magnitude on the top panel of Figure 6.For each cycle, we calculate the slope of linear fitting and the average magnitude, and plot them on the bottom panel of Figure 6.It can be noticed that ON 231 exhibits a BWB trend at both middle-term (one cycle) and long-term timescales.The slope of CI versus magnitude becomes smaller when the source becomes brighter (see the bottom panel of Figure 6).This phenomenon is quite similar to the behavior in Figure 5.The decreasing slope with increasing flux indicates that there is a negative correlation between the peak frequency ν p and spectral curvature b, i.e., as more relativistic electrons are accelerated to the high-energy regine by the shock in the jet, the electron spectrum becomes harder and the curvature parameter b becomes smaller.This ν p -b trend can be explained in the framework of energydependent acceleration mechanisms where the acceleration probability decreases with the particle energy (Massaro et al. 2004) or through a stochastic acceleration process of relativistic particles undergoing momentum diffusion (Tramacere et al. 2011).According to the equation slope of CI versus magnitude into the slope of α V −R versus log F ν .The slope of the whole period is 0.521 which closely matches the K obs = 0.57 ± 0.08 in Section 3.3.1.The resulted slopes of the six cycles range from 0.399 to 0.889, which are smaller than the slopes in Table 2.The spectral index of the photometric data measures flux integrated on a range of wavelengths, and it is flatter than that of the spectroscopy data.If we consider that the variation is modulated by ν ′ p (K ν ′ p ∈ [0.30, 8.82]), the resulted slopes can still be understood in the framework of the shock-in-jet model.

X-ray Photon Index and Hardness Ratio
The X-ray PI and HR are effective indicators for manifesting the spectral variability of the source along with its brightness.In Figure 7, we plot the light curve of X-ray in the soft band, hard band, and HR.We select two typical periods (a), and (b), which correspond to a active state from MJD 54450 to 54650, and a quiescent state but with a flare in hard band from MJD 56300 to 56900.In Figure 8, the upper left panel exhibits the behaviors of X-ray PI versus the flux in 0.3-10 keV while the upper right panel exhibits the correlation between flux in soft band versus flux in hard band.The behaviors of X-ray PI with the soft band flux and the hard band flux are plotted on the two bottom panels.We notice that the variation behavior on each panel has a turning point and two branches.The variation in period (a) and in period (b) coincide on different branches.Combined with the bottom panel of Figure 7, the HR <1 in period (a) indicates that the variation is more significant at the soft band.Conversely, the HR >1 in period (b) indicates the variation is more significant at the hard band.The different behaviors during period (a) and period (b) demonstrate that the X-ray radiation in different energy bands could originate from different processes.
Based on the above analysis, we consider that the hard or soft X-ray emission includes both the synchrotron and Synchrotron Self-Compton (hereafter SSC) components, i.e., F obs ν h,s = F syn ν h,s + F ssc ν h,s .Previous studies show the fact that the log-parabolic model has advantages in modeling the spectral distributions of the synchrotron and the SSC components (Massaro et al. 2006).Therefore, we adopt two log-parabolic functions to represent the synchrotron component and the SSC component, respectively.The specific formulas are given as follows here, the ν p1 , F ′ In the framework of the shock-in-jet model, we consider that the variation process is dominated by the variation of the intrinsic peak frequency ν ′ p .In theory, the X-ray PI is calculated by α X = (log F ν h − log F νs )/(log ν h − log ν s ).The hard and soft frequencies corresponding to 5 keV and 0.67 keV (roughly the geometric median of the hard and soft band) are log ν h = 18.08 and log ν s = 17.21, respectively.To model the behaviors of PI versus log F X and log F H versus log F S , we first set the parameters as m = 3, δ = 10, γ = 1300.Then, we manually vary the peak frequency ν ′ p [Hz] ranging from 12.0 to 13.5, and adjust the other to match the data as closely as possible.The simulation results, marked by the red asterisked lines, are plotted in Figure 8.The corresponding free parameter values are shown in the caption of Figure 8. Considering that the free parameters are coupled, the fitted results still have a certain degree of degeneracy.It is clear that the theoretical model can well reproduce the above variation on all panels of Figure 8,  especially for the turning point and two branches.Taking panel (a) for example, while the synchrotron component dominates the X-ray emisson, the target exhibit a HWB trend with a smaller slope.While the SSC component dominates the X-ray emisson, the target exhibit a HWB trend with a higher slope.If the contribution of the synchrotron component and SSC component about equal, the target could exhibit a SWB trend.In summary, we tend to conclude that the X-ray emission originates from the superposition of the synchrotron component and the SSC component, and the variation is modulated by the shift of the intrinsic peak-frequency ν ′ p caused by the shock accelerating in the jet.

Broadband Spectral Energy Distribution
To understand the diverse behaviour of the source during different periods, we select two epochs from period (A) (high luminosity state), and the other two epochs from period (D) (low luminosity state).Epochs (A1), (A2), (D1), and (D2) are around MJD 54863, 54884, 56716, and 56740, respectively.The time ranges for these epochs is 1 ∼ 3 days, and a longer time range is used in rare cases due to the lack of data.In Figure 9, we construct SEDs and model them using the theoretical framework mentioned in Section 3.3.3.Specifically, we use two log-parabolic functions to represent the synchrotron component and SSC component, respectively.The specific formulas are given in Equation 6and Equation 7. The fitting results of the broadband SEDs in different epochs are listed in Table 3.The SEDs during these four epochs exhibit different behaviors.In epoch (A1), the optical band is around the synchrotron peak frequency ν p , and the X-ray emission includes both the synchrotron and SSC components.In epoch (A2), the optical band falls on the left side of the synchrotron peak frequency ν p , and the X-ray emission is from the synchrotron tail.In epoch (D1), the optical band is located on the right side of the synchrotron peak frequency ν p , and the X-ray emission includes both the synchrotron and SSC components but the SSC component dominates.In epoch (D2), the optical band falls on the right side of the synchrotron peak frequency ν p , and the X-ray emission is from SSC component.Combined with the fitting results in Table 3, we find that the variation of this target is dominanted by the shift of the intrinsic peak-frequency ν ′ p .As ν ′ p increases, the optical band varies from being on the right side of synchrotron peak frequency ν p to around ν p , and then to the left side of ν p .Meanwhile, the dominant component contributing to the X-ray emission varies from the SSC component to the synchrotron component.In most cases, the two components are in competition with each other and superimposed to compose the X-ray emission.In conclusion, the broadband SEDs visually demonstrate that the variation is dominated by ν ′ p , and X-ray emission is the superposition of both the synchrotron and SSC components.

Optical Polarization and Optical Flux
The long-term monitoring of optical polarization has shown that polarization behavior is complex, and it serves as a significant window to figure out the jet structure and radiation mechanism in BL Lacs.Generally, the shock-in-jet scenario is favored for the correlation between optical flux and PD, whereas the anticorrelation can be explained by the multi-component model (Otero-Santos et al. 2022).
We plot logPD versus log F ν (V -band flux) in Figure 10.The linear fitting results were listed in Table 4. Four segments are selected and represented by different colors in Figure 10, and each segment is divided into two subperiods.Noteworthy, three segments (marked by the blue, green, and red colors) show the V-shaped behav- ior, one segment (marked by the yellow color) shows the inverted V-shaped behavior in the correlation between PD and flux, i.e., for each segment, one sub-period exhibits the positive correlation while the other sub-period exhibits the negative correlation.Marscher (2014) proposed a Turbulent Extreme Multi-Zone (TEMZ) model to explain the variability of the flux and polarization.
The model divides the emitting region into a large number of emission zones ("cells") with different magnetic field orientations, and the superposition of synchrotron radiation from these cells can lead to a decreasing polarization with increasing flux.In this work, we consider the two-component model.Both components are of jet.One is the polarized background component, representing the post-flaring blobs, and the other component represents a newly appearing flare with different magnetic field orientations.With the increasing luminosity caused by the flaring component, the polarization from the background component is diluted, and the observed PD decreases.When the flaring component becomes dominant, it exhibits a positive correlation between PD and flux, which could be interpreted by the shock-in-jet model.In Figure 10, the three V-shaped behaviors have a translation relation, and it seems that the lower limit of the logPD of the V-shaped behaviors for different segments decreases with increasing flux.A similar translation relation is also found in the plot of α o versus log F ν in Figure 4.Both the PD and spectral behaviors could be explained if there exists a long-term slightly varying background component.For PD, the V-shaped behavior shifts towards the high flux segment as the background component brightens.For the spectral behaviors, the Vshaped behavior (see Figure 4, the green and light blue points, from MJD 57930 to 58260) could be understood if the background component has a bluer spectrum than the new flare.Thus, The translation phenomenon can be understood well by the two-component model, and the variation is dominated by the shock in the jet.
F l u x < 2 2 < F l u x < 2 .52 .5 < F l u x < 3 3 < F l u x < 4 4 < F l u x 0 .0 2 .0 4 .0 6 .0 8 .0 Figure 12.Each dot in the top three panels corresponds to one short period in Table 5. Panel (a), (b), and (c) exhibit the variation of the PD-slope vs. the average flux, the CI-slope vs. the average flux, and the PD-slope vs. the CI-slope, respectively.The dots with blue color represent the period with the CW rotation behavior while red corresponds to the ACW rotation behavior.The bottom three panels exhibit the statistical results of Table 5.The bar charts exhibit the number of periods with different behaviors while the line charts exhibit the proportion of periods with different behaviors.Panel (d), (e), and (f ) exhibit the variation of the correlation of logPD vs. log Fν , the CI behavior, and the rotation behavior on qu-plane with increasing flux, respectively.

Rotation Behavior of Polarization Angle
The PA, χ, is defined within a 180 • -interval, thus the rotation behavior of PA can't be directly presented because of the nπ-ambiguity: χ = χ±nπ, n ∈ N (Marscher et al. 2008;Abdo et al. 2010b;Kiehlmann et al. 2016).To escape the nπ-ambiguity problem, we study the variation behavior on the qu-plane.Here, q = Q/I, u = U/I, Q, U , and I are the Stokes parameters.We plot the q versus u for different periods in Figure 11.The left panel shows the one-year binned qu data, where the circle size represents the flux value and the colors represent the observation date.The average of q and u of all data points (marked by a blue asterisk) deviates significantly from the origin (q = u = 0).This translates a constant polarization degree P ave = 10.7% ± 5.1% and polarization angle Θ ave = 69.4 In the left panel of Figure 11, we also find an interesting phenomenon that the rotation of PA is first clockwise at higher luminosity and then anticlockwise or oscillation at lower luminosity.We think that the rotation and oscillation behaviors of PA originate from the flaring activities.Additionally, we select a middle-term timescale period with a typical anticlockwise (hereafter ACW) rotation behavior and plot the qu data on the middle panel of Figure 11.Interestingly, the ACW rotation of PA coincides in time with multiwavelength flares.We plot the simultaneous light curves of γ-ray, X-ray, and optical V -band on the right panel of Figure 11, among which the γ-ray peak at MJD 56377.5, the X-ray peak at MJD 56380.9, the optical V -band peak at MJD 56388.4.There appear to be time lags between multiwavelength flares, i.e., the γ-ray band leads the X-ray band while the X-ray band leads the V -band.Generally, the above phenomenon suggests that the polarimetric variation is a joint contribution of a stable polarization background component and a variable flaring component.

Statistics of Behavior at Short-term Timescale
To investigate the mechanism of the short-term variation, we divide the light curves into many short-term periods of about one week.The division principle is mainly based on the optical data sampling of the SO.Some periods include a longer time range due to the lack of data points or the guarantee for the integrity of variation.We calculate the parameters (Slope, Pearson's r, P-value) of linear fitting logPD versus log F ν and V − R versus V during every short period.For simplicity, the slope of logPD versus log F ν is named as PD-slope, and the slope of V − R versus V is named as CI-slope.These parameters together with the rotation behavior of PA, and the V -band flux are listed in Table 5.In Figure 12, panel (a) exhibits that about half of the periods show a positive correlation of logPD versus log F ν while the other half showed a negative correlation.It seems that the correlation appears to be independent of the rotation behavior of PA.On panel (b), most of the periods show the BWB trend, among which the periods with high CI-slopes show the ACW rotation behavior.Additionally, we also notice that the CI-slope becomes flatter as the flux increases, and this phenomenon also seems to be present at middle-term timescales in Figure 6.It can be explained by the shock-in-jet model with a varying curvature parameter, as discussed in Section 3.3.2.For panel (c), there seems to be a weak positive correlation between PD-slope and CI-slope if those periods with high PD-slope or CI-slope are excluded.
On the three bottom panels of Figure 12, we sort out the number and proportion of the periods with different behaviors.The bar charts represent the number of different behaviors while the line charts represent the per-centage of different behaviors.As readily seen, the periods with the positive correlation of logPD versus log F ν and the CW rotation behavior gradually dominate with the increasing flux.However, there are always more periods with BWB trends than that with RWB trends no matter whether in the active or quiet state.According to the statistical results at short-term timescales, the twocomponent model could also be valid in explaining the transitions of different behaviors with increasing flux.

CONCLUSIONS
In this work, we gather the multiwavelength light curves of ON 231, including γ-ray, X-ray, optical band, optical PD, optical PA, and 15 GHz radio.We perform the correlation analysis and calculate time lags among them.We discern the variation mechanism by analyzing the spectral behaviors at γ-ray, X-ray, and optical band.To further constrain the variation mechanism, we investigate the correlation of the optical flux versus PD and the rotation behavior of PA.In addition, we also investigate the variation of the above behaviors at short-term timescales.The principal findings of this work are summarized as follows: • According to the cross-correlation analysis results, we find that no correlations are beyond the 3σ significance.However, there seem to be some 2σ correlations among the light curves of γ-ray, V -band, and radio with large time lags.Specifically, γ-ray leads V -band 165.4 +56.8  −123.3 days, V -band leads the radio 246.6 +26.4   −38.4   days, and γ-ray leads the radio 387.8 +13.7 −38.7 days.Additionally, there is a negative correlation between PD and optical V -band at the near-zero time lag.The above phenomenon implies that the one-zone emission model is not favored.
• The behaviors of the multi-band spectral index are investigated.For the γ-ray PI, it shows a SWB trend.The variation mechanism could be understood by the two-component model, in which a stable hard background component reserves further study.For the optical spectral index, it usually exhibits a HWB trend.But, in one period, it shows a trend transition from SWB to HWB.Additionally, we find that the slope of α o versus log F ν decreases as the flux increases.Combined with the behavior of CI, it shows a BWB trend, and the CI-slope also becomes smaller when the source becomes brighter.We present a theoretical frame to quantitatively explain these spectral behaviors.We find that both the geometric (δ) and intrinsic (ν ′ p ) variations can produce a HWB trend.Considering the slope of α o versus log F ν , the shock-in-jet model is more promising than the helical jet model as the dominant mechanism of variation.The decreasing slope indicates the energy-dependent acceleration processes of the radiative particles.The most complex variation is in the X-ray spectrum, it shows a trend transition from HWB to SWB to HWB.We consider that the X-ray spectrum is the superposition of both the synchrotron and SSC components described by two log-parabolic functions.In the case of the variation modulated by ν ′ p , the theoretical framework could well describe the spectral behavior of X-ray.Additionally, We constructed SEDs in different epochs, and the variation also mainly modulated by ν ′ p .It again confirms that the variation mechanism tends to be explained by the shock-in-jet model.
• We find that the correlation between the optical PD and flux exhibits the typical V-shaped behavior.We also notice that there is a rotation behavior of PA that coincides in time with the multiwavelength flares.In addition, the statistics results of short-term behavior exhibit a similar trend transition when the source becomes brighter.These polarization phenomena could also be explained by the two-component model in the shock-in-jet scenario, in which one is a highly polarized background component, and the other is a significantly varied flaring component.
Figure1.From top to bottom, panels exhibit the light curves of γ-ray, X-ray, optical V -band, optical R-band, PD, PA, and radio 15GHz, respectively.Unfortunately, we are unable to show radio data obtained from the OVRO directly in the publication due to some data usage policy.C1∼ C10 represent the periods divided by the SO sampling interval.Different background colors also represent different periods.The purple zone is named period (A), the blue zone is named period (B), the green zone is named period (C), the red zone is named period (D), and the orange zone is named period (E).

Figure 2 .
Figure2.The LCCF analysis results among light curves of different wavebands are plotted.The significance levels of 1σ, 2σ, and 3σ are indicated by the orange, green, and red lines, respectively.The blue dots in some panels represent the LCCF results based on the corrected radio and optical light curves.

Figure 3 .
Figure 3.The variation of the theoretical slope of αo versus log Fν with log νp in different cases are plotted.The red color represents cases where ν ′ p dominates the variation and the blue color represents cases where δ dominates the variation.Solid, dashed, and dotted lines correspond to the cases b = 0.2, b = 0.4, and b = 0.6 respectively.

Figure 4 .
Figure 4.The γ-ray PI vs. log Fγ is plotted on the left panel while the optical spectral index αo vs. log Fν is plotted on the right panel.Different colors correspond to different periods.

Figure 5 .
Figure 5.The slopes of αo vs. log Fν (see Table2) against the average flux on the corresponding period are plotted.Different colors represent different periods, consistent with those in Figure4.In this panel, period (E) is eliminated due to its V-shaped behavior.

Figure 6 .
Figure6.The CI V − R vs. V magnitude is plotted on the top panel and the slope vs. average magnitude on each period is plotted on the bottom panel.Different colors correspond to different periods.

Figure 7 .
Figure 7. From top to bottom, panels exhibit the light curves of FX in the energy ranges of 0.3 − 1.5 keV, 2.5 − 10.0 keV, and the HR.Different colors represent different periods.The dashed line in the bottom panel represents HR=1.

Figure 8 .
Figure 8.The upper left panel exhibits the behaviors of X-ray PI with flux in 0.3 − 10 keV.The upper right panel exhibits the correlation between log νH FH vs. log νSFS.The bottom left panel exhibits the behaviors of X-ray PI with flux in 0.3 − 1.5 keV while the bottom right panel exhibits the behaviors of X-ray PI with flux in 2.5 − 10 keV.Different colors represent different periods.The red asterisked lines represent the numerical results of the model.On all panels, we set the parameters as b1 = 0.08, b2 = 0.55, F ′ ν ′ p 1 = 10 −28 mJy, F ′ ν ′ p 2 = 10 −32.8 mJy.log ν ′ p [Hz] ranges from 12.0 to 13.5.
represent the observed peak frequency, the intrinsic flux at the peak frequency, and the curvature parameter of the synchrotron component, respectively.The ν p2 , F ′ ν ′ p 2 , b 2 represent those of the SSC component.The observed peak frequency of the synchrotron component and SSC component can be related by ν p2 /ν p1 = 4γ 2 , where γ is the Lorentz factor of electrons.

Figure 9 .
Figure 9.The broadband SED fitted by two log-parabolic models at four different epochs.The dashed line represents the synchrotron component and the dotted line represents the SSC component.

Figure 10 .
Figure10.Different colors correspond to different segments.Each segment is divided into two sub-periods.The sub-periods represented by the solid dots exhibit a positive correlation while the sub-periods represented by the hollow dots exhibit a negative correlation.

Figure 11 .
Figure 11.The rotation behavior on qu-plane based on one-year binned qu data is shown on the left panel.The circle size represents the averaged flux for the corresponding period.The colors of the circles represent the observation date, where values are listed in the color bar.The middle panel exhibits a typical ACW rotation behavior while the corresponding multiwavelength flares are plotted on the right panel.The blue asterisk represents the center of each rotation behavior.
s r P-value Correlation Slope Pearson's r P-value Tendency Min Max Average 5

Table 1 .
The results of time lag analysis

Table 2 .
Linear fitting results of PI versus log Fγ and αo versus log Fν

Table 3 .
The fitting results of the broadband SED in different epochs

Table 4 .
Linear fitting results of logPD versus log Fν

Table 5 .
The statistical results of the behaviors at short-term timescale.Here, PC or NC denotes the positive or negative correlation, respectively.