Analytical ray-tracing in planetary atmospheres

TL;DR

This paper develops an analytical framework for modeling photon paths and light times in planetary atmospheres, including non-spherical cases, with high accuracy validated against numerical methods.

Contribution

It introduces a comprehensive mathematical formalism for analytical integration of photon trajectories in atmospheres of any symmetry, extending beyond traditional spherical models.

Findings

Analytical solutions achieve relative errors of 10^-8 for light time and 10^-5 for refractive bending.

The method accurately accounts for atmospheric non-sphericity, such as quadrupolar moments.

Validation against numerical integration confirms the high precision of the analytical approach.

Abstract

Ground-based astro-geodetic observations and atmospheric occultations, are two examples of observational techniques requiring a scrutiny analysis of atmospheric refraction. In both cases, the measured changes in observables are geometrically related to changes in the photon path and the light time of the received electromagnetic signal. In the context of geometrical optics, the change in the physical properties of the signal are related to the refractive profile of the crossed medium. Therefore, having a clear knowledge of how the refractivity governs the photon path and the light time evolution is of prime importance to clearly understand observational features. Analytical studies usually focused on spherically symmetric atmospheres and only few aimed at exploring the effect of the non-spherical symmetry on the observables. In this paper, we analytically perform the integration of the photon path and the light time of rays traveling across a planetary atmosphere. We do not restrict our attention to spherically symmetric atmospheres and introduce a comprehensive mathematical framework which allows to handle any kind of analytical studies in the context of geometrical optics. To highlight the capabilities of this new formalism, we carry out five realistic applications for which we derive analytical solutions. The accuracy of the method of integration is assessed by comparing our results to a numerical integration of the equations of geometrical optics in the presence of a quadrupolar moment $J_2$. This shows that the analytical solution leads to the determination of the light time and the refractive bending with relative errors at the level of one part in $10^8$ and one part in $10^5$, for typical values of the refractivity and the $J_2$ parameter at levels of $10^{-4}$ and $10^{-2}$, respectively.

Tables

Table 1. Table 1: Definition of the trigonometric functions 𝒜 l ​ ( f ) subscript 𝒜 𝑙 𝑓 \mathcal{A}_{l}(f) which are used in this section. We do not specify the dependencies in ι 𝜄 \iota and ω 𝜔 \omega , since they are not expected to vary as quickly as f 𝑓 f at first order in the components of the perturbation.

Equations

$$\mathbf{s}\equiv \frac{\mathrm{d}\mathbf{h}}{\mathrm{d}s}\text{.}$$

$$\nabla \mathcal{S}\equiv n\text{s}$$

$$\nabla \mathcal{S}\cdot \nabla \mathcal{S}=n^2\text{,}$$

$$\frac{\mathrm{d}}{\mathrm{d}s}\big{(}\bm{\nabla}\mathscr{S}\big{)}=\bm{\nabla}n\text{.}$$

$$\frac{\mathrm{d}t}{\mathrm{d}s}=\frac{n}{c}\text{.}$$

$${\tau }_{\text{atm}}=\frac{1}{c}{\int }_{s_1}^{s_2}n\mathrm{d}s-t_{\text{vac}}\text{.}$$

$$\frac{{\nu }_2}{{\nu }_1}=\frac{(u^{\alpha }{k}_{\alpha }{)}_2}{(u^{\alpha }{k}_{\alpha }{)}_1}=\frac{(u^0{k}_0{)}_2}{(u^0{k}_0{)}_1}\frac{(1+{\hat{u}}^i{\hat{k}}_i{)}_2}{(1+{\hat{u}}^i{\hat{k}}_i{)}_1}\text{,}$$

$${\hat{k}}^i={\partial }^i\mathcal{S}\text{.}$$

$$\frac{\nu_{2}}{\nu_{1}}=\frac{\Gamma_{2}}{\Gamma_{1}}\frac{\Big{(}1-\hat{\textbf{u}}\cdot\bm{\nabla}\mathscr{S}\Big{)}_{2}}{\Big{(}1-\hat{\textbf{u}}\cdot\bm{\nabla}\mathscr{S}\Big{)}_{1}}\text{.}$$

$$A_{\text{def}}=10\log \left(1-S\frac{\mathrm{d}ϵ}{\mathrm{d}K}\right)\text{,}$$

$$\alpha \equiv -\frac{\mathrm{d}n}{\mathrm{d}\Phi }$$

$$\mathbf{K}=\mathbf{h}\times \nabla \mathcal{S}$$

$$K=nh\sin \varphi \text{,}$$

$$\mathbf{h}=h\cos \theta \hat{\mathbf{x}}+h\sin \theta \hat{\mathbf{y}}\text{.}$$

$$\hat{\mathbf{\rho }}$$

$$U_0=-\frac{\mu }{h}\text{,}$$

$$n(h)=\eta +\frac{\alpha \mu }{h}\text{.}$$

$$=\frac{p}{1+e\cos \kappa f}\hat{\mathbf{\rho }}\text{,}$$

$$\nabla \mathcal{S}$$

$$\hat{\mathbf{\rho }}=$$

$$\hat{\mathbf{\tau }}=$$

$$\hat{\mathbf{\sigma }}=$$

$$\hat{\mathbf{e}}=$$

$$K\cos \iota$$

$$-K\sin \iota \cos \Omega$$

$$K\sin \iota \sin \Omega$$

$$K\sin \iota \cos \omega$$

$$\sin \iota \sin \omega$$

$$K=\alpha \mu e\text{,}$$

$${\kappa }^2=\frac{e^2-1}{e^2}\text{.}$$

$$p=\frac{\alpha \mu }{\eta }(e^2-1)\text{,}$$

$$=\frac{{\kappa }^{2}K}{\eta }{\left(\frac{\alpha \mu }{K}+\cos \kappa f\right)}^{-1}\hat{\mathbf{\rho }}\text{,}$$

$$\nabla \mathcal{S}$$

$$-a=\frac{\alpha \mu }{\eta }\text{.}$$

Full text

*Aims. *In this paper, we analytically perform the integration of the photon path and the light time of rays traveling across a planetary atmosphere. We do not restrict our attention to spherically symmetric atmospheres and introduce a comprehensive mathematical framework which allows to handle any kind of analytical studies in the context of geometrical optics.

*Methods. *Assuming that the index of refraction of the medium is a linear function of the Newtonian potential, we derive an exact solution to the equations of geometrical optics. The solution’s arbitrary constants of integration are parametrized by the refractive response of the medium to the gravitational potential, and by the first integrals of the problem. Varying the constants, we are able to reformulate the equation of geometrical optics into a new set of osculating equations describing the constants’ evolution for any arbitrary changes in the index of refraction profile.

*Results. *The osculating equations are identical to the equation of geometrical optics, however, they offer a comprehensive framework to handle analytical studies which aim at exploring non-radial dependencies in the refractive profile. To highlight the capabilities of this new formalism, we carry out five realistic applications for which we derive analytical solutions. The accuracy of the method of integration is assessed by comparing our results to a numerical integration of the equations of geometrical optics in the presence of a quadrupolar moment $J_{2}$. This shows that the analytical solution leads to the determination of the light time and the refractive bending with relative errors at the level of one part in $10^{8}$ and one part in $10^{5}$, for typical values of the refractivity and the $J_{2}$ parameter at levels of $10^{-4}$ and $10^{-2}$, respectively.

I Introduction

In the context of Maxwell’s theory of electromagnetism, the refraction phenomenon is understood as the superposition of an incident electromagnetic wave with wavelets generated by some electric particles being accelerated in response to a local electromagnetic field. If the medium containing the particles is sufficiently tenuous (like for a gas), the local field is proportional to the incident field itself. The waves’ superposition generates a resulting signal which exhibits a phase shift with respect to the incident one. This apparent change in the primary signal’s phase is parametrized thanks to a dimensionless parameter called the index of refraction of the medium. It encloses the physical properties related to the interaction between the electromagnetic signal and the material filling the medium. Therefore, any variation of the physical properties of the medium will change the index of refraction which in turn will modify the features of the incident signal. In Astronomy, those changes in the signal’s characteristics are the observables, thus, it is of prime importance to understand the relationship between refraction and light propagation.

Ground-based (GB) astro-geodetic observations and atmospheric occultations (AO) are two examples of observational techniques requiring a careful analysis of refraction. For GB observations, refraction is an accuracy-limiting factor and it has to be modeled in the analysis software; for AO, it is thoroughly modeled to infer the physical properties of the refractive medium. In both cases, the refraction usually occurs at characteristic lengths which are much larger than the wavelength of the electromagnetic signal and the problem can be studied in the approximation of geometrical optics. In that context, the different observables (range, Doppler, or attenuation) can be related to the geometry of the photon path and the light time which are both controlled by the index of refraction. Thus, a clear knowledge of how the geometry of the light propagation is related to the index of refraction is important to understand the changes in observables.

For GB observations, especially for techniques operating for the realization of the International Earth Rotation and Reference Systems Service (IERS), the accuracy of the analysis is largely affected by errors in modeling the group delay during propagation of the signal through the atmosphere Petit et al. (2010). Usually, the atmospheric delay is evaluated at zenith and is computed with an analytical model of Marini and Murray (1973) or with the more recent model (also analytic) of Mendes and Pavlis (2004). Then, the projection to a given elevation angle is operated with the help of the mapping function developed by Mendes et al. (2002). Those models are formulated under the assumption that the atmosphere is spherically symmetric and assume a radial dependency for the index of refraction. Consequently, they leave aside some possible effects due to the presence of horizontal gradients in the atmosphere111Those horizontal gradients can arise at two different levels. Firstly, they can be generated by an atmosphere with a non-spherical shape, and secondly, they can be due to a local horizontal variation of the physical quantities entering in the computation of the index of refraction e.g. the temperature., which might be somehow inaccurate. This issue has only been partially overcome by Chen and Herring (1997) for Very Long Baseline Interferometry (VLBI), and by Hulley and Pavlis (2007) for Satellite Laser Ranging (SLR), by performing numerical integration (numerical ray-tracing) across multi-layered spherical atmosphere, assuming constant refractivity inside each layer. Even if in those cases the index of refraction may possess non-radial dependencies, the shape of the atmosphere is still considered for being spheric. So, it is worth mentioning that beyond the radial dependency of the refractive index, only numerical integrations are usually used for such GB observations.

In AO experiments, the basic idea is to establish radio links between a transmitter and a receiver when the latter is being occulted by a planetary atmosphere as seen from the transmitter. Conversely to GB observations, AO experiments use the effect of the refraction by the atmosphere to determine the physical properties of the medium. This is often referred to as the inverse problem. In the literature, there exists different approaches devoted to the resolution of the inverse problem. For instance, the initial data processing of the Mariner IV occultations was based primarily on model fitting Kliore et al. (1965). The model used for determining the index of refraction profile was approximate and only valid for a planet with a thin spherically symmetric atmosphere Fjeldbo and Eshleman (1965, 1968). Later, it as been shown Phinney and Anderson (1968) that determining directly the refractivity profile from the Doppler frequencies (or from the bending angle) is a special case of Abelian integral inversion when the index of refraction is assumed to be driven by spherical symmetry. In the past, this method has been successfully applied to the dense atmosphere of Venus Fjeldbo et al. (1971) using the data of the occultation of Mariner V. More recently, Abelian integral inversion was used to process radio data of Cassini occultation by Titan Schinder et al. (2011, 2012) or GPS/MET (Global Positioning System Meteorology Experiment) occultation by the Earth atmosphere Kursinski et al. (1997); Steiner et al. (1999). Later, with the Voyager 2 occultation by Neptune in 1989 and the more recent Cassini mission orbiting Saturn, the problem of studying non-spherically symmetric atmosphere was raised. A first analytical model was proposed by Lindal (1992). This solution based on a Taylor series expansion is suited for a numerical evaluation, but does not provide a clear understanding of the effect of non-spherical symmetry on the photon path and the time-delay. More recently, in order to account for additional effects such the light dragging effect by the winds, Schinder et al. (2015) proposed a purely numerical ray-tracing technique to solve the inverse problem. This method is of course more general than analytical models but necessitates a much higher computational cost.

In this context, we propose a fully analytical and comprehensive description of the light path inside planetary atmospheres. The solution is expressed in terms of a free parameter $\alpha$, which characterizes the refractive response of the medium to the gravitational pull of the planet. This parameter is either an input for GB observations or an output for AO experiments. The geometric properties of the ray (as the photon path and the time-delay) are expressed in term of $\alpha$ and in return the observables, which are related to the light geometry, can also be expressed in terms of this parameter. In this paper, we present a useful mathematical framework which allows to perform analytical studies beyond the usual spherically symmetric case. The validity of the derived solution is assessed by comparing it with numerical ray-tracing performed on a simulated atmosphere. We try to keep the discussion as general as possible in order to let the possibility of applying the incoming description to any situation involving light propagation in planetary atmospheres e.g. GB observations or AO experiments.

II Definitions and hypothesis

In this paper, we use the geometrical optics approximation which is characterized by neglecting the wavelength of electromagnetic waves Landau and Lifshitz (1975). In geometrical optics, the concept of rays is introduced as curves whose tangents at each point coincide with the direction of propagation of the electromagnetic waves. In the following, we will mainly focus on the evolution of the separation vector, locating a point along the light ray, and the unit vector tangent to the ray at that point. In this chapter, we introduce the basic principles for determining the evolution of these quantities.

II.1 Equations of geometrical optics

Let us introduce h, the separation vector between the center of a reference frame and a point along the light ray trajectory. Also, we define s, the unit vector which is tangent to the light ray at h’s location

[TABLE]

$\mathrm{d}\mathbf{h}$ is an elemental displacement vector along the ray and $\mathrm{d}s$ represents its length, which is defined as $\mathrm{d}s^{2}=\mathrm{d}\mathbf{h}\cdot\mathrm{d}\mathbf{h}$. Because of this quadratic relation, one can infer that s is a unit vector, and this reduces the number of independent constants of integration from six to five. Together, h and s provide the description of the evolution of the separation vector as well as the direction of the ray.

For convenience, we introduce the new following quantity

[TABLE]

where $\mathscr{S}(\textbf{h})$ is a real scalar function of the position and is usually defined as the optical path Born and Wolf (1999). Here, the gradient is computed with respect to h, and $n(\textbf{h})$ is an arbitrary scalar field called the index of refraction of the medium. Remembering that s is a unit vector, the square of the gradient of the optical path leads to the well-known eikonal equation

[TABLE]

which is the fundamental equation of geometrical optics. The function $\mathscr{S}$ is also referred to as the eikonal.

The height of the ray, its direction of propagation or the eikonal, changes in space according to the spatial variation of the index of refraction of the medium in which the light ray propagates. Thus, it is possible to directly specify the ray with the index of refraction by differentiating Eq. (2) with respect to $s$, which leads to Born and Wolf (1999)

[TABLE]

This expression is given in a reference frame at rest with respect to the refractive medium. Eq. (4) can only be integrated after having introduced some hypotheses about the dependencies of $n$. In the following, we will refer to Eq. (4) as the equation of geometrical optics.

In order to assess a time description along the path, we introduce an additional well-known formula describing the rate of change of the light time with respect to the geometrical length

[TABLE]

The introduction of this additional equation, leads the total number of independent constants of integration from five to six.

Eqs. (1), (4), and (5) are the equations that we have to deal with in order to solve simultaneously the photon path and the light time which in turn can be used to determine the observables, as shown in the next section.

II.2 Observables given by the geometry

We recall how the observables like the range, the range-rate (Doppler shift), or also the defocussing can be determined once $\bm{\nabla}\mathscr{S}$ and $t$ are known. Therefore, this discussion highlights the great benefit of having analytical expressions describing the evolution of the tangent to the ray or the light time.

The range is directly given by integration of Eq. (5). From it, we deduce the group delay expression which is the time-delay due to atmospheric effects

[TABLE]

