© The Authors 2023

1 Introduction

Electromagnetic signals propagating through the atmosphere of Earth are used for time and frequency transfer between distant atomic clocks located on the surface of the Earth or at satellites with applications in fundamental science (Delva et al. 2017, 2018; Safronova et al. 2018; Beloy et al. 2021), geodesy (Müller et al. 2018; Mehlstäubler et al. 2018; Delva et al. 2019), and metrology (Fujieda et al. 2014; Hachisu et al. 2014; Riedel et al. 2020; Pizzocaro et al. 2021). Future perspectives in this field are described, for example, by Alonso et al. (2022) and Derevianko et al. (2022).

Recent developments in optical atomic clocks that reach fractional systematic uncertainty in frequency on the order of 10⁻¹⁸ in case of stationary clocks (Ludlow et al. 2015; McGrew et al. 2018; Lisdat et al. 2021) or 10⁻¹⁷ in case of transportable clocks (Koller et al. 2017; Cao et al. 2017; Origlia et al. 2018), together with developments in free-space optical links (see Bodine et al. 2020 and references therein, and also Djerroud et al. 2010, Gozzard et al. 2018, Kang et al. 2019, Dix-Matthews et al. 2021 and Shen et al. 2021) enable increasing accuracies of the experiments. Therefore, the atmospheric effects originating from spatial and temporal variations of the refractive index of air or from winds start to play a more important role.

In case of the ground-to-satellite time and frequency transfer, two projects are going to be realized in the near future that will place an atomic clock on the International Space Station (ISS), namely the Atomic Clock Ensemble in Space (ACES), with a Cs clock and H-maser (Cacciapuoti et al. 2020; Lilley et al. 2021), and the Space Optical Clock on the ISS (I-SOC), with an optical lattice clock (Cacciapuoti & Schiller 2017; Origlia et al. 2018). In the ACES mission, the absolute frequency accuracy of the on-board clock is about 10⁻¹⁶ in fractional frequency, whereas the same parameter for the I-SOC mission is about 10⁻¹⁷. This gives an uncertainty limit to one of the main objectives of the experiments, that is, the test of the gravitational redshift. To compare time and frequency between various clocks, ACES will use a microwave link (MWL) and an optical link of the European Laser Timing (ELT) experiment (Schreiber et al. 2010). When used for time comparisons between ground clocks, the MWL is expected to provide absolute synchronization of the ground clock timescales with an uncertainty of 100ps. For the ELT, the overall planned accuracy of the time transfer is 50 ps (Turyshev et al. 2016). The I-SOC proposal uses an improved version of the MWL as well as the optical link of ELT, aiming at an accuracy of a few ps in the clock synchronization (Cacciapuoti & Schiller 2017). Therefore, based on the planned accuracies of these experiments, a relativistic model of one- and two-way time and frequency transfer is needed that includes the atmospheric effects and has an accuracy of approximately 1 ps for the time transfer and a relative accuracy of the lower multiples of 10⁻¹⁸ for the frequency transfer.

Blanchet et al. (2001) developed the relativistic theory of one- and two-way time and frequency transfer in vacuum, including terms up to the order of c⁻³ to cover contributions relevant for experiments with an uncertainty of 5 × 10⁻¹⁷ in frequency transfer and 5 ps in time transfer. This theory is based on a solution of the null geodesic equation for a light ray in vacuum, static, spherically symmetric spacetime.

Relativistic theory of the propagation of electromagnetic signals in the broader context of astrometric measurements in the Solar System is discussed, for example, in Kopeikin et al. (2011) and references therein. The propagation time and frequency shift of a signal in the gravitational field of multiple moving bodies were studied especially in Kopeikin & Schäfer (1999) and Kopeikin & Mashhoon (2002), where the Liénard-Wiechert representation of the metric tensor up to the first post-Minkowskian order was used.

Linet & Teyssandier (2002) applied the Synge world function formalism to the gravitational field of an axisymmetric rotating body, and they showed that certain c⁻⁴ terms in frequency shift approach the 10⁻¹⁸ value for frequency transfer from a ground station to the ISS. They also discussed the influence of the quadrupole moment of the Earth J₂ at the c⁻³ order. The theory of Linet & Teyssandier (2002) was further developed by Le Poncin-Lafitte et al. (2004) and Teyssandier & Le Poncin-Lafitte (2008), where the formalism of time transfer functions was introduced. Hees et al. (2014a,b) employed this formalism to compute various observables relevant for actual space missions in the Solar System, including the coordinate propagation time and the frequency shift of light.

The relativistic theoretical works mentioned above considered the vacuum case alone and did not take the atmospheric effects into account. The theory of Teyssandier & Le Poncin-Lafitte (2008) was extended to include the atmospheric effects by Bourgoin (2020). Here, a general formalism of time transfer functions was developed using the Gordon optical metric (Gordon 1923), and it was applied to a case of stationary optical spacetime. The effects of refractivity and motion of a neutral atmosphere were added to the gravitational effects, and the light-dragging effect due to the steady rotation of the atmosphere of Earth was discussed. The models of Linet & Teyssandier (2002) and Bourgoin (2020) focused on one-way time and frequency transfer only. The results of Bourgoin (2020) were later used to model the atmospheric occultation experiments (Bourgoin et al. 2021). In this work, explicit formulas for time and frequency transfer were derived for a steadily rotating spherically symmetric atmosphere up to the first post-Minkowskian order, with the remark that the method described in the paper easily allows extending the results to higher orders and beyond the spherical symmetry. Feng et al. (2022) used the formalism of Gordon’s optical metric in the context of global navigation satellite systems to include the atmospheric effects in the relativistic framework. The null geodesics of the Gordon metric were solved numerically. The effect of motion of the atmosphere of Earth on the speed of light (the Fresnel-Fizeau effect) and the related propagation time delay were discussed also by Kopeikin & Han (2015) in the context of geodetic very long baseline interferometry.

Several effects originating from the atmospheric refractivity have been addressed in the literature in a nonrelativistic framework. The influence of the refractive index fluctuations given by atmospheric turbulence has been discussed by Sinclair et al. (2014, 2016), Robert et al. (2016), Belmonte et al. (2017), Swann et al. (2019), and Taylor et al. (2020). Nonrelativistic analytical solutions of light rays in planetary atmospheres were discussed by Bourgoin et al. (2019). Propagation-time nonreciprocity in two-way ground-to-satellite time transfer due to the distribution of the atmospheric refractive index was discussed by Stuhl (2021). Earlier findings on tropospheric and ionospheric corrections were summarized, for example, in Petit & Luzum (2010). A two-way time transfer model including the effect of ionospheric dispersion was discussed by Duchayne et al. (2009), who also provided uncertainty requirements on the orbit determination of space clocks.

The aim of this paper is to develop a relativistic model of one- and two-way time and frequency transfer that includes gravitational and atmospheric effects for realistic fields of the gravitational potential, the refractive index, and the wind speed, assuming neither stationarity nor spherical or axial symmetry. The intended accuracy of the model is given by the experimental needs described above: approximately 1 ps for time transfer, and lower multiples of 10⁻¹⁸ for frequency transfer. The model is based on the solution of the equation of motion for a light ray, that is, the null geodesic equation in Gordon’s optical metric, which is given by the fields of the gravitational potential, the atmospheric refractive index, and the wind speed in the vicinity of the Earth. Dispersion of the medium and free electric charges or currents are not included in the model.

The paper is organized as follows. After defining the notation and conventions in Sect. 2, we introduce the coordinates, the spacetime metric, and the fields describing the medium in Sect. 3. Then, Sect. 4 summarizes the main results of the paper for easy reference. In the Sect. 4, the level of approximation of the model is explained and the solution for the light rays in the spherically symmetric, static case is derived because this solution is needed to express the results for the general case without the symmetries, the formulas for one- and two-way time and frequency transfer are given, and several examples are shown that quantify the atmospheric effects for some typical situations whose parameters are similar to the experiments at the ISS. A complete derivation of the results is then given in Sect. 5. In the appendix, some basic facts about the refractive index of air and its distribution in the atmosphere of Earth are summarized.

2 Notation and conventions

In this paper, c is the speed of light in a vacuum, and G is the Newtonian gravitational constant. The spacetime metric signature (− + + +) is used. Small Latin indices run from 1 to 3, and small Greek indices run from 0 to 3. Partial derivatives with respect to spacetime coordinates are denoted as ∂_α = ∂/∂x^α and also ∂_t = ∂/∂t. Einstein’s summation convention is used. The Euclidean norm of a triplet of components aⁱ is denoted as $\left|a^i\right|\equiv \sqrt{{\delta }_{ij}{a}^i{a}^j}$, with δ_ij being the Kronecker delta. Quantities ωⁱ, ${\upsilon }_R^i$, υⁱ, Vⁱ, Aⁱ, hⁱ ,χⁱ, ςⁱ, and γⁱ in this paper are treated as Euclidean three-vectors when their index is lowered or raised, for example, V_i = δ_ijV^j The three-dimensional antisymmetric Levi-Civita symbol is denoted as ϵ_ijk, eventually ${ϵ}^i{}_{jk}={\delta }^{im}{ϵ}_{mjk}$.

3 Coordinates, fields, and observers

3.1 Background metric and coordinate system

The Geocentric Celestial Reference System (GCRS) is a coordinate system in the surroundings of the Earth that is centered in the center of mass of the Earth and is nonrotating with respect to distant stars. We denote the GCRS coordinates as ${\tilde{x}}^{\alpha }=\left({\tilde{x}}^0,{\tilde{x}}^i\right)$, where ${\tilde{x}}^0/c=t$ is the Geocentric Coordinate Time (TCG), and ${\tilde{x}}^i$ are the spatial coordinates of the system. The components of the spacetime metric in the GCRS coordinates up to the order that is sufficient for the analysis further are (Soffel et al. 2003) $${\tilde{g}}_{00}=-1+\frac{2W}{c^2}+O\left(c^{-4}\right)$$(1a) $${\tilde{g}}_{0i}=O\left(c^{-3}\right)$$(1b) $${\tilde{g}}_{ij}={\delta }_{ij}\left(1+\frac{2W}{c^2}\right)+O\left(c^{-4}\right),$$(1c)

with W being the scalar gravitational potential defined by Soffel et al. (2003) with the convention W ≥ 0.

To evaluate the time and frequency transfer, however, we use coordinates that rotate with respect to the GCRS system, with a constant angular velocity vector approximating the real angular velocity vector of the rotation of the Earth around its axis, such that the surface of the Earth is nearly at rest in the new coordinates, and we can use the advantage of nearly vanishing velocities of stationary observers on the ground. The use of this corotating system leads to more complicated formulas for oneway time and frequency transfer, but, on the other hand, the fact that the position of the stationary observer on the ground is almost constant makes the calculation of two-way corrections more straightforward.

The rotating coordinate system x^α = (x⁰, xⁱ) is related to the GCRS system as $x^0={\tilde{x}}^0$, $x^i=R^i{}_j{\tilde{x}}^j$, where Rⁱ_j(t) ∊ SO(3), or inversely, ${\tilde{x}}^i=R_j^i{x}^j$, where $R_j^i={\delta }_{jk}{\delta }^{il}{R}^k{}_l$. The time-dependent rotation matrix Rⁱ_j(t) corresponds to the rigid uniform rotation with a constant angular velocity vector that approximates the real angular velocity vector of the Earth. Therefore, we have $$\frac{\text{d}{R}_n^i}{\text{d}t}={ϵ}_{ijk}{\tilde{\omega }}^j{R}_n^k,$$(2)

where ${\tilde{\omega }}^j$ are constant components of the angular velocity vector of the coordinate rotation expressed in the ${\tilde{x}}^i$ frame. Since ${\tilde{\omega }}^j$ is defined to be constant and the coordinate rotation is rigid, the irregularities in the rotation of the Earth and the deformations of the Earth given mainly by the solid Earth tides are observed as tiny movements of the surface of the Earth in the coordinates xⁱ. We denote ${\tilde{\upsilon }}_R^i=\text{d}{\tilde{x}}^i/\text{d}t={ϵ}_{ijk}{\tilde{\omega }}^j{\tilde{x}}^k$ the velocity of a point that is fixed in the rotating frame with respect to the GCRS frame. In terms of components in the rotating frame, we have ${\upsilon }_R^i=R^i{}_j{\tilde{\upsilon }}_R^j={ϵ}_{ijk}{\omega }^j{x}^k$ with ${\omega }^i=R^i{}_j{\tilde{\omega }}^j$.

The components of the metric (1) in the rotating coordinates are $$g_{00}=-1+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{c^2}+O\left(c^{-4}\right)$$(3a) $$g_{0i}=\frac{1}{c}{\upsilon }_R^i+O\left(c^{-3}\right)$$(3b) $$g_{ij}={\delta }_{ij}\left(1+\frac{2W}{c^2}\right)+O\left(c^{-4}\right),$$(3c)

where ${\upsilon }_R^2={\delta }_{ij}{\upsilon }_R^i{\upsilon }_R^j$. The monopole part of the scalar gravitational potential is denoted as $\widehat{W}\left(r\right)$, and the remaining part containing the higher multipole moments and the time-dependent components as ∆W(x^α) such that $$\begin{array}{cc}W\left(x^{\alpha }\right)=\widehat{W}\left(r\right)+\text{Δ}W\left(x^{\alpha }\right)& \text{With}\end{array}$$(4) $$\widehat{W}\left(r\right)=\frac{GM}{r},$$(5)

where r = ∣xⁱ∣ is the radial coordinate, and M is the mass of the gravitating object. For details of the gravitational potential W(x^α) in the vicinity of the Earth, see, for example, Wolf & Petit (1995) and Voigt et al. (2016).

3.2 Medium

In our model, the atmosphere of Earth is represented by a nondis-persive medium, and we assume that its optical properties are given by a refractive index field n(x^α). This simplification can be applied, for example, in the case of a monochromatic signal by setting the frequency (or vacuum wavelength) variable in the refractive index of real air to a fixed value, assuming that the slight frequency variations along the signal trajectory, as observed from the local rest frames of the medium, have a negligible dispersion effect. We also assume that the medium is electrically neutral and has a vanishing density of the electric current.

The refractive index of air as a function of temperature, pressure, composition, and vacuum wavelength is discussed, for example, by Owens (1967) and Ciddor (1996; see also Appendix A.1). When the fields of the atmospheric temperature, the pressure, and the composition are known, the refractive index field n(x^α) for a given vacuum wavelength can be determined.

The atmospheric temperature and pressure as functions of altitude show a systematic behavior at large scales that is discussed, for example, in ICAO (1993) and NOAA (1976; see also Appendix A.2). Spatial and temporal local fluctuations are present as well. Therefore, it is reasonable to decompose the atmospheric refractive index field as $$n\left(x^{\alpha }\right)=\widehat{n}\left(r\right)\left(1+\alpha \left(x^{\alpha }\right)\right),$$(6)

where $\widehat{n}\left(r\right)$ is the part of the refractive index field that depends on the radial coordinate alone (i.e., it is static and spherically symmetric), and α(x^α) is its relative correction, which can be time dependent and in general has no spatial symmetry.

In our model, we consider general fields $\widehat{n}\left(r\right)$, α(x^α). When the model is to be applied to a specific case, these fields must be specified. A detailed discussion of the possible definitions of these fields is beyond scope of this paper, but some guidelines for determining $\widehat{n}\left(r\right)$ for altitudes below 80 km are given in the appendix. The correction α(x^α) then consists of various more or less predictable components, such as the deviation from spherical symmetry caused by the centrifugal potential, the effects of the changing position of Sun as a heat source, the effects of large-scale meteorological structures, or components coming from variations in local temperature, pressure, and composition, and from atmospheric turbulence.

The total refractivity N(x^α) and the static, spherically symmetric refractivity $\widehat{N}\left(r\right)$ are then defined as deviations of the corresponding refractive index from its vacuum value, $$N=n-1,$$(7) $$\widehat{N}=\widehat{n}-1.$$(8)

Denoting $x_M^i\left(t,x_0^{\alpha }\right)$ the trajectory of the medium element in the rotating coordinates that passes through a spacetime point $x_0^{\alpha }=\left(c{t}_0,x_0^i\right)$, the winds in the atmosphere are described by their velocity field Vⁱ(x^α), which is defined by $$V^i\left(x_0^{\alpha }\right)=\frac{\partial x_M^i}{\partial t}\left(t,x_0^{\alpha }\right)|{}_{t=t_0}$$(9)

for any $x_0^{\alpha }$ where the medium is present. We also define the following quantity that is useful when the optical effects of the winds are evaluated: $$A^i=\left(1-n^2\right)V^i.$$(10)

3.3 Observers

For a given trajectory xⁱ (t) of an observer in the rotating coordinates, we denote υⁱ = dxⁱ/dt the observer’s velocity in these coordinates.

4 Time and frequency transfer through the atmosphere of Earth. Summary of the results

This section summarizes the main results of the paper: the corrections for one- and two-way time and frequency transfer through the atmosphere of Earth. It also provides numerical examples that give an idea about the magnitude of the atmospheric effects. The details of the derivation of the results are given in Sect. 5.

4.1 Approximation level of the model

As mentioned in the introduction, a model accuracy of approximately 1 ps for the time transfer and relative accuracy of lower multiples of 10⁻¹⁸ for the frequency transfer is needed in order to meet the requirements of the forthcoming clock-on-satellite experiments. The corresponding approximation level of the model should be defined.

The influencing quantities entering the equation of motion and the evaluation of the time and frequency transfer in the rotating frame originate from gravitation given by the Newtonian potential W/c² (higher-order gravitational effects such as the Lense-Thirring effect are negligible considering our intended accuracy), from inertial effects given by the rotation velocity ${\upsilon }_R^i/c$, from atmospheric effects given by the air refractivity $\widehat{N}$ together with the correction α and by the wind speed Vⁱ /c and from the observers’ velocities υⁱ/c at the emission and reception side.

We define the “order” of a term in an expansion according to the following procedure: (a) we express all formulas in terms of the refractivity $\widehat{N}\left(r\right)$ and its correction α(x^α), (b) we denote ${\widehat{N}}_M$ the maximum value of $\widehat{N}$ and α_M the maximum value of ∣α∣, (c) we define normalized functions ${\widehat{N}}_{*}=\widehat{N}/{\widehat{N}}_M$ and α_* = α/α_M and in all formulas we express $\widehat{N}={\widehat{N}}_M{\widehat{N}}_{*}$ and α = α_Mα_* (d) the order of a term in an expansion is then given by its factor ${\left(c^{-1}\right)}^{p_{c^{-1}}}{\left({\widehat{N}}_M\right)}^{p_{\widehat{N}}}{\left({\alpha }_M\right)}^{p_{\alpha }}$, where the powers $p_{c^{-1}}$, $p_{\widehat{N}}$, and p_α are nonnegative integers, and finally, (e) the derivative of the fields ∂₀ = c⁻¹ ∂_t increases the power of c⁻¹ by one. The approximation level of an expansion is then given by a selection of terms defined by a subset in the grid of possible triplets of the exponents $\left(p_{c^{-1}},p_{\widehat{N}},p_{\alpha }\right)$. The expansion of n in terms of ${\widehat{N}}_M$, ${\widehat{N}}_{*}$ α_M, and α_* serves for the purpose of defining the order, but is not shown explicitly in formulas where it is not practical.

First, considering terms on the order of ${\left(c^{-1}\right)}^{p_{c^{-1}}}{\left({\widehat{N}}_M\right)}^{p_{\widehat{N}}}$ (i.e., p_α = 0) in the time and frequency transfer formulas, we take $p_{c^{-1}}\le 3$ corresponding to the vacuum model of Blanchet et al. (2001), and we take an arbitrary value of $p_{\widehat{N}}$ for $p_{c^{-1}}\le 2$. Terms on the order of $c^{-3}{\left({\widehat{N}}_M\right)}^{p_{\widehat{N}}}$ with $p_{\widehat{N}}\ge 1$ are neglected. Some of the terms on the order $c^{-3}{\widehat{N}}_M$ can slightly exceed the 10⁻¹⁸ level in case of ground-to-ISS frequency transfer with a near horizon position of the satellite. The c⁻⁴ terms, which we neglect as well, can approach the 10⁻¹⁸ level for ground-to-ISS frequency transfer, as was shown by Linet & Teyssandier (2002). Then, taking the field α into account, its effect is considered to be minor compared to the spherical, static part of the refractivity $\widehat{N}$. We therefore only include the field α to a lower order. Namely, for p_α ≥ 1, we take $p_{c^{-1}}+p_{\widehat{N}}+p_{\alpha }\le 3$. Possible magnitudes of the neglected terms with p_α ≥ 1 have not been evaluated, however.

Next, we denote the total $p_{c^{-1}}+p_{\widehat{N}}+p_{\alpha }=p_T$, and we can summarize the approximation level of the model as follows: the expressions for the signal propagation time, the relative frequency shift, and the related two-way time and frequency transfer corrections in this paper include all terms with the order from a set 𝒫(3), which is given as $$𝒫\left(3\right)=\left\{\left(p_{c^{-1}},p_{\widehat{N}},p_{\alpha }\right):p_T\le 3\vee \left(p_{c^{-1}}\le 2,p_{\alpha }=0\right)\right\}.$$

The only exception in which we express the results up to a lower order only are the terms in the two-way corrections that originate from the motion of a stationary observer on the ground in the corotating frame during the back and forth propagation of the signal. This motion is given by the deviations of the surface motion of Earth from rigid uniform rotation and therefore is very small. The approximation level of these terms is discussed together with the particular formulas.

Next we introduce the notation O(q) for expressions containing terms with p_T ≥ q only. Expressions containing terms with $p_{\widehat{N}}=p_{\alpha }=0$ and $p_{c^{-1}}\ge k$ only are denoted O(c^−k) as usual. It is also useful to define the following subset of the exponent triplets: $$𝒫\left(2\right)=\left\{\left(p_{c^{-1}},p_{\widehat{N}},p_{\alpha }\right):p_T\le 2\vee \left(p_{c^{-1}}\le 1,p_{\alpha }=0\right)\right\}.$$

4.2 Strategy of solving the problem

The approach used in this paper is based on the formulation and on an approximate solution of the equation of motion for a light ray in a flowing medium and a gravitational field. This equation is given as the null geodesic equation of the corresponding Gordon optical metric (Gordon 1923). Details of the solution are given in Sect. 5.

In summary, our strategy of solving the equation of motion is the following. We formulate the equation in the Newtonian form in the corotating coordinate system with an Euclidean length of the light ray (in the same system) as a parameter, and we split the equation into its static, spherically symmetric part plus a correction. First, an exact solution of the static, spherically symmetric part is found with the given spatial positions of the emission and reception events. We denote this solution as ${\widehat{x}}^i\left(\widehat{l}\right)$, with $\widehat{l}$ being a parameter given by the Euclidean length of this path from its initial point. Then, an approximate linearized equation is solved for the path correction $\text{Δ}{x}^i\left(\widehat{l}\right)$ caused by the deviations of the refractive index and the gravitational potential from sphericity and staticity, by winds and by the Coriolis and centrifugal forces. Finally, a separate first-order differential equation for the coordinate time $t\left(\widehat{l}\right)$ is solved.

The main results of this paper are expressed in terms of the integrals along the path ${\widehat{x}}^i\left(\widehat{l}\right)$. We therefore discuss the solution of the static, spherically symmetric part of the equation of motion in detail in this summary.

4.3 Light rays in a static, spherically symmetric atmosphere

The equation of motion for a light ray in a static, spherically symmetric distribution of the refractive index $\widehat{n}\left(r\right)$ and gravitational potential Ŵ(r) is given as (see Eq. (132)) $$\frac{{\text{d}}^2{\widehat{x}}^i}{\text{d}{\widehat{l}}^2}=\left(\frac{1}{\widehat{n}}{\partial }_j\widehat{n}+\frac{2}{c^2}{\partial }_j\widehat{W}\right)\left({\delta }^{ij}-\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\right),$$(11)

