Radio Jet Proper-motion Analysis of Nine Distant Quasars above Redshift 3.5

Up to now, jet kinematic studies of radio quasars have barely reached beyond the redshift range at $z>3.5$. This significantly limits our knowledge of high-redshift jets, which can provide key information for understanding the jet nature and the growth of the black holes in the early Universe. In this paper, we selected 9 radio-loud quasars at $z>3.5$ which display milliarcsec-scale jet morphology. We provided evidence on the source nature by presenting high-resolution very long baseline interferometry (VLBI) images of the sample at 8.4~GHz frequency and making spectral index maps. We also consider Gaia optical positions that are available for 7 out of the 9 quasars, for a better identification of the jet components within the radio structures. We find that 6 sources can be classified as core--jet blazars. The remaining 3 objects are more likely young, jetted radio sources, compact symmetric objects. By including multi-epoch archival VLBI data, we also obtained jet component proper motions of the sample and estimated the jet kinematic and geometric parameters (Doppler factor, Lorentz factor, viewing angle). Our results show that at $z>3.5$, the jet apparent transverse speeds do not exceed 20 times the speed of light ($c$). This is consistent with earlier high-redshift quasar measurements in the literature and the tendency derived from low-redshift blazars that fast jet speeds ($>40\,c$) only occur at low redshifts. The results from this paper contribute to the understanding of the cosmological evolution of radio AGN.


INTRODUCTION
Active galactic nuclei (AGNs) are the most powerful persistent objects in the Universe, luminous across the whole range of the electromagnetic spectrum. They can be found from the local Universe up to very high redshifts. AGNs have been observed throughout almost the entire history of the Universe, providing important information for studying the co-evolution of supermassive black holes (SMBHs) and their host galaxies. The most distant quasars discovered so far are ULAS J134208.10+092838.61 at redshift z = 7.54 (Bañados et al. 2018b) and J0313−1806 at z = 7.64 (Wang et al. 2021). Representing the end era of the cosmic reionization, they offer important clues to explore this episode when the neutral hydrogen gas in the Universe became ionized. More than 300 quasars with z > 5.7 are known to date, most of which were discovered in optical observations (Bañados et al. 2015;Jiang et al. 2016;Shen et al. 2019). In recent years, in addition to traditional near-infrared spectroscopic observations, X-ray and mm-wavelength radio observations have also contributed to the discovery of quasars with z > 6 (e.g., Maiolino et al. 2005;Medvedev et al. 2020). Among these high-redshift AGNs, SMBHs with masses exceeding 10 9 M are found. Although there is an observational bias that the most massive and highly accreting black holes are the easiest to observe due to sensitivity limitations, at least the existence of these SMBHs indicates that they have completed their rapid growth in less than one tenth of the current age of the Universe. The formation and rapid growth of the first SMBHs is one of the mysteries of current cosmology and AGN astrophysics.
Relativistic jets are thought to play an important role in promoting rapid accretion of early AGNs. If a significant fraction of the gravitational energy can be channeled into forming and maintaining relativistic jets, then their black holes can remain at high-rate accretion for a long time (Jolley & Kuncic 2008;Ghisellini et al. 2013). In this sense, the dynamical properties and the lifetime of the high-redshift jets are key to testing this scenario. Although the number of high-z AGNs discovered has been increasing over the years (Jiang et al. 2016;Wang et al. 2019), their radio counterparts are less intensely studied. An important reason is that only ∼ 10% of optically detected quasars have active radio counterparts. The ratio of radio-loud and radio-quiet AGNs seems to be consistent at high redshifts and in the local Universe (Bañados et al. 2015), but it has also been suggested that the fraction of radio-loud AGNs tends to evolve with redshift (Diana et al. 2022). Producing relativistic jets is not only one of the ways of AGN energy release (Blandford et al. 2019), but also has an important impact on the feedback to the environment in the host galaxy (Fabian 2012). To date, the role of jets in the formation and evolution of first SMBHs remains an open question.
The extreme subclass of AGNs are blazars, remarkable for their high bulk Lorentz factors in the jet plasma, the small viewing angles of the jets with respect to the line of sight, and thus the large Doppler-boosting factors. Blazars make up nearly half of the radio-loud AGNs with flux density S 1.4GHz > 100 mJy in the high-redshift population (Sotnikova et al. 2021). The jets of high-redshift blazars are typically short, possibly due to the projection effect caused by the small viewing angle. Perhaps also the environment of the early galaxies is not favorable for the development of largescale jets. In the latter case, the relativistic electrons in the extended jets are heavily depleted via inverse-Compton scattering by the enhanced number density of cosmic microwave background (CMB) photons in the high-redshift Universe (Ghisellini et al. 2014). Also, the interaction between the jet and the dense interstellar medium may cause a large amount of the jet's mechanical energy to be consumed and the growth of the jet to be stalled (An et al. 2022a). Typical radio images of the distant jetted AGNs are characterized by compact flat-spectrum cores or core-jet structures. Very long baseline interferometry (VLBI), with its unique high-resolution capability, provides a direct way to probe the pc-scale compact structure and jet properties of high-z AGNs.
To date, there are only six radio-loud quasars imaged with VLBI at extremely high redshifts (z > 6). Their VLBI morphologies indicate different types of radio sources. They do not show a strong tendency towards highly-beamed sources (i.e. blazars) that would be expected from a selection effect caused by relativistic beaming. Here we list the radio-loud quasars at z > 5 with VLBI imaging available in the literature (Table 1). As one can see, most of the objects are marked as a compact quasar. These sources show a compact structure or minor extension on milliarcsec (mas) angular scale (down to a few tens of pc in linear scale). They typically do not have very high brightness temperatures or show flat radio spectrum. Their VLBI studies suggest that these compact high-redshift quasars could be young jetted AGNs in the early Universe (see references in Table 1). The observing wavelengths of current VLBI facilities are typically at a few cm, and the corresponding rest-frame frequencies of high-redshift quasars are much higher than those of their local counterparts studied at the same observing frequency. The emitted frequency is (1 + z) times the observed frequency. This makes the steep-spectrum jets relatively fainter and could cause the lack of radio quasars observed with prominent extended jets at high redshifts (e.g. Gurvits et al. 2015;Krezinger et al. 2022).
The small sample significantly limits our knowledge of high-redshift jets. It was found that there are two distinct cosmological epochs of SMBH formation: the number density of non-jetted AGNs hosting SMBHs peaks at z ∼ 2−2.5, while the formation of SMBHs in jetted AGNs took place earlier, at z ∼ 4 (Sbarrato et al. 2015;Sbarrato 2021, and references therein). This seems to indicate that black holes with powerful jets grow faster than those of the same mass but without jets. Thus, observational studies of jetted AGNs at z ∼ 4 can bring us crucial information about the rapid growth of the first black holes. Moreover, the sample size of radio-loud AGNs in this cosmic period is much larger than during the cosmic reionization (z 7), facilitating statistical studies.
The kinematics of high-redshift jets provides key information for understanding the jet nature. In particular, VLBI can measure the jet proper motion and independently calculate the Doppler-boosting factor, which can be used to estimate the bulk Lorentz factor and viewing angle of the jet flow. These place strong observational constraints on the modeling of the spectral energy distribution (SED) of the jet (Sbarrato et al. 2022). Current detections of jet proper motions of high-redshift (z 4) AGNs are limited to a few sources (Veres et al. 2010;Frey et al. 2015;Perger et al. 2018;An et al. 2020b;Zhang et al. 2020). All these are blazars, with apparent superluminal speeds distributed in a wide range from 2 c to about 20 c, where c denotes the speed of light. In contrast, the jets of high-redshift galaxies move at mildly relativistic speeds (An et al. 2022a), consistent with the traditional notion of their large viewing angle,  (2022) without strong beaming effects. More proper motion measurements help extend our understanding of the nature and geometry of high-redshift jets.
In the present work, we provide new high-quality VLBI images of nine high-redshift radio-loud quasars. By combining archival multi-frequency, multi-epoch VLBI data as well as optical data, we classify their source nature and determine for the first time their jet proper motions and other radio properties on pc scales. This work increases the size of the proper motion sample of z > 3.5 jets by a factor of two. The subsequent sections of the paper are organized as follows. In Section 2, we describe the method of sample selection, the new Very Long Baseline Array (VLBA) observations, and the data processing. Section 3 presents the new VLBI images which are used as the reference for source classification. Section 4 contains the analysis and discussion of the jet kinematics. Section 5 summarizes the main results and conclusions of this paper. Comments on individual target sources are given in Appendix A, extended tables and figures can be found in Appendix B. Throughout this paper, we use the cosmological parameters derived from a flat Λ Cold Dark Matter (ΛCDM) model (Komatsu et al. 2011) with Ω m = 0.27, Ω Λ = 0.73, and H 0 = 70 km s −1 Mpc −1 .

