# Kinematics of Parsec-Scale Jets of Gamma-Ray Blazars at 43 GHz during Ten
Years of the VLBA-BU-BLAZAR Program

Zachary R. Weaver Institute for Astrophysical Research, Boston University, 725
Commonwealth Avenue, Boston, MA 02215, USA  Svetlana G. Jorstad Institute for
Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA
02215, USA  Astronomical Institute, St. Petersburg State University,
Universitetskij, Pr. 28, Petrodvorets, St. Petersburg 198504, Russia  Alan P.
Marscher Institute for Astrophysical Research, Boston University, 725
Commonwealth Avenue, Boston, MA 02215, USA  Daria A. Morozova Astronomical
Institute, St. Petersburg State University, Universitetskij, Pr. 28,
Petrodvorets, St. Petersburg 198504, Russia  Ivan S. Troitsky Astronomical
Institute, St. Petersburg State University, Universitetskij, Pr. 28,
Petrodvorets, St. Petersburg 198504, Russia  Iván Agudo  Instituto de
Astrofísica de Andalucía (IAA), CSIC, Glorieta de la Astronomía s/n, E-18008,
Granada, Spain  José L. Gómez  Instituto de Astrofísica de Andalucía (IAA),
CSIC, Apartado 3004, E-18080, Granada, Spain  Anne Lähteenmäki  Aalto
University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä,
Finland  Aalto University Department of Electronics and Nanoengineering, P.O.
Box 15500, FI-00076, AALTO, Finland  Joni Tammi  Aalto University Metsähovi
Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland  Merja
Tornikoski  Aalto University Metsähovi Radio Observatory, Metsähovintie 114,
FI-02540 Kylmälä, Finland

###### Abstract

We analyze the parsec-scale jet kinematics from 2007 June to 2018 December of
a sample of $\gamma$-ray bright blazars monitored roughly monthly with the
Very Long Baseline Array at 43 GHz under the VLBA-BU-BLAZAR program. We
implement a novel piece-wise linear fitting method to derive the kinematics of
521 distinct emission knots from a total of 3705 total intensity images in 22
quasars, 13 BL Lacertae objects, and 3 radio galaxies. Apparent speeds of
these components range from $0.01c$ to $78c$, and 18.6% of knots (other than
the “core”) are quasi-stationary. One-fifth of moving knots exhibit non-
ballistic motion, with acceleration along the jet within 5 pc of the core
(projected) and deceleration farther out. These accelerations occur mainly at
locations coincident with quasi-stationary features. We calculate the physical
parameters of 273 knots with statistically significant motion, including their
Doppler factors, Lorentz factors, and viewing angles. We determine the typical
values of these parameters for each jet and the average for each subclass of
active galactic nuclei. We investigate the variability of the position angle
of each jet over the ten years of monitoring. The fluctuations in position of
the quasi-stationary components in radio galaxies tend to be parallel to the
jet, while no directional preference is seen in the components of quasars and
BL Lacertae objects. We find a connection between $\gamma$-ray states of
blazars and their parsec-scale jet properties, with blazars with brighter 43
GHz cores typically reaching higher $\gamma$-ray maxima during flares.

galaxies: active — galaxies: jets — techniques: interferometric

††facilities: VLBA, Fermi-LAT, Metsähovi Radio Observatory††software: AIPS
(van Moorsel et al., 1996), astropy (Astropy Collaboration et al., 2013,
2018), Difmap (Shepherd, 1997), lifelines (Davidson-Pilon, 2019), pwlf (Jekel
& Venter, 2019), statsmodels (Seabold & Perktold, 2010)

Draft 3in

##  1 Introduction

The ultra-relativistic flow speeds of the jets of some active galactic nuclei
(AGNs) strongly affect their observed characteristics through beaming,
aberration, and light-travel delays. These physical effects cause the illusion
of apparent superluminal motion, greatly enhanced flux densities, and reduced
timescales of variability of flux and polarization (e.g., Blandford et al.,
1977). Blazars (Angel & Stockman, 1980), with their jet directions oriented
nearly along our line of sight (e.g., Lister et al., 2016), exhibit the most
extreme characteristics of all AGNs. The extraordinary resolution afforded by
very-long baseline interferometry (VLBI) enables detailed studies of the inner
jets of AGNs on distance scales from $<1$ to $>100$ pc from the central
supermassive black hole. Monitoring programs over at least a decade are able
to (1) systematically characterize the long-term variability of the parsec-
scale jet, including knot motions, jet wobbling, and precession; (2) reveal
the connection between multi-wavelength outbursts and multi-messenger events
in the jet; and (3) provide an observational basis for theoretical models of
jet formation, acceleration, and collimation.

Most blazar jets feature “knots” of emission that travel away from a (presumed
stationary) “core” at apparent speeds that range from $\lesssim 1c$ to
$\gtrsim 50c$ (e.g., Lister et al., 2016; Jorstad et al., 2017). In a number
of blazars, the jets are observed (in projection) to undergo rapid swings of
position angle in their innermost regions, but remain collimated on kiloparsec
scales (for a review of several sources, see Agudo, 2009). The projected jet
position angles can vary considerably (by up to 150$\arcdeg$ over a 12-16 yr
time span; Lister et al., 2013). Monotonic swings with no apparent periodicity
(“wobbles”) have been reported in some objects (e.g., Jorstad et al., 2004;
Agudo et al., 2007, 2012), while precession (sinusoidal variations) has been
claimed in others (e.g., Stirling et al., 2003; Bach et al., 2005; Lobanov &
Roland, 2005; Savolainen et al., 2006; Marti-Vidal et al., 2013; Britzen et
al., 2018, 2019). The origin of the wobbles and precession is not well
understood, but the presence of such changes in only the innermost (parsec-
scale) regions of the jet implies that the causal mechanisms are related to
the origin of the jet, e.g., accretion disk precession (e.g., McHardy et al.,
1990; Conway & Wrobel, 1995; Villata & Raiteri, 1999; Britzen et al., 2005;
Perucho et al., 2012; Fromm et al., 2013; Larionov et al., 2013; Britzen et
al., 2017; Raiteri et al., 2017, 2021), instabilities in the jet flow (e.g.,
Nakamura et al., 2001; Moll et al., 2008; Mignone et al., 2010; Liska et al.,
2018), or interactions in a binary supermassive black hole system (e.g.,
Britzen et al., 2018; Qian et al., 2019).

In order to understand the physical processes in blazars, it is important to
relate changes observed in the radio jet to multi-wavelength variability and
outbursts. Connections between the high-energy emission and radio jets of AGNs
have been theorized for $\sim$40 years (e.g., Königl, 1981; Marscher, 1987),
with observational evidence provided by joint radio and $\gamma$-ray studies,
the latter based on data from the _Compton Gamma Ray Observatory_ (Hartman et
al., 1999). The results have highlighted the dichotomy between blazar
subclasses. In the traditional AGN unification scheme (Urry & Padovani, 1995),
flat-spectrum radio quasars (FSRQs) and BL Lacertae objects (BLs) are
differentiated by their optical properties and radio morphology (e.g., Wardle
et al., 1984). During the _Compton_ era, FSRQs were found to have higher
$\gamma$-ray fluxes when they also exhibited high optical polarization and a
high-radio-frequency outburst (Valtaoja & Terasranta, 1995). There was no
similarly clear correlation between the $\gamma$-ray emission in BLs and their
radio variations (Lähteenmäki & Valtaoja, 2003). For both subtypes, however,
there existed a correlation between the brightness of the radio “core” and
$\gamma$-ray flux (Jorstad et al., 2001b), as well as a connection between the
ejections of superluminal radio knots and $\gamma$-ray outbursts (Jorstad et
al., 2001a).

The _Fermi_ Large Area Telescope (LAT) has ushered in a new era of
$\gamma$-ray-radio studies. Individual sources show a wide variety of behavior
during major $\gamma$-ray events. The typical duration of a high $\gamma$-ray
state in blazars is several months (Williamson et al., 2014), although
shorter-timescale fluctuations are often observed (e.g., Agudo et al., 2011a;
Weaver et al., 2019, 2020). The majority of these $\gamma$-ray and very-high-
energy (VHE, $\geq 100$ GeV) outbursts can be connected with the propagation
of a knot through the core or other quasi-stationary component in the parsec-
scale radio jet (e.g., Marscher et al., 2010; Agudo et al., 2011a; Casadio et
al., 2015; Liodakis et al., 2020; H. E. S. S. Collaboration et al., 2021;
MAGIC Collaboration et al., 2021). Two major issues exist with this model for
blazar flares: the presence of “orphan” $\gamma$-ray flares (e.g., Krawczynski
et al., 2004; MacDonald et al., 2015) with no counterparts at longer
wavelengths, and the need for an intense external photon field at these
locations in the jet (e.g., Aleksić et al., 2011; Joshi et al., 2014), the
source of which is unknown. Emission-line flares with short time delays
relative to the $\gamma$-ray outburst have also been observed in several
blazars (León-Tavares et al., 2013; Larionov et al., 2020; Hallum et al.,
2022), which could provide additional seed photons (Isler et al., 2015; León-
Tavares et al., 2015; Larionov et al., 2020). Detection of a 290 TeV neutrino
with a direction of origin coincident with the flaring blazar TXS 0506+056
(IceCube Collaboration et al., 2018a, b), suggests an association between
flaring blazars and high-energy neutrinos (e.g., Rodrigues et al., 2018;
Murase et al., 2018). Thus, blazar jets may produce high-energy hadronic
cosmic-rays (e.g., Keivani et al., 2018; Cerruti et al., 2019) alongside
electromagnetic outbursts and changes in the parsec-scale structure (for
discussion of the parsec-scale jet of TXS 0506+056, see Li et al., 2020).

VLBI observations can also provide constraints on theoretical models of AGN
jet launching, acceleration, and collimation, which are challenging to model.
The launching mechanism in AGNs can be ascribed to one (or both) of two
processes. The Blandford-Znajek mechanism (Blandford & Znajek, 1977) utilizes
the rotational energy of the black hole, with the spin of the black hole
winding up the magnetic field lines anchored to the event horizon, resulting
in a relativistic Poynting flux dominated jet (e.g., Beskin & Kuznetsova,
2000; Komissarov, 2001, 2005; Tchekhovskoy et al., 2010). Alternatively,
magnetic field lines anchored to the accretion disk can launch mildly
relativistic mass-dominated winds (Blandford & Payne, 1982). One of the
primary goals of VLBI monitoring programs of blazar jets is to provide an
observational framework for comparing numerical simulations of jets (e.g.,
Vlahakis & Königl, 2003; Komissarov, 2005; Lyubarsky, 2009; Meier, 2012; Pu et
al., 2015; Chatterjee et al., 2019) with the actual parsec-scale behavior.

A number of observing programs follow time changes in the jets of a large
number of blazars. The largest has monitored several hundreds of the brightest
AGNs in the northern sky, first under the 2 cm VLBA survey (Kellermann et al.,
1998), and later as part of the
MOJAVE111https://www.physics.purdue.edu/MOJAVE/ program. Lister et al. (2019)
find that $60\%$ of the 409 monitored sources exhibit acceleration of emission
features in the jet. Other than a few outliers, the 15 GHz emission features
have a distribution of apparent velocities that fall within $\pm 40\%$ of the
median speed of a given source. No jet features are reported to have a bulk
Lorentz factor $\Gamma>50$. Through Monte-Carlo modeling of the parent
population, the authors show that the data are consistent with a jet
population that has a simple unbeamed power-law luminosity function
incorporating pure luminosity evolution, and a power-law Lorentz factor
distribution ranging from 1.25 to 50, with a slope $-1.4\pm 0.2$.

The Boston University group and its collaborators have been monitoring blazar
jets with the VLBA at 43 GHz since 1993. The research focuses on the
connection between jet properties and multi-wavelength behavior. This includes
analysis of events in the jets of $\gamma$-ray blazars during the _Compton_
era (Jorstad et al., 2001b), accretion disk-jet associations (Marscher et al.,
2002; Chatterjee et al., 2011), and relations between optical and millimeter-
wave linear polarization (Jorstad et al., 2005). Over the past 10 years, we
have been monitoring the jet kinematics of a sample of 38 blazars at a roughly
monthly cadence at 43 GHz under the VLBA-BU-BLAZAR
program222http://www.bu.edu/blazars/BEAM-ME.html. For comparison with the
longer-wavelength data, we obtain well-sampled $\gamma$-ray light curves from
the _Fermi_ LAT. The higher observing frequency relative to other VLBA
monitoring programs allows emission features to be tracked closer to the
“core” of the sources. Also, the more frequent monitoring allows faster,
shorter-lived features to be observed. Based on the first 5 years of data, we
found that nearly one-third of moving knots show evidence of acceleration in
the jet (Jorstad et al., 2017). Furthermore, knots in FSRQs tend to have
higher Doppler and Lorentz factors, with smaller viewing and opening angles,
than the knots in BLs and radio galaxies (RGs). A number of radio jet -
$\gamma$-ray event associations have been detected using these data (Marscher
et al., 2010; Agudo et al., 2011a, b; Wehrle et al., 2012; Jorstad et al.,
2013; Casadio et al., 2015; Jorstad & Marscher, 2016; Larionov et al., 2020;
H. E. S. S. Collaboration et al., 2021; MAGIC Collaboration et al., 2021).

In this paper, we provide the data and detailed analyses of the parsec-scale
jet behavior of a sample of $\gamma$-ray bright AGNs over our decade-long
monitoring program, VLBA-BU-BLAZAR. We implement a new piece-wise linear
fitting method to determine the jet kinematics. For consistency, we reanalyze
previously-published data to produce a uniform data set that covers the first
10.5 years of VLBA monitoring associated with _Fermi_ LAT observations, 2007
June to 2018 December. The paper is structured as follows. In $\S$2, we detail
the observations, their calibration, and analysis from 2013 January through
2018 December, and describe the piece-wise fitting of the core-knot
separations versus time that we apply to all of the data. The data for
individual knots and sources are available in Tables 2 and 3. The parsec-scale
jet structure and the behavior of the jet position angle with time (Table 7)
of each source are discussed in $\S$3. Section 4 describes the
moving/stationary features (Tables 8 and 9) and accelerations (Table 10) of
all knots in our sample. We determine the statistically robust physical
parameters of knots in our sample in $\S$5 (Table 12), and provide estimates
for the “typical” physical parameters of each jet (Table 13). We briefly
discuss connections between millimeter-wave and $\gamma$-ray states of the
AGNs in $\S$6. A summary of our findings is given in $\S$7. Throughout this
work, we assume a standard $\Lambda$CDM cosmology with rounded values of
cosmological parameters based on those obtained by the Planck collaboration
(Planck Collaboration et al., 2020), the results from SNIa studies (Riess et
al., 2016), and VLT-KMOS data (González-Morán et al., 2021): $H_{\circ}=70$ km
s-1 Mpc-1, $\Omega_{\text{m}}=0.3$, and $\Omega_{\Lambda}=0.7$.

##  2 Observations and Data Analysis

