Astrometric Apparent Motion of High-redshift Radio Sources

Radio-loud quasars at high redshift (z>4) are rare objects in the Universe and rarely observed with Very Long Baseline Interferometry (VLBI). But some of them have flux density sufficiently high for monitoring of their apparent position. The instability of the astrometric positions could be linked to the astrophysical process in the jetted active galactic nuclei in the early Universe. Regular observations of the high-redshift quasars are used for estimating their apparent proper motion over several years. We have undertaken regular VLBI observations of several high-redshift quasars at 2.3 GHz (S band) and 8.4 GHz (X band) with a network of five radio telescopes: 40-m Yebes (Spain), 25-m Sheshan (China), and three 32-m telescopes of the Quasar VLBI Network (Russia) -- Svetloe, Zelenchukskaya, and Badary. Additional facilities joined this network occasionally. The sources have also been observed in three sessions with the European VLBI Network (EVN) in 2018--2019 and one Long Baseline Array (LBA) experiment in 2018. In addition, several experiments conducted with the Very Long Baseline Array (VLBA) in 2017--2018were used to improve the time sampling and the statistics. Based on these 37 astrometric VLBI experiments between 2017 and 2021, we estimated the apparent proper motions of four quasars: 0901+697, 1428+422, 1508+572, and 2101+600.


INTRODUCTION
The third version of the International Celestial Reference Frame (ICRF3) is defined by the positions of 4536 extragalactic radio sources (Charlot et al. 2020) measured by the geodetic Very Long Baseline Interferometry (VLBI) technique in dual-frequency mode (S/X), in the framework of the observing campaigns coordinated by the International VLBI Service for Geodesy and Astrometry (IVS, Nothnagel et al. 2017). These extragalactic radio sources are active galactic nuclei (AGN), presumably with a supermassive black hole in the nucleus, mostly appearing in VLBI images as a set of very compact radio-bright knots (a core and jet components). Due to the physical processes in the nucleus area of the sources, the structure of the AGN on milliarcsec (mas) scales probed by VLBI at GHz frequencies is variable on the timescales down to months, leading to either physical motion of the compact knots, or variations of their brightness, or both. As a result, highly accurate positions of the reference radio sources are unstable on milliarcsecond scales. If a bright jet component is moving persistently along a particular direction for several years, then the effect might be detected as apparent motion of the radio source with a typical magnitude of ∼ 10 − 100 µas yr −1 (e.g. Fey et al. 1997;Feissel-Vernier 2003;Titov 2007;Moór et al. 2011).
High-redshift AGN are tracers of the properties of early cosmological epochs. Thus, the astrometric observing project described in this paper has astrophysical relevance, too. Only a small fraction of high-redshift AGN are strong enough at radio frequencies to be available for VLBI observations using medium-size radio telescopes. Therefore, a multi-epoch astrometry of the high-redshift AGN is a subject of special interest.
Historically, in astrometry the term "proper motion" refers to the physical motion of stars in our Galaxy. This effect is observed as a coordinate shift of the comparatively nearby stars with respect to more distant stars or extragalactic objects. As the quasars are at cosmological distances from the Earth-based observer, the corresponding proper motion induced by the same physical reasons is negligible (Bartel et al. 1986). For instance, even for the closest galaxy with compact radio core (M 81) at a distance of 3 Mpc, the proper motion does not exceed ∼ 10 µas yr −1 (Kardashev 1986) and has not been detected so far. The motion of one or more compact components of the AGN jet, or variations of the brightness that lead to large apparent change of the radio source position are not linked to a bulk physical motion of the whole AGN with respect to the observer. Therefore, this astrometric effect is used to be called "apparent proper motion" to separate it from the stellar proper motion that have a real physical nature.
Apart from the apparent proper motions induced by the relativistic jets that are randomly distributed in direction for a sample of objects, small systematic effects may also appear. In particular, the dipole systematic of several µas yr −1 discloses an impact of the Galactocentric acceleration of the Earth (e.g. Gwinn et al. 1997;Sovers et al. 1998;Titov et al. 2011;MacMillan et al. 2019;Charlot et al. 2020). On the cosmological scale of events, theoretical studies suggest the existence of primordial gravitational waves. The primordial gravitational waves, if they are strong enough to be detected, would cause a quadrupole systematic in the apparent proper motion of the reference radio sources (e.g. Kristian & Sachs 1966;Gwinn et al. 1997;Titov et al. 2011). Moreover, in some cosmological scenarios, the magnitude of the proper motions could grow with redshift. If the cosmological signal is too weak, it might not be detected with the low-redshift radio sources. However, as this signal grows with distance, it might be detected with radio sources at higher redshift. As the number of high-redshift radio sources regularly observed in the geodetic VLBI program is limited, it is not possible to extract a specific systematic signal from a large set of proper motions. The only possible way is to check whether the observed high-redshift sources demonstrate the proper motion that cannot be related to the usual relativistic jet activity.
Therefore, in our observational project, we are focusing on a set of radio sources with redshift z ≥ 4, trying to detect unusually large proper motions on a timescale of several years. Detection of any astrometric anomaly would trigger a follow-up cosmological research.

