Correlation analysis between several observables of OJ 287 radio jet

We collected the archival data of blazar OJ~287 from heterogeneous very long baseline interferometry (VLBI) monitoring programs at 2.3 GHz, 8.6 GHz, 15 GHz and 43 GHz. The data reduction and observable extraction of those multi-band multi-epoch observations are batch-processed consistently with our automated pipeline. We present the multivariate correlation analysis on the observables at each band. We employ the cross-correlation function to search the correlations and the Monte Carlo (MC) technique to verify the certainty of correlations. Several correlations are found. The foremost findings are the correlations between the core flux density and the jet position angles on different scales, which validated the plausible predictions of the jet with precession characteristics. Meanwhile, there is a variation in the offset between the core EVPA and the inner-jet position angle over time at 15~GHz and 43~GHz.


INTRODUCTION
Blazars are a particular class of active galactic nucleus (AGN) with their jet oriented at a very small angle with respect to the line of sight (Urry & Padovani 1995).The relativistic jet propagates from the vicinity of the actively accreting supermassive black hole (SMBH) in the center.The blazar subclass consists of BL Lacertae objects (BL Lacs) and flat-spectrum radio quasars (FS-RQs) identified by their optical spectra.Blazars show flux variability at almost all wavelengths of the electromagnetic (EM) spectrum and significant polarization variability in both optical and radio bands.
The VLBI is a powerful technique to increase the radio telescope's resolution by forming a synthetic aperture with long baselines.More and more high-resolution VLBI observations have revealed the swinging innermost jet position angle and complex innermost jet structure in the plane of the sky (Lister et al. 2013;Jorstad et al. 2017).Lister et al. (2021) analyzed a large data set and concluded that most of the jets show variations between Corresponding author: M. Zhang mgnahz@gmail.com10 • and 50 • in their inner-jet position angle over time.The essential reason for this jet bending phenomenon in the innermost region may be the wobbling flow instability near the jet base, which includes Lense-Thirring precession of the black hole and the accretion disk (Caproni & Abraham 2004), jet precession caused by the binary black hole system (Britzen et al. 2018;Dey et al. 2021) or the warping of the accretion disk (Lai 2003), magnetohydrodynamic (MHD) flow instabilities (Matveyenko & Seleznev 2015).
Long-term observed flux density variations can be comodulated by multiple mechanisms induced by intrinsic and extrinsic factors, which are often difficult to distinguish (Böttcher 2019).Extrinsic factors mainly include the geometrical effects caused by non-radial jet motion accompanied by the viewing angle variations and consequent Doppler factors variations.Meanwhile, the nonradial jet motion essentially forms the intrinsic helical or curved structure, and the projection on the sky plane is S-shaped or some complex morphology.If the above geometric effect is the main reason for the variability in radio emission and jet morphology, a correlation between the flux density and the jet position angle could be predicted in that case.Precession is a special and impor-tant mechanism that causes geometric effects.So, the source OJ 287 with long-term multi-band VLBI monitoring and well-studied precession characteristics is selected for correlation analysis.
Through the time-varying properties of the jet direction caused by the jet precession, we can study whether the rotation of the core EVPA is consistent with that of the jet position angle at the core scale.That is, whether both maintain relatively constant directions during longterm monitoring, as the theory studies suggest (Laing 1980;Lyutikov et al. 2005).It should be mentioned that the results of various VLBI studies on radio core EVPA and local position angle alignment in AGN jets are observationally controversial (Agudo et al. 2018;Hodge et al. 2018).
OJ 287 is a low synchrotron peaked (LSP) BL Lac object which is highly active at all wavelengths and famous for its repetitive optical flares with roughly a 12year period (Sillanpaa et al. 1988) and variable jet morphology (Gómez et al. 2022).As seen in many other BL Lac objects, OJ 287 has a one-sided core-jet structure.Up to now, the jet of OJ 287 has been extensively investigated at several wavelengths and with different long-term VLBI monitoring programs, like the MOJAVE1 , the 2 cm survey (Kellermann et al. 1998) and the VLBA-BU-BLAZAR program2 (Jorstad et al. 2017).Tateyama & Kingham (2004) reported the position angle of the more collimated unresolved jet rotated clockwise by ˜30 • from 1994 to 2002 due to the ballistic precession of the jet.Agudo et al. (2012) found a sharp innermost jet position angle swing during 2004 and 2006.In recent years, Cohen (2017) suggested that the jet is rotating on a 30-year period by analyzing the jet ridge lines.And Britzen et al. (2018) found that the OJ 287 jet is precessing on a time-scale of ˜22 years.The radio emission is highly polarized, the time-domain variations of the degree of the linear polarization and the EVPA have been extensively studied in radio bands (Roberts et al. 1987;Cohen et al. 2018).
The structure of this paper is as follows: In Section 2, we describe the observations, the data reduction techniques, and the methods to derive observables.In section 3, we describe the method used in the correlation analysis in detail.Section 4 shows the results of correlations between observables.Section 5 includes some discussions of current correlation analysis results and the significance and merits of correlation research.The final section gives a summary of our findings.All position angles are expressed in degrees from north to south, between 0 • to −180 • in the clockwise and 0 • to 180 • in the counterclockwise direction.The RDV program utilizes simultaneous 2.3 GHz/8.6 GHz global-VLBI observations to monitor the radio source structures with geodetic antennas around the world.We use the data observed from July 1994 to December 2003 until the source was dropped from the astrometric monitoring due to the structural complexity.The RDV data were acquired through collaboration with the Bordeaux VLBI Image Database (BVID).For this source, the database maintained by the MOJAVE team comprises 127 observed epochs over a span of 25 years from April 1995 to June 2020.The MOJAVE database consists of two parts: the 2 cm Survey and the MOJAVE program initiated in 2002 as the successor to the former 3 .MOJAVE program observations were carried out in dual polarization mode using frequencies centered at 15.4 GHz, with a bandwidth of 32 MHz.We use the self-calibrated visibility data from the streamed release of the MOJAVE database in our analyzing pipeline.In addition, we also use the 43 GHz data from the VLBA-BU-BLAZAR program, which includes data from June 2007 and onwards.The VLBA-BU-BLAZAR observations were made in continuum mode, with both left and right circular polarizations at a central frequency of 43.1 GHz.It also provided self-calibrated visibility data.In this paper, we analyze their released data onwards until January 2020.

