# Astrometric Apparent Motion of High-redshift Radio Sources

Oleg Titov Geoscience Australia, PO Box 378, Canberra 2601, Australia  Sándor
Frey Konkoly Observatory, ELKH Research Centre for Astronomy and Earth
Sciences, Konkoly Thege Miklós út 15-17, 1121 Budapest, Hungary  CSFK, MTA
Centre of Excellence, Konkoly Thege Miklós út 15-17, 1121 Budapest, Hungary
Institute of Physics, ELTE Eötvös Loránd University, Pázmány Péter sétány 1/A,
H-1117 Budapest, Hungary  Alexey Melnikov Institute of Applied Astronomy,
Russian Academy of Sciences, Kutuzova Embankment 10, St. Petersburg, 191187,
Russia  Fengchun Shu Shanghai Astronomical Observatory of the Chinese Academy
of Science, 80, Nandan Road, Shanghai, China  Bo Xia Shanghai Astronomical
Observatory of the Chinese Academy of Science, 80, Nandan Road, Shanghai,
China  Javier González  Observatorio de Yebes (IGN), Apartado 148, 19180
Yebes, Spain  Belén Tercero  Observatorio de Yebes (IGN), Apartado 148, 19180
Yebes, Spain  Leonid I. Gurvits Joint Institute of VLBI ERIC, Oude
Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands  Faculty of Aerospace
Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The
Netherlands  Aletha de Witt South African Radio Astronomy Observatory, 2 Fir
Street, Black River Park, Cape Town, 7925, South Africa  Jamie McCallum
University of Tasmania, Hobart, Australia  Mikhail Kharinov Institute of
Applied Astronomy, Russian Academy of Sciences, Kutuzova Embankment 10, St.
Petersburg, 191187, Russia  Vladimir Zimovsky Institute of Applied Astronomy,
Russian Academy of Sciences, Kutuzova Embankment 10, St. Petersburg, 191187,
Russia  Máté Krezinger Department of Astronomy, Institute of Geography and
Earth Sciences, ELTE Eötvös Loránd University, Pázmány Péter sétány 1/A,
H-1117 Budapest, Hungary  Konkoly Observatory, ELKH Research Centre for
Astronomy and Earth Sciences, Konkoly Thege Miklós út 15-17, 1121 Budapest,
Hungary

###### Abstract

