© The Authors 2026

1. Introduction

Active galactic nuclei (AGNs) are central compact regions of the galaxies that emit electromagnetic radiation ranging from radio to γ-rays. AGNs are among the most powerful sources in the universe, containing supermassive black holes (SMBHs) with masses ranging from ∼10⁶ to 10¹⁰ M_⊙ at their centers, surrounded by an accretion disk, and relativistic jets perpendicular to the disk. Magnetic fields are believed to be responsible for the collimation and acceleration of relativistic jets throughout their propagation. These jets emit nonthermal synchrotron radiation, resulting from the interaction of relativistic particles with magnetic fields. The emission radiated from the jets is highly polarized and variable with characteristic timescales of hours to years (Krolik 1999). Very long baseline interferometry (VLBI) has been extensively used to study the innermost regions in AGNs, providing high-resolution imaging capabilities. In particular, utilizing data from the VLBA-BU-BLAZAR monitoring program (Jorstad & Marscher 2016), several studies have investigated the temporal evolution of jets or the relationship between γ-ray variability and jet kinematics (e.g., Jorstad & Marscher 2010; Hodgson et al. 2021; Rani et al. 2015). The interferometric monitoring of gamma-ray bright AGNs (iMOGABA) employs the Korean VLBI network (KVN) at 22, 43, 86, and 129 GHz (Lee et al. 2016), and has contributed to our understanding of various physical properties of AGNs with bright γ-rays (e.g., Lee et al. 2020; Kim et al. 2022; Park et al. 2019).

Blazars, a subclass of AGNs, belong to Type 0 (Urry & Padovani 1995), which have very unusual spectral characteristics (including high luminosity and strong variability) with a small viewing angle (e.g., less than a few degrees) of jets relative to the line of the sight, causing substantial Doppler boosting. A typical spectral energy distribution (SED) exhibits two humps: one attributed to synchrotron radiation at relatively low frequencies, and the other to inverse Compton (IC) radiation at higher frequencies, arising from the scattering of photons, by relativistic electrons, from synchrotron emission or external components (e.g., accretion disk, broad-line region (BLR), or dusty torus). Blazars are divided into two subclasses according to their spectral features: flat-spectrum radio quasars (FSRQs) and BL Lac objects (BL Lacs). FSRQs are highly luminous and emit strong radiation, with broad emission lines, indicating the presence of a dense BLR, while BL Lacs show weak or absent emission lines. Magnetic fields in jets play a crucial role in shaping the physical properties of blazars, influencing both jet formation and high-energy emission processes (Marscher et al. 2008). Sokolovsky et al. (2010) showed that the upper limits of the magnetic field strength in blazars range from 10⁻² to 10² Gauss. Many studies have measured the magnetic field strengths of blazars using various methods, including those based on the synchrotron self-absorption (SSA) spectrum (Lee et al. 2017; Kang et al. 2021; Jeong et al. 2023; Kim et al. 2022), the core shift effect (O’Sullivan & Gabuzda 2009; Kutkin et al. 2014), and the synchrotron luminosity (Kim et al. 2017). Blazars exhibit γ-ray flares when superluminal jet components pass through the core, or stationary features, which are often associated with recollimation shocks within the jet (Marscher et al. 2008; Kovalev et al. 2009; Li et al. 2024). The studies have shown that these γ-ray flares occur at locations far from the central engine (e.g., ≳1 parsec), while the origin of the seed photons responsible for the emission remains under debate.

PKS 1222+216 (also known as 4C +21.35) is a FSRQ with redshift z = 0.435 (Wang et al. 2004), located at α = 12^h24^m54.458^s,  δ = +21° 22′46.388′ in the J2000 position. Its luminosity distance is D_L = 2.4 Gpc (Jorstad et al. 2015), and the central SMBH of the source has a mass of 6 × 10⁸ M_⊙ (Farina et al. 2012). The source has been monitored by Fermi-LAT, and γ-ray flaring events were detected in 2010, 2014, and 2019 (Kramarenko et al. 2022). Among these events, the 2010 flare was identified at very high energies (VHE; E > 100 GeV) by MAGIC (Aleksić et al. 2011). Jorstad et al. (2015) asserted that the γ-ray maxima coincide with interactions between moving jet knots and stationary structures, such as the core and A1, which are interpreted as recollimation shocks. They considered that the absence of a spectral cutoff constrains the γ-ray emission region, causing it to lie outside the BLR to avoid the pair-production absorption of VHE γ-rays. This supports the study presented in Tavecchio et al. (2011), where the emission from a compact blob and a large region located beyond the BLR adequately explains the SED during the 2010 γ-ray flaring period. Lee et al. (2019) analyzed the trajectory of a newly identified component, C, at 22 GHz using the KVN and VLBI Exploration of Radio Astrometry (VERA) array (KaVA) from 2014 to 2016, and suggested that component C may be associated with the γ-ray flare observed in 2014. Kinematic analyses of the jet components in PKS 1222+216 have been extensively conducted using VLBI imaging techniques. Troitskiy et al. (2016) identified five superluminal knots (9–22 c) and a stationary component A1 (0.14 ± 0.04 mas from the core) using 43 GHz VLBA data from 2008 to 2015, including the emergence of knot K6 in 2015. Lee et al. (2019) analyzed a subsequent epoch (2014–2016) using KaVA, presenting the kinematics of two jet knots with apparent speeds of 3.5 c to 6.8 c at 43 GHz. Weaver et al. (2022) analyzed the motion of jet components at 43 GHz for PKS 1222+216 using data collected from 2007 to 2018. The study showed physical parameters for each component, including the timescale of variability, Doppler factor (δ), Lorentz factor (Γ), and viewing angle. Fourteen components were identified, with the parameters for components B9, B10, and B11 estimated as δ ∼ 3–6, Γ ∼ 3–13, and with a viewing angle of approximately 6°. In addition, several studies have reported estimates of the magnetic field strength. Pushkarev et al. (2012) applied core-shift measurements between 8 GHz and 15 GHz to probe the physical conditions in AGN jets. For PKS 1222+216, they estimated the magnetic field strength to be 0.9 G at a distance of 1 pc from the central engine, and inferred the location of the 15 GHz core to be 23.41 pc from the central engine. Lee et al. (2019) found that the decay timescales of the jet flux densities at 22 and 43 GHz ranged from ∼0.5 yr to 3 yr, resulting in the magnetic field strength B of about 10–30 mG for PKS 1222+216.