Data reduction
The data are reduced using the Astronomical Image Processing System (AIPS) package and our SAND pipeline (Zhang et al. 2018) was used to automate the procedures of mapping and model-fitting, as well as the observable extraction and post-processing.

Imaging
The CLEANed images are produced with the AIPS task IMAGR using the BGC CLEAN algorithm (Clark 1980).We choose to use uniform weighting and a proportional-to-beam cell size of 0.5, 0.15, 0.08 and 0.02 mas pixels −1 from the lowest frequency (2.3 GHz) to the highest frequency (43 GHz), respectively.We use restoring beam smaller than the nominal CLEAN beam obtained by fitting a Gaussian to the central lobe of the dirty beam to create slightly super-resolved images.In most cases, our 15 GHz and 43 GHz observations have full Stokes parameters successfully correlated, thence the pipeline will automatically produce polarization maps.Fig. 1 displays a sample of quad-band CLEANed images contours with the Gaussian modelfitting components superimposed at each frequency.In the same band, all images have been convolved with the same restoring beam, i.e., the FWHM dimensions of circular restoring beam size are 3, 0.7, 0.4, 0.2 mas from the lowest frequency to the highest frequency, respectively.

Model fitting
Our pipeline utilizes the task SAD to extract sources in the CLEANed images, and fit Gaussian components to those.The task can search for potential sources with bright peaks above a certain flux threshold.We use the CPARM parameters to set the threshold level in descending order of powers of the root mean square (RMS) noise.The background noise level is obtained by examining pixels well away from the sources in each Stokes I image with the task IMSTAT.The components are fitted within rectangular islands, and it may be possible to get bad solutions, so the spurious components are rejected by setting DPARM parameters.

Properties
Usually, the most upstream visible component is brighter than downstream components, which can be assigned as the 'core'.However, there can be confused situations when the image configuration is complex as in the high-frequency bands.In this paper, we usually defined the brightest and compact emission feature as the VLBI core, but not always.For special cases of complex situations in high frequency, we carefully identify the core by taking flux density, degree of compact, and position of components in adjacent epochs into consideration.

Core measurements
In Table 1, we list the general properties of the fitted core component.A previous study by Liu et al. (2012) used the position angle of the major axis of the fitted core component to determine and define the inner-jet position angle.In our research, we also follow this definition.Obviously, if the ratio between the major axis and the minor axis of the fitted core component is close to 1, the major axis position angle can be an arbitrary value.In this case, the outlier will be substituted with the mean value of two adjacent epochs' inner-jet position angles.As we know, the non-directional major axis position angle is centrosymmetric, so any position angle is equivalent to itself with an nπ flip.To ensure that the correct inner-jet position angle is chosen, we select the multi-epoch position angle based on the jet direction matching criteria.