Radio-loud quasars at high redshift ($z\geq 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-2018$
were 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.

galaxies: active — galaxies: high-redshift — radio continuum: galaxies —
astrometry

††software: sked (Gipson, 2010), SFXC (Keimpema et al., 2015), PIMA (Petrov et
al., 2011), OCCAM (v6.3, Titov et al., 2004), AIPS (Greisen, 2003), Difmap
(Shepherd, 1997)

##  1 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 $\sim 10-100$ $\mu$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 $\sim 10$ $\mu$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 $\mu$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\geq 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.

##  2 The sample and observations

We started composing the sample from a list of previously observed quasars at
redshift $z\geq 4$ with total flux density at 8.4 GHz $S_{8.4}\geq 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 $\geq 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 $\mu$as. The Very Long
Baseline Array (VLBA) experiments had used short snapshot scans for all radio
sources.

Table 1: List of the high-redshift radio sources discussed in this paper. Positions are from the aus2021b.crf solution (see footnote [3]) produced by the Geoscience Australia IVS Analysis Centre.  Source | $z$ | Ref. | Right ascension (2000.0) | Declination (2000.0) | No. of experiments
---|---|---|---|---|---
0901$+$697 | 5.47 | 1 | 09 06 30.74876 |  $+$69 30 30.8287 | 30
1428$+$422 | 4.715 | 2 | 14 30 23.74163 |  $+$42 04 36.4911 | 25
1508$+$572 | 4.301 | 3 | 15 10 02.92237 |  $+$57 02 43.3759 | 104
2101$+$600 | 4.575 | 4 | 21 02 40.21905 |  $+$60 15 09.8366 | 21
Table 2:  The list of 37 experiments in which the high-redshift radio sources from Table LABEL:table_sources were observed. Telescope codes listed in the column “Array” are explained in Table LABEL:table_stations. The maximum baseline length is shown in the last column.  Session | Date | Array | High-$z$ sources (on-source time, h) | Maximum baseline (km)
---|---|---|---|---
UF001A | 2017 Jan 16 | VLBA | 0901$+$697 (0.08) | 8611
RUA011 | 2017 Jan 28 | Sv-Zc-Bd | 0901$+$697 (3.3) | 4405
RUA012 | 2017 Feb 04 | Sv-Zc-Bd-Ys | 0901$+$697 (2.48), 2101$+$600 (2.87) | 7080
RUA015 | 2017 Apr 08 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (3.89), 2101$+$600 (5.06) | 9170
RUA016 | 2017 Apr 22 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (3.28), 2101$+$600 (4.08) | 9170
UF001G | 2017 Apr 28 | VLBA | 1508$+$572 (0.02) | 8611
UF001H | 2017 May 01 | VLBA | 1508$+$572 (0.02) | 8611
UF001O | 2017 Aug 05 | VLBA | 1508$+$572 (0.02) | 8611
RUA017 | 2017 Aug 19 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (3.1), 1428$+$422 (1.12), 2101$+$600 (4.08) | 9170
RUA018 | 2017 Aug 26 | Sv-Bd-Ys-Sh | 0901$+$697 (4.84), 1428$+$422 (1.57), 2101$+$600 (4.67) | 9170
RUA019 | 2017 Nov 18 | Sv-Bd-Ys-Sh | 0901$+$697 (5.23), 1428$+$422 (0.07), 2101$+$600 (4.86) | 9170
RUA020 | 2017 Nov 25 | Sv-Zc-Bd-Ys | 0901$+$697 (5.04), 1428$+$422 (0.07), 2101$+$600 (4.67) | 7080
UG002B | 2018 Jan 25 | VLBA | 0901$+$697 (0.07), 1508$+$572 (0.03) | 8611
RUA022 | 2018 Feb 03 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (3.88), 2101$+$600 (5.05) | 9170
RUA023 | 2018 Mar 31 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (3.68), 2101$+$600 (5.59) | 9170
UG002H | 2018 May 19 | VLBA | 1508$+$572 (0.02) | 8611
UF002I | 2018 Jun 03 | VLBA | 1508$+$572 (0.02) | 8611
ET036A | 2018 Jun 04 | Wn-Mc-Nt-On-Ef- | 0901$+$697 (3.75), 1428$+$422 (3.48), |
|  | Sv-Zc-Bd-Ys-Hh | 2101$+$600 (3.87) | 7080
RUA024 | 2018 Jun 23 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (4.26), 2101$+$600 (7.18) | 9170
V515C | 2018 Jul 19 | Sv-Zc-Bd-Sh-Km- | 0901$+$697 (1.25), |
|  | Ho-Hh | 2101$+$600 (2.0) | 6711
RUA025 | 2018 Jul 21 | Sv-Zc-Bd-Ys | 0901$+$697 (3.68), 2101$+$600 (5.38) | 7080
RUA026 | 2018 Sep 08 | Sv-Zc-Bd-Ys | 0901$+$697 (3.88), 2101$+$600 (5.39) | 7080
ET036B | 2018 Oct 27 | Wz-Mc-On-Ef-Sv- | 0901$+$697 (3.73), 1428$+$422 (1.45), |
|  | Zc-Bd-Ys-Hh-T6 | 2101$+$600 (4.28) | 9170
UG002T | 2018 Nov 18 | VLBA | 1508$+$572 (0.02) | 8611
RUA028 | 2019 Jan 05 | Sv-Zc-Bd-Ys | 0901$+$697 (4.07), 2101$+$600 (4.49) | 7080
ET036C | 2019 Feb 25 | Wn-Mc-On-Ef-Sv- | 0901$+$697 (3.38), 1428$+$422 (1.65), |
|  | Zc-Bd-Ys-Hh-T6-Ur | 2101$+$600 (3.62) | 9170
RUA030 | 2019 Jun 29 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (3.29), 1508$+$572 (2.21), 2101$+$600 (5.69) | 9170
RUA031 | 2019 Jul 06 | Sv-Zc-Bd-Km-Ho | 0901$+$697 (2.52), 1508$+$572 (2.08), |
|  |  | 2101$+$600 (3.89) | 6761
RUA032 | 2019 Jul 13 | Sv-Zc-Bd-Ys-Sh-Hh | 0901$+$697 (1.94), |
|  |  | 1508$+$572 (2.34) | 9170
RUA033 | 2019 Sep 28 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (2.91), 1508$+$572 (2.62), 2101$+$600 (4.79) | 9170
RUA037 | 2020 Jan 25 | Sv-Zc-Bd-Ys | 0901$+$697 (3.26), 1508$+$572 (2.21), 2101$+$600 (5.03) | 7080
RUA038 | 2020 Jun 27 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (4.56), 1508$+$572 (1.88), 2101$+$600 (5.04) | 9170
RUA039 | 2020 Nov 28 | Sv-Zc-Bd-Ys | 0901$+$697 (2.13), 1508$+$572 (2.21), 2101$+$600 (3.29) | 7080
RUA040 | 2020 Dec 26 | Sv-Zc-Bd-Ys-Sh | 0901$+$697 (2.71), 1508$+$572 (1.77), 2101$+$600 (5.39) | 9170
RUA041 | 2021 Mar 20 | Sv-Zc-Bd-Sh | 0901$+$697 (2.99), 1508$+$572 (0.69), 2101$+$600 (3.97) | 6761
RUA042 | 2021 Mar 27 | Sv-Zc-Bd-Sh | 0901$+$697 (1.34), 1508$+$572 (0.24), 2101$+$600 (2.59) | 6761
RUA043 | 2021 Apr 10 | Sv-Zc-Bd-Sh | 0901$+$697 (0.74), 2101$+$600 (1.37) | 6761
Table 3: Radio telescopes involved in the study presented in this paper, in the alphabetical order of their codes used in Table LABEL:table_sessions.  Telescope | Telescope | Telescope | Institute or observatory | Country
---|---|---|---|---
code | name | diameter [m] |  |
Bd | Badary | 32 | Institute of Applied Astronomy | Russia
Ef | Effelsberg | 100 | Max Planck Institute for Radio Astronomy | Germany
Hh | Hartebeesthoek | 26 | South African Radio Astronomy Observatory | South Africa
Ho | Hobart | 26 | Mt. Pleasant Radio Observatory, Tasmania | Australia
Km | Kunming | 40 | Yunnan Astronomical Observatory | China
Mc | Medicina | 32 | INAF Institute of Radio Astronomy | Italy
Nt | Noto | 32 | INAF Institute of Radio Astronomy | Italy
On | Onsala | 20 | Onsala Space Observatory | Sweden
Sh | Sheshan | 25 | Shanghai Astronomical Observatory | China
Sv | Svetloe | 32 | Institute of Applied Astronomy | Russia
T6 | Tianma | 65 | Shanghai Astronomical Observatory | China
Ur | Nanshan (Urumqi) | 25 | Xinjiang Astronomical Observatory | China
VLBA | VLBA | $10\times 25$ | National Radio Astronomy Observatory | USA
Wn | Wettzel-North | 13 | Geodetic Observatory Wettzell | Germany
Wz | Wettzel | 20 | Geodetic Observatory Wettzell | Germany
Ys | Yebes | 40 | Yebes Observatory | Spain
Zc | Zelenchukskaya | 32 | Institute of Applied Astronomy | Russia

Table LABEL:table_sessions lists the observing sessions, their dates, radio
telescopes involved, observed target sources and on-target times. Table
LABEL:table_stations explains the codes of radio telescope shown in Table
LABEL:table_sessions. 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 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.

Table 4: Recording setups of the sessions Session | Sample rate | Sampling | 3.5 cm | 13 cm
---|---|---|---|---
names | Mbit s-1 | bits | setup | setup
RUA011 | 256 | 1 |  $10\times 8$ MHz |  $6\times 8$ MHz
ET036A–ET036C | 512 | 2 |  $10\times 8$ MHz |  $6\times 8$ MHz
RUA012–RUA041 | 1024 | 2 |  $10\times 16$ MHz |  $6\times 16$ MHz
V515C | 1024 | 2 |  $10\times 16$ MHz |  $6\times 16$ MHz
RUA042 | 1024 | 2 |  $16\times 16$ MHz | —
RUA043 | 2048 | 2 |  $16\times 32$ MHz | —
UF001A–UG002T | 2048 | 2 |  $12\times 32$ MHz |  $4\times 32$ MHz

RUA and EVN sessions listed in Table LABEL:table_sessions 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 LABEL:table_setups. 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
LABEL:table_sources) 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-$\alpha$ peak
and other characteristic properties of high-redshift objects, e.g. Ly-break at
$\lambda=912$ Å rest-frame wavelength, Ly-$\beta$ and CIV emission lines. The
only known spectrum of the fourth object, 2101$+$600, has unreliable
identification of the Ly-$\alpha$ 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-$\alpha$ line, the blueshift effect of
CIV, and low intensity in Ly-$\beta$. Therefore, a conservative estimate of
the redshift formal errors at $z\sim~{}3$ is about $\sigma_{z}=0.005$ Titov et
al. (e.g. 2013, 2017).

