Revisiting the light time correction in gravimetric missions like GRACE and GRACE follow-on

Abstract

The gravity field maps of the satellite gravimetry missions Gravity Recovery and Climate Experiment (GRACE ) and GRACE Follow-On are derived by means of precise orbit determination. The key observation is the biased inter-satellite range, which is measured primarily by a K-Band Ranging system (KBR) in GRACE and GRACE Follow-On. The GRACE Follow-On satellites are additionally equipped with a Laser Ranging Interferometer (LRI), which provides measurements with lower noise compared to the KBR. The biased range of KBR and LRI needs to be converted for gravity field recovery into an instantaneous range, i.e. the biased Euclidean distance between the satellites’ center-of-mass at the same time. One contributor to the difference between measured and instantaneous range arises due to the nonzero travel time of electro-magnetic waves between the spacecraft. We revisit the calculation of the light time correction (LTC) from first principles considering general relativistic effects and state-of-the-art models of Earth’s potential field. The novel analytical expressions for the LTC of KBR and LRI can circumvent numerical limitations of the classical approach. The dependency of the LTC on geopotential models and on the parameterization is studied, and afterwards the results are compared against the LTC provided in the official datasets of GRACE and GRACE Follow-On. It is shown that the new approach has a significantly lower noise, well below the instrument noise of current instruments, especially relevant for the LRI, and even if used with kinematic orbit products. This allows calculating the LTC accurate enough even for the next generation of gravimetric missions.

1 Introduction

The twin GRACE satellites observed Earth’s gravity field and, more importantly, the monthly time variations of it from the launch in 2002 until their reentry in 2017. These variations reflect the mass transport on large scale in and on Earth. The measurement principle is based on low-low satellite-satellite tracking (LL-SST), i.e. measuring distance variations between the orbiters, which are separated on the same polar orbit by approx. 200 km Tapley et al. (2004). The inter-satellite range variations were measured by the K-Band Ranging system (KBR) with a noise level of approx. 1 $\upmu \hbox {m}/\sqrt{\text {Hz}}$ at a Fourier frequency of 0.1 Hz, and with elevated noise towards lower frequencies.

Due to the enormous success of GRACE, a successor mission called GRACE Follow-On (GFO) was launched on May 22, 2018. Its payload, an evolved version of the original GRACE, is comprised of, among others, GNSS receivers for precise orbit determination, accelerometers for the measurement of non-gravitational accelerations, star cameras and inertial measurement units for attitude determination and the aforementioned KBR system Kornfeld et al. (2019). In addition, GRACE Follow-On hosts the novel Laser Ranging Interformeter (LRI) which is a technology demonstrator, and it is the first inter-satellite laser interferometer in space. It has demonstrated an excellent performance and reliability of all subsystems and exhibits a noise level of approx. 1 $\hbox {nanometer/}\sqrt{\text {Hz}}$ at a Fourier frequency of 0.1 Hz, well below the requirements Abich et al. (2019). The novel LRI and conventional KBR are operated in parallel and, since both should measure the same Euclidean distance variations after some post-processing corrections that are described below, inter-comparisons and cross-calibrations can be performed in order to characterize the instruments and their behavior.

Both instruments rely on the transmission of electro-magnetic radiation, back and forth, between the satellites. The LRI operates at an optical frequency of $\approx$ 281 THz in a so-called active transponder configuration Sheard et al. (2012), while the KBR, often also called the microwave ranging instrument (MWI), uses two microwave frequencies, one in the K and one in the Ka band, in the so-called dual one-way ranging (DOWR) scheme Kim (2000); Thomas (1999). Both instruments rely on tracking the phase of a beatnote signal at low radio frequencies ($\le 18$ MHz). The tracked phase is—up to an unknown offset—proportional to the travel time of the radiation between the orbiters, thus, proportional to the inter-satellite distance variations from an initial epoch where phase tracking started. When the phase measurements are rescaled to a displacement, they are usually referred to as biased range observations in the official data products.

The gravity field recovery algorithms usually are based on the corrected (i.e.  instantaneous Euclidean) biased range or on its time derivative called range rate Naeimi and Flury (2017). The former one means the Euclidean biased distance between both satellites’ center-of-mass at the same epoch, which differs from the measured biased range due to effects from the finite speed of light and due to the fact that the measurements are not referred to the center-of-mass. The difference between biased and corrected instantaneous range is usually expressed as the sum of three terms: the light-time correction, the ionospheric correction, and the antenna phase center correction—often called tilt-to-length coupling in the context of laser interferometry.

The LRI was designed to have a minimal tilt-to-length coupling, which has been confirmed by in-flight measurements to be below 150 $\upmu \hbox {m/rad}$ Wegener et al. (2020). The coupling is significantly lower than for the KBR Wu et al. (2006), where the reference point for the range measurement is offset by approx. 1.4 m from the center-of-mass. The ionospheric effect is also insignificant in the case of the LRI due to the shorter wavelength of the optical radiation. The ionospheric correction for the KBR is briefly addressed in this paper, mainly to show that there is a cross-coupling of ionospheric effect and light time correction (LTC), but it is highly suppressed to a level below picometers in the employed two-way measurements. The main focus of this paper lies on the LTC, which is relevant for KBR and LRI and which was mentioned first for the GRACE satellites in Thomas (1999). Later, Kim Kim (2000) described a method to analytically calculate the light time correction based on absolute spacecraft velocities, i.e. only the special relativistic contribution. Turyshev et al. Turyshev et al. (2014) established an extensive description of general relativistic observables in GRACE-like missions, which includes an analytical model for the light time correction, among others. However, in our opinion, it is not straightforward to apply the formalism to actual flight data, because the relevant LTC expressions are derived under the assumption of nearly-matched Keplerian orbits for the satellites and approximations are used to derive closed-form expressions for the LTC. This enables the authors to understand and discuss the individual terms, but is also a restriction with regard to generality.

Thus, we derive the light time correction from first principles, and stay close to the data products and processing strategy in gravimetric mission, such that the results are easily applicable. The potential of Earth’s gravity field is expressed in terms of Stokes coefficients of a spherical harmonic (SH) expansion and the equations are formulated with quantities available from the official public data of the missions. In the following Sect. 2, the equations of motion are introduced in the general relativistic context, which are needed to describe the propagation of electro-magnetic waves. The propagation time of light between satellites is derived and split into the contributions from relativity (Sect. 3) and atmosphere (Sect. 4). However, actual calculations require a solution of an implicit equation (Sect. 5), which can be solved iteratively or by means of an analytical approximation. The analytical approach offers some advantages, since it allows us to replace some orbit product quantities that drive the numerical precision with more precise ranging observations. The analytical solution is combined in Sect. 6 into the dual-way light time corrections for KBR and in Sect. 7 into the round-trip LTC for LRI. Section 8 addresses the sensitivity of the ranging instruments and sketches a potential goal for the precision of the analytical equations and background models for the LTC. In the subsequent section 9, the analytical expressions for the one-way LTC are verified against numerical results and a parameter study is performed regarding background model accuracy and degree of approximations. We compare our results for the LTC against the results from official datasets for GRACE and GRACE Follow-On in Sect. 10, while Sect. 11 addresses further potential improvements in the light time correction calculation.

2 Equations of motion in general relativity

In order to derive a precise light time correction, the travel time of light between satellites is needed in a general relativistic context. For this, it is convenient to describe the light or microwaves in terms of mass-less particles, the photons, which move on geodesics according to the equations of motion in general relativity. We denote the coordinates of an object in the Geocentric Celestial Reference System (GCRS) as:

$$\begin{aligned} x^{\alpha }=\left( c_{0} \cdot t, x, y, z\right) =\left( c_{0} \cdot t, \mathbf {r}\right) =\left( x^{0}, x^{1}, x^{2}, x^{3}\right) ^{\top } \end{aligned}$$

(1)

where the common four-vector notation from relativity is used, and $c_{0}$ is the proper speed of light for vacuum in a local Lorentz frame with a numerical value of 299 792 458 m/s, $\mathbf {r}$ is the three dimensional spatial vector.

We employ the sign convention $\gamma _{\alpha \beta } = \text {diag} \{-1, +1, +1, +1\}$ for the Minkowski metric as used, for instance, by Kopeikin et al. Kopeikin et al. (2011). The Greek indices such as $\alpha$ and $\beta$ range from 0 to 3, while Latin letters like m and n denote spatial components and range from 1..3. $c_{n}$ is the coordinate speed of light.

The metric tensor $g_{\alpha \beta }$ of the Earth in the GCRS is approximated by a Post-Newtonian expansion as Soffel and Langhans (2012) and Turyshev et al. (2014):

$$\begin{aligned} \begin{aligned} g_{00}&=\gamma _{00} + \frac{2 W}{c_{0}^{2}}-\frac{2 W^{2}}{c_{0}^{4}}+{\mathcal {O}}\left( c_{0}^{-6}\right) \\ g_{0m}&=g_{m0}=-\frac{4 \mathbf {V}_m}{c_{0}^{3}}+{\mathcal {O}}\left( c_{0}^{-5}\right) \\ g_{m m}&=\gamma _{mm}+\frac{2 W}{c_{0}^{2}}+{\mathcal {O}}\left( c_{0}^{-4}\right) \end{aligned} \end{aligned}$$

(2)

with

$$\begin{aligned} W = W_{e} + \sum _i W_{\text {cb},i} \end{aligned}$$

(3)

where $W_e$ is the classical Newtonian potential due to the mass distribution of the Earth. Moreover, W contains a sum of potentials $W_{\text {cb},i}$ giving rise to the direct tidal acceleration towards other celestial bodies, in particular the Sun and the Moon. The vector potential $\mathbf {V}$ in Eq. (2) accounts for Earth’s spin moment with $\mathbf {V}_m$ denoting the mth component of $\mathbf {V}$.

We describe the potential $W_{e}$ as the sum of a central term $W_\text {PM} = GM_\mathrm{e}/\mathrm{r}$ and of higher moments of the gravity field $W_\text {HM}$, i.e. 

$$\begin{aligned} W_{e} = W_\text {PM} + W_\text {HM} = W_\text {PM} + W_\text {G} + W_\text {tidal} +W_\text {non-tidal}, \end{aligned}$$

(4)

whereby $W_\text {HM}$ is formed by the higher moment of static mass distribution potential $W_G$, by the potential $W_\text {tidal}$ describing the distortion of the mass distribution due celestial bodies such as Moon and Sun, and by the non-tidal potential $W_\text {non-tidal}$ describing small variations in the atmosphere, oceans, hydrology, ice and solid earth (AOHIS). These non-tidal variations contain highly interesting information for Earth sciences and the measurement of them is the main objective of GRACE-like missions.

The potentials describing higher moments of the gravity field are usually expressed in terms of a SH expansion Heiskanen and Moritz (1967):

$$\begin{aligned} W_\text {HM}(r, \varTheta ,\lambda )= & {} \frac{G M_{e}}{R_{e}} \sum _{l=1}^{\infty }\left( \frac{R_{e}}{r}\right) ^{(l+1)} \sum _{k=0}^{l}\left( {\overline{C}}_{l k} \cos (k \lambda )\right. \nonumber \\&\left. +{\overline{S}}_{l k} \sin (k \lambda ) \right) {\overline{P}}_{l k}(\cos \varTheta ) \end{aligned}$$

(5)

where G is the gravitational constant, $M_{e}$ is the mass of the Earth, $R_{e}$ is Earth’s average radius, $(r, \varTheta ,\lambda )$ are the spherical position coordinates, ${\overline{P}}_{l k}$ are the normalized Legendre functions of the second kind, l and k are the degree and order of the series expansion, and ${\overline{C}}_{l k}$ and ${\overline{S}}_{l k}$ are the normalized dimensionless Stokes coefficients. The Stokes coefficients of the static, tidal and non-tidal models given in Table 1 can be summed up in order to yield the total field $W_\text {HM}$.

Table 1 List of background models used in calculations

The direct acceleration towards a celestial body, which is often called direct tidal acceleration, has in the Earth-centered frame the potential $W_{\text {cb},i}$ Montenbruck et al. (2002):

$$\begin{aligned} W_{\text {cb},i}=\frac{G M_{\text {cb},i}}{R_{\text {cb},i}} \sum _{l=2}^{\infty }\left( \frac{r}{R_{\text {cb},i}}\right) ^{l} {\overline{P}}_{l}(\cos \varsigma _i) \end{aligned}$$

(6)

where G is the gravitational constant, $M_{\text {cb},i}$ is the mass the of i-th celestial body, $R_{\text {cb},i}$ is the distance between Earth and celestial body, r is the distance between Earth center and the satellite, ${\overline{P}}_{l}$ are the normalized Legendre functions of the first kind, $\varsigma _i$ is the angle between $\mathbf {R}_{\text {cb},i}$ and the satellite position vector $\mathbf {r}_{s}$, and l is the degree of the series expansion. In this paper we consider only the Sun and the Moon, since they are dominating the direct tidal acceleration.

