Proper motion of the radio jets in two blazars at redshift above 3

There is still a limited number of high-redshift ( 𝑧 > 3) active galactic nuclei (AGN) whose jet kinematics have been studied with very long baseline interferometry (VLBI). Without a dedicated proper motion survey, regularly conducted astrometric VLBI observations of bright radio-emitting AGN with sensitive arrays can be utilized to follow changes in the jets, by means of high-resolution imaging and brightness distribution modeling. Here we present a first-time VLBI jet kinematic study of NVSS J080518 + 614423 ( 𝑧 = 3 . 033) and NVSS J165844 − 073918 ( 𝑧 = 3 . 742), two flat-spectrum radio quasars that display milliarcsecond-scale jet morphology. Archival astrometric observations carried out mainly with the Very Long Baseline Array, supplemented by recent data taken with the European VLBI Network, allowed us to monitor changes in their radio structure in the 7 . 6 − 8 . 6 GHz frequency band, covering almost two decades. By identifying individual jet components at each epoch, we were able to determine the apparent proper motion for multiple features in both sources. Apparent superluminal motions range between ( 1 − 14 ) 𝑐 , and are found to be consistent with studies of other high-redshift AGN targets. Using the physical parameters derived from the brightness distribution modeling, we estimate the Doppler-boosting factors ( 𝛿 ≈ 11 . 2 and 𝛿 ≈ 2 . 7), the Lorentz factors ( Γ ≈ 7 . 4 and Γ ≈ 36 . 6) and the jet viewing angles ( 𝜃 ≈ 4 . ◦ 4 and 𝜃 ≈ 8 . ◦ 0), for NVSS J080518 + 614423 and NVSS J165844 − 073918, respectively. The data revealed a stationary jet component with negligible apparent proper motion in NVSS J165844 − 073918.