Core polarization properties
If full Stokes parameters are correlated for a data record, SAND will generate the Q and U images.And then, the Stokes Q and U images are combined using the AIPS task COMB to create the polarized intensity and position angle maps.The EVPA is calculated as EVPA = (1/2)tan −1 (U/Q).We averaged the EVPAs of all pixels within an area of 3 × 3 pixels centered on the core of the Stokes I image in the position angle map as the nominal core EVPA.We used the following equations to estimate the errors in EVPA: where p = Q 2 + U 2 , Q and U are the Stokes parameters in the given pixels, and σ Q and σ U uncertainties in the Stokes Q, and U data.σ Q and σ U should include contributions from the rms error, D-term error and CLEAN error, defined as follows: (2) where σ rms , σ Dterm , and σ CLEAN denote rms noise, Dterm error, and CLEAN errors, respectively.According to Hovatta et al. (2012)'s simulations that considering the instrument polarization, the additional error σ Dterm is defined as (3) where N ant is the number of antennas, N IF is the number of IFs, N scan is the number of scans with independent parallactic angles, and N peak is the peak total intensity map.We used AIPS task PRTAN and task LISTR  with optype='scan' to read N ant , N IF , N peak value from the actual recorded data epoch by epoch for σ Dterm calculation.We adopted the error from imperfect EVPA calibration of 3 • at 15 GHz (Hovatta et al. 2012) and of 7 • at 43 GHz (Kravchenko et al. 2020;Jorstad et al. 2005), and added this error to Equation (1)(σ EVPA ) in quadrature.We also averaged the error over a small region of 9 pixels as the core EVPA error.The core EVPAs measured and its error parameters are listed in Table 1.
Polarized waves are affected by Faraday rotation when propagating through non-relativistic plasma within or external to the source.In the case of external Faraday rotation, the observed EVPA(χ obs ) and wavelength squared appear a linear dependence and follow Equation 4 where χ 0 is the intrinsic EVPA and RM is the rotation measure.
The rotating medium located outside the jet mainly includes a sheath surrounding the jet or the broad/narrow-line regions (BLRs/NLRs), even intergalactic and galactic plasma.It has been reported that galactic RM of OJ 287 is 30 rad m −2 (Rudnick & Jones 1983), which rotates the 15 GHz EVPA and 43 GHz EVPA far less than 1 • .For this source, Hovatta et al. (2012) used VLBA observations that carried out four frequencies between 8 GHz and 15 GHz in 2006 April and gave the median RM over the core is −307.9rad m −2 .The rotation value calculated by RM = −307.9rad m −2 is about 7 • at 15 GHz.Assuming that the RM is also available at 43 GHz, the rotation value will be smaller there.Previous studies have shown that the opaque core of OJ 287 not show signatures of internal Faraday rotation.For these reasons, we do not need to correct the Faraday rotation of the external to the source.A significant focus question in analyzing polarization data is the nπ ambiguity (Marscher et al. 2008;Larionov et al. 2008;Abdo et al. 2010) of the observed EVPA.Several different methods (Kiehlmann et al. 2016) were proposed to work out this issue.These methods are based on the assumption of minimal variation between adjacent data points or between the current and a couple of previous data points.For our sparsely-sampled data, we choose the final core EVPA time series, which fall within the narrowest EVPA range, in the set of data points and their ±nπ flips.Admittedly, it is difficult to capture the rapid rotation of the EVPA (≥180 • ) in the interval of two adjacent epochs due to the sampling rate limitation.At radio wavelengths, the rapid rotation phenomenon is rarer than that at optical wavelengths (Aller et al. 2003;

Jet measurement
In Table 2, we list the modeled properties of the recognized jet component.Various methods have been utilized to identify the jet position angle, such as using the position angle of the brightest jet component with respect to the core or taking a flux-density-weighted position angle average of all fitted jet components from the core (Lister et al. 2013).We follow the prescription of Valtonen & Wiik (2012) and Moór et al. (2011) to fit a regression line that is forced to go through the location of the core component.Slightly different from previous methods, we use the orthogonal distance regression (ODR) algorithm to get the best linear fitting rather than the ordinary least-squares (OLS) fit (Moór et al. 2011) or the ordinary least-squares bisector (OLSbisector) regression (Valtonen & Wiik 2012).This method can be used to take into account the errors in both relative coordinates.The data points are weighted by the inverse square of their positional uncertainty.The uncertainty of the position of the fitted Gaussian components in coordinates along their major and minor axes are estimated by the following formula: (5) µ contains information in two directions: µ(x 0 ) and µ(y 0 ).µ(x 0 ) and µ(y 0 ) are the rms lengths of the major and minor axes of the error ellipse composed of position uncertainty.SNR is the signal-to-noise of the component, which is determined by dividing the peak flux density of the fitted Gaussian component by the off-source noise.Θ image size is the deconvolved size (the full-width at half-maximum, FWHM) of major or minor axis which deconvolves the Clean beam Θ beam from the fitted component size.The error ellipse major-axis position angle φ is measured from north to east, then the rightascension and declination errors, i.e., µ(α) and µ(δ), are calculated by (Condon 1997) The ODR result outputs the slope of the best linear regression line (k) and the standard error (σ k ); we calculate the jet position angle and its error with k and σ k .Fig. 2 illustrates our fitting method for three contiguous epochs at 15 GHz and 43 GHz.If only one jet component is fitted within a reasonable distance from the core, the position angle of this component is assumed as the jet position angle.All calculated jet position angles are listed in column (6) in Table 2.

Time-domain observables
We extracted some observables, namely the core flux density, the inner-jet position angle, the jet position angle, and the core EVPA, their evolutions over time are shown in Fig. 3.We follow Homan et al. (2002)'s work and use the dominant calibration error, estimated to be ˜5% of the core flux density, as the core flux density error.The inner-jet position angle error is estimated by AIPS using the formulae given from (Fomalont 1999).The error of the jet position angle and the core EVPA is as described above.

CORRELATION ANALYSIS METHOD
The cross-correlation function (CCF) is commonly employed in the study of AGN, including probing the structure of the broad emission-line region by reverberation mapping (Peterson et al. 2004), studying the continuum emission mechanism by correlating multi-band light curves (Max-Moerbeck et al. 2014), and finding correlations between the radiation variability and other AGN properties (Rani et al. 2014).We also employed the CCF to seek correlations among variations of the flux density, jet position angle, and EVPA.
When calculating the cross-correlation, each point in one time series of observables must be paired with a point in the other time series, so the data points should be regularly spaced, otherwise the pairing is only possible when the two-time series are aligned.The datasets we collected were unevenly sampled due to observation schedules, project time allocations, and other constraints such as equipment failures and severe weather.First, we expanded the amount of data three times by linear interpolation to make the data set even at each band.The time interval for our interpolation seems to be a reasonable choice since it is shorter than the majority of the real data intervals, but not so short that a large number of artificial data points interfere with the analysis.We also compared another interpolation scheme (e.g., cubic spline interpolation) but found that this interpolation sometimes overfits, not always reliable.We divided the entire time series into two segments for correlation analysis at the time of the large changes of the inner-jet position angle of 15 GHz or the jet position angle of 43 GHz.Moreover, three data points after the 15 GHz inner-jet position angle change were removed, as well as several data points before the 43 GHz jet position angle change.Our purpose in doing this is to make the time series of each segment have quasi-stationary properties as much as possible.Then we calculated the normalized CCF between each pair of observables, with the maximum and secondary peaks obtained from the correlation coefficient sequence.We analyze the significance of cross-correlation results by Monte Carlo (MC) trials.At the same time, we should not ignore the fact that since the overlap between the observables is smaller for large values of time delay, we deem these peaks relatively less credible.
We constructed 10,000 entirely simulated time series for each observable that span the approximate range of variation measured in the real data but have their own variations independent of each other in time domain.We repeated the step of the previous standard cross-correlation analysis for 10,000 simulated observable's time series pairs and accumulated the acquired cross-correlation coefficients for each time delay.And for each time delay, the significance level of the real data cross-correlation coefficient is estimated based on the distribution of the simulated cross-correlation coefficients.
We also employed MC sampling to estimate the certainty of these correlations affected by all errors on the observables.The correlation coefficient errors were de- termined from 10,000 MC trials by adding sampled normally distributed errors to the expected value of the observables and then performing correlation calculations for each trial.We subtracted the correlation coefficient calculated from the expected value of the observables (ρ expt.,i ) from the MC trial's correlation coefficient (ρ M C,i,n ) to get a series of correlation coefficient error values (∆ ρ i,n ), i.e., ∆ ρ i,n = ρ M C,i,n − ρ expt.,i, i refers to the index of the time delay slice and n denotes the nth MC trial.We verified ∆ ρ i of each time delay slice following the normal distribution at a significance level of 5% by the Shapiro Wilk tests (SHAPIRO & WILK 1965).When each ρ M C,i,n and the ρ expt.,i are very close, i.e., the mean and the standard deviation of the ∆ρ i distribution close to 0, it can be considered that the error of the observables does not affect the final correlation result.

RESULTS OF CORRELATIONS BETWEEN OBSERVABLES
The correlation results which have physical meaning between each pair of observables are shown in Fig. 4 and Table 3.
It is very likely to make wrong physical interpretations based solely on the maximum peak of CCF without considering its significance.If and only if there is no secondary peak near zero time delay falling within plus or minus 0.8 times the maximum peak's value of the CCF and this peak with a 95% significance level, we claim that the maximum peak is significant and the two observables are correlated.If there is a secondary peak which also has a 95% significance level, then it indicates that the two observables are more likely to correlate at that peak.
The estimation of the significance of the crosscorrelation achieved with the MC method is shown in Fig. 5.Meanwhile, the MC trials involving the effect of the observational error on correlation indicate that the maximum peak of the n-th trial's correlation coefficient sequence (ρ M C,n ) and the corresponding time delay have slight changes, but the error of the observables does not affect the correlation results for each pair of observables (Fig. 6).
The core flux density and the inner-jet position angle have a significant negative correlation with each other at each band.In segment I, at 8.6 GHz and 15 GHz, the core flux density appears to be anti-correlated with the jet position angle; and at 2.3 GHz, it shows a positive correlation between them (Fig. 4 (b)).In segment II, the core flux density appears to be correlated with the jet position angle at 15 GHz and 43 GHz (Fig. 4 (b')).What can be clearly seen in Fig. 4 (c), the inner-jet position angle and the jet position angle have a significant negative correlation at 2.3 GHz; a significant positive correlation at 8.6 GHz and 15 GHz.In segment II, There is still a correlation between the two at a significance level greater than 95%.Fig. 4 (d/d') show the extremely weak or even indistinguishable correlations between the core EVPA and the inner-jet position angle at 15 GHz and 43 GHz; and the MC results show that we have no evidence of a correlation between the two as their peak significance is so low.Inevitably, due to the insufficient sampling rate, interpolation of uneven sample data can affect the above statistical results.

DISCUSSION
In this section, we present interpretations of the results of correlations between observables through two themes.The first theme discusses the correlations between the core flux density, the inner-jet position angle, and the jet position angle.Another theme discusses the correlation between the core EVPA and the inner-jet position angle.The higher the frequency of observation, the closer the core position is to the jet base, i.e., the 'core shift' effect (Blandford & Königl 1979;Konigl 1981;Sokolovsky et al. 2011;Pushkarev et al. 2012).It makes sense that the radio core identified at different resolutions is seen as the different part of the body of the 'true' jet and the jet position angle at the core scale (i.e., inner-jet position angle) changes with the jet swings rather than standing still.
The results that the core flux density and the innerjet position angle have a clear correlation conform to the prediction of the jet precession model or jet precessionnutation model proposed by Tateyama & Kingham (2004) and Britzen et al. (2018), respectively.This is because both models consider that the radiation variability in the radio jet could be attributed to the geometrical effects, that is, the precessing jet motion is accompanied by the viewing angle variations and consequent Doppler beaming changes.Meanwhile, OJ 287 has the intrinsic helical shape in the source frame due to the jet precession mechanism; and this source exhibits the bending jet morphology in the observer's frame, especially at higher observing frequencies.In our visual inspection, we can see apparent bent jets at 43 GHz, while at lower observing frequencies, the visual bending effect is not obvious (Fig. 1).
Moreover, geometric effects affect the lower frequency radio emission more strongly than higher frequency radiation.This is because the high-frequency radiations arise from higher-energy electrons with shorter cooling time scales, which can naturally lead to stronger intrinsic variability than at lower frequencies.When all data have the same excellent data quality, our correlation analysis can confirm the influence of geometric effects in different bands.Currently, under the premise that 2.3 GHz and 8.6 GHz data have basically the same time coverage and sampling rate, the correlation results show that the correlation between the core flux density and the inner-jet position angle at 8.6 GHz is significantly weaker than their correlation at 2.3 GHz; so we can conclude for the time being that the weakened correlation at high frequencies is most likely due to the differences in the geometric effect of different bands.

The core flux density and the jet position angle
As mentioned in the previous subsection, we can infer that a correlation between the total flux density and the jet position angle for all bands if the source has a precession mechanism.To verify this inference, we calculated the correlation between the total flux density and the jet position angle in four bands.In this calculation, for individual epoch observations, the total flux density is estimated as the sum of the flux densities of all fitted components.Looking at Fig. 7, it is apparent that the total flux density and the jet position angle show a significant positive correlation at 2.3 GHz and a clear negative correlation at 8.6 GHz and 15 GHz.The correlation coefficient values of the total flux density and the jet position angle at the cross-correlation peak of the core flux density and the jet position angle are 0.709, −0.685 and −0.546 (segment I)/0.624 (segment II)/ at the 2.3 GHz, is the result from segment-II data.Arrows located at time delays corresponding to cross-correlation peaks.The arrow with a question mark on the left side of the third panel in Panel (a) identifies the position of the maximum peak.However, based on our principle of determining the significant peak, in this case, the significant peak is considered to be at the arrow close to 0. The shaded area indicates parts with a correlation coefficient greater than (or less than) the 0.8 times maximum peak of the cross-correlation function.
At 2.3 GHz, 8.6 GHz, and 15 GHz, we notice that the direction of the correlation between the core flux density and the jet position angle is consistent with that between the total flux density and the jet position angle.At 2.3 GHz, the amplitude of the correlation between the core flux density and the jet position angle is basically the same as that between the total flux density and the jet position angle.However, at 8.6 GHz and 15 GHz, the amplitude of the correlation between the core flux density and the jet position angle differs from that between the total flux density and the jet position angle.
We interpret these results as the emission of core dominance increases towards lower frequencies.The main reason for this core dominance property is the limited jet length detected at 43 GHz, 15 GHz, 8.6 GHz, and 2.3 GHz and progressively decreasing angular resolution that results in blending in the core region swallowing more and more jet emissions.In other words, the con-tribution of jet components flux density to the total flux density is minimal at 2.3 GHz and becomes greater at higher frequencies.Two pieces of evidence support this contribution characteristic, i.e., the core dominance gradually increases with decreasing frequency.The first is a near-perfect linear relationship between the core flux density and the total flux density appearing at 2.3 GHz, not at relatively higher frequencies (see Fig. 8).The second is the source compactness factor, defined as the ratio of core flux density to total flux density, which decreases with increasing frequency.The mean values of the compactness factor are 0.964, 0.875, 0.905, 0.755 for 2.3 GHz, 8.6 GHz, 15 GHz and 43 GHz, respectively.
For the core dominance gradually decreases as the frequency increases, the correlation between the total flux density and the jet position angle at 8.6 GHz and 15 GHz is not transferred to the correlation between the core flux density and the jet position angle.Mathematically, suppose the correlation coefficients among three random variables A, B, C are ρ AB , ρ AC , ρ BC , respectively.If two correlation coefficients are known (e.g., ρ AB and ρ AC ), the value range of the third correlation coefficient (ρ BC ) is: ρ BC ρ AB ρ AC + 1 − ρ 2 AB 1 − ρ 2 AC .So through the correlation coefficient sequence of the total flux density and the jet position angle and the correlation coefficient of the core flux density and the total flux density at the time delay equal to 0, the value range of the correlation coefficient sequence of the core flux density and the jet position angle can be limited.We verify that our correlation coefficient sequence of the core flux density and the jet position angle is entirely within the limited range.
Interestingly, the correlation between the core flux density and the jet position angle shows the reverse correlation direction of 2.3 GHz and 8.6 GHz.We propose a plausible jet configuration based on the precession mechanism to explain this reverse correlation (see Fig. 9).Because the resolution of each observing frequency is different, the positions and scales of the fitted components for each band are also different, and the emerging components observed in some bands may not be resolved at lower resolutions.Generally speaking, the precession model considered that the superluminal fluid elements are ejected from a jet 'nozzle' that precesses around a fixed axis.Stirling et al. (2003) have found that some bright components evolve in straight paths, and their ejection angles are consistent with the 'nozzle' direction at the ejection time.
Obviously, the curvature of the upstream jet is greater, while that of the downstream jet is less.It is conceivable that the curvature of the downstream jet is decreased by the previous components moving in the ejection direction.As shown in Fig. 9, the geometric effect caused by the above physical reasons can lead to a configuration where the jet position angle changes in the opposite direction when observing the 2.3 GHz and 8.6 GHz.However, the angle variation between the whole jet and our line of sight is consistent in each band.Therefore, we can capture that the correlation direction between the core flux density and the jet position angle is exactly opposite at 2.3 GHz and 8.6 GHz, i.e., one positive correlation and one negative correlation.Distinctive features of the precessing jet include a gradually-curved helical jet.The definition of the innerjet and main jet direction depends on the resolution and thus on the frequency of observation.The fitted components for different epochs on different observing frequencies are located at different positions of the S-like shaped apparent jet with complex curvature.Hence the evolution of the jet position angle is also complicated; this situation does apply to the inside of the fitted core.The ad hoc configuration of the unresolved core components may lead to an assembled kinematics of the inner-jet, which is either in-phase or out-of-phase with those of the main jet.For this reason, if the jet has a direction change near the base, the correlations with changing direction of the jet may be positive or negative depending on the viewing angle and the effective resolution, as new features emerge from the core region.

Correlation analysis between the core EVPA and the inner-jet position angle
The observed polarization could be attributed to the nature of the locally ordered magnetic field.If the magnetic field is tangled and isotropic, the linear polarization integrated and averaged in some parts of the jet would be absent, and the direction of the polarization is also disorganized.In theory, the EVPA is orthogonal to the synchrotron magnetic field in optical thin regions.Virtually, all polarized emissions we observe in radio jet images arise in the optically thin regions, including polarization from the core region (Gabuzda 2021).
Usually, the AGN core EVPA is moderately alignment with the local jet direction, especially the BL Lac objects (Gabuzda & Cawthorne 2000;Hodge et al. 2018).However, it should be mentioned that some research findings on whether the EVPA and jet orientations in the core region tend to align are inconsistent, OJ 287 is no exception.For OJ 287, some studies investigated the offset between core EVPA and local jet position angle for individual epoch (D'arcangelo et al. 2009), while others explored the average results of the evolution of their offset over time (Sasada et al. 2018;Hodge et al. 2018).Sasada et al. (2018) reported that the differences between the core EVPA and the local jet position angle are scattered between 0 • to 90 • from June 2007 to March 2016 at 43 GHz.Hodge et al. (2018) reported that the median absolute value and variance of the differences between EVPA and jet position angle in the core region is 70 • and 30 • , respectively.These two statistics were calculated from a sample of 97 observation epochs over 11 years at 15 GHz.The results above consistently conclude that the angle between the core EVPA and the local jet position angle of OJ 287 is not a relatively stable value for OJ 287.
We re-investigated the long-term variation of offset between the core EVPA and the inner-jet position an-gle4 , which is shown as a function of time (left panels of Fig. 10).These results show that the offset is not always a stable value in either band, which is consistent with the previous studies.Moreover, it can be noticed that the sudden jump in the offset occurred in early 2009 at 15 GHz.The inner-jet position angle caused this sudden jump in the offset because it underwent rapid rotation in January 2009, the core EVPA only shows erratic rather than abrupt variation over time, especially during 2009.The offset distributions of the two frequencies are visualized by histograms (right panels of Fig. 10), where the 15 GHz offset dataset is divided into two groups.At 43 GHz, the offset's probability distribution exhibits Jshaped distribution; and at 15 GHz, it consists of two J-shaped distributions that are mirror images of each other.
We have not found a strong evidence correlation between these two observables because the variations of the core EVPA and the inner-jet position angle cannot keep constant during long-term monitoring.This sign shows that EVPAs in the 15 GHz and 43 GHz unresolved core region are less ordered.

Significance of correlation analysis
The correlation analysis is crucial for further understanding the radio radiation variability causes and the jet magnetic field configurations.
This fact that the changing flux density and jet position angle are correlated may be of geometric origin.Namely, the jet direction change can explain the variability of the radio flux density via viewing angle changes and the Doppler beaming.The precession is a special and important mechanism that causes a change in the jet direction that brings out this geometric connection.A common method for researchers to determine the presence of jet precession is to extract the jet position angles from the VLBI radio observations and use the obtained data to estimate precession model parameters (Britzen et al. 2018;Dominik et al. 2021).The better and more convenient way to search for AGNs with geometrically connected characteristics is to search for correlations between multiple observables, which can greatly improve the efficiency of searching for candidate sources with the jet precession.
The previous study has shown that the jet position angles of MOJAVE AGNs typically do not vary on the level of the core EVPAs, OJ 287 is an exceptional source that shows jet position angle change at the core scale in the plane of the sky (Lister et al. 2013;Hodge et  2018).In our study, the inner-jet position in OJ 287 shows time-varying characteristics at both 15 GHz and 43 GHz.In the long-term multi-band observations, rotation in the inner-jet and core polarization directions provides unique opportunities to gain insight into the polarization and magnetic properties of the innermost jet.

SUMMARY
This paper carries out a multivariate correlation analysis of observables of the OJ 287 which is a candidate with precessing jet.We analyzed the correlation between the core flux density, the inner-jet position angle, and the jet position angle for four frequency bands and that between the core EVPA and the inner-jet position angle for 15 and 43 frequency bands.The primary conclusions are summarized as follows.
In the four bands, there is a clear correlation between the core flux density, the inner-jet position angle and the jet position angle.The correlations that have been found can be explained by changes in the jet direction possibly due to precession.We emphasize that the correlations between observables found in OJ 287 are not an accidental property of this particular dataset, and the above correlations can be found if a source has precession characteristics.Admittedly, interference with correlation from relatively poor data quality is inevitable.
The offset of the core EVPA and the inner-jet position angle varies with time; at 15 GHz, the offset jumps rapidly at this moment because of the sudden rotation of the inner-jet position angle in 2009.In the 15 GHz and 43 GHz bands, there are no significant correlation between the core EVPA and the inner-jet position angle.These results may not support the assertion that the magnetic field becomes more ordered down the jet, even though this property is clearly demonstrated by Pushkarev et al. (2017), who reported an increase in fractional polarization with core separation.Nevertheless, we emphasize that longer-term VLBI monitoring is needed to confirm or refute this assertion.
of West China Program of the Chinese Academy of Sciences (2017-XBQNXZ-A-008).
This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018).
This study makes use of VLBA data from the VLBA-BU BLAZAR Monitoring Program (BEAM-ME and VLBA-BU-BLAZAR; http://www.bu.edu/blazars/BEAM-ME.html), funded by NASA through the Fermi Guest Investigator Program.The VLBA is an instrument of the National Radio Astronomy Observatory.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc.
Monitoring observations A collection of archived multi-band multi-epoch VLBI data from several geodetic/astrometric and astrophysical VLBI programs are analyzed in this study.It includes the S/X-band (2.3 GHz/8.6 GHz) data from Research and Development VLBA (RDV) program, the K u -band (15 GHz) data from the MOJAVE program and the Q-band (43 GHz) data from the VLBA-BU-BLAZAR program.

Figure 1 .Figure 2 .
Figure1.CLEANed images of OJ 287 for the first, median, and last epoch of multi-epoch data at different frequencies.The lowest I contours are plotted at the ± 3 σ (σ is the off-source noise) level of all images and the positive contour levels increase by a factor of 2. The linearly polarized intensity contours in blue overlayed in I contours and the lowest P contours are plotted at the ± 6 σ level of all images and the positive contour levels also increase by a factor of 2. The center and semi-major and semi-minor axes of the red ellipse represent the position and size (FWHM) of the Gaussian components fitted in the image plane, respectively.The solid green lines indicate the electric polarization vector directions, the length of this vector is proportional to the intensity of the linear polarization with 1 mas corresponding to 0.0625 Jy/beam at 15 GHz and 0.25 Jy/beam at 43 GHz.The FWHM restoring beam is represented by the circle in the bottom left corner of each panel.

Figure 3 .
Figure 3. Evolution of the four observables of four bands with time.Panel (a): Multi-band VLBI core flux density plotted as a function of time.Panel (b): Multi-band inner-jet position angle plotted as a function of time.Panel (c): Multi-band jet position angle plotted as a function of time.Panel (d): Multi-band core EVPA plotted as a function of time.Two double-headed arrows with the text 'SEGMENT I' or 'SEGMENT II', and a background of different colors identify the segments considered for analysis in our work.

Figure 4 .
Figure 4. Cross-correlation analysis of each pair observables at 2.3 GHz, 8.6 GHz, 15 GHz and 43 GHz.The left panes are the results of segment I, and the right panels are the results of segment II.Panel (a/a'): core flux density vs. inner-jet position angle.Panel (b/b'): core flux density vs. jet position angle.Panel (c/c'): inner-jet position angle vs. jet position angle.Panel (d/d'): core EVPA vs. inner-jet position angle.The positive time delay corresponds to the latter observable's time series leading the former observable when doing the cross-correlation, and vice versa.Arrows located at time delays corresponding to cross-correlation peaks.

5. 1 .
Correlation analysis between the core flux density and the jet position angle at different spatial scales 5.1.1.The core flux density and the inner-jet position angle

Figure 5 .
Figure 5. Cross-correlation significance results.The color lines represent the cross-correlation for each pair observables, while the color contours show the distribution of random cross-correlations obtained by the Monte Carlo simulation with 1σ (dark color/inner layer near horizontal axis), 2σ (light color).Panel (a/a'): core flux density vs. inner-jet position angle.Panel (b/b'): core flux density vs. jet position angle.Panel (c/c'): inner-jet position angle vs. jet position angle.Panel (d/d'): core EVPA vs. inner-jet position angle.Panel (a/b/c/d) is the result from segment-I data, and Panel (a'/b'/c'/d') is the result from segment-II data.Arrows located at time delays corresponding to cross-correlation peaks.The arrow with a question mark on the left side of the third panel in Panel (a) identifies the position of the maximum peak.However, based on our principle of determining the significant peak, in this case, the significant peak is considered to be at the arrow close to 0. The shaded area indicates parts with a correlation coefficient greater than (or less than) the 0.8 times maximum peak of the cross-correlation function.

Figure 6 .
Figure 6.MC error analysis results of correlations between each pair observables.Panel (a/a'): The ∆ ρCF,IA represents the correlation coefficient error values between each band's core flux density and inner-jet position angle.Panel (b/b'): The ∆ ρCF,JA represents the correlation coefficient error values between each band's core flux density and jet position angle.Panel (c/c'): The ∆ ρIA,JA represents the correlation coefficient error values between each band's inner-jet position angle and jet position angle.Panel (d/d'): The ∆ ρEA,IA represents the correlation coefficient error values between each band's core EVPA and inner-jet position angle.The black line is connected by the mean values of the correlation coefficient error values of each time delay slice.The red lines is connected by the above mean values plus and minus 3 times standard deviation (σ), where σ is determined by a Gaussian fit of the distribution of correlation coefficient error values over each time delay slice.Panel (a/b/c/d) is the result from segment-I data, and Panel (a'/b'/c'/d') is the result from segment-II data.

Figure 7 .Figure 8 .
Figure 7. Cross-correlation analysis of the total flux density and the jet position angle.The short line marks the location of the cross-correlation peak of the core flux density and the jet position angle.

Figure 9 .
Figure 9. Sketch of a projection of the sky plane of a helical jet caused by precession.The 8.6 GHz (X-band) and 2.3 GHz (S-band) jet configurations are located in the left and right panels, respectively.The black and gray circles represent the individual fluid elements in epoch 1 and epoch 2, respectively.The yellow circle represents the emerging fluid element, we only draw an emerging fluid element for clear visualization.Each black arrow connects each black circle with the corresponding grey circle, indicating the projected direction of the ballistic motion of each fluid element in space.The blue oval indicates the fitted core component, and the green oval and pink oval indicate the fitted jet component at epochs 1 and 2. The red line roughly represents the direction of the jet, i.e., jet position angle.
Figure 10.Time evolution of the offset of the core EVPA and the inner-jet position angle for OJ 287 at 15 GHz and 43 GHz.At 15 GHz, before January 2009, the offset values are marked in black triangles, and the offset values after January 2009 are marked in red triangles.Distributions of the offset of the EVPA and the inner-jet position angle at 15 GHz (Panel (c)) and 43 GHz (Panel (d)).At 15 GHz, the black border histogram shows the offset distribution of all data between January 1996 and January 2009 and the red border histogram shows the offset distribution of all data between February 2009 and June 2020.

Table 2 .
Jet components properties.Columns are as follows: (1) date of observation; (2) observing frequency in GHz; (3) fitted jet component number; (4) the jet component's relative right ascension shift to the core component in mas; (5) the jet component's relative declination shift to the core component in mas; (6) jet position angle derived by fitting regression line.(This table is available in its entirety in machine-readable form.) Kikuchi et al. (1988)here are few research reports on the EVPA day-level rapid rotation event of OJ 287 in the radio band.To date, merelyKikuchi et al. (1988)'s study reported the EVPA varied by 80 • over-5 days in the radio region for OJ 287.To sum up the above, for the observation cadence of the datasets in this study, the method we used can better extract the correct core EVPA values.

Table 3 .
Correlation results between observables for each observing frequency.Columns are as follow: (1) two observables for calculating the correlation; (2) observing frequency; (3-4) the correlation coefficient value of significant peaks; their corresponding time delays are in parentheses, in years.In the first column, 'CF' represents core flux density, 'IA' represents inner-jet position angle, 'JA' represents jet position angle, 'EA' represents core EVPA.