The vector potential $\mathbf {V}$ in Eq. (2) is usually approximated as Turyshev et al. (2014):

$$\begin{aligned} \mathbf {V}(t, \mathbf {r}) \approx \frac{G M_{e}}{2 \cdot r^{3}} \cdot \mathbf {S} \times \mathbf {r}+{\mathcal {O}}\left( x^{-4}, c^{-2}\right) \end{aligned}$$

(7)

where $\mathbf {S}$ is Earth’s spin moment, or its angular momentum per unit of mass. It can be approximated by the angular momentum of a homogeneous sphere:

$$\begin{aligned} \mathbf {S} \approx \frac{2}{5} \cdot R_{e}^{2} \cdot {\varvec{\omega }}_{e} \end{aligned}$$

(8)

where ${\varvec{\omega }}_{e}$ is Earth’s angular velocity vector.

The equations of motions of a point particle, e.g. satellites or light read in the context of General Relativity (GR) as Kopeikin et al. (2011):

$$\begin{aligned} \frac{ \text {d}^2 x^k }{ \text {d} t^2 } = -\varGamma ^{k}_{~\alpha \beta } \cdot \frac{\text {d} x^\alpha }{\text {d} t} \cdot \frac{\text {d} x^\beta }{\text {d} t} + \frac{1}{c_0} \varGamma ^{0}_{~\alpha \beta } \cdot \frac{\text {d} x^\alpha }{\text {d} t} \cdot \frac{\text {d} x^\beta }{\text {d} t} \cdot \frac{\text {d} x^k}{\text {d} t} \qquad \text {with} \quad k = 1..3, \end{aligned}$$

(9)

where t is the coordinate time, and $\varGamma ^{k}_{~\alpha \beta }$ are the Christoffel symbols, which depend on derivatives of the metric tensor $g_{\alpha \beta }$. It is straightforward to numerically integrate these differential equations in order to obtain a trajectory for a given set of initial conditions. For a photon, the trajectory appears bent with approximately twice the classical Newtonian acceleration towards Earth’s center, consistent with one of the very early results of GR Einstein (1911); Soldner (1921). The selection of the initial velocity of a photon requires the coordinate speed of light, which depends on the metric tensor and on the propagation direction. It can be derived from the following ansatz for the four velocity:

$$\begin{aligned} \frac{\mathrm {d} x^{\alpha }}{\mathrm {d} t}=\left( c_{0}, \mathbf {d_{0}} .c_{n}\right) ^{\mathrm {T}} \end{aligned}$$

(10)

where $c_{n}$ is the coordinate speed of light in a vacuum in the GCRS, $\mathbf {d_{0}}$ is the normalized propagation direction of the photon and t is the coordinate time of the GCRS.

The interval $\text {d}s^2$ of a world line or trajectory of a massless particle vanishes Kopeikin et al. (2011):

$$\begin{aligned} \mathrm {d} s^{2}=g_{\alpha \beta }(t, \mathbf {r}) \cdot \partial x^{\alpha } \cdot \partial x^{\beta }=0. \end{aligned}$$

(11)

After dividing by $\text {d}t^{2}$ and plugging Eq. (2) into Eq. (11), one obtains a quadratic equation for $c_n$

$$\begin{aligned} \begin{aligned} 0&=g_{\alpha \beta }(t, \mathbf {r}) \cdot \mathrm {d} x^{\alpha } / \mathrm {d} t \cdot \mathrm {d} x^{\beta } / \mathrm {d} t \\&=c_{0}^{2} \cdot g_{00}+\mathbf {G} .\mathbf {d_{0}} \cdot c_{n} \cdot c_{0}+c_{n}^{2} \cdot g_{mm}, \end{aligned} \end{aligned}$$

(12)

where $\mathbf {G}=2(g_{01},g_{02},g_{03})^{T}=-8\mathbf {V}/c_{0}^{3}$ and $g_{mm} = g_{11}= g_{22}= g_{33}$. The post-Newtonian effect is very small, such that $g_{00}$ and $g_{mm}$ are close to unity. The quadratic equation can be solved and the solution with positive propagation velocity is taken for the coordinate speed of light:

$$\begin{aligned} c_{n}= & {} c_{0} \cdot \sqrt{-\frac{g_{00}}{g_{mm}}+\frac{\left( \mathbf {G} .\mathbf {d}_{0}\right) ^{2}}{4 \cdot \left( g_{mm}\right) ^{2}}}-c_{0} \cdot \frac{\mathbf {G} .\mathbf {d}_{0}}{2 \cdot \left( g_{mm}\right) }\nonumber \\= & {} \sqrt{\frac{c_{0}^{6}-2 \cdot c_{0}^{2} \cdot W^{2}+4 \cdot W^{3}+16 \cdot \left( \mathbf {V} .\mathbf {d}_{0}\right) ^{2}}{\left( c_{0}^{2}+2 \cdot W\right) ^{2}}}\nonumber \\&+\frac{4 \cdot \mathbf {V} .\mathbf {d}_{0}}{c_{0}^{2}+2 \cdot W}. \end{aligned}$$

(13)

The infinitesimal propagation time $\text {d}t$ of a photon is related to the coordinate pathlength $\text {d}s$ through

$$\begin{aligned} \text {d}t = \frac{n}{c_{n} } \cdot \text {d}s = \frac{1+2 \cdot W / c_{0}^{2}-4 \cdot \mathbf {V} .\mathbf {d}_{0} / c_{0}^{3}}{c_{0}} \cdot n \cdot \text {d}s+{\mathcal {O}}\left( c_{0}^{-5}\right) , \end{aligned}$$

(14)

where n denotes the refractive index at the location of the photon.

For a one-way ranging measurement, the propagation time $\varDelta t$ of a photon traveling along path ${\mathcal {P}}$ can be written as

$$\begin{aligned} \varDelta t=\int _{{\mathcal {P}}} \frac{n}{c_{n}\left( t, \mathbf {r}_{\mathrm {ph}}\right) } \mathrm {d} s \approx \underbrace{\int _{{\mathcal {P}}} \frac{1}{c_{n}\left( t, \mathbf {r}_{\mathrm {ph}}\right) } \mathrm {d} s}_{\varDelta t_\text {rel}} + \underbrace{\frac{1}{c_0} \int _{{\mathcal {P}}} (n-1)~\mathrm {d} s}_{\varDelta t_\text {media}}, \end{aligned}$$

(15)

where $\mathbf {r}_{\mathrm {ph}}$ is the position of the photon on the path ${\mathcal {P}}$ and t is the coordinate time. Since $c_n$ is close to $c_0$ and since the effect of the refractive index due to the ionospheric and neutral atmosphere is small, such that $(n-1)$ is close to zero, it is possible to approximate the integral as the sum of the relativistic effect ($\varDelta t_\text {rel}$) and a contribution from the refractive index of the media ($\varDelta t_\text {media}$). Both effects are analyzed in more detail in the next two subsections.

3 Light time correction $\varDelta t_\text {rel}$ due to relativity

The light path ${\mathcal {P}}$ between satellites in a gravimetric mission can be assumed as a straight line in the three-dimensional coordinate system, which can be parameterized by a parameter $\lambda \in [0,1]$:

$$\begin{aligned} \mathbf {r}_{\mathrm {ph}}(\lambda )=\mathbf {r}_{e}+\left( \mathbf {r}_{r}-\mathbf {r}_{e}\right) \cdot \lambda \end{aligned}$$

(16)

where $\mathbf {r}_{r}$ is the three-dimensional position of the photon reception and $\mathbf {r}_{e}$ is the three-dimensional position of the photon emission.

This neglects the relativistic light bending, which arises from an apparent acceleration $a_c$ of the photons towards the geocenter with twice the Newtonian acceleration Kopeikin et al. (2011), i.e. $a_c = 2GM_{e}/r^2$. The displacement of a photon in radial direction w.r.t. a straight line is of the order of $a_c \cdot (\varDelta t)^2 / 2 \approx 4~\upmu \text {m}$, where a propagation time of $\varDelta t = 200~\text {km} / c_0 \approx 0.66~\text {msec}$ and a satellite position of $r = 6731$ km was assumed. Temporal variations of the displacement due to higher moments of the gravity field are much smaller. In the domain of phasefronts, the light-bending yields a negligible static phasefront tilt of the order of $4~\upmu \text {m}/200~\text {km} \approx 2 \cdot 10^{-11} \text {rad}$.

Thus, one can anticipate that the light-time correction derived from the bent light path will differ only insignificantly from a correction derived on the straight line. The approximation is further justified in Sect. 9, where our simplified analytical results are compared to results obtained via numerically integrating Eq. (9) and thus, accounting for the full GR effects.

Evaluating the propagation time $\varDelta t_\text {rel}$ in Eq. (15) with the photon path ${\mathcal {P}}$ yields:

$$\begin{aligned} \varDelta t_\text {rel}&=\int _{\lambda =0}^{1} c_{n}^{-1}\left( t, \mathbf {r}_{\mathrm {ph}}\right) \cdot \left| \frac{\mathrm {d} \mathbf {r}_{\mathrm {ph}}}{\mathrm {d} \lambda }\right| \mathrm {d} \lambda =\left| \mathbf {r}_{e}-\mathbf {r}_{r}\right| \cdot \int _{\lambda =0}^{1} c_{n}^{-1}\left( t, \mathbf {r}_{\mathrm {ph}}\right) \mathrm {d} \lambda \end{aligned}$$

(17)

$$\begin{aligned}&\approx \underbrace{\frac{\left| \mathbf {r}_{r}-\mathbf {r}_{e}\right| }{c_{0}}}_{\varDelta t_\text {SR}}+ \underbrace{2 \cdot \varDelta t_\text {SR} \cdot \int _{\lambda =0}^{1} \frac{G M_{e}}{c_{0}^{2} \cdot \left| \mathbf {r}_{\mathrm {ph}}(\lambda )\right| } d \lambda }_{{\mathcal {T}}_{{\mathrm {PM}}}} \nonumber \\&\quad +\underbrace{2 \cdot \varDelta t_{{\mathrm {SR}}} \cdot \int _{\lambda =0}^{1} \frac{W_{{\mathrm {HM}}}\left( t(\lambda ), \mathbf {r}_{\mathrm {ph}}(\lambda )\right) }{c_{0}^{2}} \mathrm {d} \lambda }_{{\mathcal {T}}_{{\mathrm {HM}}}}\nonumber \\&\quad +\underbrace{\varDelta t_{{\mathrm {SR}}} \cdot \int _{\lambda =0}^{1} \frac{-4 \cdot \mathbf {V}\left( \mathbf {r}_{\mathrm {ph}}(\lambda )\right) .\mathbf {d}_{0}}{c_{0}^{3}} \mathrm {d} \lambda }_{{\mathcal {T}}_{ {\mathrm {SM}}}}, \end{aligned}$$

(18)

where terms with the order of $c_{0}^{-4}$ and smaller were omitted and where the normalized propagation direction of the photon $\mathbf {d_{0}}$ was abbreviated by

$$\begin{aligned} \mathbf {d}_{0}=\frac{\mathbf {r}_{r}-\mathbf {r}_{e}}{\left| \mathbf {r}_{r}-\mathbf {r}_{e}\right| }. \end{aligned}$$

(19)

In upper Eq. (18), the first term $\varDelta t_{\mathrm {SR}}$ is the propagation time from special relativity in flat space-time, the second term ${\mathcal {T}}_{\mathrm {PM}}$ is the time delay due to Earth’s central field, the third term ${\mathcal {T}}_{\mathrm {HM}}$ is the time delay from higher moments of the gravitational potential due to Earth’s mass distribution and due to other celestial bodies, and the fourth term ${\mathcal {T}}_{\mathrm {SM}}$ is the time delay due to Earth’s spin moment.

The term ${\mathcal {T}}_{{\mathrm {PM}}}$ is commonly called Shapiro time delay and it has a closed analytical form Turyshev et al. (2014)

$$\begin{aligned} {\mathcal {T}}_{{\mathrm {PM}}}= & {} \frac{2 \cdot G M_{e}}{c_{0}^{3}} \cdot \ln \left( \frac{\left| \mathbf {r}_{r}\right| +\mathbf {d_{0}} .\mathbf {r}_{r}}{\left| \mathbf {r}_{e}\right| +\mathbf {d}_{0} .\mathbf {r}_{e}}\right) \nonumber \\= & {} \frac{2 \cdot G M_{e}}{c_{0}^{3}} \cdot \ln \left( \frac{\left| \mathbf {r}_{r}\right| +\left| \mathbf {r}_{e}\right| +\left| \mathbf {r}_{r}-\mathbf {r}_{e}\right| }{\left| \mathbf {r}_{r}\right| +\left| \mathbf {r}_{e}\right| -\left| \mathbf {r}_{r}-\mathbf {r}_{e}\right| }\right) . \end{aligned}$$