INTRODUCTION
Active galactic nuclei (AGN) can be found and studied through almost the entire history of the Universe.These are the most powerful non-transient objects we currently know and can be observed even at extremely high redshifts ( > 7, e.g.Bañados et al. 2018;Wang et al. 2021;Castellano et al. 2023;Larson et al. 2023).AGN provide essential information about the behaviour of supermassive black holes (SMBHs) and the connection to their host galaxies.They are powered by accretion onto their central SMBHs, thus releasing a large amount of energy across the whole electromagnetic spectrum.The most massive high-redshift SMBHs known reach masses of ∼ 10 10 M ⊙ (e.g.Wang et al. 2013;Shen et al. 2019;Sbarrato 2021).The rapid growth of SMBHs in the early history of the Universe is one of the most intriguing fields of research in modern astrophysics and cosmology ★ E-mail: krezinger.mate@csfk.org(see e.g.Fan et al. 2022 and references therein).The early accretion of AGN can be explained by maintaining powerful relativistic jets that carry off a significant fraction of the released gravitational energy (Jolley & Kuncic 2008;Ghisellini et al. 2013).This way a black hole can keep a high accretion rate for a long time, which means faster mass growth.It is also possible for the jets to trigger an infall of galactic matter in the vicinity of the black hole (Fabian 2012).Such feedback has already been observed in local active galaxies (e.g Silk & Rees 1998;Cresci & Maiolino 2018;Yuan et al. 2018).High-redshift jetted AGN more likely harbour black holes with mass exceeding 10 9 M ⊙ (Sbarrato et al. 2022), suggesting that jets indeed play an important role in the fast formation of the first SMBHs.They also contribute to the AGN-host galaxy feedback (Fabian 2012).
The number of the AGN found at high redshifts has been increasing steadily over the years (e.g.Jiang et al. 2016;Wang et al. 2019;Onken et al. 2022;DESI Collaboration et al. 2023;Yue et al. 2023).Only a small fraction (∼ 10%) of these quasars are radio-loud, i.e., have powerful relativistic jets producing synchrotron radio emission.Perger et al. (2017) showed that above  ≳ 4, there are only about 300 AGN with radio emission above a typical sensitivity threshold of modern radio telescopes.The majority of them are weak radio sources with mJy-level flux densities at GHz frequencies.At high redshifts, radio-loud AGN are mainly represented by blazars (e.g.Sbarrato et al. 2022).
Blazars are a subclass of AGN in which the jet viewing angle is close ( ≲ 10 • ) to the line of sight (Urry & Padovani 1995).Due to relativistic beaming effects, they tend to have Doppler-boosted radio emission and extremely high brightness temperatures.These objects are often luminous X-ray and/or -ray sources, but less frequently found with increasing redshift (e.g.Sbarrato et al. 2013;Ghisellini & Tavecchio 2015).At high redshifts and with high angular resolution, blazars are usually characterised by a compact, flat-spectrum core or a core-dominated, moderately extended jet structure (e.g.Coppejans et al. 2016;Cao et al. 2017;Krezinger et al. 2022).Apart from the jet base (the core), jet components have steep spectrum, i.e. their brightness rapidly decreases as the frequency increases.The restframe frequency ( 0 ), where the radiation is emitted, is a function of the observed frequency ( obs ) and the redshift as  0 = (1 + )  obs .Thus, at a given  obs , the higher the redshift, the fainter is the steepspectrum jet (e.g.Gurvits 2000;Gurvits et al. 2015).This effect causes the steep-spectrum jets with extended radio emission to remain largely undetected, therefore it is the cores of flat-spectrum radio quasars (FSRQs, Healey et al. 2007) with strong, compact, enhanced synchrotron emission which are typically observed at high redshifts.Milliarcsecond (mas) scale jets of high-redshift blazars tend to appear shorter than jets in the low-redshift Universe.Apart from the observational effect explained above, this is possibly also because of an intrinsic phenomenon: the dense circumnuclear matter that hinders the development of large-scale jets, as only the most powerful ones are able to penetrate it (Ghisellini & Sbarrato 2016;An et al. 2022).Large-scale emission structures are also dimmer due to the interaction of their electrons with photons of the cosmic microwave background (Ghisellini et al. 2015).To detect the radio emission coming from the innermost regions of these distant sources, high angular resolution and high sensitivity are needed.The technique of very long baseline interferometry (VLBI) is capable of achieving mas-scale angular resolution which corresponds to pc scales in linear resolution.
To study the kinematics of the jets, and to acquire key information for understanding the nature of their function, VLBI is essential.With long-term VLBI monitoring, the proper motion of the jet components can be directly measured.From the derived apparent brightness temperature values, we can calculate the properties of the relativistically enhanced plasma, i.e. the Doppler-boosting factor, the bulk Lorentz factor, and the jet inclination angle.This information can be used to further refine the modeling of the spectral energy distribution (SED, e.g.Sbarrato et al. 2022).Existing high-redshift jet proper motion studies are still scarce and usually only deal with a few sources at a time (e.g.Veres et al. 2010;Frey et al. 2015;Perger et al. 2018;An et al. 2020;Zhang et al. 2020Zhang et al. , 2022)).Jet component proper motion is measurable for core-jet type sources, and these are relatively rare at high redshifts.Without a dedicated long-term VLBI observing programme for determining the proper motion of high-redshift jets, a straightforward way is to look for frequently observed sources with core-jet characteristics suitable for component identifications across the epochs.
While only about 20 sources with jet proper motion analysis exist at  > 3, they provide important constraints on early SMBH growth and relativistic jet properties in the early Universe.Even incremental increases in sample size can improve the statistical constraints on the physics of high-redshift jets.To expand this limited sample in order to work towards improving the statistical constraints, we specifically targeted two FSRQs, NVSS J080518+614423 (J0805+6144 in short,  = 3.033, Sowards-Emmerd et al. 2005) and NVSS J165844−073918 (J1658−0739,  = 3.742, Healey et al. 2008) because they are among the brightest high-redshift blazars accessible to VLBI imaging.Their flat radio spectra and compact mas-scale core-jet structures made them suitable for multi-epoch model-fitting analysis to measure jet kinematics.Both can be found among the targets of geodetic and astrometric VLBI observations (Nothnagel et al. 2017), but they have not been the subjects of dedicated VLBI monitoring.
J0805+6144, with coordinates in the 3rd realization of the International Celestial Reference Frame (ICRF3, Charlot et al. 2020) right ascension  VLBI = 08 h 05 m 18. s 179562 and declination  VLBI = 61 • 44 ′ 23.′′ 70056, is also known as a -ray source in the Fermi Large Area Telescope (LAT) catalogues (Abdo et al. 2010;Acero et al. 2015).Significant -ray variability is found by Nolan et al. (2012), andLi et al. (2018) identified two flaring events, one in 2009 January, showing intra-day -ray variability, and another one in 2010 January.Contrary to the -ray variability of the object, the Swift X-ray Telescope (XRT) measured constant X-ray flux during the first ten years of Fermi-LAT operation, from 2008 August to 2018 August (Sahakyan et al. 2020).Several Doppler-factor estimates can be found in the literature, based on multiwavelength SED fitting.Chen ( 2018) derived  = 6.4,also giving an estimate on the central black hole mass 10 9.1 M ⊙ .Later, Sahakyan et al. (2020) found  = 14.0 ± 0.6, while Tan et al. (2020) obtained a value higher by a factor of ∼ 2,  = 28.8.Paliya et al. (2020) estimated  = 22.6, while also giving a Lorentz-factor Γ = 14.In contrast, Sahakyan et al. (2020) found Γ = 1.13 ± 0.10, however, based only on a single X-ray observation.J0805+6144 is bright enough in optical to be detected by the Gaia astrometric space telescope (Gaia Collaboration et al. 2016Collaboration et al. , 2021)).
J1658−0739 (ICRF3 position:  VLBI = 16 h 58 m 44.s 061966,  VLBI = −07 • 39 ′ 17. ′′ 6945) is also a regularly observed radio source.Its NVSS 1.4-GHz flux density is (0.83±0.03)Jy.Its FSRQ classification was confirmed by Healey et al. (2007).The source is considered as a -ray-quiet blazar (Paliya et al. 2017).At high energies, it has only a single X-ray detection by the Swift-XRT instrument (Paliya et al. 2017).With SED fitting, the Doppler factor and the bulk Lorentz factor are estimated as  = 15.7 and Γ = 10, respectively.They found that the jet has an inclination angle of 3. • 0 and the black hole mass is estimated to be 10 9.48 M ⊙ .In optical, the source is detected by Gaia.J1658−0739 was imaged in the VSOP (VLBI Space Observatory Programme) Prelaunch Survey with the VLBA at 5 GHz (Fomalont et al. 2000).The image shows a bright core with jet extending to the northeastern direction.However, later the source was not observed on ground-space interferometer baselines in the VSOP AGN Survey (Hirabayashi et al. 2000).
In this work, we present jet proper motion studies of the two high-redshift blazars, J0805+6144 and J1658−0739, using VLBI measurements spanning almost two decades, based on data from new and archival observations.Such analysis is carried out for the first time for these two radio quasars, increasing the size of the  > 3 proper motion sample by ∼ 10 per cent.Combining the multi-epoch VLBI observations with the available radio spectral information and the Gaia Data Release 3 (DR3) positions, we investigate the nature of the sources by determining their jet proper motions and constraining key physical parameters such as Doppler and Lorentz factors.We aim to add new data points to the sparsely sampled  > 3 region of the apparent proper motion-redshift relation.The results will help refine our understanding of the typical velocities and orientations of the most relativistic jets during the era of early black hole growth.Throughout this paper, we assume a standard ΛCDM cosmological model with Ω m = 0.3, Ω Λ = 0.7, and  0 = 70 km s −1 Mpc −1 .To determine the projected linear sizes and luminosity distances, we used the cosmology calculator of Wright (2006).