The Selected High-redshift Sample
Our focus is on the jet kinematics of a sample of high-redshift blazars at z ∼ 4. To measure the proper motion of the jet, at least two epochs of VLBI observations at the same observing frequency are required, preferably with similar (u, v) coverage at each epoch. The time interval between epochs should also be taken into account when measuring jet proper motions at high redshifts. The minimum time gap between epochs can be estimated as t gap = Dmin µcomp , where D min is the minimum distance that can be distinguished between available epochs, usually depending on the observing resolution and frequency. Here µ comp denotes the proper motion of a well-identified component. In the case of VLBI measurements, D min can be measured in mas, µ comp in mas yr −1 , and thus t gap in yr. For high-redshift radio jets, their µ comp could be much smaller than their low-redshift counterparts due to the cosmological time dilation, thus a large t gap is usually needed to obtain reliable measurements. This is one of the reasons why only a few z > 4 jet proper motions have been measured so far.
To select suitable AGNs, we first checked the VLBI archive database for radio sources with 3.5 z 4.5. We used the Astrogeo database 1 to build our sample, which is currently the largest collection of VLBI imaging data available. It accumulates data from a series of astrometric and geodetic VLBI surveys such as the VLBA Calibrator Surveys (VCS, Beasley et al. 2002;Fomalont et al. 2003;Petrov et al. 2005Petrov et al. , 2006Kovalev et al. 2007;Petrov et al. 2008;Petrov 2016), the United States Naval Observatory surveys (USNO; Hunt et al. 2021), and various other projects. The Astrogeo database now provides more than 100,000 VLBI images of more than 17,000 AGNs, mostly observed in snapshot mode at 2.3 and 8.4 GHz. Although the quality of the snapshot images is not very high, we only focus on compact radio structures from bright radio sources with flux densities above 50 mJy, so these images are sufficient for our research purposes. By cross-matching the list with the NASA/IPAC Extragalactic Database 2 (NASA/IPAC Extragalactic Database (NED) 2019) and some large optical spectroscopic projects (e.g., Sloan Digital Sky Survey Data Releases 9 and 12, Ahn et al. 2012;Pâris et al. 2017), over 60 z > 3.5 quasars are found (the full sample will be given in our subsequent study). Since VLBI is sensitive to compact sources with high brightness temperatures, the high-redshift extragalactic sources detected in the flux density-limited VLBI surveys are all radio-loud jetted AGNs. From this parent sample, we then selected a sub-sample for the jet proper motion study, according to the following criteria: • The redshift is 3.5 z 4.5.
• The flux density at 8.4 GHz is S 8.4 > 50 mJy, to enable high dynamic range imaging.
• The source shows resolved jet structure in the existing 8.4-GHz archival images, a premise for being able to measure jet component motion.
• Multiple epochs of VLBI imaging are available at 8.4 GHz. A minimum time gap t gap ≈ 7 yr can be estimated, by assuming a moderate jet component speed 10 c and D min from 3-σ position error. A rough estimate of ≈ 0.2 for the position error is based on the 10% restoring beam size of typical 8.4-GHz VLBA observations. The corresponding µ comp is 0.08 − 0.09 mas yr −1 for the given redshift range.
Based on the above criteria, we selected nine targets whose detailed information can be found in Table 2. There is no other bias in this sample, except for the limitation on the flux density and the declination (Dec. > −40 • ). The latter is because the archival data were mostly obtained with the VLBA located on the northern hemisphere. In particular, our sample selection did not restrict the AGN optical properties. At the time of selecting the present sample, J1606+3124 met our criteria. However, a subsequent study (An et al. 2022a) found that its optical classification and redshift are debatable and need further confirmation. An et al. (2022a) suggest that the source is more likely a high-redshift radio galaxy rather than a quasar.
We have collected images from the Astrogeo archive for a total of 39 epochs for these nine sources, with 3-6 epochs per source and a maximum time span of ∼ 24 yr (which equals approximately 4 yr in the source's rest frame). In snapshot mode, a target is observed in a few scans (observation time segments), each lasting for a few minutes. Scans of multiple sources are interleaved to ensure that each source has good (u, v) coverage. Our selected target sources are routinely monitored in the VLBA calibrator surveys, the (u, v) coverage at 8.4 GHz are similar from epoch to epoch, and the prominent components can be clearly imaged with similar quality, even with limited sensitivity. In order to obtain high-quality images of high-z AGNs for better identification of the radio components and classification of their radio structures, we initiated a new epoch of VLBA observations of the selected 9 sources. The new data are used together with the archival data for jet kinematic studies. Details of the new VLBA observations are presented in the next subsection.

VLBA Observations and Data Reduction
The VLBA observations (project code: BZ064, PI: Y. Zhang) were conducted on 2017 February 5 and 2017 March 19 at 8.4 GHz. Nine out of the ten VLBA antennas participated in the observations (see Table B1. The observations were performed using the DDC (digital downconverter) system with a baseband bandwidth of 128 MHz, divided into 256 spectral channels of 500 kHz each, using four baseband channels (IFs) of left and right circular polarization and 2-bit sampling. This configuration results in a total data rate of 2048 Mbit s −1 . Since the target sources are all compact and bright (S 8.4 > 50 mJy), fringe search can be done using themselves as calibrators, so a standard continuum observation mode was used. A total of 6 h of observation time was divided into two sessions: BZ064A (2 h) and BZ064B (4 h), according to the distribution of the right ascensions of the 9 target sources (see Table 2). During the observations, pointing on sources in the same group were alternated, each obtaining about 0.5 h on-source time. This allowed us to reach a fairly good (u, v) coverage and high image quality. We refer to Table B1 for the details of the observations. After observations, the data were correlated with 2-s integration time using the DiFX correlator (Deller et al. 2011) in Socorro (New Mexico, USA). The correlated data were then downloaded to the China Square Kilometre Array (SKA) Regional Centre prototype ) via the internet for further calibration and imaging. We calibrated the data using the VLBI data processing pipeline developed by our group (An et al. 2022b). The main steps of the procedure are detailed below. We used the US National Radio Astronomy Observatory (NRAO) Astronomical Image Processing System (AIPS) software package (Greisen 2003) to calibrate the amplitude and phase of the visibility data. Our pipeline is written in Python language 3 , using ParselTongue as an interface to convert AIPS tasks into Python scripts (Kettenis et al. 2006), so that most of the operations can be executed automatically.
The pipeline first loaded the data into AIPS and performed a couple of tasks to assist the user in inspecting the data quality. Input parameters (names for fringe finders, bandpass calibrators and target sources, reference antenna, solution intervals, etc.) were determined manually. Then the pipeline automatically conducted the amplitude, phase, and bandpass calibrations, and splitted the visibility data into single-source files. The procedures and parameters used in our experiment are described as follows. The visibility amplitudes were calibrated using the antenna gain curves and system temperatures measured at each station during the observations. The AIPS procedure VLBATECR was then used to correct for atmospheric opacity as well. Next, the bright calibrator sources (J2148+0657 for BZ064A and 3C 273 for BZ064B; Table B1) were used to calibrate the delays and global phase errors of the instruments using the FRING task. This operation calibrated and aligned the delays and phases between different sub-bands. Next, global fringe fitting was performed on all sources and the resulting gain solutions were interpolated and applied to calculate and remove the phase errors. We checked and found that over 98% of the phase solutions were successful. Finally, the antenna-based bandpass functions were solved by using data from the calibrators and applied to correct the target sources' visibility data. At this point, the initial calibration was complete. The calibrated data were exported to an external single-source FITS file by averaging over each sub-band (128 MHz each) and a time interval of 2 s.
After calibration, together with the archival calibrated visibility data obtained from the Astrogeo archive, all singlesource data files exported from AIPS were loaded into the Caltech Difmap software package (Shepherd 1997) for final imaging and model-fitting. The hybrid mapping process consists of several iterations of CLEAN and self-calibration. Final CLEAN images were obtained after a few iterations of phase and amplitude self-calibration, repeated by gradually reducing the solution intervals from a few hours to 1 min. The data from the Astrogeo archive have been calibrated, so we just conducted the imaging and model-fitting following similar steps in Difmap as above.
To parameterize the core and the jet components, we conducted model fitting on the self-calibrated visibility data in Difmap. This utilizes the Levenberg-Marquardt non-linear least-squares minimization technique in the visibility domain (Shepherd 1997). During the processing, we followed the rules of continuity and simplicity, assuming that most components move downstream of the jet and the number of jet components does not change much over the successive epochs. Since we care more about the change in the position of the jet components, most components are fitted with circular Gaussian brightness distribution models. Very few components whose structure was too compact were fitted with a point source model. In a few cases, the jet components showed extended structure and could not be adequately fitted by a single circular Gaussian model. These features could result from diffuse emission or newly ejected components along the jet trajectory. Regarding this, several point-source models which were not treated as physical components were applied to account for the extra flux density. These components are also listed in Table B5. Fomalont (1999) introduced a method to estimate the uncertainty of the fitted model component parameters. This method considers only the statistical error of the image and is therefore closely related to the image quality, while the actual error in the observed data is often higher than the statistical error. For the peak brightness and integrated flux density of each component, we considered an additional 5% calibration uncertainty originating from the measured gain curves and system temperatures. For the position of each component, we found that the Fomalont (1999) method underestimates the actual error of the snapshot observations. Due to the short observing time, the sparse (u, v) coverage of the archival Astrogeo data could result in side lobes. In such conditions, we chose to identify components based on our own long-track VLBA observations, and used the synthesized beam size instead of the fitted component size to calculate the position error, i.e., θ beam SNR , where SNR denotes the signal-to-noise ratio. In such snapshot observations, our conservative estimation of the errors provided a more reasonable assessment of the observed properties of the jet, and similar applications are common in large VLBA surveys (e.g. the MOJAVE survey; Lister et al. 2009Lister et al. , 2016Lister et al. , 2019.

RESULTS
Here we show mas-resolution images of the nine high-z radio-loud quasars and classify their radio structures by combining their radio morphology, radio spectrum, and Gaia optical astrometric data. Figure 1 displays the highquality images obtained from our VLBA observations in 2017. The lowest noise level of the images is ∼ 0.1 mJy beam −1 (Table B2), close to the thermal noise level. The components marked with labels are from the model fitting results from Section 2.2; their parameters can be found in Table B5. Our new VLBA observations are more sensitive than the archival astrometric snapshot data from Astrogeo, thus the models fitted to our own 2017 VLBA data are used as a reference. Certain model components are too weak or too diffuse to be detected in some lower-quality images. Those components are not used for proper motion determination. Since these sources were intentionally chosen to have a resolved structure in VLBI images (see Section 2), all images demonstrate a rich jet structure but present two different types of characteristics: • The other three sources (J0048+0640, J0753+4231, and J1606+3124) show two compact components with com-parable brightness at either end of the radio structure.
The archival Astrogeo images show similar structures, with minor differences because faint components are not always detected in the less sensitive snapshot observations. An example is J1 component in J1939−1002. Basic information on the archival VLBI observations is given in Table B3. A conclusive classification of radio structures from VLBI images alone is not always possible, and for this reason, we searched optical astrometric data of these sources from the Gaia (Gaia Collaboration et al. 2016) Early Data Release 3 (Gaia Collaboration et al. 2021) via the Gaia archive (Gaia collaboration 2020). For quasars, the optical nucleus detected by Gaia corresponds to the accretion disk, with some contribution from the innermost part of the synchrotron-emitting jet (e.g. Plavin et al. 2019). It better represents the central black hole than the radio core, which is the synchrotron self-absorbed jet base at the given frequency. We found Gaia data for 7 of our 9 sources. The remaining two (J1606+3124 and J2102+6015) appear too weak in optical, which is consistent with their identification as a galaxy (see below for a detailed discussion).
In Figure 1, we marked the Gaia positions with crosses whose size represents the 3-sigma positional errors. By comparing optical positions with the mas-scale radio structure of the 7 sources, we found the following: • J0048+0640: the Gaia position lies between the two components, NE and SW. Considering the comparable brightness of NE and SW components and their almost equal distances from the optical nucleus, we can classify J0048+0640 as a compact symmetric object (CSO; Phillips & Mutel 1982).
• J0753+4231: the optical position is close to the northernmost component NE2, so this AGN may be a core-jet source. The brightest VLBI component is NE1, with a separation of about 2 mas downstream along the jet. NE1 may be a moving knot or a standing shock in the jet.
• J1230+1139: the optical position is close to the brightest component C, but with a significant deviation. The situation is similar to that of J0753+4231. Moreover, there is a clear extension from component C toward the optical position, indicating that the radio core is weaker than the bright jet knot.
• J1316+6726: the Gaia position corresponds to the brightest radio component C within the positional uncertainty. This AGN is classified as a core-jet source.
Since the astrometric/geodetic snapshot observations found in the Astrogeo database are carried out in dual frequency bands (either at around 2 and 8 GHz, or at 4 and 8 GHz), the simultaneous dual-frequency observations with the same pointing centers enable us to produce spectral index maps for each target source. This can assist in the classification of their radio structure. In making the spectral index maps, we first selected a set of 2/8-GHz (S/X-band) data with relatively high-quality images for each source, and performed the same hybrid mapping procedure. We set the same image size, pixel size, and restoring beam size for both bands. Model fitting was also performed on the self-calibrated 2-GHz visibility data to identify the optically thin components that were used for calculating the offset between the S/X image brightness peaks, in order to properly align them. The spectral index maps were finally made using our Python script which follows similar calculations as done in the AIPS task COMB. Features with brightness lower than 3 times the rms noise in the 2-GHz maps were omitted to make the spectral index images clearer and more reliable.
The resulting spectral index images are shown in Figure B1 (The spectral index α is defined as S ∝ ν α , where S is the flux density and ν the frequency). J0048+0640 shows a steep spectrum throughout the whole emission structure (mainly NE and SW components), reinforcing its classification as a CSO. The southern component of J0753+4231 has a steep spectrum and the northern component has a flat spectrum, which is consistent with its core-jet classification. In J1230+1139, J1316+6726, J1421−0643, J1445+0958, and J1939−1002, the flattest spectrum (or rising spectrum with positive spectral index) appears at one end of the radio structure, corresponding to the position of the optical nucleus, and the rest of the structure shows a steep spectrum. This information supports that the flat-spectrum components are associated with the radio cores and these sources are core-jet quasars. J1606+3124 lacks a distinct flat-spectrum component, and it is most likely a CSO (see the discussion in An et al. 2022a). Based on its VLBI images, J2102+6015 is a rather peculiar object. Both its eastern and western parts have relatively flat spectra. In the higher-resolution images, the eastern component is resolved into three sub-components . However, none of the VLBI components have detectable proper motion (Zhang et al. 2021). The spectral index map supports that J2102+6015 could be a CSO candidate (see Zhang et al. 2021).  Table B2. The green crosses mark the optical positions detected by Gaia, the bar length are 3 times the Gaia positional errors reported in https://gea.esac.esa.int/archive.
In summary, the mas-scale morphology, the radio spectral index map, and, where Gaia data are available, the cross-match between the optical nucleus and radio components all together suggest the CSO classification of the radio structure in J0048+0640, J1606+3124, and J2102+6015. The other 6 high-redshift quasars are core-jet sources.

DISCUSSION
In this section, we focus on the jet proper motion and discuss the properties of the relativistic jets in our sample of high-redshift jetted AGNs.

Core Brightness Temperature and Doppler Boosting
Extremely high brightness temperatures observed in AGN cores, close to or above the inverse-Compton limit (Kellermann & Pauliny-Toth 1969), are usually considered as due to the beaming effect of relativistic jets pointing close to the observer's line of sight. VLBI observations can be used directly to measure the brightness temperature: where S ν is the flux density of the VLBI component (in Jy), ν the observing frequency (in GHz), and θ comp is the diameter (full width at half-maximum, FWHM) of the circular Gaussian component (in mas). The estimated T B,obs values for the VLBI components in each source are listed in Table B5. The radio core is conventionally defined as the optically-thick section of the jet base in the vicinity of the central SMBH. For radio-loud quasars with a core-jet structure, the core is usually the brightest and most compact component in the VLBI image. Assuming that the magnetic field energy density and the particle energy density are in an equipartition state, the brightness temperature will have a maximum value called equipartition brightness temperature, T B,eq (Readhead 1994). This can be considered as the intrinsic brightness temperature of the relativistic jet. Observed brightness temperatures above the limit (T B,eq ≈ 5 × 10 10 K) are generally considered to be caused by the Doppler boosting effect of beamed jets. The Doppler factor can be estimated from the observed T B,obs values as δ = T B,obs /T B,eq .
We must be aware that the equipartition estimate of the Doppler factor is valid for the optically thick parts of the jet, where the frequency of the T B,obs measurement should be close to the spectral peak of the source. For the highredshift quasars in our sample, the observed 8.4 GHz corresponds to a rest-frame frequency of ∼ 40 GHz, which most likely exceeds the spectral peak frequencies. To estimate the Doppler factors of our sample, we adopted the following approach to obtain appropriate values. Firstly, we used the T B,obs values measured from all the 8-GHz (X-band) epochs to calculate an average brightness temperature for each source. During the averaging calculations, the extreme T B,obs values were omitted, reducing the impact of possible model-fitting biases due to some potentially poor-quality data. The average T B,obs is likely more characteristic to the quiescent states of the AGN core. The resulting average T B,obs for each source can be found among the individual comments on the sources (Appendix A).
Recently, Cheng et al. (2020) estimated the correlation between the observed T B,obs and the frequencies, based on large VLBA surveys of compact radio AGNs at multiple frequencies. They found that the spectral peaks are at around 7 GHz in the rest frame of the sources. Based on the empirical correlation from Cheng et al. (2020), and using the equipartition brightness temperature as the intrinsic brightness temperature that is valid at the spectral peak frequency, we can then estimate the intrinsic brightness temperature T B,int for each of our observed sources, taking the actual rest-frame frequencies into account. The extrapolated T B,int values are in the range of (1.1 − 1.4) × 10 10 K for the sources in our sample.
From the estimates above, we find that all the six core-jet sources in our sample have high core brightness temperatures exceeding the extrapolated T B,int , confirming that they contain highly relativistic jets with Doppler-boosted radio emission. We calculated their Doppler factors that can be found in Table 3. The brightness temperatures of J1316+6726, J1445+0958, and J1939−1002 are high and these values have prominent temporal variations. For J0753+4231, the brightest jet component NE1 has the highest brightness temperature. The three CSOs, J0048+0640, J1606+3124, and J2102+6015, do not have identifiable cores, but the brightness temperatures of their hot spots exceed T B,eq , suggesting that the equipartition assumption probably does not apply in these regions. The brightness temperatures of J0048+0640 hot spots show a significant decrease at epochs 2017 September 18 and 2018 January 18, while a significant increase in the component size occurs. This can be explained by the adiabatic expansion of the hot spots.

Jet Proper Motion
Since the cosmological time dilation is proportional to (1 + z), the jets in the high-redshift sample appear changing slower. Thus reliable component proper motion measurements require VLBI observations spanning a longer time interval in the observer's frame (e.g. Frey et al. 2015;Perger et al. 2018;An et al. 2020b;Zhang et al. 2020).
Using our fitted Gaussian model components, we measured apparent proper motions based on the time evolution of the separation of core and jet components. For components that are too close to the core, their fitted positions may be affected by the finite restoring beam size and perhaps newly ejected features. On the other hand, components too far away from the core ( 10 mas) are usually weak and diffuse, their positional uncertainties are too large. Therefore, we excluded these components from our proper motion determination.
Finally, we derived proper motions for 18 jet features from 9 target sources. The results are shown in Table B4. We fitted the linear proper motion along the RA (µ x ) and Dec (µ y ) directions separately, using a least-squares method, and then calculated the total proper motion as µ 2 x + µ 2 y . For the core-jet sources, we calculated the rate of change of the jet component position relative to the core with time. For CSO sources, we calculated the separation speed of the two opposite hot spots using a particular terminal hot spot as a reference. Assuming that both hot spots advance with the same speed, the advance speed of a particular hot spot is half of the calculated separation speed. Figure B2 demonstrates the jet component trajectories and proper motion fitting results of the VLBI components in the sample. For each source, we selected the most appropriate jet components for proper motion measurements: they have data from the most epochs, are clearly distinguishable in VLBI images, and have the smallest positional errors.
The detailed radio properties of two CSOs, J1606+3124 and J2102+6015, have been presented before (An et al. 2022a;Zhang et al. 2021), so we only briefly summarize the results here. The separation speed between terminal hot spots S and N in J1606+3124 is 0.013 ± 0.002 mas yr −1 (1.60 ± 0.25 c), and between the hot spot S and the inner jet knot C is 0.006 ± 0.002 mas yr −1 (0.74 ± 0.25 c). This leads to a hot spot advance speed of 0.8 c (An et al. 2022a). The separation speed between the eastern and western components in J2102+6015 is 0.023 ± 0.011 mas yr −1 (2.8 ± 1.4 c) (Zhang et al. 2021). In the third CSO in our sample, J0048+0640, we obtained a separation speed of 0.005 ± 0.002 mas yr −1 (0.6 ± 0.2 c), giving a hot spot advance speed of ∼ 0.3 c.
The apparent jet component speeds in the core-jet sources are in the range of 1.4 − 17.5 c, consistent with the known proper motions in low-redshift radio-loud quasars (e.g. Piner et al. 2012;Lister et al. 2019). The maximum speed is close to that of the fastest high-z jet observed before (Zhang et al. 2020), but lower than the maximum value found in low-redshift AGNs. The quasars J1230−1139, J1316+6726, and J1445+0958 contain fast-moving components in the outer part of the jets. These large proper motions could be caused by the projection effect of a large jet bending.

Lorentz Factors and Viewing Angles
High-redshift blazar sources are potentially valuable for studying the cosmological evolution of radio source number density and SMBH accretion. Understanding the distribution of the jet Lorentz factors is fundamental for assessing the number density of AGNs with jets misaligned with respect to the line of sight. These constitute the mostly hidden parent population of highly-beamed jetted sources (blazars). The latter objects are more easily detectable in flux density-limited observations, because of their Doppler-boosted emission. Kinematic measurements of high-redshift AGN jets and estimates of jet Lorentz factors such as presented here are still very rare. For sources at lower redshifts, Lister et al. (2019) found that the distribution of Lorentz factors peaks between Γ = 5 − 15, with a shallow tail reaching Γ ≈ 50. The misaligned parent population of low-Γ jetted sources is larger because those jets require very small viewing angles to be detected as blazars.
The Lorentz factor can be obtained by fitting the broad-band spectral energy distribution (SED) of a blazar (e.g. Boettcher et al. 1997), or directly from VLBI observations. Based on the Doppler factor and the apparent superluminal speed of the jet derived from the VLBI observations, we can estimate the bulk Lorentz factor and viewing angle (see Eqs. B5 and B7 in Ghisellini et al. 1993). The bulk Lorentz factors and jet viewing angles are calculated for five core-jet sources (with the exception of J1316+6726, due to the lack of valid proper motion measurements). The results are presented in Table 3, along with values for other high-z radio quasars taken from the literature. Two sources, J1230−1139 and J1421−0643, show relatively large Lorentz factors, which exceed the typical maximum values (Γ > 20, e.g. Kellermann et al. 2004) based on the β app,max for the studied high-redshift sources. In the relativistic beaming model (see Appendix A in Urry & Padovani 1995), this can happen when the viewing angle is very small (i.e. θ < θ crit , where the critical angle θ crit = arcsin(Γ −1 )). In such cases, a minor change in the viewing angle would greatly increase the jet apparent speed. An alternative possibility could be that the sizes of the VLBI cores are overestimated, which could be caused by the overlap of the emitting components or variability. This will lead to a decrease in the estimated Doppler factor. Both cases above are common in blazars with core-jet structures. These observational evidences, including higher Lorentz factors and smaller viewing angles, support that these sources belong to the blazar class.

High-redshift Quasars in the Apparent Proper Motion-Redshift Diagram
Except for J0906+6930 (Zhang et al. 2017;An et al. 2020b), J1026+2542 , J1430+4204 (Zhang et al. 2020), and J2134−0419 , there are no robust proper motion measurements in other quasar jets at z > 3.5. The main reason is the lack of known sources that are sufficiently bright with prominent well-resolved mas-scale radio jet structure. Also, because of the cosmological time dilation, reliably detecting positional changes in jet components requires decades-long history of VLBI monitoring observations. Thanks to the accumulated astrometric snapshot VLBI data, here we could derive jet proper motions for another 6 high-redshift blazars. This increases the available proper motion sample in the early Universe substantially (Table 3).
An even larger sample of high-redshift jet proper motions could eventually become useful to constrain cosmological model parameters through the apparent proper motion-redshift (µ − z) relation (e.g. Cohen et al. 1988;Vermeulen & Cohen 1994;Kellermann et al. 1999). Earlier studies based on large but lower-redshift source samples (Figure 2)  found that the upper bound on the µ − z relation is consistent with the predictions of the ΛCDM cosmology and a distribution of jet Lorentz factors where the vast majority of the jets have Γ 25. For a given Lorentz factor, the apparent jet component speed cannot exceed √ Γ 2 − 1 ≈ Γ (e.g. Urry & Padovani 1995). To better populate the high-redshift region, we added our own proper motion measurements to the µ − z diagram constructed from literature data (Figure 2). Large VLBI surveys (Britzen et al. 2008;Piner et al. 2012;Lister et al. 2021) at lower redshifts, as well as measurements for individual high-redshift quasars are considered. From our sample studied here, only the highest jet component speeds of each source are plotted. Note that VLBI observations made at different frequencies (ranging from 5 to 15 GHz) are collected in Figure 2. This may result in systematically different apparent proper motion values (e.g. Kellermann et al. 2004). Nevertheless, the general trend in which our new measurements also fit is clearly seen: jets with the fastest apparent proper motion, β app 40, are only found at z 2. Given that our estimated Lorentz factors reach about 40 (Table 3), it is expected that with larger samples of high-redshift quasars, their apparent proper motions could be as high as µ ≈ 0.4 mas yr −1 (Figure 2), without violating the current cosmological paradigm and without requiring extremely high bulk Lorentz factors in the jet. Our study helps populating the high-redshift end of the apparent proper motion-redshift diagram with reliable jet proper motions measured with VLBI.

SUMMARY
In this paper, we reported new 8.4-GHz VLBA observations of 9 high-redshift (z > 3.5) jetted quasars. Based on archival dual-band (usually 2-and 8-GHz) astrometric VLBI data and our high-quality VLBA data taken in 2017, we presented high-resolution radio images and spectral index maps of each source in the sample. Accurate optical positions from Gaia were also considered for 7 out of the 9 sources, in the context of their mas-scale radio structure. By fitting the source visibility data with circular Gaussian brightness distribution model components, and using appropriate component identifications across multiple observing epochs, kinematic properties of the jets were derived.
In our sample, six sources (J0753+4231, J1230−1139, J1316+6726, J1421−0643, J1445+0958, and J1939−1002) were classified as core-jet blazars. These sources have flat-spectrum (or inverted-spectrum) features at the brighter end of the radio structure (the so-called cores), which also closely correspond to the position of the optical nucleus where Gaia measurements are available. The brightness temperatures of their VLBI cores all exceed the estimated equipartition threshold. In the framework of the relativistic beaming model, we also estimated the jet kinematic and geometric parameters (Doppler factor, Lorentz factor, and jet viewing angle) of these core-jet blazars (except for J1316+6726 due to the lack of valid proper motion measurements). The results are presented in Table 3.
Three sources in our sample (J0048+0640, J1606+3124, and J2102+6015) were classified as CSOs or CSO candidates. The spectral index maps of J0048+0640 and J1606+3124 show steep-spectrum emission at both sides of their jet ends. For J2102+6015, its eastern and western features have relatively flat (or inverted) spectra, and none of the VLBI components have detectable proper motion (see also Zhang et al. 2021). The maximum jet proper motion detected in our sample is ∼ 0.15 mas yr −1 , corresponding to an apparent jet speed of ∼ 18 c. The range of our jet proper motions shows good consistency with low-redshift quasars, where large values of β app only appear at low redshifts (e.g. Piner et al. 2012;Lister et al. 2019). Our study substantially increases the sample of high-redshift radio quasars with reliable jet proper motions measured with VLBI. It may serve as an important starting point for accumulating data for future studies of high-redshift AGN jets.  The first VLBI image of this source was made at 5 GHz from a European VLBI Network (EVN) observation in 1996. It showed a double-component source extended along the northeast-southwest direction (Paragi et al. 1999). The two components were separated by ∼ 3.8 mas, consistently with our 8.4-GHz imaging results. The source is marginally resolved at 2.3 GHz and clearly resolved into two components at 8.4 GHz (Figure 1). In producing the spectral index image, we used the midpoints of the 2.3-and 8.4-GHz images for the alignment (the two components can be extracted from the best-quality 2.3-GHz data with model fitting). Both the NE and SW components exhibit flat spectra and high brightness temperatures. The average T B is 7.0 ± 1.3 × 10 10 K for NE and 4.6 ± 3.9 × 10 10 K for SW, with the Gaia nucleus positioned in between. The possibility of gravitational lensing can be ruled out because the two components would violate the preservation of surface brightness (see e.g. Spingola et al. 2019;Casadio et al. 2021). Considering the source compactness, it is very unlikely to be a dual AGN, but this scenario is difficult to rule out with observations at other wavebands because of the insufficient resolution. The most likely explanation is that J0048+0640 is a CSO. No significant proper motion was detected between the two components over a 14-yr time baseline. Recent multi-frequency total flux density measurements of z > 3 AGNs with the RATAN-600 radio telescope show a peaked spectrum for J0048+0640 with ν peak ≈ 4 GHz in the observer's frame (Sotnikova et al. 2021), reinforcing that it is a young CSO source.

A.2. J0753+4231
The source was included in the CSO candidate sample of CSOs observed in the northern sky (COINS) and VLBA Imaging and Polarimetry Survey (VIPS) (Peck & Taylor 2000;Helmboldt et al. 2007;Tremblay et al. 2016). In recent studies of the radio spectra of high-z quasars, this source was also identified as a candidate MHz-peaked spectrum source with ν peak < 1 GHz in the observer's frame (e.g. Mingaliev et al. 2013;Sotnikova et al. 2021). In our study, J0753+4231 shows a double structure in both the 2.3-and 8.4-GHz images. The Gaia position and the spectral index map support that the source is a one-sided core-jet quasar rather than a CSO. The core should be located in the outermost NE2 component. The average brightness temperatures T B of NE1 and NE2 are 8.4 ± 1.7 × 10 10 K and 1.4 ± 0.4 × 10 10 K, respectively. We estimate the maximum proper motion of the jet to be 0.018 ± 0.001 mas yr −1 (1.8 ± 0.1 c) using NE2 (core) as the reference point. In a previous study, Britzen et al. (2008) also detected four components at 5 GHz, and conducted proper motion measurements based on three VLBI epochs, leading to a maximum speed ∼ 0.1 mas yr −1 . Compared to their proper motions, our measured value is much smaller, but it is based on a longer time period (22 yr) and has higher accuracy.

A.3. J1230−1139
The 8.4-GHz image of this source shows a prominent curved jet that bends from the west to the southwest at a projected distance of about 45 pc. Two components, J1 and J2, are detected in the curved southwest jet section, and J3 is located where the jet bends. The brightest component is C, but it is not a radio core. The extension east of C shows a flat spectrum and should host the core. The maximum proper motion is ∼ 0.10 mas yr −1 , corresponding to an apparent transverse speed of ∼ 12 c. The mean T B of this source is 3.4 ± 2.7 × 10 10 K. The relatively lower brightness temperature may be due to the presence of self-absorption in the core, resulting in an underestimate of the flux density and an overestimate of its size. The estimated viewing angle of the innermost jet is around 10 • , derived from the Doppler factor and maximum apparent jet speed, classifying it as a blazar. The source also shows significant variability, consistent with its blazar identification. It shares similar high-resolution jet properties with other high-z blazars (e.g. Zhang et al. 2020;An et al. 2020b).

A.4. J1316+6726
The redshift of this source is from photometric measurements, but with a very high probability (∼ 95%, Richards et al. 2009). Since the jet is well-separated from the core, we included it in our sample and tried to check its high-z nature by determining its jet proper motions following the method introduced by An et al. (2020a). In the VLBI images of this source, a weak and extended jet can be found at around 8 mas southeast of the bright core component, with a sharp bending at a projected distance of 60 pc from the core. The jet is only marginally detected at 3 available epochs, and we did not detect a significant (i.e. ≥ 3σ) outward motion due to the significant positional errors and the short time period. The mean T B of the core component is 13.3 ± 5.5 × 10 10 K, indicating strong non-thermal emission. The apparent inward motion needs to be verified by future observations. If this proper motion is caused by the motion of a newly-formed jet component in the core, then the new jet speed along the jet direction is 14.2 c, which does not exceed the expected proper motion upper limit for a high-z jet (see an example of the opposite situation in An et al. 2020a), and thus its photometric redshift can be considered plausible.

A.5. J1421−0643
This source shows a typical core-jet radio structure. The mean T B is 4.2 ± 2.1 × 10 10 K. This source is a prominent high-z blazar. Previous studies have discovered large kpc-scale radio and X-ray jets extending towards the northeast of this nucleus (Worrall et al. 2020). Our VLBI images reveal the pc-scale jet structure, showing three jet components moving away from the core radially toward the northeast direction. The position angle of the pc-scale jet is consistent with that of the kpc-scale jet (about 30 • ). The estimated maximum jet proper motion is ∼ 0.15 mas yr −1 (∼ 16 c).

A.6. J1445+0958
This object, also known as OQ 172, has a rich jet structure starting toward the west and turning the south. It was once identified as a GHz-peaked spectrum (GPS) galaxy due to the less variable radio and optical emission, the non-flat radio spectrum and the lack of pc-scale structural evolution (see Punsly et al. 2015, and the references therein). Previous multi-band VLBI observations revealed a clockwise jet trajectory from high-frequency to low-frequency images, which was attributed to the interaction between the relativistic jet and the narrow-line region medium (Liu et al. 2017). In our image at 8.4 GHz, the emission of the jet that bends from west to south was modeled and analyzed. Since the jet is becoming weak and diffuse in the southern tail, barely detected in the less sensitive archival snapshot observations, we only fit the brighter jet section in the northern part (Figure 1). From seven epochs spanning 21 yr, we were able to obtain jet proper motion in J1445+0958 for the first time. Among the three jet components, J1 stands out as the fastest one, with a proper motion of 0.13 ± 0.01 mas yr −1 . The mean T B of this source is 15.5±6.4×10 10 K, suggesting a highly beamed jet. From the spectral index map, the high brightness temperature, and the Gaia optical position, the source could be a blazar. Although its GPS-type radio spectrum is inconsistent with the typically flat spectrum of low-redshift blazars, it is not unprecedented in high-redshift quasars (e.g. Sotnikova et al. 2021).

A.7. J1606+3124
This source was previously identified as a flat-spectrum radio quasar (Healey et al. 2007;Torniainen et al. 2007;Coppejans et al. 2016). However, further radio spectral studies based on simultaneous multi-frequency observations suggested that the source is a GPS radio galaxy in the early Universe (RATAN-600, Mingaliev et al. 2012;Sotnikova et al. 2019). In our study, by analyzing the spectral index map and the archival radio spectra, we conclude that the source could be the a high-redshift CSO source. Additional results and a detailed discussion of the CSO identification have been presented elsewhere (An et al. 2022a).

A.8. J1939−1002
The VLBI morphology resembles a core-jet radio source. The two compact jet components lie close to the bright core, and a distant jet component (J1) is in the northeastern direction at a distance of about 22 mas. This distant feature shows a slight change in the position angle. It could possibly be explained by the rise of the Doppler-boosted region caused by the helical jet path, which is usually seen in flat-spectrum radio quasars (FSRQs) (e.g. Alberdi et al. 2000;Hong et al. 2004). However, the possibility of jet interaction with the surrounding interstellar medium (ISM) cannot be ruled out either (e.g. Gómez et al. 2000;An et al. 2020b). The mean T B of this source is 21.4 ± 6.7 × 10 10 K, which well exceeds the equipartition value and indicates strong Doppler boosting with δ ≥ 4.
Since the jet component J1 is weak and diffuse at most available epochs, we just measured proper motions of J2 and J3. The resulting jet proper motion is −0.014 ± 0.001 mas yr −1 for J2 and 0.005 ± 0.001 mas yr −1 for J3. The apparent inward jet motion is physically meaningless. Despite position errors, the curved jet motion across the line of sight and the newly-emerging features that move beyond the time separations among the available epochs could be the possible reasons for the negative proper motion of J2 (see e.g. Lister et al. 2013Lister et al. , 2016.

A.9. J2102+6015
Previous studies of this source claim it as an FSRQ with moderate Doppler boosting effects (Coppejans et al. 2016), but recent works prefer J2102+6015 to be a GPS radio source (e.g. Coppejans et al. 2017;Frey et al. 2018). We find further interesting characteristics of this source. Considering the high-resolution images, the spectral index behavior, and the component separation speeds, we believe this high-z AGN is most likely a CSO. From the deep imaging and model-fitting results, we found that neither the E nor the W features are unresolved. To further constrain the source nature, we collected more high-resolution VLBI data from further epochs, and used somewhat more conservative positional error estimates (i.e. one-tenth of the restoring beam size) to conduct component identification and proper motion estimates. A more robust proper motion estimate of 0.023 ± 0.011 mas yr −1 was obtained. The details are presented in a separate paper (Zhang et al. 2021).   Table B1 lists the basic observing information for the proposed VLBA session BZ064. Table B2 presents the parameters of the VLBI CLEAN images for each target source in our VLBA observation BZ064. Table B3 lists the 8-GHz VLBI imaging parameters of the target sources based on their archival VLBA experiments. Table B4 is showing the jet proper motion parameters of the major jet components within each source in our sample, which were estimated based on our multi-epoch 8.4 GHz VLBA observations. Table B5 catalogued the parameters of the brightness distribution models (Gaussian or point-source models) that are used to fit the target sources in the visibility domain. Figure B1 shows the spectral index maps of the sources in our sample, made from their simultaneous 2.3 and 8,4 GHz VLBA images. Figure B2 demonstrates how the jet components move with respect to their "cores" (optically thick radio cores or referencing components) in our sample and exhibits the fitted results of their radial motions and position angle changes.  Note-Columns: (1) J2000 source name; (2) date of the observation; (3) peak intensity of the image; (4)-(5) major and minor axes from the FWHM size of the elliptical Gaussian restoring beam; (6) position angle of the beam major axis, measured from north to east; (7) the rms noise of the post-fit image.  (3) project code; (4) observing frequency; (5) peak intensity of the image; (6) rms noise, (7)-(8) major and minor FWHM axes of the elliptical Gaussian restoring beam; (9) position angle of the beam major axis, measured from north to east; (10) observing bandwidth of each IF × number of IFs.  Note-Columns: (1) J2000 source name; (2) component identifier; (3)-(6) fitted proper motion in mas yr −1 : µ rad -the jet radial motion; µx, µy -proper motion components along the right ascension and declination axes, respectively; µtot: total proper motion calculated from µx and µy; (7)-(8) apparent transverse speed of jet radial (β rad ) and total (βtot) motions, respectively, in units of c; (9) change of the jet position angle, in • yr −1 .