Here, the integral represents the total light time accounting for the atmospheric effects and $t_{\text{vac}}$ is the light time in vacuum. The integration is performed along the light path, from the transmission point (labeled 1) to the receiving point (labeled 2).

The instantaneous relativistic Doppler shift derives from the fact that the ratio between the received and the transmitted frequencies can be expressed as Synge (1960); Teyssandier et al. (2008)

[TABLE]

where $\nu_{1}$ and $\nu_{2}$ are the transmitted and the received frequencies respectively, $u^{\alpha}$ denotes the 4-velocity vector, $k_{\alpha}$ is the covariant form of the 4-wave vector which is tangent to light ray, and $\hat{u}^{i}$ and $\hat{k}_{i}$ are respectively given by $\hat{u}^{i}=u^{i}/u^{0}$, and $\hat{k}_{i}=k_{i}/k_{0}$ (see Teyssandier et al. (2008); Hees et al. (2014) for definitions and notations). In this expression, the subscripts 1 and 2 specify that the terms between parenthesis must be evaluated at the position of the transmitter and the receiver respectively.

We can specify the frequency ratio in term of $\bm{\nabla}\mathscr{S}$ by following e.g. Leonhardt and Piwnicki (1999). Indeed, for a static medium, it can be shown that

[TABLE]

Using the flat Minkowski metric222Without loss of generality for the discussion, we neglect general relativity effects and focus our attention on the refraction only. If needed, a more generic treatment of the Doppler shift in general relativity can be found in Synge (1960). with the signature $(+,-,-,-)$, we have $k_{i}=-k_{0}\hat{k}^{i}$ and $u^{0}=\Gamma$, where $\Gamma$ is the Lorentz factor of special relativity which is defined as $\Gamma\equiv(1-\hat{\textbf{u}}\cdot\hat{\textbf{u}})^{-1/2}$. In addition, it turns out that $k_{0}$ is constant along null geodesic for that space-time metric, then Eq. (7) becomes

[TABLE]

Finally, one can introduce another observable, namely the attenuation of the signal. The expression for the change in intensity, in dB (also called the attenuation), due to the refractive defocussing, is given by Kliore et al. (2004) as

[TABLE]

where $S$ is the distance from the transmission till the closest approach, $K$ is the impact parameter of the light beam trajectory, and $\epsilon$ is the refractive bending experienced by the ray between the transmission and the reception333If $\epsilon$ can be derived from the cross product between $\bm{\nabla}\mathscr{S}_{2}$ and $\bm{\nabla}\mathscr{S}_{1}$ as $\epsilon=\arcsin(|\bm{\nabla}\mathscr{S}_{2}\times\bm{\nabla}\mathscr{S}_{1}|/n_{1}n_{2})$, in the following, we will use an alternative way for computing it..

Those three quantities [cf. Eqs. (6), (9), and (10)] are three examples of observables that can directly be determined knowing the evolution of $\bm{\nabla}\mathscr{S}$ and $t$ between the transmission and the reception points. However, as seen in Eqs. (1)–(5), they can only be determined after specifying the index of refraction’s dependencies.

II.3 Refractive profile dependency

Except in Sec. IV where we remain as general as possible, in this work we will consider that the index of refraction coincides with the surface of constant $\Phi$, where $\Phi$ is a generalized potential which will be defined later in Sec. V.

This hypothesis is motivated by the hydrostatic equilibrium assumption applied to a body which is made up with fluid elements. Indeed, assuming a non-rotating fluid body, it can be demonstrated (cf. p. 91–92 of Poisson and Will (2014)) that the surfaces of constant density, pressure, and gravitational potential all coincide. This result can even be extended to a steady rotating body as shown in Poisson and Will (2014). However in that case, if surfaces of constant density and pressure still coincide together they do not coincide anymore with surfaces of constant gravitational potential, $U$. Instead, they coincide with the level surfaces of a more generalized potential $\Phi=U-\Phi_{C}$, where $\Phi_{C}$ is known as the centrifugal potential.

So, if we now assume that the index of refraction is related to the spatial distribution of the fluid elements (e.g. in the context of elementary theory of dispersion Born and Wolf (1999)), it becomes obvious that the surface of constant $n$ also coincide with the surfaces of constant $\Phi$. Such a dependence between $n$ and $\Phi$ has also been used recently in Schinder et al. (2015) to analyze the radioscience data of Cassini in the context of occultations by Saturn’s atmosphere.

In order to simplify future notations, we introduce

[TABLE]

which possesses the same dimensions as the inverse of a potential ($L^{-2}T^{2}$).

In the context of GB observations, $\alpha$ must be seen as an input parameter describing the variation of the index of refraction in different layers of the atmosphere. Conversely, in the context of AO experiments it is seen as an output parameter which must be determined following an accurate modeling of the observable. In general, $\alpha$ is not a constant parameter and several measurements are needed through the atmosphere in order to integrate Eq. (11), as it is discussed in Schinder et al. (2015).

II.4 Presentation of the work

In this work, we introduce a mathematical framework which provides a comprehensive description of the light ray trajectory inside planetary atmospheres using geometrical optics.

In Chap. III, we derive a simple solution to Eq. (4) assuming that the index of refraction of a planetary atmosphere is a linear function of the monopole term of the gravitational potential of the central planet. In other words, we assume that $\alpha$ [see Eq. (11)] is constant, which is a reasonable assumption between each layer of a spherically symmetric atmosphere. This solution is then referred to as the reference solution. It is given in terms of constant parameters called hyperbolic elements (in reference to elliptic elements of celestial mechanics) which describe together the shape and the spatial orientation of the photon path as well as the time-delay.

In Chap. IV, we derive a comprehensive mathematical framework which allows to study generic refractive profiles. It is obtained using the reference solution derived in Sec. III as an osculating solution in order to turn Eqs. (1), (4), and (5) into a set of oscultating equations describing the rate of change of the hyperbolic elements following a change in the index of refraction profile. The osculating equations (also called the perturbation equations) are perfectly equivalent to the equation of optics, therefore, in that sense, they can used in place of Eqs. (1), (4), and (5). The main interest is that they are well suited for analytical study of light propagation inside planetary atmospheres.

In Chap. V, we find approximate solutions to the perturbation equations for a large class of additional contributions in the refractive profile of the central planet’s atmosphere. First, we study the effects on the hyperbolic elements due to non-linear dependencies of the refractive profile to the generalized potential. We focus on the limiting case where the generalized potential reduces to the monopole term of the gravitational potential of the central planet (this reduces to the spherically symmetric assumption). Then, we assume that besides the monopole term there exists an additional contribution to the refractive profile, which is simultaneously the centrifugal potential, the axisymmetric part of the gravitational potential of the central planet, and finally an external gravitational potential raised by a perturbing body. We also spotlights the versatility offered by the formalism by studying the light dragging effect caused by a rotating medium. All those perturbations are mainly related to the shape of the atmosphere. The last application concerns the effect due to horizontal gradients444Physically, such gradients may be generated e.g. by horizontal changes in temperature. inside a spherically symmetric atmosphere for a specific application to GB observations.

Finally in Chap. VI, we compare the analytical results obtained for the axisymmetric part of the gravitational potential with a solution obtained from a numerical integration of the equations of geometrical optics (numerical ray-tracing).

III Reference solution

The purpose of this chapter is to build a reference solution, for an idealized problem. As we will see in the next chapter, the main advantage of this reference solution lays in the fact that it is exact and simple, providing a good understanding of the effect of a central refractive profile on the photon path, and also a good starting point to acquire insights in more complex problems.

The reference solution is found by means of successive approximations made on the refractive profile dependencies. The least stringent is the monopole approximation which implies spherical symmetry. The most stringent is the constancy of $\alpha$ which imposes a particular type of spherical symmetry. Each are discussed in turn.

III.1 Spherical symmetry

In this section, we remind some general results which arise assuming a spherically symmetric refractive profile. The spherical symmetry implies that the index of refraction is only function of the magnitude of the separation vector, i.e. $n=n(h)$ with $h=|\mathbf{h}|$. Therefore, it is shown from Eqs. (1) and (4) (see also Born and Wolf (1999)) that the quantity which is defined as

[TABLE]

is constant all along the path of the light ray. Regarding the magnitude, we deduce the following relation which is known as the Bouguer’s rule

[TABLE]

in which $\phi$ is the angle between the separation vector h and the unit tangent vector to the ray s (see Fig. 1). $K=|\mathbf{K}|$ is a constant value called the impact parameter, hence, in the following $\mathbf{K}$ will be referred to as the impact parameter vector.

When the index of refraction is a function of the height only, the impact parameter vector is constant in direction as well as in magnitude and the light ray propagates into a plane which remains orthogonal to K. In such a case, it is helpful to describe the light path with the help of the polar coordinates $h$ and $\theta$ which are defined such that the components of the separation vector are given by

[TABLE]

We have imposed that the $(\hat{\mathbf{x}},\hat{\mathbf{y}})$ plane coincides with the propagation plane. The $\hat{\mathbf{z}}$-direction is collinear to K, the $\hat{\mathbf{x}}$-axis is directed along the line of nodes (defined in Sec. III.3), and the $\hat{\mathbf{y}}$-axis is orthogonal to $\hat{\mathbf{x}}$, and $\hat{\mathbf{z}}$. We take the triad of vector $(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}})$ to form a right-handed vectorial basis that we will refer to as the propagation basis since $(\hat{\mathbf{x}},\hat{\mathbf{y}})$ coincides with the propagation plane of the ray. The geometry is depicted in Fig. 1.

Let us introduce a new rotating vectorial basis which is well suited for the polar description of the light propagation. Later, this frame will be referred to as the polar basis. The first member of the basis is $\hat{\bm{\rho}}\equiv\mathbf{h}/h$. The third member is $\hat{\bm{\sigma}}\equiv\mathbf{K}/K$, which is necessarily orthogonal to $\hat{\bm{\rho}}$ and collinear to $\hat{\mathbf{z}}$. We notice that in the specific case where the index of refraction is a function of the height only, $\hat{\bm{\sigma}}$ is a constant vector. The second member of the basis is $\hat{\bm{\tau}}$ which is defined to be orthogonal to $\hat{\bm{\rho}}$ and $\hat{\bm{\sigma}}$. We take the triad of vectors $(\hat{\bm{\rho}},\hat{\bm{\tau}},\hat{\bm{\sigma}})$ to form a right-handed vectorial basis. In the propagation plane, we have the relations

[TABLE]

III.2 Monopole approximation and constancy of $\alpha$

In Sec. II.3, we have discussed that assuming hydrostatic equilibrium together with elementary theory of dispersion, one can infer that $n=n(\Phi)$ where $\Phi$ is a generalized potential containing the gravitational potential $U$ of the central planet. The first approximation of the gravitational potential is the so-called monopole term (point-mass limit), which is given by the well-known Newtonian potential

[TABLE]

in which $\mu=GM$. Here, $G$ is the Newtonian gravitational constant and $M$ is the mass of the central planet.

Therefore, assuming that the generalized potential contains only the monopole term ($\Phi=U_{0}$), we deduce that the refractivity profile becomes spherically symmetric (see Sec. III.1 for the implications) since it reduces to $n=n(h)$.

We derive the reference solution, for a very peculiar type of spherical symmetry. Indeed, we are to consider that the refractive profile evolves linearly with the monopole term of the gravitational potential. This is the case when $\alpha$ is constant. Therefore, the expression of the index of refraction is explicitly inferred from Eqs. (11) and (16) assuming $\alpha$ constant

[TABLE]

In this expression, $\eta$ is a constant representing the value of the index of refraction at infinity555If the value of the index of refraction is known for a given height, $H$ [let us call $n_{H}=n(H)$ this value], we have the relation $\eta=n_{H}-\alpha\mu/H$.. The quantity which is defined as $N(h)\equiv n(h)-\eta$ is usually referred to as the refractivity of the medium. From Eq. (17), it is seen that the refractivity profile is spherically symmetric and evolves as $1/h$.

We can now determine the quantities related to the description of the geometrical path of photons and the light time of the ray in the propagation plane. The complete derivation is given in appendix A, we just summarize the results here.

From the conservation of $K$ and from Eq. (17), we obtain the evolution of the height of the ray as well as the tangent to the ray, both expressed in the propagation plane [making use of Eqs. (15)]

[TABLE]

The light path, which is described by $h$, is an hyperbola. The parameter $\kappa$ is a constant which is a simple function of $K$ [cf. Eq. (71)] and $p$ and $e$ are two constants of integration being respectively the semi-latus rectum and the eccentricity of the conic-section. The constancy of $p$ and $e$ [cf. Eqs. (73a) and (84)] are both linked to the conservation of the magnitude of the impact parameter. Simultaneously, $e$ is also linked to the conservation of $E$ which is the dimensionless energy [cf. Eq. (82)]. In Eqs. $\eqref{eq:h_ns_sol}$, $f$ is the true anomaly and is defined in Eq. (73b) by means of an other constant of integration, $\omega$, known as the argument of the closest approach. The constancy of $\omega$ is linked to the conservation of the eccentricity vector defined in Eq. (85).

The planar description is completed by using the fact that the index of refraction does not depend explicitly on the length or the time. This gives rise to an additional constant of integration, $S$, being the geometrical length till the closest approach. With that constant we now assess a length description of the trajectory given the true anomaly by means of a Kepler-like equation [cf. Eq. (93)]. In order to directly assess a time description instead of a geometrical length, we derived a second Kepler-like equation for the time [cf. Eq. (94)]. For that last case, the arbitrary constant of integration is, $T$, the elapsed time till the closest approach.

The solution is consequently described in the propagation plane by means of five constants of integration $(p,e,\omega,S,T)$. Its validity has been assessed by comparison with a numerical integration of Eqs. (1), (4), and (5) for the index of refraction given by Eq. (17). The hyperbolic elements remain constants at the level of the numerical noise.

III.3 Hyperbolic path in space

The solution which is described in Sec. III.2 and appendix A achieves a remarkable degree of simplicity. By assuming that the index of refraction is a linear function of the Newtonian potential of the planet, we have shown that the light path is totally contained inside a plane which remains orthogonal to the parameter vector. We have taken advantage of this particularity by choosing for convenience the polar coordinates $(h,\theta)$ attached to the propagation plane. Then, we have demonstrated that the equation of light propagation possesses an exact solution which is an hyperbola contained inside that plane. However, a more general description involving a more generic frame can be preferred. Indeed, when the reference solution is perturbed, the propagation plane is no longer fixed in space, and the introduction of a fixed frame is required to describe the changes. The projection of the reference solution within a more generic frame is the purpose of the present section.

Now, let us introduce the new fundamental reference frame $(\hat{\mathbf{X}},\hat{\mathbf{Y}},\hat{\mathbf{Z}})$. We adopt the $(\hat{\mathbf{X}},\hat{\mathbf{Y}})$ plane as a reference plane and the $\hat{\mathbf{Z}}$-direction as a reference direction. As before, the $(\hat{\mathbf{x}},\hat{\mathbf{y}})$ plane is the propagation plane while the $\hat{\mathbf{z}}$-direction is the orthogonal direction to the propagation plane. The reference plane could be superimposed with the inertial equator of the central planet at a given time and the reference direction would be the direction of the pole. In that case, we impose that the two frames share the same origin. To achieve a complete description of the light path in the reference frame, we have to introduce two additional angles as it can be seen from Fig. 2. The first one, labeled $\Omega$, is called the longitude of the node and is the angle between the $\hat{\mathbf{X}}$-direction and the intersection between the propagation and the reference planes. The second one, labeled $\iota$, is called the inclination and represents the angle between the $\hat{\mathbf{z}}$-direction orthogonal to the propagation plane and the $\hat{\mathbf{Z}}$-direction orthogonal to the reference plane. The line forming the intersection between the two planes is called the line of nodes. From Fig. 2, we can easily determine the components of the unit vectors $\hat{\bm{\rho}}$, $\hat{\bm{\tau}}$, and $\hat{\bm{\sigma}}$ into the reference frame. The radial, tangent, and normal unit vectors are given by

[TABLE]

Regarding the direction of the closest approach, we have

