© The Authors 2025

1 Introduction

Geodetic Very Long Baseline Interferometry (VLBI) measures the difference in arrival times of radio waves – referred to as ‘group delay’ – at two radio telescopes receiving signals from distant sources, and a relative fringe frequency (also known as delay rate) which compensates for the geocentric motion of each radio telescope due to the Earth’s rotation. This measurement achieves remarkable precision, with accuracies of 20–40 picoseconds (ps) for the group delay and 1–100 frad for the fringe frequency (Schuh & Behrend 2012). VLBI observations are conducted in two frequency bands: 2.3 GHz (S-band) and 8.4 GHz (X-band), where the lower frequency serves to calibrate ionospheric fluctuations present in the X-band data. By combining data from a network of radio telescopes distributed globally, this technique can determine Earth orientation parameters (EOP) with an accuracy of approximately 0.1 milliarcseconds (mas) during a standard 24-hour geodetic VLBI experiment (Bizouard et al. 2019). The positions of the radio telescopes, in Cartesian coordinates, are provided by the International Terrestrial Reference Frame (ITRF), while the coordinates of reference radio sources are taken from the International Celestial Reference Frame (ICRF3) catalogue. Historically, the group delay measured by geodetic VLBI has been used to derive polar motion relative to the Celestial Intermediate Pole (CIP) within the ITRF, utilising three parameters: the two components of polar motion (x – and y-pole motion) and Universal Coordinated Time (UTC). The CIP is linked to the Geocentric Celestial Reference Frame (GCRF) pole through a precession-nutation model. This traditional concept relies on the position of the equinox, which is steadily moving along the ecliptic due to the motion of the Earth’s rotational axis with respect to the distant celestial objects, with a period of about 26 000 years. The alternative method proposes an introduction of the non-rotating origin (NRO) instead of the equinox point (Guinot 1979). To connect the instantaneous rotation pole and the celestial reference frame, new angles (X and Y) must be introduced; these are linked to the precession-nutation parameters (Capitaine 1990). The theory links the position of the so-called celestial intermediate pole (CIP) with the celestial pole offsets X and Y, and the corrections (dX and dY) should be estimated. Fukushima (2001) showed that the NRO adopted by the International Astronomical Union (IAU) leads to the appearance of a global secular rotation of the celestial ephemeris origin with a rate of −4.15″/year, and suggested additional theoretical options (Fukushima 2003). However, this NRO approach has not been widely used to date due to a lack of observations that could directly provide the X and Y values. Here, we show that VLBI delay rate observations produce three components of Earth’s angular rotation vector, which can be used to calculate the Ω magnitude. The conversion of these components to the celestial pole positions (X and Y angles) will be presented in future publications.

2 Basic VLBI delay rate equations

The fringe frequency measured with VLBI reflects the generalised Doppler effect, which describes the frequency shift between an ‘emitter’ and a ‘receiver’ moving relative to each other and to the direction of signal propagation. The resultant frequency shift is the difference between the two Doppler shifts observed at each radio telescope. Since the signal originates from the same radio source, the contribution from the ‘emitter’ cancels out, leaving only the contributions from the two ‘receivers’. It should be noted that the fringe frequency is calculated simultaneously with the group delay during the correlation process, rather than derived from the group delay (Sovers et al. 1998). The primary term of the fringe frequency can be expressed as shown by Cohen & Shaffer (1971): $$\frac{\mathrm{\Delta }{f}_{21}}{f}=-\frac{(w_2\cdot s)-(w_1\cdot s)}{c},$$(1)