In this paper we present observational results based on long-term, multi-frequency (22–86 GHz) KVN data obtained from 2013 to 2020, supplemented by VLBA radio data, with a focus on constraining the magnetic field strength under synchrotron cooling conditions and investigating the relationship between γ-ray flares and radio emissions. In Section 2 we describe the observation details and how to reduce the data. Our analysis results are shown in Section 3, and our results are discussed in Section 4. Finally, we summarize our study in Section 5.

2. Observations and multi-wavelength data

2.1. KVN observations and data reduction

We obtained multi-frequency simultaneous KVN data for the target source PKS 1222+216 as part of the KVN key science program (iMOGABA). The observations were conducted at 22 (K), 43 (Q), 86 (W), and 129 GHz (D) bands utilizing the KVN Yonsei, Ulsan, and Tamna stations. The observations were conducted monthly, except during maintenance seasons from June to August every year. The full bandwidth of 256 MHz was divided into four frequencies, providing a bandwidth of 64 MHz for each frequency in left circular polarization. The source was observed for one to eight scans of 5 minutes in length, distributed over several hours. Further detailed observations can be found in Lee et al. (2016). The data presented here cover the period from December 4, 2012 (MJD 56265), to March 6, 2020 (MJD 58914). A total of 60 epochs were obtained at K and Q bands, while 63 and 41 epochs were available for W and D bands, respectively. The KVN data were processed by the DiFX software correlator at the Korea Astronomy and Space Science Institute (KASI) (Lee et al. 2014). The data were also calibrated with the AIPS package utilizing an automatic pipeline for KVN data (Hodgson et al. 2016), including the procedures of phase and amplitude calibrations, fringe fitting, and bandpass correction.

We used the pipelined KVN data to obtain CLEAN images via the DIFMAP software (Shepherd et al. 1994). Amplitude outliers were flagged, and the imaging process began with the self-calibrating of the visibility data based on a point-source model of 1 Jy using the DIFMAP command STARTMOD. Subsequently, CLEANing and phase self-calibration were performed iteratively until the residual maps were left with only Gaussian random noise, i.e., a uniform distribution of residual noise with comparable positive and negative maxima. We obtained CLEAN images of the target source in 54 epochs on the K and Q bands, in 53 epochs on the W band, and 23 epochs on the D band. The number of images on the D band is less than half that on the K, Q, and W bands due to the lack of baselines; accordingly, we used only the K, Q, and W bands for the data analysis. The typical amplitude calibration errors applied here are 10% for K and Q bands and 20% for the W band. Fig. 1 shows the CLEAN image from the KVN data reduced for the 2016 January epoch. The CLEAN beam sizes are typically 5.66 × 3.35 mas, 2.77 × 1.74 mas, and 1.52 × 0.79 mas at 22, 43, and 86 GHz, respectively.

Fig. 1. CLEAN images of PKS 1222+216 obtained by KVN observations at 22, 43, and 86 GHz on January 13, 2016. The radio cores are aligned to the (0,0) position. The color and contours indicate the intensity of each map. The contours start at three times the root mean square (RMS) of the residual map and increase by a factor of 1.4. The RMS levels are 0.012, 0.019, and 0.013 Jy/beam at 22, 43, and 86 GHz, respectively. Synthesized beams are plotted in the lower right corner of each map as a gray ellipse.

2.2. VLBA data

PKS 1222+216 has been observed at 43 GHz with the VLBA under the VLBA-BU-BLAZAR program (Jorstad & Marscher 2016), which monitors 34 blazars and three radio galaxies. We selected 15 VLBA epochs from the period between April 12, 2015, and September 5, 2016 (MJD 57124–57636), to complement the 43 GHz flux density measurements obtained with the KVN. The VLBA 43 GHz data were pre-calibrated with CLEAN models available on the program website¹. We obtained the CLEAN total flux densities of PKS 1222+216 from these CLEAN models. Typical amplitude calibration uncertainty is estimated to be 5%. For kinematic analysis, we also selected 21 VLBA epochs spanning from July 28, 2014, to September 5, 2016 (MJD 56866–MJD 57636), to investigate the inner region of the source, taking advantage of the superior angular resolution of 0.1 mas (Jorstad et al. 2017). We fitted circular Gaussian models to the data using the MODELFIT task of DIFMAP, taking into account the findings of previous studies (Jorstad et al. 2015; Troitskiy et al. 2016; Weaver et al. 2022). We iteratively added Gaussian components to the model, accepting each new component only if it yielded a χ² reduction of at least 10%. To maintain continuity, the final model from the previous epoch was used as the initial guess for the subsequent epoch. Fig. 2 shows the results of the component model fitting, where the yellow circles represent the full width at half maximum (FWHM) size of each component.

Fig. 2. CLEAN maps at VLBA 43 GHz from July 2014 to September 2016. The contour starts from three times the RMS noise level and increases by a scale factor of 1.4. The yellow circles indicate the sizes (FWHM) of the individual components at each observation. The gray ellipses indicate the synthesized beam sizes for each epoch.

The 15 GHz flux density of PKS 1222+216 was obtained from the Monitoring Of Jets in AGN with the VLBA Experiments (MOJAVE) program, referring to the long-term and systematic monitoring of relativistic jets in AGNs on parsec scales using VLBA (Lister et al. 2018). We selected eight VLBA observations taken from June 16, 2015, to March 11, 2017 (MJD 57189–57823). We utilized the total CLEANed Stokes I flux density in the VLBA image at 15 GHz, publicly available on the program website². The typical amplitude calibration error was estimated to be 5%.

2.3. [γ]-ray data

The Fermi Large-Area Telescope (LAT) Light Curve Repository (LCR) provides γ-ray light curves binned on timescales of three days, one week, and one month for 1525 sources, including PKS 1222+216 (Abdollahi et al. 2023)³. From this repository, we retrieved the γ-ray light curve of PKS 1222+216 with a one-week cadence and energy range of 0.1–100 GeV from December 7, 2012 (MJD 56268), to March 6, 2020 (MJD 58914).

3. Results

3.1. Multi-wavelength light curves

We obtained CLEAN total flux density light curves at 15, 22, 43, and 86 GHz from MOJAVE and iMOGABA, covering a time span of approximately seven years from December 2012 to March 2020 (MJD 56265–58914), as shown in Fig. 3. The source appears compact within the KVN data (see Fig. 1) and as such it is most likely that the CLEAN total flux is dominated by the emission from the compact region. To enhance the light curve cadence at 43 GHz, we incorporated the CLEAN total flux densities of the BU 43 GHz data. We found that the multi-wavelength CLEAN total flux densities exhibited a declining trend at all frequencies from early 2015 to early 2017 during the observation period. Following the methodology described in Section 3.4, we determined the decay periods of the flux densities for each wavelength (gray boxes in Fig. 3). Table 1 summarizes the range of total flux densities, the decay periods, and fractional flux density decreases for each frequency. In the following sections, we focus on a detailed examination of the physical properties during the cooling period.