[TABLE]

Thanks to these expressions, the height and the direction of the ray [cf. Eqs. (18)] can now be expressed into the reference frame. We can also provide some simple definitions for the angles of the problem using the solution in the reference frame. For instance, from the previous equations, we deduce the following useful relationships

[TABLE]

We have used the previous expressions of $\hat{\bm{\rho}}$, $\hat{\bm{\tau}}$, $\hat{\bm{\sigma}}$, and $\hat{\mathbf{e}}$. These definitions are fundamental since they are expressed with vectors of the problem. So, an alternative description of the light path consists in using the set of the hyperbolic elements666If we do not specify that s is a unit vector (like in appendix A), we count seven hyperbolic elements. Otherwise, one of the hyperbolic elements is not independent anymore and we end up with six hyperbolic elements. This point is further developed in Sec. III.4. $(p,e,\iota,\Omega,\omega,S,T)$ instead of the h, s, and $t$.

The complete solution described here and in the previous section is the analogous to the solution of the two body problem in celestial mechanics Goldstein (1950); Brouwer and Clemence (1961); Murray and Dermott (2000); Poisson and Will (2014). However, as in celestial mechanics, the usefulness of this solution is not in its direct application, but rather in the fact that it is exact and simple, and can therefore be used as a reference solution (or osculating solution) for studying much more complex problems, as discussed in Sec. IV.

III.4 Limiting case

In this section and the next one, we close the topic related to the hyperbolic solution i) by showing that the solution can be further simplified by making use of the unit vector property of s, and ii) by discussing its direct application to realistic cases.

In the derivation proposed so far, we never made use of the fact that s is a unit vector. This property of s was left aside on purpose, and the reason will become obvious in Sec. IV.2. Here, we are about to further simplify the previous expressions derived in appendix A by using $\mathrm{d}s=(\mathrm{d}\mathbf{h}\cdot\mathrm{d}\mathbf{h})^{1/2}$, which implies that s is a unit vector. Thus, all formulas which are presented in this section are determined under this assumption and should only be applied when the index of refraction is exactly given by Eq. (17). Although they should not be used in general, they provide a simplified and compact version of the formulas derived in appendix A.

If the tangent to the ray is a unit vector, by definition [cf. Eq. (2)] we end up with the eikonal equation [cf. Eq. (3)], and the total number of independent constants of integration becomes six. After inserting the eikonal equation into Eq. (82), we deduce $E=\eta^{2}/2$, and thus, from Eq. (84), one infers

[TABLE]

which shows that the eccentricity reduces to the inverse of the refractivity at $K$. This expression allows to express $\kappa$ [cf. Eq. (71)] as a unique function of the eccentricity

[TABLE]

From these two last equations, it is straightforward to show that the semi-latus rectum [cf. (73a)] is expressed as

[TABLE]

and is therefore a function of the eccentricity only. The meaning is that $p$ and $e$ are not independent constants anymore when s is a unit vector.

Inserting these simplifications into Eqs. (18) provides the evolution of the separation vector as well as the direction of the tangent to the ray

[TABLE]

From Eqs. (22), (24), and (76), it is easily shown that the semi-major axis, $a$, is now given by

[TABLE]

Insertion of Eqs. (22) and (26) into the two Kepler-like equations [cf. Eqs. (93) and (94)], let one to deduce

[TABLE]

as well as

[TABLE]

$F$ is the hyperbolic anomaly [cf. Eqs. (90)]. From these two relationships, we deduce that the two events $(ct,s)$ and $(cT,S)$ are linked by the expression

[TABLE]

This last expression can be useful to convert directly the length in time or vice-versa.

In addition, after having inserted the simplifications into the expression of $\bar{e}$ [cf. Eq. (100)], we obtain $\bar{e}=0$, and we deduce the evolution of the argument of the bending angle

[TABLE]

We can now examine the vacuum limiting case. As already mentioned, the parameter $\alpha$ characterizes the coupling between the geometrical propagation of the light ray and the gravitational attraction of the central planet through the refractive properties of the medium. Thus, it is understood that if $\alpha\rightarrow 0$, the light ray should become insensitive to the presence of the medium and should therefore propagate along a straight line trajectory. This can be easily checked by inserting the expression of $p$ [cf. Eq. (24)] into Eq. (72), and by taking the limit when $\alpha\rightarrow 0$. In this case, it is seen from Eqs. (22) and (26) that

[TABLE]

also, we end up with

[TABLE]

which represents, indeed, the equation of a straight line in polar coordinates (let us note that the vacuum limit is obtained for $\eta=1$). In addition, from Eq. (29), it is seen that $\alpha\rightarrow 0$ implies $c(t-T)-\eta(s-S)=0$, which just reflects the fact that the two events $(ct,s)$ and $(cT,S)$ are linked by a signal traveling along a straight line at the phase velocity $c/\eta$ (which reduces to $c$ when $\eta=1$ in vacuum). The straight line propagation behavior can also be inferred from Eqs. (30) and (31)

[TABLE]

which shows that $\psi$ remains constant whatever is the value of the true anomaly, meaning that the bending angle is null and thus that the trajectory is a straight line.

III.5 Applications

At this level, the formulas presented in Sec. III.4 can already be used either to model atmospheric delays or to process AO data. In this section, we describe succinctly the procedure to simulate an atmospheric time-delay or to analyze Doppler data using the hyperbolic solution described so far.

The solutions in Eqs. (25), (28), and (30) allow us to compute explicit expressions of the observables [cf. Eqs. (6), (9), and (10)] using directly the reference solution.

The expression of the atmospheric time-delay is immediate, since we only have to subtract the light time in a vacuum from the total light time which is explicitly given in Eq. (28). The light time in vacuum is simply given by taking the limit of Eq. (28) when $\alpha\rightarrow 0$ and by using Eq. (90a). We end up with

[TABLE]

where $F_{1}=F(f_{1})$ and $F_{2}=F(f_{2})$.

The expression for the Doppler shift [cf. Eq. (9)] is also immediate since it is given by the scalar product between Eq. (25b) and $\hat{\textbf{u}}$, where $\hat{\textbf{u}}$ is the coordinate velocity vector of the transmitter (subscript 1) and the receiver (subscript 2), divided by the speed of light in vacuum. The substitution leads to

[TABLE]

where $\hat{\bm{\rho}}_{1}=\hat{\bm{\rho}}(f_{1})$ and $\hat{\bm{\rho}}_{2}=\hat{\bm{\rho}}(f_{2})$, likewise for $\hat{\bm{\tau}}$. If Eq. (35) is the most general form that can be obtained for the hyperbolic solution, it is also interesting to express the Doppler frequency shift (see appendix C for the complete derivation) in term of the total refractive bending angle $\epsilon\equiv\psi(f_{2})-\psi(f_{1})$, which is easily inferred from Eq. (30).

Finally, we can also find an expression for the refractive defocussing effect [cf. Eq. (10)] by performing directly the differentiation of Eq. (30) with respect to $K$, remembering that the constants of integration, e.g. the eccentricity, are function of the impact parameter. Such an expression may be tedious to derive, thus a numerical evaluation of the derivative $\mathrm{d}\epsilon/\mathrm{d}K$ from finite differences may be preferred.

From the evolution of the observables, we can now think of an easy procedure to apply to simulate delays and to analyze range-rate measurements. However, we remind once again that all the analytical computations made in Secs. III.4 and III.5 have been carried out assuming that the index of refraction is exactly given by Eq. (17). Furthermore, we have seen in Sec. III.2 that this refractive profile corresponds to a very peculiar type of spherical symmetry ($\alpha=\text{cst}$). So, in order to extend the reference solution to any kind of spherical symmetries ($\alpha\neq\text{cst}$), we are to assume (only for this discussion) that the atmosphere under study is made with $m$ concentric layers, where each layer possesses its own value of $\alpha_{k}$ with $k=1,\ldots,m$. In addition, inside each layer the index of refraction follows Eq. (17) where the value of $\eta_{k}$ must be chosen in such a way that it connects the values of the index of refraction at the interfaces between two layers.

For GB observations, values of $\alpha_{k}$ can be assumed or deduced from different sources. For instance, we can use the profile of temperature and pressure of the International Standard Atmosphere Vaughan (2010), or we can apply the same procedure as the one described in Marini and Murray (1973) or Mendes and Pavlis (2004). For a given geometry, we can thus determine the value of the hyperbolic elements inside each layer. The atmospheric time-delay inside each layer is simply given by Eq. (34), and finally the total delay is the sum of all the contributions.

Conversely, the purpose of AO experiments is to deduce the value of $\alpha_{k}$ at different heights inside the atmosphere. The observables are the Doppler frequency shift measurements, for which the analytical counterpart is given in Eq. (35). It depends on the evolution of $\bm{\nabla}\mathscr{S}$ which relies on $\alpha_{k}$. Therefore, the difference between the observed and the computed Doppler shift can be minimized by varying the value of $\alpha_{k}$. At the end, each layer provides its own determination of $\alpha_{k}$. Therefore, $n$ can be retrieved overall the vertical profile by integrating Eq. (11).

IV Perturbation equations

In the case of a spherically symmetric refractivity, we saw that the direction of K [see Eq. (12)] remains the same while the beam ray propagates through the atmosphere. However, if the gradient of the index of refraction have non-null transverse or normal components, the vector K does not remain constant anymore and the previous light path which used to occur in a fixed plane is no longer a solution.

In this chapter, we address the problem of finding a new solution to a more complex situation where the index of refraction is as much general as possible.

IV.1 Perturbing gradient of the index of refraction

Let us introduce an additional variable contribution besides the hyperbolic term of the index of refraction. We will use the following notation

[TABLE]

with $n_{0}$ being the index of refraction given in Eq. (17) and $\delta n$ being a new non-constant contribution777If $\delta n$ is a pure constant, it must be treated as a new limit value at infinity, such $\eta^{\prime}=\eta+\delta n$.. Therefore, this new contribution will modify the second member of Eq. (4) which is now given by a much more general expression than in Eq. (66)

[TABLE]

Here, $g$ still refers to the magnitude of the local spherical acceleration [cf. Eq. (67)], and the first term on the right-hand side is the gradient of the index of refraction which produces the hyperbolic solution for the light path. The additional term, $\mathbf{f}_{\text{pert}}$, represents a supplementary spatial variation of the index of refraction and is given by $\mathbf{f}_{\text{pert}}=\bm{\nabla}\delta n$. We do not make any assumption about its magnitude with respect to the hyperbolic part, and we do not specify its dependency. In addition, we release the isopotential surfaces hypothesis. We only assume that the perturbing gradient can be decomposed in a certain rotating basis such that its radial, transverse, and normal components are given by

[TABLE]

Here it must be pointed out that the basis $(\hat{\bm{\rho}},\hat{\bm{\tau}},\hat{\bm{\sigma}})$ has a different meaning with respect to the one introduced in Chap. III. Indeed, as discussed above, if $\mathbf{f}_{\text{pert}}=\bm{0}$, only the monopole approximation remains and thus the direction of K remains constant. However, when the perturbing gradient is no longer null, K is changing in direction (also in norm) and then $\hat{\bm{\sigma}}\equiv\mathbf{K}/K$ does not remain fixed in space as before. Good insight can be obtained on the effects of the perturbing gradient on the impact parameter vector by differentiating Eq. (12) with respect to $s$. Making use of Eqs. (4), (1), and (37), we quickly arrive to $\mathrm{d}\mathbf{K}/\mathrm{d}s=\mathbf{h}\times\mathbf{f}_{\text{pert}}$, which can be expressed in terms of the perturbing gradient’s components

[TABLE]

Then, a comparison between Eq. (39) and the derivative of $\hat{\bm{\sigma}}$ with respect to $s$ reveals that the magnitude of K changes according to

[TABLE]

From these two last equations, we predict that a normal perturbation will only affect the direction of the impact parameter vector while a transverse one will also change its magnitude. These equations are perfectly equivalents to the ones of celestial mechanics where only the non-radial components of the perturbing acceleration change the angular momentum vector Burns (1976); Poisson and Will (2014).

We can also demonstrate that K is not the only first integral being affected by the perturbing gradient. Indeed, as seen from the form of the relationships in Eqs. (79) and (85), we can expect that changes in the index of refraction are going to impact the dimensionless energy parameter and the eccentricity vector. Then, since all the hyperbolic elements $(p,e,\iota,\Omega,\omega,S,T)$ are determined from K, $E$, and e (see appendix. A), they are also expected for not being constant anymore. The determination of their evolution is the subject of the next two sections.

IV.2 Method of variation of arbitrary constants

In this section, we will refer to the general solution of the unperturbed problem as

[TABLE]

where $\mathbf{h}_{0}$, $\bm{\nabla}\mathscr{S}_{0}$, and $t_{0}$ are the solutions respectively for the separation vector, the direction, and the light time along the ray for an hyperbolic path, i.e. for $n_{0}$. The $\mathbf{C}$-vector represents the hyperbolic elements which are constants for the unperturbed trajectory. These solutions have been fully described in Chap. III.

We are now interested in studying the very general case where the perturbing gradient is non-null. In this context, the method of variation of arbitrary constants is very well adapted (see e.g. Poisson and Will (2014) for application of the method of variation of arbitrary constants in the context of celestial mechanics). The core of the method is to consider that $\mathbf{h}_{0}$ and $\bm{\nabla}\mathscr{S}_{0}$ [cf. Eqs. (18)] are still solutions of the perturbed problem, avoiding the apparent contradiction by allowing the constants of integration (hyperbolic elements) to change as $s$ evolves. The physical meaning is that at any length $s$ along the ray, the trajectory is taken to be hyperbolic with the parameters $\mathbf{C}(s)$. The general solution must be expressed as

[TABLE]

A subtle consequence of this choice, is that s and $\textbf{s}_{0}$ cannot be two unit vectors at the same time. Indeed, from the definition of $\bm{\nabla}\mathscr{S}$ in Eq. (2) and from (42), we notice that $\textbf{s}_{0}\cdot\textbf{s}_{0}=(n/n_{0})^{2}\ \textbf{s}\cdot\textbf{s}$. Moreover, considering that the eikonal equation (3) must hold true for the index of refraction under study (in this case $n$), we infer that $\mathbf{s}$ is a unit vector, thus $\textbf{s}_{0}$ is not. Consequently, from now, we are not allowed to use the hyperbolic expressions derived in Secs. III.4 and III.5, which are based on the assumption that $\textbf{s}_{0}$ is a unit vector. However, we remind that the solution which is derived in appendix A is free of this assumption, then all the expressions in there can be used.

Another important consequence is that, saving the expression of $\bm{\nabla}\mathscr{S}$ [cf. Eq. (18b)], allows us to fix the expression of the Doppler frequency shift once and for all. Indeed, Eq. (9) depends only of $\bm{\nabla}\mathscr{S}$, and if Eq. (18b) must hold by virtue of the method of variation of arbitrary constant, then the expression of the Doppler frequency shift [as inferred using Eq. (18b)] must also hold true for any arbitrary index of refraction. However, since the hyperbolic elements are not constants anymore, their values at the reception point might be different from the ones at the transmission point.

Applying the proposition (42), we can determine an explicit expression for the variation of the dimensionless hyperbolic energy $E$. Differentiating Eq. (79) and making use of Eqs. (38), and (18), we get the following expression

[TABLE]

in which

[TABLE]

is a new perturbing component. The term $\delta n$ has been introduced in Eq. (36) and refers to the difference between the actual and the hyperbolic index of refraction. We recall that $\delta n$ must contain only the variable part of the difference as discussed in footnote 7. We see from Eq. (43) that only the components of the perturbing gradient which are contained inside the propagation plane can impact the dimensionless parameter $E$.

Considering that all the hyperbolic elements can be expressed in terms of the first integrals $E$, $\mathbf{K}$, and e, we can now infer their length rate of change as functions of the perturbing gradient’s components using Eqs. (39), (40), and (43). However, we first need to express the differentials of the hyperbolic elements as functions of $\mathrm{d}E/\mathrm{d}s$, $\mathrm{d}\mathbf{K}/\mathrm{d}s$, and $\mathrm{d}\mathbf{e}/\mathrm{d}s$.