Table 1: The VLBA-BU-BLAZAR Sample. Source | Name | Type | $z$ | $\langle S_{37}\rangle$ | $\langle S_{\gamma}\rangle$ | $\langle S_{R}\rangle$
---|---|---|---|---|---|---
|  |  |  | [Jy] | [$10^{-8}$ ph cm-2 s-1] | [mag]
0219+428aaSource added to the sample in 2009. | 3C 66A | BL | 0.444ffRedshift has not been confirmed. | $\phantom{0}0.71\pm 0.24$ | $\phantom{0}10.81\pm\phantom{0}12.86$ | $14.40\pm 0.42$ |
0235+164 | A0 0235+16 | BL | 0.940ffRedshift has not been confirmed. | $\phantom{0}1.99\pm 1.16$ | $\phantom{0}11.77\pm\phantom{0}14.91$ | $17.45\pm 0.98$ |
0316+413bbSource added to the sample in 2010. | 3C 84 | RG | 0.0176 | $26.11\pm 4.57$ | $\phantom{0}35.54\pm\phantom{0}20.40$ | $12.82\pm 0.14$ |
0336–019 | CTA 26 | FSRQ | 0.852 | $\phantom{0}2.24\pm 0.56$ | $\phantom{0}12.34\pm\phantom{0}10.13$ | $16.86\pm 0.52$ |
0415+379 | 3C 111 | RG | 0.0485 | $\phantom{0}4.06\pm 2.13$ | $\phantom{0}\phantom{0}7.43\pm\phantom{0}\phantom{0}4.68$ | $17.38\pm 0.22$ |
0420–014 | OA 129 | FSRQ | 0.916 | $\phantom{0}4.11\pm 2.14$ | $\phantom{0}\phantom{0}7.68\pm\phantom{0}\phantom{0}5.18$ | $17.50\pm 0.78$ |
0430+052ccSource added to the sample in 2012. | 3C 120 | RG | 0.033 | $\phantom{0}3.40\pm 1.64$ | $\phantom{0}\phantom{0}5.43\pm\phantom{0}\phantom{0}4.84$ | $13.90\pm 0.13$ |
0528+134 | PKS 0528+13 | FSRQ | 2.060 | $\phantom{0}1.84\pm 1.20$ | $\phantom{0}\phantom{0}8.85\pm\phantom{0}\phantom{0}9.93$ | $19.38\pm 0.27$ |
0716+714 | S5 0716+71 | BL | 0.3ffRedshift has not been confirmed. | $\phantom{0}2.44\pm 1.19$ | $\phantom{0}20.78\pm\phantom{0}14.67$ | $13.19\pm 0.63$ |
0735+178 | PKS 0735+17 | BL | 0.424 | $\phantom{0}0.77\pm 0.23$ | $\phantom{0}\phantom{0}6.21\pm\phantom{0}11.77$ | $15.98\pm 0.47$ |
0827+243 | OJ 248 | FSRQ | 0.939 | $\phantom{0}1.23\pm 0.61$ | $\phantom{0}\phantom{0}6.87\pm\phantom{0}11.26$ | $16.91\pm 0.39$ |
0829+046 | OJ 049 | BL | 0.182 | $\phantom{0}0.83\pm 0.37$ | $\phantom{0}\phantom{0}6.53\pm\phantom{0}\phantom{0}3.93$ | $15.52\pm 0.38$ |
0836+710 | 4C+71.07 | FSRQ | 2.172 | $\phantom{0}2.12\pm 0.69$ | $\phantom{0}13.91\pm\phantom{0}19.40$ | $16.60\pm 0.09$ |
0851+202 | OJ 287 | BL | 0.306 | $\phantom{0}5.94\pm 2.20$ | $\phantom{0}\phantom{0}9.08\pm\phantom{0}\phantom{0}7.46$ | $14.37\pm 0.49$ |
0954+658 | S4 0954+65 | BL | 0.368ffRedshift has not been confirmed. | $\phantom{0}1.19\pm 0.40$ | $\phantom{0}10.87\pm\phantom{0}20.24$ | $15.58\pm 0.96$ |
1055+018aaSource added to the sample in 2009. | 4C+01.28 | FSRQggClassified as a BL object in Jorstad et al. (2017). | 0.890 | $\phantom{0}5.22\pm 1.36$ | $\phantom{0}\phantom{0}9.91\pm\phantom{0}\phantom{0}5.72$ | $15.92\pm 0.83$ |
1101+384bbSource added to the sample in 2010. | Mkn 421 | BL | 0.030 | $\phantom{0}0.49\pm 0.16$ | $\phantom{0}17.83\pm\phantom{0}\phantom{0}9.50$ | $12.41\pm 0.39$ |
1127–145 | PKS 1127–14 | FSRQ | 1.184 | $\ldots$ | $\phantom{0}\phantom{0}6.80\pm\phantom{0}\phantom{0}5.05$ | $16.74\pm 0.18$ |
1156+295 | 4C+29.45 | FSRQ | 0.729 | $\phantom{0}2.08\pm 0.83$ | $\phantom{0}16.51\pm\phantom{0}21.60$ | $15.71\pm 1.00$ |
1219+285 | WCom | BL | 0.102 | $\phantom{0}0.41\pm 0.16$ | $\phantom{0}\phantom{0}6.29\pm\phantom{0}\phantom{0}2.65$ | $15.26\pm 0.37$ |
1222+216 | 4C+21.35 | FSRQ | 0.435 | $\phantom{0}1.98\pm 0.61$ | $\phantom{0}33.62\pm\phantom{0}55.61$ | $15.10\pm 0.43$ |
1226+023 | 3C 273 | FSRQ | 0.158 | $17.76\pm 4.58$ | $\phantom{0}29.03\pm\phantom{0}40.00$ | $12.63\pm 0.13$ |
1253–055 | 3C 279 | FSRQ | 0.538 | $19.00\pm 5.31$ | $\phantom{0}66.30\pm 117.31$ | $15.02\pm 0.66$ |
1308+326 | B2 1308+32 | FSRQ | 0.998 | $\phantom{0}1.74\pm 0.72$ | $\phantom{0}\phantom{0}5.25\pm\phantom{0}\phantom{0}4.24$ | $18.03\pm 0.83$ |
1406–076 | PKS 1406–07 | FRSQ | 1.494 | $\phantom{0}0.94\pm 0.26$ | $\phantom{0}\phantom{0}6.54\pm\phantom{0}10.19$ | $18.28\pm 0.42$ |
1510–089 | PKS 1510–08 | FSRQ | 0.361 | $\phantom{0}3.31\pm 1.25$ | $\phantom{0}96.60\pm 120.47$ | $15.88\pm 0.53$ |
1611+343 | DA 406 | FSRQ | 1.40 | $\phantom{0}3.11\pm 0.71$ | $\phantom{0}\phantom{0}5.25\pm\phantom{0}\phantom{0}2.44$ | $17.31\pm 0.21$ |
1622–297hhRemoved from the sample after 2017 July. | PKS 1622–29 | FSRQ | 0.815 | $\ldots$ | $\phantom{0}12.09\pm\phantom{0}\phantom{0}7.57$ | $18.18\pm 0.25$ |
1633+382 | 4C+38.41 | FSRQ | 1.814 | $\phantom{0}3.78\pm 1.39$ | $\phantom{0}30.31\pm\phantom{0}30.63$ | $16.68\pm 0.63$ |
1641+399 | 3C 345 | FSRQ | 0.595 | $\phantom{0}5.83\pm 1.16$ | $\phantom{0}\phantom{0}9.22\pm\phantom{0}\phantom{0}5.83$ | $17.15\pm 0.46$ |
1652+398ddSource added to the sample in 2014. | Mkn 501 | BL | 0.034 | $\phantom{0}0.93\pm 0.17$ | $\phantom{0}\phantom{0}5.27\pm\phantom{0}\phantom{0}2.38$ | $13.26\pm 0.09$ |
1730–130 | NRAO 530 | FSRQ | 0.902 | $\phantom{0}3.47\pm 0.94$ | $\phantom{0}11.96\pm\phantom{0}\phantom{0}9.99$ | $17.63\pm 0.38$ |
1749+096aaSource added to the sample in 2009. | OT 081 | BL | 0.322 | $\phantom{0}3.82\pm 1.21$ | $\phantom{0}\phantom{0}8.89\pm\phantom{0}12.39$ | $16.49\pm 0.52$ |
1959+650eeSource added to the sample in 2017. | 1ES 1959+65 | BL | 0.047 | $\phantom{0}0.32\pm 0.14$ | $\phantom{0}\phantom{0}5.13\pm\phantom{0}\phantom{0}3.29$ | $14.24\pm 0.25$ |
2200+420 | BL Lac | BL | 0.069 | $\phantom{0}4.46\pm 2.47$ | $\phantom{0}40.48\pm\phantom{0}25.08$ | $13.64\pm 0.56$ |
2223–052 | 3C 446 | FSRQ | 1.404 | $\phantom{0}3.80\pm 2.00$ | $\phantom{0}\phantom{0}5.02\pm\phantom{0}\phantom{0}2.89$ | $18.09\pm 0.35$ |
2230+114 | CTA 102 | FSRQ | 1.037 | $\phantom{0}3.50\pm 0.89$ | $118.81\pm 185.75$ | $15.65\pm 1.04$ |
2251+158 | 3C 454.3 | FSRQ | 0.859 | $17.01\pm 7.72$ | $206.62\pm 337.67$ | $15.02\pm 0.79$ |

As part of the VLBA-BU-BLAZAR monitoring program, we have obtained roughly
monthly observations with the VLBA at 43 GHz (7 mm wavelength) of a sample of
AGNs detected as $\gamma$-ray sources. The results of the observations from
2007 June to 2013 January have been presented in Jorstad et al. (2017). In
this work, we extend the analyzed period with observations obtained from 2013
January through 2018 December, bringing the total monitoring period to 10.5
years. In order to interpret the results of our study more consistently, we
apply a new formalism (described below) to the entire data set to calculate
the knot motions and kinematic parameters.

The sample studied here consists of a total of 38 sources, of which 22 are
FSRQs, 13 are BLs, and 3 are RGs. Each source has been detected at
$\gamma$-ray energies by the _Fermi_ LAT (Abdo et al., 2010; Acero et al.,
2015), with an average flux density at 43 GHz exceeding 0.5 Jy, a declination
north of $-30\arcdeg$, and an optical magnitude in _R_ -band brighter than
$18.5^{\text{m}}$. The original sample is described in Jorstad et al. (2017).
Several additional sources that have been revealed to have occasionally strong
$\gamma$-ray flux (as observed with the _Fermi_ -LAT) and sufficiently high
optical flux to meet the selection criteria have been included in the current
sample.

The observations of the sample were performed roughly monthly via dynamical
scheduling, over a 24-hr period at each epoch, using all available VLBA
antennas. (At some epochs 1-2 antennas could not collect data owing to
technical problems or weather conditions.) A total of 32 or 33 objects were
observed during each epoch, with 40-45 min integration time per source (8-9
scans of $\sim 5$-min duration), spread over the 6-10 hr time span with the
optimal visibility of the object. Eight sources were observed during
alternating sessions; these objects are in crowded right ascension ranges and
exhibit relatively long timescales of variability of the jet structure at 43
GHz. Although the observations and data processing were carried out in a
similar manner as in Jorstad et al. (2017), with the latter using the
Astronomical Image Processing System (AIPS; van Moorsel et al., 1996), some
details were different. A new data recording system (Mark5C) was installed,
allowing a data recording rate of 2048 Mb per second, as opposed to the
previous rate of 512 Mb per second. The resulting increase in data quality
allows a sufficient signal-to-noise ratio to be achieved without averaging
intermediate frequency bands (IFs). Thus, since 2014 February the final
calibrated data available on the project website include visibilities at 4 IFs
(centered at 43.0075, 43.0875, 43.1515, and 43.2155 GHz). We have calibrated
the instrumental polarization separately for each IF.

We have continued to monitor the optical _BVRI_ photometric brightness and _R_
-band linear polarization of the sample at the 1.83-m Perkins Telescope (PTO,
Boston University333Operated by Lowell Observatory prior to 2019, Flagstaff,
Arizona) over $\sim 1$ week time intervals during dark/gray lunar phases. One
such session occurred during most months.

Table 1 gives a list of the monitored sources, along with their
characteristics: 1—name based on B1950 coordinates; 2—common name; 3—subclass
of blazar; 4—redshift, $z$, according to the NASA Extragalactic
Database444https://ned.ipac.caltech.edu/; 5—average total intensity at 37 GHz,
$\langle S_{37}\rangle$, and its standard deviation as observed with the
Metsähovi Radio Observatory of Aalto University, Finland, between 2007 January
and 2018 December; 6—average $\gamma$-ray photon flux at 0.1-200 GeV, $\langle
S_{\gamma}\rangle$, and its standard deviation based on the $\gamma$-ray light
curve of each source, calculated from the photon and spacecraft data provided
by the _Fermi_ Space Science Center (for details, see Williamson et al., 2014)
between 2007 January and 2018 December; 7—average _R_ -band magnitude,
$\langle S_{R}\rangle$ and its standard deviation according to our
observations at the Perkins Telescope between 2007 January and 2018 December.

###  2.1 Total Intensity Image Modeling

The parsec-scale jet morphology of blazars features a “core” (e.g., Jorstad et
al., 2001b; Kellermann et al., 2004; Jorstad et al., 2005; Lister et al.,
2013; Jorstad et al., 2017; Lister et al., 2018), an assumed stationary
feature located at the upstream end of the jet seen in the VLBI images. For
most sources, the core is usually the brightest feature in the jet (Jorstad et
al., 2017), but sometimes in some blazars an emission component downstream of
the core can be more prominent (e.g., 2251+158; Jorstad et al., 2001b). The
exact location of this core relative to the central black hole is frequency
dependent as the core represents the region where the jet becomes optically
thick at a given frequency, with VLBI at higher frequencies being able to
image closer to the black hole (see, e.g., O’Sullivan & Gabuzda, 2009; Event
Horizon Telescope Collaboration et al., 2019). One or more “knots” of emission
are often observed downstream of the core (e.g., Jorstad et al., 2001b). Some
of these are quasi-stationary, without systematic motions relative to the
core, while others move down the jet, in most cases at apparent superluminal
speeds (e.g., Cohen et al., 1971, 1977; Jorstad et al., 2001b, 2017; Lister et
al., 2018).

In order to represent the total intensity structure of each source at each
epoch, we model the frequency-averaged visibility data with a series of
circular Gaussian components that best fit the data. The model fitting
involves iterations of the MODELFIT task within the Difmap (Shepherd, 1997)
software package. We use the term “knot” to refer to these Gaussian
components, which correspond to a (usually) compact feature of enhanced
brightness in the jet. Initially, a single circular Gaussian component
approximating the brightness distribution of the core is used, then knots are
added at the approximate locations of bright features identified in the image.
Each addition of a knot is followed by hundreds of iterations with MODELFIT to
determine the parameters of the knot that yield the best agreement between the
model and $uv$ data, according to a $\chi^{2}$ test. The knot parameters that
we fit are: $S$—flux density; $R_{\text{obs}}$—distance with respect to the
core; $\Theta_{\text{obs}}$—relative position angle (P.A.) with respect to the
core, measured north through east; and $a$—the angular size of the knot,
corresponding to the FWHM of the circular Gaussian component. The iterative
process of adding knots is ended when the addition of a new knot does not
significantly improve the $\chi^{2}$ value of the model. Often, the model of a
previous epoch is used as the starting model for the next epoch for a given
source, as we assume that the jet does not drastically change structure
between our roughly monthly observations. We utilize the previously-published
Difmap models of each source at each epoch from 2007 June to 2012 December
(Jorstad et al., 2017). In this work we provide models for each source from
2013 January to 2018 December.

Table 2: Modeling of Jets by Gaussian Components 2013–2018

Source | $\chi^{2}$ | RMS | Epoch | MJD | $S\pm\sigma_{S}$ | $R\pm\sigma_{R}$ | $\Theta\pm\sigma_{\Theta}$ | $a\pm\sigma_{a}$ | $T_{b,\text{obs}}$
---|---|---|---|---|---|---|---|---|---
|  |  | [yr] |  | [Jy] | [mas] | [deg] | [mas] | [$\times 10^{8}$ K]
0219+428 | 1.131 | 0.050 | 2013.038 | 56307 | $0.129\pm 0.009$ | $0.000\phantom{\pm}\phantom{0}\phantom{.}\phantom{0}\phantom{0}\phantom{0}$ | $\phantom{-}\phantom{.}\phantom{.}0.0\phantom{\pm}\phantom{0}\phantom{.}\phantom{0}\phantom{0}\phantom{0}$ | $0.020\pm 0.009$ | 2418.8L
|  |  | 2013.038 | 56307 | $0.071\pm 0.010$ | $0.084\pm 0.030$ | $-165.0\pm\phantom{0}5.8$ | $0.081\pm 0.022$ | 81.2
|  |  | 2013.038 | 56307 | $0.044\pm 0.012$ | $0.391\pm 0.099$ | $-157.1\pm\phantom{0}4.6$ | $0.183\pm 0.037$ | 9.9
|  |  | 2013.038 | 56307 | $0.025\pm 0.014$ | $0.737\pm 0.296$ | $-176.0\pm\phantom{0}5.8$ | $0.329\pm 0.057$ | 1.7
|  |  | 2013.038 | 56307 | $0.031\pm 0.015$ | $2.396\pm 0.443$ | $-177.8\pm\phantom{0}2.7$ | $0.512\pm 0.067$ | 0.9
0219+428 | 0.792 | 0.035 | 2013.288 | 56398 | $0.203\pm 0.013$ | $0.000\phantom{\pm}\phantom{0}\phantom{.}\phantom{0}\phantom{0}\phantom{0}$ | $\phantom{-}\phantom{.}\phantom{.}0.0\phantom{\pm}\phantom{0}\phantom{.}\phantom{0}\phantom{0}\phantom{0}$ | $0.053\pm 0.013$ | 542.0
|  |  | 2013.288 | 56398 | $0.047\pm 0.011$ | $0.149\pm 0.060$ | $-161.5\pm\phantom{0}6.8$ | $0.122\pm 0.029$ | 23.7
|  |  | 2013.288 | 56398 | $0.030\pm 0.011$ | $0.448\pm 0.098$ | $-160.4\pm\phantom{0}3.8$ | $0.148\pm 0.036$ | 10.3
|  |  | 2013.288 | 56398 | $0.021\pm 0.015$ | $0.843\pm 0.613$ | $-169.4\pm 11.1$ | $0.558\pm 0.077$ | 0.5
|  |  | 2013.288 | 56398 | $0.025\pm 0.014$ | $2.474\pm 0.372$ | $-177.3\pm\phantom{0}2.2$ | $0.398\pm 0.062$ | 1.2

Note. — Values of $T_{\text{b,obs}}$ denoted by the letter “L” represent the
lower limits to the brightness temperature; see $\S$2.1 for details.
(Table 2 (10192 lines) is published in its entirety in the machine-readable
format. A portion is shown here for guidance regarding its form and content.)

The uncertainties of the model parameters were calculated using the formalism
described in Jorstad et al. (2017), which is based on an empirical relation
between the uncertainties and the brightness temperatures of knots,
$T_{b,\text{obs}}=7.5\times 10^{8}S/a^{2}$ K (e.g., Jorstad et al., 2005;
Casadio et al., 2015). Estimates of the errors in each model parameter are as
follows: $\sigma_{\text{X}}\approx 1.3\times 10^{4}\ T_{b,\text{obs}}^{-0.6}$,
$\sigma_{\text{Y}}\approx 2\sigma_{\text{X}}$, $\sigma_{S}\approx 0.09\
T_{b,\text{obs}}^{-0.1}$, and $\sigma_{a}\approx 6.5\
T_{b,\text{obs}}^{-0.25}$. Here, $\sigma_{\text{X}}$ and $\sigma_{\text{Y}}$
are the $1\sigma$ uncertainties in right ascension and declination in mas,
$\sigma_{S}$ the uncertainty in the flux density in Jy, and $\sigma_{a}$ the
uncertainty in component size in mas. Note that the typical north-south
synthesized beam is $\sim 2\times$ the typical east-west beam size, which
increases the uncertainty in declination. As in Jorstad et al. (2017), we have
also added a minimum positional error of 0.005 mas (related to the resolution
of the observations) and a typical amplitude calibration error of 5% to these
uncertainties.