where w₁ and w₂ are vectors of the geocentric velocities of the radio telescopes with geocentric positions assigned by vectors r₁ and r₂, s is the unit vector in the direction of the radio source, c is the speed of light, and (·) denotes the dot product. This analytical equation Eq. (1) can be considered as the derivative $\frac{\mathrm{\Delta }{f}_{21}}{f}=\frac{d\tau }{dt}$ of the group delay, $\tau =-\frac{(b_{21},s)}{c}=-\frac{(r_2\cdot s)-(r_1\cdot s)}{c}$, followed by the complete relativistic equation provided in Appendix A (again noting that the observed value $\frac{\mathrm{\Delta }{f}_{21}}{f}$ is not a derivative of the observed group delay). For rotational motion, the linear velocity vector is calculated as w = [Ω × r], where Ω is the vector of the Earth’s angular rotation (Ω = (ω_x, ω_y, ω_z), |Ω| ≈ 7.2921151467 · 10⁻⁵rad/s), r is the vector of the radio telescope geocentric position, and [×] denotes the vector cross product. Equation (1) can then be written as $$\frac{\mathrm{\Delta }{f}_{21}}{f}=-\frac{([\mathrm{\Omega }\times b_{21}]\cdot s)}{c}=-\frac{([s\times b_{21}]\cdot \mathrm{\Omega })}{c},$$(2)

where the baseline vector b₂₁ = r₂ − r₁ now appears explicitly. The magnitude of the effect $\frac{\mathrm{\Delta }{f}_{21}}{f}$ for the Earth-based baseline is approximately 0.1 ÷ 1 μrad. A full relativistic equation for the fringe frequency includes additional terms describing the effects of special and general relativity Eqs. (A.1) and (A.3) to the level of 1 frad =10⁻¹⁵ rad. The corresponding expected formal error of the vector Ω components is 10⁻¹²−10⁻¹³ rad/s, or 10⁻⁸−10⁻⁹ in relative units, $\frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{\Omega }}$.

2.1 Solution setup

The positions of the celestial intermediate pole in the X and Y directions are predicted by the IAU2000A/IAU2006 nutation and precession models (Mathews et al. 2002; Capitaine et al. 2003). The standard approach focusses on analysing group delay and applying corrections to the rotation angle values (dX, dY) to calculate minor adjustments to precession and nutation amplitudes, e.g. (Petit & Luzum 2010; Belda et al. 2017). The least-squares method is then employed to solve the system of equations after incorporating all relativistic corrections. Typically, this method requires a priori information about the estimated EOP. While it is often asserted that geodetic techniques cannot measure the X and Y angles directly (Cerveira et al. 2009; Schreiber et al. 2023), this paper demonstrates that the fringe frequency measured with geodetic VLBI, as described in Eq. (2), can effectively fill this gap. Eq. (2) encompasses the complete vector of angular rotation Ω, rather than merely small LOD corrections.

The Earth’s angular rotation vector Ω defines the positions of the IRP. Thus, the IRP direction could be extracted using the VLBI delay rate and Eq. (2). However, while the two poles, CIP and IRP, are close enough, their definitions are not equivalent, and the transition of Ω = (ω_x, ω_y, ω_z) to the IRP is not discussed in this paper.

The system of equations for the least-squares method accounts for frequency fringe variations resulting from the non-hydrostatic component of the lower troposphere. The model for tropospheric delay rates is detailed in Eq. (B.2). The positions of reference radio sources, geodetic site coordinates, and the EOP were fixed; moreover, a simple parametrical model to adjust the difference in hydrogen maser frequencies, including the constant frequency offset, was used. Consequently, the model comprises seven parameters per radio telescope (six for the tropospheric delay rate and one for the difference between the two hydrogen maser frequencies) alongside three components of the vector Ω.

It is important to note that the group delays (used in the standard method) are affected by hydrogen maser phase fluctuation offsets and the local troposphere wet component instability. This necessitates the development of complex algorithms to mitigate the effects of phase variations when using standard group delay for data analysis. In contrast, the fringe frequency in Eq. (2) is not affected by such phase fluctuations. Therefore, a simple least-squares method is effective, yielding daily estimates of the three components of the vector Ω for each 24-hour VLBI experiment.