Let us start with the semi-major axis. Differentiation of Eq. (81) yields

[TABLE]

The equation for the rate of change of the eccentricity is given by differentiating Eq. (84). After some little algebra we find

[TABLE]

The variation of the semi-latus rectum is determined either from Eqs. (45) and (46) making use of Eq. (76), or directly by differentiating Eq. (73a). Both ways lead to

[TABLE]

The variation of the inclination is given after differentiating Eq. (21a)

[TABLE]

Eqs. (21b) and (21c) let to determine the tangent of the longitude of the node. After differentiation, we arrive to

[TABLE]

For the derivation of the rate of change of the argument of the closest approach, we do not differentiate the eccentricity vector as it is usually done in celestial mechanics Poisson and Will (2014). Instead, we are to simplify the algebra by using a geometrical relation Burns (1976). The computation starts by differentiating Eq. (72) and by making use of Eqs. (73), (71), (47), and (46), which leads to

[TABLE]

On the right-hand side we have substituted the following relationship

[TABLE]

which has been determined geometrically from Fig. 3. The left-hand side is obtained by making use of the method of variation of arbitrary constants [cf. Eqs. (42)]. The right-hand side contribution is determined remembering that for a certain location along the ray, h remains unchanged, so we conclude that a tilt along $\Omega$ will decrease $\theta$ by an amount equal to the projection of the magnitude of the tilt on the propagation plane. So, the geometrical contribution to $\theta$ changes is exactly given by the term $-\delta\Omega\cos\iota$. This point is illustrated in Fig. 3 and is also discussed in Burns (1976). With this geometrical relation, the expression of the change in $\omega$ has been determined without using the change in the eccentricity vector.

Finally, we consider the rate of change of the true anomaly for closing the system, even if it is not constant along the photon path. The relation is obtained by differentiating Eq. (73b) and by substituting it into Eq. (51)

[TABLE]

The left-hand side represents the hyperbolic contribution to the change in the true anomaly. The right-hand side represents the change in the direction of the closest approach relative to the $\hat{\textbf{X}}$-direction and measured inside the propagation plane (see Fig. 3). The same is also true in celestial mechanics as mentioned by Poisson and Will (2014). Substituting Eq. (50) into Eq. (52) provides the following expression

[TABLE]

The right-hand side reduces to $\hat{\bm{\tau}}\cdot\bm{\nabla}\mathscr{S}/nh$ when the perturbation is null (i.e. when $\mathrm{d}E/\mathrm{d}s$ and $\mathrm{d}K/\mathrm{d}s$ vanish).

IV.3 Osculating equations

We are now ready to derive the final expressions for the evolution of the hyperbolic elements as functions of the components of the perturbing gradient of the index of refraction. However, for a matter of compactness, we do not provide the differentials of $(S,T)$, but rather the differentials of $(f,t)$, even if $f$ or $t$ are not constants along the light path for the hyperbolic trajectory.

Substituting Eqs. (39), (40), and (43) into all the previous equations describing the rate of change of the hyperbolic elements with $\mathrm{d}E/\mathrm{d}s$, $\mathrm{d}\textbf{K}/\mathrm{d}s$, and $\mathrm{d}K/\mathrm{d}s$, it is possible to infer the following relationships

[TABLE]

These equations provide a good understanding of how the components of the perturbing gradient govern the length rate of change of the hyperbolic elements. For instance, it can be seen that $p$ is only affected by the transverse component of the perturbing gradient, where $\iota$ and $\Omega$ only change because of the normal component. We also see that $e$ and $f$ are affected by the radial and the transverse components, therefore any perturbation which is contained inside the propagation plane will induce variations of these two elements. In addition, they are also expected to vary because of a term in $\delta n$ (through $\mathcal{N}$), which represents a non-constant change in the expression of the index of refraction. Finally, we see that $\omega$ is the only one being affected by all the components of the perturbing gradient.

This set of equations is similar to the osculating equations of celestial mechanics which are also called perturbation equations. However, one important difference stands in the terms in $\mathcal{N}$ which does not possess any equivalent in celestial mechanics. The reason is that the motion of planets is not directly sensitive to the gravitational potential itself but rather to its gradient. In the case of geometrical optics, the equations of light propagation not only contain the gradient of the index of refraction, but also the value of that index, see e.g. Eq. (5). This can also be inferred once Eq. (4) is written as $\mathrm{d}\textbf{s}/\mathrm{d}s=n^{-1}[\bm{\nabla}n-\textbf{s}(\textbf{s}\cdot\bm{\nabla})n]$, which represents the rate of change of the curvature of the ray Born and Wolf (1999).

As in celestial mechanics, when the inclination is null, $\omega$ and $\Omega$ are not defined since the direction of the line of nodes is not determined as seen in Fig 2. These singularities are related to the choice we have made in the definition of angles. They can be removed Brouwer and Clemence (1961); Murray and Dermott (2000); Morbidelli (2002); Poisson and Will (2014) by introducing new angles, such the longitude of the periapsis which is defined as $\varpi=\omega+\Omega$. In the next, we will continue with $\omega$ and $\Omega$, keeping in mind that a null inclination induces coordinate singularities.

The perturbation equations provide an alternative description to the fundamental equations of optics [cf. Eq. (4)], also in that sense, they are generic and no approximations have been introduced in the transcription. Their validity has been assessed by comparing numerical integration of Eqs. (54) with numerical integration of Eqs. (1), (4), and (5) showing that the differences remain at the level of the numerical noise. The main advantage of the perturbation equations is to provide an easy geometrical picture to tackle the effects of a perturbing gradient with respect to an hyperbolic path.

It is worth noticing that any exact solution to the fundamental equations of optics could have been used as a reference solution, or osculating solution. For instance, we could have chosen the straight line as the osculating solution instead of the hyperbolic path. Actually, if we insert Eqs. (22)–(24) into Eqs. (54), and if we take the limit $\alpha\rightarrow 0$ (which implies $e\rightarrow\infty$), the perturbation equations reduce to a set of equations describing the evolution of an osculating straight line.

In the following, we will see how to use the set of perturbation equations when the perturbing gradient can be considered small with respect to the spherical contribution of the gravitational potential of the central planet.

IV.4 First-order approximation

The great benefit of having written the equation of fundamental optics in terms of the components of the perturbing gradient will become obvious in this section.

Indeed, in the case where the perturbing gradient is small with respect to the hyperbolic contribution (i.e. $\alpha g\gg|\mathbf{f}_{\text{pert}}|$), the changes in the hyperbolic elements are expected to be small, too. Thus, a good understanding can be reached by inserting the constant values of the hyperbolic elements in the right-hand side of Eqs. (54) in order to get the first-order approximation. For this, it is convenient to use the true anomaly as independent variable instead of the geometrical path length. This requires to invert Eq. (54f) keeping only the linear order in the components of the perturbing gradient. We find the following relations

[TABLE]

We can also add an expression for the evolution of the argument of the refractive bending which is needed to determine the defocussing and which can also be used to approximate the Doppler frequency shift (cf. discussion in Sec. III.5). To do so, one can differentiate Eqs. (13) and (95), and then makes use of Eqs. (96) to deduce

[TABLE]

which is more general than Eq. (97) since it includes, besides the hyperbolic term, the variations of $K$ and $\omega$. Substituting the differentials into that equation and keeping only the first order in the components of the perturbation, we end up with

[TABLE]

where we have removed the hyperbolic contribution, $\psi_{0}$, which is given exactly in Eq. (99). We have defined $\delta\psi=\psi-\psi_{0}$ as the non-hyperbolic contribution to the argument of the refractive bending. Once integrated, that last equation will give us the net change in the refractive bending due to a perturbation and then, Eq. (148) gives us directly the first order effect on the Doppler frequency shift.

IV.5 Method of integration

Eqs. (55) are general enough to be applied to different geometries and different kind of problems as AO experiments or GB observations. In both cases, the integration can be performed along the hyperbolic light path to get the first order effects on the hyperbolic elements. Calling $\mathbf{C}$ the vector of the following elements, $\mathbf{C}=(p,e,\iota,\Omega,\omega,\delta s,\delta t,\delta\psi)$, we have

[TABLE]

where $f_{1}$ and $f_{2}$ are respectively the true anomaly at the transmission and at the reception point, with $f_{2}\geq f_{1}$. The term $\mathrm{d}\mathbf{C}/\mathrm{d}f^{\prime}$ represents the rate of change of the elements with respect to the true anomaly and is obviously given by Eqs. (55).

In general, direct integration, as in Eq. (56), may lead to complex solutions. So, if small quantities show up in Eqs. (55), a Taylor expansion in term of these quantities can first of all let to simplify the integration and secondly lead to a more user-friendly solution.

For instance, in the context of AO experiments and GB observations, we can think of the inverse of the eccentricity as a small quantity. Indeed, in most experiments, we have to deal with a small refractivity which implies $\alpha\ll 1$ and then, from Eq. (22), we deduce $e\gg 1$ (as long as $K\gg\alpha\mu$). The meaning is that the shape of the hyperbola departs slightly from a straight-line trajectory. However, if this approximation is largely relevant for AO, it can become somehow inaccurate for GB experiments, in particular at very high elevations (close to the zenith-direction of the site). Indeed, $K$ might become the same order as $\alpha\mu$, which implies $e\sim 1$. Hence, the solutions at first order in $1/e$ are not valid for GB experiments at very high elevations (i.e. when $K\sim\alpha\mu$). For instance in the Earth’s atmosphere 888Following Aparicio and Laroche (2011), for an average parcel of air at sea level, the refractivity is approximately $(3\pm 1)\times 10^{-4}$, which represents a value of $\alpha=(5\pm 2)\times 10^{-6}\,\mathrm{km^{-2}\cdot s^{2}}$. So, the region around the zenith-direction, in which the eccentricity is $e\lesssim 10$, is enclosed inside a cone with a half top angle which remains $\lesssim 0.3^{\circ}$. All the light path trajectories with $K\lesssim 30\,\mathrm{km}$, are enclosed inside that region., if we consider that the expansion in $1/e$ does not hold for $e\lesssim 10$, we must consider only trajectories with $K\gtrsim 30\,\mathrm{km}$, which corresponds to elevations angles $\lesssim 89.7^{\circ}$.

For GB observations at very high elevations, the light ray trajectory is close to be radial, that is to say that the change in the true anomaly between the transmitter and the receiver $(f_{2}-f_{1})$ is a small quantity. For that particular case, we do not need to integrate Eqs. (56), instead an evaluation of Eqs. (55) at the level of the transmitter is sufficient to determine the change in the hyperbolic elements at first order in $(f_{2}-f_{1})$

[TABLE]

Therefore, for all the following applications (see Chap. V), we consider the leading order in $1/e$ since: i) it encompasses almost all experimental cases (usually, at zenith angles less than few degrees, source tracking can be difficult), ii) the transposition to GB observations at very high elevations (zenith angles less than few degrees) reduces to an evaluation of Eqs. (55) at the level of the transmitter.

V Applications

In this chapter, different types of perturbations are analyzed within the osculating equations formalism. In each application, we aim at determining the variations of the photon path and the light time due to these perturbations. Then, the change in the observables can be easily inferred as discussed in Sec. II.2.

We first consider that the index of refraction possesses non-linear dependences to the monopole term of the gravitational potential. Then, we explore the impact of the centrifugal potential as well as the drag effect on the light propagation due to the moving medium. We also apply the perturbation equations considering the non-spherical part of the self-gravitational potential of the central planet. Then, we study the perturbation on the light beam trajectory caused by a perturbing body raising tides on the central planet’s atmosphere. Finally, we close the section by studying the effects due to the presence of horizontal gradients in a spherically shaped atmosphere.

V.1 Definitions and hypothesis

In this section, we introduce the basic ideas which will help us to formulate the internal problem999In the Solar System, the internal problem (study of the internal structure of a self-gravitating bodies) can be mainly decoupled from the external one (inter-body dynamics). This is true as long as the inter-body distance is much more higher than a characteristic length scale within the extended bodies. Consequently, an external stress can be considered as a perturbation in the internal problem and conversely, an internal stress may be seen as a perturbation in the external problem.. The idea is to determine a generic refractive profile assuming the influence of different sources of stress.

Lately, we have made use of three reference frames, all centered at the planet’s center of mass and referred to as the propagation frame $(\hat{\mathbf{x}},\hat{\mathbf{y}},\hat{\mathbf{z}})$, the polar basis $(\hat{\bm{\rho}},\hat{\bm{\tau}},\hat{\bm{\sigma}})$, and the reference frame $(\hat{\mathbf{X}},\hat{\mathbf{Y}},\hat{\mathbf{Z}})$. The first one was helpful to derive the solution of reference (hyperbolic path) by making use of the second one to introduce polar coordinates. The last one was introduced to infer the perturbation equations in an arbitrary oriented frame. It was assumed that the medium was at rest simultaneously in the propagation frame and in the frame of reference, also the equations of geometrical optics [cf. Eq. (4)] were consequently available inside these two frames.

When considering the internal motion of the medium with respect to the reference frame, the equation of geometrical optics [cf. Eq. (4)] no longer stands into the propagation or the reference frames since it is expressed in a frame comoving with the medium. Therefore, in order to maintain the coherence with the reference solution which was expressed in the reference frame, we deal with the theory of light propagation in moving medium Fresnel (1818). Following Leonhardt and Piwnicki (1999) and Rozanov and Sochilin (2005), the equation of geometrical optics expressed in the medium comoving frame is

[TABLE]

where v is the velocity vector of the medium, with $v=|\textbf{v}|$. $\gamma$ is the Fresnel’s dragging coefficient introduced in Eq. (83). It might be seen from this expression, that the corrections to Eq. (4) are of the order of $v/c$. For gas giant planets in the Solar System (e.g. Jupiter), the typical upper bound value of the velocity for the solid rotation at the equator is $5\times 10^{4}\,\mathrm{km/h}$ with zonal winds asymmetric velocity ranging between $\pm 550\,\mathrm{km/h}$ Kaspi et al. (2018), also $v/c\lesssim 10^{-4}$. Therefore, we will consider the effect of the moving medium as a perturbation to the hyperbolic path,

[TABLE]

This contribution will be considered later in the context of a steady rotating atmosphere.

To handle the internal motion let us introduce a new reference frame $(\hat{\mathbf{X}}^{\prime},\hat{\mathbf{Y}}^{\prime},\hat{\mathbf{Z}}^{\prime})$ centered at the planet’s center of mass and rotating with it (we call $P$ the central planet). The $\hat{\mathbf{Z}}^{\prime}$-axis is chosen to be orthogonal to the equatorial plane of $P$. The angular velocity vector is considered to be constant and given by $\textbf{w}=w\hat{\mathbf{Z}}^{\prime}$. Since w is independent of time, the $\hat{\mathbf{Z}}^{\prime}$-axis is spatially fixed. Therefore, we choose the $\hat{\mathbf{Z}}$-axis of the reference frame to be collinear with the $\hat{\mathbf{Z}}^{\prime}$-axis. The frame $(\hat{\mathbf{X}}^{\prime},\hat{\mathbf{Y}}^{\prime},\hat{\mathbf{Z}}^{\prime})$ is well suited for the study of internal motions of $P$ and will be referred to as the fluid rotating frame.

It is shown in Poisson and Will (2014) (cf. p. 106) that Euler’s equation describing the time evolution of the velocity of a fluid element makes appear the following generalized potential101010We took the convention that the Newtonian gravitational potential is a negative quantity, also the signs in that expression differs from Poisson and Will (2014).

[TABLE]

