Polarized Optical Emission of the Blazar PKS 1222+216: Discovery of A 420-day Quasi-Periodic Signal

We report our search for quasi-periodic signals in long-term optical and $\gamma$-ray data for the blazar PKS~1222+216, where the data are from the Steward Observatory blazar monitoring program and the all-sky survey with the Large Area Telescope onboard the {\it Fermi Gamma-ray Space Telescope}, respectively. A quasi-periodic signal, with a period of $\simeq$420\,days and a significance of $>5\sigma$, is found in the measurements of the optical linear polarization degree for the source, while no similar signals are found in the optical and $\gamma$-ray light curves covering approximately the same time period of $\sim$10\,yr. We study the quasi-periodic variations by applying a helical jet model and find that the model can provide a good explanation. This work shows that polarimetry can be a powerful tool for revealing the physical properties, in particular the configuration of the magnetic fields of jets from galactic supermassive black holes.


INTRODUCTION
In recent years, the availability of rich data over nearly the whole electromagnetic spectrum has greatly enabled studies of variable and transient phenomena in our universe. A particular and interesting one is quasi-periodic oscillation (QPO) seen in Active Galactic Nuclei (AGN). The QPO phenomenon in the past has been frequently seen in X-ray observations of Galactic X-ray binaries (e.g., Ingram & Motta 2019 and references therein). Now as long-term monitoring data at multiple wavelength bands are available, allowing searches for different time-scale variability, quite a few long-periodicity QPOs have been reported (see, e.g., Zhang & Wang 2021 and references therein).
Given the periodicties and likely transient nature (i.e., only appearing during certain time periods), AGN QPOs are discussed to reflect the accretion activity in the innermost stable circular orbit around an SMBH (Gierliński et al. 2008), the quasi-periodic variations of the accretion disk surrounding an SMBH, or the binary nature of the SMBH system in the center (e.g., King et al. 2013;Ackermann et al. 2015, and references therein). As approximately 10% of AGN have jets and the jets can contribute with a significant fraction of the observed emission, particularly from those so-called blazars (whose jets point close to our line of sight), the observed flux variations including QPOs can be predominately caused by the variations of the jets. Since the jets are considered to be coupled with the accretion disk, their variations still reflect the physical activity in an AGN.
The AGN jets are relativistic, with the bulk Lorentz factor Γ b ∼ 10. For blazars, the large and dominant flux variations are seen arising from their jets, boosted by the Doppler beaming effect. Thus if a jet has a helical structure, an emitting blob in this jet that is moving helically would induce periodic flux modulations, as our viewing angle to the blob is changing periodically (e.g., Rieger 2004;Sobacchi et al. 2017). This geometric model has been applied to a few QPO cases such as those in the blazar PKS 2247−131 (Zhou et al. 2018) and in the narrow-line Seyfert 1 Galaxy J0849+5108 (Zhang & Wang 2021). Accompanying the flux modulation, it is conceivable to think that the polarization of the emission would be changing accordingly, as the magnetic field is considered to have a helical configuration and the synchrotron radiation, which is responsible for emission in approximately from infrared to X-ray wavelengths from a jet, is highly polarized (see e.g., discussion in Chen & Zhang 2021). Indeed, there have been studies and discussions for the related polarization changes (e.g., Lyutikov et al. 2005;Raiteri et al. 2013; Zhang 2019, and references therein). Previously Otero-Santos et al. (2020) have even reported a 4σ QPO signal (with a period of ∼1.9 yr) in both the optical broadband photometric and polarimetric measurements for the blazar source B2 1633+38. Thus polarimetry can provide insights into jets' physical properties . Moreover if QPOs are detected, they could help revealing the configuration of the magnetic fields of jets.
In this paper we report another QPO case found in the polarized optical emission from a blazar, PKS 1222+216 (or 4C+21.35). This source is a flat spectrum radio quasar (FSRQ) subtype blazar, with redshift z = 0.434 (Ahn et al. 2012). It has been detected with the Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope (Fermi) with name J1224.7+2121 in the Fermi-LAT first source catalog (1FGL; Abdo et al. 2010), and it was selected as a target in the Steward Observatory blazar monitoring program (Smith et al. 2009). Using the data from the optical spectropolarime-try and photometry carried out by the program, we conducted the analysis of the source's long-term variations. The Fermi-LAT 0.1-500 GeV data for the source were also analyzed, since the γ-ray band might provide additional information for its overall multi-band variations. Below we describe the data and analysis in Section 2, and present the QPO search results in Section 3. The results are discussed in Section 4.