(20)

The ${\mathcal {T}}_{{\mathrm {HM}}}$ integral can be readily approximated using the N-point trapezoidal rule,

$$\begin{aligned} {\mathcal {T}}_{{\mathrm {HM}}}^{(N-1)}&\approx \frac{2}{c_{0}^{2}} \cdot \sum _{n=0}^{N-1} \frac{W_{{\mathrm {HM}}}\left( {\tilde{t}}_{n}, \mathbf {r}_{\mathrm {ph}}\left( \lambda _{n}\right) \right) +W_{{\mathrm {HM}}}\left( {\tilde{t}}_{n+1}, \mathbf {r}_{\mathrm {ph}}\left( \lambda _{n+1}\right) \right) }{2}\nonumber \\&\quad \cdot \left( {\tilde{t}}_{n+1}-{\tilde{t}}_{n}\right) \end{aligned}$$

(21)

$$\begin{aligned}&= \frac{2 \cdot \varDelta t_{{\mathrm {SR}}}}{c_{0}^{2} \cdot N} \cdot \left( \sum _{n=1}^{N} W_{{\mathrm {HM}}}\left( {\tilde{t}}_{n}, \mathbf {r}_{\mathrm {ph}}\left( \lambda _{n}\right) \right) \right. \nonumber \\&\quad \left. + \frac{W_{{\mathrm {HM}}}\left( {\tilde{t}}_{N}, \mathbf {r}_{\mathrm {ph}}\left( \lambda _{N}\right) \right) + W_{{\mathrm {HM}}}\left( {\tilde{t}}_{0}, \mathbf {r}_{\mathrm {ph}}\left( \lambda _{0}\right) \right) }{2} \right) \end{aligned}$$

(22)

$$\begin{aligned} \text {with time } \quad {\tilde{t}}_{n}&= t(\lambda _0) + \varDelta t_{{\mathrm {SR}}} \cdot \lambda _{n} = t(\lambda _0) + \varDelta t_{{\mathrm {SR}}} \cdot \frac{n}{N}, \qquad 0 \le n \le N, \end{aligned}$$

(23)

with $(N-1)$ being the number of segments in the uniform grid sampling of the light path ${\mathcal {P}}$. Finally, the gravito-magnetic effect, the ${\mathcal {T}}_{{\mathrm {SM}}}$ term, can be approximated with a two-point trapezoidal rule as

$$\begin{aligned} {\mathcal {T}}_{{\mathrm {SM}}} \approx -\frac{2 GM_{e} R_{e}^{2}}{5 c_{0}^{3}} \cdot \left( {\varvec{\omega }}_{e} \times \mathbf {r}_{e}\right) .\mathbf {d}_{0} \cdot \left( \frac{1}{\left| \mathbf {r}_{e}\right| ^{3}}+\frac{1}{\left| \mathbf {r}_{r}\right| ^{3}}\right) \cdot \varDelta t_{{\mathrm {SR}}}. \end{aligned}$$

(24)

Anticipating the result, it is beneficial to separate the special relativistic contribution into a delay $\varDelta t_\text {inst}$ from the instantaneous inter-satellite range at the reception time $t_\mathrm{r}$ and into a special relativistic correction ${\mathcal {T}}_\text {SR}$, i.e. 

$$\begin{aligned} \varDelta t_{{\mathrm {SR}}}&= \frac{\left| \mathbf {r}_{r}-\mathbf {r}_{e}\right| }{c_{0}} = \frac{\left| \mathbf {r}_{B}(t_\mathrm{r})-\mathbf {r}_{A}(t_e)\right| }{c_{0}} = \frac{\left| \mathbf {r}_{B}(t_\mathrm{r})-\mathbf {r}_{A}(t_\mathrm{r})\right| }{c_{0}}\nonumber \\&\quad + {\mathcal {T}}_\text {SR} = \varDelta t_\text {inst} + {\mathcal {T}}_\text {SR}, \end{aligned}$$

(25)

where it was assumed without loss of generality that the light is received by satellite B after being emitted by satellite A at time $t_e = t_\mathrm{r} - \varDelta t$. In summary, the light propagation time $\varDelta t_\text {rel}$ can be written as

$$\begin{aligned} {\varDelta t_{\mathrm {rel}}} = \varDelta t_{\mathrm {inst}} + {\mathcal {T}} \end{aligned}$$

(26)

with the light-time correction ${\mathcal {T}}$ containing special and general relativistic contributions

$$\begin{aligned} {\mathcal {T}} = {\mathcal {T}}_{{\mathrm {SR}}} + {\mathcal {T}}_{\mathrm {GR}} = {\mathcal {T}}_{{\mathrm {SR}}} + {\mathcal {T}}_{{\mathrm {PM}}} + {\mathcal {T}}_{{\mathrm {HM}}} + {\mathcal {T}}_{{\mathrm {SM}}}. \end{aligned}$$

(27)

In order to compute all these terms, the emission position and emission time of the photon is needed, which depend on the light-time corrections. This yields an implicit light-time equation, which is solved in Sect. 5, after discussing the remaining correction for the atmosphere.

4 Light time correction $\varDelta t_\text {media}$ due to atmosphere

At orbit heights below approx. 500 km, such as the low Earth orbits of the GRACE and GRACE Follow-On satellites, the residual atmosphere may alter the speed of light due to refraction. A deviation of the refractive index n from unity arises due to the neutral atmosphere and due to free electrons in the ionosphere. The former effect is negligible for interferometric range measurements, i.e. for the time-delay $\varDelta t_\text {media}$, since the fluctuations are estimated to be below $2~\text {nm}/\sqrt{\text {Hz}}/c_0$ for mHz frequencies, and with sinusoidal variations below 1 nm/$c_0$ amplitude at once and twice the orbital frequency Müller (2017).

However, the propagation of electromagnetic waves needs to be modeled according to propagation laws in plasma due to the charged particles in the ionosphere between 75..1000 km height. The main correction to the propagation time is the first-order ionospheric delay, which is commonly expressed as Petit and Luzum (2010); Hernández-Pajares et al. (2007)

$$\begin{aligned} \varDelta t_\text {media}&= \frac{1}{c_0} \int _{{\mathcal {P}}} (n-1)~\mathrm {d} s \approx - \frac{40.3\,\text {Hz}^2 /\text {m}}{c_0 \cdot f_\text {em}^2} \cdot \frac{\text {TEC}}{1\,e^-/\text {m}^2} \nonumber \\&= - \frac{40.3\,\text {Hz}^2}{f_\text {em}^2} \cdot \frac{\langle \eta \rangle }{1\,e^-/\text {m}^3}\cdot \varDelta t_\text {SR}, \end{aligned}$$

(28)

where $f_\text {em}$ is the frequency of the electromagnetic wave and $\text {TEC}$ is the total electron content on the photon path with units of electrons per square-meter. The ionospheric delay is actually an advancement, since the correction is always negative, which is known from GNSS, where the code delay is positive, while the phase delay is negative. Due to the frequency dependence, it is possible to estimate variations of the TEC with interferometric range measurements at two different frequencies, but the absolute value of the TEC, and hence, the absolute value of $\varDelta t_\text {media}$ is not measurable, because the ranging instruments can determine only a biased range.

However, in order to simulate the effect, the TEC can be expressed as the product of the mean electron density $\langle \eta \rangle$ between the satellites and the geometrical inter-satellite distance $\varDelta t_\text {SR} \cdot c_0$. For satellites at a height of 400 km, the electron density can reach values of up to $\langle \eta \rangle =10^{12}~e^{-}/\text {m}^3$ Kelley (2009), which translate in worst-case to an absolute delay of $-13\,\text {mm}/c_0$ for a microwave frequency of $f = 24.5\,\text {GHz}$ and $\varDelta t_\text {SR} \approx 200\, \text {km} /c_0$. The effect of such a non-measurable absolute delay onto the instantaneous biased KBR range is assessed through the LTC in Sect. 6. On this occasion, we point out that ionospheric effects are negligible for laser ranging with an optical frequency of $281\,\text {THz}$, since the contributions in propagation time or biased range are reduced by the factor

$$\begin{aligned} \left( \frac{24.5\,\text {GHz}}{281\,\text {THz}} \right) ^2 \approx 7.6 \cdot 10^{-9} \end{aligned}$$

(29)

compared to the microwave K-band.

5 Solving the light-time equation

The propagation time $\varDelta t$ of electromagnetic waves or photons between the two satellites has been described so far as a function of the photon path, or more precisely, as a function of the emission time $t_e$, emission position $\mathbf {r}_e$, reception time $t_\mathrm{r}$ and reception position $\mathbf {r}_\mathrm{r}$.

We may assume that the satellite trajectories are known, in particular, the satellite position $\mathbf {r}_{A/B}$, velocity $\dot{\mathbf {r}}_{A/B}$ and acceleration $\ddot{\mathbf {r}}_{A/B}$ at the time of reception $t_\mathrm{r}$. The acceleration can be derived with a kinematic approach as time-derivative or by dynamic means using force models. Without loss of generality, we may assume that satellite B is the receiver such that the reception position becomes $\mathbf {r}_\mathrm{r} = \mathbf {r}_B(t_\mathrm{r})$ and that satellite A is the emitter.

Using Taylor expansion, the satellite’s trajectory can be approximated in the vicinity of $t_\mathrm{r}$ as

$$\begin{aligned} \mathbf {r}_A(t_\mathrm{r} - \epsilon ) \approx \mathbf {r}_A(t_\mathrm{r}) - \dot{\mathbf {r}}_A(t_\mathrm{r}) \cdot \epsilon + \ddot{\mathbf {r}}_A(t_\mathrm{r}) \cdot \epsilon ^2/2, \end{aligned}$$

(30)

which allows us to write the position at the event of photon emission as $\mathbf {r}_e = \mathbf {r}_A(t_\mathrm{r} - \varDelta t)$. In order to solve for $\varDelta t$ one has to solve the implicit equation

$$\begin{aligned} \varDelta t(t_\mathrm{r})&= \frac{| \mathbf {r}_B(t_\mathrm{r}) - \mathbf {r}_A(t_\mathrm{r} - \varDelta t) |}{c_0} + {\mathcal {T}}_\text {GR}(\mathbf {r}_e = \mathbf {r}_A(t_\mathrm{r} - \varDelta t)) \nonumber \\&\quad + \varDelta t_\text {media}(\mathbf {r}_e = \mathbf {r}_A(t_\mathrm{r} - \varDelta t)) \end{aligned}$$

(31)

A solution can be obtained by iterative means using

$$\begin{aligned} \varDelta t^{(n+1)}(t_\mathrm{r})&= \frac{| \mathbf {r}_B(t_\mathrm{r}) - \mathbf {r}_A(t_\mathrm{r} - \varDelta t^{(n)}) |}{c_0} + {\mathcal {T}}_\text {GR}(\mathbf {r}_e = \mathbf {r}_A(t_\mathrm{r} - \varDelta t^{(n)})) \nonumber \\&\quad + \varDelta t_\text {media}(\mathbf {r}_e = \mathbf {r}_A(t_\mathrm{r} - \varDelta t^{(n)})) \end{aligned}$$

(32)

with start value $\varDelta t^{(0)} = \varDelta t_\text {inst} = |\mathbf {r}_A(t_\mathrm{r}) -\mathbf {r}_B(t_\mathrm{r})|/c_0$. The three summands on the right hand side have an amplitude of approximately 200 $\hbox {km}/c_0$, 300 $\upmu \hbox {m}/c_0$ and in case of the K-band -13 $\hbox {mm}/c_0$, respectively.

The vectors in the first term have typically a magnitude of $7 \cdot 10^6$ meters, which limits direct numerical solutions of $\mathbf {r}_A - \mathbf {r}_B$ in Eq. (32) to a precision of the order of $\hbox {nanometer/}c_0$ due to the $\approx$ 15 digits precision of 64-bit (double) floating-point arithmetic. One way to overcome this limitation is to derive an analytical closed-form solution for $\varDelta t$. This can be achieved by substituting Eq. (30) into Eq. (32), taking into account the relation in Eq. (25), and evaluating the first few iterations using an algebraic manipulation software. If terms with negligible magnitude are omitted in the lengthy expression¹, the solution for $\varDelta t^{(2)}$ (and higher iteration numbers) reads