Table 2 gives the total intensity jet models and brightness temperature values
for all components in the jet of each source from 2013 January to 2018
December as follows: 1—source name according to its B1950 coordinates;
2—$\chi^{2}$ value of the model fit; 3—root-mean-square difference between the
observed and model visibilities; 4—epoch of the start of the observation;
5—MJD of the start of the observation; 6—the flux density, $S$, and its
uncertainty $\sigma_{S}$, in Jy; 7—distance with respect to the core,
$R_{\text{obs}}$, and its uncertainty $\sigma_{R}$, in mas; 8—Position angle
of the knot with respect to the core, $\Theta_{\text{obs}}$, and its
uncertainty $\sigma_{\Theta}$, measured north through east in degrees;
9—angular size of the knot, $a$, and its uncertainty $\sigma_{a}$, in mas; and
10—observed brightness temperature $T_{b,\text{obs}}$, in units of $10^{8}$ K.
The core at each epoch is located at position $(X,Y)=(0,0)$, where $X$ is the
relative right ascension and $Y$ is the relative declination, and $X$ and $Y$
are related to $R_{\text{obs}}$ and $\Theta_{\text{obs}}$ through
$R_{\text{obs}}=\sqrt{X^{2}+Y^{2}}$ and $\tan{\Theta_{\text{obs}}}=X/Y$. For
9.8% of knots, the model fit yielded a size less than 0.02 mas. This angular
size is too small to be resolved on the longest VLBA baselines. Thus, lower
limits to $T_{b,\text{obs}}$ for these knots are calculated using $a=0.02$
($\sim 20\%$ of the resolution of the longest baselines) and denoted by “L” in
column 10 of Table 2. The total intensity jet models and brightness
temperature values for all components in the jet of each source from 2007 June
to 2012 December can be found in Table 2 of Jorstad et al. (2017).

###  2.2 Identification of Components and Calculation of Apparent Speed

As discussed above, the VLBI “core” is defined as the bright, compact (either
unresolved or partially resolved) emission feature at the upstream end of the
jet, relative to which the positions of all other components are measured. To
identify components, we make use of the fact that all four model parameters
$(S,R_{\text{obs}},\Theta_{\text{obs}}$, and $a$) should not change abruptly
with time given the regularity of our observations. If a knot is identified at
$N\geq 4$ epochs, it is assigned a kinematic classification and identification
number (ID).

Previous studies (e.g., Homan et al., 2015; Jorstad et al., 2017) have found
that a given knot can exhibit both acceleration and deceleration as it
propagates down the jet. In order to describe the locations along the jet of
such regions of changing velocity, we have created a new formalism for the
calculation of kinematic parameters of a knot, utilizing piece-wise linear
fits to separation (from the core) versus time data instead of polynomial fits
as used in Jorstad et al. (2017). This simplifies the calculation of kinematic
parameters by describing the motion in terms of the mean velocity of each
segment of time within which the velocity is roughly constant and
significantly different from that of other segments. This allows us to
estimate the distance down the jet at which acceleration occurs. The new
procedure consists of the following steps:

  1. 1.

_Piece-wise Linear Fits to $X$ and $Y$:_ We start with the assumption that
knots detected at $N\geq 4$ epochs can have their $X$ and $Y$ positions fit by
a linear trend of the form

$X(t_{i})=a_{0}+a_{1}(t_{i}-t_{\text{mid}})$ (1)

and

$Y(t_{i})=b_{0}+b_{1}(t_{i}-t_{\text{mid}}),$ (2)

where $t_{i}$ is the epoch of observation, $i=1,\ldots,N$, and
$t_{\text{mid}}\equiv(t_{N}+t_{1})/2$. We use the statsmodels
package555https://www.statsmodels.org/stable/index.html (Seabold & Perktold,
2010) weighted least-squares program to calculate the best-fit parameters and
reduced $\chi^{2}$ value of the fit. These fit $\chi^{2}$ values are then
compared to the critical $\chi^{2}_{\text{crit}}$, corresponding to a
significance level of $\zeta=0.05$ for $f=N-2$ degrees of freedom (Bowker &
Lieberman, 1972). If $\chi^{2}<\chi^{2}_{\text{crit}}$, our assumption was
valid. However, if the opposite is true, the knot may be experiencing
acceleration in the jet and a more complicated fit is necessary. In such
cases, we use the pwlf
package666https://jekel.me/piecewise_linear_fit_py/index.html to perform a
weighted least-squares piece-wise linear fit to the data (Jekel & Venter,
2019). The advantage of such a fit over higher-order polynomial fits is that
it relies on simple linear fits for segments of the data, so that there are no
concerns about over-fitting the data or estimating kinematic parameters from
higher-order fits.

We define the piece-wise fits in this work in terms of the number of line
segments present for each coordinate of motion $X$ and $Y$. A single-segment
fit describes the simple linear case as defined above. A two-segment fit
indicates that the motion can be described by two different speeds with a
single acceleration region at a particular break point, $t_{b}$, and a three-
segment fit indicates three speeds and two acceleration regions and
breakpoints. In a two-segment fit, there are 5 free parameters, while in a
three-segment fit there are 8 free parameters. A piece-wise fit with the
fewest number of segments for which $\chi^{2}\leq\chi^{2}_{\text{crit}}$ is
used to fit the data, where $\chi^{2}_{\text{crit}}$ is the $\chi^{2}$ value
corresponding to a significance level of $\zeta=0.05$ for $f=N-\ell$, where
$\ell$ is the number of fit parameters. However, a two-segment fit is only
used if $N>8$, and a three-segment fit is only used if $N>15$, given the
number of degrees of freedom. In cases for which the $\chi^{2}$ criterion is
not reached for even a three-segment piece-wise fit, we use the simplest fit
that adequately describes the data ($\chi^{2}$ closest to
$\chi^{2}_{\text{crit}}$). There are only several such cases. In some cases
the $X$ and $Y$ motion breakpoints can overlap in time, generally defined as
$|t_{b,X}-t_{b,Y}|\leq 0.25$ yr. In these situations we average $t_{b,X}$ and
$t_{b,Y}$ and use this average time as a single break-point. The model fit is
then re-run with the break times frozen to the averaged values. The overlap
value of 0.25 yr was chosen because in any time interval $\Delta t=0.25$ yr,
our roughly monthly observations will lead to only 2 or 3 observations within
the time interval so that a calculation of the speed of such a segment is
unreliable. The uncertainties of all best-fit parameters are calculated from
the diagonals of the co-variance matrices of the fits.

  2. 2.

_Calculation of the Epoch of Ejection, $t_{\circ}$:_ The epoch of ejection is
defined as the time when the centroid of a knot coincides with the centroid of
the core, with the knot motion extrapolated back to $(X,Y)=(0,0)$. In general,
all observations of the motion of a knot are important for calculating its
kinematic parameters. However, for determining $t_{\circ}$, the most important
observations are when the knot is closest to the core, before any potential
acceleration. We thus extrapolate the linear fits to $X$ and $Y$ back to
$(0,0)$ from the first identified break point (or use the first ten
observations of a knot in the case of no acceleration). This extrapolation
provides the time of ejection along each coordinate, $t_{x\circ}$ and
$t_{y\circ}$. We then calculate the true time of ejection as:

$t_{\circ}=\frac{t_{\mathrm{x\circ}}/\sigma_{t_{\mathrm{x\circ}}}^{2}+t_{\mathrm{y\circ}}/\sigma_{t_{\mathrm{y\circ}}}^{2}}{1/\sigma_{t_{\mathrm{x\circ}}}^{2}+1/\sigma_{t_{\mathrm{y\circ}}}^{2}}$
(3)

Here, the $1\sigma$ propagated uncertainties in $t_{x\circ}$ and $t_{y\circ}$
are defined as $\sigma_{t_{x\circ}}$ and $\sigma_{t_{y\circ}}$, respectively.
The uncertainty in $t_{\circ}$ is then calculated as

$\sigma_{t_{\circ}}=\sqrt{\frac{(t_{\circ}-t_{\mathrm{x\circ}})^{2}/\sigma_{t_{\mathrm{x\circ}}}^{2}+(t_{\circ}-t_{\mathrm{y\circ}})^{2}/\sigma_{t_{\mathrm{y\circ}}}^{2}}{1/\sigma_{t_{\mathrm{x\circ}}}^{2}+1/\sigma_{t_{\mathrm{y\circ}}}^{2}}}$
(4)

For a few cases, weighting $t_{x\circ}$ and $t_{y\circ}$ by uncertainties does
not provide a robust solution since they are significantly different from each
other, a situation that we define as $|t_{x\circ}-t_{y\circ}|\geq 0.5$ yr.
This often occurs in the case of knot motion directly along the $X$ or $Y$
axis or for knots whose ejections occurred before our monitoring. In such
cases, we construct a best-fit linear trend to the $R$ coordinate (either the
first ten epochs of observation in the case of purely linear motion or the
first segment of motion for knots showing acceleration) and find $t_{\circ}$
as its root, with its corresponding uncertainty.

To each uncertainty in $t_{\circ}$, regardless of calculation method, we add a
minimum uncertainty of $0.005/\mu$ yr, where 0.005 mas is the size of half a
pixel in the VLBA images and $\mu$ is the proper motion of the component
(defined below), in mas yr-1. This minimum error is typically $\sim 3$-$4$
days and is generally much less than the error computed as indicated above.

  3. 3.

_Calculation of Apparent Speeds:_ Based on the piece-wise linear fits that
best fit the $X$ and $Y$ data, we calculate the proper motion $\mu$, direction
of motion $\Phi$, and apparent speed $\beta_{\mathrm{app}}$ for each time
segment. We continue to define the speed and direction of the motion as in
Jorstad et al. (2005), with
$\langle\mu_{j}\rangle=\sqrt{\langle\mu_{x,j}\rangle^{2}+\langle\mu_{y,j}\rangle^{2}}$
and
$\langle\Phi_{j}\rangle=\tan^{-1}(\langle\mu_{x,j}\rangle/\langle\mu_{y,j}\rangle)$,
where the subscript $j$ indicates the number of the line segment. As the
motion in each coordinate $X$ and $Y$ can be described independently by up to
3 linear segments, the maximum range of $j=1,2,\ldots,5$. We arrange the
break-points chronologically.

  4. 4.

_Calculation of Apparent Accelerations:_ If $j\neq 1$, there is a change in
the speed and direction of motion of the knot over time as it moves down the
jet. In the case of two or more line segments, we define a region in time
between the mid-points of the line segments as the “acceleration region,”
$\Delta t$. While the acceleration may occur over a shorter or longer time
period than this mid-point definition, using the mid-points of a segment to
determine $\Delta t$ characterizes the acceleration without _a priori_
knowledge of the cause of acceleration. The acceleration along each coordinate
$X$ and $Y$ is defined as $\dot{\mu}_{x}=\Delta\mu_{x}/\Delta t$ and
$\dot{\mu}_{y}=\Delta\mu_{y}/\Delta t$. Then the total acceleration parallel,
$\dot{\mu}_{\parallel}$, and perpendicular, $\dot{\mu}_{\perp}$, to the jet
axis can be calculated through
$\dot{\mu}_{\parallel}=\dot{\mu}_{x}\sin{\langle\Theta_{\text{jet}}\rangle}+\dot{\mu}_{y}\cos{\langle\Theta_{\text{jet}}\rangle}$
and
$\dot{\mu}_{\perp}=\dot{\mu}_{x}\cos{\langle\Theta_{\text{jet}}\rangle}-\dot{\mu}_{y}\sin{\langle\Theta_{\text{jet}}\rangle}$,
where $\langle\Theta_{\text{jet}}\rangle$ is the average jet direction. How we
calculate the average jet direction for each source is discussed in $\S$3.2.

To interpret the results of the jet kinematics in a uniform manner across ten
years of observation, we have applied the above method to all of the knots
presented in Jorstad et al. (2017), in addition to the new knots observed
since 2012 December. As an example of this fitting procedure, Figure 1 shows
the result for knot _C35_ of the FSRQ 1253$-$055\. Panel (b) and (c) presents
the linear fits to the separation of the knot versus time along the $X$ and
$Y$ dimensions, respectively. In this example, 2 segments are needed to
represent both the $X$ and $Y$ motions accurately. Panel (a) shows the motion
of the knot relative to the core along the radial $R$ dimension, combining the
$X$ and $Y$ models. For this example, there are three measurements of the
speed of the knot and two estimates for the acceleration.

Figure 1: An example of the piece-wise linear fitting procedure described in
$\S$2.2 for the knot _C35_ of the FSRQ 1253–055. The three panels represent
the total distances from the core along the (a) radial _R_ -axis, (b) _X_
-axis (relative right ascension), and (c) _Y_ -axis (relative declination).
The red lines represent the piece-wise linear fits to the data, while gray
dotted lines indicate break-point dates. The fit line in (a) is constructed
from the combination of the fits in (b) and (c). The estimated epoch of
ejection, $t_{\circ}$, is shown in (a) as a horizontal red line. See text for
details.

##  3 Parsec-Scale Jet Structure

Given the above determination of the kinematic parameters of the jet
components, we have separated the knots into different types according to
their properties. Images of all sources at all epochs possess a core
(described above), which we label as _A0_. All remaining components are
classified based on their proper motion. Those with significant proper motion
($\mu\geq 2\sigma_{\mu}$) are designated as moving features of type _M_.
Components without this proper motion are termed “quasi-stationary” (often
shortened to “stationary” below) features of type _St_. We continue the naming
scheme introduced in Jorstad et al. (2017) for all previously analyzed
sources. For the two sources added to the sample later (the BLs 1652+398 and
1959+650), designations start at 1 for both moving and quasi-stationary
features.

We have calculated the average parameters of all features in all of the
blazars, and list the results in Table 3 in the following format: 1—B1950 name
of the source; 2—ID of the knot; 3—the number of epochs, $N$, at which the
knot is detected; 4—mean flux density, $\langle S\rangle$, of the knot and its
standard deviation, in Jy; 5—average distance of the knot, $\langle R\rangle$,
and its standard deviation, in mas777This average represents the typical
distance at which the knot was observed and is not necessarily the mid-point
of its trajectory.; 6—average position angle of the knot with respect to the
core, $\langle\Theta\rangle$, and its standard deviation, in degrees. For the
core, $\langle\Theta\rangle$ is calculated using the weighted average of the
$\Theta$ of the all knots in the jet over all epochs at which $R\leq 0.7$ mas,
and thus represents the average projected direction of the inner jet at a
given epoch. Counter-jet features are not included in this calculation; 7—mean
size of the knot, $\langle a\rangle$, and its standard deviation, in mas; and
8—kinematic type of the knot (see above).

Table 3: Jet Structure

Source | Knot | $N$ | $\langle S\rangle$ | $\langle R\rangle$ | $\langle\Theta\rangle$ | $\langle a\rangle$ | Type
---|---|---|---|---|---|---|---
|  |  | [Jy] | [mas] | [deg] | [mas] |
0219+428 | A0 | 62 | $0.24\pm 0.09$ | $0.00\pm 0.02$ | $-166.0\pm 10.0$ | $0.04\pm 0.02$ | Core
| A1 | 52 | $0.08\pm 0.04$ | $0.17\pm 0.04$ | $-167.9\pm 10.9$ | $0.09\pm 0.03$ | St
| A2 | 36 | $0.04\pm 0.02$ | $0.37\pm 0.06$ | $-168.5\pm\phantom{0}9.5$ | $0.14\pm 0.04$ | St
| A3 | 50 | $0.03\pm 0.01$ | $0.65\pm 0.09$ | $-173.2\pm\phantom{0}5.8$ | $0.21\pm 0.09$ | St
| A4 | 56 | $0.03\pm 0.01$ | $2.31\pm 0.09$ | $-177.1\pm\phantom{0}1.4$ | $0.46\pm 0.14$ | St
| B1 | 6 | $0.02\pm 0.01$ | $1.03\pm 0.19$ | $-175.7\pm\phantom{0}1.3$ | $0.27\pm 0.06$ | M
| B2 | 10 | $0.02\pm 0.01$ | $1.04\pm 0.22$ | $-172.7\pm\phantom{0}6.0$ | $0.40\pm 0.10$ | M
| B3 | 6 | $0.04\pm 0.03$ | $0.48\pm 0.26$ | $-160.7\pm 18.3$ | $0.15\pm 0.11$ | M
| B4 | 7 | $0.01\pm 0.01$ | $1.35\pm 0.45$ | $-168.6\pm\phantom{0}1.3$ | $0.41\pm 0.10$ | M
| B5 | 11 | $0.02\pm 0.02$ | $0.94\pm 0.56$ | $-169.5\pm\phantom{0}3.2$ | $0.28\pm 0.14$ | M
| B6 | 4 | $0.01\pm 0.01$ | $1.40\pm 0.24$ | $-177.5\pm\phantom{0}3.2$ | $0.43\pm 0.06$ | M
0235+164aaThe direction of the inner jet for 0235+164 is uncertain, since features detected in the jet have a range of P.A. exceeding 180$\arcdeg$. | A0 | 126 | $1.30\pm 0.83$ | $0.00\pm 0.02$ | $\phantom{-000.0}\ldots$ | $0.06\pm 0.03$ | Core
| B1 | 6 | $0.16\pm 0.23$ | $0.40\pm 0.17$ | $-\phantom{0}17.7\pm\phantom{0}5.8$ | $0.26\pm 0.13$ | M
| B2 | 15 | $1.03\pm 0.88$ | $0.22\pm 0.08$ | $\phantom{-}163.1\pm 12.2$ | $0.15\pm 0.06$ | M
| B3 | 40 | $0.25\pm 0.14$ | $0.21\pm 0.07$ | $\phantom{-}164.9\pm\phantom{0}6.4$ | $0.13\pm 0.05$ | M
| B4 | 16 | $0.15\pm 0.11$ | $0.19\pm 0.06$ | $-\phantom{0}63.4\pm 21.0$ | $0.13\pm 0.10$ | M
| B5 | 17 | $0.41\pm 0.29$ | $0.28\pm 0.10$ | $\phantom{-0}59.1\pm\phantom{0}4.0$ | $0.20\pm 0.07$ | M
| B6 | 19 | $0.19\pm 0.16$ | $0.27\pm 0.11$ | $\phantom{-0}51.2\pm\phantom{0}6.2$ | $0.23\pm 0.14$ | M
| B7 | 6 | $0.04\pm 0.01$ | $0.26\pm 0.03$ | $\phantom{-00}9.6\pm\phantom{0}5.5$ | $0.11\pm 0.05$ | M

