Relativistic theory for time and frequency transfer through flowing media with an application to the atmosphere of Earth

Context. Several space missions that will use atomic clocks on board of an Earth-orbiting 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 (I-SOC). 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 I-SOC project is about 10 − 17 . Aims. We develop a relativistic model of one-and two-way 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 two-way 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 two-way ground-to-satellite 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 two-way time transfer and from 10 − 17 to 10 − 13 for two-way 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 two-way time transfer for normal atmospheric conditions, but for the two-way 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 ms − 1 . Conclusions. The atmospheric effects including the effect of wind should be considered in the forthcoming clock-on-satellite experiments such as ACES or I-SOC.


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;Safronova et al. 2018;Delva 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. (2021).
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 free-space optical links (see Bodine et al. (2020) and references therein and also Djerroud et al. (2010), Gozzard et al. (2018), Kang et al. (2019), Dix-Matthews et al. (2021) and Shen et al. (2021)) enable increasing accuracies of the experiments.Therefore, the atmospheric effects originating from spatial and temporal variations of the refractive index of air or from winds start to play a more important role.
In case of the ground-to-satellite time and frequency transfer, two projects are going to be realized in the near future that will place an atomic clock on the International Space Station (ISS), namely the Atomic Clock Ensemble in Space (ACES), with a Cs clock and H-maser (Lilley et al. 2021;Cacciapuoti et al. 2020), and the Space Optical Clock on the ISS (I-SOC), with an optical lattice clock (Cacciapuoti & Schiller 2017;Origlia et al. 2018).In the ACES mission, the absolute frequency accuracy of the on-board clock is about 10 −16 in fractional frequency, whereas the same parameter for the I-SOC 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 100 ps.For the ELT, the overall planned accuracy of the time transfer is 50 ps (Turyshev et al. 2016).The I-SOC proposal uses an improved version of the MWL as well as the optical link of ELT, aiming at an accuracy of a few ps in the clock synchronization (Cacciapuoti & Schiller 2017).Therefore, based on the planned accuracies of these ex-periments, a relativistic model of one-and two-way time and frequency transfer is needed that includes the atmospheric effects and has an accuracy of approximately 1 ps for the time transfer and a relative accuracy of the lower multiples of 10 −18 for the frequency transfer.Blanchet et al. (2001) developed the relativistic theory of one-and two-way time and frequency transfer in vacuum, including terms up to the order of c −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énard-Wiechert representation of the metric tensor up to the first post-Minkowskian order was used.Linet & Teyssandier (2002) applied the Synge world function formalism to the gravitational field of an axisymmetric rotating body, and they showed that certain c −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 Poncin-Lafitte et al. (2004) and Teyssandier & Le Poncin-Lafitte (2008), where the formalism of time transfer functions was introduced.Hees et al. (2014b) and Hees et al. (2014a) employed this formalism to compute various observables relevant for actual space missions in the Solar System, including the coordinate propagation time and the frequency shift of light.
The relativistic theoretical works mentioned above considered the vacuum case alone and did not take the atmospheric effects into account.The theory of Teyssandier & Le Poncin-Lafitte (2008) was extended to include the atmospheric effects by Bourgoin (2020).Here, a general formalism of time transfer functions was developed using the Gordon optical metric (Gordon 1923), and it was applied to a case of stationary optical spacetime.The effects of refractivity and motion of a neutral atmosphere were added to the gravitational effects, and the 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 one-way time and frequency transfer only.The results of Bourgoin (2020) were later used to model the atmospheric occultation experiments (Bourgoin et al. 2021).In this work, explicit formulas for time and frequency transfer were derived for a steadily rotating spherically symmetric atmosphere up to the first post-Minkowskian order, with the remark that the method described in the paper easily allows extending the results to higher orders and beyond the spherical symmetry.Feng et al. (2022) used the formalism of Gordon's optical metric in the context of global navigation satellite systems to include the atmospheric effects in the relativistic framework.The null geodesics of the Gordon metric were solved numerically.The effect of motion of the atmosphere of Earth on the speed of light (the Fresnel-Fizeau effect) and the related propagation time delay were discussed also by Kopeikin & Han (2015) in the context of geodetic very long baseline interferometry.
Several effects originating from the atmospheric refractivity have been addressed in the literature in a nonrelativistic framework.The influence of the refractive index fluctuations given by atmospheric turbulence has been discussed by Sinclair et al. (2014), Robert et al. (2016), Sinclair 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 two-way ground-to-satellite time transfer due to the distribution of the atmospheric refractive index was discussed by Stuhl (2021).Earlier findings on tropospheric and ionospheric corrections were summarized, for example, in Petit & Luzum (2010).A two-way time transfer model including the effect of ionospheric dispersion was discussed by Duchayne et al. (2009), who also provided uncertainty requirements on the orbit determination of space clocks.
The aim of this paper is to develop a relativistic model of one-and two-way time and frequency transfer that includes gravitational and atmospheric effects for realistic fields of the gravitational potential, the refractive index, and the wind speed, assuming neither stationarity nor spherical or axial symmetry.The intended accuracy of the model is given by the experimental needs described above: approximately 1 ps for time transfer, and lower multiples of 10 −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 two-way time and frequency transfer are given, and several examples are shown that quantify the atmospheric effects for some typical situations whose parameters are similar to the experiments at the ISS.A complete derivation of the results is then given in Sect. 5.In the appendix, some basic facts about the refractive index of air and its distribution in the atmosphere of Earth are summarized.

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 |a i | ≡ δ i j a i a j , with δ i j being the Kronecker delta.Quantities , and γ i in this paper are treated as Euclidean three-vectors when their index is lowered or raised, for example, V i = δ i j V j .The three-dimensional antisymmetric Levi-Civita symbol is denoted as ϵ i jk , eventually ϵ i jk = δ im ϵ m jk .
the center of mass of the Earth and is nonrotating with respect to distant stars.We denote the GCRS coordinates as xα = ( x0 , xi ), where x0 /c = t is the Geocentric Coordinate Time (TCG), and xi 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) with W being the scalar gravitational potential defined by Soffel et al. (2003) with the convention W ≥ 0.
To evaluate the time and frequency transfer, however, we use coordinates that rotate with respect to the GCRS system, with a constant angular velocity vector approximating the real angular velocity vector of the rotation of the Earth around its axis, such that the surface of the Earth is nearly at rest in the new coordinates, and we can use the advantage of nearly vanishing velocities of stationary observers on the ground.The use of this corotating system leads to more complicated formulas for oneway time and frequency transfer, but, on the other hand, the fact that the position of the stationary observer on the ground is almost constant makes the calculation of two-way corrections more straightforward.The rotating coordinate system x α = (x 0 , x i ) is related to the GCRS system as The time-dependent 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 where ω j are constant components of the angular velocity vector of the coordinate rotation expressed in the xi frame.Since ω j is defined to be constant and the coordinate rotation is rigid, the irregularities in the rotation of the Earth and the deformations of the Earth given mainly by the solid Earth tides are observed as tiny movements of the surface of the Earth in the coordinates x i .We denote ṽi R = d xi /dt = ϵ i jk ω j xk 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 The components of the metric (1) in the rotating coordinates are where The monopole part of the scalar gravitational potential is denoted as Ŵ(r), and the remaining part containing the higher multipole moments and the time-dependent components as ∆W(x α ) such that 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).

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 part A.1 of the appendix).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 part A.2 of the appendix).Spatial and temporal local fluctuations are present as well.Therefore, it is reasonable to decompose the atmospheric refractive index field as where n(r) 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 n(r), α(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 n(r) for altitudes below 80 km are given in the appendix.The correction α(x α ) then consists of various more or less predictable components, such as the deviation from spherical symmetry caused by the centrifugal potential, the effects of the changing position of Sun as a heat source, the effects of large-scale meteorological structures, or components coming from variations in local temperature, pressure, and composition, and from atmospheric turbulence.
The total refractivity N(x α ) and the static, spherically symmetric refractivity N(r) are then defined as deviations of the corresponding refractive index from its vacuum value, Denoting x i M (t, x α 0 ) the trajectory of the medium element in the rotating coordinates that passes through a spacetime point x α 0 = (ct 0 , x i 0 ), the winds in the atmosphere are described by their velocity field V i (x α ), which is defined by for any x α 0 where the medium is present.We also define the following quantity that is useful when the optical effects of the winds are evaluated: (10)

Observers
For a given trajectory x i (t) of an observer in the rotating coordinates, we denote v i = dx i /dt the observer's velocity in these coordinates.

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

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 clock-on-satellite experiments.The corresponding approximation level of the model should be defined.
The influencing quantities entering the equation of motion and the evaluation of the time and frequency transfer in the rotating frame originate from gravitation given by the Newtonian potential W/c 2 (higher-order gravitational effects such as the Lense-Thirring effect are negligible considering our intended accuracy), from inertial effects given by the rotation velocity v i R /c, from atmospheric effects given by the air refractivity N together with the correction α and by the wind speed V i /c and from the observers' velocities v 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 N(r) and its correction α(x α ), b) we denote NM the maximum value of N and α M the maximum value of |α|, c) we define normalized functions N * = N/ NM and α * = α/α M and in all formulas we express N = NM N * and α = α M α * , d) the order of a term in an expansion is then given by its factor (c −1 ) p c −1 ( NM ) p N (α M ) p α , where the powers p c −1 , p N , 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 (p c −1 , p N , p α ).The expansion of n in terms of NM , N * , α M , and α * serves for the purpose of defining the order, but is not shown explicitly in formulas where it is not practical.
First, considering terms on the order of (c −1 ) p c −1 ( NM ) p N (i.e., p α = 0) in the time and frequency transfer formulas, we take p c −1 ≤ 3 corresponding to the vacuum model of Blanchet et al. (2001), and we take an arbitrary value of p N for p c −1 ≤ 2. Terms on the order of c −3 ( NM ) p N with p N ≥ 1 are neglected.Some of the terms on the order c −3 NM can slightly exceed the 10 −18 level in case of ground-to-ISS 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 ground-to-ISS frequency transfer, as was shown by Linet & Teyssandier (2002).Then, taking the field α into account, its effect is considered to be minor compared to the spherical, static part of the refractivity N. We therefore only include the field α to a lower order.Namely, for p α ≥ 1, we take p c −1 + p N + p α ≤ 3. Possible magnitudes of the neglected terms with p α ≥ 1 have not been evaluated, however.
Next, we denote the total p c −1 + p N + p α = p T , and we can summarize the approximation level of the model as follows: the expressions for the signal propagation time, the relative frequency shift, and the related two-way time and frequency transfer corrections in this paper include all terms with the order from a set P(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 two-way corrections that originate from the motion of a stationary observer on the ground in the corotating frame during the back and forth propagation of the signal.This motion is given by the deviations of the surface motion of Earth from rigid uniform rotation and therefore is very small.The approximation level of these terms is discussed together with the particular formulas.
Next we introduce the notation O(q) for expressions containing terms with p T ≥ q only.Expressions containing terms with p N = p α = 0 and p c −1 ≥ k only are denoted O(c −k ) as usual.It is also useful to define the following subset of the exponent triplets:

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 xi ( l), with l being a parameter given by the Euclidean length of this path from its initial point.Then, an approximate linearized equation is solved for the path correction ∆x i ( l) caused by the deviations of the refractive index and the gravitational potential from sphericity and staticity, by winds and by the Coriolis and centrifugal forces.Finally, a separate first-order differential equation for the coordinate time t( l) is solved.
The main results of this paper are expressed in terms of the integrals along the path xi ( l).We therefore discuss the solution of the static, spherically symmetric part of the equation of motion in detail in this summary.

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 n(r) and gravitational potential Ŵ(r) is given as (see Eq. ( 132)) where the fields n, ∂ j n, and ∂ j Ŵ are evaluated in xk ( l), and the solution has a unit tangent satisfying corresponding to the fact that l is its Euclidean length parameter.We assume that the initial and final points of the path are known as input parameters, and we later relate them to spatial coordinates of the emission and reception events of the one-or two-way time and frequency transfer setups.We denote the total path length between the initial and final points as L, such that the parameter has the range l ∈ [0, L].Defining a field w as Eq. ( 11) can be written as where we briefly denote xi (ζ) ≡ xi ( l(ζ)).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 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 xi (ζ) and using the assumption of the central field, which leads to 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 ψ( l) the angle between the tangent d xi ( l)/d l and the position vector xi ( l) and taking the magnitude of Eq. ( 17) with use of Eq. ( 16), we obtain where r( l) = | xi ( l)| 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 l(r) being the inverse function to r( l) (assuming dr/d l 0 all along the trajectory), Eq. ( 18) can be written as Selecting a right-handed Cartesian coordinate system y i with the origin at r = 0, the y 1 -axis in xi (0) direction, and the y 3 -axis in the h i direction (for h = 0, we can take any direction perpendicular to xi (0)), defining polar coordinates by y 1 = r cos φ and y 2 = r sin φ and expressing the trajectory xi ( l) in the polar coordinates as φ(r), with the radius as the parameter, it follows from the definition of ψ that Integrating this equation with tan ψ expressed from the Snell law (19), we obtain with ± = + for dr/d l > 0 (ascending trajectory) and ± = − for dr/d l < 0 (descending trajectory).Thus, up to the final integration, we obtain the solution of Eq. ( 11) in polar coordinates parameterized by r.The remaining step is to express the constant h in terms of the input parameters, which, in this case, are the boundary values of the trajectory xi ( l).
In the following, we denote quantities related to the initial point xi (0) by the index I and to the final point xi ( L) by the index F, for example, xi , and ψ F = ψ( L).We also define the angle φ between xi  1).The following relations hold: 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 d xi /d l| 0 and d xi /d l| L, respectively, from the vector xi F − xi I (see Fig. 1), so that we have ε Evaluating the Snell law ( 19) in the boundary points and using Eqs.( 24) and ( 25), we obtain 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 d xi /d l| 0 and d xi /d l| L, which satisfies Solving Eq. ( 26) together with Eq. ( 27) for the angles ε I and ε F , we obtain The total bending angle also satisfies where the angle φ can be expressed as φ = φ(r F ) using Eq. ( 21), and we can also write where on the right-hand side (RHS), we use Eq. ( 19).Thus, we obtain 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 ground-to-satellite or satellite-toground 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 φ 0, we have The solution parameterized by r is then given as For φ = 0, it is simply xi ( l(r)) = re i 1 .
The reparameterization function l(r) can be expressed with the use of where the ± sign has the same meaning as in Eq. ( 21).Integrating this formula using the Snell law ( 19), we obtain For L we then get L = l(r F ).
It is also useful to define an orthonormal right-handed basis (χ i , ς i , γ i ) with χ i pointing in direction from xi I to xi F , γ 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) Defining the angle ε( l) between the path tangent d xi /d l and the vector χ i satisfying ε(0) = ε I (negative value) and ε( L) = ε F (positive value), we can write In some formulas, the lowest-order approximation of the tangent d xi /d l and of the length L is sufficient.Expanding Eq. ( 40) in ε and using |ε| ≤ ε T = O(1), which follows from Eq. ( 32), we obtain Next we obtain where the expansion of cos ε was used in the last step.
In the following text, we refer to the solution xi ( l) of Eq. ( 11) with the given boundary conditions as the base path.

One-way time transfer
In one-way time transfer, we wish to express a coordinate time t F − t I of the propagation of the signal emitted from a position x i I at a time t I and received in a position x i F at a time t F .We consider a base path xi ( l) with boundary points xi (0) = x i I and xi ( L) = x i F .Based on the analysis given in Sect.5, we obtain the following formula for the propagation time, including all terms on the order of P(3): where the fields n, α, W, v Ri , A i and their partial derivatives with respect to the space-time coordinates are evaluated in t = t I and , respectively, and the function D( l, l′ ) 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 Fresnel-Fizeau effect of dragging of light by a medium.The remaining terms containing α are the effects of nonsphericity and nonstationarity of the refractive index field.Numerical integration of the terms in Eq. ( 43) can be performed, for example, by transforming the integration variable l into the radial coordinate r using Eq. ( 36) and by expressing the integrated functions in spherical coordinates with the base path given by Eq. ( 21).

Two-way time transfer
In two-way transfer, we consider the signal emitted from a stationary observer on the ground A at the coordinate time t A1 , which corresponds to the proper time τ A (t A1 ) of clock A. The spatial coordinates of this event are denoted x i A1 .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 ) is obtained from the measurements.Then a signal is sent from observer B at the coordinate time t B and is received by observer A at the coordinate time t A2 , which corresponds to the proper time τ A (t A2 ) of clock A. The spatial coordinates of this event are denoted x i A2 .This signal can either be the same signal reflected or another signal that is synchronously sent when the one-way signal is received by observer B. We call this set-up the Λ-configuration.

