Open Access
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/0004-6361/202345994
Published online 23 May 2023

© The Authors 2023

Licence Creative CommonsOpen 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 free-space optical links (see Bodine et al. 2020 and references therein, and also Djerroud et al. 2010, Gozzard et al. 2018, Kang et al. 2019, Dix-Matthews et al. 2021 and Shen et al. 2021) enable increasing accuracies of the experiments. Therefore, the atmospheric effects originating from spatial and temporal variations of the refractive index of air or from winds start to play a more important role.

In case of the ground-to-satellite time and frequency transfer, two projects are going to be realized in the near future that will place an atomic clock on the International Space Station (ISS), namely the Atomic Clock Ensemble in Space (ACES), with a Cs clock and H-maser (Cacciapuoti et al. 2020; Lilley et al. 2021), and the Space Optical Clock on the ISS (I-SOC), with an optical lattice clock (Cacciapuoti & Schiller 2017; Origlia et al. 2018). In the ACES mission, the absolute frequency accuracy of the on-board clock is about 10−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 100ps. For the ELT, the overall planned accuracy of the time transfer is 50 ps (Turyshev et al. 2016). The I-SOC proposal uses an improved version of the MWL as well as the optical link of ELT, aiming at an accuracy of a few ps in the clock synchronization (Cacciapuoti & Schiller 2017). Therefore, based on the planned accuracies of these experiments, a relativistic model of one- and two-way time and frequency transfer is needed that includes the atmospheric effects and has an accuracy of approximately 1 ps for the time transfer and a relative accuracy of the lower multiples of 10−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 J2 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. (2014a,b) employed this formalism to compute various observables relevant for actual space missions in the Solar System, including the coordinate propagation time and the frequency shift of light.

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

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

The aim of this paper is to develop a relativistic model of one- and two-way time and frequency transfer that includes gravitational and atmospheric effects for realistic fields of the gravitational potential, the refractive index, and the wind speed, assuming neither stationarity nor spherical or axial symmetry. The intended accuracy of the model is given by the experimental needs described above: approximately 1 ps for time transfer, and lower multiples of 10−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.

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 ai is denoted as , with δij being the Kronecker delta. Quantities ωi, , υi, Vi, Ai, hi i, ςi, and γi in this paper are treated as Euclidean three-vectors when their index is lowered or raised, for example, Vi = δijVj The three-dimensional antisymmetric Levi-Civita symbol is denoted as ϵijk, 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 two-way corrections more straightforward.