Fig. 3. Light curves of PKS 1222+216 observed from December 2012 to March 2020 (MJD 56265–58914) at 15, 22, 43, and 86 GHz. The total CLEAN flux densities are shown for each frequency. In the third panel, the green circles represent the iMOGABA data, while the olive circles correspond to the BU data. The error bars used typical values of 10% for K and Q, and 20% for W. BU and MOJAVE were applied at 5% (Jorstad et al. 2017). The bottom panel shows the 0.1–100 GeV γ-ray light curve from the Fermi-LAT LCR. The upper limit of the flux density is indicated by the gray triangles. The gray vertical dashed line refers to the γ-ray peak at MJD 56975.5 (Ezhikode et al. 2022). The gray boxes are the periods that fit the exponential decay profiles at each frequency.

We also utilized γ-ray data ranging from MJD 56268 to MJD 58914, a period coinciding with the iMOGABA observations (including the upper limits). The light curve peaked on MJD 56975.5 with a peak flux density of 3.2 × 10⁻⁶ ± 1.20 × 10⁻⁸ photon cm⁻² s⁻¹ and a MJD flaring period of 56974.6–56977.8, as defined in Ezhikode et al. (2022). We found that after the γ-ray flare, the radio light curves peaked at MJD 57189 (15 GHz), MJD 57107 (22 GHz), MJD 57124 (43 GHz), and MJD 57142 (86 GHz), followed by a monotonic decline, as described above. After a significant flare in 2014, the flux densities began to rise again during 2018–2019, peaking initially at 86 GHz, in agreement with the γ-ray data. This intriguing behavior could be explored further in future studies.

3.2. Spectral indices of the iMOGABA flux density

We analyzed multi-wavelength iMOGABA data to investigate the spectral properties of PKS 1222+216 and derived spectral index values from simultaneous KVN observations. We obtained the spectral index α at two adjacent frequencies (frequency ν₁ and ν₂) as

$$\begin{array}{c}\hfill \alpha =\frac{log(F_{{\nu }_2}/F_{{\nu }_1})}{log({\nu }_2/{\nu }_1)},\end{array}$$(1)

where F_ν is the CLEAN flux density at the observing frequency at ν and ν₂, ν₁ (ν₂ > ν₁) are the observing frequencies.

The top panel of Fig. 4 shows the spectral indices for 53 epochs with corresponding average values of −0.60 ± 0.15 and −0.59 ± 0.23 at 22–43 GHz and 43–86 GHz, indicating that the source is optically thin at 22–86 GHz on the mas-scale. To explore the spectral properties of the variable emission regions, we analyzed the source in both quiescent and active states. The quiescent state refers to periods without significant flux enhancements, and the quiescent flux was defined as the minimum CLEAN total flux density observed during the monitoring period: 1.21, 0.69, and 0.45 Jy at the K, Q, and W bands, respectively, on MJD 58850 (K and Q) and 58888 (W). Fig. 5 shows the quiescent spectral index. To account for the effect of quiescent states, we subtracted the quiescent flux from the total CLEAN flux to obtain the active flux, from which the active spectral index was subsequently derived. The active spectral indices were measured and found to be −0.19 ± 0.71 at 22–43 GHz and −0.37 ± 0.97 at 43–86 GHz (see the bottom panel of Fig. 4). We found that the active spectral index at 43–86 GHz increased between MJD 57076 and MJD 57289, transitioning from optically thin to optically thick conditions, before returning to optically thin conditions (α = −0.78) on MJD 57319. This dynamic behavior may be due to the kinematic interaction of the 43 GHz components, but a detailed investigation should be conducted in future studies.

Fig. 4. Radio spectral indices at 22–43 GHz and 43–86 GHz. The top panel represents the values derived from the CLEAN total flux, while the bottom panel corresponds to the active flux. The orange scatters represent the spectral index at 22–43 GHz, while the green scatters indicate the spectral index at 43–86 GHz. The horizontal dashed line indicates that the spectral index is zero, and the dotted line in the top panel is the mean value with α = −0.6. In the bottom panel, the orange dotted line represents the average active spectral index of −0.19 at 22–43 GHz, while the green dotted line represents the average active spectral index of −0.37 at 43–86 GHz. The gray box indicates the decreasing period from 57124 to 57636 at 43 GHz.

Fig. 5. Spectral index of PKS 1222+216 in the quiescent state. The black points represent the quiescent-state flux densities at 22, 43, and 86 GHz, with error bars indicating flux uncertainties. The orange slope corresponds to the quiescent spectral index between 22 and 43 GHz (αKQ = −0.81), while the green slope represents the quiescent spectral index between 43 and 86 GHz (αQW = −0.628).

3.3. BU 43 GHz component analysis

In order to investigate the inner jet region (less than 1 mas) and its physical parameters, we conducted a detailed component analysis using the BU program data. The components were identified based on the results in Weaver et al. (2022), following the assumption that the geometric distance from the core and flux density of each component does not vary rapidly over time (see Appendix B for a detailed comparison and identification process). Table A.1 presents the model fitting parameters for each component. In Fig. 6, the results of the component analysis are shown, including the distances from the core and the flux densities for each component.

Fig. 6. Upper panel: Radial distance of 43 GHz jet components relative to the core, which is shifted to the phase center at (0,0) position, from April 2015 to September 2016. The black arrows represent their position angles, and the solid lines represent their fitted trajectories. The colors correspond to the newly ejected components: C9 (blue), C10 (magenta), C11 (green), and C12 (cyan). The gray box highlights the interaction period between components C9 and A1. The horizontal yellow and pink lines indicate the mean positions of stationary components A1 and A2, respectively. Middle panel: Flux density of the core (red) and the newly ejected components, C9, C10, C11, and C12. Bottom panel: Flux density of components with relatively low flux density: C5, C6, C8, C12, A1, and A2. The black vertical line at MJD 56975 represents the epoch of the γ-ray flare.