$$\begin{aligned}&\varDelta t(t_\mathrm{r}) = \varDelta t_\text {inst}(t_\mathrm{r}) + {\mathcal {T}}_\text {SR}(t_\mathrm{r}) + {\mathcal {T}}_\text {GR}(t_\mathrm{r}) + \varDelta t_\text {media}(t_\mathrm{r}) \end{aligned}$$

(33)

$$\begin{aligned}&{\mathcal {T}}_\text {SR} = \varDelta t_\text {inst} \frac{\mathbf {d}_0 .\dot{\mathbf {r}}_A }{c_0} + \varDelta t_\text {inst}^2 \frac{\mathbf {d}_0 .\ddot{\mathbf {r}}_A }{2 \cdot c_0} \nonumber \\&\quad + \frac{ \varDelta t_\text {inst}^2 \cdot (- \mathbf {d}_0 .\ddot{\mathbf {r}}_A \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_A - \dot{\mathbf {r}}_A .\ddot{\mathbf { r}}_A/2) + \varDelta t_\text {inst}/2 \cdot ( (\mathbf {d}_0 .\dot{\mathbf {r}}_A)^2 + |\dot{\mathbf {r}}_A|^2 )}{c_0^2} \nonumber \\&\quad + \frac{\varDelta t_\text {inst} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_A \cdot |\dot{\mathbf {r}}_A|^2 }{c_0^3} + ({\mathcal {T}}_\text {GR}+\varDelta t_\text {media}) \frac{\mathbf {d}_0 .\dot{\mathbf {r}}_A}{c_0} \nonumber \\&\quad + {\mathcal {O}} \left( 10^{-12}\,\text {m}/c_0 \right) \end{aligned}$$

(34)

$$\begin{aligned}&{\mathcal {T}}_\text {GR} = {\mathcal {T}}_\text {GR}(\mathbf {r}_e = \mathbf {r}_A(t-\varDelta t_\text {inst} - \varDelta t_\text {inst} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_A / c_0 ))\nonumber \\&~~~~\quad \approx {\mathcal {T}}_\text {GR}(\mathbf {r}_e = \mathbf {r}_A(t-\varDelta t)) \end{aligned}$$

(35)

where all quantities, also the one used for calculating $\mathbf {d}_0$ with Eq. (19), are evaluated at the photon reception time $t_\mathrm{r}$. Thus, the equation can be directly applied with orbit data from GRACE or GRACE Follow-On.

The overall light-time correction ${\mathcal {T}} = {\mathcal {T}}_\text {SR} + {\mathcal {T}}_\text {GR}$ is dominated by the first term in Eq. (34), which has an amplitude of the order of $-5$ $\hbox {m/}c_0$ for $\varDelta t_\text {inst} \approx 200\,\text {km}/c_0$ and $\mathbf {d}_0 .\dot{\mathbf {r}}_A \approx -7.6\,\text {km/s}$. The derivation of ${\mathcal {T}}$ assumed so far a single path of a photon from one satellite to the other, i.e. an one-way ranging approach. However, the ranging systems in GRACE and GRACE-Follow-On exchange light in both directions and the light-time correction becomes a linear combination of two (LRI) or four (MWI) one-way corrections (${\mathcal {T}}$). As will turn out subsequently, these linear combinations have a significantly lower magnitude due to a high common-mode rejection.

6 Light time correction in dual one-way ranging (DOWR)

Fig. 1

Minkowski diagram of the light path (red arrows) in a dual-one way ranging (DOWR) scheme at a particular frequency (left plot) and in the two-way ranging (TWR) scheme (right plot). For the DOWR, the emission (e) and reception (r) events are located at the antenna phase centers (grey trajectories) of the two satellites (A and B). In the TWR case, these events occur at the center-of-mass (solid black lines) of the master (M) and transponder (T) satellite. The reflection event on the transponder side is denoted as Tp

The dual-one way ranging concept is used by the microwave ranging systems in GRACE and GRACE Follow-On Wen et al. (2019), where the ionospheric effect needs to be removed using measurements at two frequencies, namely at the K-band with 24.5 GHz and at the Ka-band with 32.7 GHz frequency. Each satellite (A and B) records two phase measurements ($\varPhi _A^{K}$, $\varPhi _A^{Ka}$, $\varPhi _B^{K}$ and $\varPhi _B^{Ka}$) using heterodyne interferometry and phase tracking, which represent the phase difference between a local (LO) and a received (RX) electromagnetic field at reception time $t_\mathrm{r}$, i.e. (Kim (2000), eq. 2.14)

$$\begin{aligned}&\varPhi _{B}^{K/Ka}(t_\mathrm{r}) = \varPhi _\mathrm{Br}^{K/Ka} = -\left( \varphi _\text {RX,B} - \varphi _\text {LO,B} \right) \nonumber \\&\quad = - \left( {\hat{f}}_{A}^{K/Ka} \cdot \tau _A^\text {USO}(t_\mathrm{r}-\varDelta t_\mathrm{AeBr}^{K/Ka}) - {\hat{f}}_{B}^{K/Ka} \cdot \tau _B^\text {USO}(t_\mathrm{r}) \right) \end{aligned}$$

(36)

$$\begin{aligned}&\quad \approx -\left( {\hat{f}}_{A}^{K/Ka} \cdot \tau _A^\text {USO}(t_\mathrm{r}) - {\hat{f}}_{B}^{K/Ka} \cdot \tau _B^\text {USO}(t_\mathrm{r})\right) \nonumber \\&\qquad + {\hat{f}}_{A}^{K/Ka} \cdot \frac{\text {d} \tau _A^\text {USO}}{\text {d}t}\cdot \varDelta t_\mathrm{AeBr}^{K/Ka} + \text {const.} \end{aligned}$$

(37)

$$\begin{aligned}&\quad = -\left( {\hat{f}}_{A}^{K/Ka} \cdot \tau _A^\text {USO}(t_\mathrm{r}) {-} {\hat{f}}_{B}^{K/Ka} \cdot \tau _B^\text {USO}(t_\mathrm{r})\right) \nonumber \\&\qquad + f_{A}^{K/Ka}(t_\mathrm{r})\cdot \varDelta t_\mathrm{AeBr}^{K/Ka} {+} \text {const.} \end{aligned}$$

(38)

$$\begin{aligned}&\varPhi _{A}^{K/Ka}(t_\mathrm{r}) = \varPhi _{Ar}^{K/Ka} = +\left( \varphi _\text {RX,A} - \varphi _\text {LO,A} \right) \nonumber \\&\quad = + \left( {\hat{f}}_{B}^{K/Ka} \cdot \tau _B^\text {USO}(t_\mathrm{r}-\varDelta t_{BeAr}^{K/Ka}) - {\hat{f}}_{A}^{K/Ka} \cdot \tau _A^\text {USO}(t_\mathrm{r}) \right) \end{aligned}$$

(39)

$$\begin{aligned}&\quad \approx +\left( {\hat{f}}_{B}^{K/Ka} \cdot \tau _B^\text {USO}(t_\mathrm{r}) - {\hat{f}}_{A}^{K/Ka} \cdot \tau _A^\text {USO}(t_\mathrm{r})\right) \nonumber \\&\qquad - {\hat{f}}_{B}^{K/Ka} \cdot \frac{\text {d} \tau _B^\text {USO}}{\text {d}t}\cdot \varDelta t_{BeAr}^{K/Ka} + \text {const.} \end{aligned}$$

(40)

$$\begin{aligned}&\quad = +\left( {\hat{f}}_{B}^{K/Ka} \cdot \tau _B^\text {USO}(t_\mathrm{r}) - {\hat{f}}_{A}^{K/Ka} \cdot \tau _A^\text {USO}(t_\mathrm{r})\right) \nonumber \\&\qquad - f_{B}^{K/Ka}(t_\mathrm{r}) \cdot \varDelta t_{BeAr}^{K/Ka} + \text {const.} \end{aligned}$$

(41)

The phases $\varphi _{\ldots }$ of the electro-magnetic fields are given as the product of a static nominal frequency ${\hat{f}}_{A/B}^{K/Ka}$ and USO time $\tau _{A/B}^\text {USO}$, which differs from the proper time $\tau _{A/B}$ due to clock errors. These clock errors account for noise and errors sources, in particular for deviations of the USO frequency from the nominal or design values: ${\hat{f}}_{A}^{K} = 5076 \cdot 4.832\,\text {MHz}$, ${\hat{f}}_{A}^{Ka} = 6768 \cdot 4.832\,\text {MHz}$, ${\hat{f}}_{B}^{K} = 5076 \cdot 4.832099\,\text {MHz}$ and ${\hat{f}}_{B}^{Ka} = 6768 \cdot 4.832099\,\text {MHz}$Wen et al. (2019). The clock errors can be estimated during precise orbit determination (see CLK1B and USO1B data products in GRACE-FO) and allow to derive the apparent frequencies $f_{A/B}(t) = {\hat{f}}_{A/B} \cdot \text {d} \tau _{A/B}^\text {USO} / \text {d}t$, which are relevant for the ranging and contain effects from relativistic time dilation and clock errors, e.g. USO frequency deviations. For the purpose of calculating the light-time-correction, which is significantly smaller than the actual ranging signal, it is usually sufficient to drop the time-dependency and use a (daily) mean value $\langle f_{A/B} \rangle$, since the deviations of $f_{A/B}(t)/\langle f_{A/B} \rangle$ from unity are below $10^{-10}$ in magnitude for both, the daily clock drifts and the relativistic modulation.²

The first part of $\varPhi$ in line (38) and (41) is proportional to ${\hat{f}}_B \cdot \tau _B^\text {USO} (t_\mathrm{r})- {\hat{f}}_A \cdot \tau _A^\text {USO}(t_\mathrm{r})$ and describes a constant positive phase ramp with a slope of approx. 500 kHz and 670 kHz for the K- and Ka-band, respectively. The frequency order is reversed between the spacecraft. Usually, phase trackers are not aware of the frequency order and return a positive slope, which means that the sign of the second term ($\varDelta t$) is reversed between both S/C. This sign convention is consistent with the usual description of phase-tracking in the laser ranging instrument (see next section). However, it is opposite to the usual literature for microwave ranging (see Wen et al. (2019)), where the phase ramps on satellite A (GFO-C) have negative slope. The term $\varDelta t_\mathrm{AeBr}$ in above eq. describes the propagation time of the microwaves from satellite A to B, while $\varDelta t_{BeAr}$ denotes the opposite path. The last summand const. represents the fact that the phase measurement always have an unknown bias, which is constant unless the phase-tracking is interrupted or cycle slips occur. The MWI measures distance variations between the antenna phase center (APC), which are offset on each satellite by approx. 1.4 m in the direction of the distant satellite.

By subtracting the two phase observations in the K- or Ka-band, and dividing with the sum of the measured apparent frequencies $f_{A/B,\text {meas}}^{K/Ka}$ (cf. Kim (2000), eq. 2.16), one can obtain a range observation at the K- and Ka-band, i.e

$$\begin{aligned} \rho _\text {DOWR}^{K/Ka}(t)&= c_0 \cdot \int ^t_0 \frac{ \text {d}~\left( \varPhi _{Br}^{K/Ka}(t^\prime ) - \varPhi _{Ar}^{K/Ka}(t^\prime ) \right) / \text {d}t^\prime }{f_{A,\text {meas}}^{K/Ka}(t^\prime ) + f_{B,\text {meas}}^{K/Ka}(t^\prime )} \text {d}t^\prime \end{aligned}$$

(42)

$$\begin{aligned}&\approx c_0 \cdot \frac{ \varPhi _{Br}^{K/Ka} - \varPhi _{Ar}^{K/Ka}}{\langle f_{A,\text {meas}}^{K/Ka}\rangle + \langle f_{B,\text {meas}}^{K/Ka} \rangle } \nonumber \\&= c_0 \cdot \frac{ f_{A}^{K/Ka}(t) \cdot \varDelta t_\mathrm{AeBr}^{K/Ka} + f_{B}^{K/Ka}(t) \cdot \varDelta t_{BeAr}^{K/Ka}}{ \langle f_{A,\text {meas}}^{K/Ka} \rangle + \langle f_{B,\text {meas}}^{K/Ka}\rangle } + \text {const.} \end{aligned}$$

(43)