where the fields $\widehat{n}$, ${\partial }_j\widehat{n}$, and ∂jŴ are evaluated in ${\widehat{x}}^k\left(\widehat{l}\right)$, and the solution has a unit tangent satisfying $${\delta }_{{}^{ij}}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}=1,$$(12)

corresponding to the fact that is its Euclidean length parameter. We assume that the initial and final points of the path are known as input parameters, and we later relate them to spatial coordinates of the emission and reception events of the one- or two-way time and frequency transfer setups. We denote the total path length between the initial and final points as $\widehat{L}$, such that the parameter has the range $\widehat{l}\in \left[0,\widehat{L}\right]$.

Defining a field w as $$w=\widehat{n}\exp \left(2\widehat{W}/c^2\right),$$(13)

which depends on the spatial coordinates through the radius r = ∣xⁱ∣ only and introducing a new parameter ζ with a function dependence $\zeta \left(\widehat{l}\right)$ satisfying $$\frac{\text{d}\zeta }{\text{d}\widehat{l}}=\frac{1}{w\left({\widehat{x}}^i\left(\widehat{l}\right)\right)},$$(14)

Eq. (11) can be written as $$\frac{{\text{d}}^2{\widehat{x}}^i\left(\zeta \right)}{\text{d}{\zeta }^2}=\frac{1}{2}{\partial }_i{w}^2\left({\widehat{x}}^j\left(\zeta \right)\right),$$(15)

where we briefly denote ${\widehat{x}}^i\left(\zeta \right)\equiv {\widehat{x}}^i\left(\widehat{l}\left(\zeta \right)\right)$. Thus, we obtain an equation of motion of the signal in the central potential field. The method of solving an analogous equation in mechanics is known from the literature (see, e.g., Landau & Lifshitz 1976). It is based on two conservation laws. One law corresponds to the total energy and can be obtained from Eqs. (12) and (14) in our case, and it is given as $${\delta }_{ij}\frac{\text{d}{\widehat{x}}^i\left(\zeta \right)}{\text{d}\zeta }\frac{\text{d}{\widehat{x}}^j\left(\zeta \right)}{\text{d}\zeta }=w^2\left({\widehat{x}}^k\left(\zeta \right)\right).$$(16)

The second law corresponds to the angular momentum and can be obtained by taking a cross product of Eq. (15) with the signal position vector ${\widehat{x}}^i\left(\zeta \right)$ and using the assumption of the central field, which leads to $${ϵ}_{jk}^i{\widehat{x}}^j\left(\zeta \right)\frac{\text{d}{\widehat{x}}^k\left(\zeta \right)}{\text{d}\zeta }=h^i,$$(17)

with hⁱ being a constant vector. Equation (17) shows that for ∣hⁱ∣ ≠ 0, the signal propagates in a plane that intersects the origin r = 0 and has a normal aligned with hⁱ (for ∣hⁱ∣ = 0, the signal path is just a straight radial line). Denoting $\psi \left(\widehat{l}\right)$ the angle between the tangent $\text{d}{\widehat{x}}^i\left(\widehat{l}\right)/\text{d}\widehat{l}$ and the position vector $\widehat{x}\left(\widehat{l}\right)$ and taking the magnitude of Eq. (17) with use of Eq. (16), we obtain $$r\left(\widehat{l}\right)w\left({\widehat{x}}^i\left(\widehat{l}\right)\right)\sin \psi \left(\widehat{l}\right)=h=\text{const.,}$$(18)

where $r\left(\widehat{l}\right)=\left|{\widehat{x}}^i\left(\widehat{l}\right)\right|$ and h = ∣hⁱ∣. This is the analog of the Snell law in static, spherically symmetric media with a gravitational field of the same symmetry. For the nongravitational case (Ŵ = 0), Eq. (18) gives the Bouguer formula of classical geometrical optics (see Eq. (7) in Sect. 3.2 of Born & Wolf 1999). Understanding r as the trajectory parameter, with (r) being the inverse function to $r\left(\widehat{l}\right)$ (assuming $\text{d}r/\text{d}\widehat{l}\ne 0$ all along the trajectory), Eq. (18) can be written as $$rw\left(r\right)\sin \psi \left(\widehat{l}\left(r\right)\right)=h.$$(19)

Selecting a right-handed Cartesian coordinate system yⁱ with the origin at r = 0, the y¹-axis in ${\widehat{x}}^i\left(0\right)$ direction, and the y³-axis in the hⁱ direction (for h = 0, we can take any direction perpendicular to ${\widehat{x}}^i\left(0\right)$), defining polar coordinates by y¹ = r cos φ and y² = r sin φ and expressing the trajectory ${\widehat{x}}^i\left(\widehat{l}\right)$ in the polar coordinates as φ(r), with the radius as the parameter, it follows from the definition of ψ that $$r\frac{\text{d}\phi }{\text{d}r}=\text{tan\hspace{0.17em}}\psi \left(\widehat{l}\left(r\right)\right).$$(20)

Integrating this equation with tan ψ expressed from the Snell law (19), we obtain $$\phi \left(r\right)=\pm {\displaystyle {\int }_{r_I}^r\text{d}{r}^{\prime }}\frac{h}{r^{\prime }}{\left({r^{\prime }}^{2}w{\left(r^{\prime }\right)}^2-h^2\right)}^{-\frac{1}{2}},$$(21)

with $\pm =+\text{\hspace{0.17em}for\hspace{0.17em}}\text{d}r/\text{d}\widehat{l}>0$ (ascending trajectory) and ± = − for $\text{d}r/\text{d}\widehat{l}<0$ (descending trajectory).

Thus, up to the final integration, we obtain the solution of Eq. (11) in polar coordinates parameterized by r. The remaining step is to express the constant h in terms of the input parameters, which, in this case, are the boundary values of the trajectory ${\widehat{x}}^i\left(\widehat{l}\right)$.

In the following, we denote quantities related to the initial point ${\widehat{x}}^i\left(0\right)$ by the index I and to the final point ${\widehat{x}}^i\left(\widehat{L}\right)$ by the index F, for example, ${\widehat{x}}_I^i={\widehat{x}}^i\left(0\right)$, $r_I=\left|{\widehat{x}}_I^i\right|$, $w_I=w\left({\widehat{x}}_I^i\right)$, ψ_I = ψ(0), and ${\widehat{x}}_F^i={\widehat{x}}^i\left(\widehat{L}\right)$, $r_F=\left|{\widehat{x}}_F^i\right|$, $w_F=w\left({\widehat{x}}_F^i\right)$, and ${\psi }_F=\psi \left(\widehat{L}\right)$. We also define the angle $\overline{\phi }$ between ${\widehat{x}}_F^i$ and ${\widehat{x}}_I^i$, the angle θ_I between ${\widehat{x}}_F^i-{\widehat{x}}_I^i$ and ${\widehat{x}}_I^i$, and the angle θ_F between ${\widehat{x}}_F^i-{\widehat{x}}_I^i$ and ${\widehat{x}}_F^i$, and we denote $D=\left|{\widehat{x}}_F^i-{\widehat{x}}_I^i\right|$ (see Fig. 1). The following relations hold: $$r_F\sin {\theta }_F=r_I\sin {\theta }_I,$$(22) $$r_F\cos {\theta }_F=r_I\cos {\theta }_I+D.$$(23)

The boundary value problem is more complicated than the initial value problem, in which case h is given simply as h = r_Iw_I sin ψ_I. In the boundary value problem, we need to express h as a function of r_I, r_F and one of the angles $\overline{\phi }$, θ_I, or θ_F. One possible approach is presented below. We denote as ε_I and ε_F the angular deviations of the path tangents ${\text{d}{\widehat{x}}^i/\text{d}\widehat{l}|}_0$ and ${\text{d}{\widehat{x}}^i/\text{d}\widehat{l}|}_{\widehat{L}}$, respectively, from the vector ${\widehat{x}}_F^i-{\widehat{x}}_I^i$ (see Fig. 1), so that we have $${\epsilon }_I={\psi }_I-{\theta }_I,$$(24) $${\epsilon }_F={\psi }_F-{\theta }_F.$$(25)

Evaluating the Snell law (19) in the boundary points and using Eqs. (24) and (25), we obtain $$h=r_I{w}_I\sin \left({\theta }_I+{\epsilon }_I\right)=r_F{w}_F\sin \left({\theta }_F+{\epsilon }_F\right).$$(26)

Now we need to express the angle ε_I or ε_F. For this purpose, we introduce a total bending angle of the path ε_T that is the angle between ${\text{d}{\widehat{x}}^i/\text{d}\widehat{l}|}_0$ and ${\text{d}{\widehat{x}}^i/\text{d}\widehat{l}|}_{\widehat{L}}$, which satisfies $${\epsilon }_F-{\epsilon }_I={\epsilon }_T.$$(27)

Solving Eq. (26) together with Eq. (27) for the angles ε_I and ε_F, we obtain $$\tan {\epsilon }_I=\frac{\frac{w_F}{w_I}\sin \left({\theta }_F+{\epsilon }_T\right)-\sin {\theta }_F}{-D/r_F-\frac{w_F}{w_I}\text{\hspace{0.17em}cos}\left({\theta }_F+{\epsilon }_T\right)+\cos {\theta }_F},$$(28) $$\tan {\epsilon }_F=\frac{\frac{w_I}{w_F}\sin \left({\theta }_I-{\epsilon }_T\right)-\sin {\theta }_I}{D/r_I-\frac{w_I}{w_F}\text{\hspace{0.17em}cos}\left({\theta }_I-{\epsilon }_T\right)+\cos {\theta }_I}.$$(29)

The total bending angle also satisfies $${\epsilon }_T=\overline{\phi }+{\psi }_F-{\psi }_I,$$(30)

where the angle $\overline{\phi }$ can be expressed as $\overline{\phi }=\phi \left(r_F\right)$ using Eq. (21), and we can also write $${\psi }_F-{\psi }_I={\displaystyle {\int }_{r_I}^{r_F}\text{d}r\frac{\text{d}\psi \left(\widehat{l}\left(r\right)\right)}{\text{d}r}},$$(31)

where on the right-hand side (RHS), we use Eq. (19). Thus, we obtain $${\epsilon }_T\left(h\right)=\mp {\displaystyle {\int }_{r_I}^{r_F}\text{d}r}\frac{h}{w\left(r\right)}\frac{\text{d}w\left(r\right)}{\text{d}r}{\left(r^{2}w{\left(r\right)}^2-h^2\right)}^{-\frac{1}{2}}$$(32)

with ∓ = − for the ascending path and ∓ = + for the descending path.

Then, the constant h can be determined by the following iterative procedure: (a) we estimate h = r_Fw_F sin θ_F, corresponding to ε_F = 0 in Eq. (26), (b) we calculate ε_T(h) according to Eq. (32), (c) we calculate ε_F according to Eq. (29), (d) we calculate h according to Eq. (26) and start another loop from (b). As a result, not only the constant h is obtained, but also the angles ε_F, ε_T, and ε_I = ε_F − ε_T. Alternatively, the procedure can be performed using the angle ε_I given by Eq. (28) instead of ε_F. We also note that for the ground-to-satellite or satellite-to-ground transfer, the ε_F/I angle at the satellite position is at least one order of magnitude smaller than the corresponding angle at the ground position. Therefore, neglecting the ε_F/I angle at the satellite position gives a better initial estimate of h according to Eq. (26) than neglecting the corresponding angle at the ground position.

To retrieve the solution in the coordinates xⁱ, we express the components of the first two basis vectors of yⁱ in the coordinate system xⁱ. Assuming $\overline{\phi }\ne 0$, we have $$\begin{array}{cc}{e}_1^i=\frac{{\widehat{x}}_I^i}{r_I},& e_2^i=\frac{1}{\sin \overline{\phi }}\left(\frac{{\widehat{x}}_F^i}{r_F}-\text{cos\hspace{0.17em}}\overline{\phi }\frac{{\widehat{x}}_I^i}{r_I}\right)\end{array}.$$(33)

The solution parameterized by r is then given as $${\widehat{x}}^i\left(\widehat{l}\left(r\right)\right)=r\text{\hspace{0.17em}cos\hspace{0.17em}}\phi \left(r\right)e_1^i+r\sin \phi \left(r\right)e_2^i.$$(34)

For $\overline{\phi }=0$, it is simply ${\widehat{x}}^i\left(\widehat{l}\left(r\right)\right)=r{e}_1^i$.

The reparameterization function $\widehat{l}\left(r\right)$ can be expressed with the use of $$\frac{\text{d}\widehat{l}}{\text{d}r}=\pm \sqrt{1+r^2{\left(\frac{\text{d}\phi }{\text{d}r}\right)}^2}=\frac{1}{\cos \psi \left(\widehat{l}\left(r\right)\right)},$$(35)

where the ± sign has the same meaning as in Eq. (21). Integrating this formula using the Snell law (19), we obtain $$\widehat{l}\left(r\right)=\pm {\displaystyle {\int }_{r_I}^r\text{d}{r}^{\prime }{\left(1-\frac{h^2}{{r^{\prime }}^{2}w{\left(r^{\prime }\right)}^2}\right)}^{-\frac{1}{2}}}.$$(36)

For $\widehat{L}$ we then get $\widehat{L}=\widehat{l}\left(r_F\right)$.

It is also useful to define an orthonormal right-handed basis (χⁱ, ςⁱ, γⁱ) with χⁱ pointing in direction from ${\widehat{x}}_I^i$ to ${\widehat{x}}_F^i$, γⁱ pointing in direction of the normal to the propagation plane hⁱ (assuming h ≠ 0), and ςⁱ completing the triplet. We have (see also Fig. 1) $${\chi }^i=\frac{1}{D}\left({\widehat{x}}_F^i-{\widehat{x}}_I^i\right),$$(37) $${\varsigma }^i=\frac{1}{D}\left({\widehat{x}}_F^i\text{\hspace{0.17em}cot\hspace{0.17em}}{\theta }_I-{\widehat{x}}_I^i\text{\hspace{0.17em}cot\hspace{0.17em}}{\theta }_F\right),$$(38) $${\gamma }^i=\frac{{ϵ}^i{}_{jk}{\widehat{x}}_I^j{\widehat{x}}_F^k}{r_I{r}_F\sin \overline{\phi }}.$$(39)

Defining the angle $\epsilon \left(\widehat{l}\right)$ between the path tangent $\text{d}{\widehat{x}}^i/\text{d}\widehat{l}$ and the vector χⁱ satisfying ε(0) = ε_I (negative value) and $\epsilon \left(\widehat{L}\right)={\epsilon }_F$ (positive value), we can write $$\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}={\chi }^i\cos \epsilon \left(\widehat{l}\right)+{\varsigma }^i\sin \epsilon \left(\widehat{l}\right).$$(40)

In some formulas, the lowest-order approximation of the tangent $\text{d}{\widehat{x}}^i/\text{d}\widehat{l}$ and of the length $\widehat{L}$ is sufficient. Expanding Eq. (40) in ε and using ∣ε∣ ≤ ε_T = O(1), which follows from Eq. (32), we obtain $$\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}={\chi }^i+O\left(1\right).$$(41)

Next we obtain $$D={\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\chi }_i\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}={\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\cos \epsilon \left(\widehat{l}\right)=\widehat{L}+O\left(2\right),$$(42)

where the expansion of cos ε was used in the last step.

In the following text, we refer to the solution ${\widehat{x}}^i\left(\widehat{l}\right)$ of Eq. (11) with the given boundary conditions as the base path.

4.4 Coordinate time transfer

4.4.1 One-way time transfer

In one-way time transfer, we wish to express a coordinate time t_F − t_I of the propagation of the signal emitted from a position $x_I^i$ at a time t_I and received in a position $x_F^i$ at a time t_F. We consider a base path ${\widehat{x}}^i\left(\widehat{l}\right)$ with boundary points ${\widehat{x}}^i\left(0\right)=x_I^i$ and ${\widehat{x}}^i\left(\widehat{L}\right)=x_F^i$.

Based on the analysis given in Sect. 5, we obtain the following formula for the propagation time, including all terms on the order of 𝒫(3): $$\begin{array}{l}{t}_F-t_I=\frac{1}{c}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}n+\frac{1}{c^3}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}2W\hfill \\ +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left({\upsilon }_{Ri}+A_i\right)\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\hfill \\ +\frac{1}{c^3}\frac{D}{2}{\upsilon }_{RI}^i{\upsilon }_{RF}^j\left({\delta }_{ij}+{\chi }_i{\chi }_j\right)\hfill \\ +\frac{1}{c^2}{ϵ}^i{}_{jk}{\chi }^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_i\alpha }\hfill \\ +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}{\partial }_t\alpha }\hfill \\ +\frac{1}{2c}\left({\delta }^{ij}-{\chi }^i{\chi }^j\right){\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}{\widehat{l}}^{\prime }}D\left(\widehat{l},{\widehat{l}}^{\prime }\right){\partial }_i\alpha \left(t_I,{\widehat{x}}^k\left(\widehat{l}\right)\right){\partial }_j\alpha \left(t_I,{\widehat{x}}^k\left({\widehat{l}}^{\prime }\right)\right),\hfill \end{array}$$(43)

where the fields n, α, W, υ_Ri, A_i and their partial derivatives with respect to the space-time coordinates are evaluated in t = t_I and $x^i={\widehat{x}}^i\left(\widehat{l}\right)$ if not stated explicitly otherwise, the velocities ${\upsilon }_{RI}^i$ and ${\upsilon }_{RF}^i$ are defined as ${\upsilon }_R^i\left({\widehat{x}}_I^k\right)$ and ${\upsilon }_F^i\left({\widehat{x}}_F^k\right)$, respectively, and the function $D\left(\widehat{l},{\widehat{l}}^{\prime }\right)$ is given by Eq. (149). Equation (43) is expressed in corotating coordinates, which means that all components of the fields are expressed as functions of the corotating coordinates, the wind speed Vⁱ in Aⁱ is the speed in the corotating frame, and the base path is defined by the boundary conditions that are given by the spatial coordinates of the emission and reception events in the corotating frame as well.

Equation (43) corresponds to Eq. (5) (or (A.39)) of Blanchet et al. (2001), transformed to the rotating frame and generalized by including the atmospheric effects. Its first term is the leading Newtonian term in the refractive medium, the second term is the Shapiro delay, the first term in the second line is the Sagnac correction, which appears due to the use of the rotating frame as well as the term in the third line, and the second term in the second line is the effect of the wind, corresponding to the Fresnel-Fizeau effect of dragging of light by a medium. The remaining terms containing α are the effects of nonsphericity and nonstationarity of the refractive index field. Numerical integration of the terms in Eq. (43) can be performed, for example, by transforming the integration variable $\widehat{l}$ into the radial coordinate r using Eq. (36) and by expressing the integrated functions in spherical coordinates with the base path given by Eq. (21).

4.4.2 Two-way time transfer

In two-way transfer, we consider the signal emitted from a stationary observer on the ground A at the coordinate time t_A1 , which corresponds to the proper time τ_A(t_A1) of clock A. The spatial coordinates of this event are denoted $x_{A1}^i$. The signal is then received by observer B (ground or satellite) at the coordinate time t_B, which corresponds to the proper time τ_B(t_B) of clock B (see Fig. 2). The spatial coordinates of this event are denoted $x_B^i$. The pseudo-time-of-flight τ_B(t_B) − τ_A(t_A1) is obtained from the measurements. Then a signal is sent from observer B at the coordinate time t_B and is received by observer A at the coordinate time t_A2, which corresponds to the proper time τ_A(t_A2) of clock A. The spatial coordinates of this event are denoted $x_{A1}^i$. This signal can either be the same signal reflected or another signal that is synchronously sent when the one-way signal is received by observer B. We call this set-up the Λ-configuration. The pseudo-time-of-flight τ_A(t_A2) − τ_B(t_B) is also obtained from measurements.

Using the coordinate-time synchronization convention, the desynchronization of the clocks is defined as the difference τ_B(t_B) − τ_A(t_B), where τ_A(t_B) is the proper time measured by observer A at the coordinate time t_B. In case of the two-way time transfer, it can be expressed as $$\begin{array}{l}{\tau }_B\left(t_B\right)-{\tau }_A\left(t_B\right)=\frac{1}{2}\left({\tau }_B\left(t_B\right)-{\tau }_A\left(t_{A1}\right)\right)\hfill \\ -\frac{1}{2}\left({\tau }_A\left(t_{A2}\right)-{\tau }_B\left(t_B\right)\right)\hfill \\ +\frac{1}{2}\left(\text{Δ}{\tau }_{-}-\text{Δ}{\tau }_{+}\right),\hfill \\ \hfill \end{array}$$(44)

where ∆τ₊ = τ_A(t_B) − τ_A(t_A1) and ∆τ₋ = τ_A(t_A2) − τ_A(t_B). The first two lines of Eq. (44) are the measured pseudo-times-of-flight, and the difference ∆τ₋ − ∆τ₊ needs to be computed.

Denoting ∆t₊ = t_B − t_A1 the coordinate time of the signal propagation from A to B and ∆t₋ = t_A2 − t_B the coordinate time of the signal propagation from B to A, the computed term can be approximated as $$\text{Δ}{\tau }_{-}-\text{Δ}{\tau }_{+}\approx \text{Δ}{t}_{-}-\text{Δ}{t}_{+}.$$(45)

The error of this approximation is several orders below the 1 ps level in terrestrial conditions.

We also note that in the corotating coordinates we use, the position of the stationary observer on the ground A is nearly fixed and changes only due to the deformations of Earth, which are given, for example, by the solid Earth tides or by the irregularity of the rotation of Earth. Thus, the difference between $x_{A1}^i$ and $x_{A2}^i$ is very small, and we include its effect up to the largest order only using the velocity ${\upsilon }_{A1}^i$ of observer A in the corotating system at the emission event.

The propagation times ∆t₊ and ∆t₋ can be calculated using Eq. (43) or its slightly extended version, Eq. (183). In their difference, the leading term in the first line of Eq. (43) and several other terms are compensated for. Denoting ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ the base path with the boundary points ${\widehat{x}}_{+}^i\left(0\right)=x_{A1}^i$ and ${\widehat{x}}_{+}^i\left(\widehat{L}\right)=x_B^i$, we arrive at the following result for the two-way time transfer correction, including all terms on the order of 𝒫(3) except the p_T = 3 order terms originating from $x_{A2}^i-x_{A1}^i$ $$\begin{array}{l}\text{Δ}{t}_{-}-\text{Δ}{t}_{+}=-\frac{2}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left({\upsilon }_{Ri+}+A_{i+}\right)\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}\hfill \\ +\frac{2}{c^2}{ϵ}^i{}_{jk}{\chi }_{+}^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\left(\widehat{L}-\widehat{l}\right){\partial }_i{\alpha }_{+}\hfill \\ +\frac{2}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left(\widehat{L}-\widehat{l}\right){\partial }_t{\alpha }_{+}\hfill \\ -\frac{2}{c^2}\widehat{L}{\chi }_{i+}{\upsilon }_{A1}^i,\hfill \end{array}$$(46)

where the fields υ_Ri, A_i, ∂_iα, and ∂_tα with the index + are evaluated in t = t_B and $x^i={\widehat{x}}_{+}^i\left(\widehat{l}\right)$, and ${\chi }_{+}^i$ is the vector χⁱ corresponding to the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ (i.e., it is the unit vector in the direction from $x_{A1}^i$ to $x_B^i$).

The main contribution to the two-way time transfer correction in Eq. (46) is given by the first term in the first line, which is the well-known Sagnac effect. We denote $${\left(\text{Δ}{t}_{-}-\text{Δ}{t}_{+}\right)}_{\text{Sag}}\text{=\hspace{0.17em}\hspace{0.17em}}-\frac{2}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\upsilon }_{Ri+}\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}.$$(47)

Then, taking into account that ${\upsilon }_{Ri+}={ϵ}_{i,jk}{\omega }^j{\widehat{x}}_{+}^k\left(\widehat{l}\right)$, the correction given by Eq. (47) can be expressed as $${\left(\text{Δ}{t}_{-}-\text{Δ}{t}_{+}\right)}_{\text{Sag}}\text{=\hspace{0.17em}\hspace{0.17em}}-\frac{4}{c^2}{\omega }^j{\text{Σ}}_j,$$(48)

where Σ_j is defined as $${\text{Σ}}_j=\frac{1}{2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{ϵ}_{jki}{\widehat{x}}_{+}^k\left(\widehat{l}\right)\frac{\text{d}{\widehat{x}}_{+}^i\left(\widehat{l}\right)}{\text{d}\widehat{l}}.$$(49)

The magnitude Σ = ∣Σ_j∣ equals the area of the plane section surrounded by three lines: 1) the straight line from the origin of the coordinate system (the Earth center) to $x_{A1}^i$, 2) the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ from $x_{A1}^i$ to $x_B^i$, and 3) the straight line from $x_B^i$ to the origin of the coordinate system (see Fig. 3). The direction of Σ_j is normal to this plane. The difference between the Sagnac effect in the Earth atmosphere and in vacuum is then given by the change in the area Σ when the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ changes from the case including atmospheric refraction and gravitation to a case that only includes gravitation.

A convenient formula for Σ_j can be obtained when the integrand of Eq. (49) is expressed using Eq. (17) together with Eq. (14) and using the definition γ_i = h_i/h. Then, transforming the integration variable to r using Eq. (36), we obtain $${\text{Σ}}_j=\frac{1}{2}{\gamma }_j{}_{+}h{\displaystyle {\int }_{r_A}^{r_B}\text{d}r{\left(w{\left(r\right)}^2-\frac{h^2}{r^2}\right)}^{-\frac{1}{2}},}$$(50)

where $r_A=\left|x_{A1}^i\right|,r_B=\left|x_B^i\right|,{\gamma }_{j+}$ is the vector γ_j corresponding to the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$, and w(r) is given by Eq. (13). This integral can be computed numerically for any refractive index profile $\widehat{n}\left(r\right)$ of interest.

The second term in the first line of Eq. (46) is the effect of the wind, and the following two lines give the effects of nonsphericity and nonstationarity of the refractive index field. The last line of Eq. (46) then describes the effect of motion of the stationary observer on the ground in the corotating frame.

Possible magnitudes of the atmospheric contribution to the Sagnac effect and of the effect of wind are discussed in the next section. A numerical evaluation of the other effects is left for future research.

4.5 Time transfer. Examples and magnitudes of the effects

4.5.1 Sagnac correction

In this section, we evaluate the Sagnac correction, Eq. (48), for the example of ground-to-satellite transfer through an isothermal atmosphere at a temperature T₀ with an unchanging composition corresponding to a molar mass of air M_a. We assume that the atmosphere rotates rigidly with the Earth and is in hydrostatic equilibrium with its gravitational and centrifugal forces. In this case, the refractivity $\widehat{N}\left(r\right)=\widehat{n}\left(r\right)-1$ is given by Eq. (A.13), which gives the following spherical refractive index field: $$\widehat{n}\left(r\right)=1+N_A\exp \left(\frac{M_a}{R{T}_0}\left(\widehat{W}-{\widehat{W}}_A\right)\right),$$(51)

where Ŵ = GM/r, Ŵ_A = GM/r_A, R is the universal gas constant, and N_A is the air refractivity in the reference position $x_{A1}^i$, which also equals $\widehat{N}\left(r_A\right)$ (for details, see the appendix). For the sake of simplicity, we use Eq. (51) in the full range of the radius between the ground station and the satellite. The function w(r) in Eq. (50) is then calculated according to Eq. (13).

We consider r_A = 6371 km, corresponding to the surface of Earth, r_B = r_A + 408 km, approximately corresponding to the altitude of the ISS, and T₀ = 288.15 K (15 °C) and M_a = 28.964 × 10⁻³ kg mol⁻¹, corresponding to dry air with a carbon dioxide molar fraction x_c = 450 ppm. For an atmospheric pressure of 101 325 Pa in the position $x_{A1}^i$ and a light signal with a vacuum wavelength of 1 µm, we obtain (see, e.g., Ciddor 1996) $$N_A\approx 2.742\times {10}^{-4}.$$(52)

In two-way transfer, we denote the θ_I angle of the ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ path simply as θ. It represents the zenith angle of the satellite at the reception/re-emission event (see Fig. 3).

Based on Eqs. (48) and (50), the Sagnac correction can be calculated as a function of θ. In Fig. 4, this function is shown (the dashed blue curve) together with the effect of the atmosphere (the red full curve), which is defined as the Sagnac correction in the atmosphere (N_A given by Eq. (52)) minus the Sagnac correction in vacuum (N_A = 0) for the same value of θ. The values in the plots in Fig. 4 are obtained for the case ωⁱγ_i+ = −ω, which corresponds to the paths ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ propagating in the equatorial plane against the rotation of the Earth. This is the case with the highest positive values of the Sagnac correction.

Figure 4 shows that the Sagnac correction increases from 0 to approximately 24 ns, with θ increasing from 0 to 90°. The effect of the atmosphere is below 0.1 ps for θ < 78°, it reaches 1 ps for θ ≈ 86.6°, and it increases to approximately 5 ps as θ approaches 90°. Therefore, for the time transfer at an accuracy level of 1 ps, this effect is only significant for large angles θ.

Since the gravitational effect to ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ is much smaller than the effect of atmospheric refraction, the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ can be approximated by a straight line for angles θ where the effect of the atmosphere is negligible. Thus, for example, for θ ∊ [0°, 78°], the following formula for the Sagnac correction has a bias below 0.1 ps for the example studied in this section, $${\left(\text{Δ}{t}_{-}-\text{Δ}{t}_{+}\right)}_{\text{Sag}}=-\frac{2}{c^2}{\omega }^i{\gamma }_i+r_A{r}_B\sin \overline{\phi }.$$(53)

4.5.2 Effect of the wind

We evaluate the effect of the wind speed in the two-way time transfer correction (46), which is given by $${\left(\text{Δ}{t}_{-}-\text{Δ}{t}_{+}\right)}_{\text{wind}}=-\frac{2}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{A}_{i+}}\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}.$$(54)

For a rough estimate of the magnitude of the effect, we assume that the signal propagates in the part of the atmosphere in which the refractive index n is approximately constant along the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$, as well as the component of the air velocity Vⁱ (in the corotating frame) in the direction of the path tangent $\text{d}{\widehat{x}}_{+}^i/\text{d}\widehat{l}$, which we denote V. The approximately constant n corresponds, for example, to the horizontal signal path from the ground to the ground or to the initial part of the signal path from the ground to the satellite when the satellite is in a position near the horizon. In this case, we obtain $${\left(\text{Δ}{t}_{-}-\text{Δ}{t}_{+}\right)}_{\text{wind}}\approx \frac{2}{c^2}\left(n^2-1\right)V\widehat{L}.$$(55)

For a 1 ps correction in the air with the refractivity value of Eq. (52), which is typical on the surface of the Earth, we need $V\widehat{L}\approx 8.2\times {10}^7$ m² s⁻¹, which, for example, for the path length of $\widehat{L}=100$ km leads to V ≈ 820 ms⁻¹. Thus, the effect is negligible for time transfer at an accuracy level of 1 ps under normal atmospheric conditions. However, it can be detected by inter-ferometric methods, as was shown, for example, by Tselikov & Blake (1998), who used the effect in a different context for flow metering applications in laboratory conditions.

4.6 Frequency transfer

4.6.1 One-way frequency transfer

In one-way frequency transfer, the goal is to express the ratio of the proper frequency ν¡ of a signal emitted by an observer with velocity ${\upsilon }_I^i$ from a position $x_I^i$ at a time t_I and a proper frequency ν_F of the same signal as is received by the observer with a velocity ${\upsilon }_F^i$ in the position $x_F^i$ at the time t_F. Again, we consider a base path ${\widehat{x}}^i\left(\widehat{l}\right)$ with the boundary points ${\widehat{x}}^i\left(0\right)=x_I^i$ and ${\widehat{x}}^i\left(\widehat{L}\right)=x_F^i$, and we denote $W_I=W\left(t_I,x_I^i\right)$, $W_F=W\left(t_F,x_F^i\right)$, ${\upsilon }_{RI}^i={\upsilon }_R^i\left(x_I^j\right)$, and ${\upsilon }_{RF}^i={\upsilon }_R^i\left(x_F^j\right)$

Based on the analysis included in Sect. 5, we obtain the following formula for the frequency ratio: $$\frac{v_I}{v_F}=\frac{1+\frac{W_I}{c^2}+\frac{1}{2c^2}{\left|{\upsilon }_{RI}^i+{\upsilon }_I^i\right|}^2}{1+\frac{W_F}{c^2}+\frac{1}{2c^2}{\left|{\upsilon }_{RF}^i+{\upsilon }_F^i\right|}^2}\frac{1-\frac{1}{c}{\upsilon }_I^i{\left({\overline{l}}_i\right)}_I}{1-\frac{1}{c}{\upsilon }_F^i{\left({\overline{l}}_i\right)}_F}\exp I,$$(56)

where ${\left({\overline{l}}_i\right)}_I$ and ${\left({\overline{l}}_i\right)}_F$ are quantities related to tangents of the light-ray trajectory at the emission and reception events, and they are given by Eqs. (196)–(197d) with η = 0 for ${\left({\overline{l}}_i\right)}_I$ and η = 1 for ${\left({\overline{l}}_i\right)}_F$. In these formulas, we set $\Delta x_F^i=0$ and ${\widehat{l}}_0=0$ corresponding to ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)=\left(c{t}_I,{\widehat{x}}^i\left(\widehat{l}\right)\right)$. The details of the expressions for ${\left({\overline{l}}_i\right)}_I$ and ${\left({\overline{l}}_i\right)}_F$ are given in the Theory section, but the corresponding contribution to two-way transfer is discussed further in this summary. The exponent I in the last term of Eq. (56) is given as $$\begin{array}{l}I={\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left(\frac{1}{c}{\partial }_{t}n+\frac{1}{c^2}\left({\partial }_t{A}_i\right)\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}+\frac{2}{c^3}{\partial }_{t}W\right)\hfill \\ +\frac{1}{c^2}{ϵ}^i{}_{jk}{\chi }^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right)}{\partial }_i{\partial }_t\alpha \hfill \\ +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}{\partial }_t{\partial }_t\alpha }\hfill \\ +\frac{1}{c}\left({\delta }^{ij}-{\chi }^i{\chi }^j\right){\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\int }_0^{\widehat{L}}\text{d}{\widehat{l}}^{\prime }D\left(\widehat{l},{\widehat{l}}^{\prime }\right){\partial }_i{\partial }_t\alpha \left(t_I,{\widehat{x}}^k\left(\widehat{l}\right)\right){\partial }_j\alpha \left(t_I,{\widehat{x}}^k\left({\widehat{l}}^{\prime }\right)\right)}\hfill \end{array}$$(57)

where the fields n, α, A_i, W, and their partial derivatives with respect to the space-time coordinates are evaluated in t = t_I and $x^i={\widehat{x}}^i\left(\widehat{l}\right)$ if not stated explicitly otherwise.

Equation (56) together with Eqs. (196)–(197d) and Eq. (57) gives the frequency ratio including all terms on the order of 𝒫(3). The first fraction in Eq. (56) contains the gravitational redshift effect and the second-order Doppler effect, the second fraction contains the first-order Doppler effect, and the exponential term is the effect of nonstationarity.

4.6.2 Two-way frequency transfer

In two-way frequency transfer, a stationary observer on the ground A emits a signal with the proper frequency ν_A1 from the position $x_{A1}^i$ at the time t_A1. The signal is received by observer B (ground or satellite) with the proper frequency ν_B in the position $x_B^i$ at the time t_B and is immediately transponded back with the same frequency to observer A, where it is received with the proper frequency ν_A2 in the position $x_{A2}^i$ at the time t_A2. The velocities of observers A and B at the particular emission and reception events are denoted as ${\upsilon }_{A1}^i$, ${\upsilon }_B^i$, and ${\upsilon }_{A2}^i$. For the stationary (fixed at the ground) observer A, the velocities ${\upsilon }_{A1}^i$ and ${\upsilon }_{A2}^i$ in the corotating coordinates and the corresponding spatial shift $x_{A2}^i-x_{A1}^i$ are given just by the deformations in the body of the Earth, for example, the solid tides, or by irregularities in the rotation of Earth. Thus, these quantities are very small, and the main contribution to the Doppler shift comes from the satellite part and not from the ground.

The frequency ratio ν_A2/ν_B, which is needed for the frequency comparison, is then expressed in terms of the ratio ν_A2 /ν_A1, which is measured by the reference clock at ground station A, and the correction, which must be computed. Ashby (1998) and Blanchet et al. (2001) defined the correction Δ by the following relation: $$\frac{v_{A2}}{v_B}=\frac{1}{2}\frac{v_{A2}}{v_{A1}}+\text{Δ}+\frac{1}{2}.$$(58)

The correction Δ, including all terms on the order of 𝒫(3) except the p_T ≥ 3 order terms originating from the motion of observer A in the corotating frame and except the terms quadratic in ${\upsilon }_{A1}^i/c$, is given by Eqs. (212)–(215) in the Theory section. In this summary, we present a simpler result for Δ that is obtained by neglecting terms of higher than the third order (p_T > 3) and terms proportional to c⁻²N_B, which are negligible for satellite applications because the atmospheric refractivity N_B in the satellite position is very low for satellite altitudes of several hundred kilometers.

The fields evaluated in the spacetime points $x_{A1}^{\alpha }$, $x_B^{\alpha }$, or $x_{A2}^{\alpha }$ of the emission and reception events are denoted by the corresponding index, for example, $W_{A1}=W\left(t_{A1},x_{A1}^j\right)$ and ${\upsilon }_{R|A1}^i={\upsilon }_R^i\left(x_{A1}^j\right)$). The base path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ is again defined by the boundary points ${\widehat{x}}_{+}^i\left(0\right)=x_{A1}^i$ and ${\widehat{x}}_{+}^i\left(\widehat{L}\right)=x_B^i$, and the fields denoted with the index + are evaluated in t = t_B and $x^i={\widehat{x}}_{+}^i\left(\widehat{l}\right)$, for example, ${\partial }_b{\alpha }_{+}=\left({\partial }_b\alpha \right)\left(t_B,{\widehat{x}}_{+}^i\left(\widehat{l}\right)\right)$. The basis; ${\chi }_{+}^i,$, ${\varsigma }_{+}^i,$, ${\gamma }_{+}^i$ is associated with the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$, and ε_B is the angle ε_F of the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$. The acceleration of observer A at the time t_A1 in the rotating frame is denoted $a_{A1}^i$ (i.e., $a_{A1}^i={\text{d}}^2{x}_A^i/{\text{d}{t}^2|}_{t_{A1}}$). Using this notation, we obtain $$\begin{array}{l}\text{Δ}=\frac{1}{2c^2}\left[W_{A1}+W_{A2}-2W_B+{\left|{\upsilon }_{R|A1}^i\right|}^2+2{\upsilon }_{Bi}{\upsilon }_{R|A1}^i-{\left|{\upsilon }_{R|B}^i+{\upsilon }_B^i\right|}^2\right]\hfill \\ \times \left(1-\frac{1}{c}{\upsilon }_B^j{\chi }_{j+}\right)\hfill \\ +\frac{1}{c^2}D{\chi }_{i+}\left(a_{A1}^i+{ϵ}^i{}_{jk}{\omega }^j{\upsilon }_{A1}^k\right)+\frac{1}{c^2}{\upsilon }_{Bi}{\upsilon }_{A1}^i\hfill \\ +\frac{1}{c}{\upsilon }_B^i{\beta }_{i\left(2\right)}+I_{+\left(c^{-2}\right)}-\frac{1}{c^2}{\upsilon }_B^i{\chi }_{i+}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\partial }_t{n}_{+},}\hfill \end{array}$$(59)

where the atmospheric effects are included in the last line, in which $$\begin{array}{l}{\beta }_{i\left(2\right)}=\frac{1}{c}{\chi }_{+}^{i}D{\omega }_j{\gamma }_{+}^j{\epsilon }_B\hfill \\ +\frac{1}{c}\left({\delta }^{ij}-{\chi }_{+}^i{\chi }_{+}^j\right){ϵ}^b{}_{jk}{\omega }^k\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_b{\widehat{n}}_{+}}\hfill \\ +\frac{1}{c}\left({\delta }^{ia}-{\chi }_{+}^i{\chi }_{+}^a\right){ϵ}^b{}_{jk}{\chi }_{+}^j{\omega }^k\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\widehat{l}}^2\left(\widehat{l}-\widehat{L}\right){\partial }_a{\partial }_b{\widehat{n}}_{+}}\hfill \\ +\frac{1}{c}{ϵ}_{ijk}{\chi }_{+}^j\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\text{\hspace{0.17em}curl}{A}_{+}^k}\hfill \\ +\frac{1}{c}\left({\delta }^{ij}+2{\chi }_{+}^i{\chi }_{+}^j\right){ϵ}^b{}_{jk}{\omega }^k\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_b{\alpha }_{+}}\hfill \\ +\frac{1}{c}{ϵ}^b{}_{jk}{\chi }_{+}^j{\omega }^k\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\widehat{l}}^2\left(\widehat{l}-\widehat{L}\right){\partial }_i{\partial }_b{\alpha }_{+}}\hfill \\ +\frac{1}{c}\left({\delta }^{ij}-{\chi }_{+}^i{\chi }_{+}^j\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_t{\partial }_j{\alpha }_{+}}\hfill \end{array}$$(60)

and $$\begin{array}{l}{I}_{+\left(c^{-2}\right)}=\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left({\partial }_t{A}_{i+}\right)}\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}\hfill \\ +\frac{1}{c^2}{ϵ}^i{}_{jk}{\chi }_{+}^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_i{\partial }_t{\alpha }_{+}}\hfill \\ +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_t{\partial }_t{\alpha }_{+}.}\hfill \end{array}$$(61)

When the atmospheric effects vanish, the correction (59) reduces to the vacuum correction of Blanchet et al. (2001), as can be verified by transforming this correction into the rotating coordinates.

The examples in the next section focus on the numerical evaluation of the atmospheric effects originating from the spherically symmetric, static part of the refractive index $\widehat{n}\left(r\right)$ (the first three lines of Eq. (60)) and from the wind (the fourth line of Eq. (60)). The estimation of the remaining atmospheric terms describing the effects of nonsphericity and nonstationarity of the refractive index and the nonstationarity of the wind field is left for future research.

4.7 Frequency transfer. Examples and magnitudes of the effects

4.7.1 Effect of the atmospheric refractivity. The spherical part

First, we evaluate the contribution of the spherical part of the atmospheric refractive index $\widehat{n}\left(r\right)$ to the two-way frequency transfer correction (59), which is given by the first three lines of Eq. (60). We denote $${\text{Δ}}_{\text{At}}^{\left(1\right)}=\frac{1}{c^2}{\upsilon }_{Bi}{\chi }_{+}^{i}D{\omega }_j{\gamma }_{+}^j{\epsilon }_B,$$(62a) $${\text{Δ}}_{\text{At}}^{\left(2\right)}=\frac{1}{c^2}{\upsilon }_{Bi}\left({\delta }^{ij}-{\chi }_{+}^i{\chi }_{+}^j\right){ϵ}^b{}_{jk}{\omega }^k\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_b{\widehat{n}}_{+}},$$(62b) $${\text{Δ}}_{\text{At}}^{\left(3\right)}=\frac{1}{c^2}{\upsilon }_{Bi}\left({\delta }^{ia}-{\chi }_{+}^i{\chi }_{+}^a\right){ϵ}^b{}_{jk}{\chi }_{+}^j{\omega }^k\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\widehat{l}}^2\left(\widehat{l}-\widehat{L}\right){\partial }_a{\partial }_b{\widehat{n}}_{+}},$$(62c) $${\text{Δ}}_{\text{At}}^{\left(S\right)}={\text{Δ}}_{\text{At}}^{\left(1\right)}+{\text{Δ}}_{\text{At}}^{\left(2\right)}+{\text{Δ}}_{\text{At}}^{\left(3\right)}.$$(62d)

As an example, we consider an observer on the ground, located at the equator with r_A = 6371 km and a satellite moving on a circular orbit in the equatorial plane at a radius of r_B = r_A+408 km in the direction of the rotation of the Earth. The corresponding satellite velocity in the co-rotating system is υ_B = 7.17 km s⁻¹. We consider the spherical refractive index field given by Eq. (51), which corresponds to the isothermal atmosphere in hydrostatic equilibrium with a constant air composition, and we use the same parameters as given in the initial part of Sect. 4.5.1 up to Eq. (52).

The resulting corrections ${\text{Δ}}_{\text{At}}^{\left(1\right)}$ to ${\text{Δ}}_{\text{At}}^{\left(3\right)}$ and their sum ${\text{Δ}}_{\text{At}}^{\left(S\right)}$ as functions of the satellite zenith angle θ are shown in Fig. 5. One value of θ corresponds to two possible positions of the satellite at the reception event B: one before the satellite zenith, and one after. The value of all the three corrections is the same for both satellite positions. The corrections ${\text{Δ}}_{\text{At}}^{\left(2\right)}$ and ${\text{Δ}}_{\text{At}}^{\left(3\right)}$ were calculated by transforming the integration variable into r using Eq. (36).

Figure 5 shows that the sum of the corrections starts at 10⁻¹⁷ for θ = 0, reaches 10⁻¹⁶ for θ ≈ 56°, and ends at 10⁻¹³ for θ approaching 90°. Therefore, the corrections should be taken into account in the forthcoming clock-on-satellite experiments.

4.7.2 Effect of the wind

In this section, we evaluate the effect of the wind, which enters the correction (59) through the fourth line of Eq. (60). In case of a nonstationary A_i, another effect is given by the first line of Eq. (61), which is not discussed here. We denote $${\text{Δ}}_{\text{At}}^{\text{wind}}=\frac{1}{c^2}{\upsilon }_B^i{\chi }_{+}^j{ϵ}_{ijk}\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\text{\hspace{0.17em}curl}{A}_{+}^k}.$$(63)

To show how large the effect can be, we evaluate Eq. (63) for the simple situation with the satellite position $x_B^i$ radially above the emission point on the ground $x_{A1}^i$, in which case, the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ is just a straight vertical line connecting $x_{A1}^i$ with $x_B^i$ (i.e., ${\widehat{x}}_{+}^i\left(\widehat{l}\right)=x_{A1}^i+{\chi }_{+}^i\widehat{l}$ and $\widehat{L}=r_B-r_A$). Next, we assume that the wind field V_j at t = t_B is perpendicular to the path tangent ${\chi }_{+}^i$ in the neighborhood of ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ (horizontal flow). In this case, using Eq. (10) and writing $${\chi }_{+}^j{ϵ}_{ijk}\text{curl}{A}_{+}^k={\chi }_{+}^j{\partial }_i{A}_{j+}-{\chi }_{+}^j{\partial }_j{A}_{i+},$$(64)

the first term on the RHS of Eq. (64) vanishes, and the second term, which equals $-\text{d}{A}_{i+}/\text{d}\widehat{l},$ can be integrated by parts to obtain $${\text{Δ}}_{\text{At}}^{\text{wind}}=\frac{1}{c^2}{\upsilon }_B^i\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{A}_{i+}},$$(65)