THE SAMPLE AND OBSERVATIONS
We started composing the sample from a list of previously observed quasars at redshift z ≥ 4 with total flux density at 8.4 GHz S 8.4 ≥ 100 mJy. This simple selection criterion gave us four radio sources: 0901+697, 1428+422, 1508+572, and 2101+600. The sources' celestial positions favoured regular observations of these quasars by a set of radio telescopes in the northern hemisphere. To achieve a good signal-to-noise ratio (SNR) for a fringe detection, we have practiced prolonged tracking of the high-redshift radio sources during the RUA, EVN and LBA experiments, up to 18-20 minutes for selected radio sources. Other 20-25 radio sources used to tie their positions to the ICRF3 were observed in the standard geodetic mode (e.g. Sovers et al. 1998) to achieve SNR ≥ 7. All the sources were selected from the list of 303 ICRF3 "defining" radio sources to provide the best reference frame around the sky with positions known to be better than 30 µas. The Very Long Baseline Array (VLBA) experiments had used short snapshot scans for all radio sources. Table 2 lists the observing sessions, their dates, radio telescopes involved, observed target sources and on-target times. Table 3 explains the codes of radio telescope shown in Table 2. RUA project code is associated with the regular 24 h astrometric sessions by the network of five radio telescopes, including Yebes (Spain), Sheshan (China), Badary, Svetloe, and Zelenchukskaya (Russia). The latter three radio telescopes are the part of the Russian National Quasar VLBI Network (Shuygina et al. 2019). The Hobart (Australia) and Kunming (China) radio telescopes participated in RUA031, the Hartebeesthoek (South Africa) radio telescope participated in RUA032. The experiment V515C was Note-Redshift references: 1 - Romani et al. (2004), 2 -Hook & McMahon (1998), 3 -Hook et al. (1995), 4 -Sowards-Emmerd et al. (2004 conducted using the Australian Long Baseline Array (LBA) -a part of the Australia Telescope National Facility. Three ET036 experiments were carried out by the European VLBI Network (EVN). We also used data from several observing sessions of the Very Long Baseline Array (VLBA) conducted by the United States Naval Observatory (USNO) and the National Aeronautics and Space Administration Goddard Space Flight Center (NASA GSFC) in the projects not related to ours. RUA and EVN sessions listed in Table 2 were scheduled with the sked software (Gipson 2010) in the Institute of Applied Astronomy of the Russian Academy of Sciences (IAA RAS) in Saint-Petersburg. The LBA session (V515C) was scheduled using sked in the Shanghai Astronomical Observatory of the Chinese Academy of Sciences (SHAO CAS). Standard geodetic/astrometric approach of scheduling the numerous scans on the target radio sources mixed with calibrator source scans was used. All RUA and EVN sessions used standard IVS 3.5/13 cm wide-band setup in right circular polarization, except for two experimental RUA sessions that used only 3.5 cm wavelength. The LBA session used modified 3.5/13 cm wide-band setup in right circular polarization. The standard IVS 3.5/13 cm setup includes 8 upper side band and 2 lower side band frequency channels (IFs) at 3.5 cm and 6 upper side band IFs at 13 cm, specifically placed over bandwidths of 720 MHz and 140 MHz, respectively. The 13 cm frequency coverage of LBA session is reduced to 90 MHz due to radio-frequency interference (RFI). Recording setups of the RUA, EVN, LBA, and VLBA sessions are listed in Table 4. Sessions differ by the IF bandwidths, bit sampling, and total sample rate. The RUA042 session used 3.5 cm frequencies of IVS geodetic wide-band setup and included 8 upper side band and 8 lower side band IFs distributed over 720 MHz bandwidth. The RUA043 session used 32-MHz bandwidth IFs forming almost continuous frequency coverage between 8188 and 8972 MHz. The detailed description of the VLBA experiments can be found in Hunt et al. (2021).
Three of the four observed radio sources (the top three in Table 1) have multiple spectroscopic measurements of their redshifts with reasonable confidence. Recently the Sloan Digital Sky Survey (SDSS, York et al. 2000) clearly displayed the dominating Ly-α peak and other characteristic properties of high-redshift objects, e.g. Ly-break at λ = 912Å rest-frame wavelength, Ly-β and CIV emission lines. The only known spectrum of the fourth object, 2101+600, has unreliable identification of the Ly-α emission line (Sowards-Emmerd et al. 2004). However, the Panoramic Survey Telescope and Rapid Response System (Chambers et al. 2016, Pan-STARRS) survey provides photometric magnitudes for 2101+600 as g = 24.9 and r = 22.6. The large colour difference g − r = 2.3 might hint on the drop in continuum between g and r, supporting the high redshift tentatively determined by Sowards-Emmerd et al. (2004).
The redshifts for the three objects (0901+697, 1428+422 and 1508+572) are compromised by damping of the Ly-α line, the blueshift effect of CIV, and low intensity in Ly-β. Therefore, a conservative estimate of the redshift formal errors at z ∼ 3 is about σ z = 0.005 Titov et al. (e.g. 2013Titov et al. (e.g. , 2017.

DATA PROCESSING AND ANALYSIS
The data of RUA sessions were transferred to and correlated at the Correlator Processing Center of the Russian Academy of Sciences at the IAA RAS in Saint-Petersburg. Correlation was performed using ARC (Surkis et al. 2008) and DiFX (Deller et al. 2011) correlators. Session RUA011 was correlated using ARC, DiFX version 2.4.0 was used for sessions RUA012 -RUA025, version 2.5.2 for RUA026 -RUA033, version 2.6.2 for RUA037 -RUA043, respectively. Three EVN sessions of the ET036 project were correlated using the EVN Software Correlator at JIVE (SFXC, Keimpema et al. 2015). The LBA session was correlated at the Pawsey Supercomputing Centre with DiFX.  Table 3. Radio telescopes involved in the study presented in this paper, in the alphabetical order of their codes used in Table 2.  Visibility data of 0.5 s integration time and of 125 kHz spectral resolution correlated using DiFX were post-processed with the pima software (Petrov et al. 2011) following the Data analysis pipeline from the PIMA User Guide 1 . The LBA and EVN output data in geodetic Mark4 format were post-processed with the Haystack Observatory Postprocessing System 2 (HOPS, Whitney et al. 2004).
Standard approach of calibrating the ionospheric delay was used for the RUA observations at 3.5/13 cm wavelength (Sovers & Jacobs 1996) (sessions RUA011 -RUA030). The ionospheric contribution was calculated using daily Global Ionospheric Maps (GIMs) produced by the Center for Orbit Determination in Europe (CODE) Analysis Center (Schaer et al. 1996) for the rest of observations (sessions RUA031 -RUA043). Various research have shown that ionospheric contribution calculated using GIMs and derived from geodetic VLBI are in good agreement. In particular, Hobiger et al. (2006) analysed available VLBI sessions from 1995 to 2006, their estimate of mean difference between the two techniques is less than 3 TECU (∼57 ps). Recently, Gordon (2010) analysed ten 3.5/13 cm IVS RDV (Research and Development) sessions and found a moderate increase of the source position uncertainties at 3.5 cm when using GIMs by ∼0.5 mas. However, his analysis of the thirteen 3.5/13 cm VLBA Calibrator Survey (VCS) sessions did not show a significant difference between the standard dual frequency calibration of the ionosphere and the newly developed GIMs approach. Moreover, several radio sources were detected at 3.5 cm only.
In some RUA experiments, 13 cm band data were seriously corrupted by RFI, making it impossible to get a solution from 3.5/13 cm data. We decided to work only with 3.5 cm data using the algorithm described in Aksim et al. (2019). Taking possible contamination by RFI into account, some of the latest RUA experiments were specially designed as a single 3.5 cm band experiment. That allows us to increase the number of observations of faint target sources due to higher recording data rate and therefore shorter on-source time. We used GIMs to calculate ionospheric contribution for the experiments from 2019.5 to 2021.5. This period corresponds to a low solar activity, in contrary to the above mentioned works where the VLBI data covered almost a full solar cycle of activity. We expect that the ionospheric contribution using GIMs during the period of low solar activity is negligible. To ensure this, we compared source position estimates for a number of RUA experiments in this period and got an agreement within 1-σ formal errors. The total delays obtained from 3.5 cm data along with the ionospheric delay from either dual-band combination or using vertical total electron content (TEC) from GIMs were then astrometrically analysed.
The data of all presented in this paper VLBI experiments were processed with the OCCAM 6.3 software (Titov et al. 2004) in the session-wise mode. The astronomical reduction of the data (precession, nutation, tidal correction, etc.) was realized in accordance with the IERS Conventions 2010 (Petit & Luzum 2010). The astrometric and geodetic parameters were estimated by the least squares collocation method (Titov & Schuh 2000). The list of parameters comprises the Earth pole position components, correction to the Universal Time (UT1-UTC), nutation offsets, and daily positions of the radio telescopes. All the parameters were estimated for each individual experiment along with the position corrections to the high redshift radio sources from Table 1. A priori positions of the radio telescopes were referred to the International Terrestrial Reference Frame (ITRF2014, Altamimi et al. 2016). We estimated the positions of the radio sources from Table 1 only, the positions of all other observed radio sources were fixed to the aus2020a catalogue values 3 . This catalogue comprises new observations during 2018-2020, after the official release of ICRF3. The most essential change with respect to ICRF3 is an improvement of positions of rarely observed radio sources in the southern hemisphere, while positions of the defining radio sources in the northern hemisphere are very close to the original ICRF3 catalogue values. Among ∼ 150 defining radio sources from the northern hemisphere, the positions of more than 100 in the aus2020a catalogue do not differ from the ICRF3 positions by more than 0.025 mas in both coordinates, so one can state that the daily positions of the four high-redshift objects studied in this paper are very close to the ICRF3 system.
A priori dry zenith delays were determined from local atmospheric pressure values (Saastamoinen 1972), which were then mapped to the elevation of the observed sources using the Vienna mapping function (VMF1, Boehm et al. 2006). The elevation cut-off angle was set to 5 • . To calibrate the troposphere instability induced by the non-equilibrium water vapour during a 24-h experiment, we estimated the so-called zenith wet delays with North-South and East-West gradients at each observational epoch (Titov & Schuh 2000). The instability of the hydrogen maser clock at all VLBI stations was also estimated for each observational epoch.
Section 4 provides detailed results for individual sources, along with an interpretation in the astrophysical context. We approximated the time series of daily right ascension and declination components of the four radio sources with linear functions to calculate the apparent proper motion by the conventional weighted least squares method. Variations of 2101 + 600 were also treated in a special way (see below).

0901+697
This is the most distant quasar (z = 5.47, Romani et al. 2004) among those available for geodetic VLBI observations, due to its comparatively high total flux density. We observed 0901+697 in 23 RUA experiments, three EVN experiments, and added two VLBA experiments conducted in independent projects (UF001A and UG002B). The VLBA data are publicly available from the U.S. National Radio Astronomy Observatory archive 4 and the IVS database.
The source has been imaged with VLBI at various frequencies, up to 43 GHz and multiple epochs since 2004 Zhang et al. 2017;Frey et al. 2018;An et al. 2020). Its pc-scale structure is dominated by a prominent compact core. The jet is extended to southwest with a component separated by ∼ 0.7 mas, then apparently bends southwards, ending in a component located within ∼ 1.5 mas from the core. By detecting polarized radio emission on pc scales and estimating jet component proper motions for the first time in this source, An et al. (2020) found that the properties of 0901+697 are the best explained with a nascent jet embedded in and interacting with a dense surrounding matter. Note that we use the term "jet component proper motion" when speaking about an apparent change in the relative separation of core and jet components in VLBI images measured between different observing epochs. This is not to be confused with the apparent proper motion of the quasar as a whole, like the values given in Table 5.
While the position of the inner jet component, possibly marking the location of a standing shock produced by the jet encountering the ambient medium, is nearly stationary, the more distant southern jet component has a small but significant proper motion, 0.019 ± 0.006 mas yr −1 , measured between 2004 and 2018. According to the 2.3-GHz image sensitive to more extended steep-spectrum radio emission (Frey et al. 2018), the jet bending continues towards the southeast, up to ∼ 2.5 mas from the core. According to occasional multi-epoch VLBI observations at various radio frequencies from 8.4 to 43 GHz, the source flux density is changing on time scales of years (see the data in Zhang et al. 2017;An et al. 2020). But the total 15-GHz flux density regularly monitored at the 40-m antenna of the Owens Valley Radio Observatory (OVRO, Richards et al. 2011), which is dominated by the compact flat-spectrum core at this high frequency, only slightly decreased in 2017-2018, the period overlapping with our astrometric observations (see the supplementary figure 1 in An et al. 2020). The comparison of OVRO single-dish and VLBI flux densities obtained with very different angular resolutions is justified by the fact that near-simultaneous 15-GHz measurements of 0901+697 indicate no appreciable difference between them , suggesting negligible contribution of any emission extended to spatial scales beyond those probed by VLBI. This applies for strongly core-dominated variable radio sources in general (e.g. Lister et al. 2011).
No trend was found in any of the coordinates after linear least squares fitting (Table 5). The potential impact of the source structure (the core fading and the outward motion of the southern jet component) needs to be verified with future observations.

1428+422
The source 1428+422 was identified as a high-redshift (z=4.715) quasar by Hook & McMahon (1998). It has been observed with VLBI as a part of the astrometric program occasionally since 1996, but the number of delays in the majority of experiments did not exceed 10, therefore these early data are not included in the current analysis. We collected 7 data points (4 from RUA experiments and 3 from EVN) between 2017 and 2019. No significant proper motion was detected from this set of experiments (Fig. 2).
The most sensitive 2.3-GHz VLBA image of the source shows a core-jet structure extending to ∼ 20 mas in the southwest direction ). The quasar displays remarkable variations in its flux density producing a major outburst in 2005. The 15-GHz VLBA imaging observations performed shortly before and after the outburst did not reveal a birth of a new component in the inner jet (Veres et al. 2010). However, based on 5 epochs of 8.4 GHz VLBI imaging data spanning 22 yr, Zhang et al. (2020) studied the jet kinematics and found two components within ∼ 3 mas separation from the core moving with apparent jet component proper motions 0.017 ± 0.002 mas yr −1 and 0.156 ± 0.015 mas yr −1 , with the quickly fading outer jet component being the fastest. The latter value corresponds to an apparent transverse speed of 19.5 ± 1.9 times the speed of light, marking the jet in 1428+422 one the fastest-moving jets in z > 4 sources among known to date ).

