MOJAVE. XVII. Jet Kinematics and Parent Population Properties of Relativistically Beamed Radio-loud Blazars

Relativistic jets from active galactic nuclei (AGNs) represent some of the most energetic known phenomena in the universe, and they played a key role in regulating galaxy formation at early epochs via feedback processes (Blandford et al. 2018). One of the most powerful tools for investigating these outflows is the Very Long Baseline Array (VLBA), which can be used to provide full polarization, submilliarcsecond scale imaging at radio wavelengths.

Since the VLBA's inauguration in 1994, we have carried out a long-term program to investigate the parsec-scale properties of several hundred of the brightest AGN jets in the northern sky. This effort started out as the 2 cm VLBA survey (Kellermann et al. 1998) and continued as the MOJAVE survey in 2002 with the addition of full polarization imaging of a complete flux density-limited sample. We have presented the results from MOJAVE in a number of papers in this series, including our most recent analysis of jet kinematics based on multiepoch data obtained between 1994 August 31 and 2013 August 20 (Lister et al. 2016).

In this paper we perform a new kinematics analysis that adds VLBA data taken up to 2016 December 26, and extends the number of AGN jets studied from 274 to 409. Most of the new AGNs were added to the MOJAVE program on the basis of their detection in GeV gamma-rays by the LAT instrument on board the Fermi observatory. We also update and expand our 1.5 Jy flux density-limited sample from 181 to 230 AGNs according to data from the RATAN 600 m telescope and Owens Valley Radio Observatory (OVRO) 40 m telescope monitoring observations at 15 GHz. This sample is now the largest and most complete radio-loud blazar sample to date, covering 75% of the entire sky. Using Monte Carlo simulations, we deconvolve the effects of Doppler boosting and Malmquist bias in this sample to uncover the intrinsic jet properties of the bright radio-loud quasar population.

The layout of the paper is as follows. In Section 2 we describe our VLBA observations and the new flux density-limited 1.5 Jy Quarter Century (1.5JyQC) AGN sample. We describe our Gaussian fitting of bright jet features and their apparent trajectories, and discuss our general findings on the parsec-scale jet kinematics of our sample in Section 3. In Section 4 we describe the best-fit parent population properties for 174 quasars in the 1.5JyQC sample based on Monte Carlo simulations. We use this best-fit simulation to describe the likely viewing angle, Lorentz factor, and Doppler factor distributions of bright radio-loud quasars.

Throughout this paper we adopt the convention ${S}_{\nu }\propto {\nu }^{\alpha }$ for spectral index α, and use the cosmological parameters Ω_m = 0.27, Ω_Λ = 0.73, and ${H}_{o}=71\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ (Komatsu et al. 2009).

The observational data set consists of 15 GHz VLBA observations of 409 AGNs obtained between 1994 August 31 and 2016 December 26 as part of the MOJAVE program, with supplementary data from the NRAO archive. These AGNs all have 15 GHz flux density ≳0.1 Jy, and have at least five VLBA epochs spaced in time. The epoch coverage and cadence varies considerably among the AGNs as they are members or candidate members of various radio and gamma-ray selected samples that have been added at various stages of the program (see Lister et al. 2018). We have previously presented the VLBA total intensity and polarization images in Lister & Homan (2005), Lister et al. (2009a, 2013, 2016), and Lister et al. (2018). These images are also available from our online data archive.¹³ We obtained observer frame values for the low energy (synchrotron) peak frequency from the literature or via the Agenzia Spaziale Italiana Data Center (ASDC) spectral energy distribution (SED) builder (Stratta et al. 2011). We list the overall properties of the AGNs in Tables 1 and 2. The latter contains 19 AGNs that have not been observed in the MOJAVE VLBA program, but are new additions to the new 1.5 Jy sample, as we describe in the next section.

2.1. The MOJAVE 1.5 Jy Quarter Century Sample

3.1. Gaussian Modeling

3.2. Jet Kinematics

As in our previous papers (Homan et al. 2009, 2015; Lister et al. 2009b, 2013, 2016), we analyze the kinematics of jet features using three methods: (i) a simple one-dimensional radial motion fit, (ii) a nonaccelerating vector fit in two (sky) dimensions, and (iii) a constant acceleration fit (for features with 10 or more epochs). We use the radial fit for diagnostic purposes only (see below), and do not tabulate those fit results here. In all cases, we assume the bright core feature (id = 0 in Table 3) to be stationary, and measure the positions of jet features at all epochs with respect to it.

We have modified our model slightly from our previous papers, and now fit for the sky position of each feature at a reference middle epoch t_mid, rather than fitting for the epoch of origin in the x (R.A.) and y (decl.) sky directions. Our new parameterization is as follows:

$$\begin{eqnarray}&&x(t)={x}_{\mathrm{mid}}+{\mu }_{x}(t-{t}_{\mathrm{mid}})+\displaystyle \frac{{\dot{\mu }}_{x}}{2}{\left(t-{t}_{\mathrm{mid}}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&y(t)={y}_{\mathrm{mid}}+{\mu }_{y}(t-{t}_{\mathrm{mid}})+\displaystyle \frac{{\dot{\mu }}_{y}}{2}{\left(t-{t}_{\mathrm{mid}}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where t_mid is the numerical mean of the first and last observation epoch dates for the feature being fitted, and μ_x and μ_y are the fitted angular speeds in each sky direction. For the vector fits, the accelerations ${\dot{\mu }}_{x}$ and ${\dot{\mu }}_{y}$ are fixed to zero, and for the radial motion fits, we used the alternate parameterization $r(t)={r}_{\mathrm{mid}}+{\mu }_{r}(t-{t}_{\mathrm{mid}})$, where r(t) is the radial distance from the core feature at time t.

We made radial and vector motion fits using all of the available data from 1994 August 31 to 2016 December 26 on 1744 robust jet features in 382 jets. There were 27 jets in which we could not identify any robust features because of a lack of sufficiently strong downstream jet flux or a suitably stable core feature, or insufficient spatial resolution. We are carrying out a follow-up 43 GHz multiepoch VLBA study on several of these jets.

In Table 4 we list the results of the vector motion fits. Because of the nature of our kinematic model, which naturally includes the possibility of accelerated motion, we did not estimate ejection epochs (Column 12) for any features where we could not confidently extrapolate their motion to the core. Jet features for which we list an ejection epoch had the following properties: (i) significant motion (μ ≥ 3σ_μ), (ii) no significant acceleration, (iii) a velocity vector direction ϕ within 15° of the outward radial direction to high confidence, i.e., $| \langle \vartheta \rangle -\phi | +2\sigma \leqslant 15^\circ$, where ϑ is the mean position angle, (iv) an extrapolated position at the ejection epoch no more than 0.2 mas from the core, and (v) a fitted ejection epoch that differed by no more than 0.5 yr from that given by the radial motion fit.

Table 4.  Vector Motion Fit Properties of Jet Features

			$\langle S\rangle$	$\langle R\rangle$	$\langle {d}_{\mathrm{proj}}\rangle$	$\langle \vartheta \rangle$	ϕ	$| \langle \vartheta \rangle -\phi |$	μ	βapp			αm	δm
Source	I.D.	N	(mJy)	(mas)	(pc)	(deg)	(deg)	(deg)	(μas yr−1)	(c)	tej	tmid	(μas)	(μas)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)
0003+380	1	8	5	4.23	15.36	120.7	96 ± 17	24 ± 17	158 ± 43	2.30 ± 0.63	⋯	2008.81	3691 ± 74	−2169 ± 80
0003+380	2	6	19	1.78	6.45	112.6	120.1 ± 3.1	7.5 ± 3.1	317 ± 25	4.61 ± 0.36	⋯	2007.71	1662 ± 29	−694 ± 11
0003+380	4	5	16	1.25	4.53	114.9	205 ± 14	90 ± 14b	39 ± 10	0.57 ± 0.15	⋯	2009.54	1130 ± 11	−527 ± 14
0003+380	5	8	40	0.75	2.71	117.5	21 ± 89	96 ± 89	2.7 ± 7.6	0.04 ± 0.11	⋯	2010.26	663 ± 20	−342 ± 10
0003+380	6	10	98	0.39	1.43	115.4	335 ± 46	141 ± 46	12.7 ± 8.4d	0.19 ± 0.12	⋯	2009.90	350 ± 22	−158 ± 19
0003−066	2	5	222	1.05	5.12	322.9	226.3 ± 4.9	96.6 ± 5.0b	191 ± 15	4.09 ± 0.33	⋯	1997.80	−585.9 ± 8.9	883 ± 37
0003−066	3	9	119	2.82	13.73	296.9	284.8 ± 4.7	12.1 ± 4.8	250 ± 39	5.36 ± 0.83	⋯	1999.33	−2375 ± 98	1237 ± 41
0003−066	4a	26	120	6.61	32.23	285.6	284 ± 11	2 ± 11	41 ± 14	0.87 ± 0.29	⋯	2004.83	−6326 ± 60	1768 ± 22
0003−066	5a	14	1031	0.70	3.40	10.7	350.9 ± 5.3	19.9 ± 5.5b	88.1 ± 4.3	1.888 ± 0.091	⋯	2004.37	138 ± 18	634.1 ± 9.0
0003−066	6a	10	97	1.01	4.92	290.2	210 ± 15	81 ± 15b	55 ± 17	1.18 ± 0.37	⋯	2003.78	−941 ± 15	359 ± 33

Notes. Columns are as follows: (1) B1950 name, (2) feature number, (3) number of fitted epochs, (4) mean flux density at 15 GHz in milliJansky, (5) mean distance from core feature in milliarcseconds, (6) mean projected distance from core feature in parsec, (7) mean position angle with respect to the core feature in degrees, (8) position angle of velocity vector in degrees, (9) offset between mean position angle and velocity vector position angle in degrees, (10) proper motion in μas yr⁻¹, (11) apparent speed in units of the speed of light, (12) estimated epoch of origin, (13) date of reference (middle) epoch used for fit, (14) fitted R.A. position with respect to the core at the middle epoch in μas, (15) fitted decl. position with respect to the core at the middle epoch in μas. A question mark indicates a feature whose motion is not consistent with outward, radial motion but for which the possibility of inward motion and its degree of non-radialness are uncertain.

^aAcceleration model fit indicates significant accelerated motion. ^bFeature has significant non-radial motion according to the vector motion fit. ^cFeature has significant inward motion according to the vector motion fit. ^dFeature has slow pattern speed.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

A total of 881 of the robust features met the ≥10 epoch criterion for an acceleration fit, and we tabulate these results in Table 5. The majority (59%) of these well-sampled features display either significant acceleration or non-radial motion, which confirms our previous finding that accelerated motions are common in parsec-scale AGN jets (Lister et al. 2016).

Table 5.  Acceleration Fit Properties of Jet Features

		ϕ	$| \langle \vartheta \rangle -\phi |$	μ	βapp	$\dot{\mu }$	ψ	${\dot{\mu }}_{\perp }$	${\dot{\mu }}_{\parallel }$	αm	δm
Source	I.D.	(deg)	(deg)	(μas yr−1)	(c)	(μas yr−2)	(deg)	(μas yr−2)	(μas yr−2)	(μas)	(μas)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
0003+380	6	333 ± 44	142 ± 44	13.4 ± 8.6	0.20 ± 0.12	9.8 ± 8.4	309 ± 53	−4.0 ± 9.4	9.0 ± 9.0	371 ± 33	−175 ± 28
0003-066	4a	277.3 ± 3.8	8.3 ± 3.8	50.9 ± 5.3	1.09 ± 0.11	28.5 ± 2.3	73.7 ± 3.1	11.4 ± 2.1	−26.1 ± 2.5	−6582 ± 32	1693 ± 20
0003-066	5a	353.9 ± 3.0	16.8 ± 3.1b	87.2 ± 4.4	1.868 ± 0.093	26.6 ± 4.9	274 ± 10	−26.3 ± 4.9	4.5 ± 4.8	199 ± 15	630 ± 14
0003-066	6a	211.3 ± 9.6	78.9 ± 9.6b	54 ± 11	1.16 ± 0.24	65 ± 16	336 ± 11	54 ± 13	−37 ± 18	−901 ± 16	268 ± 35
0003-066	8a	290.7 ± 1.6	3.5 ± 1.6	330.4 ± 9.7	7.08 ± 0.21	67 ± 12	127 ± 10	−19 ± 12	−64 ± 12	−2444 ± 30	1121 ± 28
0003-066	9	295.2 ± 4.1	7.5 ± 4.3	278 ± 20	5.96 ± 0.42	99 ± 35	110 ± 22	9 ± 37	−99 ± 35	−1769 ± 52	582 ± 53
0010+405	1	340.7 ± 4.4	11.9 ± 4.4	432 ± 42	6.99 ± 0.68	44 ± 83	147 ± 76	11 ± 53	−43 ± 70	−4259 ± 76	6991 ± 107
0010+405	2	9 ± 123	41 ± 123	2 ± 14	0.04 ± 0.23	4 ± 22	152 ± 123	2 ± 21	−3 ± 23	−898 ± 30	1470 ± 48
0010+405	3	138 ± 83	170 ± 83	2.5 ± 5.4	0.041 ± 0.088	6.9 ± 6.1	99 ± 57	−4.3 ± 8.8	5.4 ± 9.0	−493.6 ± 9.5	783 ± 15
0010+405	4	113 ± 98	145 ± 98	1.4 ± 4.5	0.022 ± 0.072	5.2 ± 8.9	318 ± 69	−2.2 ± 7.1	−4.7 ± 8.1	−240.6 ± 9.0	382 ± 14

Notes. Columns are as follows: (1) B1950 name, (2) feature number, (3) proper motion position angle in degrees, (4) offset between mean position angle and proper motion position angle in degrees, (5) proper motion in μas yr⁻¹, (6) apparent speed in units of the speed of light, (7) acceleration in μas yr⁻², (8) acceleration vector position angle in degrees, (9) acceleration perpendicular to velocity direction in μas yr⁻², (10) acceleration parallel to velocity direction in μas yr⁻², (11) fitted R.A. position with respect to the core at the middle epoch in μas, (12) fitted decl. position with respect to the core at the middle epoch in μas. A question mark indicates a feature whose motion is not consistent with outward, radial motion but for which the possibility of inward motion and its degree of non-radialness are uncertain.

^aFeature shows significant accelerated motion. ^bFeature shows significant non-radial motion according to the acceleration fit. ^cFeature shows significant inward motion according to the acceleration fit.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

In Figure Set 1 we plot the angular separation of features from the core in each jet versus time. The robust features are plotted with filled colored symbols and solid lines representing the fit. The feature identification number is overlined if the acceleration model was fit and yielded a >3σ acceleration. An underlined identification number indicates a feature with non-radial motion, i.e., its velocity vector did not point back to the core location within the errors. We plot the individual trajectories and fits on the sky for all the robust features in Figure Set 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (409 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (1743 panels)

Tarball of figure set images

3.2.1. Pattern Speeds

In a previous kinematic study (Lister et al. 2013), we found that in many individual AGN jets, there is no single apparent speed β_app = v_app/c at which bright features propagate downstream. Instead, there is typically a single characteristic speed with a modest spread around this value. Because trackable features emerge only every few years in most bright blazar jets, continuous monitoring periods of a decade or more are needed to establish the characteristic speed of a jet, and whether any individual feature may have an atypically low pattern speed (see also a recent analysis of MOJAVE kinematics results by Plavin et al. 2018).

In Figure 3 we show the distribution of speed differences from the jet's median speed for 436 features in 26 jets that have 10 or more robust features. This plot contains nearly twice as many jet features as our previous kinematic study, and is qualitatively similar. Most features lie within ±40% of the jet's median speed. There is also a small tail consisting of atypically fast features. The jet with the largest range of speeds is 4C +15.05 (0202+149), which has 10 features with apparent speeds ranging from 0.1 c to 16 c.

Zoom In Zoom Out Reset image size

We have identified 55 features in 42 AGN jets that have appreciably slower speeds than other features in the same jet. Our specific criteria are that the feature (i) does not have a >3σ acceleration, (ii) has an angular speed smaller than 20 $\mu \mathrm{as}\,{\mathrm{yr}}^{-1}$, and (iii) has a speed at least 10 times slower than the fastest feature in the same jet. Figure 4 shows the distribution of projected distance from the core for 53 slow pattern speed features in 40 AGNs with known redshifts. The vast majority are located within 4 pc of the core feature (∼100 pc deprojected, given typical viewing angles <2°). This is consistent with the 43 GHz VLBA survey of 36 AGNs by Jorstad et al. (2017), who found 21% of jet features to be quasi-stationary, with most located at projected core distances below 3 pc.

Zoom In Zoom Out Reset image size

Of the 1744 robust jet features that we have studied, only 44 (2.5%) have velocity vectors that are directed inward toward the core feature. We might expect to see rare instances of apparent inward motion when a feature moving along a curved trajectory crosses our line of sight (e.g., as in the case of 4C +39.25; Alberdi et al. 2000). It is also possible that small changes in the brightness distribution of a large diffuse feature may alter its best-fit Gaussian centroid location, creating apparent inward motion. We note two instances (feature id = 1 in 87 GB 061258.1+570222 and id = 1 in 8C 1944+838) where this may be the case. Inward motion can also result from incorrect identification of the core feature, or variable structure near the core that is below the angular resolution of our observations that may alter the fitted core location. We note that in 16 of 33 AGN jets with inward motion, the inward-moving feature is the closest feature to the core, and four of these jets (all associated with BL Lac objects: UGC 00773, 3C 66A, Mrk 421, and ON 325) have more than one close-in inward-moving feature.

3.2.2. Speed Distributions

We have calculated maximum and median speed statistics for the jets in our sample using the method described in Lister et al. (2013). For accelerating features, we note that the speeds are determined at the middle epoch, and thus may not represent the maximum speed attained by the feature. In the case of two AGNs for which we could not identify any robust features (AO 0235+164 and 1ES 1959+650), we adopted maximum speeds from the literature according to VLBA observations made at other wavelengths. We plot the distributions of these statistics in Figure 5. Slow apparent speeds are common, with very few measured speeds above 30 c. As discussed by Vermeulen & Cohen (1994) and Lister & Marscher (1997), the shape of the distributions is incompatible with all jets having the same bulk Lorentz factor, and instead suggests a power-law parent distribution that is weighted toward slow speeds. Single-valued parent Γ distributions predict an excess of high apparent jet speeds, while Gaussian Γ distributions do not reproduce the gradual fall-off in the number of jets with higher apparent speeds.

Zoom In Zoom Out Reset image size

3.2.3. Statistical Trends

In Figure 6 we plot maximum apparent jet speed versus rest-frame synchrotron SED peak frequency. The plot includes AGNs from our survey, as well as those from Piner & Edwards (2018) and Jorstad et al. (2017). AGNs with $\lt 3\sigma$ maximum speeds are indicated with upper limit symbols. The crosses indicate BL Lacs with no known redshift, and their extents correspond to lower and upper redshift limits published in the literature. For clarity, we have omitted BL Lacs for which the redshift limits give a possible range of β_app greater than 20. There is a clear upper envelope to the distribution, with the highest jet speeds being found only in AGNs with low synchrotron peak frequencies.

Zoom In Zoom Out Reset image size

The filled symbols indicate AGNs that have been detected at TeV gamma-ray energies with the airshower telescopes VERITAS, HESS, or MAGIC. The large fraction of high synchrotron peaked (HSP) AGNs that are TeV-detected in this plot is a selection effect because these have been specifically targeted for long-term VLBA kinematic study by Piner & Edwards (2018). The ISP AGNs have been targeted in MOJAVE on the basis of their detection at GeV energies by Fermi, while most of the low synchrotron peaked (LSP) AGNs are from the radio-selected MOJAVE sample (Lister et al. 2016). Although a fast jet speed does not guarantee a TeV detection, it does appear to be a minimum requirement for the intermediate- and low-synchrotron peaked AGNs. This implies a direct connection between the bulk jet speed measured on parsec scales and the Doppler boosting level of the TeV emission. Of the 14 non-HSP TeV detected AGNs in Figure 6, only three have maximum jet speeds below 6 c. Two of these (3C 84 and M 87) are very nearby (<75 Mpc) radio galaxies, and the third (TXS 0506+056) is an unusual ISP BL Lac with a measured maximum speed of 0.98c ± 0.3c that lies within the sky error circle of a high-energy neutrino event detected in 2017 (IceCube Collaboration et al. 2018).

The interpretation of parsec-scale AGN jet kinematic studies presents a challenge in the sense that the individual objects that are most easily studied (i.e., high flux density, with proper motions observable on time periods of approximately a few years) are blazars, whose selection is highly affected by Doppler bias (Scheuer & Readhead 1979). In principle, the observed redshift, luminosity, and apparent speed distributions of a complete flux density-limited jet sample can be used to recover the intrinsic properties of the blazar parent population, but the Doppler selection effects need to be carefully accounted for. Vermeulen & Cohen (1994) and Lister & Marscher (1997) have shown that this can be done analytically only in the case of very simplistic, non-realistic assumptions. These include a non-evolving single power-law luminosity function (LF) and a single-valued or uniform distribution of bulk Lorentz factors, neither of which provides satisfactory fits to the data. The typical approach (e.g., Lister & Marscher 1997; Bloom 2008; Lister et al. 2009b; Giommi et al. 2012; Liodakis & Pavlidou 2015) has been to generate simulated flux density-limited samples from jet parent populations whose properties are drawn from specified probability distributions, and find the set of distribution parameters that best fit the data. In this section we carry out this type of Monte Carlo analysis on our MOJAVE data, based on the method of Lister & Marscher (1997).

4.1. Simulated Jet Properties

The observed flux density S_ν from a spherical optically thick source of radiation with an isotropically emitted rest-frame luminosity L_ν ∝ ν^α, moving with bulk Lorentz factor Γ at an angle θ to the line of sight and located at redshift z (with corresponding luminosity distance D_L(z)) is (e.g., Blandford & Königl 1979; Condon & Matthews 2018)

$$\begin{eqnarray}&&{S}_{\nu }=\displaystyle \frac{{L}_{\nu }{\delta }^{p}\,{\left(1+z\right)}^{(1+\alpha )}}{4\pi {D}_{L}^{2}(z)},\end{eqnarray} \tag{ 3 }$$

where ν is the observing frequency and L_ν is the luminosity emitted in the jet frame at that same frequency. The exponent p of the Doppler factor

$$\begin{eqnarray}&&\delta ={\left[{\rm{\Gamma }}-\sqrt{({{\rm{\Gamma }}}^{2}-1)}\cos \theta \right]}^{-1}\end{eqnarray} \tag{ 4 }$$

is p = 3 − α in the scenario described above. However, in the case of a continuous jet made up of many such spheres, one cannot distinguish the lifetimes of the individual emitting particles, and a time dilation factor of δ is no longer applicable, hence $p=2-\alpha$ (Cawthorne 1991).

The minimum properties required to simulate the observed flux density of an AGN jet are therefore z, L_ν, Γ, θ, α, and p. Actual AGN jets present complications in terms of the geometry of their emitting regions, optical depth variations, and flow accelerations, but the highest contributions to the observed flux density will come from regions where the synchrotron emission coefficient is highest, and where the Doppler factor is largest (e.g., $\theta \lesssim {\cos }^{-1}\beta$, where β is the flow velocity in units of the speed of light). Stacked-epoch MOJAVE VLBA images of blazars show mainly conical jet profiles (Pushkarev et al. 2017) in which adiabatic expansion and synchrotron losses exponentially reduce the electron energies and magnetic field strength with distance down the jet (e.g., Konigl 1981). The bulk of the synchrotron emission therefore originates near the base of the jet, as confirmed by very long baseline interferometry (VLBI) morphologies that typically consist of a bright optically thick core feature accompanied by a much weaker jet. The exceptions to this are (i) young AGN jets of the compact symmetric object/gigahertz peaked spectrum class, which have high-luminosity radio lobes that are interacting with the interstellar medium of the host galaxy (O'Dea 1998), and (ii) rare instances where a bent downstream jet flow crosses the line of sight and experiences maximum Doppler boosting (e.g., 4C +39.35, Alberdi et al. 2000).

There are therefore good reasons to expect that a simulated population where each jet consists of a single (core) emitting region can provide a good representation of a suitably chosen blazar sample. The 1.5JyQC sample is well-suited in several respects, as it is a complete flux density-limited sample selected at high radio frequency, where the relative flux density contribution of the steep-spectrum downstream jet emission is low compared to the (typically flat-spectrum) core. It is also selected on the basis of VLBI flux density, which includes no contribution from any large kiloparsec-scale emission. Any contaminating CSO/GPS sources can be rejected on the basis of available spectral and morphological information, and most importantly, the sample is large enough to statistically constrain the best-fit parameters of the Monte Carlo simulations. After dropping two GPS quasars (PKS B0742+103 and OI −072) and six AGNs with no optical spectral information, there are 174 1.5JyQC quasars with redshift z ≥ 0.15 suitable for comparison with our simulations.

4.2. Simulation Parameters

Our simulation method is to generate a parent population of jets drawn from specified redshift, Lorentz factor, radio luminosity, and viewing angle distributions, calculate their predicted flux densities, and retain those jets that exceed the specified 1.5 Jy flux density limit. Because the 1.5JyQC sample includes all AGNs above decl. −30° known to have exceeded 1.5 Jy over a 25 yr period, we do not include any flux variability in our simulations, but instead compare our simulated jet flux densities to the maximum jet flux density for each AGN measured during the 1.5JyQC selection period (column 5 of Table 2).

4.2.1. Luminosity Function

Despite many studies on the radio LFs of AGNs, there is still no consensus on whether radio-loud AGN LFs evolve with lookback time in a manner consistent with increasing number density, increasing luminosity, or a mixture of both (Best et al. 2014; Smolčić et al. 2017; Yuan et al. 2018). There are also indications that lower-power (i.e., FR I) AGNs may evolve differently than the high-power (FR II) population (Rigby et al. 2008). Given these uncertainties, we have adopted a simple pure luminosity evolution parameterization for flat spectrum radio quasars used by Ajello et al. (2012) and Mao et al. (2017):

$$\begin{eqnarray}&&{\rm{\Phi }}(L,z)\propto {\rm{\Phi }}(L/e(z)),\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&e(z)={\left(1+z\right)}^{k}{e}^{z/\eta },\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{\rm{\Phi }}(L/e(z=0))\propto {L}^{\gamma }.\end{eqnarray} \tag{ 7 }$$

Our approach is to find the best-fit values of γ, η, and k using the MOJAVE data. We restrict our comparisons to quasars in the 1.5JyQC sample only, given the possibility that the BL Lac objects may be drawn from a different (i.e., lower power, or FR I) parent population (Urry & Padovani 1995). We set the lower limit on the parent LF at 10²⁴ W Hz⁻¹ based on the least powerful known FR II radio galaxies (e.g., Antognini et al. 2012).

4.2.2. Redshift Distribution

By adopting a pure luminosity evolution model, we assume that the parent jet population has a constant comoving density with redshift. All of the 1.5JyQC quasars have redshifts greater than 0.15, with the exception of TXS 0241+622 (z = 0.045). In order to avoid small number statistics in this nearby volume of space, we drop this AGN from our data comparisons and set the lower redshift limit of our simulation to z = 0.15. Because the form of LF evolution is not well known at very high redshift, we set the upper redshift limit in our simulations to that of the highest redshift 1.5JyQC quasar: OH 471 (z = 3.4).

4.2.3. Bulk Lorentz Factor Distribution and Doppler Boosting Index

Because of the strong selection biases associated with Doppler boosting, any large flux density-limited jet sample should contain some jets with the maximum Lorentz factor in the population (viewed at small θ). In the MOJAVE sample the fastest instantaneous measured jet speed is approximately 50 c for an accelerating feature in the jet of PKS 0805−07 (Lister et al. 2016), which corresponds to a Γ_max ≃ 50. In light of our discussion of the observed apparent velocity distributions in Section 3.2.2, we adopt a power-law Lorentz factor distribution for our simulated jets of the form $N({\rm{\Gamma }})\propto {{\rm{\Gamma }}}^{b}$, where b is a free parameter with values less than zero and Γ ranges from 1.25 to 50. The lower limit on Γ ($\beta \simeq 0.6c$) is based on a Bayesian analysis of the relative prominence of radio cores and kiloparsec-scale jets in FR II radio sources by Mullin & Hardcastle (2009). We assume no evolution of the jet Lorentz factor distribution with redshift.

The brightest radio-loud AGN cores are known to have a range of spectral indices with a mean value α = 0.22 ± 0.03 (Hovatta et al. 2014); however, the intrinsic distribution is not well known because of the difficulty of deconvolving relativistic beaming and projection effects. The spectral index enters into the simulated jet flux density via the small (1 + z)^1+α k-correction, and, more importantly, the Doppler boost index p. Any spread of α in the parent population will be effectively smoothed out in the observed LF, so we fix α = 0 for all our simulated jets and assume continuous jet emission such that p = 2. We discuss other fixed values of α in Section 4.4.

The Monte Carlo analysis of Lister & Marscher (1997) included the possibility of an intrinsic correlation between jet Lorentz factor and synchrotron luminosity of the form L ∝ Γ^ξ. They found that both ξ = 0 and $\xi \ne 0$ models produced very similar fits to the Caltech-Jodrell Flat-Spectrum AGN sample data. As we will show in Section 4.4, we are able to obtain good fits to the 1.5JyQC quasar sample assuming no L − Γ correlation, so we explore only the ξ = 0 case in this paper.

4.3. Simulation Procedure

In order to search for the best-fit parent population parameters, we constructed a grid of simulations with equally spaced parameter values spanning the ranges listed in Table 6. The procedure used to create each simulation in the grid is as follows:

Table 6.  Monte Carlo Jet Model Parameters

Jet Property	Distribution	Fixed Parameters	Free Parameter Ranges
Lorentz Factor	N(Γ)dΓ ∝ Γb	Γmin = 1.25	−1.8 ≤ b ≤ −0.2, step = 0.2
		Γmax = 50
Luminosity Function	Φ(L, z) ∝ Φ(L/e(z))	Lmin = 1024 W Hz−1	−0.65 ≤ η ≤ −0.25, step = 0.05
	e(z) = (1 + z)kez/η	Lmax = 1031 W Hz−1	4.5 ≤ k ≤ 8.5, step = 0.5
	Φ(L/e(z = 0)) ∝ Lγ		−3.2 ≤ γ ≤ −2.4, step = 0.1
Beamed Luminosity	P = Lδp	p = 2 + α	α = −0.5, 0, 0.22
Viewing Angle	$p(\theta )d\theta =\sin \theta$	θmin = 0°	⋯
		θmax = 90°

Download table as:  ASCIITypeset image

(i) Select values for b, γ, k, and η.

(ii) Generate z, L_ν, Γ, and θ values for a single jet from the probability distributions listed in Table 6.

(iii) Calculate the observed flux density of the jet according to Equation (3). We ignore any contribution from the counter-jet because it will be negligible for AGNs in a highly Doppler-biased sample (see Section 4.4.1).

(iv) If S_ν ≥ 1.5 Jy, keep the simulated jet.

(v) Repeat steps (ii) through (iv) until a sample of jets 10 times larger than the 1.5JyQC comparison sample is obtained, and record the total size of the parent population needed to produce this sample.

By creating samples larger than the data sample in step (v), we reduce the amount of statistical fluctuations associated with selecting a relatively small number of bright AGNs from a very large parent population. In doing so, we are effectively creating simulated jet samples from 10 universes and are comparing the mean properties of these samples to the data.

4.4. Comparisons to MOJAVE Data

For each simulation in the four-dimensional parameter grid (b, γ, k, and η) we compared the simulated flux density, redshift, radio luminosity (P_ν), and apparent velocity distributions to the 1.5JyQC sample of 174 quasars using the Anderson–Darling (A-D) test. The latter is a non-parametric test that assesses whether two samples are drawn from different parent populations and is sensitive to a wider variety of possible distribution differences than the frequently used Kolmogorov–Smirnov test (Engmann & Cousineau 2011). Our method was to randomly select a sample of 174 jets from the simulation and perform the A-D tests against the 1.5JyQC sample. We repeated this process 10 times and recorded the median A-D test probabilities p_S, p_z, p_P, and p_βapp, corresponding to the probability of the null hypothesis that the simulated and 1.5JyQC distributions are drawn from the same parent population. Because the completeness of the observational data is high (100% for S, z, and P, and 87% for β_app), we did not use any bootstrapping procedures to simulate the missing data.

In Figure 7 we plot two-dimensional projections of the grid parameter space, where the false color corresponds to the maximum value of p_z for any simulation having that particular parameter combination. The 1.5JyQC redshift distribution serves to constrain the parameter space to a limited number of k − η (LF evolution parameter) combinations, as seen in the lower right panel. The top row of plots in Figure 7 indicates, however, that the 1.5JyQC redshifts can be well-reproduced with many different combinations of the parent LF and the Lorentz factor distribution parameters.

Zoom In Zoom Out Reset image size

The two-dimensional projections in Figure 8, in which the false-color corresponds to the maximum values of p_βapp, serve to further constrain the region of viable parameter space for the simulations. The β_app distribution is best fit with simulations with k > 5.5. Also, the values of k and η that provide the best fits to the 1.5JyQC apparent velocity distribution (lower right panel) yield relatively poor fits to the observed luminosity distribution.

Zoom In Zoom Out Reset image size

Within the full grid, the simulation with the highest A-D probability summed over all four observable quantities has b = −1.4, γ = −3.1, k = 8.0, and η = −0.35 (model A). There are no other simulations in the grid that have an A-D probability greater than 0.4 in all four quantities. We investigated the effect of random statistical outliers on the A-D probability values for this best-fit simulation by first creating a simulated flux density-limited sample of 174,000 jets (i.e., 1000 universes), then selecting a random subset of 174 jets to compare with the 1.5JyQC data. After repeating the random subset selection 1000 times, the standard deviations on p_S, p_z, p_P, and p_βapp were 0.25, 0.3, 0.25, and 0.2 respectively. We therefore consider any simulation that has all four A-D probabilities within 1σ of those of the best-fit simulation to also be an acceptable fit to the data.

In Figure 9 we show a corner plot with false color indicating the number of acceptable best-fit simulations having particular parameter combinations. On the basis of the plot, we find acceptable fits for the parameter ranges −1.6 ≤ a ≤ −1.2, $-3.2\leqslant \gamma \leqslant -2.8$, 7.5 ≤ k ≤ 8, and $-0.35\leqslant \eta \leqslant -0.30$.

Zoom In Zoom Out Reset image size

We constructed two additional simulation grids to investigate whether better fits could be obtained using a fixed value of α = +0.22 (corresponding to a Doppler boosting index of p = 1.78), and α = −0.5 (corresponding to p = 2.5). The best-fit simulation in the p = 1.78 case (model B in Table 7) gave acceptable fits to the flux density, redshift, and luminosity distributions, but provided a relatively poor fit to the apparent speed distribution. Although the best-fit simulation in the p = 2.5 grid (model C) provided a good fit to the apparent speed distribution, none of the simulations in the grid gave A-D probabilities greater than 0.03 in all four observable parameters simultaneously.

Table 7.  Best-Fit Monte Carlo Grid Simulations

Model	Parameter						Anderson–Darling Test Probabilities
Name	b	γ	k	η	α	p	pS	pz	pP	pβapp
A	−1.40	−3.1	8.0	−0.35	0	2	0.65	0.43	0.52	0.40
B	−1.00	−2.6	8.5	−0.30	0.22	1.78	0.70	0.67	0.24	0.10
C	−1.40	−3.2	7.0	−0.35	−0.5	2.5	0.13	0.044	0.026	0.46

Note. The simulation parameters are defined in Table 6. Model A has the highest overall Anderson–Darling test probability sum p_S+p_z+p_P+p_βapp of any grid simulation.

Download table as:  ASCIITypeset image

4.4.1. Best-fit Parent Population Properties

In Figure 10 we show the distributions of observable quantities for the 1.5JyQC quasar sample (red lines), as well as our best-fit (model A) simulation. The blue bands represent 1σ ranges on the bin values that we derived by producing a simulation 1000 times the size of the 1.5JyQC, and then randomly choosing a subsample of 174 jets from it, repeating the latter step 10,000 times. We note that the simulation plotted in Figure 10 provides the best overall fit to the data; however, other combinations of fit parameters gave better fits to individual observable quantities. We have scaled the simulated apparent speed distribution in the top right panel by a factor of 151/174 = 0.87 to take into account the 23 missing jet speeds in the 1.5JyQC quasar sample.

Zoom In Zoom Out Reset image size

We plot the distributions of several intrinsic (indirectly observable) quantities from our best-fit simulation A in Figure 11. As expected from Doppler orientation bias, nearly all of the quasar jets in the 1.5JyQC sample are predicted to have viewing angles less than ∼10° from the line of sight, with the distribution peaking at 2°. The bottom left panel shows the distribution in terms of the critical angle ${\theta }_{\mathrm{cr}}={\sin }^{-1}(1/{\rm{\Gamma }})$, and indicates that the most likely viewing angle is not θ_crit as commonly cited in the literature, but approximately half of this value (e.g., Vermeulen & Cohen 1994; Lister & Marscher 1997; Cohen et al. 2007). The top middle panel shows the Lorentz factor distribution, which is broadly peaked between Γ ≃ 5 and Γ ≃ 15, with a rapid falloff past Γ = 20. The breadth of the Γ distribution indicates that adopting a single value of Γ = 10 for all blazars is not well supported by the observational data. The Doppler factor distribution has a similar shape to the Γ distribution, and peaks at δ ≃ 10, declining rapidly past δ ≃ 30.

Zoom In Zoom Out Reset image size

Liodakis et al. (2018) recently carried out a Bayesian light-curve analysis of OVRO 15 GHz monitoring data on the original 1.5 Jy sample and calculated variability Doppler factors and distributions of Lorentz factor and viewing angle. We find a high degree of consistency between these distributions for the 1.5 Jy quasars and those of our best-fit Monte Carlo simulation in Figure 11.

In the bottom middle panel we plot the distribution of intrinsic (unbeamed) luminosity L. Although the intrinsic parent LF peaks at L ≃ 10²⁴ W Hz⁻¹, most of the jets in the simulated flux density-limited sample have intrinsic (unbeamed) luminosities roughly an order of magnitude higher due to the combined effects of Doppler and Malmquist bias. This implies that the parent population of the brightest radio quasars consists of powerful FR II radio galaxies with a relatively narrow range of unbeamed 15 GHz radio luminosity between ∼10²⁵ and ∼10²⁶ W Hz⁻¹.

Liodakis et al. (2017) used Monte Carlo simulations to investigate the predicted distribution of jet-counterjet flux density ratios due to relativistic beaming in flux density-limited blazar samples. In the bottom right panel of Figure 11 we plot the distribution of this quantity for our best-fit model. We obtain very similar results, with most jets having ratios of 10⁴–10⁷. These are much higher than can be probed in our snapshot MOJAVE VLBA images, given their image rms levels of ∼0.1 mJy beam⁻¹ and typical jet brightnesses of <100 mJy beam⁻¹ downstream from the core.

In Figure 12 we plot the distribution of parent population sizes for the 10,000 subsamples, which is approximately Gaussian. For our best-fit simulation parameters, typically (3.5 ± 0.3) × 10⁵ parent jets are needed to reproduce the 174 quasar jets in the MOJAVE 1.5JyQC sample. Given the comoving simulated volume of 1334 Gpc³, this implies a parent space density of 261 ± 19 Gpc⁻³, which is comparable to the value of 200 Gpc⁻³ obtained for FR II radio galaxies by Snellen & Best (2001) using the LF of Dunlop & Peacock (1990).

Zoom In Zoom Out Reset image size

A rule of thumb sometimes used in the literature is that for every blazar jet found in a survey with Lorentz factor Γ there are Γ² parent jets (e.g., Mutel 1990; Ghisellini 2000; Berton et al. 2016). This assumption is based on the ratio of solid angle subtended by blazar jets viewed within the critical angle 1/Γ and the full range of jet viewing angle in the parent population, but fails to properly take into account the biases of flux density-limited sampling.

In Figure 13 we plot for our best-fit Monte Carlo simulation the number of parent jets divided by the number of simulated S_ν > 1.5 Jy jets in binned intervals of Lorentz factor between 1 and 50. For Γ ≳ 15, there is a shallow increase in the predicted number of parent jets for each jet found with a particular Lorentz factor in the 1.5JyQC sample, from N ∼ 750 to N ∼ 1300 at Γ = 50. This is much shallower than the rule of thumb Γ² dependence, and is a result of the fact that (i) very high Γ jets are rare in the parent population, and (ii) most of these high Γ jets do not exceed the 1.5 Jy flux density cutoff not only because of their viewing angle, but also their redshift and/or unbeamed luminosity. The large range of possible parent sizes for Γ > 40 reflects the statistical fluctuations associated with selecting from this small cohort of jets in a flux density-limited sample.

Zoom In Zoom Out Reset image size

A different behavior is seen below Γ ≃ 15. These jets are abundant in the parent population, yet most have low unbeamed luminosities and require either substantial Doppler boosting or a low redshift to exceed the flux density cutoff. Every low Γ jet in the 1.5JyQC requires significantly more parent objects, because its maximum possible Doppler boost is only (2Γ)^p. This is the exact opposite of the N ∝ Γ² prediction.

We have carried out a study of the parsec-scale jet kinematics of 409 bright radio-loud AGNs above decl. −30°, based on 15 GHz VLBA data obtained between 1994 August 31 and 2016 December 26. These AGNs have been part of the 2 cm VLBA survey or MOJAVE programs and have ≳0.1 Jy of correlated flux density at 15 GHz. By modeling the jet emission with a series of Gaussians in the interferometric visibility plane, we identified and tracked 1744 individual features in 382 jets over at least five epochs. We fitted their sky trajectories with simple radial and vector motion models, and additionally carried out a constant acceleration fit for 881 features that had 10 or more epochs.

A primary goal of the MOJAVE program is to characterize the jet properties of a well-defined flux density-limited sample in order to better understand the blazar parent population. Using the extensive OVRO and UMRAO single-dish monitoring databases, as well as the MOJAVE VLBA archive, we constructed the MOJAVE 1.5 Jy Quarter Century sample, which consists of all 232 AGNs north of J2000 decl. −30° that are known to have exceeded 1.5 Jy in 15 GHz VLBA correlated flux density between 1994.0 and 2019.0. We carried out Monte Carlo simulations to determine the best-fit parent population parameters that reproduced the redshift, radio luminosity, and apparent velocity distributions of the 174 quasars with z ≥ 0.15 in the 1.5JyQC sample.

We summarize our conclusions as follows:

1. A total of 382 of 409 jets had at least one robust bright feature that could be tracked for five or more epochs. A majority (59%) of the well-sampled jet features showed evidence of accelerated motion at the >3σ level.

2. We examined the distribution of apparent speeds within 26 individual jets that had ten or more robust features, and confirmed that each jet tends to have a characteristic speed that is likely related to the underlying flow. Other than a few fast outliers and some slow pattern speeds, the speeds of features in a jet typically lie within ∼±40% of the characteristic speed.

3. We were able to identify 55 features in 42 jets that had unusually slow pattern speeds (μ < 20 $\mu \mathrm{as}\,{\mathrm{yr}}^{-1}$ and at least 10 times slower than the fastest feature in the jet). We confirm the 43 GHz VLBA results of Jorstad et al. (2017) that the vast majority of these lie within 4 pc (projected) of the core feature, and may represent quasi-stationary standing shocks near the jet base.

4. Only 2.5% of the features we studied had velocity vectors directed inward toward the core. In some cases, these are likely due to brightness variations affecting the fitted centroid position of a large diffuse feature, or a feature on a bent trajectory that is crossing the line of sight. In other cases there may be a misidentification of the true core position. We find that in 16 of the 32 jets with apparent inward motion, the inward-moving feature is the closest feature to the core, and that four BL Lac jets have more than one close-in inward-moving feature.

5. We examined the distribution of maximum apparent jet speed for the AGNs in our full sample and the 1.5JyQC sample, and find that it is peaked at low values, with very few speeds above 30 c. Given the fact that large Doppler-biased jet samples should contain examples of the fastest jets in the parent population, and that our survey has not measured any instantaneous speeds above 50 c, this implies that the parent distribution of jet Lorentz factors is not single-valued, but is weighted toward low Γ, with decreasing numbers of jets up to Γ_max = 50.

6. We find a strong correlation between apparent jet speed and synchrotron peak frequency, with the highest jet speeds being found only in AGNs with low ν_p values. Although a fast jet speed does not guarantee that a jet will be detected at TeV gamma-ray energies, it appears to be a minimum requirement for LSP and ISP AGNs. The exceptions to date are the two very nearby radio galaxies 3C 84 and M87, and the BL Lac TXS 0506+056 that has been associated with a high energy neutrino detection event.

7. Our large grid of Monte Carlo parent population simulations yielded several parameter combinations that could adequately reproduce the flux density, redshift, radio luminosity, and apparent velocity distributions of the 174 quasars in the 1.5JyQC sample. These simulations have an unbeamed LF above 10²⁴ Hz with power-law slope −3.2 ≤ γ ≤ −2.8, and pure luminosity evolution of the form e(z) = (1 + z)^ke^z/η, where $7.5\leqslant k\leqslant 8$ and −0.35 ≤ η ≤ −0.30. The parent jet population has a power-law distribution of Lorentz factors with slope $-1.6\leqslant b\leqslant -1.2$, ranging from Γ = 1.25 to Γ = 50, and a Doppler boosting index p = 2. The best-fit parent population (with b = −1.4, γ = −3.1, k = 8.0, and η = −0.35) has a space density of 261 ± 19 Gpc⁻³, which is consistent with that of FR II radio galaxies. Most of the quasars in the 1.5JyQC have a relatively narrow range of intrinsic (unbeamed) parsec-scale 15 GHz radio luminosity between ∼10^24.5 and $\sim {10}^{26.5}$ W Hz⁻¹.

8. Our best-fit simulation indicates that nearly all of the 1.5JyQC quasar jets are viewed at less than ∼10° from the line of sight, with the distribution peaking at 2°. As previously discussed by Vermeulen & Cohen (1994), Lister & Marscher (1997), and Cohen et al. (2007), the most probable jet viewing angle is ∼0.5 times the critical angle ${\theta }_{\mathrm{cr}}={\sin }^{-1}(1/{\rm{\Gamma }})$ where β_app = Γβ.

9. The Lorentz factor distribution of the 174 bright radio quasars in the flux density-limited 1.5JyQC sample peaks between Γ = 5 and Γ = 15, with a rapid falloff past Γ = 20. The breadth of the Γ distribution indicates that adopting a single value of Γ = 10 for all blazars is not well supported by the observational data. The Doppler factor distribution has a similar shape to the Γ distribution, and peaks at δ ≃ 10, declining rapidly past δ ≃ 30. Both distributions are similar to those inferred from variability Doppler factor estimates using OVRO 15 GHz monitoring data by Liodakis et al. (2018).

10. We find that the oft-cited rule of thumb that for every jet found in a survey with Lorentz factor Γ there are Γ² parent jets is incorrect for flux density-limited blazar samples. Above ${\rm{\Gamma }}\simeq 15$, there is only a shallow increase in the expected number of parent jets per source with Γ, while for lower Lorentz factors, the number of parent jets increases rapidly with decreasing Γ.

The MOJAVE project was supported by NASA-Fermi grants Nos. NNX08AV67G, NNX12A087G, and NNX15AU76G. Y.Y.K. and A.B.P. were supported by the Basic Research Program P-28 of the Presidium of the Russian Academy of Sciences and the government of the Russian Federation (agreement 05.Y09.21.0018). T.S. was supported by the Academy of Finland projects Nos. 274477, 284495, and 312496. T.H. was supported by the Academy of Finland project 317383. The Very Long Baseline Array and the National Radio Astronomy Observatory are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research has used observations with RATAN-600 of the Special Astrophysical Observatory of the Russian Academy of Sciences. This work made use of the Swinburne University of Technology software correlator (Deller et al. 2011), developed as part of the Australian Major National Research Facilities Programme and operated under licence. This research has made use of data from the OVRO 40 m monitoring program Richards et al. (2011), which is supported in part by NASA grants Nos. NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants Nos. AST-0808050 and AST-1109911. This research has made use of data from the University of Michigan Radio Astronomy Observatory which has been supported by the University of Michigan and by a series of grants from the National Science Foundation, most recently No. AST-0607523. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Here we provide comments on individual AGNs supplementing those given in Lister et al. (2013, 2016).

0111+021 (UGC 00773): All five features closest to the core in this nearby BL Lac jet (z = 0.047) have inward or possibly inward-directed motions.

0118−272 (OC −230.4): New Gaussian fitting to the epoch data indicates that feature id = 3 no longer has inward motion at the >3σ level as was reported by Lister et al. (2016).

0256+075 (OD 94.7): The large time gaps between epochs made it impossible to reliably cross-identify any robust features in this quasar.

0300+470 (4C +47.08): Additional new epochs and a reanalysis of the data indicates that the previously reported inward-moving feature (id = 2) in this AGN is not robust. The redshift for this BL Lac is unknown, with the NED value of z = 0.475 being an arbitrary assignment. Shaw et al. (2013b) find 0.37 < z < 1.63 based on an optical spectrum.

0518+211 (RGB J0521+212): The two innermost jet features in this BL Lac object have statistically significant inward motion.

0710+196 (WB92 0711+1940): The jet features in this quasar were too weak (<10 mJy) to identify as robust.

1101+384 (Mrk 421): All three innermost jet features of this nearby BL Lac object show inward motion.

1118+073 (MG1 J112039+0704): The location of the core in this jet is uncertain. We assumed the core to lie at the northeasternmost point in the jet.

1148−001 (4C −00.47): We identified the core as the most compact feature of the jet, with a 2 mas feature (id = 5) being located upstream.

1215+303 (ON 325): All three innermost features of this low redshift BL Lac object (z = 0.131) show small but significant inward motion of approximately 25 $\mu \mathrm{as}\,{\mathrm{yr}}^{-1}$ (0.2 c).

1224−132 (PMN J1226−1328): We were unable to identify any robust jet features in this BL Lac object.

PG 1246+586: None of the jet features in this BL Lac object were sufficiently bright or compact enough to be identified as robust.

1253−055 (3C 279): The VLBA epochs in 2013–2014are affected by the emergence of two very bright features (>10 Jy). The most consistent fits during 2014–2015 required fitting an upstream feature (id = 17) that could be the true core that is only strong enough to be visible during these epochs. The reference "core" position that we use in all of our fits may thus be a strong quasi-stationary feature in the flow.

1300+248 (VIPS 0623): There was no jet feature in this BL Lac object that was sufficiently bright or compact enough to be identified as robust.

PKS 1402+044: The innermost jet feature (id = 5) of this quasar shows statistically significant inward motion.

1458+718 (3C 309.1): A reanalysis of the complex located 23 mas south of the core now indicates no significant inward motion.

PG 1553+113: We were unable to identify any robust jet features in this BL Lac object.

1557+565 (VIPS 0926): The NED redshift of z = 0.3 is not confirmed by Shaw et al. (2013a), who find a lower limit of z > 1.049. The innermost jet feature (id = 4) of this BL Lac object shows statistically significant inward motion.

1656+482 (4C +48.41): The innermost jet feature (id = 4) of this BL Lac object shows statistically significant inward motion.

TXS 1811+062: We were unable to identify any robust jet features in this BL Lac object.

1928+738 (4C +73.18): In 2012 a counter-jet feature with an apparent speed of 0.8 c emerged in this quasar jet.

8C 1944+838: The outermost feature (id = 1) in this jet has a statistically significant inward speed, but may not represent true motion due to the large, diffuse nature of the emission.

1ES 1959+650: The jet structure was too weak and compact at 15 GHz to reliably measure any robust features. Piner et al. (2010) obtained a maximum speed measurement of 0.0322 ± 0.0064 $\mathrm{mas}\,{\mathrm{yr}}^{-1}$ at 43 GHz.

2028+492 (MG4 J202932+4925): We were unable to identify any robust jet features in this BL Lac object.

2234+282 (CTD 135): An unpublished 43 GHz VLBA image by Tao An suggests that core is located in the southwest portion of the jet.

TXS 2308+341: The brightest feature in this jet does not appear to be a stable reference point, with correlated positional changes seen in the positions of downstream features seen at several epochs.

S5 2353+816: No robust jet features could be identified in this BL Lac object because of the large time gap in the data set.