We identified ten components, including new knots, C9, C10, C11, and C12, which appeared on November 15, 2014 (MJD 56976), April 12, 2015 (MJD 57124), September 22, 2015 (MJD 57287), and July 31, 2016 (MJD 57600), respectively. We found that the γ-ray flare peak (see Section 3.1) nearly coincided with the appearance of C9 with a one-day offset. C9 appears to correspond to a new component, K6, that was ejected in 2015, as mentioned in Troitskiy et al. (2016), as K6 appeared after the high γ-ray peak. Interestingly, the position angle of the inner region (within 1 mas) significantly changed from −3.95° (Epoch 2, A1) to 26.23° (Epoch 3, C9) when component C9 was ejected. As shown in Figs. 2 and 6, the jet in the inner region appears to bend eastward after the third epoch, which is associated with the γ-ray flare. This agrees with the 43 GHz KaVA results from 2016 by Lee et al. (2019), which also revealed a ∼20° westward tilt in the inner region (see Figure 3). We also discovered stationary components A1 and A2 at 0.15 ± 0.01 mas and 0.35 ± 0.05 mas from the core, respectively, and which are denoted by the yellow and light red dotted horizontal lines in Fig. 6. These findings are consistent with previous studies (Troitskiy et al. 2016; Weaver et al. 2022; Jorstad et al. 2015). Additionally, the average geometric distances from the core for the components C5, C6, C8, C9, C10, and C11 are 2.32 ± 0.35, 1.78 ± 0.36, 1.05 ± 0.44, 0.65 ± 0.39, 0.63 ± 0.36 and 0.52 ± 0.28 mas, respectively. The fluxes of C9, C10, and C11 decreased rapidly during the observation period, with corresponding reductions of 4.67%, 11.11%, and 11.29%, during the decay period described in Section 3.4

3.4. Decay timescales

We found that all light curves at each observation frequency decreased during the observation period. Following Terasranta & Valtaoja (1994) and Burbidge et al. (1974), we undertook exponential profile fitting to the light curves to determine the logarithmic variability timescale. We used the following exponential decay model,

$$\begin{array}{c}\hfill F_{\nu }=A_{\nu }exp(-B_{\nu }t)+C_{\nu },\end{array}$$(2)

where F_ν is flux density; A_ν, B_ν, and C_ν are constants; and B_ν = 1/τ_ν, with τ_ν the decay timescale at frequency ν. The nonlinear least-squares method was utilized in curve_fit of SciPy, which is a fundamental algorithm for scientific computing in Python (Vugrin et al. 2007).

It is crucial to select the precise observation period during the exponential fitting to obtain accurate physical quantities and ensure reliable fitting results. We determined the decay period at each frequency using the method described below. We selected data only from the period of the overall large decrease in the flux densities found during the iMOGABA observation period. Starting from the peak of the light curves, the data were grouped, including four epochs subsequently, after which the group’s standard deviation and mean value were calculated. Data grouping was stopped when the mean value of the current group exceeded the mean value of the previous group. The decay period ended at the epoch with the minimum value in the final group.

For the iMOGABA data, the selected decaying periods are MJD 57107–57680 (11 epochs) at 22 GHz, MJD 57142–57623 (9 epochs) at 43 GHz, and MJD 57142–57719 (10 epochs) at 86 GHz. For the VLBA 43 GHz, we selected the data from MJD 57124–57680 (15 epochs) for the decaying period. For the 15 GHz observations of MOJAVE, we used the data covering MJD 57189–57823 (8 epochs). For the model-fitting results of the 43 GHz BU data, the selected epochs for the C9, C10, and C11 components are MJD 56976–57549 (16 epochs), MJD 57123–57573 (13 epochs), and MJD 57205–57465 (7 epochs), respectively. In the second panel in Fig. 6, the solid lines represent the results of exponential fits to the flux decay of each component.

The best-fitting decay timescales obtained using the above method are summarized in Table 2, yielding 404 ± 130 days at 15 GHz, 401 ± 88 days at 22 GHz, 227 ± 62 days at 43 GHz, and 150 ± 42 days at 86 GHz. The timescales for the C9, C10, and C11 components are 108 ± 6, 43 ± 4, and 222 ± 88 days, respectively, as presented in Table 3.

3.5. Physical parameters of the jet components

We identified three newly emerging superluminal knots (C9, C10, and C11) in the source at 43 GHz. The estimated physical parameters and their corresponding uncertainties for these components are summarized in Table 3. By fitting a linear function to the angular distance from the core over time, we determined the proper motion, μ_ang for each component. The corresponding apparent speed, β_app (in units of the speed of light), was found to be approximately 13.4 for all of the newly ejected components, with a typical uncertainty of ∼0.5. We assumed that the decay in the light curve is attributable to synchrotron cooling (e.g., Lee et al. 2019, and see Appendix C). We regard the decay timescale τ as the synchrotron cooling timescale τ_cool. Based on the assumption of the synchrotron cooling knots, we can use the cooling timescales τ_cool to determine the variability Doppler factor δ_var of the components as (Jorstad et al. 2005, 2017)

$$\begin{array}{c}\hfill {\delta }_{\mathrm{var}}=\frac{15.8s{D}_{\text{L}}\text{L}}{{\tau }_{\mathrm{cool}}(1+z)},\end{array}$$(3)

where s is the angular size of the component in mas, taken as 1.8 times the FWHM for an optically thin spherical brightness distribution (Pearson 1995); D_L is the luminosity distance to the source in Gpc; and τ_cool is the cooling timescale in years. The variability Doppler factors of the superluminal knots range from 4 ± 1.8 to 22 ± 2.4.

In addition to the variability Doppler factor obtained from the Equation (3), we also sought to estimate the Doppler factor using the brightness temperature T_b of the components as

$$\begin{array}{c}\hfill T_{\text{b}}=1.22\times {10}^{12}\frac{S}{s^2{\nu }^2}(1+z)\mathrm{K},\end{array}$$(4)

where S is the flux density at the peak in Jy, s is the angular size in mas, ν is the observing frequency in GHz (i.e., 43 GHz), and z is the redshift of the source (z = 0.435). The Doppler factor δ is then estimated as the ratio of the maximum intrinsic brightness temperature T_int to the observed brightness temperature as follows (Lähteenmäki & Valtaoja 1999):

$$\begin{array}{c}\hfill \delta =\frac{T_{\text{b}}}{T_{\text{int}}}·\end{array}$$(5)

We assume that T_int ≈ T_eq and adopt the equipartition value of 5 × 10¹⁰ K suggested in Readhead (1994). The brightness temperature varies between (0.57 ± 0.02)×10¹¹ K and (0.99 ± 0.04)×10¹¹ K, with the Doppler factors of the superluminal knots ranging from 1.1 ± 0.1 to 2.0 ± 0.1. The viewing angles Θ are calculated with the apparent speeds and Doppler factors from Equation (5) as follows:

