---
abstract: 'We study the linear polarization of the radio cores of eight blazars simultaneously at 22, 43, and 86 GHz with observations obtained by the Korean VLBI Network (KVN) in three epochs between late 2016 and early 2017 in the frame of the Plasma-physics of Active Galactic Nuclei (PAGaN) project. We investigate the Faraday rotation measure (RM) of the cores; the RM is expected to increase with observing frequency if core positions depend on frequency due to synchrotron self-absorption. We find a systematic increase of RMs at higher observing frequencies in our targets. The RM–$\nu$ relations follow power-laws with indices distributed around 2, indicating conically expanding outflows serving as Faraday rotating media. Comparing our KVN data with contemporaneous optical polarization data from the Steward Observatory for a few sources, we find indication that the increase of RM with frequency saturates at frequencies of a few hundreds GHz. This suggests that blazar cores are physical structures rather than simple $\tau=1$ surfaces. A single region, e.g. a recollimation shock, might dominate the jet emission downstream of the jet launching region. We detect a sign change in the observed RMs of CTA 102 on a time scale of $\approx$1 month, which might be related to new superluminal components emerging from its core undergoing acceleration/deceleration and/or bending. We see indication for quasars having higher core RMs than BL Lac objects, which could be due to denser inflows/outflows in quasars.'
author:
- 'Jongho Park, Minchul Kam, Sascha Trippe, Sincheol Kang, Do-Young Byun, Dae-Won Kim, Juan-Carlos Algaba, Sang-Sung Lee, Guang-Yao Zhao, Motoki Kino, Naeun Shin, Kazuhiro Hada, Taeseok Lee, Junghwan Oh, Jeffrey A. Hodgson, and Bong Won Sohn'
title: 'Revealing the Nature of Blazar Radio Cores through Multi-Frequency Polarization Observations with the Korean VLBI Network'
---

=1

Introduction \[sect1\]
======================

Blazars, characterized by violent flux variability across the entire electromagnetic spectrum, are a sub-class of active galactic nuclei (AGNs) which show highly collimated, one-sided relativistic jets (see @UP1995 for a review). Large-scale magnetic fields which are strongly twisted in the inner part of the accretion disc or the black hole’s ergosphere play a crucial role in launching and powering of relativistic jets [@BZ1977; @BP1982]. Jets appear to be gradually accelerated and collimated magneto-hydrodynamically [@VK2004; @komissarov2007; @komissarov2009; @AN2012; @TT2013; @hada2013; @asada2014; @mertens2016; @hada2017; @walker2018] and they are directly linked to accretion process onto supermassive black holes [@marscher2002; @chatterjee2009; @chatterjee2011; @ghisellini2014; @PT2017]. Their parsec-scale radio morphology is characterized by (a) the ‘VLBI core’, a (radio) bright, optically thick, compact feature, and (b) an extended, optically thin, jet (e.g., @fromm2013).

The nature of the core is a matter of ongoing debate. The standard Blandford & Königl jet model describes the core as the upstream region where the conical jet becomes optically thin, i.e., at unity optical depth (e.g., @BK1979). In this scenario, the observed core position shifts closer to the physical location of the jet base at higher observing frequencies – the well-known ‘core shift effect’ [@lobanov1998]. Core shift has been observed in blazars (e.g., @OG2009b [@sokolovsky2011; @algaba2012; @pushkarev2012; @fromm2013; @hovatta2014]) as well as in nearby radio galaxies (e.g., @hada2011 [@martividal2011]), supporting the idea that the radio core marks the transition between the optically thick and thin jet regimes.

However, the plain ‘optical depth interpretation’ of the radio core ignores the physical structure of AGN jets. Especially, a standing conical shock, located at the end of the jet acceleration and collimation zone (e.g., @marscher2008), is expected (see also @PC2013a [@PC2013b] for a discussion of the transition region from a parabolic to a conical jet shape that is dominating jet synchrotron emission of blazars). Such a (quasi-)stationary feature – a ‘recollimation shock’ – may appear when there is a mismatch between the gas pressures in the jet and the confining medium (e.g., @DM1988 [@gomez1995; @gomez1997; @agudo2001; @mizuno2015; @marti2016]). Observations of the nearby radio galaxy M87 indeed reveal a stationary feature (known as HST-1) at the end of the jet collimation region [@AN2012], showing blazar-like activity such as rapid variability and high energy emission [@cheung2007]. In addition, recent studies have discovered that most $\gamma$-ray flares in blazars occur when new (apparently) superluminal jet components pass through the core (@JM2016, see also @ramakrishnan2014 [@casadio2015; @rani2015] for the case of individual sources and @jorstad2001 [@leontavares2011] for investigation of statistical significance between the two phenomena). This indicates that the core supplies the jet plasma electrons with large amounts of energy, with the possible formation of a shock as the source of high energy emission.

At first glance, these two models and corresponding observational support seem to be in contradiction. This conflict is resolved if the core consists of (a) a standing shock which is optically thin only at (sub)-mm wavelengths, plus (b) extended jet flows downstream of the shock. In this case, there is no core shift expected at millimeter wavelengths where the core becomes transparent. Interestingly, a recent study which used a bona-fide astrometric technique showed that the core shift between 22 and 43 GHz for BL Lacertae is significantly smaller than the expected one from lower frequency data, indicating that the core at these frequencies might be identified with a recollimation shock [@dodson2017]. However, a number of previous studies did not find such a trend at the same frequencies (e.g., @OG2009b [@algaba2012; @fromm2013]). This might be because the core position accuracy of previous VLBI observations is comparable to the expected amount of core shift at those frequencies.

[cccccccc]{}\[t!\]

& & 21.650 GHz & OJ 287 & 1.22 & 3C 279, OJ 287 & -28.79 & 1.04\
& & 86.600 GHz & 3C 84 & 2.25 & 3C 279, OJ 287 & 112.76 & 0.90\
& & 43.300 GHz & OJ 287 & 0.70 & 3C 279, OJ 287 & -19.52 & 2.37\
& & 86.600 GHz & OJ 287 & 1.86 & 3C 279, OJ 287 & 125.86 & 2.65\
& & 21.650 GHz & 3C 84 & 0.57 & 3C 279, OJ 287 & 11.53 & 0.94\
& & 86.600 GHz & OJ 287 & 1.25 & 3C 279, OJ 287 & 5.55 & 1.80\
& & 43.100 GHz & 3C 84 & 1.07 & 3C 279, OJ 287 & 54.29 & 0.33\
& & 129.300 GHz & $-$ & $-$ & $-$ & $-$ & $-$\
& & 21.650 GHz & 3C 84 & 1.04 & 3C 279, OJ 287 & -28.72 & 0.87\
& & 86.600 GHz & OJ 287 & 1.59 & 3C 279, OJ 287, 3C 454.3 & 32.42 & 1.04\
& & 43.100 GHz & OJ 287 & 0.90 & 3C 279, OJ 287, BL LAC & 9.07 & 1.37\
& & 129.300 GHz & $-$ & $-$ & $-$ & $-$ & $-$ \[Information\]

An alternative route is provided by multi-frequency polarimetric observations of the core that provide RMs, defined as $\rm EVPA_{obs} = EVPA_{int} + RM\lambda^2$, where $\rm EVPA_{obs}$ and $\rm EVPA_{int}$ are observed and intrinsic electric vector position angles (EVPAs) of linearly polarized emission and $\lambda$ is observing wavelength. If the core is the $\tau = 1$ surface of a continuous conical jet and the jet is in a state of energy equipartition, then the core RM obeys the relation $|{\rm RM_{\rm core, \nu}}| \propto \nu^a$, where $a$ is the power-law index of the electron density distribution given by $N_e \propto d^{-a}$, with $d$ being the distance from the jet base [@jorstad2007]. In this scenario, we observe polarized emission from regions closer to the jet base at higher frequencies due to the core shift effect, where one may expect higher particle densities and stronger magnetic fields. Looking at this argument the other way around, we would expect *no increase* in RM as a function of frequency at millimeter wavelengths *if* the core is indeed a standing recollimation shock. This provides the opportunity to uncover the nature of blazar VLBI cores, and thus the intrinsic structure of blazar jets, through multi-frequency polarimetric observations at millimeter wavelengths.

At centimeter wavelengths, many studies showed that the power law index $a$ is usually distributed around $a = 2$ (e.g., @OG2009a [@algaba2013; @kravchenko2017]), corresponding to a spherical or conical outflow [@jorstad2007]. A conical outflow is more likely than a spherical one because @pushkarev2017 showed that conical jet geometries are common in blazars. $a\approx2$ found in many blazars is in agreement with the fact that many blazars show core shift at these wavelengths. However, to the best of our knowledge, there are only a few studies of the core RM of blazars at (sub-)mm frequencies. [@jorstad2007] analyzed 7, 3, and 1 mm polarization data and obtained an average $\langle a \rangle = 1.8 \pm 0.5$ by comparing with other studies done at cm wavelengths. This result indicates that the dependence of RM on observing frequency might continue up to mm wavelengths. Some of their sources are not fitted well by $\lambda^2$ laws even at the highest frequencies, indicating that a frequency dependence of RM exists even at around 1 mm. Another study using the IRAM 30-m telescope at 3 and 1 mm found RMs (a few times $10^4 \rm\ rad/m^2$) that are much larger than those at cm wavelengths (a few hundred $\rm rad/m^2$, @hovatta2012), albeit within large errors (@agudo2014, see also @agudo2018a [@agudo2018b; @thum2018]).

A recent observation with the Atacama Large Millimeter Array (ALMA) at 1 mm has revealed a very high rotation measure of $(3.6\pm0.3)\times10^5\rm \ rad/m^2$ in 3C 273 with the core RM scaling with frequency like $\rm|RM|\propto\nu^{1.9\pm0.2}$ from cm to mm wavelengths. [@martividal2015] observed even larger RMs ($\approx 10^8 \rm\ rad/m^2$ in the rest frame) in the gravitationally lensed quasar PKS 1830–211 through ALMA observations at up to 300 GHz (about 1 THz in the rest frame). These results may suggest that (i) blazar core RMs rapidly increase as a function of frequency, as predicted by [@jorstad2007]; (ii) polarized (sub-)mm radiation might originate near the jet base, not from a recollimation shock (which presumably is located quite far from the jet base). However, it is uncertain whether this is a common behaviour of blazars or if these quasars are special. Therefore, a systematic study of blazar core RMs with multi-frequency polarimetric observations at (sub)mm wavelengths is necessary.