once written in the fluid rotating frame. In the following, we consider that the dominant term of that expression is the monopole term $U_{0}$ of the gravitational potential of the fluid planet. The other terms will be considered as perturbations before $U_{0}$. The $U_{l}$ terms are the non-spherical parts of the central planet self gravitational potential, $\Phi_{C}$ is the centrifugal potential due to planet’s proper rotation, and $U_{\text{tidal}}(t)$ is the external potential due to the presence of other massive bodies in the surrounding of the central planet. This last term is dynamic since it evolves with the positions of the perturbing bodies. Each contribution will be studied in turn.

We now have to define an expression for the index of refraction. If we assume that the refractive index is directly related to the density of the fluid, we can consider that Eq. (60) constitutes the generalized potential in Eq. (11). In addition, successive integrations by part of Eq. (11) leads to the following infinite series

[TABLE]

where $\Phi$ is given by Eq. (60), and $\alpha_{(k)}\equiv\mathrm{d}^{k}n/\mathrm{d}\Phi^{k}$ for $k\geq 2$. We keep the same notation as in previous chapters, i.e. $\alpha\equiv\alpha_{(1)}$ is still the linear term with the generalized potential. This series can represent any function of the gravitational potential, and the choice is governed by the numerical values which are assigned to the $\alpha_{k}$ coefficients.

Eqs. (58)–(61) are all we need to formulate the internal problem and study first order effect on the light propagation due to the different pieces in Eqs. (60) and (61).

V.2 Non-linearity

Previously in Chap. III, we have assumed $\mathrm{d}n/\mathrm{d}\Phi=\text{cst}$ [cf. Eq. (11)] in order to simplify the integration of Eq. (4). However, this could be too restrictive in some applications where it is not possible to build a multi-layers modeling of the atmosphere. If a multi-layers modeling can be achieved as in Schinder et al. (2015) or in Hulley and Pavlis (2007), a constant value of $\alpha$ can be associated to each layer and the total profile can be recovered by integrating Eq. (11).

In this section, we propose to explore the first order effect of non-linearity with isopotentials in the refractive profile [see Eq. (61)]. We focus on the simple case where the generalized potential is given by the monopole term of the Newtonian potential, $\Phi=U_{0}$. All the mathematical details are exposed in Sec. B.1 and the change in hyperbolic elements for any degree $k$ is given at leading order in $1/e$ in Eq. (103). The application to the quadratic contribution $(k=2)$ is given in Eq. (104).

After integration, it is shown that whatever function of the Newtonian potential the index of refraction is, $p$, $\iota$, and $\Omega$ remain always constants regardless degree $k$. This is due to the fact that the index of refraction acts like a central field with no transverse or normal components.

We copy here the expressions of the non-hyperbolic contribution to the light time and the refractive bending for $f_{2}\geq f_{1}$ and $k=2$

[TABLE]

Obviously, from an observational point of view, these expressions are of a great interest. The first one provide directly the additional time-delay with respect to the hyperbolic path, and the second one provides the supplementary contribution to the hyperbolic refractive bending. Both are function of the geometry of the problem only. In addition, from Eq. (148), it can be seen that the later provides directly the change in the Doppler measurement due to $\alpha_{(2)}$.

V.3 Steady rotating atmosphere

The second application concerns the rotation of the atmosphere. All bodies in the Solar System are rotating. The rotation motion can be decomposed into two parts. The first concerns the proper rotation of the body around a certain axis of rotation, while the second concerns the spatial orientation of that axis. The generalized potential in Eq. (60) neglects the second part and considers a steady rotating planet.

In the frame rotating with the fluid, the media experiences the effect of the proper rotation via the simplified following generalized potential $\Phi=U_{0}-\Phi_{C}$, where $\Phi_{C}$ represents the centrifugal potential. In addition, the dragging effect due to the rotational motion of the medium is described by the additional contribution $\textbf{f}_{\text{drag}}$ in Eqs. (58). The solutions describing the change in hyperbolic elements due to a steady rotating atmosphere are given in Eqs. (108) and (111). All the mathematical details are presented in Sec. B.2.

On one hand, Eqs. (108) describe the change in hyperbolic elements due to a centrifugal potential induced by the proper rotation of the planet. The expressions of the non-hyperbolic contribution to the light time and the refractive bending are

[TABLE]

with $f_{2}\geq f_{1}$.

On the other hand, Eqs. (111) describe the change in hyperbolic elements induced by the light dragging effect due to the rotation of the fluid material which composes the atmosphere of the central planet. The expressions of the non-hyperbolic contribution to the light time and the refractive bending are

[TABLE]

with $f_{2}\geq f_{1}$.

These relationships show that the first set is quadratic with the magnitude of the angular velocity vector while the latter evolves linearly with it. Therefore, the change in hyperbolic elements due to the centrifugal potential is independent of the orientation of w and remains the same for either a direct or an indirect rotation. On the contrary, the change in hyperbolic elements due to the dragging of light is dependent of the direction of the medium’s rotation. Because of a cosine of the inclination, we can also see that the effect is maximum in the planetary equator, where the distance with respect to the axis of rotation is maximum.

Another important difference between the effects due to the centrifugal potential and those due to the dragging of light can be inferred by comparing the expressions of the non-hyperbolic contribution to the light time and the path length. In the case of light dragging, it is seen in Eq. (111g) that the non-hyperbolic contribution to the light time reduces to the non-hyperbolic contribution to the geometrical length divided by the speed of light in vacuum. This comes from the fact that the $\mathcal{N}$ perturbing component is null for the light dragging effect while it is not for the centrifugal potential contribution.

V.4 Axisymmetric gravitational potential

Because of their proper rotation, non-rigid bodies tend to be flattened due to centrifugal forces. Then, the mass repartition becomes slightly different from what would be expected in spherical symmetry, and consequently, the gravitational potential also changes. This mass redistribution due to centrifugal forces is usually the most important departure from the monopole contribution to the total self-gravitational potential of the planet.

The non-spheric components of the gravitational potential of the central planet can be modeled using the simplified following generalized potential $\Phi=U_{0}+\sum_{l}U_{l}$ with $l\geq 2$. The solutions for the change in the hyperbolic elements are derived in Sec. B.3 and the results for $l=2$ are given in Eqs. (121). We copy here the expressions of the non-hyperbolic contribution to the light time and the refractive bending for $f_{2}\geq f_{1}$

[TABLE]

Considering that the centrifugal potential and the quadrupole moment of the gravitational potential of the central planet are both axisymmetric fields, one might expect similar signatures in the changes induced on the hyperbolic elements. However, because they evolve differently with $h$ (the centrifugal effect tends to grow with $h$, where the $J_{2}$ effect decreases with $h$), a comparison shows that the signatures produced on the hyperbolic elements differ. However, one can see that the ratio between Eqs. (121) and (108) is always proportional to

[TABLE]

where C represents the set of hyperbolic elements. The subscript $J_{2}$ refers to the effect due to the quadrupolar moment of the gravitational potential, where the subscript $C$ refers to the contribution due to the centrifugal potential.

For a planet like Jupiter, considering that the impact parameter of the ray remains at the level of the equatorial radius, the ratio is $0.16$. For Saturn, the Earth, and Titan we find respectively $0.11$, $0.31$, and $0.75$. In each case the centrifugal contribution is the most important. However, it must be noticed that this ratio can grow while the light ray goes deeper inside the atmosphere, since the impact parameter decreases. This can make the $J_{2}$ effect being more important than the centrifugal one, especially for Titan for which the rotation rate is slow. The ratio can also exceed unity for GB observations carried out at high elevation because $K$ might become much smaller than the equatorial radius.

V.5 Tidal potential

In systems containing close orbiting bodies, e.g. a planet and its satellites, the question of knowing whether the tidal effect due to a perturbing body on the planetary atmosphere can affect the light propagation or not can be raised.

To study this possibility we focus on the tidal contribution in the generalized potential and ignore the body’s rotational deformation and the non-spherical gravitational potential contribution which were treated previously. Hence, we consider the simplified following generalized potential $\Phi=U_{0}+U_{\text{tidal}}(t)$. In addition, we will assume that the characteristic time of variation of the external tidal field is so slow that it never takes the central planet’s atmosphere out of hydrostatic equilibrium. This assumption is known as the static tides approximation.

All the computations are detailed in Sec. B.4, and the evolution of the hyperbolic elements are given in Eqs. (135). We remind that these equations are derived considering a tide raising body moving on an equatorial circular orbit around the central planet, and a photon path lying inside the equatorial plane (i.e. $\iota=0$). We copy here the expressions of the non-hyperbolic contribution to the light time and the refractive bending for $f_{2}\geq f_{1}$

[TABLE]

Comparison between the magnitude of the change in hyperbolic elements in Eqs. (135) and (121) reveals

[TABLE]

For the tides raised by the Moon on Earth, the ratio is of the order of $10^{-6}$. For tides raised by Titan on Saturn, the value is $10^{-7}$, and for tides raised by Io on Jupiter the ratio is of the order of $10^{-5}$. For all those cases, the tides effects are at least 5 orders of magnitude smaller than the one of the oblateness. However, they can become important when the light ray passes through the atmosphere of a satellite orbiting a central massive planet. In such a configuration the tidal effects are expected to be more important. For instance, let us consider the case of a light ray crossing Titan’s atmosphere; substitution of numerical values (we took the Titan’s $J_{2}$ value from Iess et al. (2012)) reveals that the changes due to Saturn’s tides are expected to be the same order of magnitude than changes due to Titan’s oblateness.

All the results in Eqs. (135) are derived under the assumption that the main planet is made with a perfect fluid medium having a non-unity index of refraction. However, such a perfect fluid model might be somewhat inaccurate in the context of tidal dynamics. The main reason is linked to the fact that perfect fluid model does not admit a mechanism to dissipate energy and the fluid’s response to a tidal stress is purely elastic. In the context of celestial mechanics, energy dissipation mechanisms are of prime importance in a large number of applications Murray and Dermott (2000). In Sec. B.4, we compute the first order effect on the hyperbolic elements caused by the non-elastic response of the atmosphere following a tidal stress. Referring to Eq. (138), we see that the non-elastic contribution is a least $w\tau$ times the elastic one, with $w$ the magnitude of the angular velocity of the central planet, and $\tau$ the time-delay which is related to dissipative phenomenon.

For the Earth, a typical value of $\tau$ corresponding to rotational semi-diurnal deformation is $4\,\mathrm{min}$ Folkner et al. (2014), so the ratio is of the order $2\times 10^{-2}$. So, if the elastic contribution is important in some circumstances, the viscoelastic response of the atmosphere represents a perturbation of a few percent of the amplitude of the elastic one. In the case of Saturn and Jupiter, we can use the fact that $w\tau\approx(2Q)^{-1}$ where $Q$ is the specific dissipation function Murray and Dermott (2000). The numerical substitution is done by taking value of $k_{2}/Q$ (with $k_{2}$ the gravitational Love number of degree 2) from Lainey et al. (2009) together with the value of $k_{2}$ provided by Iess et al. (2018) for Jupiter and the one provided by Lainey et al. (2017) for Saturn. We find the ratios to be respectively of the order of $8\times 10^{-6}$ and $7\times 10^{-5}$, which is well beyond the elastic response.

V.6 Horizontal gradients

The last study is a direct application for GB observations operating for the realization of IERS reference systems e.g. SLR, VLBI, or the Global Positioning System (GPS). In Chen and Herring (1997) or in Hulley and Pavlis (2007) it is shown that the group delay due to the horizontal gradients in the Earth atmosphere can overpass the centimetric level while observational techniques like SLR or VLBI currently operate with sub-centimeter accuracy measurements. Therefore, not considering horizontal gradients may lead to systematic errors in the estimation of the station coordinates, which can have repercussions on the determination of the International Terrestrial Reference Frame (ITRF). Consequently, delays due to horizontal gradients must be taken into account, especially at low elevations where the effect is maximum.

The computations are carried out in Sec. B.5, and the equations describing the change in the hyperbolic elements following horizontal gradients are given in Eqs. (146). We remind that these equations are derived for a simplified geometry. Indeed, we have assumed that the transmitting source is observed at its highest elevation as seen from the receiving site on Earth. This situation is encountered when the azimuth of the transmitter equals the azimuth of the observer, that is to say when the source is at the meridian of the receiving site. In such a case, the propagation plane is aligned with the meridian and intersects Earth’s equator for $\iota=\pi/2$ and $\lambda=\Omega$. For this type of inclination, it is seen from Eqs. (146) that only a West-East horizontal variation of the refractivity can change the inclination or the longitude of the node. Conversely, only a North-South horizontal gradient is expected to change the other hyperbolic elements.

We copy here the results for the non-hyperbolic contribution to the light time

[TABLE]

For applications, the true anomaly is not a common way of expressing the results. Instead, the colatitude (referred from the direction of the North pole) may be preferred. In order to pass from true anomalies to colatitudes, we make use of hyperbolic relations defined in appendix A, in particular Eqs. (73b) and (74). Because, this transformation is only a matter of rewriting the boundary conditions of Eq. (62), we can safely approximate the photon path assuming straight line between the source and the receiver. In this case, the hyperbolic relations in appendix A can be further simplified taking the limit $\alpha\rightarrow 0$, which corresponds to light propagation in vacuum.

As an example, we give the transformation rule which allows us to pass from true anomaly to colatitude. For the photon path $1\rightarrow 2^{\prime}$, which is depicted on the left panel of Fig. 4, we get

[TABLE]

Finally, insertion of Eqs. (63) into (62) allows one to infer the atmospheric time-delay along the photon path $[\delta t(\bar{\varphi})]$ and for the following boundary conditions $(q,\bar{\varphi}_{1},\bar{\varphi}_{2})$.

In order to emphasize the behavior of $\delta t(\bar{\varphi};q,\bar{\varphi}_{1},\bar{\varphi}_{2})$, we have represented in Fig. 4, the evolution of the North-South horizontal gradient contribution to the atmospheric time-delay. The computation has been carried out assuming three stations (labeled $2$, $2^{\prime}$, and $2^{\prime\prime}$) which are simultaneously ground tracking a unique source (labeled $1$). We have assumed that the stations are located at different colatitudes along the same meridian.

For a given ratio $q$, it is seen that the largest effect is reached for station $2$, which is ground tracking the source at the minimal elevation. This is due to three cumulative effects. Firstly, the projection of the horizontal gradient along the direction of the ray is maximum for the path $1\rightarrow 2$. Secondly, the geometrical length inside the atmosphere is the most important for the path $1\rightarrow 2$. Finally, the horizontal contribution in $(\varphi-\varphi_{2})n_{\text{NS}}$ into the expression of the index of refraction [see Eq. (139)] is the most important for $\varphi=\varphi_{1}$ and for the path $1\rightarrow 2$ which presents the most important value of $\varphi_{1}-\varphi_{2}$. Consequently, for this path, the phase velocity is minimum and consequently the time-delay is maximum.

VI Comparison with numerical integration

The validity of the approximated solutions which are derived in previous section can be assessed by a direct comparison with the output of a numerical integration of the light path.

As discussed so far, the solutions can indifferently be applied to AO experiments or astro-geodetic GB observations. The differences remain in the role of $\alpha$, which is either an output or an input, and also in the geometry through the values of the true anomaly at the transmission and at the reception (or at the entrance and the exit of the region of refractive influence). However, we have seen that the method of computation for AO experiments may differ for GB observations, especially at very high elevations for which an expansion in $1/e$ might not be accurate enough.

In this chapter, we test the validity of the derived solutions in the context of AO experiments. We focus on Eqs. (121) which solve, in a non-ubiquitous way, the problem of determining analytically the light path in the presence of the atmosphere’s oblateness. As discussed in Chap. I, this problem has never been solved analytically in a complete and satisfactory way.