1508+572
The source 1508+572 was identified as a radio-loud quasar at z = 4.30 by Hook et al. (1995). Due to an insufficient angular resolution, its first 5-GHz VLBI EVN image could only provide a hint on the core-jet structre pointing to south (Frey et al. 1997). Since then, the compact mas-scale structure of 1508+572 extending to ∼ 3 mas has been imaged with VLBI at multiple frequencies at several epochs (e.g. Beasley et al. 2002;Petrov et al. 2005;Helmboldt et al. 2007;O'Sullivan et al. 2011;Pushkarev & Kovalev 2012;Gordon et al. 2016). As an example, the 8.6-GHz image displayed in Fig. 3 was made from the RUA032 experiment data obtained on 2019 Jul 13. The Astronomical Image Processing System (AIPS, Greisen 2003) software package was used for calibrating the visibility data in a standard way, and hybrid mapping was performed with the Difmap program (Shepherd 1997). We followed the same data analysis procedure as described in Frey et al. (2018).
The source 1508+572 was monitored frequently with astrometric VLBI in the 1990s but, similarly to the case of 1428+422, the number of experiments and delays dropped dramatically after 2000. However, since 2017, 1508+572 has been observed in a set of astrometric VLBI experiments (Hunt et al. 2021). We have started our own observations of this radio source only from the middle of 2019.
The coordinate time series shown in Fig. 4 displays a steady linear trend with rates −0.121 ± 0.051 mas yr −1 in right ascension and −0.147 ± 0.040 mas yr −1 in declination. The latter estimate is statistically significant (i.e. > 3σ) and it may be related to the jet structure positioned along the north-south direction.
Since jet proper motion studies are not yet available for 1508+572 in the literature, we downloaded calibrated Xband (8.6/8.7 GHz) VLBI visibility data from the Astrogeo archive 5 that overlap in time with our astrometric study (Table 2) (Shepherd 1997), where the brightness distribution of 1508+572 was also modeled with elliptical and circular Gaussian components for the core and the southern jet feature, respectively. Model fits allow us to quantitatively characterise how the source changes with time. For error estimates, we followed Punsly et al. (2021), and took 1/5 Figure 3. A 8.6-GHz VLBI image of 1508+572 obtained with a 5-element global network of radio telescopes (Sv, Zc, Bd, Yb, Sh) on 2019 Jul 13 (project code RUA032). The peak brightness is 341 mJy beam −1 , the lowest contours are at 0.812 mJy beam −1 , the positive contour levels increase by a factor of 2. The elliptical Gaussian restoring beam indicated in the lower left corner has major and minor axes 0.837 mas and 0.707 mas (half-power beam width), respectively, with a major axis position angle −0. • 4, measured from north through east. The mas-scale structure of the source is dominated by the bright core, a jet component is seen towards south.
of the projection of the elliptical Gaussian synthesized beam half-power width along the core-jet direction as the uncertainty of the component separations. The core flux density errors are assumed to be 10%, as they are dominated by the VLBI absolute flux density calibration uncertainty (Punsly et al. 2021). In Fig. 5, the separation of the jet from the core is displayed as a function of time. We fit a linear function to the radial proper motion of the jet component and obtain 0.117 ± 0.078 mas yr −1 . The light curve constructed from the fitted flux densities of the core component is shown in Fig. 6. It seems that the source underwent a period of remarkable brightening, by more than a factor of 2, in 2017-2019, signaling the rising part of an outburst.
Coincidentally, the fitted displacement rate of the southern jet component (0.117 mas yr −1 ) is similar in magnitude to the declination component of the apparent astrometric proper motion (0.147 mas yr −1 southward) within the uncertainties. However, such a weak jet component with fitted flux density of ∼ 10 − 20 mJy, compared to the compact core with ∼ 100 − 300 mJy (Fig. 6), is unlikely to significantly influence the absolute astrometric position of the entire radio source because the jet component accounts for less than ∼ 10% of the total emission. Therefore the displacement rate of the centroid of mas-scale radio emission cannot be larger than about one-tenth of the jet component proper  motion. The small changes in right ascension and declination of 1508+572 we observed in 2017-2020 (Fig. 4) could rather be related to the on-going flux density outburst (Fig. 6) which might be caused by a newly ejected bright inner jet component that is, however, still blended with the core because of the limited angular resolution. Continuing VLBI imaging of 1508+572 could be able to resolve this new component within a couple of years. Here we restricted the study of the mas-scale radio properties of this source to the time interval covered by our astrometric monitoring only. A more detailed, astrophysically motivated analysis extending to longer periods is beyond the scope of this paper.

2101+600
The quasar 2101+600 (z = 4.575) was identified by Sowards-Emmerd et al. (2004), but was not studied in detail at optical wavelengths. In contrast, it was actively observed in radio with VLBI (Beasley et al. 2002;Petrov et al. 2008;Frey et al. 2018;Zhang et al. 2021), and was found to be a gigahertz-peaked spectrum (GPS) source with its radio spectrum peaking near 1 GHz (Coppejans et al. 2017).
Our VLBI coordinate time series in right ascension and declination are shown in Fig. 7. A quasi-periodic variation along the right ascension with period about 3 years was found for 2101+600. Fitting the time series in the right  ascension yields the amplitude 0.28 ± 0.05 mas. Taking this into account, a new estimate of the linear proper motion component along right ascension is −0.55±0.25 mas yr −1 , i.e. about 20% smaller than the estimate without accounting for the quasi-harmonic variation (Table 5). Longer-term observations to cover at least two suspected periods (≥ 6 yr) would be required to verify if this apparent periodicity persists. Although the time series in declination does not show the same systematic variations, we approximated both components with the 3-year harmonic. The corresponding approximation in sky plane is shown in Fig. 8.
VLBI imaging (e.g. Frey et al. 2018;Zhang et al. 2021) shows a structure extended to ∼ 10 mas in the east-west direction, i.e. along the right ascension axis. The overall source emission is dominated by the eastern feature that can, however, be adequately modeled with not a single but 3 circular Gaussian brightness distribution components.
The western feature is comparably weak, and the faintest component in between the eastern and western features could only be recently revealed with sensitive VLBA imaging at 8.4 GHz (Zhang et al. 2021). There is no change (< 0.04 mas yr −1 ) detected in the separation of the eastern and western features over a time range of ∼ 13 yr. In contrast to a core-jet blazar with a relativistically beamed radio jet pointing towards the observer, Zhang et al. (2021) interpret the VLBI imaging data and the peaked overall radio spectrum as characteristic to a young compact  (Fig. 7) were fitted by the linear trend and harmonic signal with the period of 3 years. This approximation is shown by the black curve. The estimates of the amplitude of variations are 0.28 ± 0.05 mas in right ascension and 0.11 ± 0.05 mas in declination. The reference point corresponds to the solution aus2020a.crf. The error bars are from Fig. 7, but referred to the fitting curve here.
symmetric object (CSO) with its bipolar jets misaligned with respect to the line of sight. In this scenario, the weak central component is the core, marking the physical centre of the AGN, similarly to the case of e.g. the bright lowredshift CSO PKS 1117+146 (Bondi et al. 1998). Therefore, the quasar 2101+600 may be morphologically different from other high-redshift sources in our sample. If the CSO interpretation holds, the brightness peak is associated with a complex hotspot in the eastern lobe, rather than a compact synchrotron self-absorbed blazar core. It is the brighter eastern lobe which dominates the astrometric positioning of the source. Figure 8 shows the apparent celestial track of the astrometric solutions obtained in our project. The pattern of this track resembles that of a precessing jet. The phenomenon of mas-scale jet precession in quasars is well known (see e.g. Britzen et al. 2018, 2019, andreferences therein). It is usually attributed to either the presence of a binary supermassive black hole (SMBH) system in the nucleus of the object or the Lense-Thirring effect in a viscous "warped disc" (e.g. Tremaine & Davis 2014).
For a binary nature of the nucleus in the source 2101+600, one might analyse a "toy model" consisting of two SMBHs, with masses m 1 and m 2 . The third Kepler's law for such the binary system can be presented as where a is the semi-major axis of the elliptical orbit of the binary system of the total mass M = m 1 + m 2 , T is its rest-frame period, and G is the gravitational constant. Simple visual analysis of Fig. 8 allows us to estimate the apparent period of precession as T app ≈ 3 yr. Taking into account the cosmological time dilation, it translates into T = T app /(1 + z) ≈ 0.54 yr in the rest-frame of the quasar at z = 4.575. Then the semi-major axis of the binary system can be estimated as a = 0.32 × 10 −2 × M 10 9 × M 1/3 × T 0.54 yr 2/3 pc. (2) For an equal distribution of masses in the binary system, m 1 = m 2 = 0.5 × M = 0.5 × 10 9 M , the semi-major axis would only be a 0.25 × 10 −2 pc, or 50 times the gravitational radius r g = 2Gm 1 c −2 ≈ 0.48 × 10 −4 pc, where c is the speed of light. We intentionally normalize the estimate of the semi-major axis a for the high value of the nucleus mass to underline that such the binary system would be a source of appreciable gravitational waves (e.g. De Rosa et al. 2019). In reality, the mass is likely to be lower. But in any case, the source 2101+600 is an interesting target for studies of the growth of SMBH in early cosmological epochs. Further astrophysical analysis of this source is outside of the scope of this paper. However, given the apparent variations of the astrometric position of 2101+600, illustrated by Fig. 7, this object should not be used as a reference source involved in the formation of reference frames.

More Detailed Analysis of the Positional Time Series
The resultant position of a radio source is related to the position of the phase center as detected on the interferometer baseline. While for a point-like radio source the interferometer response is the same for all baseline lengths and orientations, the response is more complicated for an extended radio source (Charlot 1990). Since all four radio sources presented in this study display an extended structure, the estimates of each radio source position depend on the set of baselines actually tracking the source. For the custom-scheduled RUA experiments, we attempted to observe the three radio sources (0901+697, 1508+572, and 2101+600) by all available stations in the same mode whenever it was possible. However, the estimates presented in Table 5 may be affected by the network geometry. The longer baselines provide better angular resolution of the extended source structure, and, consequently, are more sensitive to the small-scale structural variability. On the other hand, the longer baselines have limited observing time compared to shorter baselines due to sky visibility limitations for non-circumpolar sources. Therefore, quantification of the impact of the long baselines on the final position estimates is a challenging task. Table 6 shows the estimates of the apparent proper motion for the three sources observed in 25 RUA experiments in two versions: for all baselines and for baselines shorter than 6000 km only. (The fourth source, 1428+422 was rarely observed in RUA experiments, therefore it is not included in this analysis.) The 6000-km threshold cuts off the longest baselines (Ys-Bd, Ys-Sh, Zc-Sh, Sv-Sh). This effectively reduces the maximum resolution capability of the total RUA network, and suppresses the effect of the structure as close components remain unresolved at the shorter baselines. This makes the daily estimates of the radio sources positions internally more consistent.
The proper motions for 0901+697 and 2101+600 in the first column of Table 6 are very close to the proper motions from Table 5 because almost all experiments were done with the RUA network. The difference in the estimates is within the reported 3-σ uncertainty. Exclusion of the long baselines from the subset of data leads to a marginal change in proper motions of 0901+697. The difference between 2101+600 proper motions in Columns 3 and 5 (Table 6) is more notable (0.114 mas) but it is still within the combined uncertainty (0.132 mas). One could conclude that baselines of different lengths are likely to produce different positions of the phase center for 2101+600, though the total effect on the proper motion based on the 4-yr period of time is negligible. We believe that the apparently precessing motion of the jet (Figure 8) plays a more important role for the analysis of the linear trend in the position of the radio source 2101+600.  Kristian & Sachs (1966) showed that as the apparent proper motion induced by anisotropic expansion of the Universe may be proportional to the redshift or may be not, the primordial gravitational waves, if they exist in the early Universe, would cause a redshift-dependent systematic effect in the apparent proper motion. The 5/6 of the systematic is to appear in second-order transverse vector spherical harmonics (quadrupole systematic pattern around the sky) equally distributed between the "electric" and "magnetic" components, and the minor fraction of the mean squared proper motion resides in the higher harmonics (Pyne et al. 1996;Gwinn et al. 1997). The second-order spherical harmonics were not found during the search among all radio sources in the past (e.g. Gwinn et al. 1997;Titov et al. 2011;Mignard & Klioner 2012), therefore, in this paper we focus on the most distant targets to detect potential anomalies. Any unusually large values of apparent proper motion of one of the high-redshift radio sources may hint that the primordial gravitational waves are strong enough to be detected.
The apparent proper motion in the positions of the reference radio sources due to the variations of the intrinsic brightness distribution is commonly within 0.1 mas yr −1 (Feissel-Vernier 2003) and rarely exceeds 0.5 mas yr −1 . In an extreme scenario, the positional evolution with amplitude up to 1 mas yr −1 was detected on a very short time scale, less than 1 yr (Titov 2007) due to, presumably, a brightening of a moving relativistic jet. However, this large value of motion in one direction is followed by a similar proper motion in opposite direction after fading out of the moving component, making the accumulated positional displacement is almost zero. The typical pattern of the positional variations over long period of time ( 20 years) is a piece-wise function made of interval length of 1-10 years (Titov 2007). As a result, the total apparent proper motion induced by the change in the intrinsic structure is suppressed over the whole period of observations. It also should be noted that the expansion of the Universe imposes a dilation factor (1 + z) −1 , which results in the apparent "slowing down" of the processes in distant sources. Therefore, one could generally consider 0.1 mas yr −1 as an upper limit on the linear drift for a typical radio source at z ≈ 1. Then, assuming that the systematic signal due to the primordial gravitational waves grows is proportional to the redshift, an apparent proper motion of 0.4 mas yr −1 , being persistent over several years, may be considered as anomalous at z ≈ 4, if no substantial structure of the given radio source is detected. The goal of this study was to monitor the astrometric positions of high-redshift radio sources and to detect any anomalous linear drift in the apparent positions. Assuming the persistent linear drift has an amplitude of 0.5 mas yr −1 , the total displacement of a radio source during a 4-yr interval would reach 2 mas.
The median redshift of sources observed with VLBI and involved in astrometric and geodetic studies is about z = 1. Only a few of the objects in this sample (less than 1%) have redshift z > 4. (Titov & Malkin 2009;Malkin 2018), We monitored the four radio sources presented here (0901+697, 1428+422, 1508+572, and 2101+600) for more than four years using large VLBI facilities around the Earth, and added VLBA observations of the source 1508+572 during 2017-2018.
We did not find any sign of unusually large apparent proper motion among the radio sources under study. However, their high-resolution VLBI images known from the literature reveal an extended structure for all the objects. The estimates of apparent proper motion (Table 5) do not exceed their statistical 3σ threshold of significance, except for the declination component of the radio source 1508+572. These minor, mostly insignificant positional variations tend to align with the intrinsic structure revealed by VLBI images.
The four frequently observed radio sources display extended, changing radio structure. The amplitudes of the estimated apparent proper motions in this study vary typically in a range between 2σ − 3σ ( Table 5). The linear trends in right ascension or declination may exceed the statistical 3σ threshold of significance, were the radio sources observed for an additional 2 − 3 yr. An apparent quasi-periodic variation was detected in the right ascension of 2101+600. Its existence could also be confirmed with longer-term astrometric observations. Finally, the apparent proper motion as well as individual variations of the intrinsic structure do not reveal extraordinary features that may be attributed to a cosmological origin.
We did not find the hypothetical anomalies of the cosmological nature in the positional change of the four radio sources studied here at 4.3 ≤ z ≤ 5.5. This is consistent with the conclusion by Makarov & Secrest (2022) based on the optical study of astrometric proper motions of 60410 quasars at redshifts 0.5 ≤ z ≤ 7.03. Our sample presented here is too small to allow us to make a definitive conclusion. However, we note that dozens of other known AGN targets at z ≥ 4.0 have been imaged more recently with VLBI (see, e.g., Krezinger et al. 2022). Those sources as well as yet to be identified even more distant AGN at redshifts beyond the current de facto frontier of the "known" Universe at around z = 7 might become legitimate targets for future astrometric studies of anomalous proper motions.

ACKNOWLEDGMENTS
The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.
The recent VLBA experiments were run by the geodetic group of the US Naval Observatory to monitor the radio reference frame sources with ten VLBA antennas. The authors acknowledge use of the VLBA under the US Naval Observatory's time allocation. This work supports USNO's ongoing research into the celestial reference frame and geodesy.
The RUA experiments were organised by the Institute of Applied Astronomy of the Russian Academy of Sciences and made use of an ad-hoc VLBI network which consists of three 32-m radio telescopes Badary, Svetloe, and Zelenchukskaya together with the 25-m radio telescope Sheshan of the Shanghai Astronomical Observatory of the Chinese Academy of Sciences, and the 40-m radio telescope of Yebes Observatory (Instituto Geográfico National, Spain). Badary, Svetloe, and Zelenchukskaya radio telescopes are operated by the Scientific Equipment Sharing Center of the Quasar VLBI Network.
The European VLBI Network is a joint facility of independent European, African, Asian and North American radio astronomy institutes. The EVN observations presented in this paper have been conducted under the project code ET036.
The Long Baseline Array is part of the Australia Telescope National Facility (grid.421683.a) which is funded by the Australian Government for operation as a National Facility managed by CSIRO. This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia.