Issue 
A&A
Volume 673, May 2023



Article Number  A144  
Number of page(s)  25  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/202345994  
Published online  23 May 2023 
Relativistic theory for time and frequency transfer through flowing media with an application to the atmosphere of Earth
Czech Metrology Institute,
Okružní 31,
63800
Brno, Czech Republic
email: jgersl@cmi.cz
Received:
25
January
2023
Accepted:
22
March
2023
Context. Several space missions that will use atomic clocks on board of an Earthorbiting satellite are planned for the near future, such as the Atomic Clock Ensemble in Space (ACES) or the Space Optical Clock on the International Space Station (ISOC). The increasing accuracies of the developed clocks and of the links connecting them with ground stations impose corresponding accuracy requirements for theoretical models of electromagnetic signal propagation through the atmosphere of Earth and for the related time and frequency transfer corrections. For example, the fractional frequency accuracy of the optical lattice clock for the ISOC project is about 10^{−17}.
Aims. We develop a relativistic model of one and twoway time and frequency transfer. In addition to the gravitational effects, it also includes the effects of atmospheric refractivity and atmospheric flows within the relativistic framework.
Methods. The model is based on an analytical solution of the equation of motion of a light ray in spacetime filled with a medium: the null geodesic equation of Gordon’s optical metric.
Results. Explicit formulas for one and twoway time and frequency transfer corrections are given using realistic fields of the gravitational potential, the refractive index, and the wind speed, taking nonstationarity and deviations from spherical symmetry into account. Numerical examples are provided that focus on twoway groundtosatellite transfer, with satellite parameters similar to those of the International Space Station. The effect of the atmospheric refractive index increases as the satellite position moves from zenith to horizon, and it is shown that the effect ranges from 0 ps to 5 ps for twoway time transfer and from 10^{−17} to 10^{−13} for twoway frequency transfer, with a steep increase as the satellite approaches the horizon. The effect of the wind contribution is well below 1 ps for the twoway time transfer for normal atmospheric conditions, but for the twoway frequency transfer, the effect can be significant: A contribution of 10^{−17} is possible for a horizontal wind field with a velocity magnitude of about 11 m s^{−1}.
Conclusions. The atmospheric effects including the effect of wind should be considered in the forthcoming clockonsatellite experiments such as ACES or ISOC.
Key words: relativistic processes / atmospheric effects / time / reference systems / methods: analytical
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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^{−18} in case of stationary clocks (Ludlow et al. 2015; McGrew et al. 2018; Lisdat et al. 2021) or 10^{−17} in case of transportable clocks (Koller et al. 2017; Cao et al. 2017; Origlia et al. 2018), together with developments in freespace optical links (see Bodine et al. 2020 and references therein, and also Djerroud et al. 2010, Gozzard et al. 2018, Kang et al. 2019, DixMatthews 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 groundtosatellite 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 Hmaser (Cacciapuoti et al. 2020; Lilley et al. 2021), and the Space Optical Clock on the ISS (ISOC), with an optical lattice clock (Cacciapuoti & Schiller 2017; Origlia et al. 2018). In the ACES mission, the absolute frequency accuracy of the onboard clock is about 10^{−16} in fractional frequency, whereas the same parameter for the ISOC mission is about 10^{−17}. 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 ISOC 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 twoway 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^{−18} for the frequency transfer.
Blanchet et al. (2001) developed the relativistic theory of one and twoway time and frequency transfer in vacuum, including terms up to the order of c^{−3} to cover contributions relevant for experiments with an uncertainty of 5 × 10^{−17} 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énardWiechert representation of the metric tensor up to the first postMinkowskian 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^{−4} terms in frequency shift approach the 10^{−18} value for frequency transfer from a ground station to the ISS. They also discussed the influence of the quadrupole moment of the Earth J_{2} at the c^{−3} order. The theory of Linet & Teyssandier (2002) was further developed by Le PoncinLafitte et al. (2004) and Teyssandier & Le PoncinLafitte (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 PoncinLafitte (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 lightdragging effect due to the steady rotation of the atmosphere of Earth was discussed. The models of Linet & Teyssandier (2002) and Bourgoin (2020) focused on oneway 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 postMinkowskian 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 FresnelFizeau 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). Propagationtime nonreciprocity in twoway groundtosatellite 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 twoway 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 twoway 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^{−18} 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 twoway 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^{i} is denoted as , with δ_{ĳ} being the Kronecker delta. Quantities ω^{i}, , υ^{i}, V^{i}, A^{i}, h^{i} ,χ^{i}, ς^{i}, and γ^{i} in this paper are treated as Euclidean threevectors when their index is lowered or raised, for example, V_{i} = δ_{ĳ}V^{j} The threedimensional antisymmetric LeviCivita symbol is denoted as ϵ_{ĳk}, eventually .
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 , where is the Geocentric Coordinate Time (TCG), and 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) (1a) (1b) (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 twoway corrections more straightforward.
The rotating coordinate system x^{α} = (x^{0}, x^{i}) is related to the GCRS system as , , where R^{i}_{j}(t) ∊ SO(3), or inversely, , where . The timedependent rotation matrix R^{i}_{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 (2)
where are constant components of the angular velocity vector of the coordinate rotation expressed in the frame. Since 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^{i}. We denote 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 with .
The components of the metric (1) in the rotating coordinates are (3a) (3b) (3c)
where . The monopole part of the scalar gravitational potential is denoted as , and the remaining part containing the higher multipole moments and the timedependent components as ∆W(x^{α}) such that (4) (5)
where r = ∣x^{i}∣ 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 nondispersive 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 (6)
where 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 , α(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 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 largescale 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 are then defined as deviations of the corresponding refractive index from its vacuum value, (7) (8)
Denoting the trajectory of the medium element in the rotating coordinates that passes through a spacetime point , the winds in the atmosphere are described by their velocity field V^{i}(x^{α}), which is defined by (9)
for any where the medium is present. We also define the following quantity that is useful when the optical effects of the winds are evaluated: (10)
3.3 Observers
For a given trajectory x^{i} (t) of an observer in the rotating coordinates, we denote υ^{i} = dx^{i}/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 twoway 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^{−18} for the frequency transfer is needed in order to meet the requirements of the forthcoming clockonsatellite 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^{2} (higherorder gravitational effects such as the LenseThirring effect are negligible considering our intended accuracy), from inertial effects given by the rotation velocity , from atmospheric effects given by the air refractivity together with the correction α and by the wind speed V^{i} /c and from the observers’ velocities υ^{i}/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 and its correction α(x^{α}), (b) we denote the maximum value of and α_{M} the maximum value of ∣α∣, (c) we define normalized functions and α_{*} = α/α_{M} and in all formulas we express and α = α_{M}α_{*} (d) the order of a term in an expansion is then given by its factor , where the powers , , and p_{α} are nonnegative integers, and finally, (e) the derivative of the fields ∂_{0} = c^{−1} ∂_{t} increases the power of c^{−1} 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 . The expansion of n in terms of , α_{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 (i.e., p_{α} = 0) in the time and frequency transfer formulas, we take corresponding to the vacuum model of Blanchet et al. (2001), and we take an arbitrary value of for . Terms on the order of with are neglected. Some of the terms on the order can slightly exceed the 10^{−18} level in case of groundtoISS frequency transfer with a near horizon position of the satellite. The c^{−4} terms, which we neglect as well, can approach the 10^{−18} level for groundtoISS 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 . We therefore only include the field α to a lower order. Namely, for p_{α} ≥ 1, we take . Possible magnitudes of the neglected terms with p_{α} ≥ 1 have not been evaluated, however.
Next, we denote the total , 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 twoway time and frequency transfer corrections in this paper include all terms with the order from a set 𝒫(3), which is given as
The only exception in which we express the results up to a lower order only are the terms in the twoway 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 and only are denoted O(c^{−k}) as usual. It is also useful to define the following subset of the exponent triplets:
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 , with 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 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 firstorder differential equation for the coordinate time is solved.
The main results of this paper are expressed in terms of the integrals along the path . 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 and gravitational potential Ŵ(r) is given as (see Eq. (132)) (11)
where the fields , , and ∂jŴ are evaluated in , and the solution has a unit tangent satisfying (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 twoway time and frequency transfer setups. We denote the total path length between the initial and final points as , such that the parameter has the range .
which depends on the spatial coordinates through the radius r = ∣x^{i}∣ only and introducing a new parameter ζ with a function dependence satisfying (14)
Eq. (11) can be written as (15)
where we briefly denote . 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 (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 and using the assumption of the central field, which leads to (17)
with h^{i} being a constant vector. Equation (17) shows that for ∣h^{i}∣ ≠ 0, the signal propagates in a plane that intersects the origin r = 0 and has a normal aligned with h^{i} (for ∣h^{i}∣ = 0, the signal path is just a straight radial line). Denoting the angle between the tangent and the position vector and taking the magnitude of Eq. (17) with use of Eq. (16), we obtain (18)
where and h = ∣h^{i}∣. 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 (assuming all along the trajectory), Eq. (18) can be written as (19)
Selecting a righthanded Cartesian coordinate system y^{i} with the origin at r = 0, the y^{1}axis in direction, and the y^{3}axis in the h^{i} direction (for h = 0, we can take any direction perpendicular to ), defining polar coordinates by y^{1} = r cos φ and y^{2} = r sin φ and expressing the trajectory in the polar coordinates as φ(r), with the radius as the parameter, it follows from the definition of ψ that (20)
Integrating this equation with tan ψ expressed from the Snell law (19), we obtain (21)
with (ascending trajectory) and ± = − for (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 .
In the following, we denote quantities related to the initial point by the index I and to the final point by the index F, for example, , , , ψ_{I} = ψ(0), and , , , and . We also define the angle between and , the angle θ_{I} between and , and the angle θ_{F} between and , and we denote (see Fig. 1). The following relations hold: (22) (23)
The boundary value problem is more complicated than the initial value problem, in which case h is given simply as h = r_{I}w_{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 , θ_{I}, or θ_{F}. One possible approach is presented below. We denote as ε_{I} and ε_{F} the angular deviations of the path tangents and , respectively, from the vector (see Fig. 1), so that we have (24) (25)
Evaluating the Snell law (19) in the boundary points and using Eqs. (24) and (25), we obtain (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 and , which satisfies (27)
Solving Eq. (26) together with Eq. (27) for the angles ε_{I} and ε_{F}, we obtain (28) (29)
The total bending angle also satisfies (30)
where the angle can be expressed as using Eq. (21), and we can also write (31)
where on the righthand side (RHS), we use Eq. (19). Thus, we obtain (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_{F}w_{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 groundtosatellite or satellitetoground 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^{i}, we express the components of the first two basis vectors of y^{i} in the coordinate system x^{i}. Assuming , we have (33)
The solution parameterized by r is then given as (34)
For , it is simply .
The reparameterization function can be expressed with the use of (35)
where the ± sign has the same meaning as in Eq. (21). Integrating this formula using the Snell law (19), we obtain (36)
For we then get .
It is also useful to define an orthonormal righthanded basis (χ^{i}, ς^{i}, γ^{i}) with χ^{i} pointing in direction from to , γ^{i} pointing in direction of the normal to the propagation plane h^{i} (assuming h ≠ 0), and ς^{i} completing the triplet. We have (see also Fig. 1) (37) (38) (39)
Defining the angle between the path tangent and the vector χ^{i} satisfying ε(0) = ε_{I} (negative value) and (positive value), we can write (40)
In some formulas, the lowestorder approximation of the tangent and of the length is sufficient. Expanding Eq. (40) in ε and using ∣ε∣ ≤ ε_{T} = O(1), which follows from Eq. (32), we obtain (41)
where the expansion of cos ε was used in the last step.
In the following text, we refer to the solution of Eq. (11) with the given boundary conditions as the base path.
Fig. 1 Definition of angles, basis vectors, and other quantities. 
4.4 Coordinate time transfer
4.4.1 Oneway time transfer
In oneway time transfer, we wish to express a coordinate time t_{F} − t_{I} of the propagation of the signal emitted from a position at a time t_{I} and received in a position at a time t_{F}. We consider a base path with boundary points and .
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): (43)
where the fields n, α, W, υ_{Ri}, A_{i} and their partial derivatives with respect to the spacetime coordinates are evaluated in t = t_{I} and if not stated explicitly otherwise, the velocities and are defined as and , respectively, and the function 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^{i} in A^{i} 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 FresnelFizeau 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 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 Twoway time transfer
In twoway 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 . 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 . The pseudotimeofflight τ_{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 . This signal can either be the same signal reflected or another signal that is synchronously sent when the oneway signal is received by observer B. We call this setup the Λconfiguration. The pseudotimeofflight τ_{A}(t_{A2}) − τ_{B}(t_{B}) is also obtained from measurements.
Using the coordinatetime 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 twoway time transfer, it can be expressed as (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 pseudotimesofflight, 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 (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 and is very small, and we include its effect up to the largest order only using the velocity 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 the base path with the boundary points and , we arrive at the following result for the twoway time transfer correction, including all terms on the order of 𝒫(3) except the p_{T} = 3 order terms originating from (46)
where the fields υ_{Ri}, A_{i}, ∂_{i}α, and ∂_{t}α with the index + are evaluated in t = t_{B} and , and is the vector χ^{i} corresponding to the path (i.e., it is the unit vector in the direction from to ).
The main contribution to the twoway time transfer correction in Eq. (46) is given by the first term in the first line, which is the wellknown Sagnac effect. We denote (47)
Then, taking into account that , the correction given by Eq. (47) can be expressed as (48)
where Σ_{j} is defined as (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 , 2) the path from to , and 3) the straight line from 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 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 (50)
where is the vector γ_{j} corresponding to the path , and w(r) is given by Eq. (13). This integral can be computed numerically for any refractive index profile 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.
Fig. 2 Scheme of the twoway time transfer. The signal leaves observer A at a coordinate time t_{A1}, corresponding to the proper time τ_{A}(t_{A1}) of clock A, is reflected from observer B at a coordinate time t_{B}, corresponding to the proper time τ_{B}(t_{B}) of clock B and τ_{A}(t_{B}) of clock A, and is finally received by observer A at the coordinate time t_{A2}, corresponding to the proper time τ_{A}(t_{A2}) of clock A. 
Fig. 3 Visualization of the Sagnac area Σ. The part of the Sagnac effect that is caused by the atmosphere is given by the change in the area Σ when the path changes from the case including atmospheric refraction and gravitation to a vacuum case that only includes gravitation ( in Eq. (13)). 
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 groundtosatellite transfer through an isothermal atmosphere at a temperature T_{0} 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 is given by Eq. (A.13), which gives the following spherical refractive index field: (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 , which also equals (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_{0} = 288.15 K (15 °C) and M_{a} = 28.964 × 10^{−3} kg mol^{−1}, 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 and a light signal with a vacuum wavelength of 1 µm, we obtain (see, e.g., Ciddor 1996) (52)
In twoway transfer, we denote the θ_{I} angle of the path simply as θ. It represents the zenith angle of the satellite at the reception/reemission 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}γ_{i+} = −ω, which corresponds to the paths 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 is much smaller than the effect of atmospheric refraction, the path 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, (53)
Fig. 4 Sagnac correction, Eq. (48), for an isothermal atmosphere and its deviation from a vacuum value denoted as effect of the atmosphere as functions of the satellite zenith angle θ. Separate plots for θ ∊ [0°, 80º] (left) and θ ∊ [80°, 90º] (right) are presented to show details of the effect of the atmosphere curve, which grows steeply for θ approaching 90°. The example was evaluated for groundtosatellite transfer in the equatorial plane, with r_{A} = 6371 km, r_{B} − r_{A} = 408 km, and inclined against the rotation of the Earth. 
4.5.2 Effect of the wind
We evaluate the effect of the wind speed in the twoway time transfer correction (46), which is given by (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 , as well as the component of the air velocity V^{i} (in the corotating frame) in the direction of the path tangent , 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 (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 m^{2} s^{−1}, which, for example, for the path length of km leads to V ≈ 820 ms^{−1}. 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 interferometric 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 Oneway frequency transfer
In oneway frequency transfer, the goal is to express the ratio of the proper frequency ν¡ of a signal emitted by an observer with velocity from a position at a time t_{I} and a proper frequency ν_{F} of the same signal as is received by the observer with a velocity in the position at the time t_{F}. Again, we consider a base path with the boundary points and , and we denote , , , and
Based on the analysis included in Sect. 5, we obtain the following formula for the frequency ratio: (56)
where and are quantities related to tangents of the lightray trajectory at the emission and reception events, and they are given by Eqs. (196)–(197d) with η = 0 for and η = 1 for . In these formulas, we set and corresponding to . The details of the expressions for and are given in the Theory section, but the corresponding contribution to twoway transfer is discussed further in this summary. The exponent I in the last term of Eq. (56) is given as (57)
where the fields n, α, A_{i}, W, and their partial derivatives with respect to the spacetime coordinates are evaluated in t = t_{I} and 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 secondorder Doppler effect, the second fraction contains the firstorder Doppler effect, and the exponential term is the effect of nonstationarity.
4.6.2 Twoway frequency transfer
In twoway frequency transfer, a stationary observer on the ground A emits a signal with the proper frequency ν_{A1} from the position at the time t_{A1}. The signal is received by observer B (ground or satellite) with the proper frequency ν_{B} in the position 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 at the time t_{A2}. The velocities of observers A and B at the particular emission and reception events are denoted as , , and . For the stationary (fixed at the ground) observer A, the velocities and in the corotating coordinates and the corresponding spatial shift 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: (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 , 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^{−2}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 , , or of the emission and reception events are denoted by the corresponding index, for example, and ). The base path is again defined by the boundary points and , and the fields denoted with the index + are evaluated in t = t_{B} and , for example, . The basis; , , is associated with the path , and ε_{B} is the angle ε_{F} of the path . The acceleration of observer A at the time t_{A1} in the rotating frame is denoted (i.e., ). Using this notation, we obtain (59)
where the atmospheric effects are included in the last line, in which (60)
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 (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 to the twoway frequency transfer correction (59), which is given by the first three lines of Eq. (60). We denote (62a) (62b) (62c) (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 corotating system is υ_{B} = 7.17 km s^{−1}. 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 to and their sum 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 and were calculated by transforming the integration variable into r using Eq. (36).
Figure 5 shows that the sum of the corrections starts at 10^{−17} for θ = 0, reaches 10^{−16} for θ ≈ 56°, and ends at 10^{−13} for θ approaching 90°. Therefore, the corrections should be taken into account in the forthcoming clockonsatellite 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 (63)
To show how large the effect can be, we evaluate Eq. (63) for the simple situation with the satellite position radially above the emission point on the ground , in which case, the path is just a straight vertical line connecting with (i.e., and ). Next, we assume that the wind field V_{j} at t = t_{B} is perpendicular to the path tangent in the neighborhood of (horizontal flow). In this case, using Eq. (10) and writing (64)
the first term on the RHS of Eq. (64) vanishes, and the second term, which equals can be integrated by parts to obtain (65)
where the term containing , which is on the order of c^{−2}N_{B}, was neglected in the integration by parts. Expansion of Eq. (10) in the refractivity N gives (66)
For the numerical evaluation of Eq. (65), we further assume that in the neighborhood of the base path , the wind speed V^{i} is constant and has a direction opposite to the satellite velocity , which is itself perpendicular to (i.e., with and V = ∣V^{i}∣). 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 (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^{−1}, which approximately corresponds to the ISS, Eq. (65) gives (68)
This means that a horizontal wind speed V > 11 m s^{−1} can cause an effect that exceeds the 10^{−17} level.
We also note that in case of negligible air composition variations, Eq. (A.3) holds, and using Eq. (66), we obtain (69)
with ρ being the air density field and N_{0}, ρ_{0} 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 (70)
with n^{i} 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 in the direction perpendicular to both and by the small Euclidean distance d. In this case, taking n^{i} in the direction of Eq. (65) leads to (71)
However, this relation between the correction 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).
Fig. 5 Three contributions of the spherically symmetric, static part of the atmospheric refractivity , and to the twoway frequency transfer correction (59) and their sum as functions of the satellite zenith angle θ. The example is evaluated for groundtosatellite transfer through an isothermal atmosphere in hydrostatic equilibrium with 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_{A} + 408 km in the direction of the rotation of the Earth. 
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 fourvelocity 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 .
According to Gordon (1923), the Maxwell equations in the medium can be reformulated using the optical metric, which is defined as (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 the covariant derivative associated with the optical metric . We introduce a notation with the bar for tensors, obtained by lowering or raising indices with the optical metric or its inverse . For example, for the electromagnetic tensor, we have . The Maxwell equations in the medium with vanishing free fourcurrent then read (73)
We use the ansatz of geometrical optics, (74)
where is a wavelength related parameter, S (x^{α}) is the socalled eikonal function, and 𝔑 denotes the real part. We define (75)
The Maxwell Eqs. (73) then give to the order (76)
Contracting the first equation with and using the second equation, we obtain (77)
therefore, is a null vector of the optical metric. At the same time, from Eq. (75), it follows that (78)
Contracting Eq. (78) with and using Eq. (77), we obtain (79)
which means that 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 . 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 to be tangent to the trajectory x^{β}(a) we search for, that is, (80)
we obtain the equation of motion with an affine parameterization (81)
where the inverse optical metric is given as (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 fourvelocity u^{α} = dx^{α}/dcτ is given as the proper time ∆τ during which the phase changes by 2π in a linear order (see, e.g., Synge 1960 and Chen & Kantowski 2008), so that we have (83)
where the minus sign is there to obtain a positive ∆τ for future oriented . Therefore, the proper frequency ν = 1/∆τ is given as (84)
We consider a signal emitted with the proper frequency ν_{I} by an observer with the fourvelocity from the spacetime point that is received with the proper frequency ν_{F} by an observer with the fourvelocity in the spacetime point lying on the signal trajectory given by the solution x^{α}(a) of the Eqs. (80) and (81), satisfying and . The ratio of the frequencies then is (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 (86)
all along the light ray. We denote κ(λ) = dλ/da the factor converting l^{α} to the affine tangent k^{α}, that is, (87)
The frequency ratio then is (88)
Now we express the ratio κ_{I}/κ_{F} in terms of the fields along the signal trajectory. The geodesic Eq. (79) gives (89)
Contracting Eq. (90) with , we obtain (91)
where in the second line we used , which follows from Eq. (86). Taking into account that the Lie derivative of with respect to the field ξ^{α} is given as (92)
Integrating Eq. (93) along the signal trajectory from λ_{I} to λ_{F}, we obtain (94)
Equation (88) then leads to the following result for the frequency ratio: (95)
Next, we select ξ^{α} to be the coordinate basis vector field corresponding to the time coordinate (i.e., considering the coordinate system (x^{0}, x^{i}) with x^{0} = ct on spacetime, we take ξ^{α} = (∂_{0})^{α}). In this case, the condition (86) leads to (96)
Denoting with υ^{i} = dx^{i} /dt the coordinate velocity of an observer, the fourvelocity of the observer is expressed as (97)
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^{0}. Equation (95) then gives (99)
This is an analogy of the formula (A.46) of Blanchet et al. (2001), generalized to the nonstationary spacetime filled with a nonstationary medium. We also note that for ξ^{α} = (∂_{0})^{α}, we have κ = −k_{0} .
The first fraction in Eq. (99) contains the gravitational redshift and the secondorder Doppler shift. The second fraction is the firstorder 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: (100) (101)
With the choice ξ^{α} = (∂_{0})^{α}, 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 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^{0}/c. For this approximation level, it is sufficient to expand the metric components on the RHS of Eqs. (100) and (101) to include all terms with the power of c^{−1} 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 fourvelocity of the medium U^{α} are (102)
where V^{i}(x^{α}) is the velocity field of the medium in the rotating coordinates given by Eq. (9), and (103)
Using the definition (10), the contravariant components of the optical metric (82) in the rotating coordinates to the approximation level described above are (104a) (104b) (104c)
The covariant components of the optical metric (72) in the rotating coordinates to the same approximation level are (105a) (105b) (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 . Expanding the RHSs of Eqs. (100) and (101) such that they include all terms on the order of 𝒫(2) (we also keep the ∂_{0}A^{i} term to determine how it affects the equation of motion even if it is on the order that we otherwise neglect), we obtain (106) (107) (108)
where we used the fact that (109)
which can be obtained when we express with δ^{ij}a_{i}a_{j} = 1, and we solve the null condition for .
Equation (107) can be inverted to give _{;} as a function of dx^{i} /dλ. Thus, we obtain the following formula including all terms on the order of 𝒫(2): (110)
Differentiating Eq. (107) with respect to λ and using Eqs. (106), (108), and (110) we arrive at the secondorder Newtontype equation of motion for a light ray, with the RHS including all terms on the order of 𝒫(2), (111)
where the antisymmetrization is defined as . 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), (112)
For the purpose of analyzing the twoway 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 (113)
In this case, the spatial part of the light ray x^{i}(l_{E}) has a unit tangent: (114)
Using l^{i} = ∣l^{i}∣dx^{i}/dl_{E} and solving the nullvector condition together with Eq. (86) for l^{0} and ∣l^{i}∣, we obtain the following formulas including all terms on the order of 𝒫(2): (115) (116)
The transformed Eqs. (111) and (112) with RHSs including all terms on the order of 𝒫(2) then are (117)
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 lightray 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_{j}being 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 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 at the event of the emission, the spatial coordinates of the receiver at the event of reception, and either the coordinate time of the emission or the coordinate time of the reception . 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: (119) (120)
For the sake of brevity, we write Eqs. (117) and (118) in the form (121) (122)
where we introduced the coefficients (123) (124) (125) (126) (127) (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 (129)
We denote by the solution of the equation (132)
with satisfying the boundary conditions (133)
where we admit certain deviations , of the boundary points from the boundary points of the complete solution given by Eq. (119), which will be useful for twoway time and frequency transfer. Practically, we use and either or , which is a spatial shift of the observer on the ground in corotating frame during the twoway back and forth propagation of the signal.
Since Eq. (132) is a special case of Eq. (117) with vanishing fields α, ∆W, , and A^{i}, the condition (114) also applies to the solution of Eq. (132), so that we have (134)
which means that is the path of a light ray between the points and for vanishing α, ∆W, , and A^{i}, parameterized by its Euclidean length, and is the total Euclidean length of the path, which can differ from L. Reparameterizing the solution x^{i}(l_{E}) of Eq. (121) to the range [0, ] by introducing a parameter on x^{i}(l_{E}), the solution of the complete Eq. (121) can be expressed as (135)
which defines the correction . Similarly, if x^{0}(l_{E}) is the solution of Eq. (122) for the time coordinate, the reparameterized solution can be expressed as (136)
where is a constant with ∊ [t_{I}, t_{F}] and is the function to be found.
The next step is to express the differential equation for the correction . To do this, we need to express the ratio in order to proceed with the derivatives . Taking a derivative of Eq. (135) with respect to and using Eqs. (114) and (134), we obtain (137)
Inserting Eqs. (135) and (136) into Eq. (121), expanding the equation up to the linear order in ∆x^{α} and and using Eqs. (129) and (132) and the linear expansion of Eq. (137) in we arrive at the equation of the form (138)
with the boundary conditions (139)
which follow from Eqs. (119), (133), and (135), and with the coefficients given by (140) (141) (142) (143) (144)
We proceed with the solution of the particular equations. The solution of the zerothorder 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 with the boundary conditions (139).
5.5.2 Solution for
Integrating Eq. (138) twice, we obtain (145)
with , being integration constants, and [ ] denoting that the function in the square brackets is evaluated in . Determining , from the boundary conditions (139), integrating by parts using (146)
and expressing the integrals over the range ∊ [0,] or ∊ as integrals over the full range ∊ [0, ] using the Heaviside step function , which is defined as for , for and for , we obtain (147)
Integrating by parts again, we can express the term containing in terms of ∆x^{j}. Thus, we obtain the equation of the form (148)
The functions have the following properties, which will be useful later: (151) (152) (153) (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 , 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 , and we define the norm on this vector space as (155)
where f^{1}, f^{2}, f^{3} are the corresponding component functions , which are all continuous. This vector space together with the norm (155) forms a Banach space that we denote B^{3} (see, e.g., Naylor & Sell 1971 for the onedimensional case). The subspace of B^{3} that is formed by functions satisfying is denoted .
The norm on the Banach space B^{3} induces an operator norm on the linear operators from B^{3} to itself. We denote by a bounded linear operator on B^{3}, and by the image of a vector f^{i} ∊ B^{3}. The operator norm of the operator is defined as (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 and that appear in Eq. (150) are continuous on , and we define the linear operators and , as (157) (158)
where the vanishing of the image value in and is ensured by the properties (151) and (152). Next, we denote by b^{i} the linear function satisfying the boundary conditions (139). This function is given as (159)
Using these definitions and denoting by the identity operator on B^{3}, Eq. (148) can be written as (160)
According to theorem 6.7.1 in Naylor & Sell (1971), for a continuous linear operator satisfying , an inverse operator to exists, and it is given by geometric series as (161)
where is the operator given by n consecutive applications of , so that , and for n > 1, we have . 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 can be proven. The assumption 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 as (162)
Using Eq. (161), the solution of Eq. (160) can be written as (163)
where the boundary conditions (139) are ensured by the definition of b^{i} and by the properties of the operators , given by Eqs. (151) and (152). For the derivative , which is also needed below, we obtain (164) (165)
and the operator is defined as (167)
5.5.3 Equation for ∆t and its solution
For further analysis, it is convenient to split the solution (163) for ∆x^{i} into a part that does not depend on ∆t and a part that does depend on it, using Eq. (140). Defining (169) (170)
Inserting Eqs. (135) and (136) into Eq. (122), using Eq. (137) for reparameterization to the parameter , expanding the equation as Taylor series in ∆x^{α} and and using Eqs. (169)–(171), we obtain the following equation for : (172)
where the fields n, W, υ_{Ri}, and A_{i} and their derivatives ∂_{i} are evaluated in . 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^{−1} . This is sufficient to reach the required approximation level for the signal propagation time, as described in Sect. 4.1. All terms containing are on a negligible order.
Now we denote the value of for which , for example, if we select , we have , or if we select , we have . Integrating Eq. (172) from to , we obtain (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^{i}. 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^{−1} only:
, which we obtain by setting and ∂_{t}n ≈ 0 in Eq. (175). Equation (175) then gives (176)
5.6 Coordinate time transfer
5.6.1 Oneway time transfer
In oneway 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 at a time t_{I} and is received in a position at a time t_{F}. Based on Eqs. (120) and (136), we obtain . Using Eq. (176), this leads to (177)
Now we proceed with the integration of ∆σ. We set , and we assume in the boundary conditions (139). We keep nonvanishing as a preparation for later use in the twoway time transfer, where we set for the backpropagating signal 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: , with being the velocity of observer A in the corotating frame in emission event A1. Thus, we consider to be on the order of c^{−1}. In the timetransfer formulas, however, we keep only the highestorder term containing , which is sufficient for an accuracy of 1 ps in case of a stationary ground clock for which the velocity is given by the deformations of Earth, such as solid Earth tides, and by the nonuniformity of the rotation of Earth. In oneway transfer, we can set at the end.
Using the properties (151) and (152) in Eq. (170), we obtain , and therefore, the boundary conditions for are the same as for , namely (178)
Integrating terms containing in ∆σ by parts with use of Eq. (178), using Eqs. (132) and (138) to express the second derivatives and using υ_{Ri} = ϵ_{ijk}ω^{j}x^{k}, we obtain the following expression including all terms on the order of 𝒫(2) except the p_{T} = 2 order terms containing , (179)
where is the firstorder (p_{T} = 1) approximation of S^{i}, which is given as (180)
The in Eq. (179) can be expressed using Eq. (169). In this expansion, only the first term with approximated by is relevant for the required accuracy. Thus, we obtain (181)
where 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 , (182)
In this formula, with given by Eq. (180), the approximation (41) is sufficient to obtain the required accuracy for the resulting timetransfer 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 : (183)
where the fields n, α, W, υ_{Ri}, and A_{i;} and their partial derivatives with respect to the spacetime coordinates are evaluated in if not stated explicitly otherwise. Selecting and , which corresponds to , approximating the terms in the third line by integration over a straight line connecting and instead of , which leads to a negligible difference in terrestrial conditions, and approximating in the fourth line (see Eq. (42)), we obtain Eq. (43) in the summary.
5.6.2 Definitions in twoway transfer
In twoway time transfer, and later also in twoway frequency transfer, we consider a signal that is emitted by a stationary observer on the ground A from the position at the coordinate time t_{A1} to an observer B, who receives the signal in the position at a coordinate time t_{B} and immediately sends it back to A, who receives it in the position 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}, , t_{F+} = t_{I−} = t_{B}, , t_{F−} = t_{A2}, and . Next, we define ) the solution of Eq. (132) with the boundary conditions , and the solution of Eq. (132) with the boundary conditions . It can be shown that (184)
From the boundary conditions for and described above, it follows that and . The shift in the position of A during the back and forth propagation of the signal we denote . For we choose (185)
which leads to and .
5.6.3 Twoway time transfer
Now we can express the time Δt_{+} := t_{F+} − t_{I+} by setting and in Eq. (183) and the time Δt_{−} := t_{F−} −t_{I−} by setting and in Eq. (183), respectively.
Denoting by 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 and expressing t_{A2} − t_{A1} = Δt_{+} + Δt_{−} using Eq. (183), we obtain up to the highest order, (186)
Transforming the integration variables in Δt_{−} as , , using (187) (188)
which follows from Eqs. (184) and (185) and using the properties (153) and (154) of the function we arrive at the following result for the twoway time transfer correction including all terms on the order of 𝒫(3) except the p_{T} = 3 order terms originating from : (189)
where the fields with the index + are evaluated in , for example, , and is a unit vector in the direction from to . This is Eq. (46) in the summary.
5.7 Frequency transfer
5.7.1 Oneway 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 υ^{i} = dx^{i}/dt in the corotating coordinates has the fourvelocity components (in the same coordinates) given by Eq. (97) with (190)
Denoting by , , , and , with , being the spacetime points of the signal emission and reception, respectively, and denoting by , the observers’ velocities at the emission and reception events, we can express the first part of Eq. (99) as (191)
Next, we proceed with the Doppler part of Eq. (99). The components 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): (192)
The normalized tangent to the light ray dx^{i}/dl_{E} can be expressed using Eqs. (135) and (137) up to the second order in as (193)
where, at the end, the terms in the third line do not contribute to the required order of . In order to obtain the boundary values and , we need to express the boundary values of . We introduce the parameter η = 0 or 1, and we evaluate Eq. (165) in . We assume , and we keep nonvanishing as a preparation for the twoway frequency transfer, similarly as in the case of the twoway time transfer. Denoting (194)
We denote by the boundary values of namely and . We write the resulting expression for as the sum (196)
including all terms on the order of 𝒫(2) except terms containing , which are given up to the largest order only. The terms in Eq. (196) are sorted according to their power of c^{−1}, as indicated by their index, and the term containing is given separately. These terms are given by Eqs. (197a)–(197d) below, (197a) (197b) (197c) (197d)
In Eqs. (197), all the fields n, , α, υ_{Ri}, A_{i}, W, and ΔW and their partial derivatives with respect to the spacetime coordinates are evaluated in , if not stated explicitly otherwise. For example, we denote , . The functions evaluated in are denoted by , for example, . The functions and in are the functions and as given by Eqs. (143), (144), and (194) with the level of approximation including all terms with and arbitrary , namely (198) (199)
with and its partial derivatives being evaluated in . The operator in is then generated by the corresponding function (200)
and is S^{i} as given by Eq. (140) with the approximation level including all terms with , p_{α} = 0 and arbitrary , namely (201)
with , and the fields and curl Â^{k} evaluated along .
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 (202)
where the fields ∂_{t}n, ∂_{t}A_{i}, and ∂_{t}W are evaluated along the signal trajectory x^{α}(l_{E}).
Next, we transform the integration variable to l using the firstorder expansion of Eq. (137) in and using Eqs. (135) and (136). Then, taking the Taylor expansion of the integrand in Δx^{α}, assuming , and using the approximations (203) (204)
we obtain the following formula including all terms on the order of 𝒫(3) except terms containing , which are given up to the largest order only: (205a) (205b) (205c) (205d) (205e)
where the fields n, α, A_{i}, and W and their partial derivatives with respect to the spacetime coordinates are evaluated in if not stated explicitly otherwise. Selecting and , which corresponds to , we obtain Eq. (57) in the summary.
5.7.2 Twoway frequency transfer
In twoway 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 and I_{+} by setting , , in Eqs. (196)–(201) and (205), and we define and I_{−} by setting , , in Eqs. (196)–(201) and (205). The proper frequencies ν and the observer velocity quantities υ^{i}, u^{0} for the particular emission and reception events are denoted by the corresponding index A1, A2, or B. Similarly, fields evaluated in the spacetime points , , or of the emission and reception events are denoted by the corresponding index, for example, ,
The goal of this section is to express the correction ∆ defined by Eq. (58). From this definition, it follows that (206)
Using Eq. (99) for the frequency shifts in the particular directions, we obtain (207) (208)
In order to proceed with the correction (206), it is useful to express the relations for the contributions to given by Eqs. (197) and for the contributions to I given by Eqs. (205) between the two directions of the lightray propagation. Transforming the integration variables in the contributions to and in the contributions to I_{−} as , , and so on, using Eqs. (187) and (188), taking into account that , using the properties (153), (154) and taking into account that , we obtain the following relations: (209a) (209b) (209c) (209d)
and (210a) (210b) (210c) (210d)
The u^{0} 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 in terms of the acceleration of the observer A as (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 : (212)
where is given by Eq. (214) below, and β_{i(1)} is its firstorder (p_{T} = 1) part given as (213)
and is given by Eq. (215). (214) (215)
In Eqs. (214) and (215), the index + for fields denotes that the field is evaluated along the path, for example, , . The functions , , and the generating function for the operator in the third and fourth line of Eq. (214) are given by Eqs. (198)–(201) with .
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^{−1} 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^{−1} N_{B}, we can write (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 endpoint tangent, we use Eq. (40) expanded up to the first order, (217)
where the index B replaces the general index for a final point F. Where possible, the approximation (41) is used. After expanding in the integral term in the second line of Eq. (214), we use the following way of integrating the term containing ω^{l}: (218)
where the second derivative of the path 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 can be processed as follows: (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 . The use of these relations leads to the cancellation of several terms in the expression for β_{i(2)}.
6 Conclusion
A relativistic model of oneway and twoway 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^{i} (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^{i}) corotating with the Earth. This equation is given as the null geodesic equation of the Gordon optical metric (Gordon 1923). First, the exact solution 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 , 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 frequencytransfer corrections are then expressed in terms of integrals of the fields along the path at a certain constant coordinate time.
In the numerical examples evaluated in this paper, we focused on the twoway groundtosatellite transfer with a satellite altitude similar to that of the ISS. In the twoway 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 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/reemission 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 twoway time transfer is given by the difference in the propagation times caused by the FresnelFizeau effect of light dragging. It is much smaller than 1 ps for normal atmospheric conditions.
In the twoway 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^{−17} for θ = 0, reaches 10^{−16} for θ ≈ 56°, and ends at 10^{−13} for θ approaching 90°.
Remarkably, in the equation of motion for a light ray, the quantity A^{i} = (1 − n^{2})V^{i} 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 twoway 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^{−17} already for a wind speed of about 11 m s^{−1} (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 frequencytransfer corrections including the effect of the wind shows that they need to be considered in the data analysis of the forthcoming clockonsatellite experiments such as ACES or ISOC. The results obtained in this paper can be used not only to calculate the time and frequencytransfer 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 (A.1)
where ρ_{1}, ρ_{2}, and ρ_{3} are the partial densities of dry, CO_{2} free air, of water vapor, and of carbon dioxide in the mixture, and R_{1}, R_{2}, and R_{3} 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 ρ_{1}(x^{α}), ρ_{2}(x^{α}), and ρ_{3}(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 (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_{2}, ρ_{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_{2} 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 (A.3)
where N_{0} and ρ_{0} 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 (A.4)
being the gravitational plus centrifugal potential. Next, we assume the perfect gas equation of state for the air, (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 (A.7)
where ϕ_{s} is the sealevel 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 (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 (A.9)
and for β = 0, it reads (A.10)
A.3 Spherically symmetric, static component of the atmospheric refractivity
The time and frequencytransfer model developed in this paper uses a spherically symmetric, static component of the refractive index n(r) (see Eq. (6)) or the related refractivity . 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_{0}. From Eqs. (A.10) and (A.3), we then obtain (A.11) (A.12)
where ρ_{0}, N_{0}, and ϕ_{0} are the air density, air refractivity, and the potential in a certain reference spacetime point. The refractivity is then given as (A.13)
where Ŵ_{0} is the value of Ŵ in the reference point.
References
 Alonso, I., Alpigiani, C., Altschul, B., et al. 2022, EPJ Quantum Technol., 9, 30 [CrossRef] [Google Scholar]
 Ashby, N. 1998, in IEEE International Frequency Control Symposium, 320 [Google Scholar]
 Belmonte, A., Taylor, M. T., Hollberg, L., & Kahn, J. M. 2017, Opt. Express, 25, 15676 [NASA ADS] [CrossRef] [Google Scholar]
 Beloy, K., Bodine, M. I., Bothwell, T., et al. 2021, Nature, 591, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Blanchet, L., Salomon, C., Teyssandier, P., & Wolf, P. 2001, A&A, 370, 320 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bodine, M. I., Deschenes, J.D., Khader, I. H., et al. 2020, Phys. Rev. Res., 2, 033395 [NASA ADS] [CrossRef] [Google Scholar]
 Born, M., & Wolf, E. 1999, Principles of Optics, 7th edn. (Cambridge University Press) [CrossRef] [Google Scholar]
 Bourgoin, A. 2020, Phys. Rev. D, 101, 064035 [Google Scholar]
 Bourgoin, A., Zannoni, M., & Tortora, P. 2019, A&A, 624, A41 [EDP Sciences] [Google Scholar]
 Bourgoin, A., Zannoni, M., Gomez Casajus, L., Tortora, P., & Teyssandier, P. 2021, A&A, 648, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cacciapuoti, L., & Schiller, S. 2017, ISOC Scientific Requirements, SCIESA HREESRISOC document [Google Scholar]
 Cacciapuoti, L., Armano, M., Much, R., et al. 2020, Eur. Phys. J. D, 74, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Cao, J., Zhang, P., Shang, J., et al. 2017, Appl. Phys. B, 123, 112 [CrossRef] [Google Scholar]
 Chen, B., & Kantowski, R. 2008, Phys. Rev. D, 78, 044040 [Google Scholar]
 Ciddor, P. E. 1996, Appl. Opt., 35, 1566 [Google Scholar]
 Delva, P., Hees, A., & Wolf, P. 2017, Space Sci. Rev., 212, 1385 [NASA ADS] [CrossRef] [Google Scholar]
 Delva, P., Puchades, N., Schönemann, E., et al. 2018, Phys. Rev. Lett., 121, 231101 [NASA ADS] [CrossRef] [Google Scholar]
 Delva, P., Denker, H., & Lion, G. 2019, Chronometric geodesy: methods and Applications, in Relativistic Geodesy: Foundations and Applications, eds. D. Pützfeld, & C. Lammerzähl (Switzerland AG: Springer Nature) [Google Scholar]
 Derevianko, A., Gibble, K., Hollberg, L., et al. 2022, Quantum Sci. Technol., 7, 044002 [NASA ADS] [CrossRef] [Google Scholar]
 DixMatthews, B. P., Schediwy, S. W., Gozzard, D. R., et al. 2021, Nat. Commun., 12, 515 [NASA ADS] [CrossRef] [Google Scholar]
 Djerroud, K., Acef, O., Clairon, A., et al. 2010, Opt. Letters, 35, 1479 [NASA ADS] [CrossRef] [Google Scholar]
 Duchayne, L., Mercier, F., & Wolf, P. 2009, A&A, 504, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feng, J. C., Hejda, F., & Carloni, S. 2022, Phys. Rev. D, 106, 044034 [NASA ADS] [CrossRef] [Google Scholar]
 Fujieda, M., Piester, D., Gotoh, T., et al. 2014, Metrologia, 51, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Gordon, W. 1923, Ann. Phys., 72, 421 [NASA ADS] [CrossRef] [Google Scholar]
 Gozzard, D. R., Schediwy, S. W., Stone, B., Messineo, M., & Tobar, M. 2018, Phys. Rev. Appl., 10, 024046 [CrossRef] [Google Scholar]
 Hachisu, H., Fujieda, M., Nagano, S., et al. 2014, Opt. Lett., 39, 4072 [NASA ADS] [CrossRef] [Google Scholar]
 Hees, A., Bertone, S., & Le PoncinLafitte, C. 2014a, Phys. Rev. D, 89, 064045 [Google Scholar]
 Hees, A., Bertone, S., & Le PoncinLafitte, C. 2014b, Phys. Rev. D, 90, 084020 [Google Scholar]
 ICAO. 1993, Manual of the ICAO Standard Atmosphere extended to 80 kilometres (262 500 feet), 3ed edn., Doc 7488/3 [Google Scholar]
 Kang, H. J., Yang, J., Chun, B. J., et al. 2019, Nat. Commun., 10, 4438 [NASA ADS] [CrossRef] [Google Scholar]
 Koller, S. B., Grotti, J., Vogt, S., et al. 2017, Phys. Rev. Lett., 118, 073601 [NASA ADS] [CrossRef] [Google Scholar]
 Kopeikin, S. M., & Han, W.B. 2015, J. Geod., 89, 829 [NASA ADS] [CrossRef] [Google Scholar]
 Kopeikin, S., & Mashhoon, B. 2002, Phys. Rev. D, 65, 064025 [Google Scholar]
 Kopeikin, S. M., & Schäfer, G. 1999, Phys. Rev. D, 60, 124002 [NASA ADS] [CrossRef] [Google Scholar]
 Kopeikin, S., Efroimsky, M., & Kaplan, G. 2011, Relativistic Celestial Mechanics of the Solar System (WileyVCH Verlag GmbH &Co. KGaA) [CrossRef] [Google Scholar]
 Landau, L. D., & Lifshitz, E. M. 1976, Mechanics, 3rd edn. (ButterworthHeinemann) [Google Scholar]
 Le PoncinLafitte, C., Linet, B., & Teyssandier, P. 2004, Class. Quant. Grav., 21, 4463 [CrossRef] [Google Scholar]
 Leonhardt, U., & Piwnicki, P. 1999, Phys. Rev. A, 60, 4301 [NASA ADS] [CrossRef] [Google Scholar]
 Lilley, M., Savalle, E., Angonin, M. C., et al. 2021, GPS Solutions, 25, 34 [CrossRef] [Google Scholar]
 Linet, B., & Teyssandier, P. 2002, Phys. Rev. D, 66, 024045 [Google Scholar]
 Lisdat, C., Dörscher, S., Nosske, I., & Sterr, U. 2021, Phys. Rev. Res., 3, L042036 [NASA ADS] [CrossRef] [Google Scholar]
 Ludlow, A. D., Boyd, M. M., Ye, J., Peik, E., & Schmidt, P. O. 2015, Rev. Mod. Phys., 87, 637 [NASA ADS] [CrossRef] [Google Scholar]
 McGrew, W. F., Zhang, X., Fasano, R. J., et al. 2018, Nature, 564, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Mehlstäubler, T. E., Grosche, G., Lisdat, C., Schmidt, P. O., & Denker, H. 2018, Rep. Prog. Phys., 81, 064401 [CrossRef] [Google Scholar]
 Müller, J., Dirkx, D., Kopeikin, S. M., et al. 2018, Space Sci. Rev., 214, 5 [CrossRef] [Google Scholar]
 Naylor, A. W., & Sell, G. R. 1971, Linear Operator Theory in Engineering and Science (Holt, Rinehart and Wilson, Inc.) [Google Scholar]
 NOAA. 1976, U. S. Standard Atmosphere, Report No. NOAAS/T 761562 [Google Scholar]
 Origlia, S., Pramod, M. S., Schiller, S., et al. 2018, Phys. Rev. A, 98, 053443 [NASA ADS] [CrossRef] [Google Scholar]
 Owens, J. C. 1967, Appl. Opt., 6, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Perlick, V. 2000, Ray Optics, Fermat’s principle, and Applications to General Relativity (SpringerVerlag Berlin Heidelberg) [Google Scholar]
 Petit, G., & Luzum, B. 2010, IERS Conventions (2010), IERS Technical Note No. 36, Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt am Main [Google Scholar]
 Pizzocaro, M., Sekido, M., Takefuji, K., et al. 2021, Nat. Phys., 17, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Riedel, F., AlMasoudi, A., Benkler, E., et al. 2020, Metrologia, 57, 045005 [NASA ADS] [CrossRef] [Google Scholar]
 Robert, C., Conan, J.M., & Wolf, P. 2016, Phys. Rev. A, 93, 033860 [Google Scholar]
 Safronova, M. S., Budker, D., DeMille, D., et al. 2018, Rev. Mod. Phys., 90, 025008 [CrossRef] [Google Scholar]
 Schreiber, K. U., Procházka, I., Lauber, P., et al. 2010, IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 57, 728 [Google Scholar]
 Shen, Q., Guan, J.Y., Zeng, T., et al. 2021, Optica, 8, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Sinclair, L. C., Giorgetta, F. R., Swann, W. C., et al. 2014, Phys. Rev. A, 89, 023805 [NASA ADS] [CrossRef] [Google Scholar]
 Sinclair, L. C., Swann, W. C., Bergeron, H., et al. 2016, Appl. Phys. Lett., 109, 151104 [CrossRef] [Google Scholar]
 Soffel, M., Klioner, S. A., Petit, G., et al. 2003, AJ, 126, 2687 [Google Scholar]
 Stuhl, B. K. 2021, Opt. Express, 29, 13706 [NASA ADS] [CrossRef] [Google Scholar]
 Swann, W. C., Bodine, M. I., Khader, I., et al. 2019, Phys. Rev. A, 99, 023855 [NASA ADS] [CrossRef] [Google Scholar]
 Synge, J. L. 1960, Relativity: The General Theory (Amsterdam: NorthHolland Publishing Company) [Google Scholar]
 Taylor, M. T., Belmonte, A., Hollberg, L., & Kahn, J. M. 2020, Phys. Rev. A, 101, 033843 [Google Scholar]
 Teyssandier, P., & Le PoncinLafitte, C. 2008, Class. Quant. Grav., 25, 145020 [NASA ADS] [CrossRef] [Google Scholar]
 Tselikov, A., & Blake, J. 1998, Appl. Opt., 37, 6690 [CrossRef] [Google Scholar]
 Turyshev, S. G., Yu, N., & Toth, V. T. 2016, Phys. Rev. D, 93, 045027 [NASA ADS] [CrossRef] [Google Scholar]
 Voigt, C., Denker, H., & Timmen, L. 2016, Metrologia, 53, 1365 [NASA ADS] [CrossRef] [Google Scholar]
 Wolf, P., & Petit, G. 1995, A&A, 304, 653 [NASA ADS] [Google Scholar]
All Figures
Fig. 1 Definition of angles, basis vectors, and other quantities. 

In the text 
Fig. 2 Scheme of the twoway time transfer. The signal leaves observer A at a coordinate time t_{A1}, corresponding to the proper time τ_{A}(t_{A1}) of clock A, is reflected from observer B at a coordinate time t_{B}, corresponding to the proper time τ_{B}(t_{B}) of clock B and τ_{A}(t_{B}) of clock A, and is finally received by observer A at the coordinate time t_{A2}, corresponding to the proper time τ_{A}(t_{A2}) of clock A. 

In the text 
Fig. 3 Visualization of the Sagnac area Σ. The part of the Sagnac effect that is caused by the atmosphere is given by the change in the area Σ when the path changes from the case including atmospheric refraction and gravitation to a vacuum case that only includes gravitation ( in Eq. (13)). 

In the text 
Fig. 4 Sagnac correction, Eq. (48), for an isothermal atmosphere and its deviation from a vacuum value denoted as effect of the atmosphere as functions of the satellite zenith angle θ. Separate plots for θ ∊ [0°, 80º] (left) and θ ∊ [80°, 90º] (right) are presented to show details of the effect of the atmosphere curve, which grows steeply for θ approaching 90°. The example was evaluated for groundtosatellite transfer in the equatorial plane, with r_{A} = 6371 km, r_{B} − r_{A} = 408 km, and inclined against the rotation of the Earth. 

In the text 
Fig. 5 Three contributions of the spherically symmetric, static part of the atmospheric refractivity , and to the twoway frequency transfer correction (59) and their sum as functions of the satellite zenith angle θ. The example is evaluated for groundtosatellite transfer through an isothermal atmosphere in hydrostatic equilibrium with 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_{A} + 408 km in the direction of the rotation of the Earth. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.