The Korean VLBI Network (KVN) has the unique capability of observing simultaneously at four frequencies, 22, 43, 86, and 129 GHz, or at two of these frequencies in dual polarization mode [@lee2011; @lee2014]. Thanks to the simultaneous observation at multiple frequencies, one can overcome rapid phase variations at high frequencies caused by tropospheric delay that reduce the coherence time by applying the fringe solutions obtained at lower frequencies to higher ones, i.e., frequency phase transfer (FPT) [@rioja2011; @rioja2014; @algaba2015; @zhao2018]. This technique increases the fringe detection rate to values larger than 80% even at 129 GHz for sources brighter than $\approx0.5$ Jy, making the KVN a powerful instrument for multi-frequency mm polarimetry of AGNs.

In early 2017, we launched a KVN large program, the Plasma-physics of Active Galactic Nuclei (PAGaN) project (see @kim2015 [@oh2015] for related studies), for monitoring about 14 AGNs at the four KVN frequencies in dual polarization mode almost every month.[^1] One of the main scientific goals of the project is a systematic study of RMs of blazars at mm wavelengths and their evolution in time. In this paper, we present the results from three observation epochs located between late 2016 and early 2017, which were performed as test observations for the initiation of the large program. We describe observations and data calibration in Section \[sect2\]. Results are shown and discussed in Section \[sect3\] and \[sect4\], respectively. In Section \[sect5\], we summarize our findings.

-100mm

[&gt;R &gt;R &gt;R &gt;R &gt;S&gt;S &gt;S &gt;S &gt;S &gt;S &gt;S &gt;S&gt;S&gt;S&gt;S&gt;S &gt;S ]{} & & & & & & & & & & & & &\
& & & & & & & $m [\%]$ & $\chi$ \[$^{\circ}$\] & $m$ & $\chi$ & $m$ & $\chi$ & & & (1) & (2) & (3) & & (4) & (5) & & & & & & & (6) & (7) & (8) & (9)

&&&&12/09-10/16&$-0.4\pm0.2$&$-0.3\pm0.2$&9.9$\pm$0.1&26$\pm$1&10.3$\pm$0.5&33$\pm$3&10.5$\pm$2.2&37$\pm$2&-0.05&(-2.0$\pm$1.0)$\times10^3$&(-4.0$\pm$4.5)$\times10^3$&1.0$\pm$1.8&&&&01/16-17/17&$-0.3\pm0.2$&$-0.7\pm0.2$&8.8$\pm$0.1&28$\pm$1&8.8$\pm$0.4&31$\pm$2&8.2$\pm$0.9&38$\pm$3&0.03&(-6.5$\pm$6.5)$\times10^2$&(-8.4$\pm$4.1)$\times10^3$&3.7$\pm$1.6&&&&03/22-24/17&$-0.1\pm0.2$&$-0.4\pm0.2$&8.6$\pm$0.1&30$\pm$1&11.2$\pm$0.2&36$\pm$2&12.4$\pm$0.4&49$\pm$3&-0.31&(-1.7$\pm$0.6)$\times10^3$&(-1.4$\pm$0.4)$\times10^4$&3.1$\pm$0.6&&&&12/09-10/16&$-0.2\pm0.2$&$-0.4\pm0.2$&4.6$\pm$0.1&-2$\pm$1&6.1$\pm$0.1&-1$\pm$2&7.4$\pm$0.5&-2$\pm$2&-0.38&(-3.8$\pm$5.6)$\times10^2$&(1.2$\pm$2.5)$\times10^3$&1.7$\pm$3.7&&&&01/16-17/17&$0.0\pm0.2$&$-0.7\pm0.2$&6.6$\pm$0.1&-26$\pm$1&8.3$\pm$0.1&-33$\pm$1&10.9$\pm$0.1&-36$\pm$2&-0.73&(1.5$\pm$0.2)$\times10^3$&(2.8$\pm$1.6)$\times10^3$&1.0$\pm$0.7[^2]&&&&03/22-24/17&$-0.2\pm0.2$&$-0.3\pm0.2$&5.0$\pm$0.1&-63$\pm$1&7.6$\pm$0.1&-63$\pm$1&8.3$\pm$0.2&-57$\pm$1&-0.43&(-0.3$\pm$3.6)$\times10^2$&(-4.8$\pm$1.5)$\times10^3$&7.1$\pm$14.6&&&&12/09-10/16&$0.3\pm0.2$&$0.3\pm0.2$&1.6$\pm$0.1&-83$\pm$3&1.5$\pm$0.2&-91$\pm$5&1.8$\pm$0.4&-117$\pm$11&0.02&(4.3$\pm$2.6)$\times10^3$&(5.2$\pm$2.4)$\times10^4$&3.6$\pm$1.1&&&&01/16-17/17&$0.4\pm0.2$&$0.1\pm0.2$&1.1$\pm$0.1&-63$\pm$2&1.4$\pm$0.2&0$\pm$3&1.4$\pm$0.1&36$\pm$8&-0.13&(-3.2$\pm$0.2)$\times10^4$&(-7.2$\pm$1.6)$\times10^4$&1.2$\pm$0.3&&&&12/09-10/16&$-0.3\pm0.2$&$-0.6\pm0.2$&2.6$\pm$0.2&-111$\pm$2&3.9$\pm$0.2&-74$\pm$3&6.6$\pm$0.8&-62$\pm$3&-0.65&(-1.2$\pm$0.1)$\times10^4$&(-1.4$\pm$0.5)$\times10^4$&0.3$\pm$0.5&&&&01/16-17/17&$-0.2\pm0.2$&$-0.7\pm0.2$&3.7$\pm$0.1&-88$\pm$1&5.8$\pm$0.2&-74$\pm$1&9.0$\pm$0.7&-61$\pm$2&-0.65&(-4.4$\pm$0.5)$\times10^3$&(-1.6$\pm$0.3)$\times10^4$&1.9$\pm$0.3&&&&12/09-10/16&$0.0\pm0.2$&$-0.2\pm0.2$&4.2$\pm$0.1&2$\pm$1&3.2$\pm$0.1&-1$\pm$3&2.2$\pm$0.5&-11$\pm$3&0.42&(4.5$\pm$6.2)$\times10^2$&(8.8$\pm$3.6)$\times10^3$&4.3$\pm$2.1&&&&03/22-24/17&$-0.1\pm0.2$&$-0.4\pm0.2$&4.6$\pm$0.2&-10$\pm$2&2.9$\pm$0.1&-22$\pm$2&2.5$\pm$0.8&-30$\pm$6&0.63&(2.6$\pm$0.5)$\times10^3$&(6.5$\pm$5.1)$\times10^3$&1.3$\pm$1.20235+164&FSRQ&0.94&25&03/22-24/17&$0.1\pm0.2$&$0.0\pm0.2$&2.8$\pm$0.1&48$\pm$1&3.3$\pm$0.1&51$\pm$2&2.4$\pm$0.3&59$\pm$3&-0.09&(-1.4$\pm$1.0)$\times10^3$&(-1.5$\pm$0.6)$\times10^4$&3.5$\pm$1.2BLLAC&BLO&0.07&4&03/22-24/17&$0.0\pm0.2$&$-0.9\pm0.2$&3.5$\pm$0.1&38$\pm$1&4.1$\pm$0.3&31$\pm$2&3.8$\pm$0.4&19$\pm$3&-0.11&(1.0$\pm$0.3)$\times10^3$&(6.4$\pm$1.9)$\times10^3$&2.7$\pm$0.61633+38&FSRQ&1.81&27&03/22-24/17&$-0.2\pm0.2$&$-0.2\pm0.2$&1.5$\pm$0.2&31$\pm$3&1.4$\pm$0.2&24$\pm$4&2.2$\pm$0.7&18$\pm$6&-0.09&(6.3$\pm$4.8)$\times10^3$&(2.4$\pm$2.7)$\times10^4$&1.9$\pm$2.0