where the term containing $A_{i+}{|}_{\widehat{L}}=A_B^i$, which is on the order of c⁻²N_B, was neglected in the integration by parts. Expansion of Eq. (10) in the refractivity N gives $$A_i=-2N{V}_i+O\left(2\right).$$(66)

For the numerical evaluation of Eq. (65), we further assume that in the neighborhood of the base path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$, the wind speed Vⁱ is constant and has a direction opposite to the satellite velocity ${\upsilon }_B^i$, which is itself perpendicular to ${\chi }_{+}^i$ (i.e., $-{\upsilon }_B^i{V}_i={\upsilon }_{B}V$ with ${\upsilon }_B=\left|{\upsilon }_B^i\right|$ and V = ∣Vⁱ∣). Next, we consider the refractive index field approximated by Eq. (51). For the position of the observer on the ground, r_A = 6371 km, and the satellite position r_B = r_A + 408 km, we obtain $${\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{N}_{+}}\approx {\displaystyle {\int }_{r_A}^{r_B}\text{d}r\widehat{N}\left(r\right)}\approx 2.32\text{\hspace{0.17em}m.}$$(67)

This value can also be used to estimate the effect for higher satellite positions because the refractivity of air is negligible for altitudes above 400 km.

Finally, considering the satellite velocity magnitude in the corotating frame υ_B = 7.36 km s⁻¹, which approximately corresponds to the ISS, Eq. (65) gives $${\text{Δ}}_{\text{At}}^{\text{wind}}=0.93\times {10}^{-18}{\left(\text{m\hspace{0.17em}}{\text{s}}^{-1}\right)}^{-1}\times V.$$(68)

This means that a horizontal wind speed V > 11 m s⁻¹ can cause an effect that exceeds the 10⁻¹⁷ level.

We also note that in case of negligible air composition variations, Eq. (A.3) holds, and using Eq. (66), we obtain $$A_i=-2\frac{N_0}{{\rho }_0}\rho V_i+O\left(2\right),$$(69)

with ρ being the air density field and N₀, ρ₀ being air refractivity and density in a certain spacetime point. At the same time, the mass of air that passes through a 2D spatial surface Ω per unit of time (mass flow rate) is given as $$\dot{m}={\displaystyle {\int }_{\text{Ω}}\rho V_i{n}^i}\text{dΩ,}$$(70)

with nⁱ being a unit normal field of the surface, and dΩ its volume element (both in the Euclidean sense). Therefore, the correction (65) is related to the mass flow rate of air through a narrow band Ω defined by dragging the vertical path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ in the direction perpendicular to both ${\upsilon }_B^i$ and ${\chi }_{+}^i$ by the small Euclidean distance d. In this case, taking nⁱ in the direction of ${\upsilon }_B^i$ Eq. (65) leads to $${\text{Δ}}_{\text{At}}^{\text{wind}}=\frac{2}{c^2}\frac{N_0}{{\rho }_0}{\upsilon }_B\underset{d\to 0}{\lim }\frac{\dot{m}}{\widehat{L}d}.$$(71)

However, this relation between the correction ${\text{Δ}}_{\text{At}}^{\text{wind}}$ and the mass flow rate ṁ is not valid for general wind fields and base paths, and it requires special conditions, such as those assumed to derive Eq. (65).

5 Theory

In this section, the underlying theory is described, and a complete derivation of the results presented in Sect. 4 is given. This section may be skipped when derivation details are not of interest.

5.1 Light rays in flowing media

We assume that an electromagnetic signal propagates in spacetime with the metric g_αβ, filled with a medium that moves with a four-velocity field U^α satisfying U^α U_α = −1. Next, we assume that the electromagnetic properties of the medium are linear, isotropic, transparent, and nondispersive, such that they are described by two scalar functions, the permittivity ϵ(x^α), and the permeability µ(x^α). The refractive index of the medium is then defined as $n\left(x^{\alpha }\right)=c\sqrt{ϵ\left(x^{\alpha }\right)\mu \left(x^{\alpha }\right)}$.

According to Gordon (1923), the Maxwell equations in the medium can be reformulated using the optical metric, which is defined as $${\overline{g}}_{\alpha \beta }\equiv g_{\alpha \beta }+\left(1-\frac{1}{n^2}\right)U_{\alpha }{U}_{\beta }.$$(72)

A summary of this reformulation and of the geometrical optics approximation is given, for example, by Chen & Kantowski (2008), and we follow this reference in this paragraph. A detailed treatment of the ray optics in relativistic media is given by Synge (1960) and Perlick (2000).

We denote by ∇_α the covariant derivative associated with the metric g_αβ and ${\overline{\nabla }}_{\alpha }$ the covariant derivative associated with the optical metric ${\overline{g}}_{\alpha \beta }$. We introduce a notation with the bar for tensors, obtained by lowering or raising indices with the optical metric ${\overline{g}}_{\alpha \beta }$ or its inverse ${\overline{g}}^{\alpha \beta }$. For example, for the electromagnetic tensor, we have ${\overline{g}}^{\alpha \mu }{\overline{g}}^{\beta v}{F}_{\mu v}={\overline{F}}^{\alpha \beta }$. The Maxwell equations in the medium with vanishing free four-current then read $${\partial }_{[\alpha }{F}_{\beta \gamma ]}=0,{\overline{\nabla }}_{\beta }\left(\sqrt{ϵ/\mu }{\overline{F}}^{\beta \alpha }\right)=0.$$(73)

We use the ansatz of geometrical optics, $$F_{\alpha \beta }=\Re \left\{{\text{e}}^{iS/\sout{\lambda }}\left(A_{\alpha \beta }+O\left(\sout{\lambda }\right)\right)\right\},$$(74)

where $\sout{\lambda }$ is a wavelength related parameter, S (x^α) is the so-called eikonal function, and 𝔑 denotes the real part. We define $$k_{\alpha }={\partial }_{\alpha }S.$$(75)

The Maxwell Eqs. (73) then give to the order ${\sout{\lambda }}^{-1}$ $$A_{[\alpha \beta }{k}_{\gamma ]}=0,{\overline{A}}^{\alpha \beta }{k}_{\beta }=0.$$(76)

Contracting the first equation with ${\overline{g}}^{\gamma v}{k}_v$ and using the second equation, we obtain $${\overline{g}}^{\alpha \beta }{k}_{\alpha }{k}_{\beta }=0,$$(77)

therefore, ${\overline{k}}^{\alpha }$ is a null vector of the optical metric. At the same time, from Eq. (75), it follows that $${\partial }_{[\alpha }{k}_{\beta ]}={\overline{\nabla }}_{[\alpha }{k}_{\beta ]}=0.$$(78)

Contracting Eq. (78) with ${\overline{k}}^{\alpha }$ and using Eq. (77), we obtain $${\overline{k}}^{\alpha }{\overline{\nabla }}_{\alpha }{\overline{k}}^{\beta }=0,$$(79)

which means that ${\overline{k}}^{\alpha }$ is tangent to the null geodesic of the optical metric. The electromagnetic signal trajectory or the light ray is defined as the integral line of the field ${\overline{k}}^{\alpha }$. Therefore, it is the null geodesic of the optical metric.

To derive the equation of motion for a light ray, we might use the geodesic Eq. (79), but a more convenient way is to directly start with Eqs. (77) and (78). Taking a partial derivative of Eq. (77) with respect to the spacetime coordinates, using Eq. (78), and assuming the field ${\overline{k}}^{\beta }$ to be tangent to the trajectory x^β(a) we search for, that is, $$\frac{\text{d}{x}^{\beta }}{\text{d}a}={\overline{k}}^{\beta }={\overline{g}}^{\alpha \beta }{k}_{\alpha },$$(80)

we obtain the equation of motion with an affine parameterization $$\frac{\text{d}{k}_{\gamma }}{\text{d}a}+\frac{1}{2}\left({\partial }_{\gamma }{\overline{g}}^{\alpha \beta }\right)k_{\alpha }{k}_{\beta }=0,$$(81)

where the inverse optical metric is given as $${\overline{g}}^{\alpha \beta }=g^{\alpha \beta }+\left(1-n^2\right)U^{\alpha }{U}^{\beta }.$$(82)

5.2 Frequency shift in nonstationary flowing media

The period of the electromagnetic wave given by Eq. (74) as observed by an observer moving along the trajectory x^α(τ) parameterized by its proper time τ and having a four-velocity u^α = dx^α/dcτ is given as the proper time ∆τ during which the phase $S\left(x^{\alpha }\left(\tau \right)\right)/\sout{\lambda }$ changes by 2π in a linear order (see, e.g., Synge 1960 and Chen & Kantowski 2008), so that we have $$2\pi =-\frac{1}{\sout{\lambda }}\frac{\text{d}S\left(x^{\alpha }\left(\tau \right)\right)}{\text{d}\tau }\text{Δ}\tau =-\frac{c\text{Δ}\tau }{\sout{\lambda }}{u}^{\alpha }{k}_{\alpha },$$(83)

where the minus sign is there to obtain a positive ∆τ for future oriented ${\overline{k}}^{\alpha }$. Therefore, the proper frequency ν = 1/∆τ is given as $$v=-\frac{c}{2\pi \sout{\lambda }}{u}^{\alpha }{k}_{\alpha }.$$(84)

We consider a signal emitted with the proper frequency ν_I by an observer with the four-velocity $u_I^{\alpha }$ from the spacetime point $x_I^{\alpha }$ that is received with the proper frequency ν_F by an observer with the four-velocity $u_I^{\alpha }$ in the spacetime point $x_F^{\alpha }$ lying on the signal trajectory given by the solution x^α(a) of the Eqs. (80) and (81), satisfying $x^{\alpha }\left(a_I\right)=x_I^{\alpha }$ and $x^{\alpha }\left(a_F\right)=x_F^{\alpha }$. The ratio of the frequencies then is $$\frac{v_I}{v_F}=\frac{{\left(u^{\alpha }{k}_{\alpha }\right)}_I}{{\left(u^{\alpha }{k}_{\alpha }\right)}_F},$$(85)

where (k_α)_I and (k_α)_F are the values of k_α in the emission and reception events, which can be calculated from Eq. (80) when the solution x^α(a) is known.

As the first step, we replace the affine parameter a that is given by solution of the system of Eqs. (80) and (81) by another parameter λ that is defined by an additional equation. Namely, we consider the vector field ξ^α with k_αξ^α ≠ 0 (specified below), and we define the parameter λ such that the corresponding geodesic tangent l^α = dx^α(a(λ))/dλ satisfies $${\overline{g}}_{\alpha \beta }{l}^{\alpha }{\xi }^{\beta }=-1$$(86)

all along the light ray. We denote κ(λ) = dλ/da the factor converting l^α to the affine tangent k^α, that is, $${\overline{k}}^{\alpha }=\kappa \left(\lambda \right)l^{\alpha }.$$(87)

The frequency ratio then is $$\frac{v_I}{v_F}=\frac{{\left(u^{\alpha }{\overline{l}}_{\alpha }\right)}_I}{{\left(u^{\alpha }{\overline{l}}_{\alpha }\right)}_F}\frac{{\kappa }_I}{{\kappa }_F}.$$(88)

Now we express the ratio κ_I/κ_F in terms of the fields along the signal trajectory. The geodesic Eq. (79) gives $${\overline{k}}^{\alpha }{\overline{\nabla }}_{\alpha }{\overline{k}}^{\beta }=\left(\kappa l^{\alpha }{\overline{\nabla }}_{\alpha }\kappa \right)l^{\beta }+{\kappa }^2{l}^{\alpha }{\overline{\nabla }}_{\alpha }{l}^{\beta }=0$$(89)

or $$\frac{\text{d}\kappa }{\text{d}\lambda }{l}^{\beta }+\kappa l^{\alpha }{\overline{\nabla }}_{\alpha }{l}^{\beta }=0.$$(90)

Contracting Eq. (90) with ${\overline{\xi }}_{\beta }$, we obtain $$\begin{array}{l}\frac{1}{\kappa }\frac{\text{d}\kappa }{\text{d}\lambda }=l^{\alpha }{\overline{\xi }}_{\beta }{\overline{\nabla }}_{\alpha }{l}^{\beta }\hfill \\ =-l^{\alpha }{l}^{\beta }{\overline{\nabla }}_{\alpha }{\overline{\xi }}_{\beta },\hfill \end{array}$$(91)

where in the second line we used $l^{\alpha }{\overline{\nabla }}_{\alpha }\left(l^{\beta }{\overline{\xi }}_{\beta }\right)=0$, which follows from Eq. (86). Taking into account that the Lie derivative of ${\overline{g}}_{\alpha \beta }$ with respect to the field ξ^α is given as $$L_{\xi }{\overline{g}}_{\alpha \beta }={\overline{\nabla }}_{\alpha }{\overline{\xi }}_{\beta }+{\overline{\nabla }}_{\beta }{\overline{\xi }}_{\alpha },$$(92)

we obtain $$\frac{1}{\kappa }\frac{\text{d}\kappa }{\text{d}\lambda }=-\frac{1}{2}{l}^{\alpha }{l}^{\beta }{L}_{\xi }{\overline{g}}_{\alpha \beta }.$$(93)

Integrating Eq. (93) along the signal trajectory from λ_I to λ_F, we obtain $$\frac{{\kappa }_I}{{\kappa }_F}=\exp {\displaystyle \underset{{\lambda }_I}{\overset{{\lambda }_F}{\int }}\frac{1}{2}{l}^{\alpha }{l}^{\beta }{L}_{\xi }{\overline{g}}_{\alpha \beta }}d\lambda .$$(94)

Equation (88) then leads to the following result for the frequency ratio: $$\frac{v_I}{v_F}=\frac{{\left(u^{\alpha }{\overline{l}}_{\alpha }\right)}_I}{{\left(u^{\alpha }{\overline{l}}_{\alpha }\right)}_F}\exp {\displaystyle \underset{{\lambda }_I}{\overset{{\lambda }_F}{\int }}\frac{1}{2}{l}^{\alpha }{l}^{\beta }{L}_{\xi }{\overline{g}}_{\alpha \beta }}\text{\hspace{0.17em}d}\lambda .$$(95)

Next, we select ξ^α to be the coordinate basis vector field corresponding to the time coordinate (i.e., considering the coordinate system (x⁰, xⁱ) with x⁰ = ct on spacetime, we take ξ^α = (∂₀)^α). In this case, the condition (86) leads to $${\overline{l}}_0=-1.$$(96)

Denoting with υⁱ = dxⁱ /dt the coordinate velocity of an observer, the four-velocity of the observer is expressed as $$u^{\alpha }=u^0\left(1,{\upsilon }^i/c\right),$$(97)

and we obtain $$u^{\alpha }{\overline{l}}_{\alpha }=u^0\left(-1+\frac{1}{c}{\upsilon }^i{\overline{l}}_i\right).$$(98)

Moreover, expressing all tensor components in the coordinate basis of the selected coordinate system, the Lie derivative in the integral reduces to the partial derivative with respect to x⁰. Equation (95) then gives $$\frac{v_I}{v_F}=\frac{u_I^0}{u_F^0}\frac{1-\frac{1}{c}{\left({\upsilon }^i{\overline{l}}_i\right)}_I}{1-\frac{1}{c}{\left({\upsilon }^i{\overline{l}}_i\right)}_F}\exp {\displaystyle \underset{{\lambda }_I}{\overset{{\lambda }_F}{\int }}\frac{1}{2}{l}^{\alpha }{l}^{\beta }{\partial }_0{\overline{g}}_{\alpha \beta }}\text{\hspace{0.17em}d}\lambda .$$(99)

This is an analogy of the formula (A.46) of Blanchet et al. (2001), generalized to the nonstationary spacetime filled with a non-stationary medium. We also note that for ξ^α = (∂₀)^α, we have κ = −k₀ ${\widehat{l}}_i=-k_i/k_0$.

The first fraction in Eq. (99) contains the gravitational red-shift and the second-order Doppler shift. The second fraction is the first-order Doppler shift in the flowing medium. Here, effects of the medium such as bending of the light ray due to the refractive index gradient or dragging of the light by the motion of the medium enter most significantly. The integral in the exponent in Eq. (99) gives the correction for nonstationarity of the optical metric, or in other words, for the time dependence of the refractive index, the flow velocity of the medium and the gravitational field along the signal path.

In order to proceed with Eq. (99) that is used for the frequency transfer and also with the calculation of the signal propagation time used for the time transfer, we now need to find the solution of the system of Eqs. (80) and (81), reparameterized to the parameter λ or another related parameter. Using Eqs. (87) and (93), we transform the system of Eqs. (80) and (81) into the following set of equations of motion in Hamiltonian form: $$\frac{\text{d}{x}^{\alpha }}{\text{d}\lambda }={\overline{g}}^{\alpha \beta }{\overline{l}}_{\beta },$$(100) $$\frac{\text{d}{\overline{l}}_{\gamma }}{\text{d}\lambda }=-\frac{1}{2}\left({\partial }_{\gamma }{\overline{g}}^{\alpha \beta }+{\overline{l}}_{\gamma }{L}_{\xi }{\overline{g}}^{\alpha \beta }\right){\overline{l}}_{\alpha }{\overline{l}}_{\beta }.$$(101)

With the choice ξ^α = (∂₀)^α, which leads to Eq. (96), the time component of Eq. (101) is satisfied identically.

In the next sections, we solve this set of equations for the conditions in the Earth atmosphere, and we express the time and frequency transfer corrections based on this solution.

5.3 Optical metric in the vicinity of the Earth

The approximation level of the model as described in Sect. 4.1 requires the Doppler terms ${\overline{l}}_i{\upsilon }^i/c$ in Eq. (99) to be evaluated including all terms on the order of 𝒫(3), and the same accuracy is required for the propagation time Δt = Δx⁰/c. For this approximation level, it is sufficient to expand the metric components ${\overline{g}}^{\alpha \beta }$ on the RHS of Eqs. (100) and (101) to include all terms with the power of c⁻¹ lower by one, which are the terms on the order of 𝒫(2). This can be deduced retrospectively from the procedure leading to the final results.

We solve the problem in the rotating coordinates, where the components of the spacetime metric are given by Eq. (3) and the components of the four-velocity of the medium U^α are $$U^{\alpha }=U^0\left(1,V^i/c\right),$$(102)

where Vⁱ(x^α) is the velocity field of the medium in the rotating coordinates given by Eq. (9), and $$U^0=1+\frac{W}{c^2}+\frac{1}{2c^2}{\left|{\upsilon }_R^i+V^i\right|}^2+O\left(c^{-4}\right).$$(103)

Using the definition (10), the contravariant components of the optical metric (82) in the rotating coordinates to the approximation level described above are $${\overline{g}}^{00}=-n^2-\frac{2W}{c^2},$$(104a) $${\overline{g}}^{0i}=\frac{1}{c}\left({\upsilon }_R^i+A^i\right),$$(104b) $${\overline{g}}^{ij}={\delta }^{ij}\left(1-\frac{2W}{c^2}\right)-\frac{1}{c^2}{\upsilon }_R^i{\upsilon }_R^j.$$(104c)

The covariant components of the optical metric (72) in the rotating coordinates to the same approximation level are $${\overline{g}}_{00}=-\frac{1}{n^2}+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{c^2},$$(105a) $${\overline{g}}_{0i}=\frac{1}{c{n}^2}\left({\upsilon }_{Ri}+A_i\right),$$(105b) $${\overline{g}}_{ij}={\delta }_{ij}\left(1+\frac{2W}{c^2}\right).$$(105c)

5.4 Equation of motion for a light ray in the atmosphere of Earth

We express Eqs. (100) and (101) using the optical metric (104) and Eq. (96). The Lie derivative in Eq. (101) is given as ${\partial }_0{\overline{g}}^{\alpha \beta }$. Expanding the RHSs of Eqs. (100) and (101) such that they include all terms on the order of 𝒫(2) (we also keep the ∂₀Aⁱ term to determine how it affects the equation of motion even if it is on the order that we otherwise neglect), we obtain $$\begin{array}{c}\frac{\text{d}{\overline{l}}_i}{\text{d}\lambda }=\frac{1}{2}{\partial }_i\left(n^2\right)+\frac{2}{c^2}{\partial }_{i}W+\frac{1}{c}{\partial }_i\left({\upsilon }_R^j\right){\overline{l}}_j\left(1+\frac{1}{c}{\upsilon }_R^k{\overline{l}}_k\right)\\ +\frac{1}{c}{\overline{l}}_j{\partial }_i{A}^j+\frac{1}{2}{\partial }_0\left(n^2\right){\overline{l}}_i+{\frac{1}{c}}_i{\overline{l}}_i{\overline{l}}_j{\partial }_0{A}^j\end{array}$$(106) $$\frac{\text{d}{x}^i}{\text{d}\lambda }={\overline{l}}_i\left(1-\frac{2W}{c^2}\right)-\frac{1}{c}\left({\upsilon }_R^i+A^i\right)-\frac{1}{c^2}{\upsilon }_R^i{\upsilon }_R^j{\overline{l}}_j$$(107) $$\frac{\text{d}{x}^0}{\text{d}\lambda }=n^2+\frac{2W}{c^2}+\frac{1}{c}\left({\upsilon }_R^i+A^i\right){\overline{l}}_i,$$(108)

where we used the fact that $${\delta }^{ij}{\overline{l}}_i{\overline{l}}_j=1+O\left(1\right),$$(109)

which can be obtained when we express ${\overline{l}}_k={\left({\delta }^{ij}{\overline{l}}_i{\overline{l}}_j\right)}^{1/2}{a}_k$ with δ^ija_ia_j = 1, and we solve the null condition ${\overline{g}}^{\alpha \beta }{\overline{l}}_{\alpha }{\overline{l}}_{\beta }=0$ for ${\left({\delta }^{ij}{\overline{l}}_i{\overline{l}}_j\right)}^{1/2}$.

Equation (107) can be inverted to give ${\overline{l}}_i$_; as a function of dxⁱ /dλ. Thus, we obtain the following formula including all terms on the order of 𝒫(2): $${\overline{l}}_i=\frac{\text{d}{x}^i}{\text{d}\lambda }\left(1+\frac{2W}{c^2}\right)+\frac{1}{c}\left({\upsilon }_R^i+A^i\right)+\frac{1}{c^2}{\upsilon }_R^i{\upsilon }_{Rj}\frac{\text{d}{x}^j}{\text{d}\lambda }.$$(110)

Differentiating Eq. (107) with respect to λ and using Eqs. (106), (108), and (110) we arrive at the second-order Newton-type equation of motion for a light ray, with the RHS including all terms on the order of 𝒫(2), $$\begin{array}{c}\frac{{\text{d}}^2{x}^i}{\text{d}{\lambda }^2}=\frac{1}{2}{\partial }_i\left(n^2\right)+\frac{2}{c^2}{\partial }_{j}W\left({\delta }^{ij}-\frac{\text{d}{x}^i}{\text{d}\lambda }\frac{\text{d}{x}^j}{\text{d}\lambda }\right)+\frac{1}{2c^2}{\partial }_i{\upsilon }_R^2\\ -\frac{1}{c}{\partial }_0{A}_j\left({\delta }^{ij}-\frac{\text{d}{x}^i}{\text{d}\lambda }\frac{\text{d}{x}^j}{\text{d}\lambda }\right)+\frac{1}{2}{\partial }_0\left(n^2\right)\frac{\text{d}{x}^i}{\text{d}\lambda }\\ +\frac{2}{c}\left({\partial }_{\left[{}_i{\upsilon }_{Rj}\right]}+{\partial }_{\left[{}_i{A}_j\right]}\right)\frac{\text{d}{x}^j}{\text{d}\lambda }\\ +\frac{1}{c^2}{\partial }_i\left({\upsilon }_{Rj}{\upsilon }_{Rk}\right)\frac{\text{d}{x}^j}{\text{d}\lambda }\frac{\text{d}{x}^k}{\text{d}\lambda },\end{array}$$(111)