NEW VLBI OBSERVATIONS AND ARCHIVAL DATA
J0805+6144 and J1658−0739 were observed as part of a sample of  > 3 blazars with the European VLBI Network (EVN) under the project code ET036 (PI: O. Titov).Altogether, three epochs of observations were made, on 2018 June 4 (segment ET036A), 2018 October 27 (ET036B), and 2019 February 25 (ET036C).The primary purpose of the project was astrometric in nature, to investigate the positional stability of prominent high-redshift radio sources (Titov et al. 2023).However, source brightness distribution maps were also made from the interferometer visibility data.The astrometric-style dual-band snapshot observations were performed at 2.3 GHz (S band) and 8.6 GHz (X band).The data were recorded in right circular polarization and included a total of 16 intermediate frequency channels (IFs), 6 in S band and 10 in X band, each 8-MHz wide.The total bandwidth was 128 MHz and the data rate was 512 Mbps.Table 1 contains the details of the EVN observations, including the dates, frequencies, target sources, and the participating telescopes in each project segment.The observed data were processed at the Joint Institute for VLBI European Research Infrastructure Consortium (JIVE, Dwingeloo, The Netherlands) with the SFXC software correlator (Keimpema et al. 2015) with 0.5 s integration time and 62.5 kHz spectral resolution.
Both J0805+6144 and J1658−0739 are frequently observed blazars and have VLBI data available in the archives covering almost two decades.We collected the archival calibrated 7.6−8.6GHz (X band) VLBI visibility data from the Astrogeo 1 database.These observations were usually performed in the framework of VLBI calibrator surveys and astrometric/geodetic VLBI experiments (Fomalont et al. 2003;Kovalev et al. 2007;Gordon et al. 2016;Schinzel et al. 2017;Hunt et al. 2021).The sources were mostly observed with the VLBA, occasionally supplemented with other radio telescopes, in snapshot mode at S/X bands.Such observations are usually scheduled with few-minute scans on the targets, rapidly alternating between objects seen in different directions (e.g.Sovers et al. 1998).For imaging purposes, snapshots widely spaced in time could provide good (, ) sampling.The time coverage of archival VLBI data is not uniform, the observations are more frequent after 2017.Table 2 summarises the archival calibrated visibility data we downloaded and imaged.Even taking the cosmological time dilation into account, the ∼ 20 years in the observer's frame still translate to 4 years in the restframe of J0805+6144 and J1658−0739, giving sufficiently long time to detect the apparent proper motion of jet components.