Observations and Data Reduction {#sect2}
===============================

We observed a total of 11 sources in the 22, 43, and 86 GHz bands with the KVN on 2016 December 9–10 and in the four bands including the 129 GHz band on 2017 January 16–17 and 2017 March 22–24 with observation time of $\approx48$ hours for each epoch. Since KVN can observe at two frequencies simultaneously in dual polarization mode, we allocated the first half of the observing time to 22/86 GHz observations and the other half to 43/129 GHz. Although we obtained the data at 129 GHz in the two epochs observations, we had a difficulty in polarization calibration of the data and thus we did not include them in this paper. More sophisticated investigation of the 129 GHz data will be presented in a forthcoming paper (Kam et al. in preparation). All sources were observed in 6–15 scans of 5–20 minutes in length, depending on source declination and brightness. We performed cross-scan observations at least twice per hour to correct antenna pointing offsets that might lead to inaccurate correlated amplitudes. The received signals were 2-bit quantized and divided into 4 sub-bands (IFs) of 16 MHz bandwidth for each polarization and each frequency. Mark 5B recorders were used at recording rates of 1024 Mbps. The data were correlated with the DiFX software correlator in the Korea-Japan Correlator Center [@lee2015a]. Table \[Information\] summarizes our observations.

A standard data post-correlation process was performed with the NRAO Astronomical Image Processing System (AIPS). Potential effects of digital sampling on the amplitudes of cross-correlation spectra were estimated by the AIPS task ACCOR. Amplitude calibration was done by using the antennas gain curves and opacity corrected system temperatures provided by the observatory. The fringe amplitudes were re-normalized by taking into account potential amplitude distortion due to quantization, and the quantization and re-quantization losses [@lee2015b].

The instrumental delay residuals were removed by using the data in a short time range of bright calibrators, either 3C 279 or 3C 454.3. To apply the FPT technique, a global fringe fitting was performed with a solution interval of 10 seconds for the lower frequency first (22 or 43 GHz), which led us to very high fringe detection rates $\gtrsim 95\%$ in most cases. Then, we transferred the obtained fringe solutions to the simultaneously observed higher frequency (86 GHz). This process corrects rapidly varying tropospheric errors in the visibility phases at high frequencies (though not the ionospheric errors that vary more slowly). Then, the residual phases have much longer coherence times, typically larger than a few minutes. Thus, we performed a global fringe fitting with a much longer solution interval of $\approx3$ minutes for the high frequency data, which resulted in quite high fringe detection rates – usually larger than 95% at 86 GHz for our sources. Bandpass calibration was performed by using scans on bright sources such as 3C 279.

The cross-hand R-L phase and delay offsets were calibrated by using the data for bright sources, such as OJ 287, 3C 84, 3C 279 and 3C 454.3, located within short time ranges, with the task RLDLY. We used the Caltech Difmap package for imaging and phase self-calibration [@shepherd1997]. Typical beam sizes are $5.6\times3.2$, $2.8\times1.6$, and $1.4\times0.8$ mas at 22, 43, and 86 GHz, respectively. We determined the feed polarization leakage (D-terms) for each antenna by using the task LPCAL [@leppanen1995] with a total intensity model of the D-Term calibrators. 3C 84 usually serves as a good D-Term calibrator thanks to its high flux density and very low degree of linear polarization ($\lesssim$0.5%) at lower frequencies but less so at high frequency (86 GHz) where its linear polarization becomes non-negligible (Kam et al., in preparation). Thus, we also used a compact, bright, and polarized source OJ 287 at 86 GHz. We chose the D-Term calibrator for each epoch and for each frequency by comparing the behaviour of observed visibility ratios on the complex plane with the D-Term models of different calibrators (see Appendix \[appendixa\]). The EVPA calibration was performed by comparing the integrated EVPAs of the VLBI maps of the EVPA calibrators after the instrumental polarization calibration with contemporaneous KVN single dish polarization observations. We performed KVN single dish observations within two days of each VLBI observations as described in [@kang2015]. For the 2016 data, we have two 86 GHz data separated by 1 day. We note that the maps for all sources after the calibration are almost identical to each other and we used the average of Stokes I, Q, and U maps of the two data for our further analysis.

Estimating errors for degree of linear polarization ($m$) and EVPA is important but not straightforward. Errors for each polarization quantity can be derived from the following relations: $$\sigma_p = \frac{\sigma_Q + \sigma_U}{2}$$ $$\sigma_{\rm EVPA} = \frac{\sigma_p}{2p}$$ $$\sigma_{m} = \frac{\sigma_p}{I}$$ where $\sigma_Q$ and $\sigma_U$ denote rms noise in the Stokes Q and U images, respectively, $p = \sqrt{Q^2 + U^2}$, and $m = p/I$ [@hovatta2012]. In most cases, random errors are quite small and systematic errors are much more dominant in the above quantities. Imperfect D-term calibration is usually the most dominant source of errors in $m$. For EVPAs, both the D-term uncertainty and EVPA correction error are important. Following [@roberts1994], errors of $m$ and EVPA caused by residual D-terms can be expressed as: $$\label{eq4}
\sigma_{m, {\rm D}} = \sigma_D(N_aN_{\rm IF}N_s)^{-1/2}$$ $$\label{eq5}
\sigma_{\rm EVPA, {\rm D}} \approx \frac{\sigma_{m, {\rm D}}}{2m}$$ where $\sigma_D$ is the D-term error, $N_a$ and $N_{\rm IF}$ are the number of antennas and IFs, respectively, and $N_s$ is the number of scans having independent parallactic angles which depend on the source declination. We estimated the D-term errors by comparing the D-terms obtained from different D-Term calibrators (see (5) in Table \[Information\] and Appendix \[appendixa\]) and estimated $\sigma_{m, {\rm D}}$ and $\sigma_{\rm EVPA, {\rm D}}$ using Equation \[eq4\] and \[eq5\]. Thanks to the number of IFs being four and the large parallactic angle coverage of our sources, we could achieve errors in $m$ (typically $0.1-0.3\%$) much smaller than the D-Term errors (typically $1-2\%$). We also assessed the EVPA correction error, $\sigma_{\Delta\chi}$, by comparing the amount of EVPA rotation calculated from different EVPA calibrators (see (6) and (8) in Table \[Information\]). Then, we added $\sigma_m$ and $\sigma_{m, {\rm D}}$ quadratically for $m$ and $\sigma_{\rm EVPA}$, $\sigma_{\rm EVPA, D}$, and $\sigma_{\Delta\chi}$ quadratically for EVPA.

In the Appendix, we show the results of D-term calibration and the temporal evolution of the D-terms. The overall D-terms are usually less than $\approx10\%$, except for Ulsan station at 86 GHz which showed D-terms as large as $\approx20\%$. The D-terms obtained from different calibrators are quite consistent with each other, showing standard deviations of $\lesssim2\%$ (see (5) in Table \[Information\]). The D-terms are more or less stable over $\approx$ 3 months, showing standard deviations of $\lesssim2\%$. We also compare our KVN 22/43 and 86 GHz data of 3C 273 observed in 2016 December with contemporaneous Very Long Baseline Array (VLBA) 15/43 GHz data, respectively. Both fractional polarization and EVPAs at a few different locations in the jet are in good agreement within errors between the data of the different instruments, considering non-negligible time gaps between the observations and the expected RM of a few hundred $\rm rad/m^2$ in the jet (e.g., @hovatta2012).

\
\
\
\
\

results {#sect3}
=======

RM at radio wavelengths
-----------------------

In Figure \[result\], we present polarization maps for multiple frequencies convolved with the KVN 22 GHz beam size (left panels), EVPAs at the core as a function of $\lambda^2$ (central panels), and RMs as a function of geometric mean observing frequency (right panels). We obtained one RM value from each adjacent data pair in the EVPA–$\lambda^2$ plots because we could not obtain good $\lambda^2$ fits across the three bands in most cases. We have to rotate the 22 and 43 GHz EVPAs by more than $720^{\circ}$ and $180^{\circ}$ to explain this behavior with $n\pi$ ambiguity, which translates into $\rm RM \gtrsim 10^5\ rad/m^2$. This RM value is too large especially at 22 GHz because almost all of the detected core RM of blazars are less than $1000\rm\ rad/m^2$ at $\lesssim15$ GHz [@hovatta2012]. Alternatively, different optical depths of the cores at different frequencies might be responsible for the non-$\lambda^2$ fits. Especially, EVPA rotations by $90^\circ$ are expected in case of a transition of the core from optically thick to optically thin [@pacholzcyzk1970]. We provide the spectral index, $\alpha$ in $S_{\nu} \propto \nu^{\alpha}$, between adjacent frequencies in columns (3) and (4) in Table \[result\] and in the center panels of Figure \[result\] (blue asterisks). We found that values of $\alpha$ measured at different frequency pairs differ from each other by more than $1\sigma$ in only four cases, 3C 279 in 2017 January, OJ 287 in 2017 January, 3C 345 in 2017 January, and BL Lac in 2017 March. When we rotate the EVPAs at the lowest frequency by $90^\circ$ for these cases, we are left with even worse $\lambda^2$ fits (higher $\chi^2$ values) compared to the case of no rotation. In addition, the degree of polarization is not much different, less than a factor of $\approx2$, at different frequencies (Section \[fracpol\]), while it should decrease by a factor of $\approx7$ if there were a $90^{\circ}$ flip [@pacholzcyzk1970]. Therefore, the progressively steeper EVPA rotations at higher frequencies are more likely to be due to the core shift effect, as shown in the numerical simulations of special relativistic magnetohydrodynamic jets [@porth2011].

We identified the origin (0,0) of the maps with the location of the cores. This might not be always the case exactly; however, the beam size is quite large and thus the effect of an offset of the core position from the origin would be insignificant. We summarize our sources’ basic information and their observed polarization quantities in Table \[table:source\]. We note that we did not consider the effect of the integrated (Galactic) RM because it is a few hundreds $\rm rad/m^2$ at most (e.g., @taylor2009); this is much smaller than the typical RM error we obtain, and at the KVN frequencies the EVPA rotation due to the integrated RM is negligible.

Since the beam sizes are quite different at different frequencies, we need to quantify the errors in polarization quantities introduced by the convolution of all maps with the 22 GHz beam. We compared $m$ and EVPA at the core found with and without using convolution for each source at each frequency and added the differences quadratically to the errors of $m$ and EVPA found from convolved maps, respectively. Using contemporaneous BU VLBA maps, We identified and excluded sources which have complex polarization structure near the core that cannot be resolved with the KVN; this leaves us with eight sources. We briefly describe the results for individual sources below.

### 3C 279

This source is characterized by longitudinal (i.e., parallel to the jet direction) EVPAs that show a smooth distribution from the core to the inner jet[^3] (e.g., @jorstad2005). Similarly, our KVN maps show basically longitudinal EVPAs but rotated by up to $\approx20^{\circ}$ as function of frequency. RMs between adjacent frequency pairs range from $\approx10^3$ to $\approx10^4\rm\ rad/m^2$. We fitted a power law function to the RMs as a function of geometrical mean observing frequency and obtained the power law index $a$ in the relation ${\rm|RM|} \propto \nu^a$. Since we calculate each RM value from only two data points, the RM errors are relatively large, which results in relatively large errors in $a$. However, $a$ values in all three epochs show a good agreement with $a=2$–3, which is quite consistent with the results of previous KVN single-dish polarization monitoring of 3C 279 [@kang2015]. We note that the core EVPAs of this source might be contaminated by polarization from the inner jet. However, EVPA rotation of the inner jet region is expected to be very small at $\gtrsim 22\rm\ GHz$ because of a relatively small RM in that region ($\lesssim250\rm\ rad/m^2$; @hovatta2012). Therefore, we conclude that the observed EVPA rotation over frequency of this source is dominated by the core polarization.

### OJ 287\[oj287\]

For this source, we could find contemporaneous VLBA data at 15 GHz from the MOJAVE program[^4], observed on 2017 January 28 and March 11, with our KVN data being obtained on 2017 January 17–18 and March 22–24, respectively. We included those data in our analysis after convolving the 15 GHz maps with the KVN 22 GHz beam (because the KVN beam is larger than the VLBA one even though it is at a higher frequency.) This source has shown slow and gradual EVPA variation in time and thus a potential variability during the time gap between the MOJAVE and our KVN observations ($\lesssim2$ weeks) would not be significant (see the AGN monitoring database of the 26-meter University of Michigan Radio Astronomy Observatory[^5]). In addition, OJ 287 has been known for relatively small RMs at cm wavelengths (e.g., @hovatta2012). Our KVN maps are consistent with zero RM (within errors) in 2016 December. However, in the 2017 January data, EVPA rotations from 15 to 86 GHz in the same direction being steeper at higher frequency pairs are observed, which results in $a = 0.98\pm0.66$. In the 2017 March data, the EVPA rotations between 15 and 43 GHz were almost zero within errors but a relatively large rotation with $\rm|RM|\approx5\times10^3\rm rad/m^2$ was detected at 43/86 GHz.

### CTA 102\[cta102\]

This source shows a high degree of polarization, up to $\approx 40\%$, in its jet at 22 GHz but becomes compact at higher frequencies. The EVPAs rotate rapidly as function of frequency, with different slopes in the EVPA–$\lambda^2$ diagram between 22/43 GHz and 43/86 GHz. RMs at 43/86 GHz are a few times $10^4\rm\ rad/m^2$ in the source rest frame in both epochs. However, $a$ value in the 2016 December data is much larger than that in the 2017 January data. The signs of RMs are different in our two epochs, while their absolute values are of the same order of magnitude. This behavior suggests that the line-of-sight component of the jet magnetic field changed its direction within $\approx1$ month, while magnetic field strength and electron density (or at least their product) did not vary substantially. This sign change might be related to a strong flare that occurred during our KVN observations [@raiteri2017]. We discuss possible reasons for the sign flip in CTA 102 in Section \[signchange\].

### 3C 345 {#3c345}

This source shows almost the same RMs at 22/43 and 43/86 GHz, $\rm |RM|\approx10^4\ rad/m^2$ in the 2016 December data, with $a$ being consistent with zero within errors. However, the RM at 43/86 GHz is much larger than that at 22/43 GHz about one month later, resulting in $a=1.86\pm0.3$. These results indicate that there is substantial time variability in this source. Similarly to the case of CTA 102, this source shows a flare during our KVN observations and the flux density in early 2017 is almost three times higher than that in mid 2016 at 1 mm.[^6]

### 1749+096

This source displays a rather compact jet geometry at all frequencies and in both epochs (2016 December and 2017 March). Interestingly, the degree of linear polarization is larger at lower frequencies, which is not usually seen in other sources (Section \[fracpol\]). This might suggest that we are looking at a mixture of different polarization components with different EVPAs and/or RMs or that internal Faraday rotation occurs in this source (Section \[sectscreen\]). The values of RM range from $\approx10^3$ to $\approx10^4\rm\ rad/m^2$ in both epochs data. The values of $a$ are consistent within errors. Therefore, the sign, absolute value, and frequency dependence of RM appear to be stable over $\approx3$ months.

### 0235+164

This source shows a rather compact jet geometry. A systematic rotation of EVPAs in the same sense as function of frequency can be seen from 22 to 86 GHz. RMs range from $\approx10^3$ to $\approx2\times10^4\rm\ rad/m^2$, with a substantially larger RM at a higher frequency pair, resulting in $a = 3.47\pm1.18$.

### BL LAC

The EVPA–$\lambda^2$ diagram shows a steeper slope between 43 and 86 GHz than between 22 and 43 GHz, providing $a = 2.65\pm 0.61$. RM values range from $\approx10^3$ to $\approx6\times10^3\rm\ rad/m^2$. This is consistent with previous measurements of the core RM between 15 and 43 GHz [@OG2009a; @gomez2016]. However, the high resolution *RadioAstron* space VLBI image shows a complex RM structure in the core region including a sign change, indicating the presence of helical magnetic fields there [@gomez2016]. Therefore, we may be looking at a blend of those structures in our KVN images (see Section \[multipleshocks\]).

### 1633+38

This source is relatively faint and shows a low degree of linear polarization ($\approx2\%$) which leads to relatively large EVPA errors. Thus, the RM at 43/86 GHz is comparable to its error and the obtained $a=1.91\pm1.96$ has also a large error. When fitting a single linear function to the EVPAs at the three frequencies available, we obtained $\rm RM = 974\pm509\ rad/m^2$. This is surprisingly low because a very high RM, $\approx2.2\times10^4\rm\ rad/m^2$, was reported previously for the core of this source using six different frequencies from 12 GHz to optical wavelength [@algaba2012]. We found that EVPAs at the core of this source in the MOJAVE program are $44^\circ$ and $32^\circ$ in 2016 November and 2017 April, respectively. If there is no substantial EVPA variability in between these two epochs, then a simple $\lambda^2$ fit can explain the data from 15 to 86 GHz, suggesting that there is no $n\pi$ ambiguity in our data and that the core RM of this source is indeed quite small. The observations in [@algaba2012] were performed in 2008 November which indicates that there is substantial temporal variability of the EVPA rotation in 1633+38. Four years before their observations, core EVPAs could not be fitted with a single $\lambda^2$ law and the obtained RM value was much smaller, $\rm RM = -570\pm430\ rad/m^2$ [@algaba2011]. This also suggests substantial temporal RM variability in 1633+38.

Optical EVPAs from the Steward observatory\[optdata\]
-----------------------------------------------------

For a few sources, we obtained quasi-contemporaneous (taken within $\lesssim$1 week) optical polarization data from the Steward Observatory blazar monitoring program[^7] (see @smith2009 for details) for some epochs. We summarize the optical data we used in Table \[Steward\]. (We excluded some additional datasets due to their large errors.) The optical polarimetry errors are usually quite small, $\lesssim1^{\circ}$, unless sources are very weakly polarized. However, optical polarization of blazars often show rapid variability on short time scales (e.g., @jorstad2013) presumably due to a smaller size of the emission region at higher frequencies [@marscher1996]. In order to take into account the uncertainty arising from the time gap between optical and radio observations we estimated errors from source variability as follows.

The Steward Observatory blazar monitoring program usually observes each source multiple times for a specified period spanning a few days in a broad optical band from 500 to 700 nm. We selected all data in the periods that are close to our KVN observations and used the mean and standard deviation of the data points as representatives of optical EVPA and typical error, respectively. We assumed that their optical EVPAs show random-walk type variations with time. (In addition to statistical variability, many blazars occasionally show smooth, systematic optical EVPA rotations that might be associated with high energy flares; e.g., @blinov2015). We multiplied the observed standard deviation by the square root of the ratio of the time gap between optical and radio observations to the duration of the period in which a set of optical data was obtained. Under this assumption, the longer the time separation, the larger the formal uncertainty of an optical EVPA at the time of the corresponding radio observation.

[clcl]{}\[t!\]

& 11/30/16–12/02/16 & $-216.5 \pm 31.4$ & $244^{+701}_{-93}$\
& 01/10/17–01/13/17 & $111.7 \pm 12.6$ & $591^{+1467}_{-276}$\
3C 279 & 01/10/17–01/14/17 & $85.2 \pm 2.7$ & $209^{+697}_{-71}$\
BL LAC & 03/29/17–03/30/17 & $-8.9 \pm 26.1$ & $138^{+298}_{-65}$

Due to lack of frequency coverage between 86 GHz and optical wavelengths, we suffered from potential $n\pi$ rotation of the optical EVPAs. Therefore, we assumed that the optical EVPAs of our sources rotate in the same direction as the ones at mm wavelengths and that the EVPA rotation between 86 GHz and the optical band does not exceed $\pi$. We present the optical EVPA values of three sources in Table \[Steward\] and plot them with the core EVPAs from our KVN observations in the left panel of Figure \[paopt\]. We also show the RMs obtained from each adjacent frequency pair in the RM–frequency diagram (the right panels of Figure \[paopt\]). Our assumption on $n\pi$ rotation appears reasonable because the optical EVPAs follow the trend of EVPA rotation established at radio frequencies, although we cannot rule out the possibility of coincidence because of the low number of sample. The RMs obtained from the EVPA difference between 86 GHz and the optical frequencies are about an order of magnitude higher than the values obtained from the frequency pair 43/86 GHz. We note that the observed RMs between 86 GHz and optical light exceed the minimum possible measurable RM by an order of magnitude except for BL Lac for which the observed RM is about two times the minimum measurable RM.

The power-law increase of RM as a function of frequency does not continue to optical wavelengths but saturates at a certain frequency (right panel of Figure \[paopt\]). We used the term *transition frequency*, $\nu_{\rm trans}$, for this frequency. We calculated asymmetric errors on $\nu_{\rm trans}$ via Monte-Carlo simulations by adding Gaussian random numbers to the best-fit parameters of the radio RM–$\nu$ power-law relation with standard deviations identical to their $1\sigma$ errors. The obtained $\nu_{\rm trans}$ are distributed from 138 to 591 GHz in the source rest frame for different sources and in different epochs (Table \[Steward\]). We note that $\nu_{\rm trans}$ for BL Lac is consistent with the observed frequency of 86 GHz within $1\sigma$ because of the relatively large minimum measurable RM.

fractional polarization\[fracpol\]
----------------------------------

We present the degree of linear polarization $m$ as function of $\lambda$ in Figure \[frac\]. Various de-polarization models are available to explain the evolution of $m$ with wavelength (see @osullivan2012 [@farnes2014] for summaries). In principle, $m-\lambda$ scalings can be used to determine whether the emitting region and the Faraday screen are co-spatial or not, whether magnetic fields in the screen are regular or turbulent, whether there are multiple components with different polarization properties on scales smaller than the spatial resolution, and so on [@burn1966; @conway1974; @tribble1991; @sokoloff1998]. However, we did not try to apply those models to our data because (i) our data provide sparse frequency sampling over a limited frequency range, (ii) the models are mostly appropriate for optically thin emitters while we are dealing with (partially) optically thick cores, and (iii) different observing frequencies might probe different physical regions, as suggested by the complicated $\chi-\lambda^2$ scalings of the EVPAs. Instead, we obtain a polarization spectral index $\beta$ by fitting $m\propto\lambda^{\beta}$ to our data (Table \[result\], see @farnes2014), which could be used for future theoretical studies (e.g., @porth2011) and for comparison with observations at lower frequencies (e.g., @farnes2014). We refer the readers to detailed studies of degree of linear polarization at different wavelengths of AGNs using broadband radio spectro-polarimetric observations (e.g., @osullivan2012 [@osullivan2017; @hovatta2018; @pasetto2018]) and investigating spatially resolved optically thin emitting regions with multi-frequency VLBI observations (e.g., @hovatta2012 [@kravchenko2017]). The median, mean, and standard deviation of $\beta$ are -0.11, -0.17, and 0.38, respectively. All sources show $\beta\lesssim0$ except for 1749+096 which showed $\beta\approx0.5$ in both epochs (2016 December and 2017 March).

discussion {#sect4}
==========

In this section, we interpret the results of the core polarization properties of eight blazars, five flat spectrum radio quasars (FSRQs) and three BL Lac objects (BLOs).

RM distributions at different frequencies {#rmdist}
-----------------------------------------

We present the distributions of the absolute RM values for different frequency pairs in the left panel of Figure \[rmhisto\]. We excluded RM values whose absolute values are smaller than their $1\sigma$ errors. The histograms show that the median RM increases with frequency. We note that the minimum possible measurable RMs are 440 and 2038 $\rm rad/m^2$ at 22/43 and 43/86 GHz, respectively, assuming a typical EVPA error of $2^\circ$, $3^\circ$, and $3^\circ$ at 22, 43, and 86 GHz, respectively. As is evident in Figure \[rmhisto\], the RM values we found are much larger than these minimum possible measurable RMs. Notably, the trend of increasing RMs with increasing observing frequencies cannot be produced artificially.

We collected the median core RMs at cm wavelengths for our sources from [@hovatta2012] and show all median RM values as function of frequency in the right panel of Figure \[rmhisto\]. As expected, RMs increase with increasing frequency (355, 2620, and 14200 $\rm rad/m^2$ for 8.1–15.4, 22–43, and 43–86 GHz, respectively). Un-weighted fitting of a power-law function returns a best-fit power-law index $a=2.42$. Although the sample size is small and the standard deviations of the RM distributions are large, the obtained power-law index is quite close to $a = 2$, indicating that Faraday rotating media of blazars core can be represented as conical outflows statistically [@jorstad2007].

Instead of comparing RM distributions of all sources at different frequencies, we collected the power-law indices $a$ obtained for each source in Figure \[aplot\]. We have 13 measurements in total, with some sources having more than one measurement. The mean and standard deviation of all $a$ values are $a = 2.25 \pm 1.28$, which is consistent with $a=2$ and the fitting results for the median values of RM distributions at different frequencies. Our results are also consistent with previous studies of blazars at both cm and mm wavelengths (e.g., @jorstad2007 [@OG2009a; @algaba2013; @kravchenko2017; @hovatta2018]). However, many $a$ values are located far from the mean value, which potentially indicates a bimodal distribution. Assuming a power-law electron density distribution as function of jet distance ($d$), $N_e \propto d^{-a}$, toroidal magnetic fields dominant in the Faraday screen, $B\propto d^{-1}$, a conical geometry of the Faraday screen, $dl\propto d$ and $a=2$, and energy equipartition, $d_{\rm core} \propto \nu^{-1}$, one obtains ${\rm RM_{core}} \propto \int N_eBdl \propto \nu^2$ [@jorstad2007]. If some of these assumptions are not satisfied, one might expect deviations from a $a=2$ scaling. For example, there is growing evidence for a parabolic geometry of the blazar cores (e.g., @algaba2017 [@pushkarev2017]). In some cases, a need for helical magnetic fields instead of dominant toroidal fields in blazars has been pointed out (e.g., @zamaninasab2013). Likewise, the assumption of energy equipartition between radiating particles and magnetic fields may not hold for some sources (e.g., @homan2006). However, it is difficult to determine accurate $a$ values for each source with the current data only due to source variability and relatively large errors in $a$. For example, the values of $a$ are likely related with source’s flaring activity as seen in the case of CTA 102 and 3C 345 (see Section \[cta102\] and  \[3c345\]). We expect that our monthly monitoring program will allow us to investigate the reason for potential difference in $a$ values for different sources and in different epochs, together with the detailed information of the compact core region provided by the ultra high resolution arrays such as the Global Millimeter-Very Long Baseline Interferometry Array (GMVA, see e.g., @kim2016) or the Radioastron space VLBI (e.g., @gomez2016).

Change of core opacities from optically thick to thin\[transition\]
-------------------------------------------------------------------

In Section \[optdata\], we show that the power-law increase of RM as a function of frequency might not continue to optical wavelengths but flatten out at a certain frequency, $\nu_{\rm trans}$. We suggest that the cores of blazars become fully transparent at $\nu > \nu_{\rm trans}$, meaning no core-shift and thus no more frequency dependence of RM at those frequencies. Accordingly, the radio core may be a standing recollimation shock at $\nu > \nu_{\rm trans}$. For CTA 102, $\nu_{\rm trans}$ increased substantially from $\approx240$ GHz to $\approx590$ GHz within one month, albeit within large errors (Table \[Steward\]). This might be related to a strong flare that occurred at the time of our observations (see Section \[cta102\]) which ejected a large amount of relativistic electrons into the core, causing it to become optically thick. We obtained $\nu_{\rm trans} \approx 210$ and 140 GHz for 3C 279 and BL Lac, respectively. This result seems to be in line with recent astrometric observations of BL Lac which found a systematic deviation of the amount of core-shift from the one expected for a Blandford–Königl type jet at 22/43 GHz [@dodson2017]. Likewise, the scaling of synchrotron cooling time with frequency in BL Lac matches a standing shock better than an optically thick jet [@kim2017]. In summary, one may expect no frequency dependence of RM and no core-shift above $\approx140$ GHz for BL Lac and above $\approx210$ GHz for 3C 279 and CTA 102. However, we stress that the conclusions presented in this section are valid only when the assumptions of (a) no $n\pi$ ambiguity and (b) EVPA rotations in the same sense from mm to optical hold. We will study core opacity evolution and RM saturation further with dedicated upcoming multi-frequency observations at mm and sub-mm, combining data from KVN and ALMA (Park et al., in preparation).

The Faraday screen\[sectscreen\]
--------------------------------

Identifying the source of Faraday rotation is very difficult. As the amount of Faraday rotation is inversely proportional to the square of the mass of charged particles, thermal electrons and/or low-energy end of radiating non-thermal electrons would be the dominant source of the observed RM. If those Faraday rotating electrons are mixed in with the emitting plasma in jets, internal Faraday rotation occurs. However, if the rotating medium is located outside the jet, e.g., in a sheath surrounding the jet or the broad/narrow line regions (BLRs/NLRs), then the observed Faraday rotation is external to the jet. Theoretical models assuming an optically thin jet with spherical or slab geometries showed that it is very difficult for internal Faraday rotation to cause EVPA rotations larger than $45^{\circ}$ without severe depolarization [@burn1966]. Multiple studies showed that many blazars indeed have EVPA rotations larger than $45^{\circ}$ without significant depolarization, indicating that the source of Faraday rotation is external to the jets usually (e.g., @ZT2003 [@ZT2004; @jorstad2007; @OG2009a; @hovatta2012]). A sheath surrounding the jet is considered to be the most viable candidate for an external Faraday rotating medium; in addition, BLRs/NLRs are unlikely sources of RM given the time variability of RMs in jets and volume filling factor arguments [@ZT2002; @ZT2004; @hovatta2012]. Nevertheless, there is indication for potential internal Faraday rotation in some sources [@hovatta2012].

We cannot identify the Faraday screen from our data because of their limitations (Section \[fracpol\]). Nevertheless, we note that the observed RM–frequency relations having $a\approx2$ (Section \[rmdist\]) and the polarization spectral indices being predominantly negative ($\beta\lesssim0$) for our sources support the conclusion of previous studies that an external jet sheath acts as Faraday screen [@ZT2002; @ZT2004; @hovatta2012]. However, for 1749+096 we observed the degree of fractional polarization at high frequencies to be smaller than the one at lower frequencies, with $\beta\approx0.5$ in both epochs (2016 December 9 and 2017 March 24) – which is almost impossible to explain with any standard external depolarization model [@hovatta2012]. Such ‘inverse depolarization’ can be due to blending of different polarized inner jet components at different frequencies [@conway1974] or internal Faraday rotation in helical or loosely tangled random magnetic field configurations [@homan2012].

RM sign change {#signchange}
--------------

We observed a RM sign change for CTA 102 within $\approx$1 month, while the absolute values of RM did not change much (Figure \[result\]). Previous studies found temporal sign reversals in RMs for other sources (e.g., @mahmud2009 [@lico2017]), sign reversals between core and jet (e.g., @mahmud2013), and sign reversals in the cores at different frequencies intervals (e.g., @OG2009a). Scenarios proposed to explain such RM sign changes include: (i) a reversal of the magnetic pole of the black hole facing the Earth; (ii) torsional oscillations of the jet; (iii) a ‘nested-helix’ magnetic field structure; and (iv) helical magnetic fields in jets seen at different orientations due to relativistic abberation, depending on whether $\theta\Gamma$ is larger or smaller than 1, where $\theta$ is the viewing angle and $\Gamma$ is the bulk Lorentz factor of jets (see @mahmud2009 [@mahmud2013] for (i)-(iii) and @OG2009a for (iv) for details). Although all scenarios are possible theoretically, we focus on the fact that CTA 102 underwent a relatively strong flare in the period of our KVN observations (Section \[cta102\]).

Evidence for the presence of helical magnetic fields in AGN jets has been provided by many studies, starting with the detection of a transverse RM gradient in the jet of 3C 273 [@asada2002] which was later confirmed by other studies [@ZT2005; @hovatta2012]. Similar behaviour has been found in many BL Lac objects (e.g., @gabuzda2004 [@gabuzda2015]), radio galaxies (e.g., @kharb2009), and quasars (e.g., @asada2008 [@algaba2013; @gabuzda2015]). Furthermore, general relativistic magnetohydrodynamic simulations of AGN jets showed that the combination of the rotation of the jet base and the outflow leads to the generation of a helical field and associated Faraday rotation gradients [@BM2010]. If helical magnetic fields pervade in the jet sheaths and if they are the main contributor of the observed RMs as speculated in Section \[sectscreen\]), the sign of RMs would be determined by whether $\theta\Gamma$ is larger or smaller than 1, as explained in [@OG2009a].

As noted in Section \[cta102\], a strong flare at multiple wavelengths occurred during our KVN observations [@raiteri2017]. Flares in blazars are usually associated with new VLBI components emerging from the cores (e.g., @savolainen2002). The flare in CTA 102 would then likewise be connected to newly ejected VLBI components. [@jorstad2005] found $\theta = 2.6^{\circ}$ and $\Gamma=17.2$ for CTA 102, which yields $\theta\Gamma = 0.78$. If there is bending and/or acceleration or deceleration of the ejected component, changes of $\theta\Gamma$ across the value $\theta\Gamma=1$ can occur and the sign of RM reverses. Assuming scenario (i) or (ii) as mechanism behind the sign reversal requires coincidence with the recent strong flaring activity of this source. Furthermore, scenario (iii) can be related to flaring since a new jet component might lead to temporal increase of the relative contribution of the inner field to the outer field in the magnetic tower model. However, in this case it is difficult to explain the observation of similar RM magnitudes in the two epochs (a few times $10^4\rm\ rad/m^2$ for CTA 102); the relative contributions by the inner and outer magnetic fields to the observed RMs must be almost exactly opposite in different epochs, which, again, would be a coincidence (but see @lico2017 for the case of Mrk 421 which supports this scenario). Therefore, we conclude that scenario (iv) provides the most natural way to explain the observed sign change in RMs of CTA 102 as it does not require substantial changes in the physical properties of the jets. Our interpretation is also consistent with modelling the multi-wavelength flare in this source in late 2016 with a twisted inhomogeneous jet [@raiteri2017]. We note that bending and acceleration/deceleration of blazar jets are quite common indeed (e.g., @lister2013).

Optical subclasses {#sect:opticalsubclasses}
------------------

Phenomenologically, blazars can be divided into two classes based on their optical properties: FSRQs and BLOs. Previous studies showed that FSRQs tend to have higher RMs than BLOs [@ZT2004; @hovatta2012]. We collected all available core RM values from all frequency pairs and present the distributions of the (logarithmic) RMs of FSRQs and BLOs in Figure \[histoclass\]. The median RM values are $1.2\times10^4\rm\ rad/m^2$ and $4.8\times10^3\rm\ rad/m^2$ for FSRQs and BLOs, respectively – the value for FSRQs is higher than that for BLOs by a factor close to three. However, a Kolmogorov-Smirnov test [@press1992] finds a probability of 7% that the FSRQ and BLO values are drawn from the same parent population. Therefore, it is possible that their RM properties are intrinsically the same.

In Section \[sectscreen\], we claimed that the observed Faraday rotation mostly originates from jet sheaths. Relatively slow, possibly non-relativistic winds launched by an accretion disk that surround and confine the highly relativistic jet spine are one of the candidates for a jet sheath (e.g., @devilliers2005). A fundamental difference between FSRQs and BLOs is their accretion luminosities relative to their Eddington luminosities, above and below $\approx1\%$, respectively (e.g., @ghisellini2011, see also @PC2015 for further discussion). This suggests that sources in high accretion states tend to have larger RMs. A simple explanation would be that high accretion rates lead to relatively larger amounts of matter in jet sheaths. There is indeed evidence for a relation between the rate of matter injection into the jet and the accretion rate (e.g., @ghisellini2014 [@PT2017]), supporting this idea. However, the strength and degree of ordering of core magnetic fields as function of blazar subclass are poorly understood yet; the difference in RM may not be solely due to the difference in particle density.

Intrinsic polarization orientation
----------------------------------

Intrinsic EVPAs (projected onto the sky plane) of AGN jets can be obtained by correcting for Faraday rotation. It has been consistently shown that BLOs have intrinsic EVPAs well aligned with their jets, while a wide range of angles between EVPAs and jet orientations, sometimes seen as double-peaked distribution of relative angles, is observed for FSRQs (e.g., @LH2005 [@jorstad2007]). The good alignment and the mis-alignment were associated with a transverse or oblique shock and a conical shock, respectively [@jorstad2007]. These results, however, used RMs obtained from a single $\lambda^2$ law description of EVPA variation between 7 and 1 mm. The RM values were of the order $10^4\rm\ rad/m^2$.

However, our results show that there is a possibility that the core RMs of blazars can increase up to $\approx10^6\rm\ rad/m^2$ at $\approx250$ GHz (Figure \[paopt\]). The possible difference in core RM between FSRQs and BLOs, discussed in Section \[sect:opticalsubclasses\], suggests that it is easier for FSRQs to have high core RMs up to $\approx10^6\rm\ rad/m^2$ than for BLOs – unless the transition frequency for BLOs is much larger than for FSRQs, which seems not to be the case (Figure \[paopt\]). Even at 1 mm, the high RM of $\approx10^6\rm\ rad/m^2$ leads to an EVPA rotation of $\approx57\deg$. Therefore, one needs to take the frequency dependence of RM into account – especially for FSRQs – when comparing the intrinsic EVPAs at mm wavelengths with the direction of the inner jet. However, FSRQs are unlikely to have intrinsic EVPAs aligned with their jets (see also @yuan2001 [@gabuzda2006; @hovatta2016]) even after correcting for Faraday rotation; even their optical EVPAs, which do not suffer from strong Faraday rotation, show bimodal distributions in the angles between jets and EVPAs (@jorstad2007, see also @LS2000).

We observe RMs to increase with increasing frequency, meaning that intrinsic core EVPAs are different for different observing frequencies. Such a frequency dependence implies that polarized emission observed at higher frequencies comes from regions closer to the jet base. This indicates that intrinsic EVPAs can vary with distance from the jet base. A similar behaviour has been observed for a few sources in other studies. @OG2009a found that the jet of BL Lac shows EVPAs well aligned to the jet direction in inter-knot regions and even when the jet bends. [@gomez2016] showed that their high resolution polarization image of the same source shows smooth but non-negligible variations of EVPA upstream and downstream from the core. Both results were interpreted as the presence of helical magnetic fields in the jet. Similarly, different intrinsic EVPAs at different frequencies might imply the presence of helical magnetic fields in the core regions. However, a firm conclusion requires confirming the frequency dependence of RM at a wide range of observing frequencies with both short and long $\lambda^2$ spacings (see Section \[multipleshocks\]).

Multiple recollimation shocks in the cores {#multipleshocks}
------------------------------------------

Theoretical studies have shown that a series of recollimation shocks can form in relativistic jets: in analytic works (e.g., @DM1988), in hydrodynamic numerical simulations (e.g., @gomez1995 [@gomez1997; @agudo2001]), and in magneto-hydrodynamic simulations (e.g., @mizuno2015 [@marti2016]). Observationally, the presence of stationary features in AGN jets in addition to their VLBI cores has been verified in many studies (e.g., @jorstad2005). Especially, high resolution images of 3C 454.3 [@jorstad2010] and BL Lac [@gomez2016] revealed that their cores may consist of multiple stationary components. For BL Lac, emission upstream the radio core, leading to multi-wavelength flares when it passes through the core, was observed [@marscher2008; @gomez2016]; this supports the idea that the core can be identified with one of a series of recollimation shocks (e.g., @marscher2009).

We found that the EVPA–$\lambda^2$ relations of our sources are usually non-linear, instead showing breaks in their slopes. We obtained the RMs for pairs of adjacent frequencies and discovered that the core RMs systematically increase with observing frequency. Based on VLBA observations at 8 different frequencies from 4.6 to 43 GHz, [@OG2009a] showed that breaks in RM appear frequently, with the best-fit lines in the EVPA–$\lambda^2$ diagram connecting smoothly over a wide range of frequencies (though not for BL Lac in their sample). In contrast, [@kravchenko2017] presented large discontinuities between the different EVPA–$\lambda^2$ fits at much lower frequencies (between 2 and 5 GHz). Therefore, one might expect that potential discontinuities in EVPA rotations might not be substantial at mm wavelengths and the assumption underlying our analysis – no RM discontinuities – might be justified.

Furthermore, these studies showed that core EVPA rotations could be fitted well by a $\lambda^2$ law in some frequency ranges, then breaks, and then shows another good $\lambda^2$ fit in other frequency ranges (e.g., @OG2009a [@kravchenko2017]). Other studies obtained good $\lambda^2$ fits for the core EVPA rotations in most cases when they used relatively small frequency intervals, e.g., 8–15 GHz (e.g., @ZT2004 [@hovatta2012]). This indicates that polarized emission from a single emission region is dominant over relatively small frequency intervals, without showing a systematic increase of RMs as a function of frequency. However, over a wide range of frequencies, the RM–frequency relations appear to show multiple breaks; this implies that $\rm |RM| \propto \nu^a$ predicted by [@jorstad2007] assuming a continuous core-shift effect might not hold for narrow frequency intervals. One possible explanation is that blazar cores actually consist of multiple recollimation shocks and we observe polarized emission from one of the shocks in a given narrow frequency interval. As one goes to higher frequencies, polarized emission from inner shocks close to the jet base becomes more dominant due to lower opacity, leading to another good $\lambda^2$ fit with higher RMs.

Conclusions {#sect5}
===========

We studied polarization properties of 8 blazars – 5 FSRQs and 3 BLOs with multi-frequency simultaneous observations with the KVN at 22, 43, and 86 GHz. We investigated the nature of blazar radio cores by means of measuring Faraday rotation measures at different observing frequencies. Our work leads us to the following principal conclusions:

1.  We found that RMs increase with frequency, with median values of $2.62\times10^3\rm\ rad/m^2$ and $1.42\times10^4\rm\ rad/m^2$ for the frequency pairs 22/43 GHz and 43/86 GHz, respectively. These values are also higher than those obtained by [@hovatta2012] at 8.1–15.4 GHz for the same sources. The median values are described well by a power-law function with $\rm|RM|\propto\nu^{a}$ with $a=2.42$. When $a$ values are obtained separately for each source, they are distributed around $a=2$ with mean and standard deviation of $a=2.25\pm1.28$. This agrees with the expectation from core-shift [@jorstad2007] for many blazars at the KVN frequencies. This finding implies that the geometry of Faraday rotating media in blazar cores can be approximated as conical.

2.  We compared our KVN data with contemporaneous (within $\approx1$ week) optical polarization data from the Steward Observatory for a few sources. When we assume that the direction of EVPA rotation at radio frequencies is the same at optical wavelengths and that there is no $n\pi$ ambiguity, the optical data show a trend of EVPA rotation similar to that of the radio data. The RM values obtained with the optical data indicate that the power-law increase of RM with frequency continues up to a certain frequency, $\nu_{\rm trans}$, and then saturates, with $\rm |RM|\approx10^{5-6}\rm\ rad/m^2$ at $\approx250\rm\  GHz$, depending on source and flaring activity. We suggest that this saturation is due to the absence of core shift above $\nu_{\rm trans}$; instead, radio cores are standing recollimation shocks. This is in agreement with other studies which concluded that the radio cores of blazars cannot purely be explained as the unity optical depth surface of a continuous conical jet but are physical structures at least in some cases.

3.  We detected a sign change in the observed RMs of CTA 102 over $\approx1$ month, while the magnitudes of RM were roughly preserved. Since this source showed strong flaring at the time of our observations, we suggest that new relativistic jet components emerging from the core undergo acceleration/deceleration and/or jet bending, thus leading to a change in the direction of the line-of-sight component of helical magnetic fields in the jet because of relativistic abberation.

4.  We found indication that the absolute values of the core RMs of FSRQs are larger than those of BLOs at 22–86 GHz, which is consistent with results found at cm wavelengths. This difference might arise from FSRQs having higher accretion rates than BLOs, resulting in larger amounts of material in the central engine.

5.  For those sources which show non-linear EVPAs–$\lambda^2$ relations, the RM-corrected (intrinsic) EVPAs might be different at different frequencies and thus at different locations of the jets. A recent ultra-high resolution image of BL Lac observed with space VLBI shows that its intrinsic EVPAs in the core region vary with different locations indeed.

6.  We suggest that the systematic increase of RM as function of observing frequency appears only when covering sufficiently large ranges in frequency, with different $\lambda^2$ laws at different frequency ranges connecting smoothly. Combining this with the fact that linear EVPA–$\lambda^2$ relations are commonly observed over narrow frequency ranges suggests that blazars cores might consist of multiple recollimation shocks such that polarized emission from one of the shocks is dominant in a given narrow frequency range.

We thank the referee for constructive comments, which helped to improve the paper. The KVN is a facility operated by the Korea Astronomy and Space Science Institute. We are grateful to the staff of the KVN who helped to operate the array and to correlate the data. The KVN and a high-performance computing cluster are facilities operated by the 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 the KREONET (Korea Research Environment Open NETwork), which is managed and operated by the KISTI (Korea Institute of Science and Technology Information). This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. This study makes use of 43 GHz VLBA data from the VLBA-BU Blazar Monitoring Program (VLBA-BU-BLAZAR; <http://www.bu.edu/blazars/VLBAproject.html>), funded by NASA through the Fermi Guest Investigator Program. The VLBA is an instrument of the Long Baseline Observatory. The Long Baseline Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team [@lister2009]. We acknowledge financial support from the Korean National Research Foundation (NRF) via Global Ph.D. Fellowship Grant 2014H1A2A1018695 (J.P.) and Basic Research Grant NRF-2015R1D1A1A01056807 (S.T., J.C.A.). S. S. Lee was supported by NRF grant NRF-2016R1C1B2006697. GYZ is supported by Korea Research Fellowship Program through the NRF (NRF-2015H1D3A1066561). Correspondence should be addressed to S.T.

D-Term calibration and evolution \[appendixa\] {#appendixa}
==============================================

We show examples of calibration of instrumental polarization in Figure \[dtermcal\]. The calibrated complex visibility ratio, LR/LL with L and R referrring to left and right hand circular polarization, respectively, from the calibrator sources used for D-term calibration of our 2017 January data is shown on the complex plane. KVN antennas are on alt-azimuth mounts and instrumental polarization does not (usually) change over time as antennas change the direction of pointing. Target parallactic angles change over time, causing the polarization signal received by the alt-az antennas to vary. Since we perform parallactic angle correction in the early stage of data pre-processing, the polarization signal intrinsic to the target remains constant in the LR/LL plane while the instrumental polarization signal rotates with parallactic angle. Therefore, the rotation of the visibility ratio on the complex plane in Figure \[dtermcal\] is (mostly) due to instrumental polarization. See [@roberts1994] and [@aaron1997] for details of instrumental polarization calibration.

Amplitudes and phases of D-terms can be derived from the rotating pattern by assuming that the center of rotation does not vary with time. This is true only when (i) the antenna D-terms do not vary during the observation, which is true in most cases, and (ii) the calibration sources are either unpolarized or are polarized but spatially un-resolved. As we have 6 independent visibility ratios from three baselines and 6 unknown D-terms for three antennas and two polarizations, solutions of D-terms can be obtained from fitting the pattern observed in the complex plane.

The green and red lines in the left and the central panels of Figure \[dtermcal\] correspond to the expected LR/LL variation caused by the D-terms of two antennas, as obtained from the AIPS task LPCAL. The blue lines are the sum of contributions of the two D-terms, which is in good agreement with the data. The right panels show the visibility ratio after the D-term correction; the data are clustered around fixed points which mark polarization signals intrinsic to the target AGNs. The scatter is mostly due to thermal noise in the data; the noise becomes larger at higher observing frequencies.

Even precise D-term measurements may be affected by substantial systematic errors like non-stationary centers of rotation in the complex plane. Our calibrator sources are not perfectly un-polarized (or are polarized with sub-structure), and the D-terms may not always be constant during an observation run. We can estimate the errors on the D-terms by comparing the values obtained from different instrumental polarization calibrators as we did in Figure \[dtermcomparison\]. We use the standard devations of the D-terms obtained using different calibrators as errors. The errors are usually less than 1–2% but sometimes up to 3%. These errors will be transferred to the polarization quantities we used in our analysis and we considered these errors as described in Section \[sect2\].

The D-terms of the VLBA antennas are known to vary only on timescales of months or longer (e.g., @gomez2002). In Figure \[dtermevolution\] we check the stability of the KVN D-terms; their amplitudes seem to be mostly stable over $\approx4$ months but sometimes show non-negligible variability. Their standard deviations are usually less than 1–2% – in agreement with the formal errors of the D-terms.

Reliability check of KVN polarimetry \[appendixb\] {#appendixb}
==================================================

Thanks to the extensive monitoring of blazars, specifically by the MOJAVE program at 15 GHz and the VLBA-BU program[^8] at 43 GHz, we could check if our KVN maps are consistent with contemporaneous VLBA images. We picked 3C 273, which shows complex jet structure in both total intensity and polarization and thus is not used for our analysis of blazar cores, as reference source and compared (a) our KVN 22 GHz data observed on 2016 December 9 with the MOJAVE 15 GHz data observed one day after, (b) our KVN 43 GHz data observed on 2016 December 10 with the BU data observed on 2016 December 24, and (c) our KVN 86 GHz data observed on 2016 December 10 with the BU data in Figure \[comparison\]. The top panels show polarization maps generated from the KVN and VLBA data next to each other. The VLBA maps are convolved with the corresponding KVN beam.

The VLBA and KVN distributions of fractional jet polarization at 15 GHz and 22 GHz, respectively, are in good agreement with each other, showing higher degrees of polarization – up to $\approx70\%$ – at the northern edge of the jet located $\approx10$ mas from the core. Likewise, the 43-GHz VLBA and KVN data are consistent with each other except that the KVN maps show more polarized emission in the outer jet, i.e., $\approx10$ mas from the core. This might be because there is a time gap of $\approx2$ weeks between the observations and/or the KVN has only short baselines (with the maximum baseline length less than 500 km) and thus is more sensitive to the extended emission. However, the 86 GHz KVN map shows an additional polarization component near the core region, while the 43 GHz VLBA map does not have such a component but has an extended polarization near the core from the inner jet polarization component (at $\approx$ 1 mas from the core) by convolution of a large beam. Interestingly, the core of this source is usually unpolarized (e.g., @jorstad2005 [@attridge2005; @hada2016]), which has been attributed to strong Faraday depolarization or intrinsically very low polarization at the core. Its inner jet components at $\lesssim1$ mas from the core show $|\rm RM| \approx a\ few \times 10^4\ rad/m^2$ [@attridge2005; @jorstad2007; @hada2016]. Therefore, our result might indicate a detection of core polarization of 3C 273 at 86 GHz presumably because of less depolarization at higher frequency. We will verify this possibility in a forthcoming paper with more data at 86 and 129 GHz (Park et al. in preparation).

The bottom panels of Figure \[comparison\] show EVPA and degree of polarization as a function of $\lambda^2$ at a few locations in the jet marked in the top panels. The 15-GHz VLBA and 22-GHz KVN results are, to first order, consistent with each other but show non-negligible differences. This might be due to Faraday rotation with a RM of a few hundred $\rm rad/m^2$ in the jet, as reported by many other studies (e.g., @hovatta2012) and possibly different polarization structure in the jet at different frequencies. The VLBA and KVN data at 43 GHz are consistent with each other within errors, especially when considering the time gap of $\approx2$ weeks. Similarly, The VLBA 43 GHz and KVN 86 GHz data are in agreement with each other for the inner jet component at $\approx$ 1 mas from the core, taking into account the time gap and the small rotation measure. Therefore, we conclude that the polarimetry mode of KVN is reliable.

Aaron, S., 1997, EVN Memo \#78, “Calibration of VLBI polarization data”, <http://pc.astro.brandeis.edu/pdfs/evnmemo78.pdf> Agudo, I., Gómez, J.-L., Martí, J.-M., et al. 2001, , 549, 183 Agudo, I., Thum, C., Gómez, J. L., & Wiesemeyer, H. 2014, , 566, 59 Agudo, I., Thum, C., Molina, S. N., et al. 2018a, , 474, 1427 Agudo, I., Thum, C., Ramakrishnan, V., et al. 2018b, , 473, 1850 Algaba, J. C. 2013, , 429, 3551 Algaba, J.-C., Gabuzda, D. C., & Smith, P. S. 2011, , 411, 85 Algaba, J.-C., Gabuzda, D. C., & Smith, P. S. 2012, , 420, 542 Algaba, J.-C., Nakamura, M., Asada, K., & Lee, S.-S. 2017, , 834, 65 Algaba, J.-C., Zhao, G.-Y., Lee, S.-S., et al. 2015, JKAS, 48, 237 Asada, K., Inoue, M., Uchida, Y., et al. 2002, , 54, 39 Asada, K., Inoue, M., Nakamura, M., et al. 2008, , 682, 798 Asada, K., & Nakamura, M. 2012, , 745, 28 Asada, K., Nakamura, M., Doi, A., et al. 2014, , 781, 2 Attridge, J. M., Wardle, J. F. C., & Homan, D. C. 2005, , 633, 85 Blandford, R. D., & Königl, A. 1979, , 232, 34 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 Blinov, D., Pavlidou, V., Papadakis, I., et al. 2015, , 453, 1669 Broderick, A. E., & McKinney, J. C. 2010, , 725, 750 Burn, B. J. 1966, MNRAS, 133, 67 Casadio, C., Gómez, J. L., Grandi, P., et al. 2015, , 808, 162 Chatterjee, R., Marscher, A. P., Jorstad, S. G., et al. 2009, , 704, 1689 Chatterjee, R., Marscher, A. P., Jorstad, S. G., et al. 2011, , 734, 43 Cheung, C. C., Harris, D. E., & Stawarz, Ł. 2007, , 663, 65 Conway, R. G., Haves, P., Kronberg, P. P., et al. 1974, , 168, 137 Daly, R. A., & Marscher, A. P. 1988, , 334, 539 De Villiers, J.-P., Hawley, J. F., Krolik, J. H., & Hirose, S. 2005, , 620, 878 Dodson, R., Rioja, M. J., Molina, S. N., & Gómez, J. L. 2017, , 834, 177 Farnes, J. S., Gaensler, B. M., & Carretti, E. 2014, , 212, 15 Fromm, C. M., Ros, E., Perucho, M., et al. 2013, , 557, 105 Gabuzda, D. C., Knuettel, S., & Reardon, B. 2015, , 450, 2441 Gabuzda, D. C., Murray, É., & Cronin, P. 2004, , 351, 89 Gabuzda, D. C., Rastorgueva, E. A., Smith, P. S., & O’Sullivan, S. P. 2006, , 369, 1596 Ghisellini, G., Tavecchio, F., Foschini, L., & Ghirlanda, G. 2011, , 414, 2674 Ghisellini, G., Tavecchio, F., Maraschi, L., et al. 2014, , 515, 376 Gómez, J. L., Marscher, A. P., Alberdi, A., et al. 2002, VLBA Scientific Memo \#30 Gómez, J. L., Lobanov, A. P., Bruni, G., et al. 2016, , 817, 96 Gómez, J. L., Martí, J. M., Marscher, A. P., et al. 1995, , 449, 19 Gómez, J. L., Martí, J. M., Marscher, A. P., et al. 1997, , 482, 33 Hada, K., Doi, A., Kino, M., et al. 2011, , 477, 185 Hada, K., Kino, M., Doi, A., et al. 2013, , 775, 70 Hada, K., Kino, M., Doi, A., et al. 2016, , 817, 131 Hada, K., Park, J. H., M. Kino, et al. 2017, , 69, 71 Homan, D. C. 2012, , 747, 24 Homan, D. C., Kovalev, Y. Y., Lister, M. L., et al. 2006, , 642, 115 Hovatta, T., Aller, M. F., Aller, H. D., et al. 2014, , 147, 143 Hovatta, T., O’Sullivan, S., Martí-Vidal, I., et al. 2018, , in press (arXiv:1803.09982) Hovatta, T., Lindfors, E., Blinov, D., et al. 2016, , 596, 78 Hovatta, T., Lister, M. L., Aller, M. F., et al. 2012, , 144, 105 Jorstad, S., & Marscher, A. 2016, Galaxies, 4, 47 Jorstad, S. G., Marscher, A. P., Larionov, V. M., et al. 2010, , 715, 362 Jorstad, S. G., Marscher, A. P., Lister, M. L., et al. 2005, , 130, 1418 Jorstad, S. G., Marscher, A. P., Mattox, J. R., et al. 2001, , 556, 738 Jorstad, S. G., Marscher, A. P., Smith, P. S., et al. 2013, , 773, 147 Jorstad, S. G., Marscher, A. P., Stevens, J. A., et al. 2007, , 134, 799 Kang, S., Lee, S.-S., & Byun, D.-Y. 2015, JKAS, 48, 257 Kharb, P., Gabuzda, D. C., O’Dea, C. P., et al. 2009, , 694, 1485 Kim, D.-W., Trippe, S., Lee, S.-S., et al. 2017, JKAS, 50, 167 Kim, J.-Y., Trippe, S., Sohn, B. W., et al. 2015, JKAS, 48, 285 Kim, J.-Y., Lu, R.-S., Krichbaum, T. P., et al. 2016, Galaxies, 4, 39 Komissarov, S. S., Barkov, M. V., Vlahakis, N., & Königl, A. 2007, , 380, 51 Komissarov, S. S., Vlahakis, N., Königl, A., & Barkov, M. V. 2009, , 394, 1182 Kravchenko, E. V., Kovalev, Y. Y., & Sokolovsky, K. V. 2017, , 467, 83 Lee, S.-S., Byun, D.-Y., Oh, C. S., et al. 2011, , 123, 1398 Lee, S.-S., Byun, D.-Y., Oh, C. S., et al. 2015, JKAS, 48, 229 Lee, S.-S., Oh, C. S., Roh, D.-G., et al. 2015, JKAS, 48, 125 Lee, S.-S., Petrov, L., Byun, D.-Y., et al. 2014, , 147, 77 León-Tavares, J., Valtaoja, E., Tornikoski, M., et al. 2011, , 532, 146 Leppänen, K. J., Zensus, J. A., & Diamond, P. J. 1995, , 110, 2479 Lico, R., Gómez, J. L., Asada, K., & Fuentes, A. 2017, , 469, 1612 Lister, M. L., Aller, H. D., Aller, M. F., et al. 2009, , 137, 3718 Lister, M. L., Aller, M. F., Aller, H. D., et al. 2013, , 146, 120 Lister, M. L., & Homan, D. C. 2005, , 130, 1389 Lister, M. L., & Smith, P. S. 2000, , 541, 66 Lobanov, A. P. 1998, , 330, 79 Mahmud, M., Coughlan, C. P., Murphy, E., et al. 2013, , 431, 695 Mahmud, M., Gabuzda, D. C., & Bezrukovs, V. 2009, , 400, 2 Marscher, A. P. 1996, in ASP Conf. Ser. 110, Blazar Continuum Variability, ed. H. R. Miller, J. R. Webb, & J. C. Noble (San Francisco, CA: ASP), 248 Marscher, A. P. 2009, in ASP Conf. Ser. 402, Approaching Micro-Arcsecond Resolution with VSOP-2: Astrophysics and Technologies, ed. Y. Hagiwar et al. (San Francisco, CA: ASP), 194 Marscher, A. P., Jorstad, S. G., D’Arcangelo, F. D., et al. 2008, , 452, 966 Marscher, A. P., Jorstad, S. G., Gómez, J.-L., et al. 2002, , 417, 625 Martí, J. M., Perucho, M., & Gomez, J. L. 2016, , 831, 163 Martí-Vidal, I., Marcaide, J. M., Alberdi, A., et al. 2011, , 533, 111 Martí-Vidal, I., Muller, S., Vlemmings, W., et al. 2015, Science, 348, 311 Mertens, F., Lobanov, A. P., Walker, R. C., & Hardee, P. E. 2016, , 595, 54 Mizuno, Y., Gómez, J. L., Nishikawa, K.-I., et al. 2015, , 809, 38 Oh, J., Trippe, S., Kang, S., et al. 2015, JKAS, 48, 299 O’Sullivan, S. P., Brown, S., Robishaw, T., et al. 2012, , 421, 3300 O’Sullivan, S. P. & Gabuzda, D. C. 2009a, , 393, 429 O’Sullivan, S. P. & Gabuzda, D. C. 2009b, , 400, 26 O’Sullivan, S. P., Purcell, C. R., Anderson, C. S., et al. 2017, , 469, 4034 Pacholczyk, A. G. 1970, Radio Astrophysics (San Francisco: W. H. Freeman) Park, J., & Trippe, S. 2017, , 834, 157 Pasetto, A., Carrasco-González, C., O’Sullivan, S., et al. 2018, , in press (arXiv:1801.09731) Porth, O., Fendt, C., Meliani, Z., & Vaidya, B. 2011, , 737, 42 Potter, W. J., & Cotter, G. 2013a, , 429, 1189 Potter, W. J., & Cotter, G. 2013b, , 431, 1840 Potter, W. J., & Cotter, G. 2015, , 453, 4070 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in FORTRAN (2nd ed.; Cambridge: Cambridge Univ. Press) Pushkarev, A. B., Hovatta, T., Kovalev, Y. Y., et al. 2012, , 545, 113 Pushkarev, A. B., Kovalev, Y. Y., Lister, M. L., & Savolainen, T. 2017, , 468, 4992 Raiteri, C. M., Villata, M., Acosta-Pulido, J. A., et al. 2017, , 552, 374 Rani, B., Krichbaum, T. P., Marscher, A. P., et al. 2015, , 578, 123 Ramakrishnan, V., León-Tavares, J., Rastorgueva-Foi, E. A., et al. 2014, , 445, 1636 Rioja, M. J., Dodson, R., Jung, T., et al. 2014, , 148, 84 Rioja, M. J., Dodson, R., Malarecki, J., & Asaki, Y. 2011, , 142, 157 Roberts, D. H., Wardle, J. F. C., & Brown, L. F. 1994, , 427, 718 Savolainen, T., Wiik, K., Valtaoja, E., et al. 2002, , 394, 851 Shepherd, M. C. 1997, in Astronomical Data Analysis Software and Systems VI, eds. G. Hunt, & H. E. Payne (San Francisco: ASP), ASP Conf. Ser., 125, 77 Smith, P.S., Montiel, E., Rightley, S., et al. 2009, arXiv:0912.3621, 1094 2009 Fermi Symposium, eConf Proceedings C091122 Sokoloff, D. D., Bykov, A. A., Shukurov, A., et al. 1998, , 299, 189 Sokolovsky, K. V., Kovalev, Y. Y., Pushkarev, A. B., & Lobanov, A. P. 2011, , 532, 38 Taylor, A. R., Stil, J. M., & Sunstrum, C. 2009, , 702, 1230 Thum, C., Agudo, I., Molina, S. N., et al. 2018, , 473, 2506 Toma, K., & Takahara, F. 2013, Progr. Theor. Exper. Phys., 2013, 083E02 Tribble, P. C. 1991, , 250, 726 Urry, C. M., & Padovani, P. 1995, , 107, 803 Vlahakis, N., & Königl, A. 2004, , 605, 656 Walker, R. C., Hardee, P. E., Davies, F. B., et al. 2018, , 855, 128 Yuan, M. S., Tran, H., Wills, B., & Wills, D. 2001, in ASP Conf. Ser. 227, Blazar Demographics and Physics, ed. P. Padovani & C. M. Urry (San Francisco: ASP), 150 Zamaninasab, M., Savolainen, T., Clausen-Brown, E., et al. 2013, , 436, 3341 Zavala, R. T., & Taylor, G. B. 2002, , 566, 9 Zavala, R. T., & Taylor, G. B. 2003, , 589, 126 Zavala, R. T., & Taylor, G. B. 2004, , 612, 749 Zavala, R. T., & Taylor, G. B. 2005, , 626, 73 Zhao, G.-Y., Algaba, J. C., Lee, S. S., et al. 2018, , 155, 26

60

[^1]: <https://radio.kasi.re.kr/kvn/ksp.php#ksp003>

[^2]: This power-law index is obtained by including the RM at 15/22 GHz; see Figure \[result\] and Section \[oj287\] for details.

[^3]: The term ‘inner jet’ denotes any polarized component in the jet that can be resolved from the core by instruments with higher angular resolution than the KVN but cannot be (well) resolved by the KVN itself, like the extended linear polarization structure of 3C 279 at $\approx1$ mas from the core seen in the BU map of 14 January 2017 (see <https://www.bu.edu/blazars/VLBA_GLAST/3c279/3C279jan17_map.jpg>)

[^4]: <http://www.physics.purdue.edu/astro/MOJAVE/index.html>

[^5]: <https://dept.astro.lsa.umich.edu/obs/radiotel/gif/0851_202.gif>

[^6]: <http://sma1.sma.hawaii.edu/callist/callist.html?plot=1642%2B398> This might be related to a substantial change in $a$ within $\approx1$ month for this source, though the sign of RM does not change during the period.

[^7]: <http://james.as.arizona.edu/~psmith/Fermi>

[^8]: https://www.bu.edu/blazars/VLBAproject.html