In order to assess the validity of Eqs. (121), we have performed a numerical integration of Eqs. (1), (4), and (5) across the refractive profile given in Eqs. (115) and (117) for $l=2$. The numerical integration has been carried out in double precision, with a numerical error tolerance of $10^{-11}$, for different values of $\alpha$ and $J_{2}$. [Actually, instead of $\alpha$, we work in the following with the dimensionless parameter $\alpha\mu/R$, which represents the hyperbolic refractivity evaluated at the level of the radius of the central planet, $N_{0}(R)=n_{0}(R)-\eta$]. The tested values of the refractivity and $J_{2}$, range from $10^{-1}$ to $10^{-5}$, and from $10^{-1}$ to $10^{-8}$ respectively. For each numerical integration, we have compared the total change in the hyperbolic elements, between the transmission and the reception, with the analytical predictions given in Eqs. (121). Let us emphasize that $\delta s$, $\delta t$, and $\delta\psi$ cannot be determined easily from the numerical integration of Eqs. (1), (4), and (5), thus, instead of working with the non-hyperbolic contributions alone, we are considering the total change in $s$, $t$, and $\psi$, which is easily inferred from results of the numerical integration. From an analytical point of view, the total change in $s$, $t$, and $\psi$, is simply given by the sum of hyperbolic [see Eqs. (93), (94), and (99)] and non-hyperbolic contributions [see Eqs. (121f)–(121h)]. In Fig. 5, we show the evolution of the relative error on the total change in the hyperbolic elements and $s$, $t$, and $\psi$.

Firstly, it is seen that the solutions, for low refractivity [$N_{0}(R)$ below $10^{-4}$] and for small values of $J_{2}$, seem to suffer from a lack of accuracy (c.f. errors on $p$, $e$, $\iota$, $\Omega$, and $\omega$ in Fig. 5). However, this loss of accuracy is only due to numerical noise, since for small values of $J_{2}$, the perturbing effect is so small that the hyperbolic elements tend to be constants. The difference between their values at the reception and at the transmission remains at the level of the numerical noise. For instance, for $N_{0}(R)=10^{-5}$, we have $p(f_{1})\sim 10^{8}$, moreover the analytical solution predicts $\Delta p_{\text{ana}}\sim 10^{-5}$ for $J_{2}=10^{-6}$, which represents a change in $\Delta p_{\text{ana}}/p(f_{1})\sim 10^{-13}$, i.e. one order of magnitude beyond the numerical double precision.

Secondly, we also notice an important feature linked to the accuracy of the solutions. Indeed, by changing the order of magnitude of the refractivity we also change the order of magnitude of the relative error on the total change in the value of all parameters. This feature reveals the approximation at leading order in $1/e$, which has been assumed in order to simplify the integration of the perturbation equations. In fact, while changing the order of magnitude of the refractivity value, which corresponds to a change in $\alpha$, we also modify the order of magnitude of the eccentricity. From Eq. (22), we notice that increasing $\alpha$ makes the eccentricity smaller, which decreases the accuracy of the solutions which are derived at leading order in $1/e$. For very high refractivity (e.g. $10^{-1}$ which corresponds here to a value of eccentricity around $10$), it is seen (cf. purple curves in Fig. 5) that the change in the hyperbolic elements, as given by the analytical expressions, is accurate at the level of $\sim 10\%$. For usual typical values of the refractivity ($\sim 10^{-4}$), the analytical solutions in Eqs. (121a)–(121e) are found to be accurate at the level of $\sim 0.01\%$.

For the geometrical length, the total light time and the refractive bending, the relative error evolves linearly with $J_{2}$ in the log-log plot. This behavior is due to the fact that the hyperbolic contribution is included within the analytical computation of $s$, $t$, and $\psi$. Indeed, for hyperbolic elements, as discussed above, when $J_{2}$ tends to be null, the computation of $\Delta\mathbf{C}_{\text{num}}$ generates a lot of numerical noise since the hyperbolic elements are constants when the perturbation vanishes. In the case of $s$, $t$, and $\psi$, when $J_{2}$ goes to zero, the computation of $\Delta\mathbf{C}_{\text{num}}$ does not generates any numerical noise since a non-null hyperbolic contribution remains when the non-hyperbolic one vanishes. In order to prove that the non-hyperbolic evolution of $s$, $t$, and $\psi$ is well described by Eqs. (121f)–(121h), one can compute the relative error considering only the hyperbolic contribution in $\Delta\mathbf{C}_{\text{ana}}$. The computation reveals that, the relative error computed without the non-hyperbolic contributions is between one and two orders of magnitude larger than the relative error computed with the non-hyperbolic contribution.

The internal accuracy of our solutions describing the evolution of the geometrical length, the light time and the refractive bending can be assessed by looking at Fig. 5. For very high refractivity (e.g. $10^{-1}$), the maximum relative error is reached for the maximum value of the oblateness ($J_{2}=10^{-1}$). For instance, considering the hyperbolic and the non-hyperbolic contributions, we are able to provide analytical solutions which are accurate at the level of $\sim 0.01\%$ for both $s$ and $t$, and at the level of $\sim 1\%$ for $\psi$. For a typical value of the refractivity ($\sim 10^{-4}$), and for a value of Jupiter or Saturn’s $J_{2}$ ($\sim 10^{-2}$), we might expect errors at the level of $\sim 10^{-6}\%$ for both $s$ and $t$, and at the level of $\sim 0.001\%$ for $\psi$. For the same values of the refractivity and $J_{2}$, if we only consider the hyperbolic contribution, the relative error grows to $\sim 10^{-5}\%$ for both $s$ and $t$, and to $0.1\%$ for $\psi$.

VII Conclusion and perspectives

In this paper we have demonstrated, in Sec. III, that the equations of geometrical optics possess an exact solution, referred to as the reference solution, when one assumes the hydrostatic equilibrium together with the constancy of the variation of the index of refraction with the Newtonian potential. Not surprisingly, since the problem reduces to a central field problem, the solution is found to be a conic-section similar to the solution of the well-known two-body problem in celestial mechanics. For the cases of interest, the reference solution is an hyperbola which is completely described with some constants of integration called hyperbolic elements. These elements describe the geometry of the light beam, namely the shape and the spatial orientation of the hyperbola. We have shown that the hyperbolic elements are related to the refractive response of the medium to the gravitational potential which is parametrized by the derivative of the index of refraction with respect to the potential.

The reference solution is found to be applicable in situation involving light propagation in planetary atmospheres, and thus, to be well suited for future applications to GB observations or AO experiments. In Sec. III.4, we provided short indications on how applying it in the context of real observations, but its practical application to real space missions’ data will be discussed in a successive paper. However, the main feature of this solution is not its direct application, but the fact that it provides a comprehensive framework that can be extended to carry out analytical studies in the case where the index of refraction has generic dependencies.

Indeed, based on the method of variation of arbitrary constants, in Chap. IV we have converted the equation of geometrical optics as length rate of change in the hyperbolic elements. These new equations are shown to be the optical counterpart of the perturbation equations in celestial mechanics. When no assumption is made on the order of magnitude of the change in the refractive profile, they are perfectly equivalent to equations of geometrical optics. If they are less compact, they present the advantage of providing a comprehensive alternative description to the path of light rays. They describe quantitatively the departure from the hyperbolic path due to additional terms besides the hyperbolic contribution into the index of refraction.

In the case where the perturbing gradient can be assumed small with respect to the hyperbolic one, we have shown in Sec. IV.5, how to approach analytically the solutions of the perturbation equations. This procedure seems to be generic enough to easily handle the first order effects due to any kind of perturbation besides the hyperbolic contribution.

To highlight the capabilities of this formalism, in Chap. V we have analyzed different sources of perturbing gradients. For instance, we have studied the effects on the light path due to, the non-linear dependencies with the Newtonian potential, the centrifugal potential and the light dragging effect due to the rigid rotation of the atmosphere, the non-spherical gravitational potential of the planet, the tidal atmospheric bulge raised by the presence of a massive perturbing body, and the horizontal gradients of refractivity inside the atmosphere. These examples of utilization are just a non-exhaustive list showing the possibility offered by the perturbation equations formalism applied to geometrical optics.

Finally, in Sec. VI, in order to assess the validity of, i) the perturbation equations, and ii) the method of integration, we have compared the analytical solution derived in the context of a quadrupolar axisymmetry in the gravitational potential of the planet with its numerical counterpart determined from a numerical integration of the equations of geometrical optics. We carried out several numerical integrations for different values of the refractivity and the $J_{2}$ parameter. For each one of them, we have compared the changes in the hyperbolic elements as provided by the numerical integration to those obtained from the analytical solution. By doing this, we have been able to assess the accuracy of the analytical solution. We have shown that for standard refractivity values ($\sim 10^{-4}$), the relative error in the total change in the hyperbolic elements is of the order of $0.01\%$ and is independent of the value of $J_{2}$. For the other elements like the geometrical length, the time-delay and the refractive bending, the relative error evolves with $J_{2}$ as shown in Fig. 5. For Jupiter or Saturn’s typical value of $J_{2}$ ($\sim 10^{-2}$), the relative error is of the order of $10^{-6}\%$ for both the length and the time, and $0.001\%$ for the bending. This represents a really good agreement considering the simplicity of the final solution. Indeed, by looking at the complete analytical expression of the light time, we notice that it is just given by the sum of the hyperbolic [cf. Eq. (93)] and the non-hyperbolic contribution [cf. Eq. (121g)].

The main original contribution of this paper is the reformulation of the equations of geometrical optics into a set of perturbation equations [cf. Eqs. (54)], describing the evolution of an osculating hyperbola along the light trajectory. This reformulation is very convenient for a use inside planetary atmospheres since the hyperbolic trajectory has proved to be an exact solution for spherical symmetry. However, an other possibility for the reformulation of the equations of geometrical optics could have been to use the straight line as an osculating solution instead of the hyperbola, offering in this way the possibility of analytically studying other symmetries involving small refractivity (e.g. for application to occultations by the Io plasma torus Bagenal and Sullivan (1981); Phipps and Withers (2017); Phipps et al. (2018)). This is also something which can can be easily achieved from Eqs. (54) by taking the limit $\alpha\rightarrow 0$.

In the context of the realization of the IERS reference systems, the perturbation equations together with the method of integration, could provide an efficient tool for improving existing models of the Earth tropospheric delay. In particular, an accurate modeling of horizontal gradients is of prime importance for the determination of station coordinates which in turn define the scale and the origin of the ITRF. However, existing models of the tropospheric delay always stick to the spherical symmetry assumption, and only few of them consider horizontal gradients (Chen and Herring (1997) for VLBI, Hulley and Pavlis (2007) for SLR). In addition, even for those which consider horizontal gradients, the integration across the atmosphere is always done numerically assuming spherical layers made of constant index of refraction. The perturbation equations could easily tackle these issues by providing e.g. the contribution to the time-delay due to the atmosphere’s oblateness and the one due to horizontal gradients, as discussed in Secs. V.4 and V.6.

In the context of AO experiments, the perturbation equations allow to assess for the first time a clear description of the photon path, the light time or the refractive bending into an oblate atmosphere. As mentioned in the introduction, this question is usually tackled using numerical ray-tracing techniques Schinder et al. (2015), and only a few studies Lindal (1992) aimed at exploring analytically this problem. Yet, the main advantage of analytical formulations is twofold. First, they help to acquire a physical understanding of the phenomenon, and secondly, they do not require high computation time while analyzing observational data. If the solution proposed by Lindal (1992) is a first analytical approach, it fails to provide a clear physical understanding of the photon path or the atmospheric time-delay due to the atmosphere’s oblateness. The reason is due to the fact that it consists in a Taylor series expansion around a local point along the photon path, and in addition it is expressed in terms of Cartesian coordinates which might be somehow abstract. Because of this Taylor expansion, the solution cannot be employed for high refractivity and is only defined in the surroundings of the expansion point. With the perturbation equations, we have solved all these difficulties since the solution is available all along the light path trajectory and can also be applied for high refractivity. In addition, the solution provides a simple geometrical interpretation by fitting at any length along the ray an hyperbolic trajectory. Another advantage which is worth mentioning is that the perturbation equations allows to express directly the light time as well as the refractive bending expressions which are needed for computing range and range-rate observables. Finally, let us mention that if analytical studies are not as precise as purely numerical ones, the formalism of the perturbation equations has been shown to be really accurate (errors of one part in $10^{8}$ and $10^{5}$ between the numerical and the analytical predictions on the light time and the refractive bending, in the presence of an oblate atmosphere characterized by $J_{2}=10^{-2}$ and $N_{0}(R)=10^{-4}$, respectively). In addition, we have demonstrated that the formalism presented in this paper offers a large flexibility since it is able to handle a wide range of perturbations including light dragging effects caused by atmosphere’s winds. Thus, a complete and purely analytical method, only based on analytical solutions derived in Chap. V, could be developed in order to process real AO data in the context of past, current, and future space missions. Our next opportunity, which actually motivated this work, is ESA’s L-class mission JUICE (JUpiter Icy moons Explorer) Grasset et al. (2013); Witasse (2016); here a radioscience experiment named 3GM (Gravity and Geophysics of Jupiter and the Galilean Moons) will take advantage of a careful Jupiter system tour design Boutonnet and Schoenmaekers (2016), offering both frequent Jupiter moon’s flybys Lari (2018); Dirkx et al. (2017) and radio occultation opportunities by Jupiter’s oblate atmosphere.

Acknowledgements.

The Authors are grateful to the University of Bologna and to Italian Space Agency (ASI) for financial support through the Agreement 2013-056-RO in the context of ESA’s JUICE mission. A.B. is grateful to F. Mignard from Observatoire de la Côte d’Azur for interesting discussions about atmospheric and astronomic refraction. A.B. is also grateful to A. Hees from Paris Observatory for valuable comments about a preliminary version of this manuscript.

Appendix A Hyperbolic path

As discussed in Chap. III, when the index of refraction is only function of the height of the ray (spherically symmetric assumption), all the light path is contained in the propagation plane. Using this fact, we derive, in this chapter, a complete solution to Eqs. (1), (4), and (5). We start by determining the shape of the light trajectory (i.e. the light path), then we focus on the time-delay.

A.1 Shape of the light path

We wish to write the equations of geometrical optics in the propagation frame and in polar coordinates by making use of the polar basis. The vectors h and $\bm{\nabla}\mathscr{S}$ are decomposed as

[TABLE]

where a dot denotes the differentiation with respect to $s$. Moreover, the left-hand side of Eq. (4) is given by

[TABLE]

The gradient of the index of refraction can also be expressed in the same coordinate system. With the help of Eq. (11) and the monopole approximation [cf. Eq. (16)], the right-hand side of Eq. (4) is given by

[TABLE]

The term $g$ represents the magnitude of the local acceleration (monopole contribution) experienced by the media which makes up the atmosphere of the central planet. It is explicitly given by

[TABLE]

It is obvious from the presence of $\alpha$ in Eq. (66), that the light ray experiences the gravitational pull via the media in which the ray propagates through.

From Eqs. (65) and (66), it is seen that the absence of a transverse component (along $\hat{\bm{\tau}}$) in the gradient of the index of refraction implies that $nh^{2}\dot{\theta}$ is a conserved quantity. From Eqs. (12) and (64) we immediately deduce

[TABLE]

Inserting this expression into the radial part of Eq. (65), we find

[TABLE]

We change the independent variable from $s$ to $\theta$ adopting the convention that a prime denotes the differentiation with respect to $\theta$. We also adopt $u\equiv 1/h$ as a convenient substitute for $h$ and derive a differential equation for it. Making all these substitutions into Eq. (69), we quickly arrive to

[TABLE]

with

[TABLE]

a constant parameter. The general solution to the simple Eq. (70) is a conic section with origin at the focus. Returning to the original radial variable the spatial solution is

[TABLE]

in which $e$ is an arbitrary constant of integration called the eccentricity of the conic. We will adopt the convention that the eccentricity is a positive quantity. The two other parameters, $p$ and $f$, are defined such that

[TABLE]

They are respectively known as the semi-latus rectum and the true anomaly. The last parameter, $\omega$, is also a constant of integration and its role is better understood when $h$ is minimal, that is to say when $\theta=\omega$. It is seen that $\omega$ defines the argument of the closest approach which is the minimal separation between the path of the ray and the center of the reference frame. This angle $\omega$ is completely determined by the initial conditions of the problem. For instance, if we call $(h_{1},\theta_{1})$ the components of the position vector of the transmission point, then from the solution in Eq. (72), we end up with