##  3 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.

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
Guide111http://astrogeo.org/pima/pima_user_guide.html, accessed on
2022.10.09.. The LBA and EVN output data in geodetic Mark4 format were post-
processed with the Haystack Observatory Postprocessing
System222https://www.haystack.mit.edu/haystack-observatory-postprocessing-
system-hops/, accessed on 2022.10.09. (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 ($\sim$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 $\sim$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-$\sigma$ 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
LABEL:table_sources. 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
LABEL:table_sources only, the positions of all other observed radio sources
were fixed to the aus2020a catalogue
values333https://cddis.nasa.gov/archive/vlbi/ivsproducts/crf/. 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 $\sim 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\arcdeg$. 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).

##  4 Observational Results

###  4.1 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 archive444https://data.nrao.edu/ and the
IVS database.

The source has been imaged with VLBI at various frequencies, up to 43 GHz and
multiple epochs since 2004 (Romani et al., 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 $\sim 0.7$ mas, then apparently bends southwards, ending in a
component located within $\sim 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
LABEL:table_pm.

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\pm 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 $\sim 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 (An et al., 2020), 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 LABEL:table_pm). 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.

Figure 1: Relative astrometric coordinates (right ascension and declination)
of 0901$+$697 as a function of time, between 2017 and 2020. The reference
point corresponds to the solution aus2020a.crf. The error bars show $\pm
1\sigma$ formal uncertainty.

###  4.2 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 $\sim 20$ mas in the southwest direction (Zhang et al., 2020).
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 $\sim 3$ mas separation from the core moving
with apparent jet component proper motions $0.017\pm 0.002$ mas yr-1 and
$0.156\pm 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\pm 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 (Zhang et al., 2020).

Figure 2: Time series of the relative right ascension (left) and declination
(right) of 1428$+$422\. The reference point corresponds to the solution
aus2020a.crf. The error bars show $\pm 1\sigma$ formal uncertainty.

###  4.3 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 $\sim 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\pm 0.051$ mas yr-1 in right ascension and $-0.147\pm 0.040$ mas
yr-1 in declination. The latter estimate is statistically significant (i.e.
$>3\sigma$) and it may be related to the jet structure positioned along the
north–south direction.

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\fdg 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.

Figure 4: Relative astrometric coordinates of 1508$+$572 as a function of
time, between 2017 and 2020. The reference point corresponds to the solution
aus2020a.crf. The error bars show $\pm 1\sigma$ formal uncertainty.

Since jet proper motion studies are not yet available for 1508$+$572 in the
literature, we downloaded calibrated $X$-band (8.6/8.7 GHz) VLBI visibility
data from the Astrogeo archive555http://astrogeo.org/cgi-
bin/imdb_get_source.csh?source_name=1508%2B572, maintained by Leonid Petrov
that overlap in time with our astrometric study (Table LABEL:table_sessions).
These data from VLBA experiments UF001G (2017 Apr 29), UF001H (2017 May 2),
UF001N (2017 Jul 17), UF001O (2017 Aug 5) (Hunt et al., 2021), UG002B (2018
Jan 25), UG002H (2018 May 20), UG002I (2018 Jun 3), UG002R (2018 Sep 25),
UG002S (2018 Nov 3), UG002T (2018 Nov 18), UG002W (2018 Dec 29), and UG002X
(2019 Jan 21) were supplemented with data from RUA experiments where we
obtained sufficient interferometric visibility points for good-quality
imaging, RUA032 (2019 Jul 13, Fig. 3) and RUA033 (2019 Sep 28). Imaging was
done in Difmap (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$ 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\pm 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.

Figure 5: Angular separation of the fitted core and jet components in
1508$+$572 at 8.6/8.7 GHz as a function of time. The line indicates the linear
least squares fit with the slope $0.117\pm 0.078$ mas yr-1. Figure 6: Fitted
core component flux density in 1508$+$572 at 8.6/8.7 GHz as a function of
time.

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
$\sim 10-20$ mJy, compared to the compact core with $\sim 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
$\sim 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.

###  4.4 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).

Figure 7: Relative astrometric coordinates (right ascension and declination)
of 2101$+$600 as a function of time, between 2017 and 2020. The reference
point corresponds to the solution aus2020a.crf. The error bars show $\pm
1\sigma$ formal uncertainty. Figure 8: Evolution of the position of the source
2101$+$600 in the sky plane over the period from 2017.0 to 2021.2. The right
ascension and declination time series (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\pm
0.05$ mas in right ascension and $0.11\pm 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.

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\pm 0.05$ mas. Taking this into account, a
new estimate of the linear proper motion component along right ascension is
$-0.55\pm 0.25$ mas yr-1, i.e. about $20\%$ smaller than the estimate without
accounting for the quasi-harmonic variation (Table LABEL:table_pm). Longer-
term observations to cover at least two suspected periods ($\geq 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 $\sim 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 $\sim 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 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 low-redshift 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, and references 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

$a=\left(\frac{GMT^{2}}{4\pi^{2}}\right)^{1/3},$ (1)

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_{\rm app}\approx 3$ yr.
Taking into account the cosmological time dilation, it translates into
$T=T_{\rm app}/(1+z)\approx 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\times 10^{-2}\times\left(\frac{M}{10^{9}\times
M_{\odot}}\right)^{1/3}\times\left(\frac{T}{0.54\,\mathrm{yr}}\right)^{2/3}\,\,\,\mathrm{pc}.$
(2)

For an equal distribution of masses in the binary system,
$m_{1}=m_{2}=0.5\times M=0.5\times 10^{9}M_{\odot}$, the semi-major axis would
only be $a\simeq 0.25\times 10^{-2}$ pc, or $\simeq 50$ times the
gravitational radius $r_{\mathrm{g}}=2Gm_{1}c^{-2}\approx 0.48\times 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.

Table 5: Proper motion estimates of the four quasars with 1-$\sigma$ formal errors for all experiments. Quasar | PM in RA | PM in Dec | Number of experiments
---|---|---|---
| mas yr-1 | mas yr-1 |
0901$+$697 | $-0.023\pm 0.023$ | $-0.066\pm 0.023$ | 28
1428$+$422 | $-0.098\pm 0.120$ | $-0.334\pm 0.173$ | 7
1508$+$572 | $-0.121\pm 0.051$ | $-0.147\pm 0.040$ | 18
2101$+$600 | $-0.066\pm 0.036$ | $+0.079\pm 0.035$ | 26

###  4.5 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 LABEL:table_pm 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 LABEL:table_pm_3 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
LABEL:table_pm_3 are very close to the proper motions from Table
LABEL:table_pm because almost all experiments were done with the RUA network.
The difference in the estimates is within the reported 3-$\sigma$ 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 LABEL:table_pm_3) 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.

Table 6: Proper motion estimates of three quasars with 1-$\sigma$ formal errors. Columns 2 and 3 – from RUA experiments only, including all baselines. Columns 4 and 5 – from RUA experiments only, for baselines shorter than $6000$ km Quasar | PM in RA | PM in Dec | PM in RA | PM in Dec | Number of
---|---|---|---|---|---
| mas yr-1 | mas yr-1 | mas yr-1 | mas yr-1 | experiments
0901$+$697 | $-0.003\pm 0.029$ | $-0.069\pm 0.027$ | $-0.028\pm 0.023$ | $-0.063\pm 0.032$ | 25
1508$+$572 | $-0.243\pm 0.086$ | $-0.341\pm 0.131$ | $-0.039\pm 0.039$ | $-0.132\pm 0.093$ | 10
2101$+$600 | $-0.080\pm 0.041$ | $+0.116\pm 0.037$ | $-0.027\pm 0.031$ | $+0.002\pm 0.024$ | 24

##  5 Summary and Conclusion

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
(Titov et al., 2011). 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 ($\gtrsim 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\approx 1$. Then, assuming that the systematic
signal due to the primordial gravitational waves grows is proportional to the
redshift, an apparent proper motion of $\gtrsim 0.4$ mas yr-1, being
persistent over several years, may be considered as anomalous at $z\approx 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 LABEL:table_pm) do not exceed their
statistical $3\sigma$ 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\sigma-3\sigma$ (Table
LABEL:table_pm). The linear trends in right ascension or declination may
exceed the statistical $3\sigma$ 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\leq z\leq
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\leq z\leq 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\geq 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.

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. We acknowledge the use of archival calibrated VLBI data from the
Astrogeo Center database maintained by Leonid Petrov. We are grateful to
Sergei Kurdubov from the Institute of Applied Astronomy RAS for a fruitful
discussion. SF was supported by the Hungarian National Research, Development
and Innovation Office (OTKA K134213). We acknowledge with gratitude very
useful comments provided by the anonymous reviewer of the manuscript of this
paper. This paper is published with the permission of the CEO of the
Geoscience Australia.

## References

  * Aksim et al. (2019) Aksim, D., Melnikov, A., Pavlov, D., & Kurdubov, S. 2019, The Astrophysical Journal, 885, 159, doi: 10.3847/1538-4357/ab499a
  * Altamimi et al. (2016) Altamimi, Z., Rebischung, P., Métivier, L., & Collilieux, X. 2016, Journal of Geophysical Research (Solid Earth), 121, 6109, doi: 10.1002/2016JB013098
  * An et al. (2020) An, T., Mohan, P., Zhang, Y., et al. 2020, Nature Communications, 11, 143, doi: 10.1038/s41467-019-14093-2
  * Bartel et al. (1986) Bartel, N., Herring, T. A., Ratner, M. I., Shapiro, I. I., & Corey, B. E. 1986, Nature, 319, 733, doi: 10.1038/319733a0
  * Beasley et al. (2002) Beasley, A. J., Gordon, D., Peck, A. B., et al. 2002, ApJS, 141, 13, doi: 10.1086/339806
  * Boehm et al. (2006) Boehm, J., Werl, B., & Schuh, H. 2006, Journal of Geophysical Research (Solid Earth), 111, B02406, doi: 10.1029/2005JB003629
  * Bondi et al. (1998) Bondi, M., Garrett, M. A., & Gurvits, L. I. 1998, MNRAS, 297, 559, doi: 10.1046/j.1365-8711.1998.01518.x
  * Britzen et al. (2018) Britzen, S., Fendt, C., Witzel, G., et al. 2018, MNRAS, 478, 3199, doi: 10.1093/mnras/sty1026
  * Britzen et al. (2019) Britzen, S., Fendt, C., Böttcher, M., et al. 2019, A&A, 630, A103, doi: 10.1051/0004-6361/201935422
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560.  https://arxiv.org/abs/1612.05560
  * Charlot (1990) Charlot, P. 1990, AJ, 99, 1309, doi: 10.1086/115419
  * Charlot et al. (2020) Charlot, P., Jacobs, C. S., Gordon, D., et al. 2020, A&A, 644, A159, doi: 10.1051/0004-6361/202038368
  * Coppejans et al. (2017) Coppejans, R., van Velzen, S., Intema, H. T., et al. 2017, MNRAS, 467, 2039, doi: 10.1093/mnras/stx215
  * De Rosa et al. (2019) De Rosa, A., Vignali, C., Bogdanović, T., et al. 2019, New A Rev., 86, 101525, doi: 10.1016/j.newar.2020.101525
  * Deller et al. (2011) Deller, A. T., Brisken, W. F., Phillips, C. J., et al. 2011, PASP, 123, 275, doi: 10.1086/658907
  * Feissel-Vernier (2003) Feissel-Vernier, M. 2003, A&A, 403, 105, doi: 10.1051/0004-6361:20030348
  * Fey et al. (1997) Fey, A. L., Eubanks, M., & Kingham, K. A. 1997, AJ, 114, 2284, doi: 10.1086/118648
  * Frey et al. (1997) Frey, S., Gurvits, L. I., Kellermann, K. I., Schilizzi, R. T., & Pauliny-Toth, I. I. K. 1997, A&A, 325, 511
  * Frey et al. (2018) Frey, S., Titov, O., Melnikov, A. E., de Vicente, P., & Shu, F. 2018, A&A, 618, A68, doi: 10.1051/0004-6361/201832771
  * Gipson (2010) Gipson, J. 2010, in Sixth International VLBI Service for Geodesy and Astronomy. Proceedings from the 2010 General Meeting, ed. R. Navarro, S. Rogstad, C. E. Goodhart, E. Sigman, M. Soriano, D. Wang, L. A. White, & C. S. Jacobs, 77–84
  * Gordon (2010) Gordon, D. 2010, 242
  * Gordon et al. (2016) Gordon, D., Jacobs, C., Beasley, A., et al. 2016, AJ, 151, 154, doi: 10.3847/0004-6256/151/6/154
  * Greisen (2003) Greisen, E. W. 2003, in Astrophysics and Space Science Library, Vol. 285, Information Handling in Astronomy - Historical Vistas, ed. A. Heck, 109, doi: 10.1007/0-306-48080-8_7
  * Gwinn et al. (1997) Gwinn, C. R., Eubanks, T. M., Pyne, T., Birkinshaw, M., & Matsakis, D. N. 1997, ApJ, 485, 87, doi: 10.1086/304424
  * Helmboldt et al. (2007) Helmboldt, J. F., Taylor, G. B., Tremblay, S., et al. 2007, ApJ, 658, 203, doi: 10.1086/511005
  * Hobiger et al. (2006) Hobiger, T., Kondo, T., & Schuh, H. 2006, Radio Science, 41, RS1006, doi: 10.1029/2005RS003297
  * Hook & McMahon (1998) Hook, I. M., & McMahon, R. G. 1998, MNRAS, 294, L7, doi: 10.1046/j.1365-8711.1998.01368.x
  * Hook et al. (1995) Hook, I. M., McMahon, R. G., Patnaik, A. R., et al. 1995, MNRAS, 273, L63, doi: 10.1093/mnras/273.1.L63
  * Hunt et al. (2021) Hunt, L. R., Johnson, M. C., Cigan, P. J., Gordon, D., & Spitzak, J. 2021, AJ, 162, 121, doi: 10.3847/1538-3881/ac135d
  * Kardashev (1986) Kardashev, N. S. 1986, Soviet Ast., 30, 501
  * Keimpema et al. (2015) Keimpema, A., Kettenis, M. M., Pogrebenko, S. V., et al. 2015, Experimental Astronomy, 39, 259, doi: 10.1007/s10686-015-9446-1
  * Krezinger et al. (2022) Krezinger, M., Perger, K., Gabányi, K. É., et al. 2022, ApJS, 260, 49, doi: 10.3847/1538-4365/ac63b8
  * Kristian & Sachs (1966) Kristian, J., & Sachs, R. K. 1966, ApJ, 143, 379, doi: 10.1086/148522
  * Lister et al. (2011) Lister, M. L., Aller, M., Aller, H., et al. 2011, ApJ, 742, 27, doi: 10.1088/0004-637X/742/1/27
  * MacMillan et al. (2019) MacMillan, D. S., Fey, A., Gipson, J. M., et al. 2019, A&A, 630, A93, doi: 10.1051/0004-6361/201935379
  * Makarov & Secrest (2022) Makarov, V. V., & Secrest, N. J. 2022, arXiv e-prints, arXiv:2202.07536.  https://arxiv.org/abs/2202.07536
  * Malkin (2018) Malkin, Z. 2018, ApJS, 239, 20, doi: 10.3847/1538-4365/aae777
  * Mignard & Klioner (2012) Mignard, F., & Klioner, S. 2012, A&A, 547, A59, doi: 10.1051/0004-6361/201219927
  * Moór et al. (2011) Moór, A., Frey, S., Lambert, S. B., Titov, O. A., & Bakos, J. 2011, AJ, 141, 178, doi: 10.1088/0004-6256/141/6/178
  * Nothnagel et al. (2017) Nothnagel, A., Artz, T., Behrend, D., & Malkin, Z. 2017, Journal of Geodesy, 91, 711, doi: 10.1007/s00190-016-0950-5
  * O’Sullivan et al. (2011) O’Sullivan, S. P., Gabuzda, D. C., & Gurvits, L. I. 2011, MNRAS, 415, 3049, doi: 10.1111/j.1365-2966.2011.18915.x
  * Petit & Luzum (2010) Petit, G., & Luzum, B., eds. 2010, IERS Conventions 2010 (IERS Technical Note No. 36)
  * Petrov et al. (2005) Petrov, L., Kovalev, Y. Y., Fomalont, E., & Gordon, D. 2005, AJ, 129, 1163, doi: 10.1086/426920
  * Petrov et al. (2008) Petrov, L., Kovalev, Y. Y., Fomalont, E. B., & Gordon, D. 2008, AJ, 136, 580, doi: 10.1088/0004-6256/136/2/580
  * Petrov et al. (2011) —. 2011, AJ, 142, 35, doi: 10.1088/0004-6256/142/2/35
  * Punsly et al. (2021) Punsly, B., Frey, S., Reynolds, C., et al. 2021, ApJ, 919, 40, doi: 10.3847/1538-4357/ac1070
  * Pushkarev & Kovalev (2012) Pushkarev, A. B., & Kovalev, Y. Y. 2012, A&A, 544, A34, doi: 10.1051/0004-6361/201219352
  * Pyne et al. (1996) Pyne, T., Gwinn, C. R., Birkinshaw, M., Eubanks, T. M., & Matsakis, D. N. 1996, ApJ, 465, 566, doi: 10.1086/177443
  * Richards et al. (2011) Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29, doi: 10.1088/0067-0049/194/2/29
  * Romani et al. (2004) Romani, R. W., Sowards-Emmerd, D., Greenhill, L., & Michelson, P. 2004, ApJ, 610, L9, doi: 10.1086/423201
  * Saastamoinen (1972) Saastamoinen, J. 1972, in The Use of Artificial Satellites for Geodesy, Vol. 15, 247, doi: 10.1029/GM015p0247
  * Schaer et al. (1996) Schaer, S., Beutler, G., Rothacher, M., & Springer, T. A. 1996, in Proceedings of the IGS Analysis Center Workshop 1996, ed. R. Neilan, P. Van Scoy, & J. Zumberge, 181–192
  * Shepherd (1997) Shepherd, M. C. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 125, Astronomical Data Analysis Software and Systems VI, ed. G. Hunt & H. Payne, 77
  * Shuygina et al. (2019) Shuygina, N., Ivanov, D., Ipatov, A., et al. 2019, Geodesy and Geodynamics, 10, 150, doi: https://doi.org/10.1016/j.geog.2018.09.008
  * Sovers et al. (1998) Sovers, O. J., Fanselow, J. L., & Jacobs, C. S. 1998, Reviews of Modern Physics, 70, 1393, doi: 10.1103/RevModPhys.70.1393
  * Sovers & Jacobs (1996) Sovers, O. J., & Jacobs, C. S. 1996, Observation model and parameter partials for the JPL VLBI parameter estimation software MODEST, 1996, Jet Propulsion Lab. Report
  * Sowards-Emmerd et al. (2004) Sowards-Emmerd, D., Romani, R. W., Michelson, P. F., & Ulvestad, J. S. 2004, ApJ, 609, 564, doi: 10.1086/421239
  * Surkis et al. (2008) Surkis, I., Bogdanov, A., Zimovsky, V., et al. 2008, in Measuring the Future, Proceedings of the Fifth IVS, ed. A. Finkelstein & D. Behrend, 124–128
  * Titov et al. (2011) Titov, O., Lambert, S. B., & Gontier, A. M. 2011, A&A, 529, A91, doi: 10.1051/0004-6361/201015718
  * Titov & Malkin (2009) Titov, O., & Malkin, Z. 2009, A&A, 506, 1477, doi: 10.1051/0004-6361/200912369
  * Titov et al. (2017) Titov, O., Pursimo, T., Johnston, H. M., et al. 2017, AJ, 153, 157, doi: 10.3847/1538-3881/aa61fd
  * Titov & Schuh (2000) Titov, O., & Schuh, H. 2000, IERS Technical Note, 28, 33
  * Titov et al. (2013) Titov, O., Stanford, L. M., Johnston, H. M., et al. 2013, AJ, 146, 10, doi: 10.1088/0004-6256/146/1/10
  * Titov et al. (2004) Titov, O., Tesmer, V., & Boehm, J. 2004, in International VLBI Service for Geodesy and Astrometry 2004 General Meeting Proceedings, ed. N. R. Vandenberg & K. D. Baver, 267
  * Titov (2007) Titov, O. A. 2007, Astronomy Letters, 33, 481, doi: 10.1134/S1063773707070055
  * Tremaine & Davis (2014) Tremaine, S., & Davis, S. W. 2014, MNRAS, 441, 1408, doi: 10.1093/mnras/stu663
  * Veres et al. (2010) Veres, P., Frey, S., Paragi, Z., & Gurvits, L. I. 2010, A&A, 521, A6, doi: 10.1051/0004-6361/201014957
  * Whitney et al. (2004) Whitney, A. R., Cappallo, R., Aldrich, W., et al. 2004, Radio Science, 39, RS1007, doi: 10.1029/2002RS002820
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513
  * Zhang et al. (2020) Zhang, Y., An, T., & Frey, S. 2020, Science Bulletin, 65, 525, doi: 10.1016/j.scib.2020.01.008
  * Zhang et al. (2017) Zhang, Y., An, T., Frey, S., et al. 2017, MNRAS, 468, 69, doi: 10.1093/mnras/stx392
  * Zhang et al. (2021) —. 2021, MNRAS, 507, 3736, doi: 10.1093/mnras/stab2289


stakeholders in the domain of CPPS engineering, we selected engineers and
stakeholders from different companies and domains.

###  A.14 Data Collection

We used _shadowing_ [63] with X participants to investigate the subjects
during performing the eIPSE approach on the _chess piece_ use case. Due to the
distributed locations of the subjects performing the evaluation, we used Zoom
to connect them to the eIPSE toolchain and provided them with remote control.
At least two of the authors shadowed the subjects during the evaluation
sessions. One author helped the evaluation subject if questions arose,
additional explanations were required, or the subjects were stuck in the
process. The other authors present took notes for later investigation.
Furthermore, one author stopped the time for each of the activities and steps
of the evaluation process. Additionally, we recorded the evaluation sessions
to replay them later during the internal result analysis. Afterward, we asked
the subjects about the perceived usefulness of the toolchain and the lessons
learned during the process.

## References

  * Ananieva et al. [2016] Ananieva, S., Kowal, M., Thüm, T., Schaefer, I., 2016\.  Implicit constraints in partial feature models, in: 7th Int. FOSD Workshop, FOSD@SPLASH 2016, Amsterdam, Netherlands, October 30, 2016, pp. 18–27.
  * Apel et al. [2013] Apel, S., Batory, D., Kästner, C., Saake, G., 2013\.  Feature-Oriented Software Development: Concepts and Implementation.  Springer.
  * Biffl et al. [2017] Biffl, S., Gerhard, D., Lüder, A., 2017.  Introduction to the Multi-Disciplinary Engineering for Cyber-Physical Production Systems, in: Multi-Disciplinary Engineering for Cyber-Physical Production Systems. Springer, pp. 1–24.
  * Campbell et al. [1990] Campbell, G.H., Faulk, S.R., Weiss, D.M., 1990.  Introduction To Synthesis.  Technical Report. INTRO_SYNTHESIS_PROCESS-90019-N, Software Productivity Consortium, Herndon, VA, USA.
  * Cañete et al. [2022] Cañete, A., Amor, M., Fuentes, L., 2022.  Supporting iot applications deployment on edge-based infrastructures using multi-layer feature models.  Journal of Systems and Software 183, 111086.
  * Clarke et al. [2015] Clarke, D., Helvensteijn, M., Schaefer, I., 2015.  Abstract delta modelling.  Mathematical Structures in Computer Science 25, 482–527.
  * Clements and Northrop [2002] Clements, P., Northrop, L., 2002\.  Software product lines.  Addison-Wesley Boston.
  * Czarnecki et al. [2012] Czarnecki, K., Grünbacher, P., Rabiser, R., Schmid, K., Wąsowski, A., 2012.  Cool Features and Tough Decisions: A Comparison of Variability Modeling Approaches, in: 6th International Workshop on Variability Modeling of Software-Intensive Systems, ACM. pp. 173–182.
  * Dhungana et al. [2011] Dhungana, D., Grünbacher, P., Rabiser, R., 2011.  The DOPLER Meta-Tool for Decision-Oriented Variability Modeling: A Multiple Case Study.  Automated Software Engineering 18, 77–114.
  * Drath [2009] Drath, R., 2009.  Datenaustausch in der Anlagenplanung mit AutomationML: Integration von CAEX, PLCopen XML und COLLADA.  Springer-Verlag.
  * Drath [2021] Drath, R. (Ed.), 2021.  AutomationML - A Practical Guide.  De Gruyter Oldenbourg.
  * Drath et al. [2008] Drath, R., Luder, A., Peschke, J., Hundt, L., 2008\.  Automationml - the glue for seamless automation engineering, in: 2008 IEEE International Conference on Emerging Technologies and Factory Automation, pp. 616–623.
  * Fadhlillah et al. [2022a] Fadhlillah, H.S., Feichtinger, K., Bauer, P., Kutsia, E., Rabiser, R., 2022a.  V4rdiac: Tooling for multidisciplinary delta-oriented variability management in cyber-physical production systems, Association for Computing Machinery, New York, NY, USA. p. 34–37.
  * Fadhlillah et al. [2022b] Fadhlillah, H.S., Feichtinger, K., Meixner, K., Sonnleithner, L., Rabiser, R., Zoitl, A., 2022b.  Towards multidisciplinary delta-oriented variability management in cyber-physical production systems, in: Proc. of the 16th International Working Conference on Variability Modelling of Software-Intensive Systems (VaMoS), ACM.
  * Fadhlillah et al. [2023a] Fadhlillah, H.S., Fernández, A.M.G., Rabiser, R., Zoitl, A., 2023a.  Managing cyber-physical production systems variability using v4rdiac: Industrial experiences, in: Proceedings of the 27th ACM International Systems and Software Product Line Conference - Volume A, Association for Computing Machinery, New York, NY, USA. p. 223–233.
  * Fadhlillah et al. [2023b] Fadhlillah, H.S., Sharma, S., Gutierrez Fernandez, A.M., Rabiser, R., Zoitl, A., 2023b.  Delta modeling in iec 61499: Expressing control software variability in cyber-physical production systems, in: 2023 IEEE 28th International Conference on Emerging Technologies and Factory Automation (ETFA), pp. 1–8.  doi:10.1109/ETFA54631.2023.10275693.
  * Fang [2019] Fang, M., 2019.  Model-Based Software Derivation for Industrial Automation Management Systems.  Ph.D. thesis. Technische Universität Kaiserslautern.
  * Fang et al. [2013] Fang, M., Leyh, G., Elsner, C., Dörr, J., 2013\.  Challenges in managing behavior variability of production control software, in: PLEASE@ICSE, IEEE Computer Society. pp. 21–24.
  * Feichtinger et al. [2022a] Feichtinger, K., Meixner, K., Biffl, S., Rabiser, R., 2022a.  Evolution support for custom variability artifacts using feature models: A study in the cyber-physical production systems domain, in: Perrouin, G., Moha, N., Seriai, A.D. (Eds.), Reuse and Software Quality, Springer International Publishing, Cham. pp. 79–84.
  * Feichtinger et al. [2020] Feichtinger, K., Meixner, K., Rabiser, R., Biffl, S., 2020\.  Variability Transformation from Industrial Engineering Artifacts: An Example in the Cyber-Physical Production Systems Domain, in: 3rd International Workshop on Variability and Evolution of Software-Intensive Systems (VariVolution), SPLC ’20: 24th ACM International Systems and Software Product Line Conference, Volume B, ACM. pp. 65–73.
  * Feichtinger et al. [2022b] Feichtinger, K., Meixner, K., Rinker, F., Koren, I., Eichelberger, H., Heinemann, T., Holtmann, J., Konersmann, M., Michael, J., Neumann, E.M., Pfeiffer, J., Rabiser, R., Riebisch, M., Schmid, K., 2022b.  Industry voices on software engineering challenges in cyber-physical production systems engineering, in: 2022 IEEE 27th International Conference on Emerging Technologies and Factory Automation (ETFA), pp. 1–8.
  * Feichtinger and Rabiser [2020a] Feichtinger, K., Rabiser, R., 2020a.  Towards transforming variability models: Usage scenarios, required capabilities and challenges, in: 24th ACM International Systems and Software Product Line Conference - Volume B, ACM, New York, NY, USA. p. 44–51.
  * Feichtinger and Rabiser [2020b] Feichtinger, K., Rabiser, R., 2020b.  Variability model transformations: Towards unifying variability modeling, in: 46th Euromicro Conference on Software Engineering and Advanced Applications, IEEE, Portoroz, Slovenia.
  * Feichtinger and Rabiser [2021] Feichtinger, K., Rabiser, R., 2021\.  How flexible must a transformation approach for variability models and custom variability representations be?, in: 4rd International Workshop on Languages for Modelling Variability (MODEVAR), co-located with SPLC 2021, ACM. pp. 69–72.
  * Feichtinger et al. [2021] Feichtinger, K., Stöbich, J., Romano, D., Rabiser, R., 2021\.  Travart: An approach for transforming variability models, in: 15th International Working Conference on Variability Modelling of Software-Intensive Systems, ACM. pp. 8:1–8:10.
  * Feldmann et al. [2015] Feldmann, S., Legat, C., Vogel-Heuser, B., 2015.  Engineering support in the machine manufacturing domain through interdisciplinary product lines: An applicability analysis.  IFAC-PapersOnLine 28, 211–218.
  * Galster et al. [2013] Galster, M., Weyns, D., Tofan, D., Michalik, B., Avgeriou, P., 2013.  Variability in software systems-a systematic literature review.  IEEE Transactions on Software Engineering 40, 282–306.
  * Geraldi et al. [2020] Geraldi, R.T., Reinehr, S., Malucelli, A., 2020.  Software product line applied to the internet of things: A systematic literature review.  Information and Software Technology 124, 106293.
  * Gunes et al. [2014] Gunes, V., Peter, S., Givargis, T., Vahid, F., 2014\.  A survey on concepts, applications, and challenges in cyber-physical systems.  KSII Transactions on Internet & Information Systems 8.
  * Hevner [2007] Hevner, A.R., 2007.  A three cycle view of design science research.  Scandinavian journal of information systems 19, 4.
  * Hevner et al. [2008] Hevner, A.R., March, S.T., Park, J., Ram, S., 2008.  Design science in information systems research.  Management Information Systems Quarterly 28, 6.
  * Hubaux et al. [2012] Hubaux, A., Xiong, Y., Czarnecki, K., 2012.  A user survey of configuration challenges in linux and ecos, in: Proceedings of the 6th International Workshop on Variability Modeling of Software-Intensive Systems, Association for Computing Machinery, New York, NY, USA. p. 149–155.
  * International Electrotechnical Commission (2014) [IEC] International Electrotechnical Commission (IEC), 2014\.  IEC 62714-1, Engineering data exchange format for use in industrial automation systems engineering - Automation markup language \- Part 1: Architecture and general requirements.
  * International Electrotechnical Commission , TC65/WG6(2012) [IEC] International Electrotechnical Commission (IEC), TC65/WG6, 2012\.  IEC 61499-1, Function Blocks - part 1: Architecture: Edition 2.0.
  * Jazdi et al. [2010] Jazdi, N., Maga, C., Göhner, P., Ehben, T., Tetzner, T., Löwen, U., 2010\.  Mehr Systematik für den Anlagenbau und das industrielle Lösungsgeschäft - Gesteigerte Effizienz durch Domain Engineering 58, 524–532.  doi:doi:10.1524/auto.2010.0867.
  * Järvenpää et al. [2019] Järvenpää, E., Siltala, N., Hylli, O., Lanz, M., 2019\.  Implementation of capability matchmaking software facilitating faster production system design and reconfiguration planning.  Journal of Manufacturing Systems 53, 261–270.
  * Kang et al. [1990] Kang, K.C., Cohen, S.G., Hess, J.A., Novak, W.E., Peterson, A.S., 1990.  Feature-oriented domain analysis (FODA) feasibility study.  Technical Report. Carnegie-Mellon Univ., Pittsburgh, Pa, Software Engineering Inst.
  * Kästner et al. [2012] Kästner, C., Ostermann, K., Erdweg, S., 2012.  A variability-aware module system, in: Proc. of the ACM Int’l Conf. on Object Oriented Programming Systems Languages and Applications, ACM, New York, NY, USA. p. 773–792.
  * Kowal et al. [2017] Kowal, M., Ananieva, S., Thüm, T., Schaefer, I., 2017\.  Supporting the Development of Interdisciplinary Product Lines in the Manufacturing Domain.  IFAC-PapersOnLine 50, 4336–4341.
  * Krüger et al. [2019] Krüger, J., Mukelabai, M., Gu, W., Shen, H., Hebig, R., Berger, T., 2019.  Where is my feature and what is it about? A case study on recovering feature facets.  Journal of Systems and Software 152, 239–253.
  * Krüger et al. [2017] Krüger, J., Nielebock, S., Krieter, S., Diedrich, C., Leich, T., Saake, G., Zug, S., Ortmeier, F., 2017\.  Beyond Software Product Lines: Variability Modeling in Cyber-Physical Systems, in: 21st International Systems and Software Product Line Conference, SPLC 2017, Volume A, Sevilla, Spain, September 25-29, 2017, ACM, New York, NY, USA. pp. 237–241.
  * Lee [1989] Lee, S., 1989.  Disassembly planning based on subassembly extraction, in: Third ORSA/TIMS Conference on Flexible Manufacturing System, pp. 383–388.
  * Meinicke et al. [2017] Meinicke, J., Thüm, T., Schröter, R., Benduhn, F., Leich, T., Saake, G., 2017\.  Mastering Software Variability with FeatureIDE.  Springer.
  * Meixner [2020] Meixner, K., 2020.  Integrating Variability Modeling of Products, Processes, and Resources in Cyber-Physical Production Systems Engineering, in: 24th ACM International Systems and Software Product Line Conference - Volume B, ACM. pp. 96–103.
  * Meixner et al. [2021a] Meixner, K., Feichtinger, K., Rabiser, R., Biffl, S., 2021a.  A reusable set of real-world product line case studies for comparing variability models in research and practice, in: 25th International Systems and Software Product Line Conference - Volume B, ACM. pp. 105–112.
  * Meixner et al. [2022] Meixner, K., Feichtinger, K., Rabiser, R., Biffl, S., 2022\.  Efficient Production Process Variability Exploration, in: Arcaini, P., Devroey, X., Fantechi, A. (Eds.), VaMoS ’22: 16th International Working Conference on Variability Modelling of Software-Intensive Systems, Florence, Italy, February 23 - 25, 2022, ACM. pp. 14:1–14:9.
  * Meixner et al. [2019] Meixner, K., Rabiser, R., Biffl, S., 2019.  Towards modeling variability of products, processes and resources in cyber-physical production systems engineering, in: 23rd International Systems and Software Product Line Conference - Volume B, ACM. pp. 68:1–68:8.
  * Meixner et al. [2020] Meixner, K., Rabiser, R., Biffl, S., 2020.  Feature Identification for Engineering Model Variants in Cyber-Physical Production Systems Engineering, in: 14th International Working Conference on Variability Modelling of Software-Intensive Systems, ACM. pp. 18:1–18:5.
  * Meixner et al. [2021b] Meixner, K., Rinker, F., Marcher, H., Decker, J., Biffl, S., 2021b.  A Domain-Specific Language for Product-Process-Resource Modeling, in: 26th IEEE International Conference on Emerging Technologies and Factory Automation, ETFA 2021, Vasteras, Sweden, September 7-10, 2021, IEEE. pp. 1–8.
  * Monostori [2014] Monostori, L., 2014.  Cyber-physical Production Systems: Roots, Expectations and R&D Challenges.  Procedia CIRP 17, 9–13.
  * Nyberg [2021] Nyberg, M., 2021.  Generating safety cases for large-scale industrial product lines.  URL: https://splc2021.net/program/keynotes. keynote at 25th ACM International Systems and Software Product Line Conference.
  * Oliinyk et al. [2017] Oliinyk, O., Petersen, K., Schoelzke, M., Becker, M., Schneickert, S., 2017.  Structuring automotive product lines and feature models: an exploratory study at opel.  Requirements Engineering 22, 105–135.
  * Paetzold [2017] Paetzold, K., 2017.  Product and systems engineering/ca* tool chains, in: Multi-Disciplinary Engineering for Cyber-Physical Production Systems. Springer, pp. 27–62.
  * Rabiser et al. [2018] Rabiser, D., Prähofer, H., Grünbacher, P., Petruzelka, M., Eder, K., Angerer, F., Kromoser, M., Grimmer, A., 2018\.  Multi-purpose, multi-level feature modeling of large-scale industrial software systems.  Software and Systems Modeling 17, 913–938.
  * Runeson et al. [2012] Runeson, P., Host, M., Rainer, A., Regnell, B., 2012\.  Case study research in software engineering: Guidelines and examples.  John Wiley & Sons.
  * Safdar et al. [2021] Safdar, S.A., Lu, H., Yue, T., Ali, S., Nie, K., 2021\.  A framework for automated multi-stage and multi-step product configuration of cyber-physical systems 20, 211–265.
  * Safdar et al. [2016] Safdar, S.A., Yue, T., Ali, S., Lu, H., 2016.  Evaluating variability modeling techniques for supporting cyber-physical system product line engineering.  Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 9959 LNCS, 1–19.
  * Schäfer et al. [2021] Schäfer, A., Becker, M., Andres, M., Kistenfeger, T., Rohlf, F., 2021.  Variability realization in model-based system engineering using software product line techniques: An industrial perspective, in: Proceedings of the 25th ACM International Systems and Software Product Line Conference - Volume A, Association for Computing Machinery, New York, NY, USA. p. 25–34.
  * Schleipen et al. [2015] Schleipen, M., Lüder, A., Sauer, O., Flatt, H., Jasperneite, J., 2015.  Requirements and concept for plug-and-work.  at-Automatisierungstechnik 63, 801–820.
  * Schmid and John [2004] Schmid, K., John, I., 2004.  A customizable approach to full lifecycle variability management.  Sci. Comput. Program. 53, 259–284.
  * Schmid et al. [2011] Schmid, K., Rabiser, R., Grünbacher, P., 2011.  A comparison of decision modeling approaches in product lines, in: 5th International Workshop on Variability Modelling of Software-Intensive Systems, ACM. pp. 119–126.
  * Shull et al. [2007] Shull, F., Singer, J., Sjøberg, D.I., 2007.  Guide to Advanced Empirical Software Engineering.  Springer.
  * Singer et al. [2008] Singer, J., Sim, S.E., Lethbridge, T.C., 2008.  Software engineering data collection for field studies.  Guide to advanced empirical software engineering , 9–34.
  * Sundermann et al. [2021] Sundermann, C., Feichtinger, K., Engelhardt, D., Rabiser, R., Thüm, T., 2021.  Yet another textual variability language? a community effort towards a unified language, in: Proceedings of the 25th ACM International Systems and Software Product Line Conference - Volume A, Association for Computing Machinery, New York, NY, USA. p. 136–147.
  * Tolio et al. [2010] Tolio, T., Ceglarek, D., ElMaraghy, H.A., Fischer, A., Hu, S.J., Laperrière, L., Newman, S.T., Váncza, J., 2010\.  Species—co-evolution of products, processes and production systems.  CIRP Annals 59, 672–693.
  * [66] VDI/VDE 3682, 2005.  VDI/VDE 3682: Formalised Process Descriptions.  Beuth Verlag.
  * [67] VDI/VDE 3695, 20010-2013.  VDI/VDE 3695: Engineering of industrial plants.  Beuth Verlag, 5 parts.
  * Wieringa [2014] Wieringa, R.J., 2014.  Design science methodology for information systems and software engineering.  Springer.
  * Wohlin et al. [2012] Wohlin, C., Runeson, P., Höst, M., Ohlsson, M.C., Regnell, B., Wesslén, A., 2012\.  Experimentation in software engineering.  Springer Science & Business Media.
  * Zhang et al. [2016] Zhang, B., Duszynski, S., Becker, M., 2016.  Variability mechanisms and lessons learned in practice, in: Proceedings of the 1st International Workshop on Variability and Complexity in Software Design, Association for Computing Machinery, New York, NY, USA. p. 14–20.
  * Zoitl et al. [2010] Zoitl, A., Strasser, T., Valentini, A., 2010.  Open source initiatives as basis for the establishment of new technologies in industrial automation: 4diac a case study, in: 2010 IEEE Int’l Symp. on Industrial Electronics, IEEE. pp. 3817–3819.