EVN observations
We calibrated the EVN data using the U.S National Radio Astronomy Observatory (NRAO) Astronomical Image Processing System (AIPS) software package (Greisen 2003), following a standard procedure (e.g.Diamond 1995).The interferometric visibility amplitudes were calibrated using the antenna gain curves and the system temperatures were measured at the telescopes.For some antennas, only nominal system temperature values were available (ET036A: Zc, T6; ET036B: Zc, Sv, Ur, Wz; ET036C: Zc, Sv, Ur, Wn; for station codes, see Table 1).Ionospheric delays were obtained from global navigation satellite systems measurements, and the phases were corrected using the measured Earth orientation parameters.As a simple bandpass correction, we flagged the first and last 4 channels in each IF.Short 1-min data scans were chosen to solve for instrumental phases and delays.Global fringe-fitting (Schwab & Cotton 1983) was performed on the target sources and the bright fringe-finder sources were also scheduled in the experiments.Finally, the calibrated data were averaged in frequency.
The integration times spent on the sources in each observing session were markedly different (Table 1), which resulted in poor image quality in 2 out of the 3 observing epochs.Since the epochs were close in time, only spanning a few months in the rest frame of the sources, it was feasible to combine the calibrated visibility data in AIPS.By doing so, it was possible to reduce the noise level in the final images, while keeping all data.We associated the date to the combined observations when the on-source time was much longer than at the other two epochs.It is 2018 Jul 4 for J0805+6144 and 2019 Feb 25 for J1658−0739.
The combined calibrated data of the target sources were exported from AIPS for further work in Difmap (Shepherd et al. 1994).There we performed standard hybrid mapping to produce the images of the sources, which includes iterations of the CLEAN algorithm (Högbom 1974) and phase-only self-calibration (Pearson & Readhead 1984) followed by a few rounds of amplitude and phase self-calibration.Gaussian brightness distribution model components were fitted to the self-calibrated visibility data (Pearson 1995), to quantitatively characterize sizes and flux densities of the core and jet components.The core of each source was fitted with elliptical Gaussian model components, while the jet components with circular Gaussians.Elliptical Gaussians fitted to the core provide directional information on the innermost, barely resolved section of the jet.This way we trace possible changes in the position angle of the major axis.

Archival data
After downloading all the available calibrated and self-calibrated X-band visibility data for J0805+6144 and J1658−0739 from the Astrogeo archive, we made images and modelfits by following the procedure as described above.Repeating self-calibration on already self-calibrated data cannot cause any bias.In turn, it may slightly improve the quality with respect to the automatically-generated Astrogeo images, especially in the case of a complex source structure and if outlier visiblity points are flagged manually.We note that, similarly to what is done with our EVN data, many of the archival Astrogeo visibility data sets are combined, containing observations made by different subarrays and at epochs typically separated by hours to days.Two imaging methods were used to ensure consistency in the resulting images and model parameters.First, the data were processed using an automatic Difmap imaging script2 , to quickly obtain an initial image of the two sources at each epoch.Then we also imaged data from all the epochs manually in Difmap.Self-calibrated  visibility data obtained after hybrid mapping served as a basis for model fitting using Gaussian brightness distribution components.Similarly, the core components were fitted with elliptical Gaussians, except for a few cases where a circular Gaussian component provided a better fit.The jet components further away from the core were fitted with circular Gaussians.At most of the epochs, two jet components were found, mainly the two closest ones to the core.At a few epochs, when the array configuration and higher sensitivity allowed imaging of some extended emission, outer jet components also appeared.For epochs close to each other in time, common starting models were used to ensure the consistency of the results.

Modelfit uncertainties
The errors of the fitted model component parameters (except for the separation parameter) were estimated following Fomalont (1999).
The method considers the statistical error in the image.We applied an additional 10 per cent absolute amplitude calibration uncertainty for the fitted flux density of each component (e.g.Middelberg et al. 2011;Pushkarev & Kovalev 2012;Mooley et al. 2018;Punsly et al. 2021).The uncertainties of the component separation were found to be underestimated (with up to ∼ 3 orders of magnitude smaller errors compared to the separation values themselves) if estimated following Fomalont (1999).To tackle this problem, we considered 20 per cent of the full width at half-maximum (FWHM) of the synthesised beam along the component position angle (e.g.Lister et al. 2009Lister et al. , 2013;;Punsly et al. 2021) as the separation uncertainty.For snapshot observations typically with short observing time and sparse (, ) coverage, this conservative estimate gives more reasonable results.Table 3 contains the modelfit results with their calculated uncertainties.2).For both sources, the core-jet structures extending up to about 10 mas are typical for blazars, although this is not unique for blazars, and has been observed in various quasars, including those at high redshifts (e.g.Frey et al. 2010;Coppejans et al. 2016;Krezinger et al. 2022).

Radio morphology of the sources
The Gaia DR3 optical coordinates are marked in the radio images with red crosses whose size represents the 1 positional errors.Broad optical emission lines used for the measurement of the quasar redshift are known to be formed in the so-called broad-line region  (BLR) in the vicinity of the central black hole (Peterson 2006).If the optical emission is dominated by the non-thermal emission from the relativistic jet, the thermal optical emission from the BLR would be faded out and the broad emission lines are not detected.This is a common situation with many blazars whose redshift is unknown so far in spite of their high brightness in the optical wavelengths.Therefore the detection of the broad emission lines in the optical spectra of both quasars favours that the thermal origin of the optical emission and, correspondingly, the Gaia position better represents the location of the central black hole than the radio core (Kovalev et al. 2017;Plavin et al. 2019).The latter is in fact the synchrotron self-absorbed base of the jet with   = 1 optical depth at the given frequency .Because absolute astrometric information is lost after fringe-fitting, we associated the position of the VLBI brightness peaks with the X-band ICRF3 position of the respective quasars.
The VLBI astrometric position generally reflects the location of the radio brightness peak (e.g.Fey et al. 1997).However, according to Porcas (2009), there could be a ∼ 0.2-mas level offset between the phase-referenced and group-delay positions along the jet direction, caused by opacity effects on the emission at the jet base.Moreovoer, its actual value may be affected by the variable core-shift effect.Up to a few mas offsets between the radio and optical AGN positions are found to statistically coincide with the VLBI jet direction (Petrov et al. 2019).Indeed, in both sources, the Gaia optical positions are offset from the radio brightness peak that indicates the location of the core within ∼ 1 mas, apparently upstream along the jet (Figs.1-2).
The 8.6-GHz image of J0805+6144 (Fig. 1) shows the core and two jet components (J1, J2) close to it, within ∼ 3 mas.The innermost jet starts pointing to the south, then seems to turn eastwards at ∼ 1 mas from the core, and continues in the southeastern direction.Further downstream, there is a third, diffuse jet component (J3) at ∼ 8 mas separation.This feature is heavily resolved, or even remains undetected at most epochs where the imaging sensitivity is insufficient.
The image of J1658−0739 (Fig. 2) shows a core and four jet components.The two inner components (J1, J2) are within ∼ 2 mas from the core.The jet points to the northeast, then gradually bends northwards.