where the antisymmetrization is defined as ${\partial }_{[i}{A}_{j]}=\frac{1}{2}\left({\partial }_i{A}_j-{\partial }_j{A}_i\right)$. Inserting Eq. (110) into Eq. (108), we obtain the time component of the equation of motion, with the RHS including all terms on the order of 𝒫(2), $$\frac{\text{d}{x}^0}{\text{d}\lambda }=n^2+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{c^2}+\frac{1}{c}\left({\upsilon }_{Ri}+A_i\right)\frac{\text{d}{x}^i}{\text{d}\lambda }.$$(112)

For the purpose of analyzing the two-way transfer, it is convenient to change the parameter λ in the equation of motion into the Euclidean length of the light ray l_E, which is defined by $$\frac{\text{d}{l}_E}{\text{d}\lambda }=\sqrt{{\delta }_{ij}\frac{\text{d}{x}^i}{\text{d}\lambda }\frac{\text{d}{x}^j}{\text{d}\lambda }}=\left|l^i\right|.$$(113)

In this case, the spatial part of the light ray xⁱ(l_E) has a unit tangent: $${\delta }_{ij}\frac{\text{d}{x}^i}{\text{d}{l}_E}\frac{\text{d}{x}^j}{\text{d}{l}_E}=1.$$(114)

Using lⁱ = ∣lⁱ∣dxⁱ/dl_E and solving the null-vector condition ${\overline{g}}_{\alpha \beta }{l}^{\alpha }{l}^{\beta }=0$ together with Eq. (86) for l⁰ and ∣lⁱ∣, we obtain the following formulas including all terms on the order of 𝒫(2): $$l^0=n^2+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{c^2}+\frac{n}{c}\left({\upsilon }_{Ri}+A_i\right)\frac{\text{d}{x}^i}{\text{d}{l}_E},$$(115) $$\left|l^i\right|=n+\frac{1}{2c^2}{\upsilon }_{Ri}{\upsilon }_{Rj}\left({\delta }^{ij}-\frac{\text{d}{x}^i}{\text{d}{l}_E}\frac{\text{d}{x}^j}{\text{d}{l}_E}\right).$$(116)

The transformed Eqs. (111) and (112) with RHSs including all terms on the order of 𝒫(2) then are $$\begin{array}{c}\frac{{\text{d}}^2{x}^i}{\text{d}{l}_E}=\left(\frac{1}{n}{\partial }_{j}n+\frac{2}{c^2}{\partial }_{j}W+\frac{1}{2c^2}{\partial }_j{\upsilon }_R^2-\frac{1}{c}{\partial }_0{A}_j\right)\left({\delta }^{ij}-\frac{\text{d}{x}^i}{\text{d}{l}_E}\frac{\text{d}{x}^j}{\text{d}{l}_E}\right)\\ +\frac{2}{nc}\left({\partial }_{\left[{}_i{\upsilon }_{Rj}\right]}+{\partial }_{\left[{}_i{A}_j\right]}\right)\frac{\text{d}{x}^j}{\text{d}{l}_E}\\ +\frac{1}{c^2}{\partial }_i\left({\upsilon }_{Rj}{\upsilon }_{Rk}\right)\frac{\text{d}{x}^j}{\text{d}{l}_E}\frac{\text{d}{x}^k}{\text{d}{l}_E}\end{array}$$(117)

and $$\frac{\text{d}{x}^0}{\text{d}{l}_E}=n+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{2c^2}+\frac{1}{c}\left({\upsilon }_{Ri}+A_i\right)\frac{\text{d}{x}^i}{\text{d}{l}_E}+\frac{1}{2c^2}{\upsilon }_{Ri}{\upsilon }_{Rj}\frac{\text{d}{x}^i}{\text{d}{l}_E}\frac{\text{d}{x}^j}{\text{d}{l}_E}.$$(118)

The terms in the first line of Eq. (117) correspond to optical acceleration due to the refractive index gradient, gravitational acceleration, centrifugal acceleration, and acceleration due to the time variation in the field A_j. These accelerations are projected onto a plane perpendicular to the light-ray direction by the projection tensor at the end of the line. This corresponds to the selected parameterization by the Euclidean length, which does not allow changes in “velocity” magnitude (see Eq. (114)). The second line is an analog of a magnetic force, with υ_Rj being the magnetic potential, leading to Coriolis acceleration, and A_jbeing an additional magnetic potential originating from the wind speed and the refractive index (see Eq. (10)). It is remarkable that the effect of the velocity field of the medium to light rays is similar to the effect of the magnetic potential to charged particles. This has been pointed out in the literature (see, e.g., Leonhardt & Piwnicki 1999).

We expressed the equation of motion in the coordinates rotating together with the Earth, but a simpler form of the equation in the GCRS coordinates can be obtained by selecting ${\upsilon }_R^i=0$ if needed.

5.5 Solving the equation of motion for a light ray

In this section, we solve the system of differential Eqs. (117) and (118) assuming that we know the spatial coordinates of an emitter $x_I^i$ at the event of the emission, the spatial coordinates $x_F^i$ of the receiver at the event of reception, and either the coordinate time of the emission $x_I^0=c{t}_I$ or the coordinate time of the reception $x_F^0=c{t}_F$. Denoting by L the Euclidean length of the light path between the emission and reception events, the range of the parameter l_E is l_E ∊ [0, L], and we assume the following boundary conditions: $$x^i\left(0\right)=x_I^i\text{\hspace{1em}and\hspace{1em}}{x}^i\left(L\right)=x_F^i,$$(119) $$x^0\left(0\right)=c{t}_I\text{\hspace{1em}or\hspace{1em}}{x}^0\left(L\right)=c{t}_F.$$(120)

For the sake of brevity, we write Eqs. (117) and (118) in the form $$\frac{{\text{d}}^2{x}^i}{\text{d}{l}_E^2}=E_j\left(x^{\alpha }\right)\left({\delta }^{ij}-\frac{\text{d}{x}^i}{\text{d}{l}_E}\frac{\text{d}{x}^j}{\text{d}{l}_E}\right)+E_{ij}\left(x^{\alpha }\right)\frac{\text{d}{x}^j}{\text{d}{l}_E}+E_{ijk}\left(x^i\right)\frac{\text{d}{x}^j}{\text{d}{l}_E}\frac{\text{d}{x}^k}{\text{d}{l}_E},$$(121) $$\frac{\text{d}{x}^0}{\text{d}{l}_E}=\text{Φ}\left(x^{\alpha }\right)+{\text{Φ}}_i\left(x^{\alpha }\right)\frac{\text{d}{x}^i}{\text{d}{l}_E}+{\text{Φ}}_{ij}\left(x^{\alpha }\right)\frac{\text{d}{x}^i}{\text{d}{l}_E}\frac{\text{d}{x}^j}{\text{d}{l}_E},$$(122)

where we introduced the coefficients $$E_i\left(x^{\alpha }\right)=\frac{1}{n}{\partial }_{i}n+\frac{2}{c^2}{\partial }_{i}W+\frac{1}{2c^2}{\partial }_i{\upsilon }_R^2-\frac{1}{c}{\partial }_0{A}_i,$$(123) $$\begin{array}{c}{E}_{ij}\left(x^{\alpha }\right)=\frac{2}{nc}\left({\partial }_{\left[{}_i{\upsilon }_{Rj}\right]}+{\partial }_{\left[{}_i{A}_j\right]}\right)\\ =\frac{1}{nc}{ϵ}_{ijk}\left(2{\omega }^k+\text{curl}{A}^k\right),\end{array}$$(124) $$E_{ijk}\left(x^i\right)=\frac{1}{c^2}{\partial }_i\left({\upsilon }_{Rj}{\upsilon }_{Rk}\right),$$(125) $$\text{Φ}\left(x^{\alpha }\right)=n+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{2c^2},$$(126) $$\text{Φ}\left(x^{\alpha }\right)=\frac{1}{c}\left({\upsilon }_{Ri}+A_i\right),$$(127) $${\text{Φ}}_{ij}\left(x^{\alpha }\right)=\frac{1}{2c^2}{\upsilon }_{Ri}{\upsilon }_{Rj}.$$(128)

5.5.1 Split of the spatial equation into a spherically symmetric static part and a correction

The strategy for the approximate solution of the spatial part of the system given by Eq. (121) with the boundary conditions (119) is to split the RHS of Eq. (121) into a spherically symmetric static part plus a correction to solve the spherically symmetric static problem and to formulate and solve the equation for a correction of this solution up to a linear order in this correction. Using Eqs. (4) and (6), we decompose Eq. (123) as $$E_i\left(x^{\alpha }\right)={\widehat{E}}_i\left(x^j\right)+\text{Δ}{E}_i\left(x^{\alpha }\right),$$(129)

with $${\widehat{E}}_i=\frac{1}{\widehat{n}}{\partial }_i\widehat{n}+\frac{2}{c^2}{\partial }_i\widehat{W},$$(130) $$\text{Δ}{\widehat{E}}_i=\frac{{\partial }_i\alpha }{1+\alpha }+\frac{2}{c^2}{\partial }_i\text{Δ}W+\frac{1}{2c^2}{\partial }_i{\upsilon }_R^2-\frac{1}{c}{\partial }_0{A}_i.$$(131)

We denote by ${\widehat{x}}^i\left(\widehat{l}\right)$ the solution of the equation $$\frac{{\text{d}}^2{\widehat{x}}^i}{d{\widehat{l}}^2}={\widehat{E}}_j\left({\widehat{x}}^k\left(\widehat{l}\right)\right)\left({\delta }^{ij}-\frac{\text{d}{\widehat{x}}^i}{d\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{d\widehat{l}}\right),$$(132)

with $\widehat{l}\in \left[0,\widehat{L}\right]$ satisfying the boundary conditions $${\widehat{x}}^i\left(0\right)=x_I^i-\text{Δ}{x}_I^i,{\widehat{x}}^i\left(\widehat{L}\right)=x_F^i-\text{Δ}{x}_F^i,$$(133)

where we admit certain deviations $\text{Δ}{x}_I^i$, $\text{Δ}{x}_F^i$ of the boundary points from the boundary points of the complete solution given by Eq. (119), which will be useful for two-way time and frequency transfer. Practically, we use $\text{Δ}{x}_I^i=0$ and either $\text{Δ}{x}_F^i=0$ or $\text{Δ}{x}_F^i=x_{A2}^i-x_{A1}^i$, which is a spatial shift of the observer on the ground in corotating frame during the two-way back and forth propagation of the signal.

Since Eq. (132) is a special case of Eq. (117) with vanishing fields α, ∆W, ${\upsilon }_R^i$, and Aⁱ, the condition (114) also applies to the solution of Eq. (132), so that we have $${\delta }_{ij}\frac{\text{d}{\widehat{x}}^i}{d\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{d\widehat{l}}=1,$$(134)

which means that ${\widehat{x}}^i\left(\widehat{l}\right)$ is the path of a light ray between the points $x_I^i-\text{Δ}{x}_I^i$ and $x_F^i-\text{Δ}{x}_F^i$ for vanishing α, ∆W, ${\upsilon }_R^i$, and Aⁱ, parameterized by its Euclidean length, and $\widehat{L}$ is the total Euclidean length of the path, which can differ from L. Reparameterizing the solution xⁱ(l_E) of Eq. (121) to the range [0, $\widehat{L}$] by introducing a parameter $\widehat{l}=\left(\widehat{L}/L\right)l_E$ on xⁱ(l_E), the solution of the complete Eq. (121) can be expressed as $$x^i\left(\widehat{l}L/\widehat{L}\right)={\widehat{x}}^i\left(\widehat{l}\right)+\text{Δ}{x}^i\left(\widehat{l}\right),$$(135)

which defines the correction $\text{Δ}{x}^i\left(\widehat{l}\right)$. Similarly, if x⁰(l_E) is the solution of Eq. (122) for the time coordinate, the reparameterized solution can be expressed as $$x^0\left(\widehat{l}L/\widehat{L}\right)={\widehat{x}}^0+\text{Δ}{x}^0\left(\widehat{l}\right)=c\widehat{t}+c\text{Δ}t\left(\widehat{l}\right),$$(136)

where ${\widehat{x}}^0=c\widehat{t}$ is a constant with $\widehat{t}$ ∊ [t_I, t_F] and $\text{Δ}t\left(\widehat{l}\right)$ is the function to be found.

The next step is to express the differential equation for the correction $\text{Δ}{x}^i\left(\widehat{l}\right)$. To do this, we need to express the ratio $L/\widehat{L}$ in order to proceed with the derivatives $\text{d}/\text{d}{l}_E=\left(\widehat{L}/L\right)\text{d}/\text{d}\widehat{l}$. Taking a derivative of Eq. (135) with respect to $\widehat{l}$ and using Eqs. (114) and (134), we obtain $$\frac{L}{\widehat{L}}=\sqrt{1+2{\delta }_{ij}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{dΔ}{x}^j}{\text{d}\widehat{l}}+{\delta }_{ij}\frac{\text{dΔ}{x}^i}{\text{d}\widehat{l}}\frac{\text{dΔ}{x}^j}{\text{d}\widehat{l}}}.$$(137)

Inserting Eqs. (135) and (136) into Eq. (121), expanding the equation up to the linear order in ∆x^α and $\text{dΔ}{x}^i/\text{d}\widehat{l}$ and using Eqs. (129) and (132) and the linear expansion of Eq. (137) in $\text{dΔ}{x}^i/\text{d}\widehat{l}$ we arrive at the equation of the form $$\frac{{\text{d}}^2\text{Δ}{x}^i}{\text{d}{\widehat{l}}^2}=S^i\left(\widehat{l}\right)+B_j^i\left(\widehat{l}\right)\text{Δ}{x}^j+C_j^i\left(\widehat{l}\right)\frac{\text{dΔ}{x}^j}{\text{d}\widehat{l}},$$(138)

with the boundary conditions $$\text{Δ}{x}^i\left(0\right)=\text{Δ}{x}_I^i,\text{\hspace{1em}Δ}{x}^i\left(\widehat{L}\right)=\text{Δ}{x}_F^i,$$(139)

which follow from Eqs. (119), (133), and (135), and with the coefficients given by $$S^i\left(\widehat{l}\right)=S_0^i\left(\widehat{l}\right)+S_t^i\left(\widehat{l}\right)\text{Δ}t\left(\widehat{l}\right),$$(140) $$\begin{array}{c}{S}_0^i\left(\widehat{l}\right)=\text{Δ}{E}_j\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\left({\delta }^{ij}-\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\right)\\ +E_{ij}\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\\ +E_{ijk}\left({\widehat{x}}^i\left(\widehat{l}\right)\right)\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}},\end{array}$$(141) $$\begin{array}{c}{S}_t^i\left(\widehat{l}\right)={\partial }_t\text{Δ}{E}_j\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\left({\delta }^{ij}-\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\right)\\ +{\partial }_t{E}_{ij}\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}},\end{array}$$(142) $$\begin{array}{c}{B}_j^i\left(\widehat{l}\right)={\partial }_j{E}_k\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\left({\delta }^{ik}-\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\right)\\ +{\partial }_j{E}_{ik}\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\\ +{\partial }_j{E}_{ikl}\left({\widehat{x}}^i\left(\widehat{l}\right)\right)\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^l}{\text{d}\widehat{l}},\end{array}$$(143) $$\begin{array}{c}{C}_j^i\left(\widehat{l}\right)=E_k\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\left(2{\delta }^{ik}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}-{\delta }_j^k\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}-{\delta }_j^i\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\right)\\ +E_{ik}\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\left({\delta }_j^k+\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\right)\\ +2E_{ijk}\left({\widehat{x}}^i\left(\widehat{l}\right)\right)\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}.\end{array}$$(144)

We proceed with the solution of the particular equations. The solution ${\widehat{x}}^i\left(\widehat{l}\right)$ of the zeroth-order Eq. (132) with Ê_j given by Eq. (130) has been discussed in Sect. 4.3. We therefore continue with the solution of Eq. (138) for the spatial path correction $\text{Δ}{x}^i\left(\widehat{l}\right)$ with the boundary conditions (139).

5.5.2 Solution for $\text{Δ}{x}^i\left(\widehat{l}\right)$

Integrating Eq. (138) twice, we obtain $$\text{Δ}{x}^i\left(\widehat{l}\right)=c_1^i+c_2^i\widehat{l}+{\displaystyle {\int }_0^{\widehat{l}}\text{d}\widehat{l}\prime {\displaystyle {\int }_0^{\widehat{l}\prime }\text{d}\widehat{l}″\left[S^i+B_j^i\text{Δ}{x}^j+C_j^i\frac{\text{dΔ}{x}^j}{\text{d}\widehat{l}}\right]\left(\widehat{l}″\right),}}$$(145)

with $c_1^i$, $c_2^i$ being integration constants, and [ ]$\left(\widehat{l}\prime \prime \right)$ denoting that the function in the square brackets is evaluated in $\widehat{l}\prime \prime$. Determining $c_1^i$, $c_2^i$ from the boundary conditions (139), integrating by parts using $${\displaystyle {\int }_0^{\widehat{l}\prime }\text{d}\widehat{l}″f^i\left(\widehat{l}″\right)=\frac{\text{d}}{\text{d}\widehat{l}\prime }\left(\widehat{l}\prime {\displaystyle {\int }_0^{\widehat{l}\prime }\text{d}\widehat{l}″f^i\left(\widehat{l}″\right)}\right)-\widehat{l}\prime f^i\left(\widehat{l}\prime \right)}$$(146)

and expressing the integrals over the range $\widehat{l}\prime$ ∊ [0,$\widehat{l}$] or $\widehat{l}\prime$ ∊ $\left[\widehat{l},\widehat{L}\right]$ as integrals over the full range $\widehat{l}\prime$ ∊ [0, $\widehat{L}$] using the Heaviside step function $H\left(\widehat{l}\prime -\widehat{l}\right)$, which is defined as $H\left(\widehat{l}\prime -\widehat{l}\right)=0$ for $\widehat{l}\prime <\widehat{l}$, $H\left(\widehat{l}\prime -\widehat{l}\right)=1$ for $\widehat{l}\prime >\widehat{l}$ and $H\left(\widehat{l}\prime -\widehat{l}\right)=1/2$ for $\widehat{l}\prime =\widehat{l}$, we obtain $$\begin{array}{c}\text{Δ}{x}^i\left(\widehat{l}\right)=\text{Δ}{x}_I^i+\left(\text{Δ}{x}_F^i-\text{Δ}{x}_I^i\right)\frac{\widehat{l}}{\widehat{L}}\\ +{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime \{\left(\left(\widehat{l}/\widehat{L}-1\right)\widehat{l}\prime +\left(\widehat{l}\prime -\widehat{l}\right)H\left(\widehat{l}\prime -\widehat{l}\right)\right)}\\ \times \left[S^i+B_j^i\text{Δ}{x}^j+C_j^i\frac{\text{dΔ}{x}^j}{\text{d}\widehat{l}}\right]\left(\widehat{l}\right)\}.\end{array}$$(147)

Integrating by parts again, we can express the term containing $\text{d}\text{Δ}{x}^j/\text{d}\widehat{l}$ in terms of ∆x^j. Thus, we obtain the equation of the form $$\begin{array}{c}\text{Δ}{x}^i\left(\widehat{l}\right)=\text{Δ}{x}_I^i+\left(\text{Δ}{x}_F^i-\text{Δ}{x}_I^i\right)\frac{\widehat{l}}{\widehat{L}}+{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime D\left(\widehat{l},\widehat{l}\prime \right)S^i\left(\widehat{l}\prime \right)}\\ +{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime K_j^i\left(\widehat{l},\widehat{l}\prime \right)\text{Δ}{x}^j\left(\widehat{l}\prime \right)}\end{array}$$(148)

with $$D\left(\widehat{l},\widehat{l}\prime \right)=\left(\widehat{l}/\widehat{L}-1\right)\widehat{l}\prime +\left(\widehat{l}\prime -\widehat{l}\right)H\left(\widehat{l}\prime -\widehat{l}\right),$$(149) $$\begin{array}{c}{K}_j^i\left(\widehat{l},\widehat{l}\prime \right)=D\left(\widehat{l},\widehat{l}\prime \right)\left[B_j^i-\frac{\text{d}}{\text{d}\widehat{l}}{C}_j^i\right]\left(\widehat{l}\prime \right)\\ +\left(-\widehat{l}/\widehat{L}+H\left(\widehat{l}-\widehat{l}\prime \right)\right)C_j^i\left(\widehat{l}\prime \right).\end{array}$$(150)

The functions have the following properties, which will be useful later: $$D\left(0,\widehat{l}\prime \right)=D\left(\widehat{L},\widehat{l}\prime \right)=0,$$(151) $$K_j^i\left(0,\widehat{l}\prime \right)=K_j^i\left(\widehat{L},\widehat{l}\prime \right)=0,$$(152) $$D\left(\widehat{l},\widehat{l}\prime \right)=D\left(\widehat{l}\prime ,\widehat{l}\right),$$(153) $$D\left(\widehat{L}-\widehat{l},\widehat{L}-\widehat{l}\prime \right)=D\left(\widehat{l}\prime ,\widehat{l}\right).$$(154)

The solution of Eq. (148) can be found in terms of infinite series corresponding to an infinite number of iterations of Eq. (148). To be able to evaluate convergence of this series, we need to introduce the notion of a norm on the space of the possible functions $\text{Δ}{x}^i\left(\widehat{l}\right)$, and then we can use some results of the linear operator theory in Banach spaces (see, e.g., Naylor & Sell 1971).

We define the vector space of all continuous functions $f^i:\left[0,\widehat{L}\right]\to {ℝ}^3$, and we define the norm on this vector space as $$\Vert f^i\Vert =\underset{\widehat{l}\in \left[0,\widehat{L}\right]}{\max }\left(\max \left\{\left|f^1\left(\widehat{l}\right)\right|,\left|f^2\left(\widehat{l}\right)\right|,\left|f^3\left(\widehat{l}\right)\right|\right\}\right),$$(155)

where f¹, f², f³ are the corresponding component functions $\left[0,\widehat{L}\right]\to ℝ$, which are all continuous. This vector space together with the norm (155) forms a Banach space that we denote B³ (see, e.g., Naylor & Sell 1971 for the one-dimensional case). The subspace of B³ that is formed by functions satisfying $f^i\left(0\right)=f^i\left(\widehat{L}\right)=0_{{ℝ}^3}$ is denoted ${\text{B}}_0^3$.

The norm on the Banach space B³ induces an operator norm on the linear operators from B³ to itself. We denote by ${𝒜}_j^i:{\text{B}}^3\to {\text{B}}^3$ a bounded linear operator on B³, and by ${𝒜}_j^i{f}^j$ the image of a vector fⁱ ∊ B³. The operator norm $\Vert A_j^i\Vert$ of the operator is defined as $$\Vert {𝒜}_j^i\Vert =\sup \left\{\frac{\Vert {𝒜}_j^i{f}^j\Vert }{\Vert f^i\Vert }:f^i\in {\text{B}}^3\text{with\hspace{0.17em}}{f}^i\ne 0_{{}^{{\text{B}}^3}}\right\},$$(156)

where on the RHS, the norm of the Banach space is used, which is given by Eq. (155) in our case (for details, see definition 5.8.1 and lemma 5.8.2 in Naylor & Sell 1971).