$$\begin{array}{c}\hfill \mathrm{\Theta }=arctan\left(\frac{2{\beta }_{\mathrm{app}}}{{\beta }_{\mathrm{app}}^2+{\delta }^2-1}\right).\end{array}$$(6)

The viewing angles between the jet axis of the superluminal knots at 43 GHz range from 8.2 ± 0.2° to 8.6 ± 0.4°.

The projected apparent speed β_app is related to the intrinsic speed β and the viewing angle Θ with respect to the line of sight as follows:

$$\begin{array}{c}\hfill \beta =\frac{{\beta }_{\mathrm{app}}}{\sin \mathrm{\Theta }+{\beta }_{\mathrm{app}}\cos \mathrm{\Theta }}·\end{array}$$(7)

The intrinsic speeds approach the speed of light for all components (β ≈ 0.999), yielding Lorentz factors of Γ ∼ 45 ± 3–78 ± 7.

3.6. Magnetic field strengths

When a gyrating electron radiates power, the synchrotron cooling time can be defined via τ_cool = E/Ė, with the electron energy E = γ_em_ec², where γ_e is the electron Lorentz factor, m_e is the electron mass, c is the speed of light, and Ė is dE/dt with time t. Using the Larmor formula and assuming an isotropic distribution of the pitch angle, Ė can be expressed in CGS units as $⟨P⟩=\frac{4}{3}{\sigma }_{T}c{U}_B{\gamma }_{\mathrm{e}}^2{\beta }_{\mathrm{e}}^2$, where σ_T is the Thomson scattering cross section, U_B is the magnetic energy density with B²/8π and B is the magnetic field strength (Rybicki & Lightman 1991). After calculating all constants in SI units, we can express the cooling time equation for the magnetic field strength in the observed frame as follows:

$$\begin{array}{c}\hfill B=\sqrt{7.74\frac{1+z}{\delta {\gamma }_{\mathrm{e}}{\tau }_{\mathrm{cool}}}}·\end{array}$$(8)

Assuming that the dominant synchrotron power is emitted at the peak frequency ν_syn of the synchrotron spectrum and that the IC spectrum is dominated by the synchrotron self Compton (SSC), the Lorentz factor of the electrons dominantly contributing to the SSC can be estimated as (Beckmann & Shrader 2012)

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{e}}=\sqrt{\frac{3{\nu }_{\mathrm{IC}}}{4{\nu }_{\mathrm{syn}}}},\end{array}$$(9)

where ν_IC is the peak frequency of the IC spectrum. Referring to the spectral analysis results for the 2014 flare of PKS 1222+216 as reported in Ezhikode et al. (2022) and Adams et al. (2022), the peak frequencies of the synchrotron and IC humps during the flaring period differ by approximately eight orders of magnitude (ν_IC/ν_syn ≈ 10⁸), yielding an estimated electron Lorentz factor γ_e of 8660. This value is close to a value located at the upper end of the interval typically covered by FSRQs. Furthermore, the magnetic field strength was determined to be 19 ± 3 mG at 15 GHz, 19 ± 2 mG at 22 GHz, 26 ± 4 mG at 43 GHz, and 32 ± 5 mG at 86 GHz, as presented in Table 2. These were calculated based on the redshift z = 0.435, the cooling timescales in Section 3.4, and the Doppler factors (δ₂₂ = 10.7 ± 0.4, δ₄₃ = 9.6 ± 0.6) derived from the inverse-variance weighted mean of jet components in Lee et al. (2019). The magnetic field strength values of C9, C10, and C11 were estimated to be 83 ± 3, 134 ± 8, and 77 ± 15 mG using the timescales determined in Section 3.4 and Doppler factor in Section 3.5, as summarized in Table 3. As a result, we found that the magnetic field strength decreases as a function of the frequency and that the region of the C10 component at 43 GHz has the most significant magnetic field strength.

4. Discussion

4.1. Connection between the [γ]-ray flare and radio components

Several investigations have indicated a correlation between the emergence of new jet components and flaring activity (Jorstad et al. 2001; Schinzel et al. 2012; Lico et al. 2022; Escudero Pedrosa et al. 2024). Motivated by this, we determined the times at which the newly emerged components C9, C10, and C11 detached from the core by fitting a linear function to their measured separations. The ejection times from the core for the components are as follows: MJD 56932 ± 8 for C9, MJD 57067 ± 9 for C10, and MJD 57167 ± 12 for C11. Additionally, we determined the interaction times with the stationary component A1 based on the trajectory fitting results, defined as the time at which the moving components passed through the centroid position of A1, located at 0.15 ± 0.01 mas from the core: MJD 57011 ± 11 for C9, MJD 57143 ± 11 for C10, and MJD 57246 ± 14 for C11. In the first panel in Fig. 6, the vertical dashed lines indicate the epochs when each component was ejected from the core, while the solid vertical lines mark the estimated interaction times between the moving components and A1. The emergence of these new components from the core is not temporally coincident with the γ-ray flare observed at MJD 56974.6–56977.8. For instance, the flare occurred 43 days after the emergence of C9 and 92 days before the emergence of C10. However, given that the interaction time between C9 and the stationary component A1 is only 36 days after the γ-ray flare, we can consider the possibility of interaction between the A1 and C9 components, affecting the γ-ray flare.

We estimated the duration of the interaction between C9 and A1 based on the angular sizes of A1 and C9, yielding an interaction period from MJD 56961 to MJD 57061 (as indicated by the gray box in Fig. 6). The γ-ray peak time (MJD 56975) coincides with this interaction period, suggesting that the γ-ray flare may be related to the interaction. This exhibits characteristics consistent with prior studies of interactions between the core or stationary components and moving components induced by a conical standing shock (e.g., Jorstad et al. 2015; Li et al. 2024). Therefore, we can suggest that the new jet component C9 passes through the stationary component A1, as illustrated in Figure 7 of Li et al. (2024), accelerating synchrotron particles for IC scattering, which up-scatters low-energy photons (e.g., BLR photons) into γ-ray photons.

Unlike the situation with C9, a γ-ray flare does not occur when C10 and C11 traverse A1 (e.g., after MJD 57123). When C10 and C11 emerged, the synchrotron particles may not have been accelerated for IC scattering (e.g., optimal shock environments are not formed) as were those for the interaction between C9 and A1. As shown in Fig. 2, the inner region of the jet appears to bend eastward starting from the third epoch, when component C9 emerges (Table A.1; position angle changes from approximately −4° for A1 to 26° for C9). In contrast, no significant change in the jet direction is observed during the ejection of components C10 and C11. While this suggests a possible connection between the γ-ray flare and the jet bending, further investigation (e.g., a high-angular-resolution polarization study) is required to confirm this relationship in the future.