Jet proper motion
To estimate the apparent proper motion of jet components, the first task was to identify the same features across the epochs spanning about two decades of observations.As observations were made with a wide range of imaging sensitivity and varying angular resolution, not all epochs provided the same number of detected jet components.The apparent proper motions of securely identified jet components were calculated based on the relative positions of the fitted circular Gaussian brightness distribution models with respect to the core.The measured separations and their corresponding errors do not allow us to estimate reliably any higher time derivative of the separation than the first one.Therefore the angular separation vs. time was fitted with linear function, and the apparent proper motion is taken as the slope of this line.Figures 1 and 2 show the core-jet component separations and position angles as a function of time, along with the best-fit linear trends for both J0805+6144 and J1658−0739.For the VLBI images shown as examples, we choose observing epochs where all fitted jet features can be seen.In Table 4, we collect the fitted kinematic parameters, as well as geometric and physical parameters of the jets calculated in Sect.4.4.Using the apparent proper motions, we can determine whether the jet is superluminal.AGN with relativistic jets pointing close to the line of sight often show superluminal motion which is an apparent relativistic effect (Rees 1966).
For J0805+6144, the apparent speeds of the innermost two jet components, J1 ( J0805 J1 = 4.5 ± 0.9) and J2 ( J0805 J2 = 6.3 ± 1.0), are superluminal and consistent with each other within the uncertainties.The trajectory of J1 seems to follow that of J2, its position angle changes from ∼ 190 • to ∼ 160 • while travelling up to ∼ 2 mas projected distance from the core (Fig. 1).From that point on, the bending jet becomes remarkably straight up to ∼ 8 mas separation.The diffuse outer component J3 was absent in the first few years and could be reliably detected only at four later epochs.The time coverage and sampling for J3 is not sufficient for a reliable proper motion determination, therefore we give an upper limit to its apparent speed in Table 4.
In the case of J1658−0739 (Fig. 2), the four jet components follow a bent trajectory, starting in the northeastern direction, and then turning northwards.The probably newly-emerging innermost components J1 and J2 were not detected at the first epoch in 2002 due to the limited angular resolution.The proper motions are apparently superluminal, except for J2 ( J1658 J2 = 0.9 ± 1.6) which is consistent with being stationary (Table 4).The apparent advance speed of the outer J3 and J4 components is about 14.It should be noted that these outer components are in the more extended parts of the jet.It is harder to detect and accurately determine their position, resulting in higher uncertainties.

Flux densities and brightness temperatures
Figures 3 and 4 show the flux density (left) and brightness temperature (middle) as a function of time for J0805+6144 and J1658−0739, respectively.
The VLBI component flux densities were determined from the Gaussian model fitting described in Section 3. We present the core flux densities, as well as the sum of the flux densities of the core and the closest jet component (J1), to better describe the innermost radio-emitting region.When the J1 jet component is still close to the core, the flux density of the compact central region may not be well represented by the core component alone.The emergence of a new jet component is often associated with a flux density outburst.For a while, C and the outward-moving J1 are blended together, as the new component cannot be distinguished from the core because of the limited angular resolution of the network.
In the case of J0805+6144 (Fig. 3), the flux density in the central region is dominated by the core, and the J1 flux density stays rather constant.From this, and the generally decreasing trend in the light curve, it is likely that the outburst possibly associated with the ejection of J1 had happened prior to the start of the VLBI observations.The VLBI flux density curve of J1658−0739 (Fig. 4) also indicates an overall decreasing trend.However, during the first half of the monitoring period, the central (C+J1) flux density is dominated by the recently ejected jet component.As J1 advances outward, its flux density decreases, while the core flux density remains stable.The brightening of the core started again around 2012.
The redshift-corrected brightness temperatures were calculated following the equation of Condon et al. (1982): where  is the redshift,   the integrated flux density of the core in Jy,  the observing frequency in GHz,  a and  b the major and minor axes (FWHM) of the fitted elliptical Gaussian in mas, respectively.Kovalev et al. (2005) give a formula for the size of the minimum resolvable source component with the interferometer.To determine   (the half-power beam width measured along an arbitrary position angle ), we followed Appendix B of Heywood et al. (2021).The calculated minimum resolvable size was substituted in Eq. 1 instead of  a and  b as an upper limit if the source was unresolved at a given epoch.This way, the brightness temperatures obtained are lower limits.
When  b exceeds the equipartition value,  b,eq ≈ 5 × 10 10 K (Readhead 1994), the emission is considered relativisically enhanced (Doppler-boosted).The Doppler factors can be calculated by assuming an intrinsic brightness temperature ( b,int ) in the source as where  b is derived from VLBI measurements at a given frequency as described above.We note that the measured brightness temperatures are variable with time (Fig. 3-4), and the intrinsic brightness temperatures may usually be lower than the equipartition value (see e.g.Homan et al. 2006Homan et al. , 2021)).For this reason, we adopted the intrinsic brightness temperature,  b,int = 4.10 × 10 10 K, from Homan et al. (2021), who found  b,int to be close to the equipartition value when taking the median core  b over many observational epochs.Before calculating the Doppler factor, we estimated the median brightness temperatures ( b,median ) based on our measured  b values of the VLBI core at 7.6 − 8.6 GHz (see Table 3), including those that are considered as lower limits.This way, the resulting  b,median values will be reasonable estimates and they are not affected much by some single  b lower limits among the measurements.The radio emission of the core is Doppler boosted in both sources in most of the epochs and the median values exceed  b,int , too.The  b,median for J0805+6144 is 4.59 × 10 11 K, while for J1658−0739,  b,median = 1.11 × 10 11 K.We used these  b,median values in the numerator of Eq. 2 to obtain a single  value for the cores.The Doppler factors determined this way (Table 4) can be considered characteristic to these jets.