Now we can proceed with Eq. (148). We assume that for any values of i, j the functions $C_j^i$ and $B_j^i-\text{d}{C}_j^i/\text{d}\widehat{l}$ that appear in Eq. (150) are continuous on $\left[0,\widehat{L}\right]$, and we define the linear operators ${𝒟}_j^i:{\text{B}}^3\to {\text{B}}_0^3$ and ${𝒦}_j^i:{\text{B}}^3\to {\text{B}}_0^3$, as $$\left({𝒟}_j^i{f}^j\right)\left(\widehat{l}\right):{\displaystyle {\int }_0^{\widehat{L}}d\widehat{l}\prime }D\left(\widehat{l},\widehat{l}\prime \right){\delta }_j^i{f}^j\left(\widehat{l}\prime \right),$$(157) $$\left({𝒦}_j^i{f}^j\right)\left(\widehat{l}\right):{\displaystyle {\int }_0^{\widehat{L}}d\widehat{l}\prime }{K}_j^i\left(\widehat{l},\widehat{l}\prime \right)f^j\left(\widehat{l}\prime \right)$$(158)

where the vanishing of the image value in $\widehat{l}=0$ and $\widehat{l}=\widehat{L}$ is ensured by the properties (151) and (152). Next, we denote by bⁱ the linear function satisfying the boundary conditions (139). This function is given as $$b^i\left(\widehat{l}\right):=\text{Δ}{x}_l^i+\left(\text{Δ}{x}_F^i-\text{Δ}{x}_I^i\right)\frac{\widehat{l}}{\widehat{L}}.$$(159)

Using these definitions and denoting by ${ℐ}_j^i$ the identity operator on B³, Eq. (148) can be written as $$\left({ℐ}_j^i-{𝒦}_j^i\right)\text{Δ}{x}^j=b^i+{𝒟}_j^i{S}^j.$$(160)

According to theorem 6.7.1 in Naylor & Sell (1971), for a continuous linear operator ${𝒦}_j^i$ satisfying $\Vert {𝒦}_j^i\Vert <1$, an inverse operator to ${ℐ}_j^i-{𝒦}_j^i$ exists, and it is given by geometric series as $${\left({ℐ}_j^i-{𝒦}_j^i\right)}^{-1}={\displaystyle \sum _{n=0}^{\infty }{\left({𝒦}_j^i\right)}^n,}$$(161)

where ${\left({𝒦}_j^i\right)}^n$ is the operator given by n consecutive applications of ${𝒦}_j^i$, so that ${\left({𝒦}_j^i\right)}^0={ℐ}_j^i,{\left({𝒦}_j^i\right)}^1={𝒦}_j^i$, and for n > 1, we have ${\left({𝒦}_j^i\right)}^n={𝒦}_{k_{n-1}}^i\dots {𝒦}_{k_2}^{k_3}{𝒦}_{k_1}^{k_2}{𝒦}_{k_j}^{k_1}$. The continuity of the linear operator is equivalent to its boundedness (see definition 5.6.3 and theorem 5.6.4 in Naylor & Sell 1971), and the boundedness of ${𝒦}_j^i$ can be proven. The assumption $\Vert {𝒦}_j^i\Vert <1$ should be verified for the physical situation of interest. For the operator (158), the operator norm (156) can be expressed in terms of the function $K_j^i\left(\widehat{l},\widehat{l}\prime \right)$ as $$\Vert {𝒦}_j^i\Vert =\max \left\{{\displaystyle \sum _{j=1}^3{\displaystyle {\int }_0^{\widehat{L}}d\widehat{l}\prime \left|K_j^i\left(\widehat{l},\widehat{l}\prime \right)\right|:\widehat{l}\in \left[0,\widehat{L}\right],i\in \left\{1,2,3\right\}}}\right\}.$$(162)

Using Eq. (161), the solution of Eq. (160) can be written as $$\begin{array}{l}\text{Δ}{x}^i={\left({ℐ}_j^i-{𝒦}_j^i\right)}^{-1}\left(b^j+{𝒟}_k^j{S}^k\right)\hfill \\ ={\displaystyle \sum _{n=0}^{\infty }{\left({𝒦}_j^i\right)}^n\left(b^j+{𝒟}_k^j{S}^k\right),}\hfill \end{array}$$(163)

where the boundary conditions (139) are ensured by the definition of bⁱ and by the properties of the operators ${𝒟}_j^i$, ${𝒦}_j^i$ given by Eqs. (151) and (152). For the derivative $\text{dΔ}{x}^i/\text{d}\widehat{l}$, which is also needed below, we obtain $$\frac{\text{dΔ}{x}^i}{\text{d}\widehat{l}}=\frac{\text{d}}{\text{d}\widehat{l}}\left(b^i+{𝒟}_j^i{S}^j\right)+\frac{\text{d}}{\text{d}\widehat{l}}\left({𝒦}_j^i{\displaystyle \sum _{n=0}^{\infty }{\left({𝒦}_k^j\right)}^n\left(b^k+{𝒟}_l^k{S}^l\right)}\right)$$(164) $$=\frac{1}{\widehat{L}}\left(\text{Δ}{x}_F^i-\text{Δ}{x}_I^i\right)+{\dot{𝒟}}_j^i{S}^j+{\dot{𝒦}}_j^i{\displaystyle \sum _{n=0}^{\infty }{\left({𝒦}_k^j\right)}^n\left(b^k+{𝒟}_l^k{S}^l\right)},$$(165)

where $$\begin{array}{l}\left({\dot{𝒟}}_j^i{S}^j\right)\left(\widehat{l}\right):=\frac{\text{d}}{\text{d}\widehat{l}}\left({𝒟}_j^i{S}^j\right)\left(\widehat{l}\right)\hfill \\ ={\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime \left(\widehat{l}\prime /\widehat{L}-H\left(\widehat{l}\prime -\widehat{l}\right)\right)S^i\left(\widehat{l}\prime \right),}\hfill \end{array}$$(166)

and the operator ${\dot{𝒦}}_j^i$ is defined as $$\begin{array}{l}\left({\dot{𝒦}}_j^i{f}^j\right)\left(\widehat{l}\right):=\frac{\text{d}}{\text{d}\widehat{l}}\left({𝒦}_j^i{f}^j\right)\left(\widehat{l}\right)\hfill \\ =C_j^i\left(\widehat{l}\right)f^j\left(\widehat{l}\right)+{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime L_j^i\left(\widehat{l},\widehat{l}\prime \right)f^j\left(\widehat{l}\prime \right),}\hfill \end{array}$$(167)

with $$L_j^i\left(\widehat{l},\widehat{l}\prime \right)=\left(\widehat{l}\prime /\widehat{L}-H\left(\widehat{l}\prime -\widehat{l}\right)\right)\left[B_j^i-\frac{d}{\text{d}\widehat{l}}{C}_j^i\right]\left(\widehat{l}\prime \right)-\frac{1}{\widehat{L}}{C}_j^i\left(\widehat{l}\prime \right).$$(168)

5.5.3 Equation for ∆t and its solution

For further analysis, it is convenient to split the solution (163) for ∆xⁱ into a part that does not depend on ∆t and a part that does depend on it, using Eq. (140). Defining $$\text{Δ}{x}_0^i={\displaystyle \sum _{n=0}^{\infty }{\left({𝒦}_j^i\right)}^n\left(b^j+{𝒟}_k^j{S}_0^k\right),}$$(169) $$\text{Δ}{x}_t^i={\displaystyle \sum _{n=0}^{\infty }{\left({𝒦}_j^i\right)}^n{𝒟}_k^j\left(S_t^k\text{Δ}t\right),}$$(170)

we can write $$\text{Δ}{x}^i=\text{Δ}{x}_0^i+\text{Δ}{x}_t^i.$$(171)

Inserting Eqs. (135) and (136) into Eq. (122), using Eq. (137) for reparameterization to the parameter $\widehat{l}$, expanding the equation as Taylor series in ∆x^α and $\text{dΔ}{x}^i/\text{d}\widehat{l}$ and using Eqs. (169)–(171), we obtain the following equation for $\text{Δ}t\left(\widehat{l}\right)$: $$c\frac{\text{dΔ}t}{\text{d}\widehat{l}}=\widehat{\sigma }\left(\widehat{l}\right)+\text{Δ}\sigma \left(\widehat{l}\right)+{\partial }_{t}n\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)\text{Δ}t,$$(172)

with $\widehat{\sigma }\left(\widehat{l}\right)$ and $\text{Δ}\sigma \left(\widehat{l}\right)$ given by $$\widehat{\sigma }\left(\widehat{l}\right)=n+\frac{1}{c}\left({\upsilon }_{Ri}+A_i\right)\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{2c^2}+\frac{1}{2c^2}{\left({\upsilon }_{Ri}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\right)}^2,$$(173) $$\begin{array}{l}\text{Δ}\sigma \left(\widehat{l}\right)={\partial }_{i}n\text{Δ}{x}_0^i+n{\delta }_{ij}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{dΔ}{x}_0^j}{\text{d}\widehat{l}}\hfill \\ +\frac{1}{2}\frac{\text{dΔ}{x}_0^i}{\text{d}\widehat{l}}\frac{\text{dΔ}{x}_0^j}{\text{d}\widehat{l}}\left({\delta }_{ij}-{\delta }_{ik}{\delta }_{jl}\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^l}{\text{d}\widehat{l}}\right)\hfill \\ +\frac{1}{c}{\upsilon }_{Ri}\frac{\text{dΔ}{x}_0^i}{\text{d}\widehat{l}}+\frac{1}{c}\left({\partial }_j{\upsilon }_{Ri}\right)\text{Δ}{x}_0^j\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}},\hfill \end{array}$$(174)

where the fields n, W, υ_Ri, and A_i and their derivatives ∂_i are evaluated in ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)$. On the RHS of Eq. (172), we kept all terms on the order of 𝒫(2), taking into account that the largest order of ∆t is c⁻¹ $\left(p_{c-1}=1,p_{\widehat{N}}=p_{\alpha }=0\right)$. This is sufficient to reach the required approximation level for the signal propagation time, as described in Sect. 4.1. All terms containing $\text{Δ}{x}_t^i$ are on a negligible order.

Now we denote ${\widehat{l}}_0$ the value of $\widehat{l}$ for which $\text{Δ}t\left({\widehat{l}}_0\right)=0$, for example, if we select $\widehat{t}=t_I$, we have ${\widehat{l}}_0=0$, or if we select $\widehat{t}=t_F$, we have ${\widehat{l}}_0=\widehat{L}$. Integrating Eq. (172) from ${\widehat{l}}_0$ to $\widehat{l}$, we obtain $$c\text{Δ}t\left(\widehat{l}\right)={\displaystyle {\int }_{{\widehat{l}}_0}^{\widehat{l}}\text{d}\widehat{l}\prime \left(\widehat{\sigma }\left(\widehat{l}\prime \right)+\text{Δ}\sigma \left(\widehat{l}\prime \right)\right)+{\displaystyle {\int }_{{\widehat{l}}_0}^{\widehat{l}}\text{d}\widehat{l}\prime {\partial }_{t}n\left({\widehat{x}}^{\alpha }\left(\widehat{l}\prime \right)\right)\text{Δ}t\left(\widehat{l}\prime \right).}}$$(175)

In principle, we could express the exact solution of Eq. (175) in terms of operator series following a similar procedure as for the spatial correction ∆xⁱ. However, an approximate solution corresponding to the approximation level of Eq. (172) itself is sufficient for our application. To obtain this solution, we express the ∆t on the RHS of Eq. (175) to its largest order c⁻¹ only:

$\text{Δ}t\left(\widehat{l}\right)\approx \left(\widehat{l}-{\widehat{l}}_0\right)/c$, which we obtain by setting $\widehat{\sigma }+\text{Δ}\sigma \approx 1$ and ∂_tn ≈ 0 in Eq. (175). Equation (175) then gives $$\text{Δ}t\left(\widehat{l}\right)=\frac{1}{c}{\displaystyle {\int }_{{\widehat{l}}_0}^{\widehat{l}}\text{d}\widehat{l}\prime \left(\widehat{\sigma }\left(\widehat{l}\prime \right)+\text{Δ}\sigma \left(\widehat{l}\prime \right)\right)+\frac{1}{c^2}{\displaystyle {\int }_{{\widehat{l}}_0}^{\widehat{l}}\text{d}\widehat{l}\prime \left(\widehat{l}\prime -{\widehat{l}}_0\right){\partial }_{t}n\left({\widehat{x}}^{\alpha }\left(\widehat{l}\prime \right)\right).}}$$(176)

5.6 Coordinate time transfer

5.6.1 One-way time transfer

In one-way time transfer, we wish to express the coordinate time t_F − t_I of the propagation of a signal that is emitted from a position $x_I^i$ at a time t_I and is received in a position $x_F^i$ at a time t_F. Based on Eqs. (120) and (136), we obtain $t_F-t_I\text{\hspace{0.17em}=\hspace{0.17em}Δ}t\left(\widehat{L}\right)-\text{Δ}t\left(0\right)$. Using Eq. (176), this leads to $$t_F-t_I=\frac{1}{c}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{\sigma }\left(\widehat{l}\right)+\text{Δ}\sigma \left(\widehat{l}\right)\right)+\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{l}-{\widehat{l}}_0\right){\partial }_{t}n\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right).}}$$(177)

Now we proceed with the integration of ∆σ. We set $\text{Δ}{x}_I^i=0$, and we assume $\text{Δ}{x}_F^i\ne 0$ in the boundary conditions (139). We keep nonvanishing $\text{Δ}{x}_F^i$ as a preparation for later use in the two-way time transfer, where we set for the back-propagating signal $\text{Δ}{x}_F^i$ equal to a shift in the position of the observer A on the ground in the corotating frame during the back and forth propagation of the signal: $\text{Δ}{x}_F^i=x_{A2}^i-x_{A1}^i\approx {\upsilon }_{A1}^i\left(2\widehat{L}/c\right)$, with ${\upsilon }_{A1}^i$ being the velocity of observer A in the corotating frame in emission event A1. Thus, we consider $\text{Δ}{x}_F^i$ to be on the order of c⁻¹. In the time-transfer formulas, however, we keep only the highest-order term containing $\text{Δ}{x}_F^i$, which is sufficient for an accuracy of 1 ps in case of a stationary ground clock for which the velocity ${\upsilon }_{A1}^i$ is given by the deformations of Earth, such as solid Earth tides, and by the nonuniformity of the rotation of Earth. In one-way transfer, we can set $\text{Δ}{x}_F^i=0$ at the end.

Using the properties (151) and (152) in Eq. (170), we obtain $\text{Δ}{x}_t^{i}0=\text{Δ}{x}_t^i\left(\widehat{L}\right)=0$, and therefore, the boundary conditions for $\text{Δ}{x}_0^i$ are the same as for $\text{Δ}{x}^i$, namely $$\text{Δ}{x}_0^i\left(0\right)=0,\text{Δ}{x}_0^i\left(\widehat{L}\right)=\text{Δ}{x}_F^i.$$(178)

Integrating terms containing $\text{dΔ}{x}_0^i/d\widehat{l}$ in ∆σ by parts with use of Eq. (178), using Eqs. (132) and (138) to express the second derivatives and using υ_Ri = ϵ_ijkω^jx^k, we obtain the following expression including all terms on the order of 𝒫(2) except the p_T = 2 order terms containing $\text{Δ}{x}_F^i$, $$\begin{array}{l}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\text{Δ}\sigma \left(\widehat{l}\right)=\frac{1}{2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\delta }_{ij}{S}_{\left(1\right)}^i\left(\widehat{l}\right)\text{Δ}{x}_0^j\left(\widehat{l}\right)}}\hfill \\ +{\delta }_{ij}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}{|}_{\widehat{L}}\text{Δ}{x}_F^j,\hfill \end{array}$$(179)

where $S_{\left(1\right)}^i$ is the first-order (p_T = 1) approximation of Sⁱ, which is given as $$S_{\left(1\right)}^i\left(\widehat{l}\right)={\partial }_j\alpha \left({\widehat{x}}^{\beta }\left(\widehat{l}\right)\right)\left({\delta }^{ij}-\frac{\text{d}{\widehat{x}}^i\left(\widehat{l}\right)}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^j\left(\widehat{l}\right)}{\text{d}\widehat{l}}\right)+\frac{2}{c}{ϵ}^i{}_{jk}{\omega }^k\frac{\text{d}{\widehat{x}}^j\left(\widehat{l}\right)}{\text{d}\widehat{l}}.$$(180)

The $\text{Δ}{x}_0^j$ in Eq. (179) can be expressed using Eq. (169). In this expansion, only the first term with $S_0^k$ approximated by $S_{\left(1\right)}^k$ is relevant for the required accuracy. Thus, we obtain $$\text{Δ}{x}_0^j\approx b^j+{𝒟}_k^j{S}_{\left(1\right)}^k,$$(181)

where $b^j=\text{Δ}{x}_F^j\widehat{l}/\widehat{L}$ does not contribute to the considered order in Eq. (179). Using the operator definition (157), we obtain the following result including all terms on the order of 𝒫(2) except the p_T = 2 order terms containing $\text{Δ}{x}_F^j$, $$\begin{array}{l}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\text{Δ}\sigma \left(\widehat{l}\right)=\frac{1}{2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\displaystyle {\int }_0^{\widehat{L}}\text{d}}\widehat{l}\prime D\left(\widehat{l},\widehat{l}\prime \right){\delta }_{ij}{S}_{\left(1\right)}^i\left(\widehat{l}\right)S_{\left(1\right)}^i\left(\widehat{l}\prime \right)}}\hfill \\ +{\delta }_{ij}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}{|}_{\widehat{L}}\text{Δ}{x}_F^j.\hfill \end{array}$$(182)

In this formula, with $S_{\left(1\right)}^i$ given by Eq. (180), the approximation (41) is sufficient to obtain the required accuracy for the resulting time-transfer formulas. In this case, using Eq. (149), some of the integrations can be performed explicitly. In the last term of Eq. (173), the approximation (41) is also sufficient. Therefore, using Eq. (177) with Eqs. (173) and (182), we obtain the following formula for the propagation time, including all terms on the order of 𝒫(3) except the p_T = 3 order terms containing $\text{Δ}{x}_F^i$ : $$\begin{array}{l}{t}_F-t_I=\frac{1}{c}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}n+\frac{1}{c^3}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}2W}\hfill \\ +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left({\upsilon }_{Ri}+A_i\right)\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}}\hfill \\ +\frac{1}{2c^3}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left({\upsilon }_R^2+{\left({\upsilon }_{Ri}{\chi }^i\right)}^2\right)}\hfill \\ -\frac{1}{6}\frac{{\widehat{L}}^3}{c^3}\left({\omega }^2-{\left({\omega }_i{\chi }^i\right)}^2\right)\hfill \\ +\frac{1}{c^2}{ϵ}^i{}_{jk}{\chi }^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_i\alpha }\hfill \\ +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{l}-{\widehat{l}}_0\right){\partial }_t\alpha }\hfill \\ +\frac{1}{2c}\left({\delta }^{ij}-{\chi }^i{\chi }^j\right){\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime D}\left(\widehat{l},\widehat{l}\prime \right){\partial }_i\alpha \left({\widehat{x}}^{\beta }\left(\widehat{l}\right)\right){\partial }_j\alpha \left({\widehat{x}}^{\beta }\left(\widehat{l}\prime \right)\right)}\hfill \\ +\frac{1}{c}{\chi }_i\text{Δ}{x}_F^i,\hfill \end{array}$$(183)

where the fields n, α, W, υ_Ri, and A_i; and their partial derivatives with respect to the spacetime coordinates are evaluated in ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)$ if not stated explicitly otherwise. Selecting $\text{Δ}{x}_F^i=0$ and ${\widehat{l}}_0=0$, which corresponds to $\widehat{t}=t_I$, approximating the terms in the third line by integration over a straight line connecting ${\widehat{x}}_I^i$ and ${\widehat{x}}_F^i$ instead of ${\widehat{x}}^i\left(\widehat{l}\right)$, which leads to a negligible difference in terrestrial conditions, and approximating $\widehat{L}=D$ in the fourth line (see Eq. (42)), we obtain Eq. (43) in the summary.

5.6.2 Definitions in two-way transfer

In two-way time transfer, and later also in two-way frequency transfer, we consider a signal that is emitted by a stationary observer on the ground A from the position $x_{A1}^i$ at the coordinate time t_A1 to an observer B, who receives the signal in the position $x_B^i$ at a coordinate time t_B and immediately sends it back to A, who receives it in the position $x_{A2}^i$ at the coordinate time t_A2. We denote the quantities related to the signal from A1 to B by the index +, whereas the quantities related to the signal from B to A2 are denoted by the index -. Using this notation, we can write t_I+ = t_A1, $x_{I+}^i=x_{A1}^i$, t_F+ = t_I− = t_B, $x_{F+}^i=x_{I-}^i=x_B^i$, t_F− = t_A2, and $x_{F-}^i=x_{A2}^i$. Next, we define ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$) the solution of Eq. (132) with the boundary conditions ${\widehat{x}}_{+}^i\left(0\right)=x_{A1}^i,{\widehat{x}}_{+}^i\left(\widehat{L}\right)=x_B^i$, and ${\widehat{x}}_{\_}^i\left(\widehat{l}\right)$ the solution of Eq. (132) with the boundary conditions ${\widehat{x}}_{-}^i\left(0\right)=x_B^i,{\widehat{x}}_{-}^i\left(\widehat{L}\right)=x_{A1}^i$. It can be shown that $${\widehat{x}}_{+}^i\left(\widehat{l}\right)={\widehat{x}}_{-}^i\left(\widehat{L}-\widehat{l}\right).$$(184)

From the boundary conditions for ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ and ${\widehat{x}}_{-}^i\left(\widehat{l}\right)$ described above, it follows that $\text{Δ}{x}_{I+}^i=\text{Δ}{x}_{F+}^i=\text{Δ}{x}_{I-}^i=0$ and $\text{Δ}{x}_{F-}^i=x_{A2}^i-x_{A1}^i$. The shift in the position of A during the back and forth propagation of the signal we denote $x_{A2}^i-x_{A1}^i:=\text{Δ}{x}_A^i$. For $\widehat{t}$ we choose $${\widehat{t}}_{+}={\widehat{t}}_{-}=t_B,$$(185)

which leads to ${\widehat{l}}_{0+}=\widehat{L}$ and ${\widehat{l}}_{0-}=0$.

5.6.3 Two-way time transfer

Now we can express the time Δt₊ := t_F+ − t_I+ by setting ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)={\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right),{\widehat{l}}_0=\widehat{L}$ and $\text{Δ}{x}_F^i=0$ in Eq. (183) and the time Δt₋ := t_F− −t_I− by setting ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)={\widehat{x}}_{-}^{\alpha }\left(\widehat{l}\right),{\widehat{l}}_0=0$ and $\text{Δ}{x}_F^i=\text{Δ}{x}_A^i$ in Eq. (183), respectively.

Denoting by $x_A^i\left(t\right)$ the trajectory of the observer A in the corotating coordinates parameterized by the coordinate time, expanding the trajectory as a Taylor series at t_A1, denoting by ${\upsilon }_{A1}^i=\text{d}{x}_A^i/\text{d}t|{}_{t_{A1}}$ and expressing t_A2 − t_A1 = Δt₊ + Δt₋ using Eq. (183), we obtain $\text{Δ}{x}_A^i$ up to the highest order, $$\begin{array}{l}\text{Δ}{x}_A^i=x_A^i\left(t_{A2}\right)-x_A^i\left(t_{A1}\right)\hfill \\ ={\upsilon }_{A1}^i\frac{2\widehat{L}}{c}+O\left(2\right).\hfill \end{array}$$(186)

Transforming the integration variables in Δt₋ as $\widehat{l}=\widehat{L}-s$, $\widehat{l}\prime =\widehat{L}-s\prime$, using $${\widehat{x}}_{-}^{\alpha }\left(\widehat{L}-S\right)={\widehat{x}}_{+}^{\alpha }\left(S\right),$$(187) $$\frac{\text{d}{\widehat{x}}_{-}^i}{\text{d}\widehat{l}}\left(\widehat{L}-S\right)=-\frac{\text{d}{\widehat{x}}_{+}^i\left(S\right)}{\text{d}s},$$(188)