We processed VLBI data used in the operational service from April 1993 to April 2024 (NEOS-A in 1993–2001 and IVS-R1, R4 in 2002–2024) with the OCCAM software platform (version 6.3) (Titov et al. 2004) to estimate the three components of the vector Ω using the aforementioned procedure. All standard reductions of geodetic VLBI data adhere to the IERS Conventions 2010 (Petit & Luzum 2010). Station positions and velocities were fixed based on ITRF2020 values (Altamimi et al. 2023), with a priori EOP sourced from IERS C 04 data and the IAU 2000A nutation model. Radio source positions were anchored by the ICRF3 catalogue coordinates (Charlot et al. 2020). A priori zenith delays were determined from local troposphere parameter values, subsequently mapped to observation elevation using the Vienna mapping function (VMF1) (Boehm et al. 2006).

There are two options for estimating the vector Ω: either add the major contribution (Eq. (1)) to the full relativistic model, or omit it. In the former case, only small corrections to the three components of the vector Ω are obtained; in the latter scenario, all three components are estimated in full. In this paper, the latter option is used.

3 Results

Three components of the vector Ω = (ω_x, ω_y, ω_z), obtained over the 30-year period of observations, are shown in Figs. 1–3, respectively. They represent approximately 1 / 1000th of the magnitude of the 26 000-year circle induced by the precession motion of the rotation pole, overlapped with the main 18.6-year nutation and the seasonal signal. The formal errors on all three graphs are approximately 1–10 prad/sec and are too small to be visible in the figures. These results are in good agreement with the time series of the vector Ω components published by the IERS(¹). Figs. 4–6 display the differences between the IERS time series and the newly calculated estimates shown in Figs. 1–3. The root mean square (rms) of the residuals in Figs. 4–6 are 1.2 prad/s, 1.2 prad/s, and 2.7 prad/s for the X, Y, and Z components, respectively. The results presented in Figs. 1–3 will be discussed in the following three subsections.

3.1 Effect of general relativity

Since the general relativity effect shifts the radio source positions, the analytical model in Eq. (A.3) should be applied for the fringe frequency reduction. Because the source positions are not estimated within the solution setup, the general relativity effect propagates into other parameters, particularly the Earth’s angular rotation velocity. This effect can optionally be switched off to compare the resultant effect of general relativity on the Ω and LOD time series. Figure 8 presents two sets of time series of the annual angular velocity values, with and without the reduction for the Sun’s gravitational field. Figure 9 shows the difference between the two graphs in Fig. 8. The offset varies from 0.25 to 0.40 prad/s for all years between 1993 and 2024, with a mean value approximately 0.33 prad/s. The corresponding offset in LOD is approximately 0.3 msec. Thus, neglecting the effects of general relativity will lead to an underestimation of the Earth’s angular velocity.

Fig. 1 Daily estimates of component ωX obtained by least-squares method using Eq. (2).

Fig. 2 Daily estimates of component ωY obtained by least-squares method using Eq. (2).

Fig. 3 Daily estimates of component ωZ obtained by least-squares method using Eq. (2).

Fig. 4 Difference in ωX between the IERS time series and the new estimates shown in Fig. 1.

Fig. 5 Difference in ωY between the IERS time series and the new estimates shown in Fig. 2.

Fig. 6 Difference in ωZ between the IERS time series and the new estimates shown in Fig. 3.

Fig. 7 Daily values of the full angular velocity Ω, derived using Eq. (3).

Fig. 8 Annual amplitude of the Earth’s angular rotation from solutions with and without the contribution of general relativity (black dots and red dots, respectively).

Fig. 9 Difference between the two time series shown in Fig. 8. The mean value is approximately 0.33 prad/s.

3.2 Angular velocity and length of day

The magnitude of the angular velocity Ω = |Ω| is calculated from Eq. (3) for each 24-hour geodetic VLBI experiment: $$\mathrm{\Omega }=\sqrt{{\omega }_x^2+{\omega }_y^2+{\omega }_z^2}.$$(3)