Note. — (Table 3 (559 lines) is published in its entirety in the machine-
readable format. A portion is shown here for guidance regarding its form and
content.)

Figure 2 presents an image of each source at a single epoch between 2013
January and 2018 December, when the most prominent features in the jet were
visible. The black circles indicate the sizes and positions of the features
according to the model parameters for that epoch. These images do not contain
labels for all observed features in the jet over the 10 years of observations.
We list the important parameters for each image in Table 4 as follows: 1—B1950
name of the source; 2—epoch of observation; 3—total intensity peak of the map,
$I_{\text{peak}}$, in mJy beam-1; 4—lowest shown contour,
$I_{\text{low}}^{\text{cnt}}$, in mJy beam-1; 5—size of the restoring beam,
used to produce images of similar fidelity across all epochs, in units of mas
$\times$ mas, with all restoring beam orientation angles set to $-10\arcdeg$;
6—antennas with no data for the source at that particular observation epoch
and 7—size and orientation angle of the synthesized beam corresponding to the
$uv$-coverage for the epoch of observation, in units of mas $\times$ mas and
degrees.

Figure 2: Total intensity images of the VLBA-BU-BLAZAR sample at 43 GHz with uniform weighting. Parameters of all images are listed in Table 4. Contours decrease by a factor of $\sqrt{2}$ from the peak flux density. Black circles on the images indicate the jet features according to the model fits. The shaded ellipse in the bottom left corner represents the restoring beam.  Figure 2: (Continued.) Figure 2: (Continued.) Table 4: Parameters of VLBA 43 GHz Images in Figure 2. Source | Epoch | $I_{\text{peak}}$ | $I_{\text{low}}^{\text{cnt}}$ | Beam Size | Unavailable | Synthesized Beam Size
---|---|---|---|---|---|---
|  |  |  |  | AntennasaaAntenna abbreviations can be found at the NRAO VLBA website (https://science.nrao.edu/facilities/vlba/docs/manuals/oss/sites). Individual antennas are sometimes unavailable due to weather, maintenance, or undetected fringes (in the case of the weakest sources, e.g., 1959+650). |
|  | [mJy beam-1] | [mJy beam-1] | [mas $\times$ mas] |  | [mas $\times$ mas, deg]
0219+428 | 2015 Dec 05 | 233 | 1.82 | $0.17\times 0.34$ | $\ldots$ | $0.16\times 0.28,\ -\phantom{0}4.1$
0235+164 | 2017 Mar 19 | 866 | 3.38 | $0.15\times 0.36$ | $\ldots$ | $0.16\times 0.33,\ -\phantom{0}3.4$
0316+413 | 2018 Aug 26 | 1572 | 12.28 | $0.15\times 0.28$ | $\ldots$ | $0.16\times 0.30,\ \phantom{-}10.0$
0336–019 | 2017 May 13 | 1214 | 6.71 | $0.16\times 0.38$ | $\ldots$ | $0.16\times 0.37,\ -\phantom{0}2.8$
0415+379 | 2016 Jan 01 | 578 | 2.26 | $0.16\times 0.32$ | $\ldots$ | $0.16\times 0.33,\ -\phantom{0}2.5$
0420–014 | 2014 Sep 23 | 1207 | 4.71 | $0.15\times 0.38$ | HN | $0.16\times 0.39,\ -\phantom{0}2.9$
0430+052 | 2018 Apr 19 | 1117 | 2.18 | $0.14\times 0.34$ | $\ldots$ | $0.15\times 0.36,\ -\phantom{0}4.4$
0528+134 | 2015 Jul 02 | 366 | 0.72 | $0.15\times 0.33$ | $\ldots$ | $0.16\times 0.37,\ -\phantom{0}0.9$
0716+714 | 2016 Jan 31 | 1498 | 8.28 | $0.15\times 0.24$ | $\ldots$ | $0.16\times 0.28,\ -\phantom{0}7.1$
0735+178 | 2017 Nov 06 | 204 | 2.25 | $0.21\times 0.35$ | BR, SC | $0.16\times 0.73,\ -20.0$
0827+243 | 2015 Aug 01 | 322 | 1.26 | $0.15\times 0.31$ | $\ldots$ | $0.16\times 0.31,\ -\phantom{0}3.1$
0829+046 | 2015 Jun 09 | 723 | 2.83 | $0.19\times 0.35$ | $\ldots$ | $0.17\times 0.39,\ -\phantom{0}7.8$
0836+710 | 2016 Apr 22 | 854 | 3.33 | $0.17\times 0.24$ | $\ldots$ | $0.16\times 0.24,\ -16.0$
0851+202 | 2015 Feb 14 | 4807 | 9.39 | $0.16\times 0.33$ | FD | $0.20\times 0.57,\ \phantom{-}26.0$
0954+658 | 2015 Sep 22 | 723 | 3.99 | $0.15\times 0.24$ | $\ldots$ | $0.17\times 0.24,\ \phantom{-}\phantom{0}8.5$
1055+018 | 2018 Jun 16 | 2195 | 8.57 | $0.14\times 0.34$ | SC | $0.15\times 0.68,\ -18.0$
1101+384 | 2014 Feb 24 | 281 | 3.10 | $0.15\times 0.24$ | FD | $0.16\times 0.28,\ -21.0$
1127–145 | 2013 Dec 16 | 1112 | 4.34 | $0.16\times 0.39$ | KP | $0.17\times 0.46,\ -\phantom{0}8.2$
1156+295 | 2014 Jul 28 | 1087 | 4.25 | $0.15\times 0.35$ | $\ldots$ | $0.16\times 0.31,\ \phantom{-}\phantom{0}0.5$
1219+285 | 2015 May 11 | 205 | 0.80 | $0.15\times 0.35$ | $\ldots$ | $0.17\times 0.30,\ -\phantom{0}8.3$
1222+216 | 2018 Dec 08 | 473 | 1.85 | $0.16\times 0.36$ | $\ldots$ | $0.16\times 0.33,\ -10.0$
1226+023 | 2014 Nov 15 | 3612 | 7.05 | $0.15\times 0.38$ | $\ldots$ | $0.18\times 0.44,\ -\phantom{0}5.6$
1253–055 | 2016 Mar 18 | 6266 | 24.47 | $0.15\times 0.36$ | $\ldots$ | $0.17\times 0.40,\ -\phantom{0}4.6$
1308+326 | 2014 Nov 15 | 571 | 2.23 | $0.15\times 0.29$ | $\ldots$ | $0.17\times 0.35,\ -15.0$
1406–076 | 2016 Oct 06 | 654 | 1.28 | $0.19\times 0.41$ | HN, MK | $0.20\times 0.82,\ \phantom{-}23.0$
1510–089 | 2016 Jan 01 | 4277 | 5.91 | $0.19\times 0.41$ | $\ldots$ | $0.14\times 0.45,\ -10.0$
1611+343 | 2018 Jun 16 | 616 | 1.70 | $0.15\times 0.36$ | SC | $0.16\times 0.37,\ -27.0$
1622–297 | 2015 May 11 | 589 | 3.25 | $0.16\times 0.58$ | HN | $0.13\times 0.40,\ -\phantom{0}4.8$
1633+382 | 2014 May 03 | 2230 | 4.36 | $0.15\times 0.33$ | NL | $0.15\times 0.30,\ -\phantom{0}9.2$
1641+399 | 2015 Dec 05 | 1800 | 1.76 | $0.14\times 0.24$ | $\ldots$ | $0.15\times 0.28,\ -17.0$
1652+398 | 2017 Feb 04 | 175 | 0.97 | $0.15\times 0.33$ | BR | $0.15\times 0.31,\ -20.0$
1730–130 | 2016 Nov 28 | 1812 | 3.54 | $0.14\times 0.45$ | $\ldots$ | $0.14\times 0.38,\ -\phantom{0}4.5$
1749+096 | 2015 Jun 09 | 3310 | 3.23 | $0.14\times 0.34$ | $\ldots$ | $0.15\times 0.36,\ -\phantom{0}7.0$
1959+650 | 2018 Jul 16 | 126 | 1.39 | $0.17\times 0.24$ | HN, PT, SC | $0.20\times 0.34,\ -\phantom{0}7.0$
2200+420 | 2018 Jul 16 | 1745 | 1.70 | $0.15\times 0.29$ | $\ldots$ | $0.16\times 0.25,\ -\phantom{0}8.1$
2223–052 | 2018 Feb 17 | 486 | 1.90 | $0.16\times 0.38$ | MK, SC | $0.27\times 0.81,\ -\phantom{0}5.0$
2230+114 | 2016 Jan 01 | 1840 | 3.59 | $0.15\times 0.35$ | $\ldots$ | $0.14\times 0.40,\ -\phantom{0}7.9$
2251+158 | 2016 Jun 10 | 5115 | 28.26 | $0.14\times 0.33$ | $\ldots$ | $0.13\times 0.33,\ -\phantom{0}6.1$

###  3.1 Observed Brightness Temperatures

The observed brightness temperatures of all jet features from 2007 June to
2012 December are listed in Table 2 of Jorstad et al. (2017), while those from
2013 January to 2018 December are listed in Table 2 of this paper. Combined,
there are a total of 19,644 estimates for the brightness temperatures of
components. Of these, 10.1% of the $T_{\text{b,obs}}$ values are lower limits
(see $\S$2.1).

Figure 3: Host-galaxy frame brightness temperature distributions of the core
(left) and non-core (knots, right) features between 2013 January and 2018
December. Unshaded regions of the distributions represent components with
lower limits to $T_{\mathrm{b,obs}}^{\mathrm{s}}$. The bin sizes are neither
linear nor logarithmic. Values in parentheses are the number of
$T_{\text{b,obs}}$ measurements of components for each subclass, including
lower limits.

For a direct comparison with the brightness temperatures analyzed in Jorstad
et al. (2017), Figure 3 (left) displays the brightness temperatures of all
cores in the host galaxy frame
($T_{\mathrm{b,obs}}^{\mathrm{s}}=T_{\mathrm{b,obs}}(1+z)$) from 2013 January
to 2018 December for each subclass of blazar. The RG distribution peaks at
$(1-5)\times 10^{11}$ K. Only a few cores have epochs when
$T_{\mathrm{b,obs}}^{\mathrm{s}}\geq 5\times 10^{11}$ K. The FSRQ distribution
is bimodal, with peaks at both $(1-5)\times 10^{11}$ K and $(1-5)\times
10^{12}$ K. The highest peak of the BL distribution is also at $(1-5)\times
10^{12}$ K, but includes a large number of lower limits. The BL distribution
also has the most significant number of cores with low brightness temperatures
$(T_{\mathrm{b,obs}}^{\mathrm{s}}<10^{11}$ K). The distributions of the
brightness temperatures for knots, other than the cores, between 2013 January
and 2018 December are shown in Figure 3 (right). In general, these knots have
lower brightness temperatures than the cores, but with much wider dispersion
about the mean. Both distributions are similar in shape to those presented in
Jorstad et al. (2017), indicating no apparent time variability in the
populations of brightness temperatures of components as a whole for the
sources in our sample.

A direct comparison of the distributions for each subclass presented in Figure
3 (as with a Kolmogorov-Smirnov (KS) Test) would ignore two crucial aspects of
the data: (1) that the brightness temperature values for the cores and the
knots of a particular object in our sample are correlated by virtue of coming
from the same source, rather than being an independent measurement from a
random blazar, and (2) the relatively large number (in some cases, $\sim
50\%$, e.g., the core of the BL 1959+650) of lower limits present in the
sample. In order to compare the brightness temperatures of the cores and knots
in our sample between each subclass, we derive “typical” values for the
brightness temperatures of components for each source by estimating the
survival function of the brightness temperature of the cores and knots to get
the median value (see Appendix A).

Figure 4 (left) shows the median brightness temperatures of core components of
blazars in our sample between 2013 January and 2018 December, while Fig. 4
(right) shows the same for the entire observing period of the program (2007
June - 2018 December). The three RGs in the sample have very similar median
core brightness temperatures, $\sim 5\times 10^{11}$ K. The distribution of
the median $T_{\mathrm{b,obs}}^{\mathrm{s}}$ for BLs is very wide-spread,
ranging from $10^{10}$ K to $\sim 10^{13}$ K. The cores of FSRQs are more
narrowly distributed around a central brightness temperature, $\sim 10^{12}$
K. As the process of calculating the median brightness temperature has already
taken into account the presence of lower-limits, we can now meaningfully
compare the distributions for each subclass through a KS test. Values of the
test-statistic, $\mathcal{D}$, and significance, $p$, values for the test
between each subclass distribution are given in Table 5. The differences
between the FSRQ and RG core distributions are statistically significant at a
$\sim 3\sigma$ level, with $p=0.009$, once all cores from 2007 June to 2018
December have been taken into account. The FSRQ and BL cores also appear to be
different, but at a less significant level ($p=0.086$). There is no
statistical difference between the BL and RG distributions, but it is
important to note that the number of sources is small (13 for BLs and 3 for
RGs).

Figure 4: Median core component host-galaxy frame brightness temperature
distributions of blazars in the VLBA-BU-BLAZAR sample, as determined through
survival analysis. The bin sizes are neither linear nor logarithmic.  Table 5:
KS test statistic, $\mathcal{D}$-, and $p$-values for the median brightness
temperatures of jet cores.

Subclass 1 | Subclass 2 | $\mathcal{D}$ | $p$
---|---|---|---
FSRQ | BL | 0.388 | 0.129
FSRQ | RG | 0.818 | 0.030
BL | RG | 0.615 | 0.229
FSRQ | BL | 0.416 | 0.086
FSRQ | RG | 0.909 | 0.009
BL | RG | 0.538 | 0.375

Note. — The data for section 1 of this table consist of observations taken
from 2013 January through 2018 December. The data for section 2 were taken
from 2007 June through 2018 December.

Similarly to the core distributions, in Figure 5 we show the median knot
brightness temperature distributions for each subclass, including knots
observed from 2013 January to 2018 December (left) and 2007 June to 2018
December (right). In general, the median knot brightness temperature
distribution for each subclass is centered on lower values than that of the
core components. The knots of FSRQs are more concentrated around the mean
value than the knots of BLs. The knots of the RG 0316+428 have a higher median
$T_{\mathrm{b,obs}}^{\mathrm{s}}$ than the other two RG in the sample. The
$\mathcal{D}$\- and $p$-values for a KS test between the knot distributions
are given in Table 6. The only distributions that show a statistically
significant difference are the knots of FSRQs compared with BLs, at a
significance $p=0.002$ for the entire observed time period.

Figure 5: Median knot component host-galaxy frame brightness temperature
distributions of blazars in the VLBA-BU-BLAZAR sample, as determined through
survival analysis. The bin sizes are neither linear nor logarithmic.  Table 6:
KS test statistic, $\mathcal{D}$-, and $p$-values for the median brightness
temperatures of jet knots.

Subclass 1 | Subclass 2 | $\mathcal{D}$ | $p$
---|---|---|---
FSRQ | BL | 0.538 | 0.010
FSRQ | RG | 0.363 | 0.750
BL | RG | 0.615 | 0.229
FSRQ | BL | 0.615 | 0.002
FSRQ | RG | 0.439 | 0.574
BL | RG | 0.615 | 0.229

Note. — The data for section 1 of this table consist of observations taken
from 2013 January through 2018 December. The data for section 2 were taken
from 2007 June through 2018 December.

Despite the low significance levels according to a KS test (likely due to the
number of sources in our sample), we see a general trend across the median
core and knot brightness temperature distributions: _FSRQs and BLs have more
intense cores than those of RGs, but the knots in the extended jets of BLs are
significantly less intense than the knots of RGs and FSRQs._ A small fraction
of analyzed knots are weak diffuse features with
$T_{\text{b,obs}}^{\text{s}}<10^{8}$ K. The highest brightness temperatures,
near $10^{13}$ K, are discussed in $\S$6.1.

###  3.2 Jet Position Angles

With over 10 years of observations, we can constrain the bulk characteristics
of the jet, such as the average position angle,
$\langle\Theta_{\text{jet}}\rangle$, and its temporal evolution. Such
constraints are important for determining whether there are sudden changes in
the direction of the jet, or if precession is common in blazars. In addition,
since the calculation of acceleration along or perpendicular to the jet axis
depends on the position angle, such information can allow us to account for
the effect of changes in $\Theta_{\text{jet}}$ on the accelerations.

In order to determine how the jet position angles vary, we have constructed
plots of the observed position angle, $\Theta$, versus time for all knots in
each source. Figure 6 shows an example of such plots for the source 0219+428.
Figure SET 6 provides a similar plot for each source in the sample.

Fig. Set 6. Jet Position Angles

Figure 6: Jet position angles of each knot, relative to the core, of the BL
0219+428 (3C 66A). The colors of each knot are the same as in Figure 7, and
are designated in the upper left corner. We classify this source as having a
constant jet position angle according to our criterion (see the text and Table
7). The black line indicates $\langle\Theta_{\text{jet}}\rangle$, and the two
gray dashed lines show the standard deviation of the PA with respect to
$\Theta_{\text{jet}}$, $\sigma_{\langle\Theta_{\text{jet}}\rangle}$, which is
listed in column 7 of Table 7. Similar plots for all sources in the sample are
available in the figure set.
(The complete figure set (38 images) is available online).

We have visually inspected the plots in Figure SET 6 to see if there are any
overall trends or abrupt changes in the position angles of knots. This leads
us to classify the average jet position angle behavior into three categories:
1) $\langle\Theta_{\text{jet}}(t)\rangle=$ constant (i.e., variations are
insignificant), labeled as _C_ ; 2) a linear dependence,
$\langle\Theta_{\text{jet}}(t)\rangle=a_{1}*(t-t_{\text{ave}})+a_{0}$, where
$t_{\text{ave}}$ is the midpoint of the observed period, labeled with _L_ ;
and 3) a cubic spline fit,
$\langle\Theta_{\text{jet}}(t)\rangle=\text{spline}(t-t_{\text{ave}})$,
labeled as _S_.