$$\begin{aligned}&\approx c_0 \cdot \varDelta t_\text {inst,APC} \nonumber \\&\quad + c_0 \cdot \frac{ \langle f_A^{K/Ka} \rangle \cdot {\mathcal {T}}_\mathrm{AeBr}^{K/Ka} + \langle f_B^{K/Ka} \rangle \cdot {\mathcal {T}}_{BeAr}^{K/Ka}}{\langle f_A^{K/Ka} \rangle + \langle f_B^{K/Ka} \rangle } \nonumber \\&\quad + c_0 \cdot \frac{ \langle f_A^{K/Ka} \rangle \cdot \varDelta t_\text {media}^{K/Ka} + \langle f_B^{K/Ka} \rangle \cdot \varDelta t_\text {media}^{K/Ka}}{ \langle f_A^{K/Ka} \rangle + \langle f_B^{K/Ka}\rangle } \nonumber \\&\quad + \text {const.} \end{aligned}$$

(44)

$$\begin{aligned}&= \rho _\text {inst,APC} + c_0 \cdot {\mathcal {T}}_\text {DOWR}^{K/Ka} + \rho _\text {media}^{K/Ka} + \text {const.}, \end{aligned}$$

(45)

which can be written as the sum of instantaneous distance between APC $\rho _\text {inst,APC}$, light time effect ${\mathcal {T}}_\text {DOWR}^{K/Ka}$ and ionospheric delay $\rho _\text {media}^{K/Ka}$. The light paths in the DOWR scheme are shown for a single frequency in the left plot of Fig. 1. Equation (42) is suited to convert the measured phases to the DOWR ranges $\rho _\text {DOWR}^{K/Ka}$. For the derivation of the much smaller light-time and ionospheric corrections, the approximations in Eqs. (43)–(45) are usually sufficient, where the distinction between true apparent and measured apparent frequency, as well as their time-dependencies, are omitted.

One can remove the ionospheric effect by a linear combination of $\rho _\text {DOWR}^{K}$ and $\rho _\text {DOWR}^{Ka}$, which yields the DOWR biased range as

$$\begin{aligned} \rho _\text {DOWR}&= a^{Ka} \cdot \rho _\text {DOWR}^{Ka} + a^{K} \cdot \rho _\text {DOWR}^{K} \nonumber \\&= \rho _\text {inst,APC} + c_0 \cdot {\mathcal {T}}_\text {DOWR} + \text {const.} \end{aligned}$$

(46)

where the light-time effect ${\mathcal {T}}_\text {DOWR}$ is, in general, a function of four ${\mathcal {T}}^{K/Ka}_{\ldots }$ terms arising from two photon paths at two frequencies:

$$\begin{aligned} {\mathcal {T}}_\text {DOWR}&= a^\text {K} \cdot {\mathcal {T}}_\text {DOWR}^{K} + a^\text {Ka} \cdot {\mathcal {T}}_\text {DOWR}^{Ka} \end{aligned}$$

(47)

$$\begin{aligned}&= b_\mathrm{AeBr}^{K} \cdot {\mathcal {T}}_\mathrm{AeBr}^K + b_\mathrm{AeBr}^{Ka} \cdot {\mathcal {T}}_\mathrm{AeBr}^{Ka} + b_{BeAr}^{K} \cdot {\mathcal {T}}_{BeAr}^{K} \nonumber \\&\quad + b_{BeAr}^{Ka} \cdot {\mathcal {T}}_{BeAr}^{Ka} \end{aligned}$$

(48)

with $a^{K/Ka}_{\ldots }$ and $b^{K/Ka}_{\ldots }$ coefficients given in Table 2.

The biased dual-one way range $\rho _\text {DOWR}$ is apportioned in Eq. (46) into the instantaneous range $\rho _\text {inst,APC}$ and an effect due to the finite speed of light $c_0 \cdot {\mathcal {T}}_\text {DOWR}$. In order to obtain the instantaneous range, one has to remove this light-time effect using an estimate or correction $\widehat{{\mathcal {T}}}_\text {DOWR}$, which can be derived from orbit data. Moreover, an antenna offset correction is applied in order to transform the biased range between APC into a biased range between the center-of-mass that is usually used for gravity field recovery.

The cross coupling of $\varDelta t_\text {media}$ into the light-time correction ${\mathcal {T}}_\text {DOWR}$ is usually omitted (cf. Eq. (34)), i.e. the K and Ka superscripts of ${\mathcal {T}}$ are dropped

$$\begin{aligned} \widehat{{\mathcal {T}}}_\text {DOWR} \approx b_\mathrm{AeBr} \cdot {\mathcal {T}}_\mathrm{AeBr} + b_{BeAr} \cdot {\mathcal {T}}_{BeAr}, \end{aligned}$$

(49)

because the absolute value of the ionospheric delay $\varDelta t_\text {media}$ is difficult to estimate and the effect on the final correction ${\mathcal {T}}_\text {DOWR}$ is well below the microwave instrument resolution. In other words, the LTC computation neglects any atmospheric effect, i.e. the photons at K- and Ka-Band have the same emission time as in vacuum and travel along the same path. However, this approximation does not affect the phase delay as determined and corrected for with the ionospheric correction (cf. Sect. (4)). The omission error in the LTC is at the sub-picometer level and can be assessed using Eqs. (28) and  (34), i.e. 

$$\begin{aligned} \left| c_0 {\mathcal {T}}_\text {DOWR,media} \right|= & {} \left| \frac{40.3\,\text {Hz}^2}{c_0} \cdot \frac{\text {TEC}}{1\,e^-/\text {m}^3}\right. \nonumber \\&\left. \cdot \left( \frac{ b_\mathrm{AeBr}^{K} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_A}{(f_A^K)^2} - \frac{ b_{BeAr}^{K} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_B}{(f_B^K)^2}\right. \right. \nonumber \\&\left. \left. \quad + \frac{b_\mathrm{AeBr}^{Ka} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_A}{(f_A^{Ka})^2} - \frac{b_{BeAr}^{{Ka}} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_B}{(f_B^{Ka})^2} \right) \right| \end{aligned}$$

(50)

$$\begin{aligned}\approx & {} \left| -2 \cdot 10^{-13}\,\text {m} - 8 \cdot 10^{-18}\,\text {m} \cdot \frac{{{\dot{\rho }}}_\text {inst}}{1\,\text {m/s}} \right| < 10^{-12}\,\text {m}\nonumber \\ \end{aligned}$$

(51)

where $\mathbf {d}_0 = (\mathbf {r}_B-\mathbf {r}_A)/|\mathbf {r}_B-\mathbf {r}_A|$, $\mathbf {d}_0 .\dot{\mathbf {r}}_A = -7.6\,\text {km/s}$, and $\mathbf {d}_0 .\dot{\mathbf {r}}_B = 7.6\,\text {km/s}+{{\dot{\rho }}}_\text {inst}$ were used as values. The range rate ${{\dot{\rho }}}_\text {inst}$ is usually below 1 m/s, hence, the modulation due to $\rho _\text {inst}$ is insignificant. The same holds for variations of the TEC, which can be expected to be well below the used upper bound estimate $\text {TEC} = 10^{12}\,\text {e}^-/\text {m}^3 \cdot 200\,\text {km}$.

The leading terms of the DOWR light-time correction in the range domain, which has to be subtracted from the measured biased range $\rho _\text {DOWR}$ to obtain the instantaneous range, reads

$$\begin{aligned} c_0 \widehat{{\mathcal {T}}}_\text {DOWR}&= \varDelta t_\text {inst} \cdot \left( b_\mathrm{AeBr} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_A-b_{BeAr} \cdot \mathbf {d}_0 .\dot{\mathbf {r}}_B \right) + \text {const.} + \ldots \nonumber \\&= -\frac{ |\mathbf {r}_B-\mathbf {r}_A| \cdot {\dot{\rho }}_\text {inst,OD}}{2 \cdot c_0} + \text {const.} + \cdots , \end{aligned}$$

(52)

where both shown terms have a typical magnitude of a few hundred micrometers (cf. Table 5). The ${\dot{\rho }}_\text {inst,OD}$ denotes the instantaneous range rate from orbit data (OD). This leading term describes approximately 99.9 % of the LTC at once and twice the orbit frequency, which may be sufficient in some cases. However, the analyses in this paper consider the full expression, not just the leading term.

Table 2 Numerical values for coefficients introduced to describe the light time correction in dual one-way ranging, which are based on carrier frequencies in the K and Ka band for the microwave ranging system

7 Light time correction in two-way ranging (TWR)

The laser ranging instrument aboard GRACE-Follow-On is based on a master-transponder scheme, which is also called a two way ranging scheme. The role of master and transponder is inter-changeable between the satellites. As shown on the right plot in Fig. 1, the master satellite emits a photon at event $M_\mathrm{e}$ using a frequency-stabilized laser source. The optical phase (in cycles) of this photon can be modelled as a function of the coordinate time t

$$\begin{aligned} \varphi _\text {M}(t) = \int _{0}^{t} {\tilde{f}}_M(t^\prime )\cdot \frac{\text {d}\tau _M(t^\prime )}{\text {d}t^\prime } ~\text {d}t^\prime \end{aligned}$$

(53)

where ${\tilde{f}}_M$ is the instantaneous optical laser frequency that would be measured in a rest-frame at the laser source and $\tau _M$ refers to the proper time of the master satellite. Imperfections of the laser or cavity, i.e. frequency variations, can be accounted for by the time-dependent ${\tilde{f}}_M$.

The photon emitted by the master satellite propagates to the transponder craft. The transponder utilizes a frequency-locked loop with 10 MHz frequency offset. This means the laser phase $\varphi _\text {LO,T}(t)$, more precisely the time-derivative of it, is controlled such that the beatnote phase $\varPhi _{T}(t)$, given as the phase difference between received (RX) and local oscillator (LO) light, becomes

$$\begin{aligned} \varPhi _{T}(t)&= \varphi _\text {LO,T}-\varphi _\text {RX,T} = \varphi _\text {LO,T}(t)-\varphi _\text {M}(t-\varDelta t_\text {MeTp}(t)) \nonumber \\&= +10\,\text {MHz} \cdot \tau _T^\text {USO}(t) + \varphi _\epsilon (t) + \text {const.}, \end{aligned}$$

(54)

where $\tau _M^\text {USO}$ is the time of the ultra-stable oscillator clock, which may differ from the proper time $\tau _M$ due to noise or errors sources. The beatnote phase $\varPhi _{T}$ implies that the optical phase of the transponder laser with units of cycles is

$$\begin{aligned} \varphi _\text {LO,T}(t) = \varphi _\text {M}(t-\varDelta t_\text {MeTp}(t)) + 10\,\text {MHz} \cdot \tau _T^\text {USO}(t) + \varphi _\epsilon (t) + \text {const.}, \end{aligned}$$

(55)

where $\varphi _\epsilon (t)$ was used to account for phase-variations that were not fully suppressed by the feedback control loop, e.g. due to finite gain and bandwidth. These are much smaller than the phase ramp with a slope of 10 MHz. The loop ensures a constant phase relation between emitted and received light on the transponder side, in other words, the transponder seems to reflect the received light at event $T_p$, however, with enhanced light power and slightly different frequency.

Eventually, the transponder photon returns to the master side at the reception event $M_\mathrm{r}$. The phase of the beatnote on the master satellite $\varPhi _{M}$ reads

$$\begin{aligned} \varPhi _{M}(t_\mathrm{r})&= \varphi _\text {RX,M} - \varphi _\text {LO,M} = \varphi _\text {LO,T}(t_\mathrm{r}-\varDelta t_\text {TpMr}) - \varphi _M(t_\mathrm{r}) \end{aligned}$$

(56)

$$\begin{aligned}&= \varphi _\text {M}(t_\mathrm{r}-\varDelta t_\text {TpMr}-\varDelta t_\text {MeTp}) - \varphi _M(t_\mathrm{r}) \nonumber \\&\quad + 10\,\text {MHz} \cdot \tau _T^\text {USO}(t_\mathrm{r}-\varDelta t_\text {TpMr}) + \varphi _\epsilon (t_\mathrm{r}-\varDelta t_\text {TpMr}) + \text {const.} \end{aligned}$$

(57)

$$\begin{aligned}&\approx -\frac{d \varphi _\text {M}}{d \tau _M} \cdot \frac{d \tau _M}{d t} \cdot (\varDelta t_\text {TpMr}+\varDelta t_\text {MeTp}) + 10\,\text {MHz} \nonumber \\&\quad \cdot \tau _T^\text {USO}(t_\mathrm{r}-\varDelta t_\text {TpMr}) + \varphi _\epsilon (t_\mathrm{r}-\varDelta t_\text {TpMr}) + \text {const.} \end{aligned}$$

(58)

$$\begin{aligned}&= -f_M(t_\mathrm{r}) \cdot (\varDelta t_\text {TpMr}+\varDelta t_\text {MeTp}) + 10\,\text {MHz}\nonumber \\&\quad \cdot \tau _T^\text {USO}(t_\mathrm{r}-\varDelta t_\text {TpMr}) + \varphi _\epsilon (t_\mathrm{r}-\varDelta t_\text {TpMr}) + \text {const.} \end{aligned}$$

(59)

The ranging information is encoded in the term containing the product of true apparent optical frequency ($f_M = {\tilde{f}}_M \cdot \text {d}\tau _M/\text {d}t$ ) and photon time of flight $\varDelta t_{\ldots }$. It can give rise to Doppler shifts of up to a few  MHz over one orbital revolution.

Subtracting both phase observations, when the transponder phase is temporally aligned to the master using an estimated one-way light travel time $\varDelta t_\text {TpMr,est}$, removes the 10 MHz phase ramp and the phase residuals $\varphi _\epsilon$. Then, the phase difference is converted to a biased range observable using an estimate of the apparent optical frequency³$f_{M,\text {est}}(t)$, as in the DOWR case (cf. Eq. (42)), i.e.

$$\begin{aligned} \rho _\text {TWR}(t)&= c_0 \cdot \int _0^t \frac{ \text {d} \left( \varPhi _{T}(t^\prime -\varDelta t_\text {TpMr,est}) - \varPhi _{M}(t^\prime ) \right) / \text {d}t^\prime }{2 \cdot f_{M,\text {est}}(t^\prime )} \text {d}t^\prime \end{aligned}$$

(60)

$$\begin{aligned}&\approx c_0 \cdot \frac{ (\langle f_{M} \rangle + \delta f_{M}(t)) \cdot (\varDelta t_\text {TpMr}+\varDelta t_\text {MeTp}) }{2 \cdot \langle f_{M,\text {est}} \rangle } + \text {const.} \end{aligned}$$

(61)

$$\begin{aligned}&= \left( 1 + \frac{\langle f_{M} \rangle -\langle f_{M,\text {est}} \rangle }{\langle f_{M,\text {est}} \rangle } + \frac{\delta f_{M}(t)}{\langle f_{M,\text {est}} \rangle } \right) \nonumber \\&\quad \cdot \frac{c_0 \cdot (2 \cdot \varDelta t_\text {inst} + {\mathcal {T}}_\text {MeTp} + {\mathcal {T}}_\text {TpMr})}{2} + \text {const.} \end{aligned}$$

(62)

$$\begin{aligned}&= \left( 1 + \kappa + \delta \kappa (t) \right) \cdot \left( \rho _\text {inst} + \frac{ {\mathcal {T}}_\text {MeTp} + {\mathcal {T}}_\text {TpMr}}{2} \right) + \text {const.} \end{aligned}$$

(63)

$$\begin{aligned}&= \rho _\text {inst}(t) + c_0 {{\mathcal {T}}}_\text {TWR}(t) + \kappa \cdot \rho _\text {inst} + \delta \kappa (t) \cdot \rho _\text {inst}(t)\nonumber \\&\quad + (\kappa + \delta \kappa (t)) \cdot c_0 {{\mathcal {T}}}_\text {TWR}(t) + \text {const.} \end{aligned}$$

(64)

The precise Eq. (60) can be used to convert the phase observables to a non-instantaneous biased range $\rho _\text {TWR}$, even with a time-dependent frequency estimate $f_{M,\text {est}}(t)$. Under the assumption of a static estimate $\langle f_{M,\text {est}} \rangle$, and with Eqs. (59), (54) and $\varDelta t_\text {TpMr,est} \approx \varDelta t_\text {TpMr}$, the expression can be approximated as Eq. (64), which illustrates the coupling of frequency errors and the light-time correction effect. The first terms are the instantaneous range $\rho _\text {inst}$ and the light-time correction ${{\mathcal {T}}}_\text {TWR} = ({\mathcal {T}}_\text {MeTp} + {\mathcal {T}}_\text {TpMr})/2$, respectively. The third term describes a static scale factor error $\kappa = (\langle f_{M} \rangle -\langle f_{M,\text {est}} \rangle )/\langle f_{M,\text {est}} \rangle$ in the conversion from phase to range, while the term proportional to $\delta \kappa = \delta f_{M}(t)/\langle f_{M,\text {est}} \rangle$ accounts for laser phase variations, commonly known as laser frequency noise Abich et al. (2019). The coupling of $\kappa$ or $\delta \kappa$ with the LTC in the fifth term is negligible compared to the same coupling with $\rho _\text {inst}$, because the magnitude of $c_0{\mathcal {T}}_\text {TWR}$ is below the millimeter level (cf. Table 5). The relevant aspect for the following sections is that the final Euclidean biased range can be computed as $\rho _\text {inst,TWR}=\rho _\text {TWR}-c_0 {{\mathcal {T}}}_\text {TWR}$.

In order to compute the propagation time $\varDelta t_\text {MeTp}$ from the master emission event (Me on right plot of Fig. 1) to the transponder reception (Tp in Fig. 1), the result of $\varDelta t_\text {TpMr}$ is needed, as apparent from the following iterative equation

$$\begin{aligned} \varDelta t_\text {MeTp}^{(n+1)}(t_\mathrm{r})&= \frac{| \mathbf {r}_T(t_\mathrm{r}-\varDelta t_\text {TpMr}) - \mathbf {r}_M(t_\mathrm{r} - \varDelta t_\text {MeTp}^{(n)} - \varDelta t_{TpMr}) |}{c_0} \nonumber \\&\quad + {\mathcal {T}}_\text {GR}(\mathbf {r}_\mathrm{r} {=} \mathbf {r}_T(t_\mathrm{r}{-}\varDelta t_\text {TpMr}), \mathbf {r}_e {=} \mathbf {r}_M(t_\mathrm{r} {-} \varDelta t_\text {MeTp}^{(n)} {-} \varDelta t_\text {TpMr})) \end{aligned}$$

(65)

which we rigorously approximate, with the same approach as utilized for Eq. (33), as

$$\begin{aligned} \varDelta t_\text {MeTp}= & {} \varDelta t_\text {inst} + {\mathcal {T}}_\text {SR,MeTp} + {\mathcal {T}}_\text {GR,MeTp} \end{aligned}$$

(66)

$$\begin{aligned} {\mathcal {T}}_\text {MeTp}= & {} \frac{\varDelta t_\text {inst} (\mathbf {d}_0.\dot{\mathbf {r}}_T-2 \mathbf {d}_0.\dot{\mathbf {r}}_M)+\varDelta t_\text {inst}^2 \left( 2 \mathbf {d}_0.\ddot{\mathbf {r}}_M-\frac{\mathbf {d}_0.\ddot{\mathbf {r}}_T}{2}\right) }{c_0}\nonumber \\&+\frac{\varDelta t_\text {inst} \left( \left| \dot{\mathbf {r}}_T-2 \dot{\mathbf {r}}_M\right| ^2+(\mathbf {d}_0.\dot{\mathbf {r}}_T)^2\right) }{2 c_0^2} \nonumber \\&+\frac{\varDelta t_\text {inst}^2 \left( -2 \mathbf {d}_0.\ddot{\mathbf {r}}_M (\mathbf {d}_0.\dot{\mathbf {r}}_M-\mathbf {d}_0.\dot{\mathbf {r}}_T)-4 \dot{\mathbf {r}}_M.\ddot{\mathbf {r}}_M+2 \dot{\mathbf {r}}_T.\ddot{\mathbf {r}}_M-\mathbf {d}_0.\ddot{\mathbf {r}}_T \cdot \mathbf {d}_0.\dot{\mathbf {r}}_T+\ddot{\mathbf {r}}_T.\dot{\mathbf {r}}_M -\dot{\mathbf {r}}_T.\ddot{\mathbf {r}}_T/2\right) }{c_0^2} \nonumber \\&+\frac{\varDelta t_\text {inst} \left( \mathbf {d}_0.\dot{\mathbf {r}}_T \left( 2 \left| \dot{\mathbf {r}}_M\right| ^2+\left| \dot{\mathbf {r}}_T\right| ^2-2 \dot{\mathbf {r}}_T.\dot{\mathbf {r}}_M\right) -2 \left| \dot{\mathbf {r}}_M\right| ^2 \mathbf {d}_0.\dot{\mathbf {r}}_M\right) }{c_0^3} \end{aligned}$$

(67)

$$\begin{aligned}&+ \frac{{\mathcal {T}}_\text {GR,TpMr} \cdot \mathbf {d}_0.\dot{\mathbf {r}}_T - ({\mathcal {T}}_\text {GR,TpMr}+{\mathcal {T}}_\text {GR,MeTp}) \cdot \mathbf {d}_0.\dot{\mathbf {r}}_M}{c_0}+{\mathcal {O}}(10^{-12}\,\text {m}/c_0). \end{aligned}$$

(68)

The satellite state vectors, $\varDelta t_\text {inst}$ and $\mathbf {d}_0 = (\mathbf {r}_M - \mathbf {r}_T)/|\mathbf {r}_M - \mathbf {r}_T|$ are evaluated at the reception time ($t_\mathrm{r}$) and are the same as those needed to compute ${\mathcal {T}}_\text {TpMr}$ with Eq. (34). The delay due to the atmosphere $\varDelta t_\text {media}$ was omitted. The general relativistic contributions ${\mathcal {T}}_\text {GR} = {\mathcal {T}}_\text {PM} + {\mathcal {T}}_\text {HM} + {\mathcal {T}}_\text {SM}$ are evaluated at

$$\begin{aligned} {\mathcal {T}}_\text {GR,TpMr}&= {\mathcal {T}}_\text {GR}(\mathbf {r}_\mathrm{r} = \mathbf {r}_M(t_\mathrm{r}),~\mathbf {r}_e \nonumber \\&= \mathbf {r}_T(t_\mathrm{r} - \varDelta t_\text {inst} - \varDelta t_\text {inst} \mathbf {d}_0.\dot{\mathbf {r}}_T/c_0) ) \end{aligned}$$

(69)

$$\begin{aligned} {\mathcal {T}}_\text {GR,MeTp}&= {\mathcal {T}}_\text {GR}\left( \mathbf {r}_\mathrm{r} = \mathbf {r}_T(t_\mathrm{r} - \varDelta t_\text {inst} - \varDelta t_\text {inst} \mathbf {d}_0.\dot{\mathbf {r}}_T/c_0 ),~\mathbf {r}_e\right. \nonumber \\&\left. = \mathbf {r}_M\left( t_\mathrm{r} - \varDelta t_\text {inst} \cdot \left( 2 c_0 + \mathbf {d}_0.\dot{\mathbf {r}}_T - \mathbf {d}_0.\dot{\mathbf {r}}_M \right) /c_0 \right) \right) , \end{aligned}$$

(70)

with the help of the Taylor expansion in Eq. (30).

It is noteworthy that the leading term in the TWR light-time correction

$$\begin{aligned} c_0 \widehat{{\mathcal {T}}}_\text {TWR}&= c_0 \frac{{\mathcal {T}}_\text {MeTp} + {\mathcal {T}}_\text {TpMr}}{2} = -\frac{ | \mathbf {r}_B - \mathbf {r}_A| \cdot {\dot{\rho }}_\text {inst,OD}}{ c_0} + \text {const.} + \ldots \end{aligned}$$

(71)

differs by a factor of two compared to the DOWR correction (cf. Eq. (52)), whereby the static part has a similar magnitude (cf. Table 5).

8 Requirements on light time correction precision

It is sensible to require that the light time corrections $c_0 {\mathcal {T}}_\text {TWR}$ and $c_0 {\mathcal {T}}_\text {DOWR}$ are precise enough to not limit the precision of the instantaneous range, which is the measured biased range with subtracted light time correction. The precision of the instantaneous range $\rho _\text {inst}$ should ideally be limited by instrument noise and errors. Noise is driven by stochastic processes and can be described with spectral densities in the frequency domain. For instance, the noise requirement for the laser ranging instrument on GRACE FO is defined in terms of the amplitude spectral density (ASD), which is the square root of the power spectral density, as Abich et al. (2019)

$$\begin{aligned}&\text {ASD}[\rho _\text {LRI,req}] = 80\,\frac{\text {nm}}{\sqrt{\text {Hz}} } \sqrt{1+\left( \frac{3\,\text {mHz}}{f}\right) ^2} \sqrt{1+\left( \frac{10\,\text {mHz}}{f}\right) ^2},\nonumber \\&\quad 2\,\text {mHz}\le f \le 100\,\text {mHz} \end{aligned}$$

(72)

while the corresponding requirement of the MWI reads Kornfeld et al. (2019)

$$\begin{aligned} \text {ASD}[\rho _\text {KBR}] = 2.62\,\frac{\upmu \text {m}}{\sqrt{\text {Hz}} } \sqrt{1+\left( \frac{3\,\text {mHz}}{f}\right) ^2}. \end{aligned}$$

(73)