The rotating coordinate system xα = (x0, xi) is related to the GCRS system as , , where Rij(t) ∊ SO(3), or inversely, , where . The time-dependent rotation matrix Rij(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 xi. 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 time-dependent components as ∆W(xα) such that (4) (5)

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

3.2 Medium

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

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

The atmospheric temperature and pressure as functions of altitude show a systematic behavior at large scales that is discussed, for example, in ICAO (1993) and NOAA (1976; see also Appendix A.2). Spatial and temporal local fluctuations are present as well. Therefore, it is reasonable to decompose the atmospheric refractive index field as (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 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 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 Vi(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 xi (t) of an observer in the rotating coordinates, we denote υi = dxi/dt the observer’s velocity in these coordinates.

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

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

4.1 Approximation level of the model

As mentioned in the introduction, a model accuracy of approximately 1 ps for the time transfer and relative accuracy of lower multiples of 10−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/c2 (higher-order gravitational effects such as the Lense-Thirring 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 Vi /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−1t 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 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 . 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 two-way 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 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 pTq only. Expressions containing terms with and only are denoted O(ck) 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 first-order 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 ∂ 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 two-way 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 .

Defining a field w as (13)

which depends on the spatial coordinates through the radius r = ∣xi∣ 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 hi being a constant vector. Equation (17) shows that for ∣hi 0, the signal propagates in a plane that intersects the origin r = 0 and has a normal aligned with hi (for ∣hi∣ = 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 = ∣hi∣. 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 right-handed Cartesian coordinate system yi with the origin at r = 0, the y1-axis in direction, and the y3-axis in the hi direction (for h = 0, we can take any direction perpendicular to ), defining polar coordinates by y1 = r cos φ and y2 = 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 = rIwI sin ψI. In the boundary value problem, we need to express h as a function of rI, rF 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 right-hand 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 = rFwF sin θF, corresponding to εF = 0 in Eq. (26), (b) we calculate εT(h) according to Eq. (32), (c) we calculate εF according to Eq. (29), (d) we calculate h according to Eq. (26) and start another loop from (b). As a result, not only the constant h is obtained, but also the angles εF, εT, and εI = εFεT. Alternatively, the procedure can be performed using the angle εI given by Eq. (28) instead of εF. We also note that for the ground-to-satellite or satellite-to-ground transfer, the εF/I angle at the satellite position is at least one order of magnitude smaller than the corresponding angle at the ground position. Therefore, neglecting the εF/I angle at the satellite position gives a better initial estimate of h according to Eq. (26) than neglecting the corresponding angle at the ground position.

To retrieve the solution in the coordinates xi, we express the components of the first two basis vectors of yi in the coordinate system xi. 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 right-handed basis (χi, ςi, γi) with χi pointing in direction from to , γi pointing in direction of the normal to the propagation plane hi (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 lowest-order 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)

Next we obtain (42)

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.

thumbnail Fig. 1

Definition of angles, basis vectors, and other quantities.

4.4 Coordinate time transfer

4.4.1 One-way time transfer

In one-way time transfer, we wish to express a coordinate time tFtI of the propagation of the signal emitted from a position at a time tI and received in a position at a time tF. 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, Ai and their partial derivatives with respect to the space-time coordinates are evaluated in t = tI 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 Vi in Ai 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 into the radial coordinate r using Eq. (36) and by expressing the integrated functions in spherical coordinates with the base path given by Eq. (21).

4.4.2 Two-way time transfer

In two-way transfer, we consider the signal emitted from a stationary observer on the ground A at the coordinate time tA1 , which corresponds to the proper time τA(tA1) 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 tB, which corresponds to the proper time τB(tB) of clock B (see Fig. 2). The spatial coordinates of this event are denoted . The pseudo-time-of-flight τB(tB) − τA(tA1) is obtained from the measurements. Then a signal is sent from observer B at the coordinate time tB and is received by observer A at the coordinate time tA2, which corresponds to the proper time τA(tA2) 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 one-way signal is received by observer B. We call this set-up the Λ-configuration. The pseudo-time-of-flight τA(tA2) − τB(tB) is also obtained from measurements.

Using the coordinate-time synchronization convention, the desynchronization of the clocks is defined as the difference τB(tB) − τA(tB), where τA(tB) is the proper time measured by observer A at the coordinate time tB. In case of the two-way time transfer, it can be expressed as (44)

where ∆τ+ = τA(tB) − τA(tA1) and ∆τ = τA(tA2) − τA(tB). The first two lines of Eq. (44) are the measured pseudo-times-of-flight, and the difference ∆τ − ∆τ+ needs to be computed.

Denoting ∆t+ = tBtA1 the coordinate time of the signal propagation from A to B and ∆t = tA2tB 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 two-way time transfer correction, including all terms on the order of 𝒫(3) except the pT = 3 order terms originating from (46)

where the fields υRi, Ai, ∂iα, and ∂tα with the index + are evaluated in t = tB 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 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 (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 = hi/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.

thumbnail Fig. 2

Scheme of the two-way time transfer. The signal leaves observer A at a coordinate time tA1, corresponding to the proper time τA(tA1) of clock A, is reflected from observer B at a coordinate time tB, corresponding to the proper time τB(tB) of clock B and τA(tB) of clock A, and is finally received by observer A at the coordinate time tA2, corresponding to the proper time τA(tA2) of clock A.

thumbnail 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 ground-to-satellite transfer through an isothermal atmosphere at a temperature T0 with an unchanging composition corresponding to a molar mass of air Ma. 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/rA, R is the universal gas constant, and NA 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 rA = 6371 km, corresponding to the surface of Earth, rB = rA + 408 km, approximately corresponding to the altitude of the ISS, and T0 = 288.15 K (15 °C) and Ma = 28.964 × 10−3 kg mol−1, corresponding to dry air with a carbon dioxide molar fraction xc = 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 two-way transfer, we denote the θI angle of the 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 (NA given by Eq. (52)) minus the Sagnac correction in vacuum (NA = 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)

thumbnail 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 rA = 6371 km, rBrA = 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 two-way 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 Vi (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 m2 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 inter-ferometric methods, as was shown, for example, by Tselikov & Blake (1998), who used the effect in a different context for flow metering applications in laboratory conditions.

4.6 Frequency transfer

4.6.1 One-way frequency transfer

In one-way frequency transfer, the goal is to express the ratio of the proper frequency ν¡ of a signal emitted by an observer with velocity from a position at a time tI and a proper frequency νF of the same signal as is received by the observer with a velocity in the position at the time tF. 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 light-ray 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 two-way 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, α, Ai, W, and their partial derivatives with respect to the space-time coordinates are evaluated in t = tI 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 second-order Doppler effect, the second fraction contains the first-order Doppler effect, and the exponential term is the effect of nonstationarity.

4.6.2 Two-way frequency transfer

In two-way frequency transfer, a stationary observer on the ground A emits a signal with the proper frequency νA1 from the position at the time tA1. The signal is received by observer B (ground or satellite) with the proper frequency νB in the position at the time tB 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 tA2. 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 pT ≥ 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 (pT > 3) and terms proportional to c−2NB, which are negligible for satellite applications because the atmospheric refractivity NB 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 = tB 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 tA1 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)

and (61)

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

The examples in the next section focus on the numerical evaluation of the atmospheric effects originating from the spherically symmetric, static part of the refractive index (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 two-way 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 rA = 6371 km and a satellite moving on a circular orbit in the equatorial plane at a radius of rB = rA+408 km in the direction of the rotation of the Earth. The corresponding satellite velocity in the co-rotating system is υB = 7.17 km s−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 clock-on-satellite experiments.

4.7.2 Effect of the wind

In this section, we evaluate the effect of the wind, which enters the correction (59) through the fourth line of Eq. (60). In case of a nonstationary Ai, 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 Vj at t = tB 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−2NB, 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 Vi is constant and has a direction opposite to the satellite velocity , which is itself perpendicular to (i.e., with and V = ∣Vi∣). Next, we consider the refractive index field approximated by Eq. (51). For the position of the observer on the ground, rA = 6371 km, and the satellite position rB = rA + 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 N0, ρ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 ni 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 ni 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).

thumbnail Fig. 5

Three contributions of the spherically symmetric, static part of the atmospheric refractivity , and to the two-way frequency transfer correction (59) and their sum 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 rA = 6371 km and a satellite moving on a circular orbit in the equatorial plane at a radius of rA + 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 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 .

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 four-current then read (73)

We use the ansatz of geometrical optics, (74)

where is a wavelength related parameter, S (xα) is the so-called 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 four-velocity uα = dxα/d 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 four-velocity from the spacetime point that is received with the proper frequency νF by an observer with the four-velocity 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)

or (90)

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)

we obtain (93)

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 (x0, xi) with x0 = ct on spacetime, we take ξα = (∂0)α). In this case, the condition (86) leads to (96)

Denoting with υi = dxi /dt the coordinate velocity of an observer, the four-velocity of the observer is expressed as (97)

and we obtain (98)

Moreover, expressing all tensor components in the coordinate basis of the selected coordinate system, the Lie derivative in the integral reduces to the partial derivative with respect to x0. 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 non-stationary medium. We also note that for ξα = (∂0)α, we have κ = −k0 .

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

In order to proceed with Eq. (99) that is used for the frequency transfer and also with the calculation of the signal propagation time used for the time transfer, we now need to find the solution of the system of Eqs. (80) and (81), reparameterized to the parameter λ or another related parameter. Using Eqs. (87) and (93), we transform the system of Eqs. (80) and (81) into the following set of equations of motion in Hamiltonian form: (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 = Δx0/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 four-velocity of the medium Uα are (102)

where Vi(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 ∂0Ai 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 δijaiaj = 1, and we solve the null condition for .

Equation (107) can be inverted to give ; as a function of dxi /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 second-order Newton-type 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 two-way transfer, it is convenient to change the parameter λ in the equation of motion into the Euclidean length of the light ray lE, which is defined by (113)

In this case, the spatial part of the light ray xi(lE) has a unit tangent: (114)

Using li = ∣li∣dxi/dlE and solving the null-vector condition together with Eq. (86) for l0 and ∣li∣, 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)

and (118)

The terms in the first line of Eq. (117) correspond to optical acceleration due to the refractive index gradient, gravitational acceleration, centrifugal acceleration, and acceleration due to the time variation in the field Aj. These accelerations are projected onto a plane perpendicular to the light-ray direction by the projection tensor at the end of the line. This corresponds to the selected parameterization by the Euclidean length, which does not allow changes in “velocity” magnitude (see Eq. (114)). The second line is an analog of a magnetic force, with υRj being the magnetic potential, leading to Coriolis acceleration, and Ajbeing 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 lE is lE ∊ [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)

with (130) (131)

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 two-way 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 two-way back and forth propagation of the signal.

Since Eq. (132) is a special case of Eq. (117) with vanishing fields α, ∆W, , and Ai, 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 Ai, parameterized by its Euclidean length, and is the total Euclidean length of the path, which can differ from L. Reparameterizing the solution xi(lE) of Eq. (121) to the range [0, ] by introducing a parameter on xi(lE), the solution of the complete Eq. (121) can be expressed as (135)

which defines the correction . Similarly, if x0(lE) is the solution of Eq. (122) for the time coordinate, the reparameterized solution can be expressed as (136)

where is a constant with ∊ [tI, tF] 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 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 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 ∆xj. Thus, we obtain the equation of the form (148)

with (149) (150)

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 f1, f2, f3 are the corresponding component functions , which are all continuous. This vector space together with the norm (155) forms a Banach space that we denote B3 (see, e.g., Naylor & Sell 1971 for the one-dimensional case). The subspace of B3 that is formed by functions satisfying is denoted .

The norm on the Banach space B3 induces an operator norm on the linear operators from B3 to itself. We denote by a bounded linear operator on B3, and by the image of a vector fi ∊ B3. 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 bi the linear function satisfying the boundary conditions (139). This function is given as (159)

Using these definitions and denoting by the identity operator on B3, 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 bi 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)

where (166)

and the operator is defined as (167)

with (168)

5.5.3 Equation for ∆t and its solution

For further analysis, it is convenient to split the solution (163) for ∆xi into a part that does not depend on ∆t and a part that does depend on it, using Eq. (140). Defining (169) (170)

we can write (171)

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)

with and given by (173) (174)

where the fields n, W, υRi, and Ai 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 ∆xi. 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 ∂tn ≈ 0 in Eq. (175). Equation (175) then gives (176)

5.6 Coordinate time transfer

5.6.1 One-way time transfer

In one-way time transfer, we wish to express the coordinate time tFtI of the propagation of a signal that is emitted from a position at a time tI and is received in a position at a time tF. 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 two-way time transfer, where we set for the back-propagating 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 time-transfer formulas, however, we keep only the highest-order 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 one-way 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ωjxk, we obtain the following expression including all terms on the order of 𝒫(2) except the pT = 2 order terms containing , (179)

where is the first-order (pT = 1) approximation of Si, 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 pT = 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 time-transfer formulas. In this case, using Eq. (149), some of the integrations can be performed explicitly. In the last term of Eq. (173), the approximation (41) is also sufficient. Therefore, using Eq. (177) with Eqs. (173) and (182), we obtain the following formula for the propagation time, including all terms on the order of 𝒫(3) except the pT = 3 order terms containing : (183)

where the fields n, α, W, υRi, and Ai; 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 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 at the coordinate time tA1 to an observer B, who receives the signal in the position at a coordinate time tB and immediately sends it back to A, who receives it in the position at the coordinate time tA2. 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 tI+ = tA1, , tF+ = tI = tB, , tF = tA2, 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 Two-way time transfer

Now we can express the time Δt+ := tF+tI+ by setting and in Eq. (183) and the time Δt := tFtI 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 tA1, denoting by and expressing tA2tA1 = Δ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 two-way time transfer correction including all terms on the order of 𝒫(3) except the pT = 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 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 υi = dxi/dt in the corotating coordinates has the four-velocity 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 dxi/dlE 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 two-way frequency transfer, similarly as in the case of the two-way time transfer. Denoting (194)

we obtain (195)

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, Ai, 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 Si 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 lE ∊ [0, L], and using Eqs. (105), (115), (116) and (118), we obtain the following expression including all terms (202)

where the fields ∂tn, ∂tAi, and ∂tW are evaluated along the signal trajectory xα(lE).

Next, we transform the integration variable to l using the first-order 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, α, Ai, 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 Two-way frequency transfer

In two-way frequency transfer, we consider the experimental setup as already described in Sect. 5.6.2, and we admit the definitions introduced there. Furthermore, we define 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, u0 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 light-ray 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 u0 terms in Eqs. (207) and (208) are expressed using Eq. (190). Quantities located at event A2 (except the potential WA2) 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 pT ≥ 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 first-order (pT = 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 (pT ≤ 2), leading to the contribution up to the third order (pT ≤ 3) in ∆. Moreover, the terms in Eq. (214) that are proportional to c−1 NB are negligible for satellite applications because the atmospheric refractivity NB 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 NB, 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 end-point 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 pT = 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 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 Vi (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, xi) 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 frequency-transfer 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 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 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 Ai = (1 − n2)Vi 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−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 frequency-transfer corrections including the effect of the wind shows that they need to be considered in the data analysis of the forthcoming clock-on-satellite experiments such as ACES or I-SOC. The results obtained in this paper can be used not only to calculate the time- and frequency-transfer corrections from given atmospheric data, but also inversely to determine atmospheric properties such as the refractive index profile or the flow rate from the time or frequency measurements. This will be a subject of a future research.

Acknowledgements

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

Appendix A Refractive index profiles

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

According to Eq. (4) of Owens (1967) the refractive index of air n depends on the densities of its components as (A.1)

where ρ1, ρ2, and ρ3 are the partial densities of dry, CO2 -free air, of water vapor, and of carbon dioxide in the mixture, and R1, R2, and R3 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 Nn − 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 xc of CO2, ρws is the density of pure water vapor at 20 °C and 1333 Pa and Naxs, and Nws 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 xc of CO2 and xw 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α), xc(xα), and xw(xα) are known.

For the special case when the air composition does not change with spacetime position (i.e., xc 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 N0 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)

with (A.5)

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 Ma 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 sea-level value of ϕ, and g is the standard acceleration due to gravity (constant). We note that H has the dimension of length. We assume that the atmosphere can be divided into intervals in H in which the temperature depends on H approximately linearly as (A.8)

where Hb is the boundary value of H between two intervals, Tb 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 frequency-transfer model developed in this paper uses a spherically symmetric, static component of the refractive index n(r) (see Eq. (6)) or the related refractivity . 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 T0. From Eqs. (A.10) and (A.3), we then obtain (A.11) (A.12)

where ρ0, N0, 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

  1. Alonso, I., Alpigiani, C., Altschul, B., et al. 2022, EPJ Quantum Technol., 9, 30 [CrossRef] [Google Scholar]
  2. Ashby, N. 1998, in IEEE International Frequency Control Symposium, 320 [Google Scholar]
  3. Belmonte, A., Taylor, M. T., Hollberg, L., & Kahn, J. M. 2017, Opt. Express, 25, 15676 [NASA ADS] [CrossRef] [Google Scholar]
  4. Beloy, K., Bodine, M. I., Bothwell, T., et al. 2021, Nature, 591, 564 [NASA ADS] [CrossRef] [Google Scholar]
  5. Blanchet, L., Salomon, C., Teyssandier, P., & Wolf, P. 2001, A&A, 370, 320 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bodine, M. I., Deschenes, J.-D., Khader, I. H., et al. 2020, Phys. Rev. Res., 2, 033395 [NASA ADS] [CrossRef] [Google Scholar]
  7. Born, M., & Wolf, E. 1999, Principles of Optics, 7th edn. (Cambridge University Press) [CrossRef] [Google Scholar]
  8. Bourgoin, A. 2020, Phys. Rev. D, 101, 064035 [Google Scholar]
  9. Bourgoin, A., Zannoni, M., & Tortora, P. 2019, A&A, 624, A41 [EDP Sciences] [Google Scholar]
  10. Bourgoin, A., Zannoni, M., Gomez Casajus, L., Tortora, P., & Teyssandier, P. 2021, A&A, 648, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Cacciapuoti, L., & Schiller, S. 2017, I-SOC Scientific Requirements, SCI-ESA- HRE-ESR-ISOC document [Google Scholar]
  12. Cacciapuoti, L., Armano, M., Much, R., et al. 2020, Eur. Phys. J. D, 74, 164 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cao, J., Zhang, P., Shang, J., et al. 2017, Appl. Phys. B, 123, 112 [CrossRef] [Google Scholar]
  14. Chen, B., & Kantowski, R. 2008, Phys. Rev. D, 78, 044040 [Google Scholar]
  15. Ciddor, P. E. 1996, Appl. Opt., 35, 1566 [Google Scholar]
  16. Delva, P., Hees, A., & Wolf, P. 2017, Space Sci. Rev., 212, 1385 [NASA ADS] [CrossRef] [Google Scholar]
  17. Delva, P., Puchades, N., Schönemann, E., et al. 2018, Phys. Rev. Lett., 121, 231101 [NASA ADS] [CrossRef] [Google Scholar]
  18. 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]
  19. Derevianko, A., Gibble, K., Hollberg, L., et al. 2022, Quantum Sci. Technol., 7, 044002 [NASA ADS] [CrossRef] [Google Scholar]
  20. Dix-Matthews, B. P., Schediwy, S. W., Gozzard, D. R., et al. 2021, Nat. Commun., 12, 515 [NASA ADS] [CrossRef] [Google Scholar]
  21. Djerroud, K., Acef, O., Clairon, A., et al. 2010, Opt. Letters, 35, 1479 [NASA ADS] [CrossRef] [Google Scholar]
  22. Duchayne, L., Mercier, F., & Wolf, P. 2009, A&A, 504, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Feng, J. C., Hejda, F., & Carloni, S. 2022, Phys. Rev. D, 106, 044034 [NASA ADS] [CrossRef] [Google Scholar]
  24. Fujieda, M., Piester, D., Gotoh, T., et al. 2014, Metrologia, 51, 253 [NASA ADS] [CrossRef] [Google Scholar]
  25. Gordon, W. 1923, Ann. Phys., 72, 421 [NASA ADS] [CrossRef] [Google Scholar]
  26. Gozzard, D. R., Schediwy, S. W., Stone, B., Messineo, M., & Tobar, M. 2018, Phys. Rev. Appl., 10, 024046 [CrossRef] [Google Scholar]
  27. Hachisu, H., Fujieda, M., Nagano, S., et al. 2014, Opt. Lett., 39, 4072 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hees, A., Bertone, S., & Le Poncin-Lafitte, C. 2014a, Phys. Rev. D, 89, 064045 [Google Scholar]
  29. Hees, A., Bertone, S., & Le Poncin-Lafitte, C. 2014b, Phys. Rev. D, 90, 084020 [Google Scholar]
  30. ICAO. 1993, Manual of the ICAO Standard Atmosphere extended to 80 kilometres (262 500 feet), 3ed edn., Doc 7488/3 [Google Scholar]
  31. Kang, H. J., Yang, J., Chun, B. J., et al. 2019, Nat. Commun., 10, 4438 [NASA ADS] [CrossRef] [Google Scholar]
  32. Koller, S. B., Grotti, J., Vogt, S., et al. 2017, Phys. Rev. Lett., 118, 073601 [NASA ADS] [CrossRef] [Google Scholar]
  33. Kopeikin, S. M., & Han, W.-B. 2015, J. Geod., 89, 829 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kopeikin, S., & Mashhoon, B. 2002, Phys. Rev. D, 65, 064025 [Google Scholar]
  35. Kopeikin, S. M., & Schäfer, G. 1999, Phys. Rev. D, 60, 124002 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kopeikin, S., Efroimsky, M., & Kaplan, G. 2011, Relativistic Celestial Mechanics of the Solar System (Wiley-VCH Verlag GmbH &Co. KGaA) [CrossRef] [Google Scholar]
  37. Landau, L. D., & Lifshitz, E. M. 1976, Mechanics, 3rd edn. (Butterworth-Heinemann) [Google Scholar]
  38. Le Poncin-Lafitte, C., Linet, B., & Teyssandier, P. 2004, Class. Quant. Grav., 21, 4463 [CrossRef] [Google Scholar]
  39. Leonhardt, U., & Piwnicki, P. 1999, Phys. Rev. A, 60, 4301 [NASA ADS] [CrossRef] [Google Scholar]
  40. Lilley, M., Savalle, E., Angonin, M. C., et al. 2021, GPS Solutions, 25, 34 [CrossRef] [Google Scholar]
  41. Linet, B., & Teyssandier, P. 2002, Phys. Rev. D, 66, 024045 [Google Scholar]
  42. Lisdat, C., Dörscher, S., Nosske, I., & Sterr, U. 2021, Phys. Rev. Res., 3, L042036 [NASA ADS] [CrossRef] [Google Scholar]
  43. Ludlow, A. D., Boyd, M. M., Ye, J., Peik, E., & Schmidt, P. O. 2015, Rev. Mod. Phys., 87, 637 [NASA ADS] [CrossRef] [Google Scholar]
  44. McGrew, W. F., Zhang, X., Fasano, R. J., et al. 2018, Nature, 564, 87 [NASA ADS] [CrossRef] [Google Scholar]
  45. Mehlstäubler, T. E., Grosche, G., Lisdat, C., Schmidt, P. O., & Denker, H. 2018, Rep. Prog. Phys., 81, 064401 [CrossRef] [Google Scholar]
  46. Müller, J., Dirkx, D., Kopeikin, S. M., et al. 2018, Space Sci. Rev., 214, 5 [CrossRef] [Google Scholar]
  47. Naylor, A. W., & Sell, G. R. 1971, Linear Operator Theory in Engineering and Science (Holt, Rinehart and Wilson, Inc.) [Google Scholar]
  48. NOAA. 1976, U. S. Standard Atmosphere, Report No. NOAA-S/T 76-1562 [Google Scholar]
  49. Origlia, S., Pramod, M. S., Schiller, S., et al. 2018, Phys. Rev. A, 98, 053443 [NASA ADS] [CrossRef] [Google Scholar]
  50. Owens, J. C. 1967, Appl. Opt., 6, 51 [NASA ADS] [CrossRef] [Google Scholar]
  51. Perlick, V. 2000, Ray Optics, Fermat’s principle, and Applications to General Relativity (Springer-Verlag Berlin Heidelberg) [Google Scholar]
  52. 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]
  53. Pizzocaro, M., Sekido, M., Takefuji, K., et al. 2021, Nat. Phys., 17, 223 [NASA ADS] [CrossRef] [Google Scholar]
  54. Riedel, F., Al-Masoudi, A., Benkler, E., et al. 2020, Metrologia, 57, 045005 [NASA ADS] [CrossRef] [Google Scholar]
  55. Robert, C., Conan, J.-M., & Wolf, P. 2016, Phys. Rev. A, 93, 033860 [Google Scholar]
  56. Safronova, M. S., Budker, D., DeMille, D., et al. 2018, Rev. Mod. Phys., 90, 025008 [CrossRef] [Google Scholar]
  57. Schreiber, K. U., Procházka, I., Lauber, P., et al. 2010, IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 57, 728 [Google Scholar]
  58. Shen, Q., Guan, J.-Y., Zeng, T., et al. 2021, Optica, 8, 471 [NASA ADS] [CrossRef] [Google Scholar]
  59. Sinclair, L. C., Giorgetta, F. R., Swann, W. C., et al. 2014, Phys. Rev. A, 89, 023805 [NASA ADS] [CrossRef] [Google Scholar]
  60. Sinclair, L. C., Swann, W. C., Bergeron, H., et al. 2016, Appl. Phys. Lett., 109, 151104 [CrossRef] [Google Scholar]
  61. Soffel, M., Klioner, S. A., Petit, G., et al. 2003, AJ, 126, 2687 [Google Scholar]
  62. Stuhl, B. K. 2021, Opt. Express, 29, 13706 [NASA ADS] [CrossRef] [Google Scholar]
  63. Swann, W. C., Bodine, M. I., Khader, I., et al. 2019, Phys. Rev. A, 99, 023855 [NASA ADS] [CrossRef] [Google Scholar]
  64. Synge, J. L. 1960, Relativity: The General Theory (Amsterdam: North-Holland Publishing Company) [Google Scholar]
  65. Taylor, M. T., Belmonte, A., Hollberg, L., & Kahn, J. M. 2020, Phys. Rev. A, 101, 033843 [Google Scholar]
  66. Teyssandier, P., & Le Poncin-Lafitte, C. 2008, Class. Quant. Grav., 25, 145020 [NASA ADS] [CrossRef] [Google Scholar]
  67. Tselikov, A., & Blake, J. 1998, Appl. Opt., 37, 6690 [CrossRef] [Google Scholar]
  68. Turyshev, S. G., Yu, N., & Toth, V. T. 2016, Phys. Rev. D, 93, 045027 [NASA ADS] [CrossRef] [Google Scholar]
  69. Voigt, C., Denker, H., & Timmen, L. 2016, Metrologia, 53, 1365 [NASA ADS] [CrossRef] [Google Scholar]
  70. Wolf, P., & Petit, G. 1995, A&A, 304, 653 [NASA ADS] [Google Scholar]

All Figures

thumbnail Fig. 1

Definition of angles, basis vectors, and other quantities.

In the text
thumbnail Fig. 2

Scheme of the two-way time transfer. The signal leaves observer A at a coordinate time tA1, corresponding to the proper time τA(tA1) of clock A, is reflected from observer B at a coordinate time tB, corresponding to the proper time τB(tB) of clock B and τA(tB) of clock A, and is finally received by observer A at the coordinate time tA2, corresponding to the proper time τA(tA2) of clock A.

In the text
thumbnail 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
thumbnail 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 rA = 6371 km, rBrA = 408 km, and inclined against the rotation of the Earth.

In the text
thumbnail Fig. 5

Three contributions of the spherically symmetric, static part of the atmospheric refractivity , and to the two-way frequency transfer correction (59) and their sum 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 rA = 6371 km and a satellite moving on a circular orbit in the equatorial plane at a radius of rA + 408 km in the direction of the rotation of the Earth.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.