In order to perform a more quantitative placement of the sources into these
categories, we have calculated the average inner jet direction,
$\langle\Theta_{\text{jet}}\rangle$, and standard deviation,
$\sigma_{\langle\Theta_{\text{jet}}\rangle}$, by averaging $\Theta$ of all
moving knots within a distance from the core of $\sim 3\times$ the observing
beam FWHM (i.e., the inner jet). For several jets, there are one or two knots
with $\Theta$ significantly different from the others, which we exclude from
this calculation. We have compared
$\sigma_{\langle\Theta_{\text{jet}}\rangle}$ with the average uncertainty of
the knot position angle measurement within the inner jet,
$\langle\sigma_{\Theta}\rangle$. If
$\sigma_{\langle\Theta_{\text{jet}}\rangle}\leq
3\times\langle\sigma_{\Theta}\rangle$, the position angle is classified as
constant, _C_. For sources where this criterion is not met, we have employed
either linear or spline fits as determined from visual inspection of the
behavior of $\Theta$ for all knots. To confirm that these fits provide a
better estimate of the jet position angle, we have compared the standard
deviation of the fit to the standard deviation of a constant
$\langle\Theta_{\text{jet}}\rangle$. In the majority of cases, the fit reduced
the standard deviation of the knot position angles to within the
$3\times\langle\sigma_{\Theta}\rangle$ criterion. This criterion remained
unmet after the fit in only four sources (0716+714, 1226+023, 1253-055, and
1749+096), but even in these cases the fits improved the standard deviation,
and so were adopted. Finally, a handful of sources have a very wide
dispersion, but no clear trend in the position angles of knots. For these
sources, we have created a fourth category (_W_ for “wide”), where we have
adopted a constant value of $\langle\Theta_{\text{jet}}\rangle$, but with a
standard deviation $>3\times\langle\sigma_{\Theta}\rangle$.

Table 7 describes the fits for each source as follows: 1—source name;
2—constant average position angle of the knots in a jet,
$\langle\Theta_{\text{jet}}\rangle$ (in degrees), regardless of the best fit;
3—average error of knot position angles, $\langle\sigma_{\Theta}\rangle$ (in
degrees); 4—categorization of the jet position angle behavior, with categories
as described above888Due to their complexity, spline fits are available upon
request. Linear fits are provided below.; 5—midpoint in time of all the knot
observations, $t_{\text{ave}}$, for the fits; 6—standard deviation of the best
fit to knot position angles, $\sigma_{\langle\Theta_{\text{jet}}\rangle}$, in
degrees; and 7—knots that were used to determine the average position angles
and fits.

The majority of sources in our sample (23/38, 60.5%) have knots in the jet
with position angles that are relatively constant over time, i.e., have a
narrow spread. For these we take the value in column 2 of Table 7 to be
$\langle\Theta_{\text{jet}}\rangle$ in the analyses that follow. Six other
sources (15.8% of the sample) are well described by a constant jet direction,
but have very wide jets in terms of the spread of the position angles of
knots. Despite the latter property, we adopt the value listed in column 2 for
the subsequent analyses. Only two sources (the FSRQs 0420$-$014 and 1226+023)
have clear linear trends of knot position angle variations over time. The
linear fit for 0420$-$014 is $\langle\Theta_{\text{jet}}(t)\rangle=(5.9\pm
0.2)*(t-2013.24)+(127.0\pm 0.7)$ degrees, while the fit for 1226+023 is
$\langle\Theta_{\text{jet}}(t)\rangle=(1.2\pm 0.1)*(t-2013.79)-(140.9\pm 0.2)$
degrees. The remaining seven sources (18.4%) display very complex behavior.
Some (e.g., the FSRQ 1222+216) visually resemble jet rotation (possibly
precession), while others exhibit abrupt changes in the jet position angle
(e.g., the FSRQ 1253$-$055). For the sources best fit with a linear trend or
spline, we use the fits to calculate the jet position angle at the time of
acceleration for the analysis in $\S$4.3. We use the average of the fit over
the observed time period for each knot for the analysis of stationary features
in $\S$4.2. For all sources, we use
$\sigma_{\langle\Theta_{\text{jet}}\rangle}$ as the uncertainty on
$\langle\Theta_{\text{jet}}\rangle$ in the subsequent analyses.

Table 7: Jet Position Angles

Source | $\langle\Theta_{\text{jet}}\rangle$ | $\langle\sigma_{\Theta}\rangle$ | TypeaaFinal fit to the jet position angle behavior, labeled as follow: C—$\langle\Theta_{\text{jet}}(t)\rangle=$ constant, for which the value in column 2 is used in subsequent analyses in this work; W—sources for which $\langle\Theta_{\text{jet}}(t)\rangle=$ constant, but $\sigma_{\langle\Theta_{\text{jet}}\rangle}>3\times\langle\sigma_{\Theta}\rangle$; L—linear dependence, $\langle\Theta_{\text{jet}}(t)\rangle=a_{1}*(t-t_{\text{ave}})+a_{0}$, where $t_{\text{ave}}$ is the midpoint of the observed period; and S—a spline fit of the 3rd order, $\langle\Theta_{\text{jet}}(t)\rangle=\text{spline}(t-t_{\text{ave}})$. | $t_{\text{ave}}$ | $\sigma_{\langle\Theta_{\text{jet}}\rangle}$ | Knots Included
---|---|---|---|---|---|---
| [deg] | [deg] |  | [yr] | [deg] |
0219+428 | $-169.6$ | $\phantom{0}5.5$ | C | $\ldots$ | $9.0$ | B1—B6
0235+164 | $54.7$ | $\phantom{0}6.8$ | C | $\ldots$ | $6.6$ | B5, B6
0316+413 | $173.3$ | $\phantom{0}0.4$ | W | $\ldots$ | $9.6$ | All except C4, C8, C12, C13
0336–019 | $78.7$ | $\phantom{0}4.6$ | C | $\ldots$ | $6.4$ | B1—B5
0415+379 | $66.1$ | $\phantom{0}2.2$ | C | $\ldots$ | $4.9$ | K1—K31
0420–014 | $121.0$ | $\phantom{0}3.7$ | L | $2013.24$ | $8.3$ | B1—B9
0430+052 | $-113.5$ | $\phantom{0}2.9$ | C | $\ldots$ | $3.2$ | C1—C15
0528+134 | $64.8$ | $\phantom{0}6.0$ | C | $\ldots$ | $13.3$ | B1—B3; C1—C5
0716+714 | $24.3$ | $\phantom{0}3.9$ | S | $2012.98$ | $18.4$ | A1, A2; B1—B14; T10
0735+178 | $78.5$ | $10.7$ | C | $\ldots$ | $13.8$ | B1—B4; C1
0827+243 | $127.1$ | $\phantom{0}5.8$ | C | $\ldots$ | $8.1$ | B1—B9
0829+046 | $63.1$ | $\phantom{0}5.7$ | C | $\ldots$ | $10.6$ | A1—A4; B1—B6
0836+710 | $-140.0$ | $\phantom{0}3.3$ | C | $\ldots$ | $9.0$ | B1—B6
0851+202 | $-24.1$ | $\phantom{0}1.2$ | S | $2012.44$ | $2.7$ | a
0954+658 | $-13.7$ | $\phantom{0}2.1$ | S | $2013.97$ | $5.5$ | A1, A2
1055+018 | $-63.8$ | $\phantom{0}6.5$ | C | $\ldots$ | $18.4$ | B1; B4; T1
1101+384 | $-14.4$ | $14.2$ | C | $\ldots$ | $13.9$ | A1; B1, B3, B4
1127–145 | $80.7$ | $\phantom{0}4.3$ | C | $\ldots$ | $10.4$ | B1, B2; C1, C2; D1
1156+295 | $10.0$ | $\phantom{0}3.8$ | C | $\ldots$ | $11.1$ | B1—B7
1219+285 | $106.2$ | $10.6$ | C | $\ldots$ | $8.0$ | B1—B4
1222+216 | $23.5$ | $\phantom{0}1.9$ | S | $2014.11$ | $4.7$ | A1, A2, A3
1226+023 | $-141.3$ | $\phantom{0}0.7$ | L | $2013.79$ | $4.3$ | B1—B22
1253–055 | $-152.7$ | $\phantom{0}1.0$ | S | $2013.37$ | $20.0$ | A1; C24—C38
1308+326 | $-64.9$ | $\phantom{0}4.7$ | C | $\ldots$ | $12.9$ | B2—B6
1406–076 | $-81.9$ | $\phantom{0}7.2$ | C | $\ldots$ | $14.5$ | B1—B5; C1
1510–089 | $-33.1$ | $\phantom{0}2.8$ | C | $\ldots$ | $8.0$ | All except B10
1611+343 | $163.3$ | $\phantom{0}3.3$ | W | $\ldots$ | $13.0$ | A2, A3; C1—C3
1622–297 | $-78.1$ | $\phantom{0}3.3$ | W | $\ldots$ | $11.7$ | B4—B6
1633+382 | $-51.6$ | $\phantom{0}6.1$ | C | $\ldots$ | $9.1$ | B1, B2, B4, B5, B7
1641+399 | $-100.5$ | $\phantom{0}2.6$ | W | $\ldots$ | $8.6$ | B1—B15; C1—C3; D1—D3
1652+398 | $119.7$ | $11.8$ | C | $\ldots$ | $26.9$ | A2—A5; C2, C3
1730–130 | $-2.8$ | $\phantom{0}1.4$ | W | $\ldots$ | $13.6$ | A1; B1—B7; C1; D; T2
1749+096 | $-1.2$ | $\phantom{0}2.3$ | S | $2014.03$ | $15.0$ | A1; B1—B10; C1; T5
1959+650 | $127.9$ | $12.9$ | C | $\ldots$ | $19.1$ | A1—A4
2200+420 | $-170.0$ | $\phantom{0}1.8$ | C | $\ldots$ | $4.5$ | B1—B16; C1—C5; D1, D2; T3
2223–052 | $114.3$ | $\phantom{0}3.2$ | W | $\ldots$ | $14.7$ | A1—A3; B1—B6; C1, C2; T5
2230+114 | $130.9$ | $\phantom{0}3.6$ | C | $\ldots$ | $10.1$ | B1—B10; D1, D2
2251+158 | $-81.8$ | $\phantom{0}1.6$ | S | $2012.67$ | $2.5$ | C

Note. — PAs are measured north through east.

##  4 Motions of Knots in the Jets

Between 2007 and 2018 we have identified 559 distinct emission features in the
parsec-scale jets of the 38 blazars in the VLBA-BU-BLAZAR sample. Of this
total number, 38 are the 43 GHz cores. Among the remaining knots, 96 have been
classified as quasi-stationary (18.4% of non-core components), leaving 425
classified as moving knots. Accelerations along the jet cause 75 of the moving
knots (17.6%) to have multiple estimates for their apparent speed as they
traverse different regions along their trajectory. The only sources that do
not exhibit apparent superluminal motion in the parsec-scale jet are the BLs
1652+398 and 1959+650. Interestingly, we find superluminal motion in the
generally subluminal jet of the RG 0316+413 (3C 84) in one knot, _C10_ , with
an apparent speed of $\beta_{\text{app}}=1.38c\pm 0.03c$. We consider this
estimate of the knot speed to be reliable, with multiple observations of its
motion close to the core and a well-defined epoch of ejection. Also, we detect
knot motions up to $\sim 2c$ in the BL 1101+384 (Mkn 421).

For each source, we have constructed a plot displaying the separation of all
knots from the core (Figure SET 7). Each knot is color-coded in the
corresponding figure in Figure SET 6 according to its motion type. In each
figure, the black vectors show the PA of knots with respect to the core at the
corresponding epoch. The solid lines represent the piece-wise linear fits to
the knot motions, the red dotted line shows the position of the core, and
dashed colored lines mark the average positions of the stationary components.
Error bars on each measurement are the approximate $1\sigma$ positional
uncertainties based on $T_{\text{b,obs}}$. Lines extrapolating the knot motion
back to the epoch of ejection $t_{\circ}$ may appear slightly bent in some
cases owing to (1) $R$ being a quadratic combination of $X$ and $Y$, or (2)
instances where $|t_{\text{x}\circ}-t_{\text{y}\circ}|\geq 0.5$ yr so that the
epoch of ejection is calculated using a linear fit to the first segment of _R_
data. These apparent bends have no effect on the apparent speeds. Figure 7
presents an example of these plots for the BL 0219+428 (3C 66A).

Fig. Set 7. Knot Separations from Core

Figure 7: Separation vs. time of the knots in the jet of the BL 0219+428 (3C
66A). The colors of the knots — the same as in Figure 6 — are coded based on
their motion designation: blue, green, and cyan represent bright moving knots,
orange and purple correspond to stationary components, and brown is used for
trailing knots (appearing to split off from other knots), relatively faint
moving features, or to differentiate between closely spaced knots. Red dotted
lines are used to show the core component _A0_. Horizontal lines along _A0_
indicate the uncertainty in time of the epoch of ejection. See the text for
more details. Plots for all sources in the sample are available in the figure
set.
(The complete figure set (38 images) is available online.) Table 8: Velocity
in the Jets

Source | Knot | $N_{\text{seg}}$ |  $t_{\text{s}}$aaKnots with a single motion segment have no error in their start and end times as these are the first and last observation of the knot, respectively. The typical error on the segment start and end times for knots fit with multiple segments is $\sigma_{T}\lesssim 0.01$ yr and are determined from the covariance matrices of fit parameters in the least-squares piece-wise linear fitting. |  $t_{\text{e}}$aaKnots with a single motion segment have no error in their start and end times as these are the first and last observation of the knot, respectively. The typical error on the segment start and end times for knots fit with multiple segments is $\sigma_{T}\lesssim 0.01$ yr and are determined from the covariance matrices of fit parameters in the least-squares piece-wise linear fitting. | $\mu$ | $\Phi$ | $\beta_{\text{app}}$ | $t_{\circ}$ | Flag
---|---|---|---|---|---|---|---|---|---
|  |  | [yr] | [yr] | [mas yr-1] | [deg] | [$c$] | [yr] | for $t_{\circ}$
0219+428 | B1 | 1 | $2009.56$ | $2010.11$ | $0.915\pm 0.178$ | $\phantom{-}179.2\pm\phantom{0}4.0$ | $24.62\pm 4.80$ | $2008.67\pm 0.24$ | R
| B2 | 1 | $2010.84$ | $2012.07$ | $0.483\pm 0.058$ | $-167.0\pm\phantom{0}4.9$ | $13.01\pm 1.56$ | $2009.24\pm 0.27$ | R
| B3 | 1 | $2012.07$ | $2013.29$ | $0.665\pm 0.070$ | $-168.5\pm 11.7$ | $17.89\pm 1.87$ | $2011.93\pm 0.01$ | XY
| B4 | 1 | $2014.15$ | $2015.93$ | $0.751\pm 0.056$ | $-170.2\pm\phantom{0}1.4$ | $20.22\pm 1.50$ | $2013.25\pm 0.08$ | XY
| B5 | 1 | $2013.96$ | $2016.31$ | $0.634\pm 0.050$ | $-168.4\pm\phantom{0}2.1$ | $17.07\pm 1.35$ | $2013.80\pm 0.03$ | XY
| B6 | 1 | $2017.85$ | $2018.94$ | $0.605\pm 0.067$ | $-165.6\pm\phantom{0}5.7$ | $16.28\pm 1.79$ | $2016.12\pm 0.28$ | R
0235+164 | B1 | 1 | $2007.66$ | $2008.62$ | $0.533\pm 0.025$ | $-\phantom{0}18.3\pm\phantom{0}3.5$ | $26.58\pm 1.27$ | $2007.45\pm 0.01$ | XY
| B2 | 1 | $2008.62$ | $2009.62$ | $0.258\pm 0.024$ | $\phantom{-}140.3\pm\phantom{0}4.9$ | $12.88\pm 1.18$ | $2008.15\pm 0.11$ | R
| B3 | 1 | $2010.58$ | $2013.96$ | $0.051\pm 0.004$ | $\phantom{-}147.1\pm\phantom{0}4.8$ | $\phantom{0}2.55\pm 0.21$ | $2008.10\pm 0.32$ | R
| B4 | 1 | $2013.49$ | $2016.31$ | $0.107\pm 0.014$ | $-124.2\pm\phantom{0}8.9$ | $\phantom{0}5.32\pm 0.68$ | $2012.62\pm 0.26$ | R
| B5 | 1 | $2016.08$ | $2017.67$ | $0.183\pm 0.017$ | $\phantom{-}\phantom{0}66.7\pm\phantom{0}4.8$ | $\phantom{0}9.13\pm 0.87$ | $2015.54\pm 0.15$ | XY
| B6 | 1 | $2016.91$ | $2018.65$ | $0.120\pm 0.013$ | $\phantom{-}\phantom{0}13.4\pm\phantom{0}8.0$ | $\phantom{0}6.00\pm 0.63$ | $2015.42\pm 0.51$ | R
| B7 | 1 | $2018.30$ | $2018.79$ | $0.191\pm 0.046$ | $\phantom{-}\phantom{0}52.3\pm 16.3$ | $\phantom{0}9.55\pm 2.28$ | $2016.55\pm 0.95$ | R
0336–019 | B1 | 1 | $2007.75$ | $2009.15$ | $0.617\pm 0.032$ | $\phantom{-}\phantom{0}81.4\pm\phantom{0}1.4$ | $28.57\pm 1.50$ | $2007.30\pm 0.02$ | XY
| B2 | 1 | $2008.81$ | $2009.50$ | $0.227\pm 0.003$ | $\phantom{-}\phantom{0}81.2\pm\phantom{0}4.3$ | $10.50\pm 0.13$ | $2008.37\pm 0.16$ | R
| B2 | 2 | $2009.50$ | $2010.45$ | $0.650\pm 0.001$ | $\phantom{-}\phantom{0}86.9\pm\phantom{0}1.5$ | $30.11\pm 0.05$ | $\ldots$ | $\ldots$