Optical data description
Nearly 100 blazars were monitored by the Steward Observatory for 10 Fermi mission cycles from 2008 October to 2018 July. PKS 1222+216 was one of the targets and for it the linear polarization measurements at a wavelength range of 5000-7000Å and photometric Vand R-band brightness measurements are available. The time range of the data extends approximately between MJD 54948-58306. We used the degree of polarization (DoP) measurements, whose minimum (at MJD 55385.18), maximum (at MJD 57127.25), and mean values in percentage are 0.25, 29.1, and 5.47, respectively, and standard deviation is 366.45%. The photometric V -and R-band brightness measurements are simultaneous, and the V -band (Rband) magnitudes have a minimum (at MJD 56715.39), maximum (at MJD 55892.51), and mean value of 13.95 (13.82), 17.79 (17.84), and 16.64 (16.58), respectively. The standard deviation is 69.69% (76.28%) in V -band (R-band). The DoP and magnitude data are shown in Figure 1, and their median and mean time intervals between two adjacent data points are ∼1 and 7 days, respectively. It can be noted that there is an approximate three-month gap in the data of each year, obviously because of the observations being ground-based. We selected the Fermi -LAT 0.1-500 GeV Pass 8 Front+Back SOURCE class photon-like events between 2008 August 4 and 2021 November 18 in a region of 20 • ×20 • (region of interest; RoI), centered at the position of PKS 1222+216. We removed the events with zenith angle > 90 • , and reduced the others with expression DATA QUAL>0 && LAT CONFIG = 1 to obtain high-quality data in the good time intervals.

Fermi-LAT data analysis
A model file was created with the script make4FGLxml.py 1 based on 4FGL, which included best-fit spectral parameters of all known 4FGL sources in the RoI. We set a few parameters free, which were the parameters of flux normalizations and spectral shapes for the sources within 5 • from the target, and the normalizations for the sources within 5 • -10 • and the ones outside 10 • but identified with significant variations (i.e., their Variablility Index 72.44). The normalizations of the Galactic and extragalactic diffuse emission components were also set free. The other parameters were fixed at the values given in 4FGL.
A binned maximum likelihood analysis was first employed to build a best-fit model file for the RoI. PKS 1222+216 has a log-parabola spectral shape in . From the likelihood analysis, we obtained a test statistic (TS) value of ∼83,470 and an average photon flux of (267.85 ± 2.12) × 10 −9 photons cm −2 s −1 in 0.1−500 GeV band from the whole data. The best-fit spectral parameters were α s = 2.324 ± 0.007, β s = 0.029 ± 0.004, and E b = 392.380±10.219 MeV. The values are in well agreement with those reported in the 4FGL data release 3 (DR3; Fermi-LAT collaboration et al. 2022). Based on the best-fit model file, a 0.1-500 GeV light curve of PKS 1222+216 was constructed by performing an unbinned likelihood analysis. In this step, we fixed the parameters of spectral shapes of the sources at their best-fit values obtained above and set the normalizations of each sources in the RoI free. We tested different time bins, such as 15, 30, and 45 day, but the obtained light curves did not have significant differences. In the bottom panel of Figure 1, the 15-day binned light curve is presented. In order to show the light curve as complete as possible, low-significant data points with 25 > TS 4 are also included.

PERIODIC VARIABILITY ANALYSIS AND RESULTS
We searched for periodicity signals in the optical and γ-ray data. Two methods were used, namely the weighted wavelet Z-transform 2 (WWZ; Foster 1996) and the generalized Lomb-Scargle Periodogram 3 (LSP; Lomb 1976;Scargle 1982;Zechmeister & Kürster 2009). The latter was used to check the results from the former as an independent method.