The pseudo-time-of-flight τ
Using the coordinate-time synchronization convention, the desynchronization of the clocks is defined as the difference τ B (t B ) − τ A (t B ), where τ A (t B ) is the proper time measured by observer A at the coordinate time t B .In case of the two-way time transfer, it can be expressed as where The first two lines of Eq. ( 44) are the measured pseudo-times-of-flight, and the difference ∆τ − − ∆τ + needs to be computed.Denoting ∆t + = t B − t A1 the coordinate time of the signal propagation from A to B and ∆t − = t A2 − t B the coordinate time of the signal propagation from B to A, the computed term can be approximated as The error of this approximation is several orders below the 1 ps level in terrestrial conditions.We also note that in the corotating coordinates we use, the position of the stationary observer on the ground A is nearly fixed and changes only due to the deformations of Earth, which are given, for example, by the solid Earth tides or by the irregularity of the rotation of Earth.Thus, the difference between x i A1 and x i A2 is very small, and we include its effect up to the largest order only using the velocity v i A1 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 xi + ( l) the base path with the boundary points xi , we arrive at the following result for the two-way time transfer correction, including all terms on the order of P(3) except the p T = 3 order terms originating from where the fields v Ri , A i , ∂ i α, and ∂ t α with the index + are evaluated in t = t B and x i = xi + ( l), and χ i + is the vector χ i corresponding to the path xi + ( l) (i.e., it is the unit vector in the direction from x i A1 to x i B ).The main contribution to the two-way time transfer correction in Eq. ( 46) is given by the first term in the first line, which is the well-known Sagnac effect.We denote Then, taking into account that v Ri+ = ϵ i jk ω j xk + ( l), the correction given by Eq. ( 47) can be expressed as where Σ j is defined as The magnitude Σ = |Σ j | equals the area of the plane section surrounded by three lines: 1) the straight line from the origin of the coordinate system (the Earth center) to x i A1 , 2) the path xi + ( l) from x i A1 to x i B , and 3) the straight line from x i B 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 xi + ( l) 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 where , γ j+ is the vector γ j corresponding to the path xi + ( l), and w(r) is given by Eq. ( 13).This integral can be computed numerically for any refractive index profile n(r) 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 The part of the Sagnac effect that is caused by the atmosphere is given by the change in the area Σ when the path xi + ( l) changes from the case including atmospheric refraction and gravitation to a vacuum case that only includes gravitation (n = 1 in Eq. ( 13)).
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.