Note. — (Table 8 (529 lines) is published in its entirety in the machine-
readable format. A portion is shown here for guidance regarding its form and
content.)

###  4.1 Moving Feature Properties

Here we describe the properties of the moving features in the jets of sources
in each of the three subclasses of blazars. Because knot _C_ of the FSRQ
2251+158 (3C 454.3) shows motion despite being a quasi-stationary feature (see
$\S$4.2), we include it in the analysis of 426 moving knots. Since any
individual knot can have multiple motions over different time periods, we
analyze a total of 529 knot speeds.

Table 8 presents the speeds of the moving knots as follows: 1—name of the
source; 2—designation of the component; 3—segment number of the fit to the
knot motion, $N_{\text{seg}}$ (1 for purely linear motion in both the $X$ and
$Y$ direction, up to 5999If best fit in a given dimension is with a 3-segment
broken linear fit, then there are 2 breakpoints in that direction of motion.
If both the _X_ and _Y_ motions are fit with a 3-segment broken linear fit,
there are 4 total breakpoints and 5 line segments.); 4—start date of the
segment, $t_{\text{s}}$, in yr; 5—end date of the segment, $t_{\text{e}}$, in
yr; 6—proper motion, $\mu$, and its associated $1\sigma$ uncertainty, in mas
yr-1; 7—direction of motion, $\Phi$, and its associated $1\sigma$ uncertainty,
in degrees; 8—apparent speed, $\beta_{\text{app}}$, and its associated
$1\sigma$ uncertainty, in units of $c$; 9—epoch of ejection, $t_{\circ}$, and
its associated $1\sigma$ uncertainty (only available for the first segment of
a knot)in yr; and 10—flag indicating whether $t_{\circ}$ was calculated to the
combined $X$ and $Y$ dimension fit (XY), or if a fit to the $R$ coordinate of
the first segment was necessary (R, see $\S$2.2). If the segment was not the
first segment for a particular knot, the flag is NA.

Figure 8 (left) shows the distributions of the apparent speeds of the features
for FSRQs, BLs, and RGs separately for the entire observing period from 2007
June to 2018 December for ease of comparison with previous literature. An
increase in the number of moving knots by a factor of $\sim 3$ has
dramatically filled in the FSRQ and BL distributions over those presented in
Jorstad et al. (2017). The main distribution of apparent speeds of features in
the jets of FSRQs covers a wide range of $\beta_{\text{app}}$, from $0.5c$ to
$\sim 50c$. Two knots (_B1_ of 0528+134 and _D2_ of 2230+114) have segments
where the speed of the knot exceeds $50c$. In order to limit the size of the
figure, these have been included in the last bin. They represent a small
percentage (0.6%) of the total number of knots and thus do not affect the
conclusions. The FSRQ distribution peaks in the $8$-$10c$ bin but has a very
long tail, extending continuously with at least one knot in each bin up to the
$38$-$40c$ range. The apparent speeds in the jets of BLs exhibit a very
different distribution, with a peak in the $0$-$2c$ bin and a similarly-sized
peak in the $8$-$10c$ bin. The distribution does not extend to
$\beta_{\text{app}}$ as high as for FSRQs, with the maximum BL speed of $\sim
27c$ occurring for a segment of the knot B3 of 1749+096. The distribution of
apparent speeds of RG knots is limited to much lower values than that of the
other two subclasses. The peak between $4c$ and $6c$ corresponds to 0415+379
and 0430+052, with the peak from $0c$ to $2c$ coming primarily from 0316+413
(3C 84).

Figure 8: Distributions of the apparent speeds of all moving knots (left), and
weighted average of the apparent speed for each source (right) detected in
FSRQs (top), BLs (middle), and RGs (bottom).

In order to more appropriately determine the differences in the distributions
of apparent speed for FSRQs, BLs, and RGs, we calculate a single average
apparent speed, $\langle\beta_{\mathrm{app}}\rangle$, for each source as the
weighted average of the apparent speeds of all knots in that source. The
distributions of these average apparent speeds are shown in Figure 8 (right).
We see that the FSRQs have a well-determined peak in the distribution with
$\langle\beta_{\mathrm{app}}\rangle\sim 12$-$14c$. The BLs have a very wide
dispersion of average apparent speeds. With just three sources, it is easy to
see (with Table 8) that the knots of 0316+413 are, on average, subluminal,
while those of 0415+379 and 0430+052 are superluminal. Despite the small
numbers of sources (especially for RGs), we have performed a KS test between
the distributions. The BL and RG distributions are the most similar, with
$\mathcal{D}=0.500$ and $p=0.525$, likely due to the small number of sources.
However, moderately-significant differences are seen between the FSRQ:BL
($\mathcal{D}=0.492,\ p=0.029$) and FSRQ:RG ($\mathcal{D}=0.955,\ p=0.003$)
distributions. Taking all of the results of Figure 8 in consideration, we see
that the knots of FSRQs have, on average, higher apparent speeds than those of
BLs, with RGs having the lowest apparent speeds.

We analyze the direction of the velocity vector of each knot by following the
same approach as in Jorstad et al. (2017), calculating the difference between
the direction $\Phi$ of the apparent velocity (given in Table 8) and the
average position angle of the jet, $\langle\Theta_{\text{jet}}\rangle$ (listed
in Table 3). Instead of determining an average value of this parameter for
each source, we consider the difference between $\Phi$ and
$\langle\Theta_{\text{jet}}\rangle$ for a given knot to be an independent
measure of the possible dispersion of the position angles for each subclass of
blazar, as all knots are required to determine the width/dispersion of a jet.
Figure 9 plots the distributions of $|\Phi-\langle\Theta_{\text{jet}}\rangle|$
for FSRQs, BLs, and RGs. The peak of each distribution lies within the first
bin ($<10\arcdeg$) and represents 34.1%, 40.3%, and 74.7% of the total number
of knots in FSRQs, BLs, and RGs respectively.

In order to measure the deviation from unidirectional motion, we first
calculate the average uncertainty in the motion direction as
$\langle\sigma(\Phi)\rangle=\sum_{i=1}^{N_{k}}\sigma(\Phi_{i})/N_{k}$, where
$N_{k}$ is the total number of knot speeds measured in each subclass. BLs have
the largest dispersion of knot motions, with $\langle\sigma(\Phi)\rangle=7\fdg
2$, while FSRQs and RGs have a dispersion $\langle\sigma(\Phi)\rangle=4\fdg 2$
and $4\fdg 8$, respectively. The much larger dispersion in this work compared
to those presented in Jorstad et al. (2017) is likely caused by a
significantly larger number of measured speeds based on the revised method for
calculating motion. In order to determine whether this dispersion is from
varying directions of the velocity vectors of knots or from knots smaller than
the jet cross-section being displaced from the jet axis by various distances,
we calculate the average position-angle dispersion of knots in a subclass as
$\langle\sigma(\langle\Theta\rangle)\rangle=\sum_{i=1}^{N_{k}}\sigma(\langle\Theta\rangle_{i})/N_{k}$.
Based on this, FSRQs have the widest range of average knot position angles,
with $\langle\sigma(\langle\Theta\rangle)\rangle=9\fdg 9$, while BLs and RGs
have smaller dispersions of $8\fdg 8$ and $5\fdg 7$, respectively. We conclude
that, while the knots of FSRQs are more inclined to move along the same
direction in the jet than are the knots in the jets of BLs, FSRQs have wider
projected opening angles. The latter can be connected with a smaller viewing
angle and/or a wider intrinsic opening angle for FSRQs. FSRQs and BLs possess
a similar percentage of knots in which the velocity vector deviates from the
jet axis by a value between $30\arcdeg$ and $60\arcdeg$, 12.5% and 14.4%
respectively. In contrast, only $2.5\%$ of the knots of RGs have motions that
fall within this deviation range.

Among superluminal knots there are 5, 13, and 19 segments in FSRQs, BLs, and
RGs, respectively, with very slow ($<1c$) or even upstream motion
corresponding to $|\Phi-\langle\Theta_{\text{jet}}\rangle|>90\arcdeg$. These
are shown in the last bin of the distribution, since they represent a small
fraction of the total number of knot segments. A KS test indicates that the
FSRQ and BL distributions are _not_ significantly different
($\mathcal{D}=0.138,\ p=0.046$). However, the FSRQ and BL distributions _are_
statistically different from the RG distribution ($\mathcal{D}=0.439,\
p=1.8\times 10^{-11}$ for FSRQs vs RGs, and $\mathcal{D}=0.369,\ p=1.2\times
10^{-6}$ for BLs vs RGs).

Figure 9: Distributions of the differences between the velocity vector,
$\Phi$, and the jet axis, $\langle\Theta_{\text{jet}}\rangle$, in FSRQs (top),
BLs (middle), and RGs
(bottom).$|\Phi-\langle\Theta_{\text{jet}}\rangle|>90\arcdeg$ indicates slow
or upstream motion of the knot segment.

###  4.2 Quasi-Stationary Features

Figure 10: _Left:_ distribution of the projected distances from the core to
the stationary features in our sample. _Right:_ trajectories of the stationary
features _A1_ in 0528+134 (FSRQ, red circles), _A2_ in 1749+096 (BL, blue
squares), and _A1_ in 0430+052 (RG, green triangles) with respect to the core,
_A0_ (the black cross at position (0,0)). The solid lines show $2\sigma$
confidence ellipses used to characterize the dispersion and orientation of the
knot motion

We have identified 96 features with no statistically significant motion (as
defined in Jorstad et al., 2017). Of these, 45 are in FSRQs, 38 in BLs, and 13
in RGs. All sources in our sample except 0235+164 (BL), 0316+413 (RG),
1156+295, and 1622$-$297 (FSRQs) continue to contain at least one — but
usually multiple — stationary feature in addition to the core. On average,
FSRQs possess $\sim 2$ and BLs $\sim 3$ such features. Our sample size of RGs
is small, however, we see 7 and 5 stationary features in 0415+379 and
0430+052, respectively. Given the sub-relativistic speeds of knots in
0316+413, it is hard to determine what features are truly stationary. Almost
half of the stationary features in our sample have appeared or disappeared
over time, indicating that while these features may be stationary in position,
they are often transient features over timescales of $\sim 10$ years. Given
that there is no way to average the properties of stationary features for a
particular source into one overarching “typical” stationary feature, we treat
each stationary feature as an independent feature when making comparisons
between the subclasses of blazar.

Table 9: Quasi-Stationary Jet Features

Source | Knot | $N$ | $\langle R\rangle$ | $X_{\text{cen}}$ | $Y_{\text{cen}}$ | $A$ | $B$ | $\varpi$
---|---|---|---|---|---|---|---|---
|  |  | [pc] | [mas] | [mas] | [mas] | [mas] | [deg]
0219+428 | A1 | 52 | $\phantom{0}0.97\pm 0.23$ | $-0.04$ | $-0.16$ | $0.08$ | $0.06$ | $21.9$
| A2 | 36 | $\phantom{0}2.11\pm 0.34$ | $-0.08$ | $-0.35$ | $0.15$ | $0.10$ | $-127.7$
| A3 | 50 | $\phantom{0}3.71\pm 0.51$ | $-0.07$ | $-0.64$ | $0.19$ | $0.11$ | $-21.7$
| A4 | 56 | $13.20\pm 0.51$ | $-0.12$ | $-2.31$ | $0.18$ | $0.11$ | $-3.3$
0336–019 | A1 | 36 | $\phantom{0}1.23\pm 0.31$ | $\phantom{-}0.14$ | $\phantom{-}0.02$ | $0.10$ | $0.08$ | $13.9$
| A2 | 7 | $14.64\pm 0.84$ | $\phantom{-}1.77$ | $\phantom{-}0.69$ | $0.28$ | $0.21$ | $16.0$
| A3 | 8 | $\phantom{0}2.15\pm 0.23$ | $\phantom{-}0.27$ | $\phantom{-}0.07$ | $0.11$ | $0.04$ | $-46.4$
| A4 | 11 | $\phantom{0}3.14\pm 0.46$ | $\phantom{-}0.39$ | $\phantom{-}0.11$ | $0.17$ | $0.09$ | $-41.8$

Note. — (Table 9 (96 lines) is published in its entirety in the machine-
readable format. A portion is shown here for guidance regarding its form and
content.)

Figure 10, left, shows the distribution of the average projected linear
distances from the core of all stationary features in our sample. Five
stationary components that appear beyond 15 pc from the core are included in
the last bin. The distribution has a prominent peak at distances $<1$ pc from
the core, as in the distribution presented in Jorstad et al. (2017). In
addition, there continues to be no statistical difference in the distributions
of FSRQs, BLs, and RGs, with a KS test yielding $\mathcal{D}=0.159$ and
$p=0.601$ for the distributions of FSRQs vs BLs, $\mathcal{D}=0.260$ and
$p=0.423$ for FSRQs vs RGs, and finally $\mathcal{D}=0.221$ and $p=0.642$ for
BLs vs RGs.

Despite the “stationary feature” classification, the brightness centroids of
many of these components do exhibit distinguishable motion. The locations of
the majority of features tend to fluctuate about their mean reported
positions. Figure 10, right, shows three examples of the trajectories of
stationary features, one in each subclass. The range of motion in both right
ascension and declination is larger than the positional errors can account
for, and the loci of $(X,Y)$ positions of a given stationary feature form a
pattern that indicates a preferred direction of motion for the fluctuations.
In order to characterize the shifts in position, we fit confidence ellipses to
the motions of each stationary feature. Table 9 gives the parameters of the
elliptical fits as follows: 1—source name; 2—knot designation; 3—number of
observations, $N$, of a knot; 4—projected linear distance from the core to
$\langle R\rangle$ of the knot, in parsecs; 5—mean position of the knot along
the $X$-axis (relative R.A.), $X_{\text{cen}}$, in mas; 6—mean position of the
knot along the $Y$-axis (relative DEC), $Y_{\text{cen}}$, in mas; 7—semimajor
axis length, $A$, in mas; 8—semiminor axis length, $B$, in mas; and
9—orientation angle of semimajor axis relative to north, $\varpi$, in degrees.
Figure 10, right shows sample $2\sigma$ confidence ellipses for the given
knots, where $2\sigma$ indicates that 95% of the observed data would lie
within the ellipse if the data were normally distributed.

Figure 11: Distributions of the differences between the major axis orientation
angles, $\varpi$, and $\langle\Theta_{\text{jet}}\rangle$, in FSRQs (top), BLs
(middle), and RGs (bottom). See the text for details.

In this characterization of the knot motion dispersion, $\varpi$ represents
the angle of the motion along the major axis. We compare this direction of
motion with the direction of the jet through
$|\varpi-\langle\Theta_{\text{jet}}\rangle|$, as defined in Table 7. For
sources with $\langle\Theta_{\text{jet}}\rangle$ fit by a linear trend or a
spline, we average the trend over the dates of observation of the particular
stationary feature. Figure 11, left, shows the distributions of this
difference angle separately for FSRQs, BLs, and RGs. A KS test between the
distributions indicates that there is no statistical difference between the
subclasses, with the FSRQ and BL distributions being the most similar
($\mathcal{D}=0.162$, $p=0.578$), while $\mathcal{D}=0.318,$ $p=0.200$ for
FSRQs versus RGs and $\mathcal{D}=0.277$, $p=0.357$ for BLs versus RGs.

