Resolving VLBI correlator ambiguity in the time delay model improves precision of geodetic measurements

The modern Very Long Baseline Interferometry (VLBI) relativistic delay model, as documented in the IERS Conventions refers to the time epoch when the signal passes one of two stations of an interferometer baseline (selected arbitrarily from the pair of stations and called the 'reference station', or 'station 1'). This model consists of the previous correlation procedure used before the year 2002. However, since 2002 a new correlation procedure that produces the VLBI group delays referring to the time epoch of signal passage at the geocenter has been used. A corresponding correction to the conventional VLBI model delay has to be introduced. However, this correction has not been thoroughly presented in peer reviewed journals, and different approaches are used at the correlators to calculate the final group delays officially published in the IVS database. This may cause an inconsistency up to 6 ps for ground-based VLBI experiments between the group delay obtained by the correlator and the geometrical model delay from the IERS Conventions used in data analysis software. Moreover, a miscalculation of the signal arrival moment to the 'reference station' could result a larger modelling error (up to 50 ps). The paper presents the justification of the correction due to transition between two epochs elaborated from the Lorentz transformation, and the approach to model the uncertainty of the calculation of the signal arrival moment. The both changes are particularly essential for upcoming broadband technology geodetic VLBI observations.


INTRODUCTION
The Very Long Baseline Interferometry (VLBI) technique measures the difference between the arrival times of a signal from a distant radio source at two radio telescopes (Schuh & Behrend (2012)). The signal is recorded at each radio telescope together with time marks from independent hydrogen masers. Due to separation of the radio telescopes by a few hundred or thousand kilometres, the plain wave front passes first telescope earlier then the second one. This difference in the arrival time of the signal at both radio telescopes is known as time delay, and the frequency shift due to the relative motion of the telescopes around the geocentre is known as delay rate.
The time delay and delay rate are found during crosscorrelation of the two independent records. There are two types of correlators (XF and FX) based on the order of the mathematical operations -cross-correlation (X) and Fourier transformation (F). Baseline-based correlators are designed as XF type correlators, and station-based correlators are FX type correlators. For the baselinebased XF-type MarkIII correlator used before 2002, the observables referred to the position of one of the two stations (station 1). For the station-based FX-type MarkIV correlator all observables for all baselines at one single multi-baseline scan are referred to the geocentre as a common reference point. As 1 ps precision is required for the time delay calculation, all first-order and secondorder effects of special relativity should be taken into account.
One of the goals of the International VLBI Service activities is to achieve 1-mm accuracy from the analysis of routine geodetic VLBI sessions. The accuracy of the daily scale factor improved dramatically in 2002 when the MarkIII correlator was replaced by MarkIV correlator (Titov & Krásná (2018)). However, so far this value varies about 3-4 mm despite technological developments since 2002.
One possible reason for the lack of improvement in accuracy is the inconsistency between the VLBI observable group delays and the relativistic delay model developed in 1980s-90s, published in the IERS Conventions 2010 (Petit & Luzum (2010)). The transition from the MarkIII to MarkIV correlator was not followed by any changes in the IERS Conventions model that still refers to the epoch of the wavefront passage of station 1. Thus, it remains consistent with the XF-type correlators. To make the output delay of the FX-correlator consistent with the IERS Conventions 2010 model (XF-type), an additional correction needs to be applied. Unfortunately, this correction has not been officially presented in explicit form. This conversion difference was called "subtle" (Whitney (2000)); however, it reaches 20 ns, which is quite significant. Corey (2000) developed a simple geometric approach under assumption of the finiteness of the speed of light to obtain this correction, but his final equation comprised a major term only, while several minor terms were not included.
In this paper, it is emphasised that the relativistic correction due to the change of the reference epoch definition should be derived from the Lorenz transformation to secure the 1-ps accuracy. Therefore, the final equation of the recommended group delay should include some minor relativistic terms due to coupling of the barycentric and geocentric velocities of the radiotelescopes to be added to the version by Corey (2000). A detailed development of the correction based on the Lorenz transformation is given in Appendix. This correction is essential from the theoretical point of view; however, its impact on the geodetic results is negligible for ground-based baselines (less than 1 mm).
A more serious problem is caused by the uncertainty in the signal arrival time as calculated by the correlators, even if the problem of the epoch calculation is fixed. Within the adopted procedure, for a single multi-station scan this time is common for all stations and is usually rounded to integer number of seconds. Meanwhile, for a multi-station scan, the factual signal arrival time is individual for each station, the output group delay is converted to the time common for all stations within one scan using a reasonable polynomial approximation. Therefore, the final output delay for each baseline is referred to the common time of scan. Theoretically, this output delay should be perfectly consistent to the delay at the time of the signal arrival to the "reference station" of each baseline. However, this is not guaranteed due to the uncertainty of the reference epoch definition (discussed in the Appendix) and hidden numerical issues during the polynomial approximation.
To estimate an additional correction, the standard parametrical model should be extended. For each scan we have a time of the signal arrival (common for N stations) and a set of N (N − 1)/2 time delays for all baselines. Instead of seeking for N (N − 1)/2 errors in the delays themselves, it would be easier to treat the signal arrival time as the parameter to be updated assuming the delays are errorless. A possible approach to model this type inconsistency is presented analytically in Finkelstein, Kreinovich, & Pandey (1983). A second order term in Equation (18) may be generalised at 1-ps accuracy in the form The case = 0 corresponds to the selection of the reference clock at station 1, and the case = ∞ corresponds to the selection of the reference clock at station 2. The relativistic group delay model from the IERS Conventions has an intrinsic assumption that = 0. A violation of this assumption results in a small deviation of the from zero. For a small value of it could be parametrized with the partial derivative By its analytical representation, this new parameter should be referred to the group of parameters to model the clock instability (offset, rate, 2nd derivative, etc). In total, (N − 1) parameters should be added to the traditional procedure of the VLBI delay modelling.
The parameters could be estimated with Equation (2) by the least squares individually for each VLBI station clock except to the clock at the network "reference station" that is assumed to be errorless. Then for two arbitrary stations (i and j) the corresponding delay is calculated as follows