Lorentz factors and jet viewing angles
There are two commonly used methods for determining the bulk Lorentz factor in a blazar jet.One is to fit the broad-band SED of a blazar (e.g.Boettcher et al. 1997).The other method, also applied in this paper, is based on VLBI observations of the radio jet, using the Doppler factor and the apparent superluminal speed.The following equations from Ghisellini et al. (1993) can be used to calculate the Lorentz factor and the jet viewing angle with respect to the line of sight: where  is the apparent superluminal speed of the jet component in the units of the speed of light .In the case of J0805+6144, the Lorentz factor is Γ ≈ 7.4 and the inclination angle is  ≈ 4.4 • (Fig. 3).The  b lower limit at the last epoch is below the calculated  b,median and its exact value could still slightly influence the resulted Doppler factor.The Lorentz factor for J1658−0739 is Γ ≈ 36.6 (Fig. 4), a rather high value, but not unprecedented at high redshifts (Zhang et al. 2022).The viewing angle is  ≈ 8.0 • .Both our targets have their jet viewing angle within  ≈ 10 • as expected for typical blazars.

Proper motion-redshift relation for jetted quasars
Distant jetted radio sources above redshift 3 are generally harder to find because of their relative weakness and the lack of bright, prominent mas-scale jet features compared to low-redshift sources.The small number of known sources (e.g.Perger et al. 2017) makes it challenging to reveal an overall picture.The few dedicated highredshift jet proper motion studies carried out with VLBI so far, either for single or multiple objects, are e.g.Frey et al. (2015); Zhang et al. (2017); Perger et al. (2018);An et al. (2020); Zhang et al. (2020Zhang et al. ( , 2022)).The cosmological time dilation requires long VLBI monitoring of the sources in the observer's frame to obtain well-

The Lorentz factors of J0805+6144 and J1658−0739
There are independent Lorentz factor estimates from the literature for both of our targets.For J0805+6144, Paliya et al. ( 2020) derived Γ = 14 from SED fitting.Our Γ ≈ 7.4 is lower by a factor of 2. On the other hand, from a single X-ray observation, Sahakyan et al. (2020) obtained Γ = 1.1 ± 0.1, which is significantly below both values.For J1658−0739, only one other Γ estimate is found in the literature, again based on SED fitting by Paliya et al. (2017).These authors obtained Γ = 10, which is well below our Γ ≈ 36.6.Paliya et al. (2017) also gave an estimate for the jet inclination angle,  = 3 • , in contrast to our  ≈ 8.0 • .In our model with  = 13.8, their jet viewing angle would correspond to  ≃ 20.
The differences between the results from the VLBI analysis and SED fitting might be because different regions of the jet are probed by radio interferometry and X-ray measurements.Also, X-ray emission can be produced at various places, like the corona, the inner jet or large, kpc-scale lobes in AGNs.Regarding VLBI imaging, it is possible that the size of the core is overestimated, leading to the underestimation of the brightness temperatures, thus the Doppler factor, too.Natarajan et al. (2017) showed that it can happen when the VLBI core component is blended with a very nearby jet component and these are unresolved by the interferometer.Moreover, variability can also be a source of uncertainty when estimating physical parameters.