Figure 12: Distributions of the relative flux densities of quasi-stationary
features for the different subclasses—FSRQs (top), BLs (middle), and RGs
(bottom). Figure 13: Distributions of the average flux of a knot relative to
the average flux of the core, $\langle S_{\text{core}}\rangle$, for the
stationary (left) and moving (right) features in FSRQs (top), BLs (middle),
and RGs (bottom).

We are able to separate the entire range of
$|\varpi-\langle\Theta_{\text{jet}}\rangle|$ into three categories: (1)
fluctuations along the jet, with
$|\varpi-\langle\Theta_{\text{jet}}\rangle|\leq 30\arcdeg$ or
$|\varpi-\langle\Theta_{\text{jet}}\rangle|>150\arcdeg$, (2) fluctuations
transverse to the jet direction,
$60\arcdeg<|\varpi-\langle\Theta_{\text{jet}}\rangle|\leq 120\arcdeg$, and (3)
oblique fluctuations with intermediate angle differences,
$30\arcdeg<|\varpi-\langle\Theta_{\text{jet}}\rangle|\leq 60\arcdeg$ or
$120\arcdeg<|\varpi-\langle\Theta_{\text{jet}}\rangle|\leq 150\arcdeg$. Figure
11 shows the divisions between parallel and oblique cases with red dotted
lines, while the division between oblique and transverse fluctuations is shown
with blue dotted lines. In FSRQs, 40% (18/45) of knots fluctuate along the jet
direction, 26.7% (12/45) oblique to the jet, and 33.3% (15/45) transverse to
the jet. The situation is similar for BLs: 44.7% (17/38), 26.3% (10/38), and
28.9% (11/38), respectively. This is supported by the KS test, which indicates
that the FSRQ and BL distributions are not significantly different. The RG
distribution, however, shows a clear preference for motions along the jet
direction, with 11 of the 13 knots (92%) moving predominately along the jet
direction. While the number of RG stationary features has increased over those
detected in Jorstad et al. (2017), there are still only a relatively small
number of such knots, so that the KS test is inconclusive.

Table 10: Acceleration in the Jets

Source | Knot | $N_{\text{acc}}$ |  $t_{\mathrm{start}}$aaWe take the typical uncertainties in the estimated start and end dates of the acceleration region to be similar to the uncertainty in the break points for the piece-wise linear fitting, which are $\sigma_{t}\lesssim 0.01$ yr. The uncertainties in the acceleration are dominated by the uncertainty in the speeds rather than the time period of acceleration. |  $t_{\mathrm{end}}$aaWe take the typical uncertainties in the estimated start and end dates of the acceleration region to be similar to the uncertainty in the break points for the piece-wise linear fitting, which are $\sigma_{t}\lesssim 0.01$ yr. The uncertainties in the acceleration are dominated by the uncertainty in the speeds rather than the time period of acceleration. | $\dot{\mu}_{\parallel}$ | $\dot{\mu}_{\perp}$ | $R_{\dot{\mu}}$ | $\eta_{\parallel}$ | $\eta_{\perp}$
---|---|---|---|---|---|---|---|---|---
|  |  | [yr] | [yr] | [mas yr-2] | [mas yr-2] | [mas] | [yr-1] | [yr-1]
0316+413 | C1a | 1 | $2014.69$ | $2017.58$ | $\phantom{-}0.127\pm 0.003$ | $\phantom{-}0.015\pm 0.023$ | 2.86 | $\phantom{-}0.538\pm 0.016$ | $\phantom{-}0.064\pm 0.098$
| C10 | 1 | $2016.32$ | $2016.89$ | $-1.778\pm 0.073$ | $-0.209\pm 0.307$ | 1.19 | $-1.558\pm 0.069$ | $-0.183\pm 0.269$
| C10 | 2 | $2016.89$ | $2018.04$ | $\phantom{-}0.404\pm 0.037$ | $\phantom{-}0.047\pm 0.077$ | 1.24 | $\phantom{-}1.827\pm 0.289$ | $\phantom{-}0.213\pm 0.349$
0336–019 | B2 | 1 | $2009.16$ | $2009.98$ | $\phantom{-}0.508\pm 0.013$ | $\phantom{-}0.102\pm 0.064$ | 0.25 | $\phantom{-}4.145\pm 0.119$ | $\phantom{-}0.832\pm 0.522$
| B3 | 1 | $2011.63$ | $2013.88$ | $\phantom{-}0.075\pm 0.002$ | $\phantom{-}0.015\pm 0.009$ | 0.28 | $\phantom{-}1.094\pm 0.039$ | $\phantom{-}0.219\pm 0.131$
| B4 | 1 | $2015.33$ | $2016.94$ | $\phantom{-}0.222\pm 0.005$ | $\phantom{-}0.044\pm 0.026$ | 0.46 | $\phantom{-}2.056\pm 0.051$ | $\phantom{-}0.407\pm 0.241$

Note. — (Table 10 (104 lines) is published in its entirety in the machine-
readable format. A portion is shown here for guidance regarding its form and
content.)

Almost all of the apparent motions of the quasi-stationary features appear to
be minor fluctuations about the average positions. A clear trend in the motion
is apparent only for knot _C_ of the FSRQ 2251+158. Over $\sim 10$ years of
observations, this feature started $\sim 0.7$ mas from the 43 GHz core in
2007, after which the separation decreased to a minimum of $\sim 0.45$ mas in
2013. Since then, the separation has been increasing more gradually. It is
unclear whether this motion is related to changes in the intensity
distribution of this extended feature, or to bright knots that were ejected
from the core and later merged with _C_. An alternative possibility is that
the change in separation between the core and _C_ in 2251+158 is due to the
motion of the core itself, contrary to the assumption made in this work that
every core is at an essentially stationary position in the jet.

Figure 12 shows distributions of the average flux density of the stationary
features, $\langle S_{\text{knot}}\rangle$, normalized by the average flux
density of the core, $\langle S_{\text{core}}\rangle$ for all three
subclasses. A KS test confirms the distributions for each subclass are
similar, with $\mathcal{D}=0.201,\ p=0.319$ for FSRQs vs. BLs,
$\mathcal{D}=0.244,\ p=0.501$ for FSRQs vs. RGs, and $\mathcal{D}=0.223,\
p=0.628$ for BLs vs. RGs.

Finally, we investigate whether the stationary and moving features in the jets
have distributions of knot:core flux density ratios that are statistically
different within a subclass. We calculate the average $\langle
S_{\text{knot}}\rangle$, median, and maximum flux density of each knot and
normalize them by the average flux density of the core, $\langle
S_{\text{core}}\rangle$. Figure 13 displays the distributions of flux ratios
for stationary features (left) and moving features (right) separately for each
subclass. There is no statistical difference between the flux ratios of the
stationary and moving features in RGs according to a K-S test
($\mathcal{D}=0.333,\ p=0.133$). The hypothesis of different distributions in
FSRQs is similarly rejected ($\mathcal{D}=0.191,\ p=0.110$). In BLs, however,
there is a moderately-significant result that the stationary features of BLs
have higher flux ratios than do the moving features ($\mathcal{D}=0.275$,
$p=0.019$). The flux distributions and KS statistics for the average and
maximum flux ratios of stationary and moving components show a similar trend,
but are not included here for brevity.

###  4.3 Knot Acceleration and Deceleration

Between 2007 and 2018, 75 of the 425 (17.6%) moving knots detected in our
sample exhibit non-ballistic apparent velocities, which we define as motions
that are best fit by multiple line segments in a piece-wise fashion. We also
include the quasi-stationary feature _C_ of the FSRQ 2251+158 in the following
analysis, as it has been reliably observed to have motion around the average
stationary value. Across the three subclasses, we have detected acceleration
in 55 features in 16 FSRQs, 9 features in 6 BLs, and 11 features in the 3 RGs.
Due to the nature of our piece-wise linear fits, the trajectories of many of
these knots contain multiple acceleration regions. This totals 104 individual
acceleration regions in the jets of 25 blazars in our sample. Given that the
multiple acceleration regions may have different properties, we consider each
region to be an independent measurement of possible acceleration within a
particular subclass of blazar, and do not calculate a “typical” acceleration
for each source. Table 10 gives the accelerations of each region in each knot
as follows: 1—name of the source; 2—designation of the component;
3—acceleration region number, $N_{\text{acc}}$; starting with 1 for the first
region present (between moving segments 1 and 2) in a knot trajectory and
reaching a maximum of 4 for a knot which is best fit by three linear segments
in both the $X$ and $Y$ direction; 4—beginning date of the acceleration,
$t_{\mathrm{start}}$; 5—end date of the acceleration, $t_{\mathrm{end}}$;
6—the acceleration parallel the jet, $\dot{\mu}_{\parallel}$, and its
$1\sigma$ uncertainty in mas yr-2; 7—the acceleration perpendicular to the
jet, $\dot{\mu}_{\perp}$, and its $1\sigma$ uncertainty in mas yr-2;
8—$R_{\dot{\mu}}$, the distance from the core to the position of the knot
based on the estimated time break point (from the piece-wise fit), in mas;
9—the normalized acceleration parallel to the jet, $\eta_{\parallel}$, and its
$1\sigma$ uncertainty in yr-1; and 10—the normalized acceleration
perpendicular to the jet, $\eta_{\perp}$, and its $1\sigma$ uncertainty in
yr-1 (see below for definitions of $\eta_{\parallel}$ and $\eta_{\perp}$). As
discussed in $\S$3.2, the accelerations parallel and perpendicular to the jet
were determined using the averaged position angle for sources with constant or
wide jets and using the linear or spline trend to calculate the jet position
angle at the time of acceleration for those with time variable jet positions.
In this way we separate the time variability of the jet as a whole from the
change in directions of knots relative to the bulk flow.

Figure 14: Left: relative accelerations parallel (top) and perpendicular
(bottom) to the jet of moving knots in FSRQs (red circles), BLs (blue
squares), and RGs (green triangles) vs. their angular distances from the core
of the jet. Right: the same accelerations vs. the average projected linear
distance downstream of the core.

The accelerations found in the jets of RGs are primarily in 0415+379 (8/11,
72.7%). Several knots in this source show very complex motions, requiring more
than two acceleration regions. Accelerations occur in only one knot in
0430+052. In contrast, only a few BLs exhibit any acceleration, and only one
BL (2200+420) contains a knot (_C2_) that shows complex enough motion to have
more than two acceleration zones. A much wider variety of knots are found to
have accelerated in FSRQs, but again a few sources dominate in terms of the
number of regions (such as 1226+023, with 21 regions for 13 knots).

Due to the fact that the 3D motion is projected onto the two-dimensional plane
of the sky, we cannot distinguish whether the acceleration in a given knot is
connected with an intrinsic change in speed (and thus in the Lorentz factor
$\Gamma$) or with a change in the intrinsic viewing angle $\Theta_{\circ}$. To
resolve this ambiguity, we follow the statistical approach suggested by Homan
et al. (2009) and followed by Jorstad et al. (2017): If, in a flux-limited
sample of beamed jets, the observed relative parallel acceleration exceeds
$\sim 60\%$ of the observed relative perpendicular acceleration (averaged over
the sample), then the observed accelerations are caused more by variations in
the Lorentz factors of the knots than by changes in the jet direction.

We compute the relative accelerations in our jets following a formalism
similar to that of Homan et al. (2009), adapted for our fitting procedure: the
acceleration parallel to the jet is
$\eta_{\parallel}=(1+z)\dot{\mu}_{\parallel}/\mu_{i}$, while acceleration
perpendicular to the jet direction is
$\eta_{\perp}=(1+z)\dot{\mu}_{\perp}/\mu_{i}$, where $\mu_{i}$ in both cases
refers to the pre-acceleration proper motion of the knot. Averaging over the
entire sample, we find $\langle\eta_{\parallel}\rangle=0.678$ yr-1 and
$\langle\eta_{\perp}\rangle=-0.056$ yr-1, thus,
$|\langle\eta_{\parallel}\rangle|>0.6|\langle\eta_{\perp}\rangle|$ and the
observed accelerations are likely caused by changes in the Lorentz factors
rather than changes in the jet direction.

One advantage of the piece-wise linear fits used in this work (instead of the
polynomial fits presented in Jorstad et al., 2017) is that, since the
accelerations can be seen to take place at specific locations in the jet, it
is straightforward to determine the distance from the 43 GHz core to the
acceleration zone. Figure 14, left, shows the computed relative accelerations
with respect to the angular distance from the core to the acceleration region,
while the right panel shows the same relative accelerations as functions of
the average projected linear distance from the core to the acceleration
region. While the accelerations in RGs tend to appear to occur closer to the
core than the accelerations in FSRQs and BLs, this is due to the observational
bias that the RG redshifts are smaller than those of the other subclasses.
Given that there are no distinct differences in relative acceleration between
the subclasses, we consider all accelerations together.

We bin both $\eta_{\parallel}$ and $\eta_{\perp}$ for all sources by distance
from the 43 GHz core of the jet, and construct histograms of the accelerations
in Figure 15. Distance bins are chosen so that a significant number of knot
accelerations would be present in each bin ($\gtrsim 10$) while also
maintaining adequate distance spacing. The left panels display
$\eta_{\parallel}$, while those on the right show $\eta_{\perp}$. In the
distributions, the red vertical lines indicate the median values, which are
listed in Table 11.

Figure 15: Histograms of relative knot accelerations (left) parallel and (right) perpendicular to the jet axis for all sources, separated by distance from the 43 GHz core. Red vertical lines in all panels indicate the median acceleration, while the gray dotted vertical lines indicate no acceleration.  Table 11: Median Values of Normalized Accelerations by Distance from the Core. Distance Bin | $\eta_{\parallel}$ | $\eta_{\perp}$
---|---|---
[pc] | [yr-1] | [yr-1]
0–1 | $\phantom{-}1.20$ | $\phantom{-}0.53$
1–2 | $\phantom{-}0.55$ | $-0.74$
2–3 | $\phantom{-}0.65$ | $\phantom{-}0.26$
3–5 | $\phantom{-}0.12$ | $-0.15$
5–10 | $-0.43$ | $-0.05$
10–35 | $-0.12$ | $-0.12$

From inspection of the distributions of $\eta_{\parallel}$, it is apparent
that accelerations close to the core are positive, increasing the speed at
which knots flow down the jet axis. As the distance from the core increases,
the magnitude of this acceleration generally decreases. Once knots reach
distances of $\sim 5$-10 pc from the core, most accelerations are negative.
While this work at 43 GHz probes projected distances of $\sim 0$-$35$ pc from
the core, data at 15 GHz (Homan et al., 2015) probes farther out, up to $\sim
100$ pc projected distance from the core. The data at longer wavelengths
supports our results here, with a general turnover from positive to
negative/no acceleration beyond $\sim 10$-20 pc projected distance. This
result indicates that, although the overall diversity in apparent speeds
likely requires a variety of different shock strengths if the knots are shock
waves (e.g., Marscher & Gear, 1985), it is common for accelerations to occur
close to the core and for constant speeds or decelerations to dominate farther
downstream ($\sim 10^{2}$ pc deprojected distance).

The distributions of $\eta_{\perp}$ indicate a different behavior of the knot
motions perpendicular to the jet axis. Beyond $\sim 2$ pc from the core, the
acceleration distributions are consistent with low, if any, acceleration
transverse to the jet axis. However, initially there appears to be positive
acceleration, resulting in knot motion away from the jet axis. This could be
due to bends in the jet close to the core or a widening of the jet flow. The
motion is reversed for distances $\sim 1$-2 pc, where there is negative
acceleration. This could be an indication of recollimation of the jet flow,
such that the jet width narrows or remains constant downstream of this point.
Stacking of the 43 GHz images of jets can be used to investigate the jet width
(e.g., Pushkarev et al., 2017; Kovalev et al., 2020; Casadio et al., 2021) and
test whether the acceleration trends are actually related to the jet width.

An important question is whether accelerations are related to the stationary
features in the jets, as one might expect if the features represent standing
shocks. To test this, Figure 16 displays the relationship between the
locations of each acceleration region and the closest stationary feature. The
gray dotted line depicts a 1:1 relationship between the two, i.e., where the
positions of acceleration and a stationary feature are the same. We have used
the standard deviation of the stationary feature position as the uncertainty
in the distance from the core to the closest stationary feature. We have
assigned a positional error of 0.02 mas ($\sim 0.2$ times our best resolution)
to the distance from the core to the acceleration region.

Figure 16: Left: Distances from the 43 GHz core to the acceleration region and
stationary feature closest to that acceleration region. The vertical segments
of data correspond to sources that contain multiple knots with acceleration
zones beyond the last stationary feature. The gray shaded region delineates
the portion of the figure that is expanded for better viewing in the right
panel. Error bars are shown in gray in both panels.

Very few pairs of acceleration region and stationary feature fall below the
1:1 relationship line. This indicates that most accelerations occur at or
beyond stationary features. An example of this trend is readily seen in the
FSRQ 0336$-$019, where knots _B2, B3_ , and _B4_ appear to accelerate after
passing through stationary features _A1_ and _A4_. Such behavior is expected
if the stationary features are structural components in the jet that can
interact with moving features, such as recollimation shocks. The
interpretation of recollimation shocks is further supported by the
acceleration distributions for projected distances $\sim 1$-4 pc, which are
predominately negative for $\eta_{\perp}$.

##  5 Physical Parameters of Jets