DATA ANALYSIS
The second term in Equation (18) (of the Appendix) is the diurnal variation of the Earth scale's factor that replaces the diurnal aberration applied for the traditional astronomical observations. This is the only term due to the Earth's rotation implemented by the FX-correlator software developers (in accordance to Corey (2000)). However, a more accurate approach based on the Lorenz transformation (12) reveals additional minor terms in Equation (18) due to coupling of the two velocities V and w 2 . The first term in Equation (18) is the coordinate term due to the transformation from the barycentric to the geocentric reference frame, and it could be ignored for the scope of this paper. We used one of the recent VGOS experiments (VT9290, 17-Oct-2019) for more detailed analysis. This 24-hour experiment included five radio telescopes (WETT13S, ONSA13NE, ONSA13SW, GGAO12M and KOKEE12M) equipped with the broad band receivers.  Standard geodetic VLBI observations operated in two frequencies, 2.3 GHz (S-band) and 8.4 GHz (X-band), are not sensitive to the effect of the time of signal arrival. Therefore, we used the new broad band VLBI data (VGOS project) to estimate the parameter . Due to the higher sample rate and the broader bandwidth of the recorded data the formal accuracy of the VGOS geodetic results is better than for standard S/X observations by an order of magnitude.
The VGOS data files were processed using the OC-CAM software (Titov, Tesmer, & Boehm (2004)) (version 6.3) in two modes. A first solution produces a standard set of parameters for estimating -(i) corrections to the positions of radio telescopes in the ITRF2014 frame (Altamimi et al. (2016)), (ii) Earth orientation parameters, (iii) wet troposphere delay and two gradients, (iv) three parameters to model the clock instability for each station except the reference one (clock offset, clock rate and second derivative), and (v) corrections to the ICRF3 positions of several radio sources that expose a high level of astrometric instability in the past. A second solution was used to estimate the parameter for all stations except for the reference one.
Estimates of the parameter for six VGOS stations operating during 2019 are shown in Table 1. About half of the estimates are statistically significant. This means that typically, the time of the radio wave arrival to the reference station is not calculated by the correlator with sufficient accuracy. The resulting group delay calculated by Equation (3) for baseline GGAO12M -KOKEE12M at the same session (17-Oct-2019, MJD = 58744) is shown in Fig 2. We selected this baseline because for both stations in this experiment the estimates of are larger than usual (−0.885 · 10 −3 for GGAO12M and 0.892 · 10 −3 for KOKEE12M). The range of the peak-to-peak variations is about 80 ps. This results in additional, hidden, source of systematic error for all other parameters.

ANALYSIS OF ASTROMETRIC RESULTS
A comprehensive analysis of geodetic parameters is beyond of the scope of this paper. Herewith we discuss only effect of the additional parameter on the astrometric positions of two well-known reference radio sources, namely 0552+398 and 1156+295. Both sources were observed in twenty 24-h broadband VLBI experiments during 2019, with a large number of observations. As a result, their formal positional errors for both components are less than 50 µas for almost all experiments. Therefore, statistical investigation of the astrometrical results would demonstrate the advantage of the new VLBI technology application and the effect of the inclusion of the additional modelling parameter.