which follows from Eqs. (184) and (185) and using the properties (153) and (154) of the function $D\left(\widehat{l},\widehat{l}\prime \right)$ we arrive at the following result for the two-way time transfer correction including all terms on the order of 𝒫(3) except the p_T = 3 order terms originating from $\text{Δ}{x}_A^i$: $$\begin{array}{l}\text{Δ}{t}_{-}-\text{Δ}{t}_{+}=-\frac{2}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left({\upsilon }_{Ri+}+A_{i+}\right)}\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}\hfill \\ +\frac{2}{c^2}{ϵ}^i{}_{jk}{\chi }_{+}^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{L}-\widehat{l}\right){\partial }_i{\alpha }_{+}}\hfill \\ +\frac{2}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{L}-\widehat{l}\right){\partial }_t{\alpha }_{+}}\hfill \\ -\frac{2}{c^2}\widehat{L}{\chi }_{i+}{\upsilon }_{A1}^i,\hfill \end{array}$$(189)

where the fields with the index + are evaluated in ${\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)$, for example, ${\upsilon }_{Ri+}={\upsilon }_{Ri}\left({\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)\right)$, and $x_{+}^i$ is a unit vector in the direction from $x_{A1}^i$ to $x_B^i$. This is Eq. (46) in the summary.

5.7 Frequency transfer

5.7.1 One-way frequency transfer

In this section, we proceed with the calculation of the frequency shift of the signal that propagates in the atmosphere of Earth according to Eq. (99).

An observer moving with the velocity υⁱ = dxⁱ/dt in the corotating coordinates has the four-velocity components (in the same coordinates) given by Eq. (97) with $$u^0=1+\frac{W}{c^2}+\frac{1}{2c^2}{\left|{\upsilon }_R^i+{\upsilon }^i\right|}^2+O\left(c^{-4}\right).$$(190)

Denoting by $W_I=W\left(x_I^{\alpha }\right)$, $W_F=W\left(x_F^{\alpha }\right)$, ${\upsilon }_{RI}^i={\upsilon }_R^i\left(x_I^{\alpha }\right)$, and ${\upsilon }_{RF}^i={\upsilon }_R^i\left(x_F^{\alpha }\right)$, with $x_I^{\alpha }$, $x_F^{\alpha }$ being the spacetime points of the signal emission and reception, respectively, and denoting by ${\upsilon }_I^i$, ${\upsilon }_F^i$ the observers’ velocities at the emission and reception events, we can express the first part of Eq. (99) as $$\frac{u_I^0}{u_F^0}=\frac{1+\frac{W_I}{c^2}+\frac{1}{2c^2}{\left|{\upsilon }_{RI}^i+{\upsilon }_I^i\right|}^2}{1+\frac{W_F}{c^2}+\frac{1}{2c^2}{\left|{\upsilon }_{RF}^i+{\upsilon }_F^i\right|}^2}+O\left(c^{-4}\right).$$(191)

Next, we proceed with the Doppler part of Eq. (99). The components ${\overline{l}}_i$ in the corotating coordinates can be obtained using Eqs. (110), (113), and (116). Thus, we obtain the following expression including all terms on the order of 𝒫(2): $$\begin{array}{l}{\overline{l}}_i=\frac{\text{d}{x}^i}{\text{d}{l}_E}\left(n+\frac{2W}{c^2}+\frac{{\upsilon }_R^2}{2c^2}\right)+\frac{1}{c}\left({\upsilon }_R^i+A^i\right)\hfill \\ +\frac{1}{c^2}{\upsilon }_{Rj}{\upsilon }_{Rk}\frac{\text{d}{x}^k}{\text{d}{l}_E}\left({\delta }^{ij}-\frac{1}{2}\frac{\text{d}{x}^i}{\text{d}{l}_E}\frac{\text{d}{x}^j}{\text{d}{l}_E}\right).\hfill \end{array}$$(192)

The normalized tangent to the light ray dxⁱ/dl_E can be expressed using Eqs. (135) and (137) up to the second order in $\text{dΔ}{x}^i/\text{d}\widehat{l}$ as $$\begin{array}{l}\frac{\text{d}{x}^i}{\text{d}{l}_E}=\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}+\frac{\text{dΔ}{x}^j}{\text{d}\widehat{l}}\left({\delta }_j^i-{\delta }_{jk}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\right)\hfill \\ -\frac{1}{2}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}{\delta }_{jk}\frac{\text{dΔ}{x}^j}{\text{d}\widehat{l}}\frac{\text{dΔ}{x}^k}{\text{d}\widehat{l}}\hfill \\ +\frac{3}{2}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\left({\delta }_{jk}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\frac{\text{dΔ}{x}^k}{\text{d}\widehat{l}}\right)-\frac{\text{dΔ}{x}^i}{\text{d}\widehat{l}}{\delta }_{jk}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\frac{\text{dΔ}{x}^k}{\text{d}\widehat{l}},\hfill \end{array}$$(193)

where, at the end, the terms in the third line do not contribute to the required order of ${\overline{l}}_i$. In order to obtain the boundary values ${\left({\overline{l}}_i\right)}_I$ and ${\left({\overline{l}}_i\right)}_F$, we need to express the boundary values of $\text{dΔ}{x}^i/\text{d}\widehat{l}$. We introduce the parameter η = 0 or 1, and we evaluate Eq. (165) in $\eta \widehat{L}$. We assume $\text{Δ}{x}_I^i=0$, and we keep nonvanishing $\text{Δ}{x}_F^i$ as a preparation for the two-way frequency transfer, similarly as in the case of the two-way time transfer. Denoting $$Q_j^i=B_j^i-\frac{\text{d}}{\text{d}\widehat{l}}{C}_j^i,$$(194)

we obtain $$\begin{array}{l}\frac{\text{dΔ}{x}^i}{\text{d}\widehat{l}}|{}_{\eta \widehat{L}}=\frac{1}{\widehat{L}}\Delta x_F^i+\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right)S^i\left(\widehat{l}\right)}\hfill \\ +\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\{(\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right)Q_j^i\left(\widehat{l}\right)-C_j^i\left(\widehat{l}\right))}\hfill \\ \times \left[{\displaystyle \sum _{n=0}^{\infty }{\left({𝒦}_k^j\right)}^n\left(b^k+{𝒟}_l^k{S}^l\right)}\right]\left(\widehat{l}\right)\}\hfill \\ +\eta C_j^i\left(\widehat{L}\right)\text{Δ}{x}_F^j.\hfill \end{array}$$(195)

We denote by ${\overline{l}}_{i|\eta }$ the boundary values of ${\overline{l}}_i$ namely ${\overline{l}}_{i|0}={\left({\overline{l}}_i\right)}_I$ and ${\overline{l}}_{i|1}={\left({\overline{l}}_i\right)}_F$. We write the resulting expression for ${\overline{l}}_{i|\eta }$ as the sum $${\overline{l}}_{i|\eta }={\overline{l}}_{i\left(c^0\right)|\eta }+{\overline{l}}_{i\left(c^{-1}\right)|\eta }+{\overline{l}}_{i\left(c^{-2}\right)|\eta }+{\overline{l}}_{i\left(\text{Δ}{x}_F\right)|\eta },$$(196)

including all terms on the order of 𝒫(2) except terms containing $\text{Δ}{x}_F^i$, which are given up to the largest order only. The terms in Eq. (196) are sorted according to their power of c⁻¹, as indicated by their index, and the term containing $\text{Δ}{x}_F^i$ is given separately. These terms are given by Eqs. (197a)–(197d) below, $$\begin{array}{ll}{\overline{l}}_{i\left(c^0\right)|\eta }=\hfill & {n|}_{\eta \widehat{L}}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}|{}_{\eta \widehat{L}}\hfill \\ \hfill & +n|{}_{\eta \widehat{L}}\left({\delta }_j^i-{\delta }_{jq}\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}|{}_{\eta \widehat{L}}\frac{\text{d}{\widehat{x}}^q}{\text{d}\widehat{l}}|{}_{\eta \widehat{L}}\right)\hfill \\ \hfill & \times \frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right){\partial }_k\alpha \left({\delta }^{jk}-\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\right)}\hfill \\ \hfill & +\left({\delta }^{ij}\left({\delta }^{kq}-{\chi }^k{\chi }^q\right)+{\chi }^j\left({\chi }^k{\delta }^{iq}-{\chi }^i{\delta }^{kq}\right)\right)\hfill \\ \hfill & \times \frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime \{D\left(\widehat{l},\widehat{l}\prime \right)\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right)}\hfill \\ \hfill & \times {\partial }_j{\partial }_k\left(\widehat{n}+\alpha \right)\left({\widehat{x}}^{\beta }\left(\widehat{l}\right)\right){\partial }_q\alpha \left({\widehat{x}}^{\beta }\left(\widehat{l}\prime \right)\right)\}\hfill \\ \hfill & +\left({\delta }^{ij}-{\chi }^i{\chi }^j\right)\hfill \\ \hfill & \times \frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\prime D\left(\widehat{l},\widehat{l}\prime \right){\chi }^k{\partial }_k\left(\widehat{n}+\alpha \right)\left({\widehat{x}}^{\beta }\left(\widehat{l}\right)\right){\partial }_j\alpha \left({\widehat{x}}^{\beta }\left(\widehat{l}\prime \right)\right)}}\hfill \\ \hfill & -\left({\delta }^{ij}-{\chi }^i{\chi }^j\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right)\alpha {\partial }_j\alpha }\hfill \\ \hfill & -\frac{1}{2}{\chi }^i{\left|\left({\delta }^{jk}-{\chi }^j{\chi }^k\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right){\partial }_k\alpha }\right|}^2\hfill \end{array}$$(197a) $$\begin{array}{lll}{\overline{l}}_{i\left(c^{-1}\right)|\eta }\hfill & =\hfill & \frac{1}{c}{\left({\upsilon }_R^i+A^i\right)|}_{\eta \widehat{L}}\hfill \\ \hfill & \hfill & +\frac{1}{c}{n|}_{\eta \widehat{L}}\left({\delta }^{ij}-{\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}|}_{\eta \widehat{L}}{\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}|}_{\eta \widehat{L}}\right)\hfill \\ \hfill & \hfill & \times \frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right)\frac{1}{n}{ϵ}_{jkl}\left(2{\omega }^l+\text{curl}{A}^l\right)\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\hfill \\ \hfill & \hfill & +{\widehat{n}|}_{\eta \widehat{L}}\left({\delta }_j^i-{\delta }_{jq}{\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}|}_{\eta \widehat{L}}{\frac{\text{d}{\widehat{x}}^q}{\text{d}\widehat{l}}|}_{\eta \widehat{L}}\right)\hfill \\ \hfill & \hfill & \times \frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\{(\left(\widehat{l}-\left(1-\eta \right)\widehat{L})Q_{k\left({\widehat{N}}^p\right)}^i\left(\widehat{l}\right)-C_{k\left({\widehat{N}}^p\right)}^j\left(\widehat{l}\right)\right)\hfill \\ \hfill & \hfill & \times \left[{\displaystyle \sum _{n=0}^{\infty }{\left({\mathcal{K}}_{l\left({\widehat{N}}^p\right)}^k\right)}^n{\mathcal{D}}_m^l{S}_{\left({\widehat{N}}^p{c}^{-1}\right)}^m}\right]\left(\widehat{l}\right)\}\hfill \\ \hfill & \hfill & +\left({\delta }^{ij}-{\chi }^i{\chi }^j\right)\frac{1}{c}\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right)\left(\widehat{l}-{\widehat{l}}_0\right){\partial }_t{\partial }_j\alpha }\hfill \\ \hfill & \hfill & +\frac{1}{c}{ϵ}_{klm}{\omega }^m{\chi }^l\left({\delta }^{ia}{\delta }^{kb}+{\chi }^a{\chi }^b{\delta }^{ik}\right)\hfill \\ \hfill & \hfill & \times \frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right)\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right){\partial }_a{\partial }_b\alpha }\hfill \\ \hfill & \hfill & +\frac{1}{c}{ϵ}_{mjk}{\omega }^j\left({\delta }^{im}{\delta }^{kl}-2{\chi }^k\left({\delta }^{im}{\chi }^l+{\delta }^{lm}{\chi }^i\right)\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_l\alpha }\hfill \\ \hfill & \hfill & +\frac{1}{c}{\chi }^i\left(\eta \widehat{L}-{\widehat{l}}_0\right){{\partial }_t\alpha |}_{\eta \widehat{L}}\hfill \end{array}$$(197b) $$\begin{array}{ll}{\overline{l}}_{i\left(c^{-2}\right)|\eta }=\hfill & \frac{1}{c^2}{\chi }^i\left(2W+\frac{1}{2}{\upsilon }_{Rj}{\upsilon }_{Rk}\left({\delta }^{jk}-{\chi }^j{\chi }^k\right)\right){\left|{}_{\eta \widehat{L}}+\frac{1}{c^2}\left({\upsilon }_R^i{\upsilon }_{Rj}\right)\right|}_{\eta \widehat{L}}{\chi }^j\hfill \\ \hfill & +\frac{{\widehat{L}}^2}{c^2}\left(\frac{1}{3}{\omega }^i{\chi }^k{\omega }_k+\frac{1}{6}{\chi }^i{\left({\omega }_j{\chi }^j\right)}^2-\frac{1}{2}{\chi }^i{\omega }^2\right)\hfill \\ \hfill & +\frac{1}{c^2}\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\{\left(\widehat{l}-\left(1-\eta \right)\widehat{L}\right)\hfill \\ \hfill & \times \left(\left({\delta }^{ij}-{\chi }^i{\chi }^j\right)\left(2{\partial }_j\text{Δ}W+\frac{1}{2}{\partial }_j{\upsilon }_R^2\right)+{\partial }_i\left({\upsilon }_{Rk}{\upsilon }_{Rl}\right){\chi }^k{\chi }^l\right)\}\hfill \end{array}$$(197c) $${\overline{l}}_{i\left(\text{Δ}{x}_F\right)|\eta }=\left({\delta }_j^i-{\chi }^i{\chi }_j\right)\frac{1}{\widehat{L}}\text{Δ}{x}_F^j$$(197d)

In Eqs. (197), all the fields n, $\widehat{n}$, α, υ_Ri, A_i, W, and ΔW and their partial derivatives with respect to the spacetime coordinates are evaluated in ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)$, if not stated explicitly otherwise. For example, we denote $\widehat{n}=\widehat{n}\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)$, ${\partial }_k\widehat{n}=\left({\partial }_k\widehat{n}\right)\left({\widehat{x}}^{\alpha }\left(\widehat{l}\right)\right)$. The functions evaluated in $\widehat{l}=\eta \widehat{L}$ are denoted by $|{}_{\eta \widehat{L}}$, for example, $\widehat{n}|{}_{\eta \widehat{L}}=\widehat{n}\left({\widehat{x}}^{\alpha }\left(\eta \widehat{L}\right)\right)$. The functions $Q_{j\left({\widehat{N}}^p\right)}^i$ and $C_{j\left({\widehat{N}}^p\right)}^i$ in ${\overline{l}}_{i\left(c^{-1}\right)|\eta }$ are the functions $Q_j^i$ and $C_j^i$ as given by Eqs. (143), (144), and (194) with the level of approximation including all terms with $p_{c^{-1}}=p_{\alpha }=0$ and arbitrary $p_{\widehat{N}}$, namely $$\begin{array}{c}{Q}_{j\left({\widehat{N}}^p\right)}^i\left(\widehat{l}\right)=\left(\frac{1}{\widehat{n}}{\partial }_k{\partial }_l\widehat{n}-\frac{2}{{\widehat{n}}^2}{\partial }_k\widehat{n}{\partial }_l\widehat{n}\right)\\ \times \left({\delta }^{ik}{\delta }_j^l-2{\delta }^{ik}\frac{\text{d}{\widehat{x}}^l}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}+{\delta }_j^i\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^l}{\text{d}\widehat{l}}\right)\\ +\frac{1}{{\widehat{n}}^2}{\partial }_k\widehat{n}{\partial }_l\widehat{n}\left({\delta }_j^i{\delta }^{kl}-{\delta }_j^l\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\right),\end{array}$$(198) $$C_{j\left({\widehat{N}}^p\right)}^i\left(\widehat{l}\right)=\frac{1}{\widehat{n}}{\partial }_k\widehat{n}\left(2{\delta }^{ik}\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}}-{\delta }_j^k\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}-{\delta }_j^i\frac{\text{d}{\widehat{x}}^k}{\text{d}\widehat{l}}\right),$$(199)

with $\widehat{n}$ and its partial derivatives being evaluated in ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)$. The operator ${𝒦}_{j\left({\widehat{N}}^p\right)}^i$ in ${\overline{l}}_{i\left(c^{-1}\right)|\eta }$ is then generated by the corresponding function $$K_{j\left({\widehat{N}}^p\right)}^i\left(\widehat{l},\widehat{l}\prime \right)=D\left(\widehat{l},\widehat{l}\prime \right)Q_{j\left({\widehat{N}}^p\right)}^i\left(\widehat{l}\prime \right)+\left(-\widehat{l}/\widehat{L}+H\left(\widehat{l}-\widehat{l}\prime \right)\right)C_{j\left({\widehat{N}}^p\right)}^i\left(\widehat{l}\prime \right),$$(200)

and $S_{\left({\widehat{N}}^p{c}^{-1}\right)}^i$ is Sⁱ as given by Eq. (140) with the approximation level including all terms with $p_{c^{-1}}=1$, p_α = 0 and arbitrary $p_{\widehat{N}}$, namely $$S_{({\widehat{N}}^p{c}^{-1})}^i(\widehat{l})=\frac{1}{\widehat{n}c}{ϵ}^i{}_{jk}(2{\omega }^k+\text{curl}{\widehat{A}}^k)\frac{\text{d}{\widehat{x}}^j}{\text{d}\widehat{l}},$$(201)

with ${\widehat{A}}^k=\left(1-{\widehat{n}}^2\right)V^k$, and the fields $\widehat{n}$ and curl Â^k evaluated along ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)$.

Now we proceed with the integral in the exponent in Eq. (99), which we denote I. We parameterize the signal trajectory by the Euclidean length l_E ∊ [0, L], and using Eqs. (105), (115), (116) and (118), we obtain the following expression including all terms $$\begin{array}{lll}I\hfill & \equiv \hfill & {\displaystyle {\int }_0^L\frac{1}{2}}{l}^{\alpha }{\partial }_0{\overline{g}}_{\alpha \beta }\frac{\text{d}{x}^{\beta }}{\text{d}{l}_E}\text{d}{l}_E\hfill \\ \hfill & =\hfill & {\displaystyle {\int }_0^L\text{d}}{l}_E\left(\frac{1}{c}{\partial }_{t}n+\frac{1}{c^2}\left({\partial }_t{A}_i\right)\frac{\text{d}{x}^i}{\text{d}{l}_E}+\frac{2}{c^3}{\partial }_{t}W\right)\hfill \end{array},$$(202)

where the fields ∂_tn, ∂_tA_i, and ∂_tW are evaluated along the signal trajectory x^α(l_E).

Next, we transform the integration variable to l using the first-order expansion of Eq. (137) in $\text{dΔ}{x}^i/\text{d}\widehat{l}$ and using Eqs. (135) and (136). Then, taking the Taylor expansion of the integrand in Δx^α, assuming $\text{Δ}{x}_I^i=0$, and using the approximations $$\text{Δ}{x}^i\approx \text{Δ}{x}_F^i\frac{\widehat{l}}{\widehat{L}}+{\mathcal{𝒟}}^i{}_j{S}_{\left(1\right)}^j,$$(203) $$\text{Δ}t\approx \frac{1}{c}\left(\widehat{l}-{\widehat{l}}_0\right),$$(204)

we obtain the following formula including all terms on the order of 𝒫(3) except terms containing $\text{Δ}{x}_F^i$, which are given up to the largest order only: $$I=I_{\left(c^{-1}\right)}+I_{\left(c^{-2}\right)}+I_{\left(c^{-3}\right)}+I_{\left(\text{Δ}{x}_F\right)},$$(205a) $$\begin{array}{lll}{I}_{\left(c^{-1}\right)}\hfill & =\hfill & \frac{1}{c}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\partial }_{t}n}\hfill \\ \hfill & \hfill & +\frac{1}{c}\left({\delta }^{ij}-{\chi }^i{\chi }^j\right){\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\prime D\left(\widehat{l},\widehat{l}\prime \right){\partial }_i{\partial }_t\alpha \left({\widehat{x}}^{\beta }\left(\widehat{l}\right)\right){\partial }_j\alpha \left({\widehat{x}}^{\beta }\left(\widehat{l}\prime \right)\right),\hfill \end{array}$$(205b) $$\begin{array}{lll}{I}_{\left(c^{-2}\right)}\hfill & =\hfill & \frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left({\partial }_t{A}_i\right)\frac{\text{d}{\widehat{x}}^i}{\text{d}\widehat{l}}\hfill \\ \hfill & \hfill & +\frac{1}{c^2}{ϵ}^i{}_{jk}{\chi }^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}}\left(\widehat{l}-\widehat{L}\right){\partial }_i{\partial }_t\alpha \hfill \\ \hfill & \hfill & +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left(\widehat{l}-{\widehat{l}}_0\right){\partial }_i{\partial }_t\alpha ,\hfill \end{array}$$(205c) $$I_{\left(c^{-3}\right)}=\frac{2}{c^3}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\partial }_{t}W,$$(205d) $$I_{\left(\text{Δ}{x}_F\right)}=\frac{1}{c\widehat{L}}\text{Δ}{x}_F^i{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left(\widehat{l}{\partial }_i{\partial }_t\alpha +{\chi }_i{\partial }_t\alpha \right),$$(205e)

where the fields n, α, A_i, and W and their partial derivatives with respect to the spacetime coordinates are evaluated in ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)$ if not stated explicitly otherwise. Selecting $\text{Δ}{x}_F^i=0$ and ${\widehat{l}}_0=0$, which corresponds to $\widehat{t}=t_I$, we obtain Eq. (57) in the summary.

5.7.2 Two-way frequency transfer

In two-way frequency transfer, we consider the experimental setup as already described in Sect. 5.6.2, and we admit the definitions introduced there. Furthermore, we define${\overline{l}}_{i+|\eta }$ and I₊ by setting ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)={\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)$, ${\widehat{l}}_0=\widehat{L}$, $\text{Δ}{x}_F^i=0$ in Eqs. (196)–(201) and (205), and we define ${\overline{l}}_{i-|\eta }$ and I₋ by setting ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)={\widehat{x}}_{-}^{\alpha }\left(\widehat{l}\right)$, ${\widehat{l}}_0=0$, $\text{Δ}{x}_F^i=\text{Δ}{x}_A^i$ in Eqs. (196)–(201) and (205). The proper frequencies ν and the observer velocity quantities υⁱ, u⁰ for the particular emission and reception events are denoted by the corresponding index A1, A2, or B. Similarly, fields evaluated in the spacetime points $x_{A1}^{\alpha }$, $x_B^{\alpha }$, or $x_{A2}^{\alpha }$ of the emission and reception events are denoted by the corresponding index, for example, $W_{A1}=W\left(x_{A1}^{\alpha }\right)$, ${\upsilon }_{\text{R}|A1}^i={\upsilon }_R^i\left(x_{A1}^{\alpha }\right)$

The goal of this section is to express the correction ∆ defined by Eq. (58). From this definition, it follows that $$\text{Δ}=\frac{{\nu }_{A2}}{{\nu }_B}\left(1-\frac{1}{2}\frac{{\nu }_B}{{\nu }_{A1}}\right)-\frac{1}{2}.$$(206)