The inner region of J1658-0739
The inner jet region of J1658−0739 showed different faces during the period covered by the VLBI observations Figs. 2,4).Notably, the J2 component has no significant proper motion in the radial direction (Table 4).It could be a standing shock (Courant & Friedrichs 1948;Daly & Marscher 1988) where a stationary knot can appear well-separated from the core.Standing shocks are belived to be formed where the jet recollimates or returns to higher density and pressure.Normally, the external pressure should drop far away from the central engine, causing the jet to become wider and less dense than at the time of the ejection.While this type of morphological feature is well represented at low redshifts based on jet monitoring VLBI programs (e.g.Jorstad et al. 2017;Lister et al. 2021), they are somewhat rare in high-redshift relativistic jets which show superluminal motion.The absence of these features might be caused by the scarce sample of high- jet kinematic studies.More stationary hot spots similar to ours are found in e.g.3C 395 (Waak et al. 1985), 4C 39.25 (Daly & Marscher 1988), also in J0753+4231 (Zhang et al. 2022) at high redshift,  = 3.595.An alternative explanation we should consider for the stationary J2 component is that it could be a projection effect rather than a physical standing shock.The gradual curvature observed in the jet supports the idea that the apparent stationarity of J2 is caused by a jet bending that slows down the flow speed when projected onto the sky plane.
There are two published VLBI images of this source in the literature, taken at different frequencies.Fomalont et al. (2000) presented a 5-GHz VLBA image (angular resolution ∼ 3 mas, dynamic range ∼ 1000 : 1) as part of the VSOP Prelaunch Survey.It shows a bright core component with an extension to the northeast, consistent with the structure seen in our 8.6-GHz map.Another, high-frequency image made at 24 GHz (angular resolution ∼ 1 mas, dynamic range ∼ 150 : 1, Charlot et al. 2010) shows a bright component at the phase centre and a much fainter emission peak towards the southwest at a separation of ∼ 1 mas and a position angle PA = 225 • .Concerning the absolute astrometric positions of J1658−0739, the X-band (8.6 GHz) and K-band (24 GHz) ICRF3 solutions (Charlot et al. 2020) agree with each other within the uncertainties, while the  optical position is separated from them by ∼ 1 mas in the position angle of about 220 • (i.e. to the southeast, Fig. 2).Notably, the Gaia position seems to coincide with the weak secondary component in the 24-GHz image (Charlot et al. 2010), suggesting that this is the actual core location.The apparent discrepancy between the source structure seen in our 8.6-GHz images (a bright compact core and a jet pointing towards the northwest, Fig. 4) and the 24-GHz image of where the brightest feature is in the northwest with a weaker secondary component towards the southeast may be reconciled in the context of the flux density variability and the structural changes observed in the source.Unfortunately, we lack 8.6-GHz observations in the period 2003−2007 when the 24-GHz image was made.However, the decreasing trend in the core+J1 flux density (Fig. 4) and the emergence of a new inner jet component in the 8.6-GHz maps starting from 2008 suggest an outburst in or shortly before 2002 when a new component was born.This could have resulted in the bright feature that appears in the 24-GHz image taken on 2007 Mar 30 (Charlot et al. 2010).Later, as J1 moved further away from the core, the two features became clearly resolved at 8.6 GHz as well.

SUMMARY
In this paper, we presented a kinematic analysis of the radio jets in two bright high-redshift (3 <  < 4) blazars, J0805+6144 and J1658−0739, for the first time.We analysed 7.6-8.6-GHzVLBI data covering nearly two decades of archival observations, supplemented by EVN imaging observations in 2018 and 2019.We used imaging data with mas-scale angular resolution to model the brightness distri-bution of the sources.We identified multiple jet components across the observing epochs in both J0805+6144 and J1658−0739, and estimated their apparent proper motion.The mas-scale jet structure of J1658−0739 contains an apparently stationary component which might be caused by a projection effect or associated with a standing recollimation shock.By interpreting the fitted component flux densities as a function of time, a higher-resolution archival 24-GHz VLBI image (Charlot et al. 2020), and the radio and optical absolute astrometric positions, we propose the explanation that the jet component marked by J1 might have been caused by a prominent outburst in this AGN that happened around 2000.
The measurements allowed us to obtain estimates of the characteristic physical and geometric parameters of these blazar jets, the bulk Lorentz factor and the inclination angle with respect to the line of sight.The results support the blazar nature of these sources.The derived apparent superluminal motions are ranging between 1 ≲  ≲ 14.Our proper motion measurements add to the sparsely sampled high-redshift part of the apparent proper motion-redshift relation.The derived Lorentz factors are consistent with values found in other  > 3 radio-loud AGNs, supporting expectations from the standard ΛCDM cosmological model without requiring unphysically high jet speeds.Regarding the future, larger samples of high- jets are needed to draw statistically robust conclusions.In addition, observations with longer time baselines would help verify stationary features, and multi-frequency VLBI studies could better identify core positions.