Another possibility involves the acceleration for component C10. Previous studies (e.g., Homan et al. 2009; Schinzel et al. 2012; Sasada et al. 2018) have shown evident acceleration of jet components associated with standing shocks. Schinzel et al. (2012), for example, reported that 43 GHz VLBI components in 3C 345 accelerate over ∼23 pc after passing the core, coinciding with γ-ray emission. Analogously, the γ-ray flare could arise during the passage of C10 through the core, followed by a ∼90 day period of acceleration toward A1, after which the component maintained a constant speed as identified in our kinematic analysis. Alternatively, the γ-ray flare could occur farther upstream, within the BLR or dusty torus, followed by the physical acceleration of the plasma, which eventually emerges as component C10. However, as described in the following Section 4.2 regarding the estimated locations of the 43 GHz core and dusty torus, such upstream regions are expected to be highly opaque to high energy photons. The dense external radiation fields from the dusty torus would lead to significant γγ pair production, effectively absorbing the γ-rays before they can escape (Donea & Protheroe 2003). Therefore, we propose the interaction between C9 and A1 as the most plausible scenario for the 2014 γ-ray flare, given that its downstream location provides a more transparent environment for γ-rays, in addition to aligning with both the shorter time lag between the flare and C9’s separation and the observed jet bending.

4.2. Estimation of the [γ]-ray flare location and the seed photon

In AGN research, γ-ray emission is likely to be generated through the IC scattering of seed photons by relativistic particles. The seed photons may be synchrotron photons originating within the jet itself for SSC scattering or may come from external sources such as BLR, dusty torus, and cosmic microwave background (CMB) radiation.

To investigate the origin of the seed photons for the IC process during the γ-ray flare, we constrained the location of the γ-ray flare caused by the interaction between the moving jet component C9 and the stationary component A1. The flare location can be estimated by determining the distance of A1 from the central engine in parsecs. Following the model outlined by Lobanov (1998), the position of the VLBI core is determined under the assumption that the core is the optical depth τ = 1 surface. The distance of the core from the central engine, r_core, depends on the observing frequency ν, as described by r_core ∝ ν^−1/k_r, where k_r is a parameter depending on the magnetic field and particle density distribution along the jet. The absolute distance of the central engine can be estimated as

$$\begin{array}{c}\hfill r_{\text{core}}(\nu )={\mathrm{\Omega }}_{r\nu }{[{\nu }^{\frac{1}{k_{\text{r}}}}sin\mathrm{\Theta }]}^{-1},\end{array}$$(10)

where Ω_rν is the core shift measure obtained from the angular offset Δr_mas in mas between core positions at two frequencies, ν₁ and ν₂ (ν₁ < ν₂):

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{r\nu }=4.85\times {10}^{-9}\frac{\mathrm{\Delta }{r}_{\text{mas}}{D}_{\text{L}}}{{(1+z)}^2}·\frac{{\nu }_1^{1/k_{\text{r}}}{\nu }_2^{1/k_{\text{r}}}}{{\nu }_2^{1/k_{\text{r}}}-{\nu }_1^{1/k_{\text{r}}}}·\end{array}$$(11)

It is expected for ideal measurements that Ω_rν is constant for all frequency pairs. We used a core shift measurement of 0.178 mas between 8.1 GHz and 15.4 GHz, as reported in Ezhikode et al. (2022), yielding Ω_r8, 15 = 7.08 pc GHz. By adopting a viewing angle of Θ = 5.3° (Pushkarev et al. 2017) and assuming k_r ≈ 1 (i.e., equipartition between the energy density of the jet particles and the magnetic field) in Equation (10), we derived the distance of the core at 43 GHz, r_core(43 GHz) = 1.765 pc. Given the averaged angular position of A1 at 0.15 ± 0.01 mas and a luminosity distance of 2.4 Gpc for the source, the projected distance of A1 from the core is estimated to be r_A1 = 9.35 ± 0.81 pc, considering the aforementioned jet viewing angle. Consequently, the absolute distance of A1 from the central engine is determined to be r_A1 = 11.11 ± 0.8 pc. The location of the γ-ray flare can also be estimated using the observed time delay, Δt_obs, between the γ-ray peak time and the moment when C9 passes through A1. Considering cosmological effects, the distance between the A1 component and the γ-ray emission region in the source frame can be expressed as follows:

$$\begin{array}{c}\hfill \mathrm{\Delta }r=\frac{{\beta }_{\text{app}}c\mathrm{\Delta }{t}_{\text{obs}}}{(1+z)sin\mathrm{\Theta }}·\end{array}$$(12)

Subsequently, the location of the γ-ray flare can be determined as r_γ = r_A1 − Δr_{C9, A1}. As the γ-ray peak time is used as MJD 56975 and the interaction time as MJD 57011, the location of the γ-ray flare is estimated to be approximately r_γ = 9.2 ± 1.0 pc. Taking into account the position and size of A1, the location of C9 at the moment when it first passes through A1 is estimated to be 10.95 ± 0.82 pc. Given the uncertainty range of the γ-ray flare location and the position where C9 intersects with A1, it is plausible to conclude that the γ-ray flare occurred as C9 passed through A1.

Furthermore, the outer radius of the BLR and the radius of the dusty torus reported by Ezhikode et al. (2022) are R_BLRout = 0.07 pc and R_DT = 2.27 pc, respectively. The estimated location of the 43 GHz core (∼1.8 pc) lies within this dusty torus region. As discussed in Section 4.1, if the γ-ray flare had originated at the core, the high energy emission would likely have been attenuated by absorption within the dusty torus. Moreover, the γ-ray flare location estimated from downstream jet components is significantly offset from the BLR and dusty torus by approximately 9 pc and 7 pc, respectively. Consequently, the seed photons participating in the IC mechanism for the γ-ray flare are less likely to have originated from the BLR or dusty torus and are most likely to have originated from the jet itself or from external CMB radiation. Alternatively, as proposed by Marscher et al. (2010), a slow sheath surrounding the jet spine could serve as a source of infrared seed photons produced by moving knots or standing shocks. This feasibility of high energy emission far from the central engine is consistent with the scenario where the 2010 γ-ray flare of PKS 1222+216 originated outside the BLR due to jet recollimation (Tavecchio et al. 2011).

4.3. Magnetic field strengths at multi-wavelengths

In Section 3.6, we determined the magnetic field strengths to be B∼ 20–30 mG in the radio core regions at frequencies of 15–86 GHz. This result is consistent with the typical blazar magnetic field of B∼ 10–100 mG (Lewis et al. 2016). We found that the magnetic field strength increases as the frequency increases.