The full magnitude of Ω is shown in Fig. 7 and displays the subtle variations associated with the length of day. Again, the formal errors on the graph are approximately 1–10 prad/sec and not visible here.

To evaluate the quality of the results, the daily values of Ω can be converted to the official LOD values published by the IERS (²). The equation linking Ω from Eq. (3) and LOD is given, for example, by (Lambeck & Cazenave 1973; Bizouard et al. 2022; Baenas et al. 2021): $$\frac{\mathrm{\Delta }LOD}{LO{D}_0}=-\frac{\mathrm{\Omega }-{\mathrm{\Omega }}_0}{{\mathrm{\Omega }}_0},$$(4)

where Ω₀ is the so-called ‘nominal’ angular velocity, Ω₀ = 72 921 151.467 prad/s, and LOD₀ = 86 400 s.

It is common to link the vector Ω components (ω_x, ω_y, ω_z) to the relative perturbations (m₁, m₂, m₃) induced by external forces (e.g. the Earth’s atmosphere) in the form (m₁, m₂, 1 + m₃), so that the angular velocity module is approximated to the first order by Ω (1 + m₃), where m₃ is the axial change. The relative variations of the length of day in Eq. (4) are then $$\frac{\mathrm{\Delta }LOD}{LO{D}_0}\approx -m_3.$$(5)

Traditionally, all fluctuations in LOD have been interpreted solely as a consequence of axial changes, i.e. 1 + m₃. A direct comparison of the vector Ω and LOD, obtained independently, reveals systematic discrepancies that can likely be attributed to the influence of minor constituents typically overlooked in standard LOD approximations (Fig. 10). However, a detailed interpretation of these systematics is beyond the scope of this paper.

Fig. 10 Comparison of LOD values. Top: yearly time series of LOD values calculated from the angular rotation amplitude Eq. (4); black dots) and IERS values (red dots). Bottom: difference between the two LOD time series.

Fig. 11 Time series of vector Ω estimates with 2-hour resolution. Top panel: data from the IVS-R1 experiment R1138 on 15 January, 2024. Bottom panel: data from the VGOS experiment VO-4059 on 28 February, 2024. Daily estimates of the vector Ω are shown in the middle of each graph and are identified by the smaller formal errors.

3.3 Higher temporal resolution results

It is often stated in the literature that the geodetic VLBI technique can only achieve Ω time series with 24-hour resolution (e.g. Schreiber et al. 2023). However, this limitation primarily applies to the standard time series of the EOP produced by the IERS. Much higher temporal resolution for Ω can be achieved using VLBI fringe frequency data.

Fig. 11 shows the 2-hour time series of the vector Ω estimates obtained for two recent geodetic VLBI experiments. The first is the standard S/X experiment conducted on 15 January 2024 IVS-R1-1138, and the second is the broadband VGOS experiment on 28 February 2024 (IVS-V0-4059). Since the number of VGOS observations (16 975) is substantially higher than that of S/X (3 940), the average number of observations in each 2-hour interval is around four times larger for the VGOS experiment. The formal errors of the 2-hour estimates vary between 4 and 16 prad/s for the S/X experiment, and between 3 and 7 prad/s for the VGOS experiment, respectively. It should be noted that two hours is not the ultimate limit for the resolution of the Ω time series. The resolution could be shortened to one hour or even smaller if a dedicated experiment is conducted for the purpose of angular rotation monitoring. Ultimately, this involves a trade-off between time series resolution and Ω uncertainty for each individual time bin.

Fig. 12 shows the estimates from all IVS-R1 and IVS-R4 experiments during 2024, compared with a set of VGOS experiments. The VGOS results display a smaller formal error and fewer range-to-range variations compared to the S/X counter-parts. However, the VGOS estimates are biased with respect to the S/X estimates by about 5 prad/s. This suggests that the nominal fringe frequency f for VGOS observations used in Eq. (1) requires some correction. Once this frequency issue is addressed, VGOS data will provide a highly effective tool for monitoring the Earth’s rotation angular velocity with record temporal resolution.