Sagnac correction
In this section, we evaluate the Sagnac correction, Eq. (48), for the example of ground-to-satellite transfer through an isothermal atmosphere at a temperature T 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 N(r) = n(r) − 1 is given by Eq. (A.13), which gives the following spherical refractive index field: where Ŵ = GM/r, ŴA = GM/r A , R is the universal gas constant, and N A is the air refractivity in the reference position x i A1 , which also equals N(r A ) (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.15K (15 • C) and M a = 28.964× 10 −3 kg/mol, corresponding to dry air with a carbon dioxide molar fraction x c = 450 ppm.For an atmospheric pressure of 101 325 Pa in the position x i A1 and a light signal with a vacuum wavelength of 1 µm, we obtain (see, e.g., Ciddor (1996)) (52)  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 ground-to-satellite transfer in the equatorial plane, with r A = 6371 km, r B − r A = 408 km, and χ i + inclined against the rotation of the Earth.
In two-way transfer, we denote the θ I angle of the xi + ( l) path simply as θ.It represents the zenith angle of the satellite at the reception/re-emission event (see Fig. 3).
Based on Eqs. ( 48) and ( 50), the Sagnac correction can be calculated as a function of θ.In Fig. 4 this function is shown (the dashed blue curve) together with the effect of the atmosphere (the red full curve), which is defined as the Sagnac correction in the atmosphere (N A given by Eq. ( 52)) minus the Sagnac correction in vacuum (N A = 0) for the same value of θ.The values in the plots in Fig. 4 are obtained for the case ω i γ i+ = −ω, which corresponds to the paths xi + ( l) propagating in the equatorial plane against the rotation of the Earth.This is the case with the highest positive values of the Sagnac correction.
Fig. 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 xi + ( l) is much smaller than the effect of atmospheric refraction, the path xi + ( l) 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,

Effect of the wind
We evaluate the effect of the wind speed in the two-way time transfer correction (46), which is given by 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 xi + ( l), as well as the component of the air velocity V i (in the corotating frame) in the direction of the path tangent d xi + /d l, which we denote V.The approximately constant n corresponds, for example, to the horizontal signal path from the ground to the ground or to the initial part of the signal path from the ground to the satellite when the satellite is in a position near the horizon.
In this case, we obtain For a 1 ps correction in the air with the refractivity value of Eq. ( 52), which is typical on the surface of the Earth, we need V L ≈ 8.2 × 10 7 m 2 /s, which, for example, for the path length of L = 100 km leads to V ≈ 820 m/s.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.One-way frequency transfer In one-way frequency transfer, the goal is to express the ratio of the proper frequency ν I of a signal emitted by an observer with velocity v i I from a position x i I at a time t I and a proper frequency ν F of the same signal as is received by the observer with a velocity v i F in the position x i F at the time t F .Again, we consider a base path xi ( l) with the boundary points xi (0) = x i I and xi ( L) = x i F , and we denote Based on the analysis included in Sect.5, we obtain the following formula for the frequency ratio: where ( li ) I and ( li ) F are quantities related to tangents of the lightray trajectory at the emission and reception events, and they are given by Eqs. ( 196) to (197d) with η = 0 for ( li ) I and η = 1 for ( li ) F .In these formulas, we set ∆x i F = 0 and l0 = 0 corresponding to xα ( l) = (ct I , xi ( l)).The details of the expressions for ( li ) I and ( li ) F are given in the Theory section, but the corresponding contribution to two-way transfer is discussed further in this summary.The exponent I in the last term of Eq. ( 56) is given as where the fields n, α, A i , W, and their partial derivatives with respect to the space-time coordinates are evaluated in t = t I and x i = xi ( l) if not stated explicitly otherwise.Equation ( 56) together with Eqs. ( 196) to (197d) and Eq. ( 57) gives the frequency ratio including all terms on the order of P(3).The first fraction in Eq. ( 56) contains the gravitational redshift effect and the second-order Doppler effect, the second fraction contains the first-order Doppler effect, and the exponential term is the effect of nonstationarity.

Two-way frequency transfer
In two-way frequency transfer, a stationary observer on the ground A emits a signal with the proper frequency ν A1 from the position x i A1 at the time t A1 .The signal is received by observer B (ground or satellite) with the proper frequency ν B in the position x i B at the time t B and is immediately transponded back with the same frequency to observer A, where it is received with the proper frequency ν A2 in the position x i A2 at the time t A2 .The velocities of observers A and B at the particular emission and reception events are denoted as v i A1 , v i B , and v i A2 .For the stationary (fixed at the ground) observer A, the velocities v i A1 and v i A2 in the corotating coordinates and the corresponding spatial shift x i A2 − x i A1 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: The correction ∆, including all terms on the order of P(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 v i A1 /c, is given by Eqs. ( 212) to (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 x α A1 , x α B , or x α A2 of the emission and reception events are denoted by the corresponding index, for example, W A1 = W(t A1 , x j A1 ) and ).The base path xi + ( l) is again defined by the boundary points xi + (0) = x i A1 and xi + ( L) = x i B , and the fields denoted with the index + are evaluated in t = t B and x i = xi + ( l), for example, . The basis χ i + , ς i + , γ i + is associated with the path xi + ( l), and ε B is the angle ε F of the path xi + ( l).The acceleration of observer A at the time t A1 in the rotating frame is denoted a i A1 (i.e., a i A1 = d 2 x i A /dt 2 | t A1 ).Using this notation, we obtain where the atmospheric effects are included in the last line, in which and 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 n(r) (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.First, we evaluate the contribution of the spherical part of the atmospheric refractive index n(r) to the two-way frequency transfer correction (59), which is given by the first three lines of Eq. ( 60).We denote As an example, we consider an observer on the ground, located at the equator with r A = 6371 km and a satellite moving on a circular orbit in the equatorial plane at a radius of r B = r A +408 km in the direction of the rotation of the Earth.The corresponding satellite velocity in the co-rotating system is v B = 7.17 km/s.We consider the spherical refractive index field given by Eq. ( 51), which corresponds to the isothermal atmosphere in hydrostatic equilibrium with a constant air composition, and we use the same parameters as given in the initial part of Sect.4.5.1 up to Eq. ( 52).
The resulting corrections ∆ (1) At to ∆ (3) At and their sum ∆ (S ) At 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 ∆ (2)  At and ∆ (3) At were calculated by transforming the integration variable into r using Eq.(36).
Fig. 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 clock-on-satellite experiments.

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 To show how large the effect can be, we evaluate Eq. ( 63) for the simple situation with the satellite position x i B radially above the emission point on the ground x i A1 , in which case, the path xi + ( l) is just a straight vertical line connecting x i A1 with x i B (i.e., xi + ( l) = x i A1 + χ i + l and L = r B − r A ). Next, we assume that the wind field V j at t = t B is perpendicular to the path tangent χ i + in the neighborhood of xi + ( l) (horizontal flow).In this case, using Eq. ( 10) and writing θ (deg) Fig. 5: Three contributions of the spherically symmetric, static part of the atmospheric refractivity ∆ (1) At , ∆ (2) At , and ∆ (3) At to the two-way frequency transfer correction (59) and their sum ∆ (S )  At as functions of the satellite zenith angle θ.The example is evaluated for ground-to-satellite 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.
the first term on the RHS of Eq. ( 64) vanishes, and the second term, which equals −dA i+ /d l, can be integrated by parts to obtain where the term containing A i+ | L = A i B , 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 xi + ( l), the wind speed V i is constant and has a direction opposite to the satellite velocity v i B , which is itself perpendicular to χ 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 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 v B = 7.36 km/s, which approximately corresponds to the ISS, Eq. ( 65) gives This means that a horizontal wind speed V > 11 m/s 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 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 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 xi + ( l) in the direction perpendicular to both v i B and χ i + by the small Euclidean distance d.In this case, taking n i in the direction of v i B , Eq. ( 65) leads to However, this relation between the correction ∆ wind At 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).

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.

Light rays in flowing media
We assume that an electromagnetic signal propagates in spacetime with the metric g αβ , filled with a medium that moves with a four-velocity field U α satisfying U α U α = −1.Next, we assume that the electromagnetic properties of the medium are linear, isotropic, transparent, and nondispersive, such that they are described by two scalar functions, the permittivity ϵ(x α ), and the permeability µ(x α ).The refractive index of the medium is then defined as n(x α ) = c ϵ(x α )µ(x α ).
According to Gordon (1923), the Maxwell equations in the medium can be reformulated using the optical metric, which is defined as 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 ḡαµ ḡβν F µν = Fαβ .The Maxwell equations in the medium with vanishing free four-current then read We use the ansatz of geometrical optics, where is a wavelength related parameter, S (x α ) is the so-called eikonal function, and R denotes the real part.We define The Maxwell equations ( 73) then give to the order −1 Contracting the first equation with ḡγν k ν and using the second equation, we obtain therefore, kα is a null vector of the optical metric.At the same time, from Eq. ( 75), it follows that Contracting Eq. ( 78) with kα and using Eq. ( 77), we obtain which means that kα 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 kα .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 equation ( 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 kβ to be tangent to the trajectory x β (a) we search for, that is, we obtain the equation of motion with an affine parameterization where the inverse optical metric is given as

Frequency shift in nonstationary flowing media
The period of the electromagnetic wave given by Eq. ( 74) as observed by an observer moving along the trajectory x α (τ) parameterized by its proper time τ and having a four-velocity u α = dx α /dcτ is given as the proper time ∆τ during which the phase S (x α (τ))/ changes by 2π in a linear order (see, e.g., Chen & Kantowski (2008) and Synge (1960)), so that we have where the minus sign is there to obtain a positive ∆τ for future oriented kα .Therefore, the proper frequency ν = 1/∆τ is given as We consider a signal emitted with the proper frequency ν I by an observer with the four-velocity u α I from the spacetime point x α I that is received with the proper frequency ν F by an observer with the four-velocity u α F in the spacetime point x α F lying on the signal trajectory given by the solution x α (a) of the Eqs.( 80) and ( 81), satisfying x α (a I ) = x α I and x α (a F ) = x α F .The ratio of the frequencies then is 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 all along the light ray.We denote κ(λ) = dλ/da the factor converting l α to the affine tangent kα , that is, The frequency ratio then is Now we express the ratio κ I /κ F in terms of the fields along the signal trajectory.The geodesic equation ( 79) gives Contracting Eq. ( 90) with ξβ , we obtain where in the second line we used l α ∇α (l β ξβ ) = 0, which follows from Eq. ( 86).Taking into account that the Lie derivative of ḡαβ with respect to the field ξ α is given as we obtain 1 κ Integrating Eq. ( 93) along the signal trajectory from λ I to λ F , we obtain Equation ( 88) then leads to the following result for the frequency ratio: 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 l0 = −1. (96) Denoting with v i = dx i /dt the coordinate velocity of an observer, the four-velocity of the observer is expressed as and we obtain 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 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 and li = −k i /k 0 .The first fraction in Eq. ( 99) contains the gravitational redshift and the second-order Doppler shift.The second fraction is the first-order Doppler shift in the flowing medium.Here, effects of the medium such as bending of the light ray due to the refractive index gradient or dragging of the light by the motion of the medium enter most significantly.The integral in the exponent in Eq. ( 99) gives the correction for nonstationarity of the optical metric, or in other words, for the time dependence of the refractive index, the flow velocity of the medium and the gravitational field along the signal path.
In order to proceed with Eq. ( 99) that is used for the frequency transfer and also with the calculation of the signal propagation time used for the time transfer, we now need to find the solution of the system of Eqs. ( 80) and ( 81), reparameterized to the parameter λ or another related parameter.Using Eqs. ( 87) and ( 93), we transform the system of Eqs. ( 80) and (81) into the following set of equations of motion in Hamiltonian form: 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.

Optical metric in the vicinity of the Earth
The approximation level of the model as described in Sect.4.1 requires the Doppler terms li v i /c in Eq. ( 99) to be evaluated including all terms on the order of P(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 P( 2).This can be deduced retrospectively from the procedure leading to the final results.
We solve the problem in the rotating coordinates, where the components of the spacetime metric are given by Eq. ( 3) and the components of the four-velocity of the medium U α are where V i (x α ) is the velocity field of the medium in the rotating coordinates given by Eq. ( 9), and Using the definition (10), the contravariant components of the optical metric (82) in the rotating coordinates to the approximation level described above are The covariant components of the optical metric ( 72) in the rotating coordinates to the same approximation level are

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 L ξ ḡαβ = ∂ 0 ḡαβ .Expanding the RHSs of Eqs. ( 100) and ( 101) such that they include all terms on the order of P(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 where we used the fact that which can be obtained when we express lk = (δ i j li l j ) 1/2 a k with δ i j a i a j = 1, and we solve the null condition ḡαβ lα lβ = 0 for (δ i j li l j ) 1/2 .Equation ( 107) can be inverted to give li as a function of dx i /dλ.Thus, we obtain the following formula including all terms on the order of P(2): Differentiating Eq. ( 107) with respect to λ and using Eqs.( 106), ( 108), and ( 110) we arrive at the second-order Newton-type equation of motion for a light ray, with the RHS including all terms on the order of P(2), where the antisymmetrization is defined as 110) into Eq.( 108), we obtain the time component of the equation of motion, with the RHS including all terms on the order of P(2), For the purpose of analyzing the two-way transfer, it is convenient to change the parameter λ in the equation of motion into the Euclidean length of the light ray l E , which is defined by In this case, the spatial part of the light ray x i (l E ) has a unit tangent: Using l i = |l i |dx i /dl E and solving the null-vector condition ḡαβ l α l β = 0 together with Eq. ( 86) for l 0 and |l i |, we obtain the following formulas including all terms on the order of P(2): The transformed Eqs. ( 111) and ( 112) with RHSs including all terms on the order of P(2) then are and The terms in the first line of Eq. ( 117) correspond to optical acceleration due to the refractive index gradient, gravitational acceleration, centrifugal acceleration, and acceleration due to the time variation in the field A j .These accelerations are projected onto a plane perpendicular to the light-ray direction by the projection tensor at the end of the line.This corresponds to the selected parameterization by the Euclidean length, which does not allow changes in "velocity" magnitude (see Eq. ( 114)).The second line is an analog of a magnetic force, with v R j 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 v i R = 0 if needed.

Solving the equation of motion for a light ray
In this section, we solve the system of differential equations ( 117) and ( 118) assuming that we know the spatial coordinates of an emitter x i I at the event of the emission, the spatial coordinates x i F of the receiver at the event of reception, and either the coordinate time of the emission x 0 I = ct I or the coordinate time of the reception x 0 F = ct F .Denoting by L the Euclidean length of the light path between the emission and reception events, the range of the parameter l E is l E ∈ [0, L], and we assume the following boundary conditions: x 0 (0 For the sake of brevity, we write Eqs. ( 117) and ( 118) in the form where we introduced the coefficients 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 with We denote by xi ( l) the solution of the equation with l ∈ [0, L] satisfying the boundary conditions where we admit certain deviations ∆x i I , ∆x i F of the boundary points from the boundary points of the complete solution given by Eq. ( 119), which will be useful for two-way time and frequency transfer.Practically, we use ∆x i I = 0 and either , which is a spatial shift of the observer on the ground in corotating frame during the two-way back and forth propagation of the signal.
Since Eq. ( 132) is a special case of Eq. ( 117) with vanishing fields α, ∆W, v i R , and A i , the condition (114) also applies to the solution of Eq. ( 132), so that we have which means that xi ( l) is the path of a light ray between the points x i I − ∆x i I and x i F − ∆x i F for vanishing α, ∆W, v i R , and A i , parameterized by its Euclidean length, and L 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, L] by introducing a parameter l = ( L/L)l E on x i (l E ), the solution of the complete Eq. ( 121) can be expressed as which defines the correction ∆x i ( l).Similarly, if x 0 (l E ) is the solution of Eq. ( 122) for the time coordinate, the reparameterized solution can be expressed as where x0 = c t is a constant with t ∈ [t I , t F ] and ∆t( l) is the function to be found.The next step is to express the differential equation for the correction ∆x i ( l).To do this, we need to express the ratio L/ L in order to proceed with the derivatives d/dl E = ( L/L)d/d l.Taking a derivative of Eq. ( 135) with respect to l and using Eqs.( 114) and ( 134), we obtain Inserting Eqs. ( 135) and ( 136) into Eq.( 121), expanding the equation up to the linear order in ∆x α and d∆x i /d l and using Eqs.( 129) and ( 132) and the linear expansion of Eq. ( 137) in d∆x i /d l, we arrive at the equation of the form with the boundary conditions which follow from Eqs. ( 119), ( 133), and ( 135), and with the coefficients given by We proceed with the solution of the particular equations.The solution xi ( l) of the zeroth-order Eq. ( 132) with Ê j given by Eq. ( 130) has been discussed in Sect.4.3.We therefore continue with the solution of Eq. ( 138) for the spatial path correction ∆x i ( l) with the boundary conditions (139).

Solution for ∆x i ( l)
Integrating Eq. ( 138) twice, we obtain with c i 1 , c i 2 being integration constants, and [ ]( l′′ ) denoting that the function in the square brackets is evaluated in l′′ .Determining c i 1 , c i 2 from the boundary conditions (139), integrating by parts using and expressing the integrals over the range l′ ∈ [0, l] or l′ ∈ [ l, L] as integrals over the full range l′ ∈ [0, L] using the Heaviside step function H( l′ − l), which is defined as H( l′ − l) = 0 for l′ < l, H( l′ − l) = 1 for l′ > l and H( l′ − l) = 1/2 for l′ = l, we obtain Integrating by parts again, we can express the term containing d∆x j /d l in terms of ∆x j .Thus, we obtain the equation of the form with The functions have the following properties, which will be useful later: 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 ∆x i ( l), and then we can use some results of the linear operator theory in Banach spaces (see, e.g., Naylor & Sell (1971)).
We define the vector space of all continuous functions f i : [0, L] → R 3 , and we define the norm on this vector space as where f 1 , f 2 , f 3 are the corresponding component functions [0, L] → R, 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 one-dimensional case).The subspace of B 3 that is formed by functions satisfying 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 i j : B 3 → B 3 a bounded linear operator on B 3 , and by A i j f j the image of a vector f i ∈ B 3 .The operator norm ∥A i j ∥ of the operator is defined as where on the RHS, the norm of the Banach space is used, which is given by Eq. ( 155) in our case (for details, see definition 5.8.1 and lemma 5.8.2 in Naylor & Sell (1971)).Now we can proceed with Eq. ( 148).We assume that for any values of i, j the functions C i j and B i j − dC i j /d l that appear in Eq. ( 150) are continuous on [0, L], and we define the linear operators D i j : B 3 → B 3 0 and K i j : B 3 → B 3 0 as where the vanishing of the image value in l = 0 and l = L 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 Using these definitions and denoting by I i j the identity operator on B 3 , Eq. ( 148) can be written as According to theorem 6.7.1 in Naylor & Sell (1971), for a continuous linear operator K i j satisfying ∥K i j ∥ < 1, an inverse operator to I i j − K i j exists, and it is given by geometric series as where (K i j ) n is the operator given by n consecutive applications of K i j , so that (K i j ) 0 = I i j , (K i j ) 1 = K i j , 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 K i j can be proven.The assumption ∥K i j ∥ < 1 should be verified for the physical situation of interest.For the operator (158), the operator norm (156) can be expressed in terms of the function K i j ( l, l′ ) as Using Eq. ( 161), the solution of Eq. ( 160) can be written as where the boundary conditions ( 139) are ensured by the definition of b i and by the properties of the operators D i j , K i j given by Eqs. ( 151) and ( 152).For the derivative d∆x i /d l, which is also needed below, we obtain where and the operator Ki j is defined as with

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 we can write Inserting Eqs. ( 135) and (136) into Eq.( 122), using Eq. ( 137) for reparameterization to the parameter l, expanding the equation as Taylor series in ∆x α and d∆x i /d l and using Eqs.( 169), ( 170) and ( 171), we obtain the following equation for ∆t( l): with σ( l) and ∆σ( l) given by σ where the fields n, W, v Ri , and A i and their derivatives ∂ i are evaluated in xα ( l).On the RHS of Eq. ( 172), we kept all terms on the order of P(2), taking into account that the largest order of ∆t is c −1 (p c −1 = 1, p N = p α = 0).This is sufficient to reach the required approximation level for the signal propagation time, as described in Sect.4.1.All terms containing ∆x i t are on a negligible order.Now we denote l0 the value of l for which ∆t( l0 ) = 0, for example, if we select t = t I , we have l0 = 0, or if we select t = t F , we have l0 = L. Integrating Eq. ( 172) from l0 to l, we obtain c∆t( l) = 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: ∆t( l) ≈ ( l − l0 )/c, which we obtain by setting σ + ∆σ ≈ 1 and ∂ t n ≈ 0 in Eq. ( 175).Equation ( 175 In one-way time transfer, we wish to express the coordinate time t F − t I of the propagation of a signal that is emitted from a position x i I at a time t I and is received in a position x i F at a time t F .Based on Eqs. ( 120) and ( 136), we obtain t F − t I = ∆t( L) − ∆t(0).Using Eq. ( 176), this leads to Now we proceed with the integration of ∆σ.We set ∆x i I = 0, and we assume ∆x i F 0 in the boundary conditions (139).We keep nonvanishing ∆x i F as a preparation for later use in the twoway time transfer, where we set for the back-propagating signal ∆x i F 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 v i A1 being the velocity of observer A in the corotating frame in emission event A1.Thus, we consider ∆x i F to be on the order of c −1 .In the time-transfer formulas, however, we keep only the highest-order term containing ∆x i F , which is sufficient for an accuracy of 1 ps in case of a stationary ground clock for which the velocity v i A1 is given by the deformations of Earth, such as solid Earth tides, and by the nonuniformity of the rotation of Earth.In one-way transfer, we can set ∆x i F = 0 at the end.Using the properties ( 151) and (152) in Eq. (170), we obtain ∆x i t (0) = ∆x i t ( L) = 0, and therefore, the boundary conditions for ∆x i 0 are the same as for ∆x i , namely Integrating terms containing d∆x i 0 /d l in ∆σ by parts with use of Eq. ( 178), using Eqs.( 132) and ( 138) to express the second derivatives and using v Ri = ϵ i jk ω j x k , we obtain the following expression including all terms on the order of P(2) except the p T = 2 order terms containing ∆x where S i (1) is the first-order (p T = 1) approximation of S i , which is given as The ∆x j 0 in Eq. ( 179) can be expressed using Eq. ( 169).In this expansion, only the first term with S k 0 approximated by S k (1) is relevant for the required accuracy.Thus, we obtain where b j = ∆x j F l/ L does not contribute to the considered order in Eq. ( 179).Using the operator definition (157), we obtain the following result including all terms on the order of P(2) except the p T = 2 order terms containing ∆x In this formula, with S i (1) given by Eq. ( 180), the approximation ( 41) is sufficient to obtain the required accuracy for the resulting time-transfer formulas.In this case, using Eq. ( 149), some of the integrations can be performed explicitly.In the last term of Eq. ( 173), the approximation (41) is also sufficient.Therefore, using Eq. ( 177) with Eqs. ( 173) and (182), we obtain the following formula for the propagation time, including all terms on the order of P(3) except the p T = 3 order terms containing ∆x i F : where the fields n, α, W, v Ri , and A i and their partial derivatives with respect to the spacetime coordinates are evaluated in xα ( l) if not stated explicitly otherwise.Selecting ∆x i F = 0 and l0 = 0, which corresponds to t = t I , approximating the terms in the third line by integration over a straight line connecting xi I and xi F instead of xi ( l), which leads to a negligible difference in terrestrial conditions, and approximating L = D in the fourth line (see Eq. ( 42)), we obtain Eq. ( 43) in the summary.

Definitions in two-way transfer
In two-way time transfer, and later also in two-way frequency transfer, we consider a signal that is emitted by a stationary observer on the ground A from the position x i A1 at the coordinate time t A1 to an observer B, who receives the signal in the position x i B at a coordinate time t B and immediately sends it back to A, who receives it in the position x i A2 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 From the boundary conditions for xi + ( l) and xi − ( l) described above, it follows that ∆x The shift in the position of A during the back and forth propagation of the signal we denote which leads to l0+ = L and l0− = 0.
Denoting by x i A (t) 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 v i A1 = dx i A /dt| t A1 and expressing t A2 − t A1 = ∆t + + ∆t − using Eq. ( 183), we obtain ∆x i A up to the highest order, Transforming the integration variables in ∆t which follows from Eqs. ( 184) and ( 185) and using the properties ( 153) and ( 154) of the function D( l, l′ ) we arrive at the following result for the two-way time transfer correction including all terms on the order of P(3) except the p T = 3 order terms originating from ∆x i A : where the fields with the index + are evaluated in xα + ( l), for example, v Ri+ = v Ri ( xα + ( l)), and χ i + is a unit vector in the direction from x i A1 to x i B .This is Eq. ( 46) in the summary.

One-way frequency transfer
In this section, we proceed with the calculation of the frequency shift of the signal that propagates in the atmosphere of Earth according to Eq. ( 99).An observer moving with the velocity v i = dx i /dt in the corotating coordinates has the four-velocity components (in the same coordinates) given by Eq. ( 97) with Denoting by , with x α I , x α F being the spacetime points of the signal emission and reception, respectively, and denoting by v i I , v i F the observers' velocities at the emission and reception events, we can express the first part of Eq. (99) as Next, we proceed with the Doppler part of Eq. ( 99).The components li 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 P(2): 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 d∆x i /d l as where, at the end, the terms in the third line do not contribute to the required order of li .In order to obtain the boundary values ( li ) I and ( li ) F , we need to express the boundary values of d∆x i /d l.
We introduce the parameter η = 0 or 1, and we evaluate Eq. ( 165) in η L. We assume ∆x i I = 0, and we keep nonvanishing ∆x i F as a preparation for the two-way frequency transfer, similarly as in the case of the two-way time transfer.Denoting we obtain We denote by li|η the boundary values of li , namely li|0 = ( li ) I and li|1 = ( li ) F .We write the resulting expression for li|η as the sum including all terms on the order of P(2) except terms containing ∆x i F , 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 ∆x i F is given separately.These terms are given by Eqs.(199) in the third and fourth line of Eq. ( 214) are given by Eqs. ( 198)-( 201) with xα ( l) = xα + ( l).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 with β i(1) given by Eq. ( 213) and β i(2) given by Eq. ( 60).Applying Eq. ( 216) in Eq. ( 212) leads to Eq. ( 59) in the summary section.Some steps in the calculation of β i up to the second order are described below.For the end-point tangent, we use Eq. ( 40) expanded up to the first order, d xi where the index B replaces the general index for a final point F.
Where possible, the approximation ( 41 where the second derivative of the path xi + ( l) 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. ( 198)-( 201) and ( 157 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 n+ .The use of these relations leads to the cancellation of several terms in the expression for β i(2) .

Conclusion
A relativistic model of one-way and two-way time and frequency transfer through flowing media was developed in this paper.The main focus was on applications to the atmosphere of Earth.The model includes gravitational and atmospheric effects given by the fields of the scalar gravitational potential W(x α ), the atmospheric refractive index n(x α ), and the wind speed V 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 xi ( l) in a spherically symmetric, static part of n(x α ) and W(x α ) is found, which is parameterized by its Euclidean length in the corotating frame l.This solution is defined by the spatial coordinates of the emission and the reception events as boundary values.Then, a complete solution x i ( l), t( l) is found using a perturbation method that takes the inertial forces, wind speed, and deviations of n(x α ) and W(x α ) from sphericity and staticity into account.The time-and frequency-transfer corrections are then expressed in terms of integrals of the fields along the path xi ( l) at a certain constant coordinate time.
In the numerical examples evaluated in this paper, we focused on the two-way ground-to-satellite transfer with a satellite altitude similar to that of the ISS.In the two-way time transfer, the main contribution of the atmosphere appears in the Sagnac effect and is given by a change in the Sagnac area that is spanned by the light path xi ( l) when the path changes from the vacuum case to the case influenced by the atmospheric refraction.In the example studied in this paper, the contribution of the atmosphere to the Sagnac effect is below 0.1 ps for the satellite zenith angle θ (the angle between the vertical line from a ground station and the line connecting the ground station with the satellite position at the reception/re-emission time) below 78 • , it reaches 1 ps for θ ≈ 87 • , and it increases to approximately 5 ps as θ approaches 90 • .Therefore, for time transfer at an accuracy level of 1 ps, this effect is only significant for large zenith angles.
The effect of the wind in the two-way time transfer is given by the difference in the propagation times caused by the Fresnel-Fizeau effect of light dragging.It is much smaller than 1 ps for normal atmospheric conditions.
In the two-way frequency transfer, the correction ∆ of Ashby (1998) and Blanchet et al. (2001) is generalized to include the atmospheric effects.In the example presented in this paper, the effect of the spherically symmetric, static part of the refractive index gives a contribution that starts at 10 −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 two-way frequency transfer, it may lead to a significant contribution to the correction ∆.
For the example of a constant horizontal wind field, the effect reaches 10 −17 already for a wind speed of about 11 m/s (the wind speed in the corotating frame is meant).
The effects of the deviation of n(x α ) from sphericity and staticity are also included in the resulting corrections.They are not evaluated numerically in this paper, however.
The evaluation of the atmospheric effects in the time-and frequency-transfer corrections including the effect of the wind shows that they need to be considered in the data analysis of the forthcoming clock-on-satellite experiments such as ACES or I-SOC.The results obtained in this paper can be used not only to calculate the time-and frequency-transfer corrections from given atmospheric data, but also inversely to determine atmospheric properties such as the refractive index profile or the flow rate from the time or frequency measurements.This will be a subject of a future research.
) which depends on the spatial coordinates through the radius r = |x i | only and introducing a new parameter ζ with a function dependence ζ( F and xi I , the angle θ I between xi F − xi I and xi I , and the angle θ F between xi F − xi I and xi F , and we denote D = | xi F − xi I | (see Fig.

Fig. 3 :
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 xi + ( l) changes from the case including atmospheric refraction and gravitation to a vacuum case that only includes gravitation (n = 1 in Eq. (13)).

Fig. 4 :
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 ground-to-satellite transfer in the equatorial plane, with r A = 6371 km, r B − r A = 408 km, and χ i + inclined against the rotation of the Earth.