The multi-wavelength magnetic field strengths and cooling timescales were fitted to a power-law function of the frequency (τ_ν ∝ ν^α and B_ν ∝ ν^β), yielding the best-fit solution of α = −0.56 ± 0.13 and β = 0.34 ± 0.04, as shown in Fig. 7. These results indicate that τ_cool decreases exponentially and the magnetic field strengths increase as a function of frequency for an optically thin plasma.

Fig. 7. Top panel: Power-law fit of the cooling timescales as a function of frequency (τν = A × να). Bottom panel: Power-law fit of the magnetic field strength as a function of frequency (Bν = B × νβ). The best-fit power-law indices are α = −0.56 ± 0.13 and β = 0.34 ± 0.04.

The relationship between the cooling timescale and frequency (τ_ν ∝ ν^{−0.56 ± 0.13}) is consistent with the theoretical prediction (τ_ν ∝ ν^−0.5) of an optically thin synchrotron plasma with critical frequency ν_c(∝γ³ω_B), an electron angular frequency ω_B(∝γ⁻¹), and cooling timescale τ(∝γ⁻¹) (see Equation (8)).

The magnetic field strength is proportional to the distance from the central engine as B ∝ r^−m, and the distance is proportional to the frequency as r ∝ ν^−1/k_r (Blandford & Königl 1979, O’Sullivan & Gabuzda 2010), yielding B ∝ ν^m/k_r. Assuming equipartition between the energy densities of jet particles and magnetic field energy densities with k_r = 1, the magnetic field strength B is expected to follow a radial profile proportional to r⁻¹ (B ∝ r⁻¹) following Konigl (1981). However, the observational data for PKS 1222+216 indicate that the magnetic field strength in the jet evolves as B ∝ ν^{0.34 ± 0.04}, indicating the quantity m/k_r ≈ 0.34. This implies that the emission region in the relativistic jet of PKS 1222+216 may not be in an equipartition condition (k_r ≳ 1) or may have a magnetic field strength profile slowly decreasing along the jet (m < 1).

5. Conclusions

In this study, we investigated the long-term decline of multi-wavelength emissions of PKS 1222+216 using radio data from the iMOGABA program, MOJAVE, and the VLBA-BU-BLAZAR program spanning 2013–2020, combined with γ-ray data from Fermi-LAT. We found that the flux densities at radio frequencies of 15, 22, 43, and 86 GHz had decreased exponentially by 37–56% over a year following a γ-ray flare in November 2014. By attributing this decay to synchrotron cooling, we estimated key physical parameters of the source during the observation period.

Although previous studies have presented model fitting results for this source, we conducted our own independent Gaussian model fitting using VLBA 43 GHz data from July 2014 to September 2016. We identified ten jet components, including three newly ejected components (C9, C10, C11) and two stationary components (A1, A2). The cooling timescales for the new components were found to range from 43 to 222 days, with viewing angles confined within approximately eight degrees. Using the electron Lorentz factor derived under the assumption of SSC emission, we estimated magnetic field strengths for each component ranging from 77 mG to 134 mG.

We found that the interaction between C9 and A1 as the most plausible scenario for the 2014 γ-ray event, based on its temporal coincidence with the γ-ray flare and the physical constraint that the core region remains opaque to high energy photons due to absorption within the dusty torus. To constrain the γ-ray emission site and the seed photon source for the IC process, we calculated the distance of component A1 from the central engine to be approximately 11.11 ± 0.8 pc, and the γ-ray flare location was estimated at 9.2 ± 1.0 pc. Given the absolute position and the size of component A1, taking into account the uncertainties, these results suggest that the flare occurred during the C9-A1 interaction, at a distance beyond the BLR and dusty torus region. Therefore, the seed photons for the IC process likely originated from the jet itself, the CMB, or a surrounding sheath, rather than from the BLR or dusty torus. Our findings support a scenario where γ-ray emission occurs farther downstream at stationary shocks, rather than close to the central engine.

Finally, we analyzed the magnetic field distribution based on optical depth variations, finding that the magnetic field scales as B ∝ r^{−0.34 ± 0.04} along the jet. This finding suggests that the emission region in the relativistic jet of PKS 1222+216 may deviate from the equipartition condition (k_r ≳ 1) or exhibit a slow decline in the magnetic field strength profile (m < 1). Our study demonstrates that analyzing cooling timescales offers a practical alternative method for determining magnetic field strengths in AGNs undergoing long-term synchrotron cooling.

Data availability

The KVN CLEAN images at 22, 43, and 86 GHz used in this study (shown in Figs. 1 and 3) are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/708/A339

Acknowledgments

We thank the anonymous reviewer for providing valuable comments that improved the manuscript. We are grateful to the staff of the KVN who helped to operate the array and to correlate the data. The KVN is a facility operated by KASI (Korea Astronomy and Space Science Institute). The KVN observations and correlations are supported through the high-speed network connections among the KVN sites provided by KREONET (Korea Research Environment Open NETwork), which is managed and operated by KISTI (Korea Institute of Science and Technology Information). This work has made use of Fermi-LAT data supplied by Abdollahi et al. (2023), https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository. This study makes use of VLBA data from the VLBA-BU Blazar Monitoring Program (BEAM-ME and VLBA-BUBLAZAR; 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. This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MIST) (2020R1A2C2009003, RS-2025-00562700).

References

Appendix A: Modelfitting parameters

We present the parameter values obtained from our independent model fitting in Table A.1. These include the total flux density, peak flux density, distance from the core, component size with respect to the FWHM, and the position angle of the major axis, each presented with the corresponding uncertainty. The uncertainties of the parameters were estimated following the method described in Section 2.4 of Lee et al. (2016). For the total flux density uncertainty, we adopted σ_tot = σ_peak (S_tot/S_peak).

Appendix B: Comparative analysis and identification of jet Components at 43 GHz