Figures 1
Figures1 and 2show examples of 8.6-GHz VLBI images taken on 2018 May 19 for J0805+6144 (experiment ug002h) and on 2017 Aug 1 for J1658−0739 (experiment rv125, Table2).For both sources, the core-jet structures extending up to about 10 mas are typical for blazars, although this is not unique for blazars, and has been observed in various quasars, including those at high redshifts (e.g.Frey et al. 2010;Coppejans et al. 2016;Krezinger et al. 2022).The Gaia DR3 optical coordinates are marked in the radio images with red crosses whose size represents the 1 positional errors.Broad optical emission lines used for the measurement of the quasar redshift are known to be formed in the so-called broad-line region 2 -Observing epoch (year-month-day); Col. 3 -Identifier of the core (C) and jet (J) components, the latter are numbered based on the increasing distance from the core; Col. 4 -Fitted flux density of the component; Col. 5 -Separation of the jet component from core; Col. 6 -FWHM size of the fitted Gaussian model component; the sizes for the core components fitted with an elliptical Gaussian are written as  × , where  and  are the major and minor axes; Col. 7 -Position angle of the component with respect to the core, measured from north to east;Col.8 -Calculated redhsift-corrected brightness temperature; Col. 9 -Position angle of the major axis of the fitted elliptical Gaussian core component, measured from north to east.
columns are as follows: Col. 1 -Source name; Col. 2 -Identifier of the core (C) and jet (J) components, the latter are numbered based on their increasing distance from the core; Col. 3 -Jet radial proper motion; Col. 4 -Rate of change of the jet position angle; Col. 5 -Apparent transverse radial speed in the units of the speed of light; Col. 6 -The median brightness temperature of the core; Col. 7 -Doppler factor calculated from the VLBI core median brightness temperature; Col. 8 -Bulk Lorentz factor; Col. 9 -Jet viewing angle.

Figure 1 .
Figure 1.Proper motion plots and 8.6-GHz VLBI image of J0805+6144.The dashed lines represent the best-fit linear model.The shaded areas represent the 1 uncertainties of each fit.The colouring of the components is the following: J1 -blue, J2 -green, J3 -red.Top left: Radial proper motion of each component.Top right: Jet component position angles as a function of time.Bottom: Locations of the jet components plotted onto the 8.6-GHz VLBI total intensity image made on 2018 May 19 (experiment ug002h).The peak intensity is 402 mJy beam −1 , with the lowest contours drawn at ±1.1 mJy beam −1 .The positive contour levels increase by a factor of 2. The size of the restoring beam is 1.2 mas × 1.0 mas (FWHM) at PA = 13.• 3 (measured from north through east), as indicated in the bottom-right corner.The shading indicates the observing time in MJD (Modified Julian Date) at the given position of each jet component.The red cross marks the Gaia DR3 optical position, its size indicates the formal uncertainties.

Figure 2 .
Figure 2. Proper motion plots and 8.6-GHz VLBI image of J1658−0739.The dashed lines represent the best-fit linear model.The shaded areas represent the 1 uncertainties of each fit.The colouring of the components is the following: J1 -blue, J2 -green, J3 -red, J4 -orange.Top left: Radial proper motion of each component.Top right: Jet component position angles as a function of time.Bottom right: Locations of the jet components plotted onto the VLBI image made on 2017 Aug 1 (experiment rv125).The peak intensity is 226 mJy beam −1 , with the the lowest contours drawn at ±1.7 mJy beam −1 .The positive contour levels increase by a factor of 2. The size of the restoring beam is 1.4 mas × 0.8 mas (FWHM) at PA = −3.• 7, as indicated in the bottom-left corner.The shading indicates the observing time in MJD at the given position of each jet component.The red cross marks the Gaia DR3 optical position, its size indicates the formal uncertainties.

Figure 3 .
Figure 3.The core flux density, redshift-corrected brightness temperature, and jet parameter plots for J0805+6144.Left: Changes in the core flux density during the time covered by the VLBI observations.The fitted flux densities of the core are plotted in black, the sum of the core and the J1 jet component flux densities in red.Middle: The brightness temperature of the core as a function of time.Only the lower limit of the  b is plotted if the source was unresolved with the interferometer.The dotted line shows the median brightness temperature while the dashed line indicates the intrinsic brightness temperature, 4.1 × 10 10 K, adopted from Homan et al. (2021).Right: The Lorentz factor (green) and the jet viewing angle (blue) as a function of the Doppler factor, based on the jet component with the fastest apparent speed.The dashed line corresponds to the estimated Doppler-factor (  ≃ 11.2) .

Figure 4 .
Figure 4.The core flux density, redshift-corrected brightness temperature, and jet parameter plots for J1658−0739.The description is the same as for Fig. 3. Right: The dashed line corresponds to the estimated Doppler-factor (  ≃ 2.7).

Table 1 .
Details of observations used from the EVN project ET036.

Table 2 .
Details of observations of the archival data obtained from the Astrogeo database Source name; Col. 2 -Project code; Col. 3 -Observing frequency; Col. 4 -Observing date; Col. 5 -Observing date in MJD; Col. 6 -Codes of participating radio telescopes.VLBA: Very Long Baseline Array.
Notes: We only present here several rows as examples, the table is available in its entirety in machine-readable form.The columns are as follows: Col. 1 -

Table 3 .
Model-fitting parameters and the calculated physical properties Notes: We only present here several rows as examples, the table is available in its entirety in machine-readable form.The 1 error is given in parentheses for the derived quantities.The columns are as follows: Col. 1 -Source name; Col.

Table 4 .
Kinematic and physical properties of the jets