Radio source 0552+398
Radio source 0552+398 is one of the most actively observed radio sources by geodetic VLBI since 1979 due to its strong flux density and good astrometric stability. It was included to the list of reference radio sources of ICRF1, (Ma et al. (1998)), ICRF2 (Fey et al. (2015)) and ICRF3 (Charlot et al. (2020)). It was also treated as a 'stable' one after independent verification by Feissel-Vernier (2003). The source 0552+398 has no apparent structure at S-and X-bands images. However, its imaging at higher frequencies (24 and 43 GHz) discloses a sub-milliarcsec jet in the east direction from the core (Charlot et al. (2010)). Recently, a second component was revealed by Bolotin et al. (2019) from the analysis of the broadband observations during the CONT17 campaign. While the daily estimates of the corrections to the declination component in Fig 3 vary around the original ICRF3 catalogue position within 0.6 mas, the estimates of the correction to right ascension (RA = 05h 55m 30s.80561207) (Charlot et al. (2020)) show a non-zero offset of approximately 0.2 mas. The available post-ICRF3 catalogues (e.g. the celestial reference frame solution aus2020a published by International VLBI Service (IVS)) including S/X observations during 2019-2020 do not detect any essential offset of the 0552+398 positions with respect to the ICRF3 catalogue coordinates. This potentially indicates that the jet observed at high frequencies (24 and 43 GHz) is also essential for frequencies between 2 and 11 GHz, even though it is not detected on the S/X images. We conclude that the broadband VLBI observations are more sensitive to the sub-milliarcsec structure than the traditional S/X VLBI observations as also hinted by Bolotin et al. (2019). Therefore, the positions of the reference radio sources observed by the new broadband technology are not necessary to be coincided with the S/X data positions.

Radio source 1156+295
Radio source 1156+295 has actively monitored for the last 30 years over a wide range of frequencies. Despite its extended structure in S-and X-bands with an elongated jet in the north direction (e.g. Kellermann et al. (1998)), the radio source 1156+295 demonstrates a moderate range of astrometric instability. At the same time, no structure was reported in 24 GHz and 43 GHz (Charlot et al. (2010)). Therefore, it was selected as one of the defining reference sources in the second ICRF realization (ICRF2) (Fey et al. (2015)), although, not included to the list of the ICRF3 reference sources. Our analysis of the broadband VLBI results highlights a higher range of astrometric instability in declination than right ascension time series (Fig 5) during 2019, presumably, induced by the jet oriented to the north direction. The average declination component is shifted approximately 0.2 mas south with respect to the ICRF3 catalogue position.
The difference between the two sets of daily estimates in Fig 4 and 6 does not reveal any noticeable astrometric signature due to inclusion of to the list of estimated parameters. For both sources the peak-to-peak variations do not exceed 0.25 mas in both components. Therefore, for radio sources 0552+398 and 1156+295, the inclusion of new parameter does not change the source position estimates essentially. However, for rare observed radio sources this difference may cause a substantial change in the final catalogue positions.

DISCUSSION AND CONCLUSION
The transition from the XF-type to FX-type correlators for processing geodetic VLBI data requires a corresponding revision of the relativistic group delay in the IERS Conventions to secure a match between the correlator output and the theoretical model. Alternatively, a special correction needs to be done at the final step of the post-correlation data processing. In Equation (18) we show in the four last terms the relativistic correction due to the time transformation from the epoch of the geocenter to the epoch of station 1. This correction is derived from the modified version of the Lorenz transformation in Equation (12). Missing of the three minor terms in Equation (18) can lead to a discrepancy of the group delay model at a level of 6 ps for long baselines. This is, in particular, pertinent for the intensive experiments for rapid estimation of Universal Time because a typical observational network consists of 2 or 3 radio telescopes separated by a long baseline (> 7000 km). We would like to recommend this equation be applied for the post-processing analysis of VLBI data at the modern FX-correlators.
Another effect, though may not be directly linked to the first one, is the uncertainty of the time of signal registration for each telescope as measured by the local clock (hydrogen maser) at the reference station and extrapolated during the process of correlation. This effect also refers to the difference of the geocentric velocities of the both radio telescopes, but it could be introduced as the extension of the clock instability model. The additional parameter describes how far the actual time of the signal arrival deviates from the time presented in the VLBI data file. Our analysis of broadband VLBI data over 2019 reveals that the parameter is statistically significant in many cases (Table 1). The corresponding systematic effect is up to 100 ps in time delays, and up to 0.25 mas in estimates of daily radio source positions.
It is not yet clear whether the source structure effect directly links to the problem of precisely determining the time of the signal arrival to the radio telescopes. The algorithm of the numerical calculation of the signal arrival time always relies on the assumption that the phase reference point of the target source is the same for all frequency bands. However, with a broadband receiver, we may have four different phase reference points at the four frequency bands. Therefore, four signals in each band may arrive to the receiver at four different times even from a point-like radio source. A standard calibration may not compensate this inconsistency perfectly, mostly due to the non-linear behaviour of the phase during the fringe-fitting process. In addition, an extended radio source may have four different phase reference points at four frequencies referring to the celestial reference frame. Thus, the actual differences between the signal arrival times for four frequency bands could change unpredictably. As a result, the signal arrival time presented in the broadband VLBI data file as a single value has some level of uncertainty making the additional parameter feasible for routine application using Equation (3). While it was not essential for the traditional S/X VLBI observations, the broadband VLBI observations are more accurate, and, more advanced parametrical model should be used to match these observations.