The trajectories, apparent velocities, and light curves derived for moving
features can be used to compute important physical parameters in the parsec-
scale jets of the blazars in our sample. These parameters are: variability
Doppler factor $\delta_{\text{var}}$, Lorentz factor $\Gamma$ of the motion,
viewing angle $\Theta_{\circ}$, and opening semi-angle $\theta_{\circ}$ of the
jet. We restrict the following analysis to the most “reliable” features, which
we define as those with a number of images in which the feature is visible
$N\geq 6$, and ejection times within our observed period of VLBA monitoring.
The resulting sample includes 191, 101, and 46 knots in FSRQs, BLs, and RGs,
respectively. Each knot can have multiple values of speed, depending on the
number of line segments used in the piece-wise linear fit. This leads to a
total of 262, 109, and 58 potential measurements of the physical parameters in
the jets of FSRQs, BLs, and RGs, respectively.

###  5.1 Variability Timescale and Doppler Factor

We follow the formalism developed by Jorstad et al. (2005) and modified in
Jorstad et al. (2017) to calculate the timescale of variability,
$\tau_{\text{var}}$, and variability Doppler factor, $\delta_{\text{var}}$,
for the moving knots in our sample. In this method, we assume that for the
majority of superluminal knots the flux density at 43 GHz declines from
radiative energy losses of the relativistic electrons rather than adiabatic
cooling. If this decay (calculated in the rest frame of the jet plasma) is
also shorter than the light-travel time, as we assume, the Doppler factor can
be determined as

$\delta_{\text{var}}\approx\frac{16sD_{\text{Gpc}}}{\tau_{\text{var}}(1+z)},$
(5)

where $D_{\text{Gpc}}$ is the luminosity distance to the source in Gpc,
$\tau_{\text{var}}$ is measured in years, $z$ is the redshift of the host
galaxy, and $s$ is the angular size of the knot in mas.

In order to estimate $\tau_{\text{var}}$, we make use of the fact that the
rise and decay portions of the millimeter-wave light curves of blazars can
both be roughly modeled by exponential profiles, so that $\ln S(t)\propto t$
(Terasranta & Valtaoja, 1994; Lister, 2001). Such millimeter-wave flares in
blazars occur roughly at the same time as the emergence of a new knot from the
core, and the decay of the flux density of a knot resembles that of the total
flux density (Savolainen et al., 2002). We thus approximate the flux decay of
the reliable knots in our sample with an exponential of the form
$\ln(S(t)/S_{\circ})=k(t-t_{\text{max}})$, where $t_{\text{max}}\leq t$ is the
epoch corresponding to the maximum flux density, $S_{\text{max}}$, $S_{\circ}$
is the flux density of a least-squares fit to the light curve decay at time
$t_{\text{max}}$, and $k$ is the slope of the fit. The timescale of
variability is then $\tau_{\text{var}}=|1/k|$ yr. While $\tau_{\text{var}}$
for most knots can be determined to within $\sim 10\%$, some light curves have
no well-defined maximum given the time coverage of the observations, which
leads to significant uncertainty in the timescale of variability. Because of
this, we remove knots from our reliable sample if
$\sigma_{\tau_{\text{var}}}>\tau_{\text{var}}/2$. We also cut from our
analysis motion segments beyond the first acceleration zone for knots with
complex motion. This reduced sample consists of estimates for the physical
parameters of 162 knots in 22 FSRQs, 80 knots in 11 BLs, and 31 knots in two
RGs. The remaining three sources either have no reliable moving knots with
defined epochs of ejections (the RG 0316+413 and the BL 1652+398) or no moving
features at all identified in this work (the BL 1959+650). While we use the
superluminal knots of the BL 1101+384 to determine its jet parameters, the low
apparent speeds suggest that we are looking at a sheath of the jet rather than
the region closer to the axis (see $\S$5.3).

Figure 17 shows the decay light curves and the fits used to calculate
$\tau_{\text{var}}$ for a knot in each subclass, shifted for clarity. The
distributions of the timescales of variability, corrected for redshift to
transform into the host galaxy frame, are similar for all three subclasses,
peaking at $\tau_{\text{var}}^{\text{s}}\lesssim 6$ months with the majority
($\gtrsim 80\%$) of all knots having a timescale
$\tau_{\text{var}}^{\text{s}}<1$ yr.

Figure 17: Light curves of moving knots B11 in the FSRQ 1510$-$089 (red
circles), B11 in the BL 0954+658 (blue squares), and K8 in the RG 0415+379
(green triangles). The dotted lines represent the least-square fits to the
flux density decays of the knots. The data and fits for B11 and K8 have been
shifted by $-1.5$ and $-3.0$ from their normalized values, respectively.

In order to test the assumption that the flux density decay of knots at 43 GHz
is due primarily to radiative losses rather than adiabatic expansion, we also
fit an exponential of the form
$\ln(a(t)/a_{\circ})=k_{\text{a}}(t-t_{\text{max}})$ to the knot sizes,
$a(t)$, over the same epochs as used to calculate the flux timescale of
variability. Here, $a_{\circ}$ is the size of the knot according to a least-
squares fit to the data at time $t_{\text{max}}$ (defined above) and
$\tau_{\text{a}}=|1/k_{\text{a}}|$ yr. Figure 18 shows the timescale of
variability of knot size compared to that of the flux density for individual
knots. We have classified knots as having reliable values of the timescale for
changes in size if $\tau_{\text{a}}\geq 2\sigma_{\tau_{\text{a}}}$, where
$\sigma_{\tau{\text{a}}}$ is the estimated uncertainty. In the figure, the
dotted gray line indicates the approximate relation between $\tau_{\text{a}}$
and $\tau_{\text{var}}$ for adiabatic losses for optically thin shocked gas
with $\alpha=0.7$ (Marscher & Gear, 1985). A substantial majority (188/287,
65.5%) of knots have $\tau_{\text{var}}<\tau_{\text{a}}$ (above the dotted
line), implying that the flux density decay is due predominantly to radiative
losses rather than adiabatic losses, and thus that our initial assumption is
valid.

Figure 18: Timescale of variability of the size vs. that of the flux density
for superluminal knots in our sample. “Reliable” estimates of
$\tau_{\text{a}}$ (for which $\tau_{\text{a}}\geq 2\sigma_{\tau_{\text{a}}}$)
are shown in black, while those do not meet the reliability criterion are
shown in red. The solid line marks the approximate relation for an optically
thin shock with $\alpha\sim 0.7$ (Marscher & Gear, 1985). Error bars are
omitted for visual clarity.

The angular size $a$ obtained in the model fits is the FWHM of a presumed
Gaussian brightness distribution. We consider a face-on disk to be a more
physically appropriate geometry, and therefore convert $a$ to the diameter $s$
of a face-on disk that has a similar interferometric visibility function,
$s=1.6a$, for use in Equation 5. We estimate $a$ as the average over the
epochs used to calculate $\tau_{\text{var}}$, so that
$a=\sum_{i=1}^{n}a_{i}/n$, where $a_{i}$ is the FWHM size of the knot at each
epoch $i$ and $n$ is the number of epochs in the decay phase of the light
curve. Table 12 presents the physical parameters of each knot segment as
follows: 1—source name; 2—knot designation; 3—timescale of variability of the
knot flux density, $\tau_{\text{var}}$, and its estimated $1\sigma$
uncertainty from the least-square fit to the decay of the light curve, in yr;
4—average size of the knot, $\langle a\rangle$, and its propagated $1\sigma$
uncertainty, in mas; 5—number of measurements used to calculate
$\tau_{\text{var}}$ and $\langle a\rangle$; $N_{t}$; 6—Doppler factor of
variability for each knot segment, $\delta_{\text{var}}$, and its $1\sigma$
uncertainty; 7—bulk Lorentz factor of each knot segment, $\Gamma$, and its
$1\sigma$ uncertainty (see $\S$5.2); and 8—intrinsic viewing angle of each
knot segment, $\Theta_{\circ}$, and its $1\sigma$ uncertainty, in degrees (see
$\S$5.2).

For comparison with past literature, Figure 19 (left) displays the
distributions of the Doppler factors for all reliable knots derived using the
above method. The Doppler factors of FSRQs and BLs are predominantly
$\delta_{\text{var}}<40$, with only eight and three segments having
$\delta_{\text{var}}>50$, respectively. These few knots are shown in the last
bins of each distribution.

To more appropriately compare between the subclasses, we show in Fig. 19
(right) the distributions of $\delta_{\mathrm{var}}$ for “typical” knots of
each source, where the “typical” knot has been determined through the method
outlined in $\S$5.3. We see that the $\delta_{\mathrm{var}}$ values in RGs are
subtantially lower than those of FSRQs and BLs. This result is statistically
significant $\sim 3\sigma$ when using a KS test to compare the FSRQ and RG
distributions ($\mathcal{D}=0.955,\ p=0.003$) but not when comparing the BL to
RG distribution ($\mathcal{D}=0.769,\ p=0.071$) due to the relatively low
number of RG sources in our sample. The FSRQ and BL distributions appear very
similar (according to a KS test; $\mathcal{D}=0.311,\ p=0.333$). However, this
simple summary statistic hides an important distinction between the
$\delta_{\mathrm{var}}$ values of FSRQs and BLs. More than half (8/13,
$61.5\%$) of the BLs have a typical $\delta_{\mathrm{var}}\lesssim 15$, while
in FSRQs more than half (12/22, $54.5\%$) of the sources have
$\delta_{\mathrm{var}}\gtrsim 15$. This higher probability for FSRQs to have
larger $\delta_{\mathrm{var}}$ values is consistent with their observed higher
brightness temperatures of jet features as compared to BLs, apparent in Figure
5.

Figure 19: Distributions of all calculated (left) and typical (defined in
$\S$5.3, right) Doppler factors of variability, $\delta_{\text{var}}$, derived
for reliable moving knots in FSRQs (top), BLs (middle), and RGs (bottom). The
typical $\delta_{\mathrm{var}}$ for the BL 1101+384 is shown in the last bin
to condense the axis.

###  5.2 Lorentz Factors and Viewing Angles

With the apparent speeds and Doppler factors derived for each knot, we can
independently calculate the Lorentz factor $\Gamma$ and viewing angle
$\Theta_{\circ}$ with the relations

$\Gamma=\frac{\beta_{\text{app}}^{2}+\delta_{\text{var}}^{2}+1}{2\delta_{\text{var}}};\hskip
7.11317pt\tan{\Theta_{\circ}}=\frac{2\beta_{\text{app}}}{\beta_{\text{app}}^{2}+\delta_{\text{var}}^{2}-1}.$
(6)

Table 12: Physical Parameters of Jet Features

Source | Knot | $\tau_{\text{var}}$ | $\langle a\rangle$ | $N_{t}$ | $\delta_{\text{var}}$ | $\Gamma$ | $\Theta_{\circ}$
---|---|---|---|---|---|---|---
|  | [yr] | [mas] |  |  |  | [$\arcdeg$]
0219+428 | B1 | $0.535\pm 0.161$ | $0.27\pm 0.02$ | 6 | $21.6\pm 6.7$ | $24.9\pm 5.6$ | $2.6\pm 0.7$
| B2 | $0.944\pm 0.300$ | $0.40\pm 0.02$ | 10 | $18.2\pm 5.9$ | $13.8\pm 1.8$ | $3.0\pm 1.3$
| B3 | $0.639\pm 0.245$ | $0.15\pm 0.02$ | 6 | $10.2\pm 4.0$ | $20.9\pm 5.4$ | $4.8\pm 1.0$
| B4 | $1.511\pm 0.330$ | $0.43\pm 0.03$ | 6 | $12.3\pm 2.8$ | $22.8\pm 3.5$ | $4.1\pm 0.5$
| B5 | $1.240\pm 0.212$ | $0.28\pm 0.02$ | 11 | $\phantom{0}9.9\pm 1.8$ | $19.8\pm 3.0$ | $5.0\pm 0.5$
0235+164 | B1 | $0.324\pm 0.068$ | $0.26\pm 0.02$ | 6 | $63.3\pm 13.9$ | $37.2\pm 5.8$ | $0.7\pm 0.2$
| B2 | $0.207\pm 0.023$ | $0.20\pm 0.01$ | 9 | $76.7\pm 9.4$ | $39.4\pm 4.5$ | $0.2\pm 0.1$
| B3 | $1.243\pm 0.119$ | $0.15\pm 0.01$ | 28 | $\phantom{0}9.3\pm 0.9$ | $\phantom{0}5.1\pm 0.4$ | $3.2\pm 0.7$
| B4 | $0.717\pm 0.102$ | $0.15\pm 0.01$ | 13 | $16.4\pm 2.5$ | $\phantom{0}9.1\pm 1.1$ | $2.1\pm 0.6$
| B5 | $0.493\pm 0.068$ | $0.23\pm 0.01$ | 13 | $37.9\pm 5.4$ | $20.0\pm 2.5$ | $0.7\pm 0.2$
| B6 | $0.361\pm 0.044$ | $0.28\pm 0.01$ | 13 | $62.5\pm 8.1$ | $31.6\pm 4.0$ | $0.2\pm 0.1$

Note. — (Table 12 (273 lines) is published in its entirety in the machine-
readable format. A portion is shown here for guidance regarding its form and
content.)

Figure 20, left, displays the distributions of $\Gamma$ for all reliable knots
of each subclass. The FSRQ and BL distributions both peak at $\Gamma\sim 10$,
but BLs have a small local peak in the $\Gamma=2$-$6$ bins, while the FSRQ
distribution has a longer high-$\Gamma$ tail. There is a single FSRQ knot for
which $\Gamma\geq 60$, _B1_ in 0528+134, with $\Gamma=81\pm 10$; it is placed
in the highest bin in the distribution. The RG distribution is comparatively
narrow and low-valued, with all but one of the knots having $\Gamma<10$. To
directly compare the Lorentz factors of FSRQs, BLs, and RGs, we utilize a KS
test on the distributions of $\Gamma$ for the “typical” knots of each source
(see $\S$5.3), shown in Figure 20 (right). Despite the low number of sources,
we continue to infer that RGs have the slowest knots ($\Gamma\lesssim 10$),
while FSRQs have, on average, faster knots, with a long tail ($\Gamma$
extending up to $\sim 60$). Both FSRQ and BL distributions are statistically
different from the RG distribution ($\sim 2\sigma$), with FSRQ vs. RG
$\mathcal{D}=0.955,\ p=0.003$, while the BL vs. RG $\mathcal{D}=0.846,\
p=0.036$). The KS test cannot distinguish between the FSRQ and BL
distributions, however, with $\mathcal{D}=0.269$ and $p=0.495$.

Figure 20: Distributions of all calculated (left) and typical (defined in
$\S$5.3, right) Lorentz factors, $\Gamma$, in FSRQs (top), BLs (middle), and
RGs (bottom), derived from the apparent speeds and variability Doppler factors
of reliable knots.

Figure 21 plots the distributions of the viewing angles for all reliable knots
in the three subclasses. It is immediately apparent that RGs have a very
different distribution than the FSRQs and BLs, peaking in the
$16\arcdeg$-$18\arcdeg$ bin, with all knots having $\Theta_{\circ}>10$. Minor
differences in can be seen between the FSRQ and BL distributions, such as the
longer high-angle tail of the BLs. The peaks are also different, with the FSRQ
distribution peaking between $0\arcdeg$ and $2\arcdeg$, while the BL
distribution peaks in the $2\arcdeg$-$4\arcdeg$ range. In both subclasses,
however, the majority ($>75\%$) of knots have $\Theta_{\circ}<10\arcdeg$.
Again to compare between the subclasses statistically, we use the viewing
angle of a “typical” knot (see $\S$5.3) for each source. The distributions for
each subclass are shown in Figure 21. The FSRQs are more centrally located
around the peak in viewing angle, $\Theta_{\circ}\sim 2-4\arcdeg$, than the
BLs, but with a KS test indicating $\mathcal{D}=0.231,\ p=0.678$. The RG
distribution is statistically different from the other two subclasses, with
$\mathcal{D}=0.955,\ p=0.003$ for the FSRQs vs RGs and $\mathcal{D}=0.769,\
p=0.071$ for the BLs vs RGs, despite the low number of sources.

Figure 21: Distributions of all calculated (left) and typical (defined in
$\S$5.3, right) viewing angles in FSRQs (top), BLs (middle), and RGs (bottom),
derived from the apparent speeds and variability Doppler factors of reliable
knots. Note the different scales in each panel.

###  5.3 Average Jet Physical Parameters and Opening Angles

Figure 22: Distributions of physical parameters $(a)\ \delta_{\text{var}}$,
$(b)\ \Gamma$, and $(c)\ \Theta_{\circ}$ of reliable knots in the RG 0415+379.
Red curves represent best-fit Gaussian distributions to the data. See the text
for details.

According to Table 12, for any individual source there is a diversity in knot
properties, which complicates the determination of the characteristic physical
conditions in the jet. The physical parameters of many jets are evidently
time-variable, perhaps stochastically determined from a distribution of
possible values. As an example, Figure 22 shows the distributions of the
physical parameters $\delta_{\text{var}},\ \Gamma$, and $\Theta_{\circ}$ of
the reliable knots of the RG 0415+379. The distributions can be approximated
by Gaussian functions, with means $\mu=2.6,\ 5.8,\ 16.3\arcdeg$ and standard
deviations $\sigma=0.7,\ 3.1,\ 5.0\arcdeg$, respectively, for
$\delta_{\text{var}},\ \Gamma$, and $\Theta_{\circ}$. However, these well-
determined distributions are the exception, since the number of knots observed
in our study is more limited for most of the sample. The majority of sources
have $N\leq 15$, too few to fit Gaussian distributions to the data accurately.
Given the knot ejection rates of the sources observed in our sample, another
$\sim 10$ years of monitoring would be required to observe enough knots to