Optical polarization
In the optical DoP data, we found a ∼420-day periodic signal. As clearly seen in the WWZ time-frequency power spectrum shown in the left panel of Figure 2, the signal was present over the whole-data time range and is significant in the time-averaged power spectrum (right panel of Figure 2). The LSP power spectrum shows a similar result, containing a sharp peak at the same period. Here we adjusted the parameter for the wavelet transform in the WWZ analysis to make the two power spectra as much as possible match with each other. We estimated the period value by fitting the power peak with a Gaussian function, which yielded 420.1±55.6 day, where the uncertainty was taken as the full width at half maximum of the peak.
A periodic signal or statistical fluctuations in a set of time series data may lead to a power peak in an LSP power spectrum. The probability (p) for the power peak arising purely from the statistical fluctuations was estimated to be < 7.8×10 −35 (with respect to the local noise floor around the peak), according to the calculation method given in Horne & Baliunas (1986) and Zechmeister & Kürster (2009). The false alarm probability (FAP) considering the period search in the frequency range is < 8.4 × 10 −33 , given by FAP = 1 − (1 − p) N , where N (=108) is the number of independent frequencies in the frequency range (i.e., the trial factor) and calculated from (f max −f min )/δ f . Since the time series data are unevenly spaced, we approximately set f max 1/(30 day) and f min 1/(1640 day) (where 1640 day ∼ half of the total length of the data), and δf is the frequency resolution, determined by the total length of the data.
We employed a light-curve simulation method to determine the significance of the periodicity signal. The LSP power spectrum was modeled with a function of a smoothly bending power law plus a constant (González-Martín & Vaughan 2012) to estimate the underlying power spectral density (PSD). The function has the form of P (f ) = Kf −α 1 + (f /f bend ) β−α −1 +C, where K, α, β, f bend , and C are the normalization, low frequency slope below bending, high frequency slope above bending, bending frequency, and Poisson noise, respectively. Their best-fit parameters were obtained with a maximum likelihood algorithm provided in Barret & Vaughan (2012), and we found K = 0.56 ± 0.80, α = 0.66±0.28, β = 1.56±0.73, f bend = 0.0003±0.0023, and C = 4.43 ± 0.69 (the model fit is shown in the right panel of Figure 2). Based on this model, 2 × 10 7 artificial light curves were generated (Emmanoulopoulos et al. 2013). For each artificial light curve, we obtained its PSD with the above LSP periodic analysis. Significance curves were calculated by counting the PSD data points at each of the frequencies. The resulting 5σ significance curve is shown in the right panel of Figure 2, and as can be seen, the LSP power peak of the 420-day periodicity exceeds this significance curve. After using the same method as the above for taking into account the trial factor, the signal is still greater than 5σ. As to a few minor peaks also seen in the power spectra, we found that none of them had a sufficiently high significance, particularly when the trial factor was taken into account.
We considered a sinusoidal modulation curve for approximately describing the periodic variations, A sin [2πf (t − t 0 )] + A 0 , where A, f, t 0 , and A 0 are the amplitude, frequency, time of zero phase, and additive offset of the sine function, respectively. The values of the best-fit parameters were provided by the LSP code, which were 2.49 ± 0.18, 0.00238 ± 0.00002 day −1 , 54935.54 ± 4.81, and 5.44 ± 0.13, respectively. This sinusoidal curve is shown in the panel C of Figure 1. Approximately 8 cycles of periodic modulation appeared in the data.
Although the 420-day value is different from a year of 365 day, we investigated whether the yearly observational gaps could induce any artifactual signals. We conducted two simulations by generating DoP data points with their times set as the observational ones. The DoP values were set as a constant (5.47, the mean of the observed data), or random ones in a range of 0-15 (the DoP values of most of the data points). In the first simulation, there were small peaks in the power spectrum, but none of them appears significant and more importantly around 420 day. In the second simulation, there were only background noises. We thus concluded that the 420-day signal could not be induced by the yearly observational gaps.

Optical and γ-ray light curves
The same analysis for periodicity search was conducted to the optical photometric V -and R-band data. The obtained power spectra are shown in Figure 3. In the power spectra resulting from both bands, a possible period of ∼500-day appears. While it can be considered to overlap with that from the DoP data, it is broad and not significantly among other nearby weak peaks.
The analysis on the Fermi-LAT γ-ray data yielded similar results to the above optical data (Figure 3). A weak signal appears, with its periodicity consistent with the error range found for the 420-day periodicity, but compared with nearby peaks, it is not significant at all. We concluded that no significant signals were found in either the optical or γ-ray light curve data.