Using Eq. (99) for the frequency shifts in the particular directions, we obtain $$\frac{{\nu }_{A1}}{{\nu }_B}=\frac{u_{A1}^0}{u_B^0}\frac{1-\frac{1}{c}{\upsilon }_{A1}^i{\overline{l}}_{i+|0}}{1-\frac{1}{c}{\upsilon }_B^i{\overline{l}}_{i+|1}}\exp I_{+},$$(207) $$\frac{{\nu }_B}{{\nu }_{A2}}=\frac{u_B^0}{u_{A2}^0}\frac{1-\frac{1}{c}{\upsilon }_B^i{\overline{l}}_{i-|0}}{1-\frac{1}{c}{\upsilon }_{A2}^i{\overline{l}}_{i-|1}}\exp I_{-},$$(208)

In order to proceed with the correction (206), it is useful to express the relations for the contributions to ${\overline{l}}_{i|\eta }$ given by Eqs. (197) and for the contributions to I given by Eqs. (205) between the two directions of the light-ray propagation. Transforming the integration variables in the contributions to ${\overline{l}}_{i-|\left(1-\eta \right)}$ and in the contributions to I₋ as $\widehat{l}=\widehat{L}-s$, $\widehat{l}\prime =\widehat{L}-s\prime$, and so on, using Eqs. (187) and (188), taking into account that ${\chi }_{+}^i=-{\chi }_{-}^i$, using the properties (153), (154) and taking into account that $\text{Δ}{x}_{F+}^i=0$, we obtain the following relations: $${\overline{l}}_{i+\left(c^0\right)|\eta }=-{\overline{l}}_{i-\left(c^0\right)|\left(1-\eta \right)}$$(209a) $${\overline{l}}_{i+\left(c^{-1}\right)|\eta }={\overline{l}}_{i-\left(c^{-1}\right)|\left(1-\eta \right)}$$(209b) $${\overline{l}}_{i+\left(c^{-2}\right)|\eta }=-{\overline{l}}_{i-\left(c^{-2}\right)|\left(1-\eta \right)}$$(209c) $${\overline{l}}_{i+\left(\text{Δ}{x}_F\right)|\eta }=0$$(209d)

and $$I_{+\left(c^{-1}\right)}=I_{-\left(c^{-1}\right)}$$(210a) $$I_{+\left(c^{-2}\right)}=-I_{-\left(c^{-2}\right)}$$(210b) $$I_{+\left(c^{-3}\right)}=I_{-\left(c^{-3}\right)}$$(210c) $$I_{+\left(\text{Δ}{x}_F\right)}=0.$$(210d)

The u⁰ terms in Eqs. (207) and (208) are expressed using Eq. (190). Quantities located at event A2 (except the potential W_A2) are expressed in terms of the quantities located at A1 using Eq. (186) and expressing the velocity difference $\text{Δ}{\upsilon }_A^i=\text{Δ}{\upsilon }_{A2}^i-\text{Δ}{\upsilon }_{A1}^i$ in terms of the acceleration $a_{A1}^i={\text{d}}^2{x}_A^i/{\text{d}{t}^2|}_{t_{A1}}$ of the observer A as $$\text{Δ}{\upsilon }_A^i=a_{A1}^i\frac{2\widehat{L}}{c}+O\left(2\right)$$(211)

Then, inserting Eqs. (207) and (208) into Eq. (206), we obtain the following result including all terms on the order of 𝒫(3) except the p_T ≥ 3 order terms originating from the motion of the observer on the ground A in the corotating frame and except the terms quadratic in ${\upsilon }_{A1}^i/c$: $$\begin{array}{lll}\text{Δ}\hfill & =\hfill & \frac{1}{2c^2}\left[W_{A1}+W_{A2}-2W_B+|{{\upsilon }_{R|A1}^i|}^2-{|{\upsilon }_{R|B}^i|}^2-{|{\upsilon }_B^i|}^2\right]\left(1-\frac{1}{c}{\upsilon }_B^j{\chi }_{j+}\right)\hfill \\ \hfill & \hfill & +\frac{1}{c}{\upsilon }_B^i{\beta }_i-\frac{1}{c^2}{\upsilon }_B^i{\beta }_{i(1)}{\upsilon }_B^j{\chi }_{j+}\hfill \\ \hfill & \hfill & +\frac{1}{c^2}D{\chi }_{i+}(a_{A1}^i+{ϵ}^i{}_{jk}{\omega }^j{\upsilon }_{A1}^k)+\frac{1}{c^2}{\upsilon }_{Bi}{\upsilon }_{A1}^i\hfill \\ \hfill & \hfill & +I_{+(c^{-2})}-\frac{1}{c^2}{\upsilon }_B^i{\chi }_{i+}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\partial }_{t}n\left({\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)\right),\hfill \end{array}$$(212)

where ${\beta }_i={\overline{l}}_{i+\left(c^{-1}\right)|1}-\frac{1}{c}{\upsilon }_{R|B}^i$ is given by Eq. (214) below, and β_i(1) is its first-order (p_T = 1) part given as $${\beta }_{i\left(1\right)}=-\frac{1}{c}{ϵ}_{ijk}{\omega }^j\left(x_B^k-x_{A1}^k\right)=-\frac{1}{c}\left({\upsilon }_{R|B}^i-{\upsilon }_{R|A1}^i\right),$$(213)

and $I_{+\left(c^{-2}\right)}$ is given by Eq. (215). $$\begin{array}{lll}{\beta }_i\hfill & =\hfill & \frac{1}{c}{A}_B^i\hfill \\ \hfill & \hfill & +\frac{1}{c}{n}_B\left({\delta }^{ij}-{\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}|}_{\widehat{L}}{\frac{\text{d}{\widehat{x}}_{+}^j}{\text{d}\widehat{l}}|}_{\widehat{L}}\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\frac{1}{n_{+}}{ϵ}_{jkl}\left(2{\omega }^l+\text{curl}{A}_{+}^l\right)\frac{\text{d}{\widehat{x}}_{+}^k}{\text{d}\widehat{l}}\hfill \\ \hfill & \hfill & +{\widehat{n}}_B\left({\delta }_j^i-{\delta }_{jq}{\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}|}_{\widehat{L}}{\frac{\text{d}{\widehat{x}}_{+}^q}{\text{d}\widehat{l}}|}_{\widehat{L}}\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\{\left(\widehat{l}{Q}_{k+\left({\widehat{N}}^p\right)}^j\left(\widehat{l}\right)-C_{k+\left({\widehat{N}}^p\right)}^j\left(\widehat{l}\right)\right)\hfill \\ \hfill & \hfill & \times \left[{\displaystyle \sum _{n=0}^{\infty }{\left({\mathcal{𝒦}}_{l+\left({\widehat{N}}^p\right)}^k\right)}^n}{\mathcal{𝒟}}_m^l{S}_{+\left({\widehat{N}}^p{c}^{-1}\right)}^m\right]\left(\widehat{l}\right)\}\hfill \\ \hfill & \hfill & +\frac{1}{c}\left({\delta }^{ij}-{\chi }_{+}^i{\chi }_{+}^j\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_t{\partial }_j{\alpha }_{+}}\hfill \\ \hfill & \hfill & +\frac{1}{c}{ϵ}_{klm}{\omega }^m{\chi }_{+}^l\left({\delta }^{ia}{\delta }^{kb}+{\chi }_{+}^a{\chi }_{+}^b{\delta }^{ik}\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\widehat{l}}^2\left(\widehat{l}-\widehat{L}\right){\partial }_a{\partial }_b{\alpha }_{+}}\hfill \\ \hfill & \hfill & +\frac{1}{c}{ϵ}_{mjk}{\omega }^j\left({\delta }^{im}{\delta }^{kl}-2{\chi }_{+}^k\left({\delta }^{im}{\chi }_{+}^l+{\delta }^{lm}{\chi }_{+}^i\right)\right)\frac{1}{\widehat{L}}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_l{\alpha }_{+}\hfill \end{array}$$(214) $$\begin{array}{lll}{I}_{+\left(c^{-2}\right)}\hfill & =\hfill & \frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\left({\partial }_t{A}_{i+}\right)\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}}\hfill \\ \hfill & \hfill & +\frac{1}{c^2}{ϵ}^i{}_{jk}{\chi }_{+}^j{\omega }^k{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\partial }_i{\partial }_t{\alpha }_{+}\hfill \\ \hfill & \hfill & +\frac{1}{c^2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left(\widehat{l}-\widehat{L}\right){\partial }_i{\partial }_t{\alpha }_{+}\hfill \end{array}$$(215)

In Eqs. (214) and (215), the index + for fields denotes that the field is evaluated along the ${\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)$ path, for example, ${\alpha }_{+}=\alpha \left({\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)\right)$, ${\partial }_k{\alpha }_{+}=\left({\partial }_k\alpha \right)\left({\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)\right)$. The functions $Q_{k+\left({\widehat{N}}^p\right)}^j$, $C_{k+\left({\widehat{N}}^p\right)}^j$, $S_{+\left({\widehat{N}}^p{c}^{-1}\right)}^m$ and the generating function $K_{l+\left({\widehat{N}}^p\right)}^k$ for the operator ${𝒦}_{l+\left({\widehat{N}}^p\right)}^k$ in the third and fourth line of Eq. (214) are given by Eqs. (198)–(201) with ${\widehat{x}}^{\alpha }\left(\widehat{l}\right)={\widehat{x}}_{+}^{\alpha }\left(\widehat{l}\right)$.

To quantify the effects in the examples in this paper and for most of the satellite applications, it is practical to express Eq. (214) explicitly up to the second order (p_T ≤ 2), leading to the contribution up to the third order (p_T ≤ 3) in ∆. Moreover, the terms in Eq. (214) that are proportional to c⁻¹ N_B are negligible for satellite applications because the atmospheric refractivity N_B in the satellite position is very low for satellite altitudes of several hundred kilometers. Thus, up to the second order and neglecting the terms proportional to c⁻¹ N_B, we can write $${\beta }_i={\beta }_{i\left(1\right)}+{\beta }_{i\left(2\right)},$$(216)

with β_i(1) given by Eq. (213) and β_i(2) given by Eq. (60). Applying Eq. (216) in Eq. (212) leads to Eq. (59) in the summary section.

Some steps in the calculation of β_i up to the second order are described below. For the end-point tangent, we use Eq. (40) expanded up to the first order, $${\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}|}_{\widehat{L}}={\chi }_{+}^i+{\varsigma }_{+}^i{\epsilon }_B+O\left(2\right),$$(217)

where the index B replaces the general index for a final point F. Where possible, the approximation (41) is used. After expanding $1/n_{+}=1-{\widehat{N}}_{+}-{\alpha }_{+}+O\left(2\right)$ in the integral term in the second line of Eq. (214), we use the following way of integrating the term containing ω^l: $$\begin{array}{lll}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}\hfill & =\hfill & \frac{1}{2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left[\frac{\text{d}}{\text{d}\widehat{l}}\left(\widehat{l}\left(\widehat{l}-\widehat{L}\right)\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}\right)+\widehat{L}\frac{\text{d}{\widehat{x}}_{+}^i}{\text{d}\widehat{l}}-\widehat{l}\left(\widehat{l}-\widehat{L}\right)\frac{{\text{d}}^2{\widehat{x}}_{+}^i}{\text{d}{\widehat{l}}^2}\right]\hfill \\ \hfill & =\hfill & \frac{\widehat{L}}{2}\left(x_B^i-x_{A1}^i\right)\hfill \\ \hfill & \hfill & -\frac{1}{2}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\left(\widehat{l}-\widehat{L}\right)\left({\delta }^{ij}-{\chi }_{+}^i{\chi }_{+}^j\right){\partial }_j{\widehat{n}}_{+}+O\left(2\right),\hfill \end{array}$$(218)

where the second derivative of the path ${\widehat{x}}_{+}^i\left(\widehat{l}\right)$ was expressed using Eq. (132) expanded up to the p_T = 1 order.

Expanding the term in the third and fourth line of Eq. (214) with the use of Eqs. (157) and (198)–(201), one of the terms coming from $Q_{k+\left({\widehat{N}}^p\right)}^j$ can be processed as follows: $$\begin{array}{l}{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\widehat{l}}^2\left(\widehat{l}-\widehat{L}\right){\chi }_{+}^i{\chi }_{+}^j{\partial }_i{\partial }_j{\widehat{n}}_{+}\hfill \\ ={\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}{\widehat{l}}^2\left(\widehat{l}-\widehat{L}\right)\frac{\text{d}}{\text{d}\widehat{l}}\left({\chi }_{+}^j{\partial }_j{\widehat{n}}_{+}\right)+O\left(2\right)\hfill \\ =-{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\left(2\widehat{l}\left(\widehat{l}-\widehat{L}\right)+{\widehat{l}}^2\right){\chi }_{+}^j{\partial }_j{\widehat{n}}_{+}+O\left(2\right)\hfill \\ =-2{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\chi }_{+}^j{\partial }_j{\widehat{n}}_{+}-{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}{\widehat{l}}^2}\frac{\text{d}}{\text{d}\widehat{l}}{\widehat{N}}_{+}+O\left(2\right)\hfill \\ =-2{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}}\widehat{l}\left(\widehat{l}-\widehat{L}\right){\chi }_{+}^j{\partial }_j{\widehat{n}}_{+}+2{\displaystyle {\int }_0^{\widehat{L}}\text{d}\widehat{l}\widehat{l}{\widehat{N}}_{+}}-{\widehat{L}}^2{\widehat{N}}_B+O\left(2\right),\hfill \end{array}$$(219)

where in the second and fourth line, the approximation (41) was used and the terms were integrated by parts. The same procedure can be applied to the analogous term in Eq. (214) with α₊ instead of ${\widehat{n}}_{+}$. The use of these relations leads to the cancellation of several terms in the expression for β_i(2).

6 Conclusion

A relativistic model of one-way and two-way time and frequency transfer through flowing media was developed in this paper. The main focus was on applications to the atmosphere of Earth. The model includes gravitational and atmospheric effects given by the fields of the scalar gravitational potential W(x^α), the atmospheric refractive index n(x^α), and the wind speed Vⁱ (x^α). The nonstationarity of the fields and deviations from spherical symmetry are also included in the model.

The method used in this paper is based on solving the equation of motion for a light ray in coordinates x^α = (ct, xⁱ) corotating with the Earth. This equation is given as the null geodesic equation of the Gordon optical metric (Gordon 1923). First, the exact solution ${\widehat{x}}^i\left(\widehat{l}\right)$ in a spherically symmetric, static part of n(x^α) and W(x^α) is found, which is parameterized by its Euclidean length in the corotating frame l̂. This solution is defined by the spatial coordinates of the emission and the reception events as boundary values. Then, a complete solution $x^i\left(\widehat{l}\right)$, $t\left(\widehat{l}\right)$ is found using a perturbation method that takes the inertial forces, wind speed, and deviations of n(x^α) and W(x^α) from sphericity and staticity into account. The time- and frequency-transfer corrections are then expressed in terms of integrals of the fields along the path ${\widehat{x}}^i\left(\widehat{l}\right)$ at a certain constant coordinate time.

In the numerical examples evaluated in this paper, we focused on the two-way ground-to-satellite transfer with a satellite altitude similar to that of the ISS. In the two-way time transfer, the main contribution of the atmosphere appears in the Sagnac effect and is given by a change in the Sagnac area that is spanned by the light path ${\widehat{x}}^i\left(\widehat{l}\right)$ when the path changes from the vacuum case to the case influenced by the atmospheric refraction. In the example studied in this paper, the contribution of the atmosphere to the Sagnac effect is below 0.1 ps for the satellite zenith angle θ (the angle between the vertical line from a ground station and the line connecting the ground station with the satellite position at the reception/re-emission time) below 78°, it reaches 1 ps for θ ≈ 87°, and it increases to approximately 5 ps as θ approaches 90°. Therefore, for time transfer at an accuracy level of 1 ps, this effect is only significant for large zenith angles.

The effect of the wind in the two-way time transfer is given by the difference in the propagation times caused by the Fresnel-Fizeau effect of light dragging. It is much smaller than 1 ps for normal atmospheric conditions.

In the two-way frequency transfer, the correction Δ of Ashby (1998) and Blanchet et al. (2001) is generalized to include the atmospheric effects. In the example presented in this paper, the effect of the spherically symmetric, static part of the refractive index gives a contribution that starts at 10⁻¹⁷ for θ = 0, reaches 10⁻¹⁶ for θ ≈ 56°, and ends at 10⁻¹³ for θ approaching 90°.

Remarkably, in the equation of motion for a light ray, the quantity Aⁱ = (1 − n²)Vⁱ related to the wind speed acts similarly as the magnetic potential in the equation of motion for a charged particle. In particular, the opposite direction of the propagation leads to opposite acceleration. In the two-way frequency transfer, it may lead to a significant contribution to the correction Δ. For the example of a constant horizontal wind field, the effect reaches 10⁻¹⁷ already for a wind speed of about 11 m s⁻¹ (the wind speed in the corotating frame is meant).

The effects of the deviation of n(x^α) from sphericity and staticity are also included in the resulting corrections. They are not evaluated numerically in this paper, however.

The evaluation of the atmospheric effects in the time- and frequency-transfer corrections including the effect of the wind shows that they need to be considered in the data analysis of the forthcoming clock-on-satellite experiments such as ACES or I-SOC. The results obtained in this paper can be used not only to calculate the time- and frequency-transfer corrections from given atmospheric data, but also inversely to determine atmospheric properties such as the refractive index profile or the flow rate from the time or frequency measurements. This will be a subject of a future research.

Acknowledgements

This research was funded by grants of Institutional Financing IF2105021501, IF2205021501, and IF2309021501 provided by Ministry of Industry and Trade of the Czech Republic to the Czech Metrology Institute. The author is also thankful to Adrien Bourgoin and Pacôme Delva for their valuable comments about the paper manuscript.

Appendix A Refractive index profiles

A.1 Refractive index as a function of the density of air

According to Eq. (4) of Owens (1967) the refractive index of air n depends on the densities of its components as $$\frac{n^2-1}{n^2+2}=R_1{\rho }_1+R_2{\rho }_2+R_3{\rho }_3,$$(A.1)

where ρ₁, ρ₂, and ρ₃ are the partial densities of dry, CO₂ -free air, of water vapor, and of carbon dioxide in the mixture, and R₁, R₂, and R₃ are the corresponding specific refractions of the components of the mixture that only depend on wavelength. Therefore, to obtain the field n(x^α) in spacetime, the fields ρ₁(x^α), ρ₂(x^α), and ρ₃(x^α) need to be known.

An approximate formula for the refractivity N ≡ n − 1 can be derived from Eq. (A.1). This formula is given, for example, by Eq. (5) of Ciddor (1996) which states $$N=\frac{{\rho }_a}{{\rho }_{axs}}{N}_{axs}+\frac{{\rho }_w}{{\rho }_{ws}}{N}_{ws},$$(A.2)

where ρ_axs is the density of pure dry air at 15°C and 101 325 Pa with certain molar fraction x_c of CO₂, ρ_ws is the density of pure water vapor at 20 °C and 1333 Pa and N_axs, and N_ws are the refractivities of the dry air and water vapor at the conditions corresponding to ρ_axs, ρ_ws. The refractivities as functions of the vacuum wavelength are known for these specific conditions, and they are given by Ciddor (1996). Similarly, ρ_a and ρ_w are the densities of the dry air component and of the water vapor component of moist air for the actual conditions. The densities ρ_a and ρ_w can be derived from the equation of state (Eq. (4) of Ciddor (1996)) as functions of the temperature T, the pressure p, and the air composition, represented, for example, by the molar fractions x_c of CO₂ and x_w of the water vapor. Therefore, the refractivity field N(x^α) for a given fixed vacuum wavelength can be determined when the fields T(x^α), p(x^α), x_c(x^α), and x_w(x^α) are known.

For the special case when the air composition does not change with spacetime position (i.e., x_c and ρ_w/ρ_α are constants in spacetime), we can divide Eq. (A.2) by the total air density ρ = ρ_a + ρ_w, and we obtain that N/ρ is constant in spacetime. In this case, we can write $$N(x^{\alpha })=N_0\frac{\rho (x^{\alpha })}{{\rho }_0},$$(A.3)

where N₀ and ρ₀ are the air refractivity and the density in a certain reference spacetime point.

A.2 Air density field of the atmosphere in hydrostatic equilibrium

The temperature of the atmosphere of Earth as a function of altitude and the related pressure and density fields derived under the assumption of hydrostatic equilibrium up to altitudes of about 80 km are discussed, for example, in ICAO (1993) and NOAA (1976) and we summarize them in this section.

Assuming that the atmosphere rigidly rotates together with the gravitating body and is in hydrostatic equilibrium with its gravitational and centrifugal forces, the atmospheric pressure p satisfies $${\partial }_{i}p=\rho {\partial }_i\varphi ,$$(A.4)

with $$\varphi =W+{\upsilon }_R^2/2$$(A.5)

being the gravitational plus centrifugal potential. Next, we assume the perfect gas equation of state for the air, $$p=\frac{RT}{M_a}\rho ,$$(A.6)

with R being the universal gas constant, T thermodynamic temperature, and M_a the molar mass of the air, which is assumed to be constant. Following ICAO (1993) we define the geopotential altitude $$H=-\left(\varphi -{\varphi }_s\right)/g,$$(A.7)

where ϕ_s is the sea-level value of ϕ, and g is the standard acceleration due to gravity (constant). We note that H has the dimension of length. We assume that the atmosphere can be divided into intervals in H in which the temperature depends on H approximately linearly as $$T=T_b+\beta \left(H-H_b\right),$$(A.8)

where H_b is the boundary value of H between two intervals, T_b is the temperature at this boundary, and β is the constant that is different for each of the intervals (see Table D of ICAO (1993)). The solution of the hydrostatic equilibrium equation (A.4) for ρ using Eqs. (A.6) and (A.8) with β ≠ 0 reads $$\rho ={\rho }_b{\left(1+\frac{\beta }{T_b}\left(H-H_b\right)\right)}^{-\left(1+\frac{M_{a}g}{\beta R}\right)},$$(A.9)

and for β = 0, it reads $$\rho ={\rho }_b\exp \left(-\frac{M_{a}g}{R{T}_b}\left(H-H_b\right)\right).$$(A.10)

A.3 Spherically symmetric, static component of the atmospheric refractivity

The time- and frequency-transfer model developed in this paper uses a spherically symmetric, static component of the refractive index n(r) (see Eq. (6)) or the related refractivity $\widehat{N}\left(r\right)=\widehat{n}\left(r\right)-1$. One option how to define this field is to use Eqs. (A.9) and (A.10) together with Eq. (A.3) and to replace the potential ϕ given by Eq. (A.5) by its monopole part Ŵ = GM/r.

In Sects. 4.5 and 4.7, where we evaluated the magnitudes of the effects for specific examples, we assumed for sake of simplicity an isothermal atmosphere with a constant temperature T₀. From Eqs. (A.10) and (A.3), we then obtain $$\rho ={\rho }_0\exp \left(\frac{M_a}{R{T}_0}\left(\varphi -{\varphi }_0\right)\right),$$(A.11) $$N=N_0\exp \left(\frac{M_a}{R{T}_0}\left(\varphi -{\varphi }_0\right)\right),$$(A.12)

where ρ₀, N₀, and ϕ₀ are the air density, air refractivity, and the potential in a certain reference spacetime point. The refractivity $\widehat{N}\left(r\right)$ is then given as $$\widehat{N}=N_0\exp \left(\frac{M_a}{R{T}_0}\left(\widehat{W}-{\widehat{W}}_0\right)\right),$$(A.13)

where Ŵ₀ is the value of Ŵ in the reference point.

References

All Figures