Deterministic or systematic errors manifest often as sinusoidal variations, so called tone errors. These should not exceed $\delta \rho = 1\,\upmu \text {m}$ peak amplitude in GRACE FO measurements. This value is specified for the MWI at twice the orbital frequency ($f = 2 f_\text {orb} \approx 0.35\,\text {mHz}$) and for the LRI between $10 f_\text {orb} \le f \le 200 f_\text {orb}$ Kornfeld et al. (2019). The $2f_\text {orb}$ KBR requirement is likely inherited and adopted from the GRACE mission (Stanton et al. (1998), p. 23), while the higher LRI requirement band ($10 f_\text {orb} .. 200 f_\text {orb}$) could be justified by the fact that other error sources like accelerometer or background model deficiencies limit the gravity field accuracy at lower frequencies. The authors recommend that both requirements are revised in future missions. Although not strictly specified by the instruments, it is reasonable to require that the LTC has no sinusoidal errors above $1\,\upmu \text {m}$ magnitude for all frequencies.

Table 3 The mean value and peak amplitudes at once and twice the orbital frequency ($f_\text {orb} =0.18$ mHz) of the difference $\text {GRA\_KBR1B\_LTC}-c_0 {\mathcal {T}}_\text {DOWR}$, where ${\mathcal {T}}_\text {DOWR}$ is computed with different accuracy levels

In the next sections, we illustrate the frequency content of time-domain signals with ASD plots, where the y-axis has units of $\mathrm{m} / \sqrt{\text {Hz}}$. These plots show the peak amplitude $\delta \rho$ of a sinusoidal variation with an amplitude of

$$\begin{aligned} \frac{\delta \rho }{\sqrt{2} \sqrt{\text {ENBW}}}, \end{aligned}$$

(74)

where ENBW is the equivalent noise bandwidth with units of Hertz. The ENBW depends on many parameters such as the length of the time-series, the sampling rate and the window function Heinzel et al. (2002). Since many gravity field recovery methods are using range rates, we recall that ASD values at a Fourier frequency f with units of $\mathrm{m} / \sqrt{\text {Hz}}$ can be converted into the range rate domain with units $\mathrm{m} / (\text {s}\sqrt{\text {Hz}})$ by a multiplication with $2 \pi f$.

The actual in-orbit ASD of the LRI is well below the 80 nm/$\sqrt{\text {Hz}}$ requirement as shown in Abich et al. (2019), i.e.

$$\begin{aligned} \text {ASD}[\rho _\text {LRI}]=\left\{ \begin{array}{@{}ll@{}} 15\,\text {nm}/\sqrt{\text {Hz}}, &{}\quad f = 35\,\text {mHz} \\ 0.3\,\text {nm}/\sqrt{\text {Hz}}, &{}\quad f = 0.85\,\text {Hz} \end{array}\right. \end{aligned}$$

(75)

which imposes stricter goals for the LTC precision at high frequencies.

9 Validation of the analytical approximations for $\varDelta t$

In order to verify the equations for the light propagation time and our implementation of the software code, we performed a closed-loop (i.e. backward-forward) simulation using reduced-dynamic orbit data of both GRACE Follow-On satellites in the Geocentric Celestial Reference Frame (GCRF) from 5th February 2019 (GNI1B Release 04). One of the two satellites is designated as receiver with position $\mathbf {r}(t_\mathrm{r})$. At each epoch $t_\mathrm{r}$ of the data, which has a sampling rate of 1 Hz, the light propagation time $\varDelta t = \varDelta t_\text {inst} + {\mathcal {T}}_\text {SR} + {\mathcal {T}}_\text {GR}$ between the satellites is computed according to Eqs. (33)–(35), which make use of Eqs. (20)–(25). With the propagation time $\varDelta t$, we compute the photon emission position $\mathbf {r}_e$ and emission time $t_\mathrm{r} - \varDelta t$. Afterwards, we determine the vectorial coordinate speed of light $c_n \cdot \mathbf {d}_0$ pointing to the receiver (Eqs. (13), (19)), which serves as the initial condition for a numerical integration of the equations of motion for photons (Eq. (9)) using the Adams–Bashforth–Moulton method Shampine and Reichelt (1997). The metric tensor used is based on a high-fidelity geopotential field, computed according to the models shown in Table 1, and takes into account the vector potential due to Earth’s spin. The integration is performed for a duration of $\varDelta t$, which yields the photon path with an end position $\mathbf {r}_\mathrm{r}^\prime$. If the analytical expressions to compute $\varDelta t$ are correct, $\mathbf {r}_\mathrm{r}^\prime$ and $\mathbf {r}_\mathrm{r}$ should be identical. Hence, we define the error $\epsilon$ in the analytically-derived $\varDelta t$ as

$$\begin{aligned} \epsilon = \left( \mathbf {r}_\mathrm{r}^\prime -\mathbf {r}_\mathrm{r}\right) .\frac{\dot{\mathbf {r}}_{r}^\prime }{|\dot{\mathbf {r}}_{r}^\prime |} \approx \left( \mathbf {r}_\mathrm{r}^\prime -\mathbf {r}_\mathrm{r}\right) .\mathbf {d}_0 \approx (\varDelta t^\prime - \varDelta t) \cdot c_0 \end{aligned}$$

(76)

which takes into consideration only the error in the propagation direction of the photon, since this contributes to the phase measurement in microwave or laser ranging. In other words, $\epsilon$ is the error of the computed $\varDelta t$ with respect to the more accurate $\varDelta t^\prime$.

The lateral part of the displacement $\mathbf {r}_\mathrm{r}^\prime -\mathbf {r}$ is of the order of 4 $\upmu \hbox {m}$ and arises due to the light bending (cf. Sect. 3), which has been omitted in our analytical approximation. By evaluating $\epsilon$, it can be shown that the bending—and omission of the bending in the analytical approximation—has a negligible effect on the phase measurement, since the longitudinal offset in propagation direction is very small and since the phasefront is, in good approximation, planar in the vicinity of $\mathbf {r}_\mathrm{r}^\prime$, i.e., the offset $\mathbf {r}_\mathrm{r}^\prime -\mathbf {r}$ vanishes when projected onto the propagation direction.

Due to the limited precision of double floating-point arithmetic, we perform the numerical integration in uniform co-moving coordinate frames, in order to have state vectors with small numerical values. This allows us to resolve even minor contributions within the light time correction.

Fig. 2

Amplitude spectral density plots of the model error $\epsilon$ and of the term ${\mathcal {T}}_\text {HM}$. a The first six traces show the model error $\epsilon$ for different contributors in the light time correction ${\mathcal {T}}$. The model error $\epsilon$ as a function of the number of sampling points N of the path integral (Eq. (22)) is shown in subfigure b, while the influence of the truncation degree for the SH expansion of the gravitational potential is illustrated in c. Subfigure d shows the ASD of a time series of ${\mathcal {T}}_\text {HM}$, where only a single gravitational potential model from Table 1 was used. All subfigures use the Nuttall4a window function. The equivalent peak height of a sinusoidal variation with 1 picometer amplitude is visualized as green dashed line in all four plots

The result of the one-way ranging validation, i.e. the time series of $\epsilon$, is shown in the spectral domain in Fig. 2a). The upper-most trace in red shows the error $\epsilon$, if special and general relativistic effects are omitted in the calculation of the light travel time $\varDelta t$, which means $\varDelta t = t_\text {inst}$. Considering ${\mathcal {T}}_\text {SR}$ yields the blue trace. The general relativistic contribution to the light propagation shows two sinusoidal variations at once and twice the orbital frequency and a continuous spectral content decaying towards higher frequencies. The peak at the orbital frequency is caused by the radially symmetric gravity field (${\mathcal {T}}_\text {PM}$), while the higher moments cause the twice per revolution peak and the continuous part.

Since the spectral plots conceal the DC component, the mean value of $\epsilon$ is provided in the legend. The figure confirms that the different contributions in the propagation time indeed reduce the error $\epsilon$ down to a mean level of $2.5 \cdot 10^{-13}\text {m}/c_0$, with fluctuations well below $1\,\text {pm}/\sqrt{\text {Hz}}/c_0$. The remaining peaks apparent at once and twice the orbital frequency from sinusoidal variations (tones) are not described properly with units of a spectral density plot (cf. Sect. 8). These variations have a peak magnitude in the time-domain of less than 1 picometer (green dashed line in Fig. 2a), if ${\mathcal {T}}_\text {PM}$ and ${\mathcal {T}}_\text {HM}$ are considered .

The contribution of the general relativistic correction ${\mathcal {T}}_\text {SM}$ due to Earth’s spin moment is present predominantly at once and twice the orbital frequency, but with a negligible magnitude (difference between brown and black trace). Hence, ${\mathcal {T}}_\text {SM}$ can be safely omitted from now on.

Table 4 The mean value and peak amplitudes at once and twice the orbital frequency ($f_\text {orb} =0.18$ mHz) of the differences $\text {GFO\_KBR1B\_LTC} - c_0 {\mathcal {T}}_\text {DOWR}$ and $\text {GFO\_LRI1B\_LTC} - c_0 {\mathcal {T}}_\text {TWR}$ for different accuracy levels of ${\mathcal {T}}_\text {DOWR}$ and ${\mathcal {T}}_\text {TWR}$

Table 5 The mean value and peak amplitudes at once and twice the orbital frequency ($f_\text {orb} =0.18$ mHz) of different constituents in the LTC ($c_0 \cdot {\mathcal {T}}$)

The dependence of the model error $\epsilon$ on the sampling point number N in Eq. (22) is shown in Fig. 2b), while Fig. 2c) visualizes the effect of the truncation degree for the SH expansion of the gravitational potential. The actual signal ${\mathcal {T}}_\text {HM}$ for different individual models of the gravitational potential (cf. Table 1) is depicted in Fig. 2d). In general, Fig. 2 can be used to decide which models and parameters are required for a particular accuracy level in the computation of the light time correction.

Although this section showed only one-way ranging results, most of the findings are also applicable for the TWR and DOWR combinations, since these are formed by the average of two one-way ranging results. Only ${\mathcal {T}}_\text {SM}$ and some terms in ${\mathcal {T}}_\text {SR}$ flip signs between the two opposite directions, which means they are canceling to a large extent in the TWR and DOWR case.

A result of this analysis is that the following parameters of ${\mathcal {T}}_\text {HM}$ are sufficient to meet the precision requirements formulated in Sects. 10 and 11: SH degree of the static gravity should be $\ge 50$, while a Solid Earth Tide (SET) model with degree 4 is sufficient; the path integral should be approximated with N$\ge 10$ and direct tidal accelerations should be taken into account at least from Sun and Moon.

10 Comparison with GRACE and GRACE FO light time correction

We compared the method to derive the light time correction presented herein with the light time correction values in the level-1b data of the GRACE and GRACE Follow-On missions. These values are provided in the KBR1B and LRI1B datasets alongside with the actual biased range. The most recent version of the GRACE data is release 03 (RL03), which is available only for the SCA1B and KBR1B data products, while for all other products RL02 is the most recent version PODAAC (2018). Details on the processing of GRACE data can be found in Wu et al. Wu et al. (2006). The GRACE Follow-On data is available in version RL04 by the time of the writing PODAAC (2019).

For the GRACE data, the GNV1B orbit data is rotated from the terrestrial to the celestial frame by a rotational matrix formed according to the IAU-2000 standard using Earth orientation parameteres Petit and Luzum (2010). The sampling rate of the orbit data is 0.2 Hz, hence, it is directly compatible with the KBR1B data. Since the LTC for microwave ranging needs to be referred to the antenna phase center (APC), the position of the phase center in the satellite frame, as provided by VKB1B,⁴ is rotated using the star camera SCA1B product into the GCRF. The COM-APC offset in the GCRF is added onto the rotated GNV1B data in order to obtain the position and velocity of the APC on each SC in the GCRF. The acceleration vector of the APC is approximated by the center-of-mass acceleration from force models, which is justified, since the angular motion of the APC on time scales of the light propagation time is negligible. The APC state vectors are used to derive the one-way LTCs ${\mathcal {T}}_\mathrm{AeBr}$ and ${\mathcal {T}}_{BeAr}$ (Eq. (34)), which are further combined using Eq. (49) into ${\mathcal {T}}_\text {DOWR}$ with K- and Ka-band frequencies from the USO1B dataset.

The difference between the light time correction from GRACE level-1b KBR data (GRA_KBR1B_LTC) and $c_0 \cdot {\mathcal {T}}_\text {DOWR}$ (Eq. (49)) with four different degrees of accuracy is shown in Fig. 3. The data used spans the GPS time between 00:00 and 06:00 on December 1st, 2008. Since the differences are minimal when only the special relativistic correction ${\mathcal {T}}_\text {SR}$ is used (red trace), it is reasonable to assume that general relativistic contributions were omitted in the GRACE level-1b light time correction. The omission error is dominated by the sinusoidal variation at the orbital frequency, however, with an amplitude of approx. 1 micrometer, i.e. close to the tone error requirement discussed in Sect. 8 for GRACE Follow-On.