Fig. 12 Daily corrections to the Ω nominal value from IVS-R1 and IVSR-4 experiments (black dots) and VGOS experiments (red dots) during 2024.

4 Conclusion

In conclusion, conventional astronomical methods do not allow the precise determination of the direction of Earth’s axis of rotation, as there is no physical marker in the sky indicating this direction. As a result, it is often assumed that it is fundamentally impossible to measure the direction of the axis using astronomical methods. However, geodetic VLBI, being a non-standard astronomical technique, is capable of measuring the motion of the ‘rotating platform’ (the Earth) using the Doppler effect, as presented in Eq. (1). This capability, previously neglected, enables the kinematic determination of the direction of the Earth’s axis of rotation using Eq. (2). This approach defines the location of an imaginary marker on the sky, which can be used for practical applications. The accuracy is comparable to the accuracy of ring laser measurements, and, moreover, the temporal resolution of the observational time series is not limited to the standard 24-hour interval. Instead, the resolution can be reduced to 1–2 hours, provided a sufficient number of observations are available. A possible application of this new method is the monitoring of the Earth’s rotational velocity in near real-time mode.

Acknowledgements

This paper is published with the permission of the CEO, Geoscience Australia. The authors acknowledge use of VLBI data obtained within the framework of the International VLBI Service for Geodesy and Astrometry (IVS). Length of Day (LOD) results published by the International Earth Rotation Service (IERS). I am grateful to Dr. Christian Bizouard from Paris Observatory and an anonymous reviewer for his comments and interesting discussion.

Appendix A Relativistic model for delay rate

Relativistic model for the observed geodetic VLBI fringe frequency Eq. (1) is calculated by differentiation of the group delay relativistic model from the IERS Conventions 2010 (Petit & Luzum (2010), Kopeikin (1990), Klioner (2003), Soffel et al. (2017)) so the final equation could be presented as follows $$\begin{array}{r}\frac{\mathrm{\Delta }{f}_{21}}{f}=-\frac{((w_2-w_1)\cdot s)}{c}\cdot (1-\frac{2GM}{c^{2}r}-\frac{(V\cdot s{)}^2}{c^2}+\frac{V^2}{c^2}-\frac{(w_2\cdot s)}{c}-\frac{(w_2\cdot V)}{c^2})+\frac{(b_{21}\cdot s)(a\cdot s)}{c^2}+\frac{(b_{21}\cdot s)(a_2\cdot s)}{c^2}-\\ -\frac{(b_{21}\cdot a)}{c^2}-\frac{((w_2-w_1)\cdot V)}{c^2}+\frac{((w_2-w_1)\cdot V)(V\cdot s)}{2c^3}\end{array}$$(A.1)

where w₁ and w₂ are the vectors of geocentric velocities of the radio telescopes, r is the geocentric distance to the Solar centre, s is the vector of direction to the radio source, b₂₁ is the baseline vector, c is the speed of light, V is the the vector of the Earth’s barycentric velocity, a is the vector of the Earth’s barycentric acceleration, a₂ is the vector of the geocentric acceleration of the second station, G is the gravitational constant, M is the Solar mass, and (·) denotes the dot product.

The effect of general relativity is calculated in the analytical form based on the model of the light deflection angle (Titov & Girdiuk (2015)) $${\tau }_{21}=\frac{2GM((b_{21}\cdot N)-(b_{21}\cdot s)(N\cdot s))}{c^3(\left|r_2\right|-(r_2\cdot s))}.$$(A.2)

Here N is the unit vector towards the Sun’s centre from the “second” telescope, i.e. $N=\frac{r_2}{\left|r_2\right|}$.