[TABLE]

where the sign is taken to be positive when $-\pi/2\leq\theta_{1}<\pi/2$ and negative elsewhere.

The minimal separation is also determined from Eq. (72) evaluated at the closest approach, such that

[TABLE]

In the following, we will consider that the signal was emitted or received (or both) from infinity, in a space region well beyond the planetary atmosphere where the index of refraction can be considered to be unity (or more generally, constant). This restriction let to avoid periodic light path. Consequently, we will restrict our attention to non-periodic solutions ($e\geq 1$). Following this, the only possibility for the eccentricity value is $e>1$, which corresponds to hyperbolic trajectories for the light path. Therefore, we can introduce a new constant known as the semi-major axis

[TABLE]

which is a negative quantity.

Non-periodic solutions, which are hyperbolic trajectories, also imply the introduction of a new important angle, $\Delta f_{\infty}$, describing the net change between the two asymptotes’ directions. Based on Eq. (72), it can be seen that $h$ goes to infinity when $\cos\kappa f$ tends to $-1/e$, which gives rise to two asymptotic solutions for $\kappa f$. Let call $f_{\text{in}}^{\infty}$ the negative solution and $f_{\text{out}}^{\infty}$ the positive one. If we introduce $\Delta f_{\infty}=f_{\text{out}}^{\infty}-f_{\text{in}}^{\infty}$, we immediately deduce

[TABLE]

Since, we focus on hyperbolic trajectories $(e>1)$, we will always have $\Delta f_{\infty}<2\pi$.

In AO experiments, the geometry is such that $\Delta f_{\infty}$ can be geometrically related to the refractive bending (indeed from Fig. 1, we deduce $\Delta f_{\infty}=\pi+\epsilon$) which is itself directly in relation with the changes in the Doppler measurements, at first order (see Eq. (148) or Fjeldbo and Eshleman (1965, 1968)). Thus, Eq. (77) can be used to infer the eccentricity from the Doppler frequency shift measurements.

Once the evolution of the height is totally determined, we can focus on the evolution of the tangent to the ray [cf. Eq. (64b)]. We thus need to determine the expressions of the radial and the transverse components; the later depends on the angular velocity. Invoking Eqs. (72) and (68) and differentiating them making use of Eqs. (73), reveals that

[TABLE]

The evolution of the tangent to the ray is provided after inserting those relationships together with Eq. (72) into Eq. (64b).

A.2 First integrals of the light path

Until now, we are missing an expression for the eccentricity and the argument of the closest approach in term of fundamental constants of the problem, such as K. We saw e.g. that the semi-latus rectum is linked to the conservation of the magnitude of the impact parameter [cf. Eq. (73a)]. In this section, we look for additional first integrals of the light path, which then could let to express the constant of integration in a fundamental way.

The first fundamental constant can be obtained by multiplying both side of Eq. (69) by $n^{2}\dot{h}$. Then, one can see that each term is a total differential with respect to $s$. In addition, one can recognized the squared components of (64b) in the integrated equation, and we end up with

[TABLE]

where $E$ is a constant of integration. That constant would be similar to the energy in classical mechanics. Indeed, by computing the square of $\bm{\nabla}\mathscr{S}$ from Eqs. (78) and (64b), we end up with

[TABLE]

where we have used Eq. (76). Then, insertion of Eq. (80) into (79), leads to

[TABLE]

which possesses the same form as the energy in the Kepler problem. From the definition of $a$, the right-hand side of the equation is directly seen for being constant and positive. Therefore, $E$ is conserved along the all light ray trajectory.

As a general remark, let us notice that Eq. (79) can be put under a more fundamental form. Indeed, by making use of the definition of $n$ [cf. Eq. (17)], we deduce

[TABLE]

where

[TABLE]

is known as the Fresnel coefficient [$\eta$ has been previously defined in Eq. (17)]. This expression of $E$ is more fundamental than Eq. (79) since it does not require a prior knowledge of the index of refraction. Indeed, the constancy of $E$ can be inferred using only Eqs. (1) and (4).

From Eq. (81), it is seen that the semi-major axis is linked to the conservation of the dimensionless energy parameter, $E$. At the same time, the constancy of the eccentricity is assured, because of the definition (76). So, in that sense the eccentricity is linked to both the conservation of energy and the conservation of the magnitude of the impact parameter. Eqs. (81), (76), and (73a) can be used to infer the eccentricity

[TABLE]

In addition to $E$, and K, we can determine an other fundamental constant of the problem. This last one, is equivalent to the eccentricity vector of celestial mechanics and comes from the very peculiar form of the index of refraction, which is linear with the gravitational potential making it evolving as $1/h$. However, because $\kappa$ is different from unity, the eccentricity vector does not show up as easily as in celestial mechanics. It is given by the following relationship

[TABLE]

where we have introduced the vector A which possesses the following components in the polar basis

[TABLE]

The vector A is expressed in terms of fundamental quantities after determining $f$ from Eq. (72), and expressing $\hat{\bm{\rho}}$ and $\hat{\bm{\tau}}$ in terms of h and K. The eccentricity in Eq. (85), must be expressed in terms of $E$ and $K$ thanks to Eq. (84).

The constancy of e can be demonstrated by substituting Eqs. (78) into (64b), then inserting the result together with Eq. (68) into (85), which finally gives with Eq. (15)

[TABLE]

Since $e$ and $\omega$ are two constants of integration, e is conserved during the propagation of light. Moreover, it is seen that the constancy of $\omega$ is linked to the conservation of the eccentricity vector. It always points towards the closest approach and its magnitude is constant and equal to the eccentricity of the light path.

The important point of this section is summarized recalling that the light path is totally contained inside a fixed plane which remains orthogonal to the impact parameter vector when the refractive profile is purely radial. In addition, we found that the shape of the trajectory is described by an hyperbola when the refractive profile is given by Eq. (17). That shape is parametrized thanks to two geometrical parameters, $p$ and $e$, and the orientation in the propagation plane is located thanks to one angle, $\omega$. Those three parameters can totally be expressed [cf. Eqs. (73a), (84), and (87)] in terms of fundamental quantities of the problem, such that $K$, $E$, and e.

A.3 Kepler-like problems

In the previous section, we have defined the shape of the trajectory as well as its orientation in the propagation plane in terms of constant parameters. In order to completely solve the problem in the plane, we need to determine a single location along the light path given a geometrical length or a light time.

The Eq. (78b) can be used to complete the description of the light path by providing a unique relationship between $f$ and $s$. This can be accomplished by integrating Eq. (78b) as

[TABLE]

where $S$ is an arbitrary constant of integration representing the traveled geometrical length to the closest approach. Additionally, we can also derive a unique relationship between $f$ and $t$, by making use of Eq. (5) in order to turn Eq. (88) into

[TABLE]

where $T$ is an arbitrary constant of integration representing the light time till the closest approach.

Eqs. (88) and (89) are the analogous of the Kepler equation of celestial mechanics. They can be exactly solved by introducing a new variable, $F$, known as the hyperbolic anomaly. The change from the true to the hyperbolic anomaly, and conversely, is given by

[TABLE]

Eqs. (90) allow us to introduce the following relations which relate the differentials of $f$ and $F$

[TABLE]

The hyperbolic anomaly can be employed instead of the true anomaly as a convenient substitute. Making all these substitutions we quickly arrive to

[TABLE]

The last equation is perfectly equivalent to Eq. (88), however it possesses the main advantage of being easily integrable. Then, making use of Eqs. (17) and (92a) in order to express the refractive index in terms of $F$, we find the following analogous to the well-known Kepler equation

[TABLE]

in which the left-hand side is the geometrical length mean anomaly.

Similarly, expressing Eq. (92c) in term of the light time instead of length path with the help of Eq. (5), let one to directly integrate Eq. (89), such that

[TABLE]

in which the left-hand side is the light time mean anomaly

From these two Kepler-like equations, it is seen that the light time [cf. Eqs. (94)] is not just the geometrical length [cf. Eqs. (93)] multiplied by the inverse of the speed of light in a vacuum. Indeed, it also involves some additional terms which account for the apparent change in the speed of propagation of the electromagnetic waves. Consequently, the geometrical length mean anomaly considers a ray traveling at the speed of light along the bent trajectory while the time mean anomaly is related to the phase path and takes into account the apparent change in the speed of the waves along the actual path of photons.

A.4 Argument of the refractive bending

In this section, we take advantage of the introduction of the hyperbolic anomaly to define a very important geometrical relationship between the true anomaly and the argument of the refractive bending, $\psi$. As seen in Fig. 1, $\psi$ is the angle between the orthogonal direction to the line of nodes and the tangent to the light ray.

In the context of AO experiments, this angle plays an important role. Indeed, from Fig. 1, it might be seen that the net change in $\psi$ between the transmitter (labeled 1) and the receiver (labeled 2) is exactly the refractive bending, $\epsilon\equiv\psi_{2}-\psi_{1}$. Moreover, as shown in Eq. (148) and in Fjeldbo et al. (1971) or Eshleman (1973), the refractive bending can be use to infer the Doppler frequency shift which is the main observable for AO experiments.

The differential relationship between $\psi$ and $f$ can be inferred from Fig. 1 which allows us deduce the following relationship

[TABLE]

We remind that the angle $\phi$ has been previously defined in Eq. (13) as the angle between $\hat{\bm{\rho}}$ and s. Hence, from Eq. (64b) we immediately deduce

[TABLE]

Then, differentiating Eqs. (13) and (95), and making use of Eqs. (96), we obtain the following relation

[TABLE]

which can be written as

[TABLE]

Making use of the relationships in Eqs. (92), we change the variable of integration from the true to the hyperbolic anomaly, and after some algebra we integrate the previous expression as

[TABLE]

in which we have introduced

[TABLE]

Since Eq. (90c) provides a unique relationship between $f$ and $F$, Eq. (99) can be used to find the argument of the bending angle knowing the true anomaly.

Appendix B Mathematical details for applications

In this chapter, we derive the change in hyperbolic elements for different contributions in the index of refraction.

B.1 Non-linearity with $U_{0}$

We first simplify the problem by assuming that the generalized potential reduces to the spherical part of the Newtonian potential ($\Phi=U_{0}$). Consequently, Eq. (61) and the perturbing gradient now read

[TABLE]

We can directly see that the non-null components of the perturbing gradient are the radial and the change in the index of refraction

[TABLE]

Eqs. (102) can then be inserted into Eqs. (55) and the integrand is expanded in power of $1/e$ (see Sec. IV.5). The change in the hyperbolic elements is given for any $k$ (with $k\geq 2$) by the following expressions

[TABLE]

where we have kept the leading order in $1/e$, and where $f_{2}\geq f_{1}$.

Considering e.g. the quadratic dependence with the gravitational potential in the index of refraction profile ($k=2$), we see that the evolution of the hyperbolic elements all along the light path is given by

[TABLE]

All those equations are defined within a constant which is obviously given by the initial condition of the problem at the level of the transmitter, so Eqs. (104) are valid for all $f\geq f_{1}$.

B.2 Steady rotating atmosphere

The term $\Phi_{C}$ in Eq. (60) is the centrifugal potential due to proper rotation of the fluid body and is given by Kaula (1966); Poisson and Will (2014)

[TABLE]

At linear order in the generalized potential, the index of refraction and its gradient are expressed in the frame comoving with the fluid as

[TABLE]

Once $\hat{\textbf{X}}^{\prime}$ and $\hat{\textbf{Y}}^{\prime}$ are projected into the polar basis, we obtain the following components

[TABLE]

Inserting Eqs. (107) into Eqs. (55), expending them in power of $1/e$ and performing the integration, allows one to find, at first order in $1/e$, that the centrifugal contribution impacts all the hyperbolic elements, such

[TABLE]

All these solutions are defined within a constant which is given by the initial condition of the problem at the level of the transmitter, so Eqs. (108) are valid for any $f\geq f_{1}$.

We remind that the previous results only account for the contribution of the centrifugal potential and does not consider the dragging effect due to the moving medium. To quantify it, let us express Eq. (58b) as

[TABLE]

in which, we have assumed a steady rotating atmosphere such that $\mathbf{v}=\textbf{w}\times\mathbf{h}$. If we express the dragging contribution in the polar basis, we obtain the following non-null components

[TABLE]

where $N$ is the refractivity defined as $N=n-\eta$. We notice the absence of the $\mathcal{N}$-component of the perturbation.

The effect of the fluid rotation on the light propagation is assessed by inserting those components into Eqs. (55), by keeping the first order in $1/e$, and then by performing the integration over the true anomaly

[TABLE]

Once again let us precise that these equations are defined within a constant which is determine with the transmitter’s position. Eqs. (111) are defined for $f\geq f_{1}$.

B.3 Axisymmetric gravitational potential

A convenient way to handle non-spherical gravitational potential is to deal with either the spherical harmonic expansion Kaula (1966); Murray and Dermott (2000) or the symmetric and trace free (STF) tensors formalism Hartmann et al. (1994); Poisson and Will (2014). Using e.g. the later one, it can be shown that the expansion of the non-spheric part of the gravitational potential is given by

[TABLE]

for $l\geq 2$. We recall that $G$ is the gravitational constant, $I^{<L>}$ is the multipole moment STF tensors of the mass distribution, and $\partial_{<L>}$ is the STF tensor made with the partial derivatives with respect to h. We refer to Hartmann et al. (1994) and Poisson and Will (2014) for notations and definitions about the use of the STF tensors.

We remind that the unit vector $\hat{\mathbf{Z}}$ has been chosen in Sec. V.1 for being aligned with the rotation axis of the extended body. We focus on the special case where the spatial changes in the rotation axis are supposed to be negligible during the light propagation event. Thus, we consider that the $\hat{\mathbf{Z}}$-direction is spatially fixed. In addition, we consider a fluid body in hydrostatic equilibrium rotating at a constant rate. Then, the resulting gravitational field is considered for being axisymmetric and independent of the longitude Poisson and Will (2014). This symmetry imposes the multipole moments STF tensors to be proportional to $\hat{Z}^{<L>}$ which is the STF tensor formed with the $\hat{\mathbf{Z}}$ unit vector. Hence, we have the relation

[TABLE]

where $M$ is the total mass of the extended body, $R$ its equatorial radius, and $J_{l}$ are its unitless multipole moments.

Making use of basic properties about STF tensors Hartmann et al. (1994); Poisson and Will (2014), such

[TABLE]

where $\chi=\hat{\bm{\rho}}\cdot\hat{\mathbf{Z}}$ and $P_{l}(\chi)$ are the well-known Legendre polynomials of degree $l$, we can rewrite the perturbing gravitational potential in Eq. (112) in terms of the multipole moments $J_{l}$. Keeping the linear part in $\Phi$ in the development of the index of refraction [cf. Eq. (61)] with $\Phi=U_{0}+\sum_{l}U_{l}$, we determine an expression for $n$ in term of the non-spherical part of the gravitational potential

[TABLE]

This profile differs from the spherical one [cf. Eq. (17)] by terms of the order $J_{l}$ when the height of the ray is at the same order of magnitude than the equatorial radius. We can now, determine the $j$th component of the gradient of the index of refraction which is given by

[TABLE]

Differentiating Eq. (112) and using the previous properties about STF tensors [cf. Eq. (114c)], we deduce

[TABLE]

In the polar basis the components of $\textbf{f}_{\text{pert}}$ are

[TABLE]

Usually, for modestly deformed bodies (e.g. the planets in the Solar System), the most important term in the developments (115) and (117) is the one proportional to the quadrupolar moment $J_{2}$. This parameter measures the rotational flattening due to the proper rotation of the extended body. From now, we restrict our attention to that single parameter.

Using the definition of Legendre polynomials

[TABLE]

with Eqs. (19), we can determine the change in the index of refraction as well as the components of the perturbing gradient

[TABLE]

Inserting those components into Eqs. (55) and applying Eq. (56), allows one to find the following estimations which are given at leading order in $1/e$