Although Weaver et al. (2022) presented the result of the Gaussian model fitting and component identifications for PKS 1222+216, we independently analyzed the 43 GHz BU data from mid-2014 to mid-2016, applying the criteria outlined in Section 2.2 and 3.3. We found fewer components in each epoch than Weaver et al. (2022). For example, they detected eight components with the core, but we detected six components on MJD 57123. We also detected C9 and C11 two epochs earlier than reported by Weaver et al. (2022). When a new component is detected in the vicinity of the stationary A1 region (r ≈ 0.14 mas), it is necessary to determine whether the component is a reappearing A1 component or a newly emerged jet component. In such cases, when a new component was detected at the typical position of A1 in the MJD 56976 epoch, we identified it as a newly emerged jet component (e.g., C9) based on its flux variation. In the previous epoch (MJD 56923), A1 had a flux density of 0.2 Jy, whereas the newly detected component exhibited a sudden increase to 1.1 Jy on MJD 56976. This abrupt flux enhancement indicates that the detected component is distinct from A1 and should be classified as a newly emerged jet component (C9). In addition to the flux variation, we also considered the distance from the core and the position angle. In the first two epochs, the component closest to the core was identified as the stationary A1, while the following component was classified as C9 due to its abrupt change in position angle (see Fig. 6). In the previous study, three components of B3, B5, and B6 were detected at distances greater than 1 mas from the core, whereas in our results, only two components, C5 and C6, were detected (i.e., non-detection of the B3 component). In Weaver et al. (2022), the stationary A1 component was detected starting from the 14th epoch (MJD 57388). However, in our results, A1 was revealed in the first and second epochs (MJD 56866 and 56923) and then reappeared in the 13th epoch (MJD 57361). The stationary A2 component was also detected one epoch later (MJD 57418) in our result compared to Weaver et al. (2022). Nonetheless, within the range of uncertainty, the average distance from the core and the flux values of A1 and A2 were consistent with the previous result (see Table B.1). We also found that the appearance epochs of C8 and C10 coincide with those reported in Weaver et al. (2022) (MJD 56866 and MJD 57124).

Appendix C: Synchrotron cooling versus adiabatic cooling

When investigating the magnetic field strength based on the cooling timescale, we assumed that the flux density decay is attributed to synchrotron cooling. However, in a relativistic jet, the particle energy density may decrease due to jet expansion, leading to adiabatic cooling. To ascertain whether this decrease was caused by adiabatic expansion or radiation, we investigated the time variation at size d (or the volume) of the jet components by fitting an exponential function to the data, yielding a jet expanding timescale of τ_d. The fitting results for each component are shown in Fig. C.1. The expansion timescale of the source size was found to be τ_d ∼ 10⁶ days, which is significantly longer than the cooling timescales of τ_cool = 43–222 days, as discussed in Section 3.4. This implies that the flux density decay is primarily attributed to synchrotron cooling rather than adiabatic cooling.

Fig. C.1. Each panel presents the results of the power-law fitting of the size evolution of each component over time. The panels, from top to bottom, show the results for C9, C10, and C11. The sizes (FWHM values) are shown as scatter points with error bars, while the best power-law fits are represented by solid black lines.

All Tables

All Figures

Fig. 1. CLEAN images of PKS 1222+216 obtained by KVN observations at 22, 43, and 86 GHz on January 13, 2016. The radio cores are aligned to the (0,0) position. The color and contours indicate the intensity of each map. The contours start at three times the root mean square (RMS) of the residual map and increase by a factor of 1.4. The RMS levels are 0.012, 0.019, and 0.013 Jy/beam at 22, 43, and 86 GHz, respectively. Synthesized beams are plotted in the lower right corner of each map as a gray ellipse.

In the text

	Fig. 2. CLEAN maps at VLBA 43 GHz from July 2014 to September 2016. The contour starts from three times the RMS noise level and increases by a scale factor of 1.4. The yellow circles indicate the sizes (FWHM) of the individual components at each observation. The gray ellipses indicate the synthesized beam sizes for each epoch.
In the text

Fig. 3. Light curves of PKS 1222+216 observed from December 2012 to March 2020 (MJD 56265–58914) at 15, 22, 43, and 86 GHz. The total CLEAN flux densities are shown for each frequency. In the third panel, the green circles represent the iMOGABA data, while the olive circles correspond to the BU data. The error bars used typical values of 10% for K and Q, and 20% for W. BU and MOJAVE were applied at 5% (Jorstad et al. 2017). The bottom panel shows the 0.1–100 GeV γ-ray light curve from the Fermi-LAT LCR. The upper limit of the flux density is indicated by the gray triangles. The gray vertical dashed line refers to the γ-ray peak at MJD 56975.5 (Ezhikode et al. 2022). The gray boxes are the periods that fit the exponential decay profiles at each frequency.

In the text

Fig. 4. Radio spectral indices at 22–43 GHz and 43–86 GHz. The top panel represents the values derived from the CLEAN total flux, while the bottom panel corresponds to the active flux. The orange scatters represent the spectral index at 22–43 GHz, while the green scatters indicate the spectral index at 43–86 GHz. The horizontal dashed line indicates that the spectral index is zero, and the dotted line in the top panel is the mean value with α = −0.6. In the bottom panel, the orange dotted line represents the average active spectral index of −0.19 at 22–43 GHz, while the green dotted line represents the average active spectral index of −0.37 at 43–86 GHz. The gray box indicates the decreasing period from 57124 to 57636 at 43 GHz.

In the text

	Fig. 5. Spectral index of PKS 1222+216 in the quiescent state. The black points represent the quiescent-state flux densities at 22, 43, and 86 GHz, with error bars indicating flux uncertainties. The orange slope corresponds to the quiescent spectral index between 22 and 43 GHz (αKQ = −0.81), while the green slope represents the quiescent spectral index between 43 and 86 GHz (αQW = −0.628).
In the text

Fig. 6. Upper panel: Radial distance of 43 GHz jet components relative to the core, which is shifted to the phase center at (0,0) position, from April 2015 to September 2016. The black arrows represent their position angles, and the solid lines represent their fitted trajectories. The colors correspond to the newly ejected components: C9 (blue), C10 (magenta), C11 (green), and C12 (cyan). The gray box highlights the interaction period between components C9 and A1. The horizontal yellow and pink lines indicate the mean positions of stationary components A1 and A2, respectively. Middle panel: Flux density of the core (red) and the newly ejected components, C9, C10, C11, and C12. Bottom panel: Flux density of components with relatively low flux density: C5, C6, C8, C12, A1, and A2. The black vertical line at MJD 56975 represents the epoch of the γ-ray flare.

In the text

	Fig. 7. Top panel: Power-law fit of the cooling timescales as a function of frequency (τν = A × να). Bottom panel: Power-law fit of the magnetic field strength as a function of frequency (Bν = B × νβ). The best-fit power-law indices are α = −0.56 ± 0.13 and β = 0.34 ± 0.04.
In the text

	Fig. C.1. Each panel presents the results of the power-law fitting of the size evolution of each component over time. The panels, from top to bottom, show the results for C9, C10, and C11. The sizes (FWHM values) are shown as scatter points with error bars, while the best power-law fits are represented by solid black lines.
In the text