ACKNOWLEDGEMENTS
We are thankful to Sergei Kopeikin (University of Missouri-Columbia), Slava Turyshev (JPL), James Anderson (GFZ), and Igor Surkis (IAA RAS) for fruitful discussions on the theoretical aspects and the technical details of the correlation process. Also we thank the PASA Editor-in-Chief, and the anonymous referee for their constructive comments and suggestions which have significantly improved the clarity of the paper. This paper is published with the permission of the CEO, Geoscience Australia.

A.1 THE CONVENTIONAL GEOCENTRIC DELAY MODEL
The equation for the relativistic group delay model has been developed in the 1980s-90s (e.g. Hellings (1986), Kopeikin (1990), Klioner (1991), Soffel et al. (1991) to approximate the observed VLBI data at the 1-ps level of accuracy. The conventional group delay model was finally adopted Petit & Luzum (2010) where b is the vector of baseline b = r 2 − r 1 , s is the barycentric unit vector of the radio source, V is the barycentric velocity of the geocenter, w 2 is the geocentric velocity of station 2, c is the speed of light, G is the gravitational constant, M is the mass of the Sun, R is the geocentric distance to the Sun, and (·) is the dotproduct operator of two vectors. The reference epoch is the UTC epoch of the wavefront passage at the reference station. In accordance with the assumption, station 1 is treated as the reference station, the geocentric velocity of station 2 is presented in (1) explicitly. A modern revision (e.g. Soffel, Kopeikin, & Han (2017)) is to add some smaller terms (less than 1 ps), but the analytical model (1) is still valid for the analysis of VLBI data.

A.2 LORENTZ TRANSFORMATION
The radio signal is received by two radio telescopes on the surface of the rotating Earth, and their coordinates are presented in the Geocentric Celestial Reference System (GCRS) comoving with the Earth. Positions of reference radio sources emitting the signals are in the Barycentric Celestial Reference System (BCRS). So, a detailed transformation of the coordinates from BCRS to GCRS is traditionally based on the metric tensor of the Solar System at the first and second post-Newtonian level (e.g. Hellings (1986), Kopeikin (1990), Klioner (1991), Soffel, Kopeikin, & Han (2017)). However, many lower order effects are not observable, therefore, a simplified approach could be developed for the relativistic model delay.
The conventional Lorenz transformation is given by is the Lorentz "gammafactor" (Mansouri & Sexl (1977), Will (1992)). It should be noted that this factor here is not the parameter γ of the Post-Newtonian formalism (PPN) used in general relativity Will (1971). Transformation (5) links the geocentric reference system S (x , t ) that is moving with velocity V around the Solar System Barycentre (SSB) and the barycentric reference system S(x, t) located at the SSB. It could be shown (Titov & Krásná (2019)) that the time delay derived from (5) may be presented in the form Whether an astronomical instrument with a reference clock were placed in the Earth's geocenter and the Solar gravitation were ignored, the equation (6) would be applied to reduction of the geodetic VLBI data. However, further complications will be discussed in two next subsections.

A.2.1 SPACE AND TIME TRANSFORMATION
INCLUDING GRAVITATIONAL POTENTIAL The relativistic model (6) does not include the term proportional to the Solar gravitational potential 2U c 2 , where U = GM R , and few terms with the geocentric velocity w 2 presented in (4). Hellings (1986) showed that the former term appears due to the Solar gravitational field (in the Schwarzschild metric) at the Earth geocentre. Therefore, Hellings (1986) has developed new equations for the relationships between intervals of physical distance and time, measured in a moving reference geocentric frame, and the intervals, given in the barycentric coordinate system including the gravitational field of the Sun is given by Transformation (7) reduces to the Lorenz transformation (5) if the Solar potential U = 0.
The corresponding equation for the relativistic group delay includes the Solar gravitational potential at the geocenter of the Earth.
Titov & Girdiuk (2015) showed that the term proportional to 2U c 2 in (8) could be unified with the general relativity effect of the gravitational delay. Therefore, we will not include it into further analysis; however, we discuss it here as it is a part of the conventional geometric part of the relativistic delay model (Petit & Luzum (2010)).

A.2.2 LORENZ TRANSFORMATION
REFERRING TO THE EPOCH OF FIRST STATION Physical clocks (hydrogen masers) used for VLBI observations are located at the Earth surface rather than at the geocenter. As two clocks separated by a long baseline are involved for a routine observational experiment, one of them should be selected as "reference" clock. This choice is completely arbitrary, though, once it is made, the geocentric velocity of the second ("no reference") clock appears explicitly in the analytical equations. The standard approach is to consider a difference between barycentric coordinates of two radio telescopes, r 1 (t 1 ) and r 2 (t 2 ), measured at the two epochs t 1 and t 2 , to expand the vector r 2 (t 2 ) as follows r 2 (t 2 ) = r 2 (t 1 ) + w 2 (t 1 )(t 2 − t 1 ), where w 2 = w 2 (t 1 ) is the geocentric velocity of the second station at epoch t 1 . Denoting B(t 1 ) a difference between two barycentric vectors at the same epoch B = B(t 1 ) = r 2 (t 1 ) − r 1 (t 1 ) one could get for the time difference (t 2 − t 1 ) It should be noted here that B is a formal threecomponent vector rather than a meaningful physical value, though it links to the physical distance between two terrestrial positions of radio telescopes on the Earth at t 1 .
Eq (10) could be obtained by alternative way. Let's introduce of a new geocentric reference frame S = S (x , t ) with the reference epoch referred to station 1 in a such way that two geocentric reference frames S and S are linked by new transformation Transformation (11) could be easily combined with the Lorentz transformation (5) It is obvious that the transformations (11) and (12) are pertinent only for an individual pair of two radio telescopes equipped with their own high precision clocks, one of which is a reference clock and the second clock is moving with instantaneous velocity w 2 . The transformation (12) is fully consistent with the special relativity postulates and reflects the situation when the position of the reference clock is not at the reference frame origin (geocentre). For a classical astronomic instrument the reference frame origin and position of the reference clock are referred to the same topocentric position of the instrument on the Earth surface. In this scenario, the geocentric velocity of the instrument is simply added to the barycentric velocity in the formulae of the Lorenz transformation, i.e. the velocity V is replaced by the sum V + w 2 in (5) followed by a substantial change in (6). From the observational point of view this results in the appearance of the classical diurnal aberration effect. In geodetic VLBI, there is no the diurnal aberration effect at all. Instead of that, as it will be shown later, the geocentric velocity w 2 contributes to the diurnal variation of the scale factor with magnitude up to 20 ns (or 6 meters in the linear scale) for a standard baseline of 6000 km in length and a geocentric velocity of 300 m/s. For calculating the group delay from (12), one needs to develop the corresponding velocity transformation. As both reference frames S and S are geocentric, the time component is only changed due to transition from (5) to (12). Traditionally, authors proceed to the equation of the relativistic time delay (6) consistent with the XF-type correlator directly (e.g. Hellings 1986, Kopeikin 1990, Soffel et al 2017. Therefore, these two transformations (5) and (11) merge together and the difference between the delays (4) and (6) is lost. However, for the FX-type correlators, this procedure must be separated into two steps to provide a proper relativistic conversion between observables produced by the XF and FX correlators.
To elaborate equation (4) (without the 2U c 2 term) from transformation (12), let's consider the velocity transfor- (14) Now apply for a standard transition to the radio source vector and, after reduction of negligible terms, This equation could be converted to the form consistent with the conventional group delay model at 1-ps level after inclusion of the Solar gravitation term (8) Development of the time delay from (17) as τ = − (b·s") c provides the conventional group delay model (4). Now it is obvious that this model is based on the modification of the Lorentz transformation (12) in which the transformation of time is presented in a non-standard way because our reference clocks are physically located at the Earth surface rather than at the geocenter.