Several vectors are variable in Eq. (A.2), but only the baseline vector changes rapidly enough to produce a measurable frequency shift for any radio source at any angular distance from the Sun. The corresponding model for the fringe frequency is developed by differentiation of the baseline vector in Eq. (A.2) $$\frac{\mathrm{\Delta }{f}_{21}}{f}=\frac{2GM((w_2-w_1)(N-s(N\cdot s))}{c^3(\left|r_2\right|-(r_2\cdot s))}=\frac{2GM(w_2-w_1)[s\times [s\times N]]}{c^3(\left|r_2\right|-(r_2\cdot s))}=\frac{2GM[\mathrm{\Omega }\times b_{12}]\cdot [s\times [s\times N]]}{c^3(\left|r_2\right|-(r_2\cdot s))}$$(A.3)

Appendix B Troposphere modeling for delay rate

The wet troposphere modeling for the group delay is given by (Boehm et al. 2006) $${\tau }_w(z)=m_w(z){\tau }_w\left(z_0\right)+m_w(z)(G_N\cos A+G_E\sin A)$$(B.1)

This model includes three parameters – zenith delay τ_w(z₀) for the zenith direction (z₀ = 0) and two gradients (North-South and East-West), G_N and G_E, and the mapping function m_w(z) to consider the actual zenith distance z of the observed radio source. The corresponding model for the delay rate includes also derivatives for the zenith delay rate and both gradients $$\begin{array}{r}{\left(\frac{\mathrm{\Delta }f}{f}\right)}_w(z)={\dot{m}}_w(z){\tau }_w\left(z_0\right)+m_w(z){\left(\frac{\mathrm{\Delta }f}{f}\right)}_w\left(z_0\right)+{\dot{m}}_w(z)\tan (z)(G_N\cos A+G_E\sin A)+\\ +m_w(z)\tan (z)({\dot{G}}_N\cos A+{\dot{G}}_E\sin A-G_N\dot{A}\sin A+G_E\dot{A}\cos A)+m_w(z)\frac{\dot{z}}{\cos {}^{2}z}(G_N\cos A+G_E\sin A)\end{array}$$(B.2)

where six parameters are estimated in total in Eq. (B.2), including three new derivatives: ${\left(\frac{\mathrm{\Delta }f}{f}\right)}_w\left(z_0\right),{\dot{G}}_N,{\dot{G}}_E$.

References

All Figures

	Fig. 1 Daily estimates of component ωX obtained by least-squares method using Eq. (2).
In the text

	Fig. 2 Daily estimates of component ωY obtained by least-squares method using Eq. (2).
In the text

	Fig. 3 Daily estimates of component ωZ obtained by least-squares method using Eq. (2).
In the text

	Fig. 4 Difference in ωX between the IERS time series and the new estimates shown in Fig. 1.
In the text

	Fig. 5 Difference in ωY between the IERS time series and the new estimates shown in Fig. 2.
In the text

	Fig. 6 Difference in ωZ between the IERS time series and the new estimates shown in Fig. 3.
In the text

	Fig. 7 Daily values of the full angular velocity Ω, derived using Eq. (3).
In the text

	Fig. 8 Annual amplitude of the Earth’s angular rotation from solutions with and without the contribution of general relativity (black dots and red dots, respectively).
In the text

	Fig. 9 Difference between the two time series shown in Fig. 8. The mean value is approximately 0.33 prad/s.
In the text

	Fig. 10 Comparison of LOD values. Top: yearly time series of LOD values calculated from the angular rotation amplitude Eq. (4); black dots) and IERS values (red dots). Bottom: difference between the two LOD time series.
In the text

	Fig. 11 Time series of vector Ω estimates with 2-hour resolution. Top panel: data from the IVS-R1 experiment R1138 on 15 January, 2024. Bottom panel: data from the VGOS experiment VO-4059 on 28 February, 2024. Daily estimates of the vector Ω are shown in the middle of each graph and are identified by the smaller formal errors.
In the text

	Fig. 12 Daily corrections to the Ω nominal value from IVS-R1 and IVSR-4 experiments (black dots) and VGOS experiments (red dots) during 2024.
In the text