[TABLE]

These equations are defined within a constant which is determined from the initial conditions at the level of the transmitter. They are valid for all $f\geq f_{1}$.

B.4 Tidal potential

In the fluid rotating frame, the perturbing tidal potential can be expressed as Poisson and Will (2014)

[TABLE]

where $\mathcal{E}_{L}(t)$ are the tidal moments being time dependent STF tensors. They are defined by

[TABLE]

where the gravitational potential created by the external bodies $U_{\neg P}$ is differentiated $l$ times with respect to $\mathbf{h}$ and the result is evaluated at the central planet’s center-of-mass (i.e. $\mathbf{h}=\bm{0}$).

Following Poisson and Will (2014), $U_{\neg P}$ is given by

[TABLE]

in which the integration must be performed in the vicinity of the external bodies $A$ where $\rho_{A}(t,\mathbf{x})$ is supposed to be non-null. Thus, the dummy variable of integration will ranges over values close to $\mathbf{r}_{A}(t)$, and we can introduce $\mathbf{x}=\mathbf{r}_{A}(t)+[\mathbf{x}-\mathbf{r}_{A}(t)]$. Then, we are left with two vectors in the denominator, the first one, $\mathbf{h}+\mathbf{r}_{AP}(t)$ with $\mathbf{r}_{AP}(t)\equiv\mathbf{r}_{P}(t)-\mathbf{r}_{A}(t)$, representing the separation between the center-of-mass of $A$ and a point $\mathbf{h}$ at which the external potential must be evaluated, and the second one, $\mathbf{x}-\mathbf{r}_{A}(t)$, representing the separation between the center-of-mass of $A$ and a point lying in the vicinity of $A$. Therefore, because $\mathbf{h}+\mathbf{r}_{AP}(t)$ is usually of the order of the inter-body separation, it is appropriate to express the denominator under the integrand as a Taylor series expansion. After some algebra, we arrive to

[TABLE]

We recall that $I_{A}^{<L^{\prime}>}$ means the multipole moment STF tensors of body $A$ and that the derivative must be taken with respect to $\mathbf{h}$ and then be evaluated at $P$’s center-of-mass (i.e. ${\mathbf{h}=\bm{0}}$). This notation is shorten noticing that $\partial_{L^{\prime}}|\mathbf{r}_{AP}+\mathbf{h}|^{-1}=\partial_{L^{\prime}}r_{AP}^{-1}$, where the derivative is taken at the extremity of the vector $\mathbf{r}_{AP}(t)$.

From Eq. (125), (122), and (123), we can define the perturbing external tidal potential

[TABLE]

That expression is very general since it takes into account the tidal effects due to the total gravitational potential (all multipole moments included) of all the external bodies acting on the total gravitational potential of the central planet.

At linear order in the generalized potential, the refractive profile and the $j$th component of the gradient of the non-constant change in the index of refraction are explicitly written as

[TABLE]

For the application that we are to consider now, we admit only one perturbing body, $A$, and we define the separation vector between $P$ and $A$ as $\mathbf{r}(t)\equiv-\mathbf{r}_{AP}(t)$. We consider that the central planet is close to spherical symmetry, also the coupling between its higher multipole moments $I^{<L>}$ and the other STF tensors can safely be neglected. In addition, we consider that the perturbing body is a satellite of the main planet $P$, and that $M_{A}\ll M$, with $M$ the mass of the central planet. Therefore we can safely consider only the monopole contribution of $A$, i.e. $l^{\prime}=0$. In addition, assuming that inter-body distance is much more larger than the typical scale of the ray closest approach to $P$ (we are interested in the light propagation in the atmosphere of $P$), it is sufficient to keep only the leading quadrupole term in the development, i.e. $l=2$. For a matter of simplicity, we also consider that the satellite follows a circular orbit of radius $r$ which lies in the equatorial plane of the central planet. Then, in the reference frame, the direction of the position vector of the satellite is given at any time $t$ by

[TABLE]

where the mean motion of the satellite around the central body is given by Kepler third law of motion $N\simeq(\mu/r^{3})^{1/2}$ (For this application, $N$ is not the refractivity anymore). Making use of Eqs. (19a)-(19c), the orbital motion of the satellite is given, in the polar frame, by

[TABLE]

in which we have defined $\lambda(t)\equiv Nt-\Omega$ for being the longitude of the perturbing body as measured from the intersection between the propagation plane and the equatorial plane (longitude of the node); this angle remains within the equatorial plane of $P$ since the perturbing body describes an equatorial circular orbit.

Following all the simplifications, the change in the index of refraction and its gradient now read

[TABLE]

Making use of properties about the STF tensors [cf. Eqs. (114)], the multi-index derivatives now read

[TABLE]

In those expressions, the components of the unit vector directed along $\mathbf{r}$ are noted $\hat{r}_{j}\equiv r_{j}/r$ and can be inferred from Eq. (129) in the polar basis.

We can now express the components of the perturbing gradient and the non-constant change in the index of refraction

[TABLE]

Here, we have introduced the unitless vector $\bm{\mathcal{A}}$ and the unitless parameter $\mathcal{A}_{\mathcal{N}}$ to simplify the notations. They are explicitly given by

[TABLE]

and once expressed in term of the fundamental angles of the problem they read

[TABLE]

One can now insert Eqs. (132) into Eqs. (55). Before integrating, we are to further simplify the problem by assuming a ray traveling in the equator of the central planet and passing through the bulge raised by the perturbing body. In such a case, $\iota=0$ and the equation can be further simplified. However, as discussed in Sec. IV.3, when the inclination is null the longitude of the node and the argument of closest approach are undefined. Therefore, instead of considering $\omega$ and $\Omega$ separately, we consider the longitude of the closest approach, $\varpi=\omega+\Omega$, which is measured from the $\hat{\textbf{X}}$-direction. Keeping the leading order in $1/e$, and integrating over the true anomaly leads to the following set of equations describing the change in hyperbolic elements following a gravitational tidal stress

[TABLE]

These equations are defined within a constant which is determined from the initial conditions at the level of the transmitter. They are valid for all $f\geq f_{1}$.

In order to account for the non-elastic response of the fluid following a tidal stress, it is shown in Poisson and Will (2014) that the first order effects can be assessed by replacing $\mathcal{E}_{L}(t)$ by $\mathcal{E}_{L}(t-\tau)$ where $\tau$ is called the time-lag and is function of the Energy dissipation. A computation of $\mathcal{E}_{L}$ at the instant $t-\tau$, changes the instantaneous position of the perturbing body [cf. Eq. (128)] into a delayed position, $\mathbf{r}(t-\tau)$. Most of the time, $\tau$ is a small quantity compared to orbital characteristic time scale. Consequently, the orbital velocity of the perturbing body is usually too small to change significantly it’s apparent position during the time interval $\tau$. In this condition, only internal process occurring on time scale faster than the orbital motion can significantly change $\mathcal{E}_{L}$. Among these process, the rigid rotation of the atmosphere appears to be a good candidate for fast rotators. The delayed position can therefore be approximated by

[TABLE]

where $w$ is the magnitude of the angular velocity which has been defined to be aligned with the normal to the equator of reference. The delayed unit vector is given by

[TABLE]

where the star refers to the non-elastic contribution which is obtained by simple identification with the previous equation.

Then, it is easily shown that the non-elastic part in the components of the perturbing gradient and the non-constant change in the index of refraction are given at first order in $\tau$ by

[TABLE]

As before, we have introduced the trigonometric coefficients $\bm{\mathcal{A}}^{*}$ and $\mathcal{A}_{\mathcal{N}}^{*}$ which are given by

[TABLE]

One can notice from the form of Eqs. (136), that

[TABLE]

Then, it is immediate to show that the change in hyperbolic elements due to the non-elastic response of the atmosphere to tidal stress is

[TABLE]

where C represents the vector of the hyperbolic elements.

B.5 Horizontal gradients

Let us consider the Earth centered frame $(\hat{\textbf{X}},\hat{\textbf{Y}},\hat{\textbf{Z}})$ with the $\hat{\textbf{Z}}$-axis directed through the North pole, the $(\hat{\textbf{X}},\hat{\textbf{Y}})$ plane superimposed within the Earth equator, and the $\hat{\textbf{X}}$-axis pointed through the intersection between the equator and the Greenwich meridian. We assume that the atmosphere is filled with a stationary medium with respect to the frame attached to the Earth, whence Eq. (4) is valid in $(\hat{\textbf{X}},\hat{\textbf{Y}},\hat{\textbf{Z}})$. A station on Earth can be located thanks to its spherical coordinates $(\varphi_{2},\lambda_{2},R)$, with $R$ the Earth radius, $\varphi_{2}$ the colatitude (measured from the North pole), and $\lambda_{2}$ the azimuth of the station121212Here, we consider the case of a receiving station but the same work can be applied for a transmitting station.. Any point along the light ray trajectory can be located through its spherical coordinates $(\varphi,\lambda,h)$.

We are going to consider a spherical atmosphere with an hyperbolic radial dependency in the refractive profile. In addition to the radial dependency, we assume an horizontal variation of the group refractivity above the transmitting site. At linear order, we can expend the index of refraction around the site of observation as it is done in Chen and Herring (1997) or Hulley and Pavlis (2007)

[TABLE]

where $n_{0}(h)$ is the hyperbolic index of refraction, and $n_{\text{NS}}$ and $n_{\text{WE}}$ are respectively given by $\partial n/\partial\varphi|_{h,\lambda}$ and $\partial n/\partial\lambda|_{h,\varphi}$. In other words, $n_{\text{NS}}$ represents the variation of the horizontal refraction along the North-South direction and $n_{\text{WE}}$ represents the variation of the horizontal refraction along West-East direction. Here, we aim at computing the effects on the hyperbolic elements due to $n_{\text{NS}}$ and $n_{\text{WE}}$.

The perturbative gradient is given by taking the gradient of $\delta n$, with $\delta n(\varphi,\lambda)=n(h,\varphi,\lambda)-n_{0}(h)$. Expressing the gradient in spherical coordinates we end up with

[TABLE]

where we have assumed that the radial variations of $n_{\text{NS}}$ and $n_{\text{WE}}$ are of second order. The two unit vectors $\hat{\textbf{u}}_{\varphi}$ and $\hat{\textbf{u}}_{\lambda}$ are part of the usual spherical vectorial basis which is given by the triad of vectors $(\hat{\textbf{u}}_{\varphi},\hat{\textbf{u}}_{\lambda},\hat{\bm{\rho}})$. It is seen that the perturbing gradient does indeed not possess any radial component if we neglect the radial variations of $n_{\text{NS}}$ and $n_{\text{WE}}$.

Now let us write the transverse and the normal components of the perturbing gradient by expressing $\hat{\textbf{u}}_{\varphi}$ and $\hat{\textbf{u}}_{\lambda}$ into the polar basis $(\hat{\bm{\rho}},\hat{\bm{\tau}},\hat{\bm{\sigma}})$.

For any arbitrary point, we can change from the spherical to the Cartesian coordinates thanks to two angles $\varphi$ and $\lambda$

[TABLE]

The radial direction can be compared to Eq. (19a) which reveals the following relationships

[TABLE]

We have introduced the trigonometric function $\mathcal{A}_{1}(f)$ which is defined in Tab. 1. From Eqs. (141) and (19) and by making use of Eqs. (142), we find

[TABLE]

where the $\mathcal{A}_{l}(f)$ coefficients are given in Tab. 1. Those equations allow us to express the perturbing gradient into the polar basis and allow us to derive its transverse and normal components.

At this level, we also need an equation for the change in the index of refraction in term of hyperbolic elements. That equation can be obtained from Eqs. (142). Indeed at first order in $f-f_{2}$ we deduce

[TABLE]

which both can be inserted into Eq. (139). We sum up the results into the following relations (the radial component is null)

[TABLE]

Finally, the effect on the hyperbolic elements is given by inserting Eqs. (145) into Eqs. (55). Before integrating we are to simplify the problem by assuming source tracking at the meridian of the site of observation. For this geometry, the propagation plane of photons is orthogonal to the plane attached to Earth’s equator ($\iota=\pi/2$). In addition, the azimuth of the source is the same as the azimuth of the receiving station, and because the inclination is $\pi/2$, the longitude of the node is the same as the azimuth of both the source and the station ($\Omega=\lambda$). Keeping the leading order in $1/e$, and integrating over the true anomaly leads to the following expressions

[TABLE]

These equations are defined within a constant which is determined from the initial conditions at the level of the transmitter. They are valid for all $f\geq f_{1}$.

Appendix C Doppler shift and refractive bending

In this chapter, we show how the Doppler frequency shift can be related to the refractive bending.

C.1 Deviation from vacuum contribution

We start by expressing Eq. (9) in an alternative way making use of the Doppler frequency shift in a vacuum, namely

[TABLE]

Inserting the definition (2) into Eq. (9), we arrive to

[TABLE]

where $\textbf{s}_{\text{vac}}$ is the light beam direction in a vacuum and $\delta\textbf{s}\equiv\textbf{s}-\textbf{s}_{\text{vac}}$ is the deviation of the actual direction of the ray with respect to $\textbf{s}_{\text{vac}}$. In vacuum, $\delta\textbf{s}=\bm{0}$, and $n-1=0$, thus $\nu_{2}/\nu_{1}$ reduces to $[\nu_{2}/\nu_{1}]_{\text{vac}}$.

For application within the Solar System, Eq. (147) can be simplified considering the small velocity approximation ($|\hat{\textbf{u}}|\ll 1$) together with the low refractivity limit ($n-1\ll 1$). The later condition imposes that the deviation with respect to the light ray direction in a vacuum is small, also $|\delta\textbf{s}|\ll 1$, and substitutions into Eq. (147) leads to the following approximation

[TABLE]

For an hyperbolic path, we can consider that $\hat{\bm{\sigma}}_{1}=\hat{\bm{\sigma}}_{2}=\hat{\bm{\sigma}}$ since K is constant in direction overall the light path. The angle $\beta$ represents the elevation of the photon path above from the straight line connecting the transmitter and the receiver (see Fig. 1).

For spherical symmetry, $\beta_{2}$ is proportional to $\beta_{1}$ by means of a trigonometric factor: $\beta_{2}=-\Lambda\beta_{1}$. This $\Lambda$-factor tends to be respectively zero or unity in the limiting cases where i) the receiver is at infinity or ii) the receiver and the transmitter are at equal distance from the center of mass ($\Lambda\rightarrow\{0;1\}$ when $h_{2}\rightarrow\{+\infty;h_{1}\}$). Moreover, from Fig. 1 it is seen that the total refractive bending can be expressed as $\epsilon\equiv\beta_{2}-\beta_{1}$, also we end up with

[TABLE]

This equation shows that in first approximation, the Doppler shift due to atmospheric effects is controlled by the total refractive bending, $\epsilon$. As discussed in appendix A, $\epsilon$ is defined as the net change in the argument of the refractive bending, $\psi$, between the transmitter and the receiver ($\epsilon\equiv\psi_{2}-\psi_{1}$), where the exact expression of $\psi$ has been given in Eq. (99).

C.2 Simplifications

For one-way Earth AO experiments Kursinski et al. (1997); Steiner et al. (1999), the index of refraction at the level of the two spacecrafts orbiting the Earth can be considered for being unity ($n_{1}=n_{2}=1$), also the previous equation reduces to

[TABLE]

when $\hat{\bm{\sigma}}=\hat{\bm{\sigma}}_{1}=\hat{\bm{\sigma}}_{2}$.

For one-way planetary AO experiments (see e.g. Fjeldbo et al. (1971); Eshleman (1973)), the index of refraction at the level of the spacecraft can be taken for being unity ($n_{1}=1$). In addition, the receiver (DSN antenna on Earth) can safely be considered at infinity, also $\Lambda\rightarrow 0$, and the simplified Doppler shift expression becomes

[TABLE]

Most of the time, the second term in the right-hand side of the equation is neglected since it can be absorbed by fitting a constant offset in the Doppler frequency shift measurements.