DISCUSSION
Taking advantage of the availability of long-term monitoring data at different wavelength bands, we have analyzed the optical and γ-ray data taken for the blazar PKS 1222+216. In the 10-year long optical photometric and polarimetric data, we have found a QPO signal of period P QPO 420 day in the DoP measurements of the latter. While there are peaks at ∼ P QPO present in the power spectra derived from the optical photometric and Fermi-LAT γ-ray data, they are too weak to be considered as possible signals. Also very-highenergy emission from the source was detected with the MAGIC telescopes in the energy range of 70 GeV up to ∼400 GeV during its active high-energy γ-ray variation period (Aleksić et al. 2011), and multi-wavelength studies of the source have clearly indicated that the major γ-ray variations, such as that around MJD 55300 (Figure 1), were flaring events (Ackermann et al. 2014;Bhattacharya et al. 2021). Non-detection of periodic signals should be expected when the γ-ray flux variations were dominated by flaring activity.
QPO activity seen from a blazar should arise from its jets, no matter what the underlying cause is, for example due to a binary SMBH system or a precessing accretion disk. This QPO in the optical polarized emission of a blazar presents an interesting case among thusfar reported QPOs, as nearly all of them were found in flux variations. Different from most of the others, this QPO possibly reveals the properties of the magnetic field in the source's jet (Raiteri et al. 2013;Marscher et al. 2017). A geometric model has been applied to explain the often-reported quasi-periodic flux variations considering a helical jet and magnetic field (Chen & Zhang 2021). An emitting blob in such a jet moves along the helical path, causing periodic changes in our view angle θ to the jet and as a result in the observed flux, since the Doppler beaming factor is a function of θ (e.g., Zhou et al. 2018).
Similar scenarios have been proposed to explain the polarization variations seen in the synchrotron emission of blazars (Lyutikov et al. 2005), where helical magnetic fields are considered. Following Raiteri et al. (2013), we assume that DoP P changes as a function of the viewing angle θ in the jet rest frame, P = P max sin 2 θ , where P max is the maximum DoP value and is shown to be ∼ 20 in Lyutikov et al. (2005). We tested to fit the DoP measurements with this function, for which sin θ = sin θ/Γ b (1 − β cos θ) (βc is the bulk speed). The viewing angle θ is derived from cos θ = sin φ sin ψ cos(2πt/P QPO ) + cos φ cos ψ, where φ is the pitch angle of an emitting blob's motion from the jet and ψ is the inclination angle of the jet to the line of sight. According to the detailed studies of the components in the jet at the center of PKS 1222+216 with the Very Long Baseline Array (VLBA) observations , the opening semi-angle of its jet is ∼ 2 • and thus we set φ = 2 • . In our fitting, we also fixed Γ b = 15, which was obtained from the radio measurements and the spectral energy distribution modelling reported in Jorstad et al. 2017 andBhattacharya et al. 2021 respectively. Since polarized emission of a jet can consist of multiple components (e.g., Marscher et al. 2017), we did not intend to fit the observed DoP values exactly and thus set P max as a free parameter. The other one set free was ψ. We obtained the best fit with P max = 9.37 ± 0.03 and ψ = 8.22 ± 0.02 deg. This fit, shown in Figure 1, can describe the overall DoP variations and is similar to the sinusoidal fit given by the LSP code. In this geometric model, the view angle changes around ψ 8 • , which is close to the average value 5.6 • ± 1.0 • estimated from the VLBA studies ). Therefore the model provides a reasonable explanation to the QPO seen in the DoP measurements, in which the magnetic field of the jet should have a helical configuration.
We note that the nearly simultaneous optical light curves, unlike what was reported for the B2 1633+38 case (Otero-Santos et al. 2020), do not show a QPO signal or similar variation patterns (cf., Figure 1). This can be explained as that the optical emission contains a significant contribution arising from the accretion disk, as the fitting results from the studies of the broadband spectral energy distribution have shown that the disk emission could be dominant (Ackermann et al. 2014;Bhattacharya et al. 2021). Alternatively, if the optical emission mainly arises from the jet, the magnetic field could consist of a turbulent component, which then explains the observed DoP values ∼5%, much lower than that expected from the synchrotron emission of highly ordered magnetic fields . This case of the QPO signal thus likely reveals a significant helical component of the magnetic field in the jet of PKS 1222+216, and shows that polarimetry can be a powerful tool for studies of blazar jets by overcoming apparently noisy brightness variations.
According to the geometric model, the jet of PKS 1222+216 was moving towards the Earth along a helical path with a cycle time of [P QPO /(1 + z)]/(1 − β cos φ cos ψ) 62 yr at the host galaxy (cf., Zhou et al. 2018). We observed ∼8 QPO cycles in 10 yr, which implies that the emitting blob moved with a distance of ∼150 pc in total or a projected distance of ∼21 pc. The luminosity distance to PKS 1222+216 is estimated to be 2480 Mpc (where the cosmological parameters from the Planck mission, H 0 67.4 km s −1 Mpc −1 , were used; Planck Collaboration et al. 2020). Thus the corresponding projected angular distance would be ∼1.8 mas. Comparing to the proper motion values of 0.01-0.5 mas yr −1 derived from the VLBA observations for several components of the jet in PKS 1222+216 , the estimated projected distance is within the range, showing the consistency of the model result with that of the high-resolution radio imaging.
As a summary, we have found a 420-day periodicity in the optical linear DoP measurements for the blazar PKS 1222+216. The QPO signal is significant at a > 5σ level, higher than that previously reported in B2 1633+38. We have applied a helical jet model to explain the quasi-periodic variations, and the parameters estimated from the model are consistent with those derived from the VLBA studies of the central radio jet. This QPO case shows that while the flux variations of a blazar may appear random, the physical properties of its jets can still be investigated and revealed through polarimetry. It is challenging to carry out polarimetry observations because of the need for detecting sufficient amount of photons, but when QPOs shown in flux vari-ations of blazars are found, close polarimetry monitoring should be warranted as the observation can potentially provide critical information for understanding the mechanism that drives the QPOs and thus jets' physical properties.