The GRACE level-1b LTC shows some artifacts above 10 mHz (magenta trace on the right subplot in Fig. 3). However, these are well below the KBR noise level and should not impede the gravity field recovery.

Fig. 3

Comparison between GRACE level-1b light time correction and ${\mathcal {T}}_\text {DOWR}$ (Eq. (49)) using different degrees of accuracy in the time (left) and spectral (right) domain. The traces on the left plot have been centered around zero by subtracting a bias shown in the legend. The difference is minimal when only the special relativistic effect is considered in ${\mathcal {T}}_\text {DOWR}$. The dominating amplitudes and the mean values are provided Table 3

For GRACE Follow-On, an additional orbit data product called GNI1B is available, which provides the satellite state in the GCRF and can be used instead of the transformed GNV1B data. The sampling rate is 1 Hz, which means that results need to be downsampled to the KBR and LRI rates of 0.2 and 0.5 Hz, respectively. A comparison with different degrees of accuracy for the light time correction is shown in Fig. 4 for February 5th, 2019. It is evident that the LTC in GRACE Follow-On takes into account the general relativistic effect ${\mathcal {T}}_\text {PM}$ due to the central field (degree 0), but not the higher moments ${\mathcal {T}}_\text {HM}$. The omission error is present predominantly at twice the orbital frequency with a peak amplitude of approx. 0.1 $\upmu \hbox {m}$ (blue trace), thus well below the discussed requirement from Sect. 8. The differences between $c_0 \cdot {\mathcal {T}}_\text {TWR}$ and the RL04 LTC in Fig. 4 are limited to a level of a few $\hbox {nm/}\sqrt{\text {Hz}}$, which is well below the LRI noise requirement.

However, the actual LRI in-orbit noise is close to 1 $\hbox {nm/}\sqrt{\text {Hz}}$ at Fourier frequencies around 0.1 Hz, hence, we study the limits of the LTC precision and propose potential improvements for the LTC in the next section.

Fig. 4

Comparison between GRACE Follow-On level-1b KBR/LRI light time correction and ${\mathcal {T}}_\text {DOWR}$ (left) and ${\mathcal {T}}_\text {TWR}$ (right) using different degrees of accuracy. The dominating amplitudes and the mean values for the different traces are provided Table 4

11 Enhancing the light time correction accuracy

In order to understand the current limit of the LTC precision of a few $\hbox {nm}/\sqrt{\text {Hz}}$, we reproduced the light time corrections provided in the GRACE Follow-On RL04 data. In a first step (step 1), the classical light time equation was solved iteratively to obtain the absolute light travel time $\varDelta t$, and, in a second step (step 2), the instantaneous contribution ($\varDelta t_\text {inst} = |\mathbf {r}_A - \mathbf {r}_B|/c_0$) was removed from $\varDelta t$ in order to obtain the one-way corrections ${\mathcal {T}}$, which are further combined into ${\mathcal {T}}_\text {DOWR}$ or ${\mathcal {T}}_\text {TWR}$.

We noted a slight inconsistency in the instantaneous Euclidean inter-satellite distance between GNI1B or GNV1B products, which shows rms differences three times higher compared to our method to rotate the GNV1B data into the GCRF (cf. left panel in Fig. 5). The precision limit of our method is the resolution of the double floating-point arithmetic, i.e. the computation error of the product of rotation matrix and position vector.

We could reproduce the light time correction of RL04 data with smallest deviations, if we used different orbit sets in step 1 and for the calculation of $\varDelta t_\text {inst}$ in step 2 (cf. green dashed trace on right subplot of Fig. 5). However, using consistent orbit sets for both steps results in a slightly lower noise for the light time correction (solid blue trace). The consistent data sets could be GNI for both steps (denoted as orbit data OD2 in the plot), or the rotated GNV data (denoted as orbit data OD3). A difference between both cases is not apparent in the spectrum, hence, the plot shows a single solid blue trace for both cases. The dashed black trace on the right plot of Fig. 5 depicts the actual in-orbit measurements of the LRI Abich et al. (2019), which contains the instrument noise but also some variations due to non-gravitational accelerations (nga) for the shown frequencies Misfeldt (2019).

Fig. 5

(Left plot:) Difference in inter-satellite distance between satellite C and D for different orbit data products. GNI and GNV are the RL04 datasets, while rotGNV denotes a dataset, which has been rotated by the authors from ITRF to GCRF. (Right plot:) Spectral density of the LTC signal (red and blue traces) and LTC differences for different orbit data sets. The LTCs have been computed in four different cases (OD1..OD4), which are based in different orbit data sets for step 1 (iter: solving for $\varDelta t$ iteratively) and step 2 (calculating $\varDelta t_\text {inst}$ and computing ${\mathcal {T}} = \varDelta t - \varDelta t_\text {inst}$). This plot was created with a log-scale amplitude spectral density (LASD) method, which produces smooth traces also at high frequencies Tröbs and Heinzel (2006)

Fig. 6

Amplitude spectral density (ASD) of the differences of the LTC with different inter-satellite range rate data for one day in August 1 2019 with logarithmic scaled frequency axis. (Left:) comparison for the LTC range of KBR, (Right:) comparison for the LTC range of LRI. In addition, the cyan blue trace in the right subplot shows the LTC from a kinematic orbit of GRACE (December 1, 2008)

The LTC accuracy can be improved further - well below the sensitivity of the LRI - by using the analytical expressions for ${\mathcal {T}}$ as discussed in Sect. 6 and 7, where the dominating terms in the single-path are proportional to $\mathbf {d}_0 .\dot{\mathbf {r}}_{A/B}$, or in the final DOWR and TWR combination ${\mathcal {T}}_\text {DOWR/TWR} \propto \mathbf {d}_0 .(\dot{\mathbf {r}}_{A} - \dot{\mathbf {r}}_{B}) \propto {\dot{\rho }}_\text {inst,OD}$ (cf. Eqs. (52) and (71)). If the satellite velocity vectors $\dot{\mathbf {r}}_{A/B}$ are derived as the time-derivative of the satellite position state vector, the accuracy of the LTC is limited to the nm/$\sqrt{\text {Hz}}$ level. However, if the velocity state vectors of GNI1B are used, the LTC noise is highly reduced as shown by the blue traces in left and right plot of Fig. 6. This results from the fact that GNI and GNV data is based on reduced-dynamic orbit determination, where the variational equations include the velocity state Bertiger et al. (2010, 2020).

It is noteworthy that the instantaneous range rate ${\dot{\rho }}_\text {inst,OD}$, which appears in the first-order approximation of the LTC (Eqs. (52) and (71)), dominates the noise in the LTC. Fortunately, the instantaneous range rate is approximately the same as the more precise measured range rate from LRI or KBR with only a minor light time correction from an orbit product, i.e.

$$\begin{aligned}&{\dot{\rho }}_\text {inst,TWR} \approx {\dot{\rho }}_\text {TWR}- \frac{\text {d}}{\text {d} t} \frac{|\mathbf {r}_A - \mathbf {r}_B| \cdot {\dot{\rho }}_\text {inst,OD}}{c_0},\nonumber \\&\quad {\dot{\rho }}_\text {inst,DOWR} \approx {\dot{\rho }}_\text {DOWR}- \frac{\text {d}}{\text {d} t} \frac{|\mathbf {r}_A - \mathbf {r}_B| \cdot {\dot{\rho }}_\text {inst,OD}}{2 c_0}. \end{aligned}$$

(77)

Thus, if ${\dot{\rho }}_\text {inst,OD}$ from the orbit product is replaced with ${\dot{\rho }}_\text {inst,TWR}$ or ${\dot{\rho }}_\text {inst,DOWR}$ in the dominating term of the LTC, the resulting LTC becomes almost independent of the orbit product. The result exhibits very low noise at high-frequencies (red trace on the right plot in Fig. 6) that is comparable to the pure GNI LTC (dashed blue trace). The deviations below 2 mHz are caused by differences between ranging and orbit data, and it is reasonable to assume that the results using Eq. (77) are more accurate than the LTC based purely on orbit data.

Moreover, the above replacement allows us to use even kinematic orbit products for the LTC calculation with acceptable high frequency noise (cyan blue trace in Fig. 6 for GRACE data Zehentner and Mayer-Gürr (2013)). For that trace, the high frequency noise above 25 mHz is driven by the KBR ranging noise. Kinematic orbits are sometimes regarded as more appropriate for gravity field recovery Naeimi and Flury (2017), since they do not rely on a-priori gravity field information.

Finally, we note that the most accurate way to determine the instantaneous biased range $\rho _\text {inst}$ is to update the LTC in the process of combined orbit determination and gravity field recovery with the most current orbit estimate in each iteration. In other words, one can consider to use the non-instantaneous biased range as observation and shift the conversion by means of the LTC into the process of precise orbit determination and gravity field recovery, where the LTC is updated iteratively.

12 Summary and conclusions

The Laser Ranging Interferometer aboard GRACE Follow-On demonstrated for the first time laser ranging between satellites in a gravimetric satellite mission. This enables inter-satellite biased range observations with an unprecedented noise level of 1 $\hbox {nm/}\sqrt{\text {Hz}}$ at the highest frequency in the level-1b data (0.25 Hz), or even 0.2 $\hbox {nm/}\sqrt{\text {Hz}}$ at the highest frequency of the level-1a data (5  Hz).

The biased range observation needs to be corrected for the effect of the finite speed of light in order to obtain the instantaneous range between the spacecraft, which is the quantity utilized in the gravity field recovery process. It is natural to seek methods to compute the light time corrections with a higher precision in order to not limit the observations of the GRACE Follow-On LRI, and potentially also of future instruments and missions.

In this paper, we revisited the calculation of the light time correction from first principles within the Post-Newtonian approximation of general relativity, taking into account state-of-the-art geopotential models. We have separated the total light time correction ${\mathcal {T}}$ into the contribution from special relativity ${\mathcal {T}}_\text {SR}$ and the general relativistic component into the effect from the scalar central field of the Earth (${\mathcal {T}}_\text {PM}$, SH degree 0), from higher moments of the gravity potential, which includes direct tidal accelerations, ${\mathcal {T}}_\text {HM}$, and from the much smaller vector potential due to Earth’s spin moment ${\mathcal {T}}_\text {SM}$. The analytical formulas were verified against the light travel time obtained by numerically integrating the equations of motion of photons.

We studied in Sect. 9 the influence of different geopotential models onto ${\mathcal {T}}_\text {HM}$, showing that to reach tone-errors below 1 pm amplitude in the LTC, one should consider the effect from the Sun and the Moon, as well as from Solid Earth tides. In order to achieve a noise level in the light time correction below 100 pm/$\sqrt{\text {Hz}}$, the SH degree of the static gravity field should be above 50 and the light path between satellites needs to be sampled with more than 10 points.

We showed that the GRACE light time correction in RL02 does not consider general relativistic effects, while GRACE Follow-On RL04 data takes into account general relativity with a radial-symmetric field (${\mathcal {T}}_\text {PM}$). The omission of ${\mathcal {T}}_\text {HM}$ causes predominantly a sinusoidal error with a peak amplitude well below 1 $\upmu \hbox {m}$ at twice the orbital frequency. The LTC in the official RL04 data is limited to a noise level of a few $\hbox {nm/}\sqrt{\text {Hz}}$ arising from numerical floating point precision and due to the fact that two slightly inconsistent orbit products (GNI and GNV) are used in each step. This level of LTC precision is comparable to the LRI instrument noise at the highest frequencies in the level-1b data.

The numerical accuracy can be easily improved to 1 $\hbox {nm/}\sqrt{\text {Hz}}$ at high frequencies by using the same orbit product in both steps. However, we recommend to use the here proposed analytical formulas as these are numerically a few orders of magnitude more accurate, as the absolute LTC accuracy depends on the models and on the orbit product quality.

In the end it was pointed out that, if the analytical formulas are employed, the dominating term of ${\mathcal {T}}_\text {TWR}$ or ${\mathcal {T}}_\text {DOWR}$ can be rewritten in terms of the measured range rate from LRI or KBR, which means the LTC becomes to first-order independent of the orbit product. Hence, kinematic orbit products that suffer higher noise can be used to compute the LTC as well.

The here presented methods to calculate the light time correction for microwave and laser ranging can be readily applied to simulated and available flight data. The analytical approximations were truncated at picometer level, which is well below the requirements for the current GRACE Follow-On mission, but may be of interest in studies for future missions.