© ESO 2019

1 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 (as 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, atmospheric radio occultations (ARO), and atmospheric stellar occultations (ASO) are three examples of observational techniques requiring a careful analysis of atmospheric refraction. For GB observations, refraction is an accuracy-limiting factor and it has to be modeled in the analysis software; for both ARO and ASO, it is thoroughly modeled to infer the physical properties of the refractive medium. For all these techniques, 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). The atmospheric delay is usually evaluated at zenith and is computed with an analytical model of Marini & Murray (1973) or with the more recent model (also analytic) of Mendes & 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 atmosphere, which might be somehow inaccurate. (These gradients could arise at two different levels. Firstly, they could be generated by an atmosphere with a non-spherical shape. Secondly, they could result from a local horizontal variation of the physical quantities entering in the computation of the index of refraction such as the temperature.) This issue has only been partially overcome by Chen & Herring (1997) for very long baseline interferometry (VLBI), and by Hulley & 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 spherical. It is also worth mentioning that beyond the radial dependency of the refractive index, only numerical integrations are usually used for such GB observations. In ARO and ASO experiments, the basic idea is to establish electromagnetic links between a transmitter and a receiver when the former is being occulted by a planetary atmosphere as seen from the receiver. (For ARO experiments, a radio link is established, while for ASO experiments the transmitter is a distant star, so that the signal’s wavelengths usually range from near infrared to visible. In the following, we mainly focus on ARO experiments, but the entire discussion applies also to ASO experiments.) Conversely to GB observations, ARO 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 & Eshleman 1965, 1968). Later, it has been shown (Phinney & 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 denseatmosphere 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 occultations by Titan (Schinder et al. 2011, 2012), GPS/MET (global positioning system meteorology experiment) occultations by the Earth’s atmosphere (Kursinski et al. 1997; Steiner et al. 1999), and also ASO data for instance by Neptune’s (Roques et al. 1994) or Titan’s atmosphere (Sicardy et al. 2006). Later, with the Voyager 2 occultation by Neptune in 1989, the more recent Cassini mission orbiting Saturn, or ASO opportunities by Titan’s distorted atmosphere (due to the presence of zonal winds, Hubbard et al. 1993; Sicardy et al. 2006), 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 as 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 α, 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 ARO and ASO experiments. The geometric properties of the ray (as the photon path and the time-delay) are expressed in term of α 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 such as GB observations, ARO or ASO experiments.

2 Definitions and hypothesis

In this paper, we use the geometrical optics approximation which is characterized by neglecting the wavelength of electromagnetic waves (Landau & 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 section, we introduce the basic principles for determining the evolution of these quantities. We provide in Table A.1 a list of the nomenclature used hereafter.

2.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 (1)

d h is an elemental displacement vector along the ray and ds represents its length, which is defined as ds² = dh ⋅dh. 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 (2)

where is a real scalar function of the position and is usually defined as the optical path (Born & Wolf 1999). Here, the gradient is computed with respect to h, and n(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 (3)

which is the fundamental equation of geometrical optics. The function is also referredto 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 & Wolf 1999) (4)

This expression is given in a reference frame at rest with respect to the refractive medium. Equation (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 (5)

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

Equations (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.

2.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 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 (6)

Here, the integral represents the total light time accounting for the atmospheric effects and t_vac is the light time invacuum. 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) (7)

where ν₁ and ν₂ are the transmitted and the received frequencies respectively, u^α denotes the four-velocity vector, k_α is the covariant form of the four-wave vector which is tangent to light ray, and ûⁱ and are respectively given by ûⁱ = uⁱ∕u⁰, and (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 by following for instance Leonhardt & Piwnicki (1999). Indeed, for a static medium, it can be shown that (8)

Without loss of generality for the discussion, we neglect general relativity effects and focus our attention on atmospheric refraction only. A more generic treatment of the Doppler frequency shift in general relativity framework including non-dispersive refractivity effects can be found in Synge (1960). Then, using the flat Minkowski metric with the signature (+, −, −, −), we have and u⁰ = Γ, where Γ is the Lorentz factor of special relativity which is defined as . In addition, it turns out that k₀ is constant along null geodesic for that space-time metric, then Eq. (7) becomes (9)

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 (10)

where S is the distance from the transmission till the closest approach, K is the impact parameter of the light beam trajectory, and ɛ is the refractive bending experienced by the ray between the transmission and the reception. If ɛ can be derived from the cross product between and such as , 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 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.

2.3 Refractive profile dependency

Except in Sect. 4 where we remain as general as possible, in this work we will consider that the index of refraction coincides with the surface of constant Φ, where Φ is a generalized potential which will be defined later in Sect. 5.

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 & 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 & 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 Φ = U − Φ_C, where Φ_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 & Wolf 1999), it becomes obvious that the surface of constant n also coincide with the surfaces of constant Φ. Such a dependence between n and Φ 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 (11)

which possesses the same dimensions as the inverse of a potential (m^-2 s² in SI).

In the context of GB observations, α 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 ARO experiments it is seen as an output parameter which must be determined following an accurate modeling of the observable. In general, α 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).

2.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 Sect. 3, 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 α (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 Sect. 4, we derive a comprehensive mathematical framework which allows to study generic refractive profiles. It is obtained using the reference solution derived in Sect. 3 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 Sect. 5, 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 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 spotlight 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 the presence of horizontal gradients inside a spherically symmetric atmosphere for a specific application to GB observations. Physically, these gradients could be generated for instance by horizontal changes in the temperature.

Finally in Sect. 6, 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).

3 Reference solution

The purpose of this section is to build a reference solution, for an idealized problem. As we will see in the next section, 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 α which imposes a particular type of spherical symmetry. Each are discussed in turn.

3.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, that is n = n(h) with h = |h|. Therefore, it is shown from Eqs. (1) and (4) (see also Born & Wolf 1999) that the quantity which is defined as (12)

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 (13)

in which ϕ is the angle between the separation vector h and the unit tangentvector to the ray s (see Fig. 1). K = |K| is a constant value called the impact parameter, hence, in the following 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 θ which are defined such that the components of the separation vector are given by (14)

We have imposed that the plane coincides with the propagation plane. The -direction is collinear to K, the -axis is directed along the line of nodes (defined in Sect. 3.3), and the -axis is orthogonal to , and . We take the triad of vector to form a right-handed vectorial basis that we will refer to as the propagation basis since 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 . The third member is , which is necessarily orthogonal to and collinear to . We notice that in the specific case where the index of refraction is a function of the height only, is a constant vector. The second member of the basis is which is defined to be orthogonal to and . We take the triad of vectors to form a right-handed vectorial basis. In the propagation plane, we have the relations (15)

3.2 Monopole approximation and constancy of α

In Sect. 2.3, we have discussed that assuming hydrostatic equilibrium together with elementary theory of dispersion, one can infer that n = n(Φ) where Φ 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 (16)

in which μ = 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 (Φ = U₀), we deduce that the refractivity profile becomes spherically symmetric (see Sect. 3.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 α is constant. Therefore, the expression of the index of refraction is explicitly inferred from Eqs. (11) and (16) assuming α constant (17)

In this expression, η is a constant representing the value of the index of refraction at infinity. For instance, if 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 η = n_H − αμ∕H. The quantity which is defined as N(h) ≡ n(h) − η 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))

The light path, which is described by h, is an hyperbola. The parameter κ is a constant which is a simple function of K (cf. Eq. (A.8)) 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. (A.10a) and (A.21)) 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. (A.19)). In Eqs. (18), f is the true anomaly and is defined in Eq. (A.10b) by means of an other constant of integration, ω, known as the argument of the closest approach. The constancy of ω is linked to the conservation of the eccentricity vector defined in Eq. (A.22).

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. (A.30)). 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. (A.31)). 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, ω, 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.

3.3 Hyperbolic path in space

The solution which is described in Sect. 3.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, θ) 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 . We adopt the plane as a reference plane and the -direction as a reference direction. As before, the plane is the propagation plane while the -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 Ω, is called the longitude of the node and is the angle between the -direction and the intersection between the propagation and the reference planes. The second one, labeled ι, is called the inclination and represents the angle between the -direction orthogonal to the propagation plane and the -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 , , and into the reference frame. The radial, tangent, and normal unit vectors are given by

Regarding the direction of the closest approach, we have (20)

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 relations

We have used the previous expressions of , , , and . 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 elements (p, e, ι, Ω, ω, S, T) instead of the h, s, and t. If we do not specify that s is a unit vector (as 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 developed later on in Sect. 3.4.

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 & Clemence 1961; Murray & Dermott 2000; Poisson & 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 Sect. 4.

3.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 Sect. 4.2. Here, we are about to further simplify the previous expressions derived in Appendix A by using , 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. (A.19), we deduceE = η²∕2, and thus, from Eq. (A.21), one infers (22)

which shows that the eccentricity reduces to the inverse of the refractivity at K. This expression allows to express κ (cf. Eq. (A.8)) as a unique function of the eccentricity (23)

From these two last equations, it is straightforward to show that the semi-latus rectum (cf. Eq. (A.10a)) is expressed as (24)

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

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

Insertion of Eqs. (22) and (26) into the two Kepler-like equations (cf. Eqs. (A.30) and (A.31)), let one to deduce (27)

as well as (28)

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

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 ē (cf. Eq. (A.37)), we obtain ē = 0, and we deduce the evolution of the argument of the bending angle (30)

We can now examine the vacuum limiting case. As already mentioned, the parameter α 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 α → 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. (A.9), and by taking the limit when α → 0. In this case, it is seen from Eqs. (22) and (26) that (31)

also, we endup with (32)

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

which shows that ψ remains constant whatever is the value of the hyperbolic anomaly, meaning that the bending angle is null and thus that the trajectory is a straight line.

3.5 Applications to realistic cases

At this level, the formulas presented in Sect. 3.4 can already be used either to model atmospheric delays or to process ARO 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. Indeed, expressions in Eqs. (25), (28), and (30) allow us to compute explicit expressions for the observables (cf. Eqs. (6), (9), and (10)) using directly the reference solution.

The expression for 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 α → 0 and by using Eq. (A.27a). We end up with (34)

where F₁ = F(f₁) and F₂ = F(f₂).

The expression for the Doppler frequency shift (cf. Eq. (9)) is also immediate since it is given by the scalar product between Eq. (25b) and , where 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 (35)

where and , likewise for . 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 ɛ ≡ ψ(f₂) − ψ(f₁), 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, such as the eccentricity, are function of the impact parameter. Such an expression may be tedious to derive, thus a numerical evaluation of the derivative d ɛ∕dK 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 Sects. 3.4 and 3.5 have been carried out assuming that the index of refraction is exactly given by Eq. (17). Furthermore, we have seen in Sect. 3.2 that this refractive profile corresponds to a very peculiar type of spherical symmetry (α = constant). So, in order to extend the reference solution to any kind of spherical symmetries (α≠ constant), 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 α_k with k = 1, …, m. In addition, inside each layer the index of refraction follows Eq. (17) where the value of η_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 α_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 & Murray (1973) or Mendes & 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 ARO experiments is to deduce the value of α_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 which relies on α_k. Therefore, the difference between the observed and the computed Doppler shift can be minimized by varying the value of α_k. At the end, each layer provides its own determination of α_k. Therefore, n can be retrieved overall the vertical profile by integrating Eq. (11).

4 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 section, we address the problem of finding a new solution to a more complex situation where the index of refraction is as general as possible.

4.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 (36)

with n₀ being the index of refraction given in Eq. (17) and δn being a new non-constant contribution¹. 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. (A.3) (37)

Here, g still refers to the magnitude of the local spherical acceleration (cf. Eq. (A.4)), 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, f_pert, represents asupplementary spatial variation of the index of refraction and is given by f_pert = ∇δn. We do not makeany 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 (38)

Here it must be pointed out that the basis has a different meaning with respect to the one introduced in Sect. 3. Indeed, as discussed above, if f_pert = 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 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 d K∕ds = h ×f_pert, which can be expressed in terms of the perturbing gradient’s components (39)

Then, a comparison between Eq. (39) and the derivative of with respect to s reveals that the magnitude of K changes according to (40)

From these two last equations, we predict that a normal perturbation will only affect the direction of the impact parameter vector whilea 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 & 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. (A.16) and (A.22), 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, ι, Ω, ω, S, T) are determined from K, E, and e (see Appendix. A), they are also not constant anymore. The determination of their evolution is the subject of the next two sections.

4.2 Method of variation of arbitrary constants

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

where h₀, , and t₀ are the solutions respectively for the separation vector, the direction, and the light time along the ray for an hyperbolic path, that is to say for n₀. The C-vector represents the hyperbolic elements which are constants for the unperturbed trajectory. These solutions have been fully described in Sect. 3.

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 & 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 h₀ and (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 C(s). The general solution must be expressed as (42)

A subtle consequence of this choice, is that s and s₀ cannot be two unit vectors at the same time. Indeed, from the definition of in Eq. (2) and from Eq. (42), we notice that . Moreover, considering that the eikonal Eq. (3) must hold true for the index of refraction under study (n), we infer that s is a unit vector, thus s₀ is not. Consequently, from now, we are not allowed to use the hyperbolic expressions derived in Sects. 3.4 and 3.5, which are based on the assumption that s₀ 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 (cf. Eq. (18b)), allows us to fix the expression of the Doppler frequency shift once and for all. Indeed, Eq. (9) depends only of , 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 Eq. (42), we can determine an explicit expression for the variation of the dimensionless hyperbolic energy E. Differentiating Eq. (A.16) and making use of Eqs. (38) and (18), we get the following expression (43)

in which (44)

is a new perturbing component. The term δn has been introduced in Eq. (36) and refers to the difference between the actual and the hyperbolic index of refraction. We recall that δn must contain only the variable part of the difference as discussed¹. 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, 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 d E∕ds, d K∕ds, and d e∕ds.

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

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

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

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

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

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 & Will 2014). Instead, we are to simplify the algebra by using a geometrical relation (Burns 1976). The computation starts by differentiating Eq. (A.9) and by making use of Eqs. (A.10), (A.8), (47), and (46), which leads to (50)

On the right-hand side we have substituted the following relationship (51)

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 Ω will decrease θ by an amount equal to the projection of the magnitude of the tilt on the propagation plane. So, the geometrical contribution to θ changes is exactly given by the term − δΩ cos ι. This point is illustrated in Fig. 3 and is also discussed in Burns (1976). With this geometrical relation, the expression of the change in ω 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. (A.10b) and by substituting it into Eq. (51) (52)

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 -direction and measured inside the propagation plane (see Fig. 3). The same is also true in celestial mechanics as mentioned by Poisson & Will (2014). Substituting Eqs. (50) into (52) provides the following expression (53)

The right-hand side reduces to when the perturbation is null (i.e. when dE∕ds and dK∕ds vanish).

4.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 dE∕ds, d K∕ds, and d K∕ds, it is possible to infer the following relationships

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 ι and Ω 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 δn (through ), which represents a non-constant change in the expression of the index of refraction. Finally, we see that ω is the only onebeing 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 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 for instance Eq. (5). This can also be inferred once Eq. (4) is written as d s∕ds = n⁻¹[∇n −s(s ⋅∇)n], which represents the rate of change of the curvature of the ray (Born & Wolf 1999).

As in celestial mechanics, when the inclination is null, ω and Ω 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 & Clemence 1961; Murray & Dermott 2000; Morbidelli 2002; Poisson & Will 2014) by introducing new angles, such the longitude of the periapsis which is defined as ϖ = ω + Ω. In the next, we will continue with ω and Ω, 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 respectto 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 α → 0 (which implies e →∞), 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.

4.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. αg ≫|f_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

In the two last equations, we have introduced the non-hyperbolic contributions to the geometrical length and the light time, such that δs = s − s₀ and δt = t − t₀, where s₀ and t₀ are the hyperbolic contributions which are given explicitly in Eqs. (A.30) and (A.31), respectively.

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 Sect. 3.5). To do so, one can differentiate Eqs. (13) and (A.32), and then makes use of Eqs. (A.33) to deduce

which ismore general than Eq. (A.34) since it includes, besides the hyperbolic term, the variations of K and ω. Substitutingthe differentials into that equation and keeping only the first order in the components of the perturbation, we end up with (55h)

where we have removed the hyperbolic contribution, ψ₀, which is given exactly in Eq. (A.36). We have defined δψ = ψ − ψ₀ 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. (B.2) gives us directly the first order effect on the Doppler frequency shift.

4.5 Method of integration

Equations (55) are general enough to be applied to different geometries and different kind of problems as ARO 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 C the vector of the following elements, C = (p, e, ι, Ω, ω, δs, δt, δψ), we have (56)

where f₁ and f₂ are respectively the true anomaly at the transmission and at the reception point, with f₂ ≥ f₁. The term dC∕df′ 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 ARO 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 α ≪ 1 and then, from Eq. (22), we deduce e ≫ 1 (as long as K ≫ αμ). The meaning is that the shape of the hyperbola departs slightly from a straight-line trajectory. However, if this approximation is largely relevant for ARO, 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 αμ, which implies e ~ 1. Hence, the solutions at first order in 1∕e are not valid for GB experiments at very high elevations (i.e. when K ~ αμ). For instance in the Earth’s atmosphere, if we consider that the expansion in 1∕e does not hold for e ≲ 10, we must consider only trajectories with K ≳ 30 km, which corresponds to elevations angles ≲ 89.7°. Indeed, following Aparicio & Laroche (2011), for an average parcel of air at sea level, the refractivity is approximately (3 ± 1) × 10⁻⁴, which represents a value of α = (5 ± 2) × 10⁻⁶ km^-2 s². So, the region around the zenith-direction, in which the eccentricity is e ≲ 10, is enclosed inside a cone with a half top angle which remains ≲0.3°. All the light path trajectories with K ≲ 30 km, are enclosed inside that region.

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₂ − f₁) 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₂ − f₁) (57)

Therefore, for all the following applications (see Sect. 5), 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.

5 Applications

In this section, 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 Sect. 2.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 thesection by studying the effects due to the presence of horizontal gradients in a spherically shaped atmosphere.

5.1 Definitions and hypothesis

In this section, we introduce the basic ideas which will help us to formulate the internal 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 as the propagation frame , the polar basis , and the reference frame . 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 & Piwnicki (1999) and Rozanov & Sochilin (2005), the equation of geometrical optics expressed in the reference frame is (58a)

where n still refers to the value of the index of refraction as measured by an observer comoving with the refractive medium. The additional term refers to a dragging contribution which is given by (58b)

where v is the velocity vector of the medium, with v = |v|. γ is the Fresnel’s drag coefficient introduced in Eq. (A.20). It might be seen from this expression, that the corrections to Eq. (4) are on 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 × 10⁴ km h^-1 with zonal winds asymmetric velocity ranging between ± 550 km h^-1 (Kaspi et al. 2018), also v∕c ≲ 10⁻⁴. Therefore, we will consider the effect of the moving medium as a perturbation to the hyperbolic path, (59)

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 centered at the planet’s center of mass and rotating with it (we call P the central planet). The -axis is chosen to be orthogonal to the equatorial plane of P. The angular velocity vector is considered to be constant and given by . Since w is independent of time, the -axis is spatially fixed. Therefore, we choose the -axis of the reference frame to be collinear with the -axis. The frame 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 & 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 potential³ (60)

once written in the fluid rotating frame. In the following, we consider that the dominant term of that expression is the monopole term U₀ of the gravitational potential of the fluid planet. The other terms will be considered as perturbations before U₀. The U_l terms are the non-spherical parts of the central planet self gravitational potential, Φ_C is the centrifugal potential due to planet’s proper rotation, and U_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 (61)

where Φ is given by Eq. (60), and α_(k) ≡−d^kn∕dΦ^k for k ≥ 2. We keep the same notation as in previous sections, meaning that α ≡ α₍₁₎ 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 α_k coefficients.

Equations (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).

5.2 Non-linearity

Previously in Sect. 3, we have assumed dn∕dΦ = constant (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 & Pavlis (2007), a constant value of α 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, Φ = U₀. All the mathematical details are exposed in Appendix B.1 and the change in hyperbolic elements for any degree k is given at leading order in 1∕e in Eq. (B.3). The application to the quadratic contribution (k = 2) is given in Eq. (B.4).

After integration, it is shown that whatever function of the Newtonian potential the index of refraction is, p, ι, and Ω 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₂ ≥ f₁ and k = 2

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. (B.2), it can be seen that the later provides directly the change in the Doppler measurement due to α₍₂₎.

5.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 Φ = U₀ − Φ_C, where Φ_C represents the centrifugal potential. In addition, the dragging effect due to the rotational motion of the medium is described by the additional contributionf_drag in Eqs. (58). The solutions describing the change in hyperbolic elements due to a steady rotating atmosphere are given in Eqs. (B.8) and (B.11). All the mathematical details are presented in Appendix B.2.

On one hand, Eqs. (B.8) 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

with f₂ ≥ f₁. On the other hand, Eqs. (B.11) 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. We copy the expressions of the non-hyperbolic contribution to the light time and the refractive bending for f₂ ≥ f₁

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. (B.11g) that the non-hyperbolic contribution to the light time reduces to the non-hyperbolic contributionto the geometrical length divided by the speed of light in vacuum. This comes from the fact that the perturbing component is null for the light dragging effect while it is not for the centrifugal potential contribution.

5.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 contributionto the total self-gravitational potential of the planet.

The non-spherical components of the gravitational potential of the central planet can be modeled using the simplified following generalized potential Φ = U₀ + ∑_lU_l where U_l ∝ R^l∕h^l+1 with R the equatorial radius of the central planet and l ≥ 2. The solutions for the change in the hyperbolic elements are derived in Appendix B.3 and the results for l = 2 are given in Eqs. (B.21). We copy here the expressions of the non-hyperbolic contribution to the light time and the refractive bending for f₂ ≥ f₁

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₂ 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. (B.21) and (B.8) is always proportional to

where C represents the set of hyperbolic elements. The subscript J₂ 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₂ 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.

5.5 Tidal potential

In systems containing close orbiting bodies, such as 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 Φ = U₀ + U_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 Appendix B.4, and the evolution of the hyperbolic elements are given in Eqs. (B.35). We remind that these equations are derived considering a tide raising body (with a standard gravitational parameter μ_A) moving on an equatorial circular orbit around the central planet, and a photon path lying inside the equatorial plane (i.e. ι = 0). We copy here the expressions of the non-hyperbolic contribution to the light time and the refractive bending for f₂ ≥ f₁

Comparison between the magnitude of the change in hyperbolic elements in Eqs. (B.35) and (B.21) reveals

For the tides raised by the Moon on Earth, the ratio is of the order of 10⁻⁶. For tides raised by Titan on Saturn, the value is 10⁻⁷, and for tides raised by Io on Jupiter the ratio is of the order of 10⁻⁵. For all those cases, the tides effects are at least five 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 thecase of a light ray crossing Titan’s atmosphere; substitution of numerical values (we took the Titan’s J₂ 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. (B.35) 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 & Dermott 2000). In Appendix 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. (B.38), we see that the non-elastic contribution is a least wτ times the elastic one, with w the magnitude of the angular velocity of the central planet, and τ the time-delay which is related to dissipative phenomenon.

For the Earth, a typical value of τ corresponding to rotational semi-diurnal deformation is 4 min (Folkner et al. 2014), so the ratio is of the order 2 × 10⁻². 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 where Q is the specific dissipation function (Murray & Dermott 2000). The numerical substitution is done by taking value of k₂ ∕Q (with k₂ the gravitational Love number of degree two) from Lainey et al. (2009) together with the value of k₂ 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 × 10⁻⁶ and 7 × 10⁻⁵, which is well beyond the elastic response.

5.6 Horizontal gradients

The last study is a direct application for GB observations operating for the realization of IERS reference systems such as SLR, VLBI, or the Global Positioning System (GPS). In Chen & Herring (1997) or in Hulley & 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 subcentimeter 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 Appendix B.5, and the equations describing the change in the hyperbolic elements following horizontal gradients are given in Eqs. (B.46). 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 isat the meridian of the receiving site. In such a case, the propagation plane is aligned with the meridian and intersects Earth’s equator for ι = π∕2 and λ = Ω. For this type of inclination, it is seen from Eqs. (B.46) 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 (62)

where n_NS represents the horizontal variation of the refractivity along the north-south direction.

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. (A.10b) and (A.11). 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. Thus, the hyperbolic relations in Appendix A can be further simplified taking the limit α → 0 (see relationships in Sect. 3.4), which corresponds to light propagation in vacuum.

As an example, we give the transformation rule which allows to pass from true anomaly to colatitude. Firstly, from the coordinates (R, f₂) of the receiving station (on the Earth’s surface) together with Eq. (32), we can insert the relation K = Rcosf₂ into Eq. (62). Then, the transformation to colatitudes is given by the following relationships (we consider the photon path 1 →2′ which is depicted on the left panel of Fig. 4) (63a)

Here, for a matter of convenience for the next, we have introduced the colatitude which is referred from ’s direction instead of north’s (see Fig. 4). The difference ω′− ω appearing on the right-hand side of Eq. (63a), is function of the location of the receiving antenna () and is also function of the ratio between the magnitude of the separation vectors of the source and the receiver (q = h₁∕R). It is given by the following expression (63b)

Finally, insertion of Eqs. (63) into (62) allows one to infer the atmospheric time-delay along the photon path and for the following boundary conditions .

In order to emphasize the behavior of , 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′, and 2″) 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 constant 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 → 2. Secondly, the geometrical length inside the atmosphere is the most important for the path 1 → 2. Finally, the horizontal contributionin (φ − φ₂)n_NS into the expression of the index of refraction (see Eq. (B.39)) is the most important for φ =φ₁ and for the path 1 → 2. Thus, for this particular path, the phase velocity is minimum and consequently the time-delay is maximum.

6 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 ARO experiments or astro-geodetic GB observations. The differences remain in the role of α, 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 ARO 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 section, we test the validity of the derived solutions in the context of ARO experiments. We focus on Eqs. (B.21) 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 Sect. 1, this problem has never been solved analytically in a complete and satisfactory way.

In order to assess the validity of Eqs. (B.21), we have performed a numerical integration of Eqs. (1), (4), and (5) across the refractive profile given in Eqs. (B.15) and (B.17) for l = 2. The numerical integration has been carried out in double precision, with a numerical error tolerance of 10⁻¹¹, for different values of α and J₂. (Actually, instead of α, we work in the following with the dimensionless parameter αμ∕R, which represents the hyperbolic refractivity evaluated at the level of the radius of the central planet, N₀ (R) = n₀(R) − η.) The tested values of the refractivity and J₂, range from 10⁻¹ to 10⁻⁵, and 10⁻¹ to 10⁻⁶ restrictively. 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. (B.21). Let us emphasize that δs, δt, and δψ 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 ψ, which is easily inferred from results of the numerical integration. From an analytical point of view, the total change in s, t, and ψ, is simply given by the sum of hyperbolic (see Eqs. (A.30), (A.31), and (A.36)) and non-hyperbolic contributions (see Eqs. (B.21f)–(B.21h)). In Fig. 5, we show the evolution of the relative error on the total change in the hyperbolic elements and s, t, and ψ.

Firstly, for values of the refractivity N₀(R) ≲ 10⁻⁵ and for values of J₂ ≲ 10⁻⁴, it is seen that the solutions seems to suffer from a lack of accuracy (cf. errors on p, e, ι, Ω, and ω in Fig. 5). However, this loss of accuracy is only due to numerical noise, since for small values of J₂, 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₀ (R) = 10⁻⁵, we obtain a value of the semi-latus rectum at p(f₁) ~ 10⁸; moreover, the analytical solution predicts Δp_ana ~ 10⁻⁵ for J₂ = 10⁻⁶, which represents a change in Δp_ana∕p(f₁) ~ 10⁻¹³, that is to say 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 ofmagnitude of the refractivity we also change the order of magnitude of the relative error on the total change in the value of all the elements. 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. Indeed, by changing the order of magnitude of the value of the refractivity, we also change the value of α which in turn changes the order of magnitude of the eccentricity. From Eq. (22), we notice that increasing α 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⁻¹ 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 ~10%. For typical values of the refractivity around 10⁻⁴, the analytical solutions in Eqs. (B.21a)–(B.21e) are found to be accurate at the level of ~ 0.01% which is a good agreement considering their simplicity.

For the geometrical length, the total light time, and the refractive bending, the relative error evolves linearly with J₂ in the log–log plot (constant gap between each curve representing a change in the value of J₂). This behavior is due to the fact that the hyperbolic contribution is included within the analytical computation of s, t, and ψ. Indeed, for the hyperbolic elements, as discussed above, when J₂ tends to be null, the computation of ΔC_num generates a lot of numerical noise since the hyperbolic elements are constants when the perturbation vanishes. In the case of s, t, and ψ, when J₂ goes to zero, the computation of ΔC_num does not generates any numerical noise since a non-null hyperbolic contribution remains besides the vanishing non-hyperbolic one. Then, it is seen that the relative error declines to finally reach the numerical noise level showing that the hyperbolic path is an exact solution when the perturbation vanishes. In order to prove that the non-hyperbolic evolution of s, t, and ψ is well described by Eqs. (B.21f)–(B.21h), one can compute the relative error considering only the hyperbolic contribution in ΔC_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 around 10⁻¹, the maximum relative error is reached for the maximum tested value of the oblateness that is J₂ = 10⁻¹. For instance, considering the hyperbolic and the non-hyperbolic contributions, we are able to provide analytical solutions which are accurate at the level of ~0.01% for both s and t, and at the level of ~ 1% for ψ. For a value of the refractivity around 10⁻⁴, and for a value of Jupiter or Saturn’s J₂ at ~ 10⁻², we might expect errors at the level of one part in 10⁸ for both s and t, and at the level of one part in 10⁵ for ψ. For the same values of the refractivity and J₂, if we only consider the hyperbolic contribution, the relative error grows to one part in 10⁷ for both the geometrical length and the light time, and to 0.1% for the argument of the refractive bending.

7 Conclusion and perspectives

In this paper we have demonstrated, in Sect. 3, 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 ARO experiments. In Sect. 3.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 Sect. 4 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 Sect. 4.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 Sect. 5 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 Sect. 6, 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₂ 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 around 10⁻⁴, the relative error on the total change in the hyperbolic elements is of the order of 0.01% and is independent of the value of J₂. For the other elements like the geometrical length, the time-delay and the refractive bending, the relative error evolves with J₂ as shown in Fig. 5. For Jupiter or Saturn’s typical value of J₂ at ~ 10⁻², the relative error is of the order of one part in 10⁸ for both the length and the time, and one part in 10⁵ 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. (A.31)) and the non-hyperbolic contribution (cf. Eq. (B.21g)).

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 & Sullivan 1981; Phipps & Withers 2017; Phipps et al. 2018). This is also something which can can be easily achieved from Eqs. (54) by taking the limit α → 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 & Herring 1997 for VLBI, Hulley & 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 tacklethese issues by providing for instance the contribution to the time-delay due to the atmosphere’s oblateness and the one due to horizontal gradients, as discussed in Sects. 5.4 and 5.6.

In the context of ARO 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⁸ and 10⁵ 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₂ = 10⁻² and N₀(R) = 10⁻⁴, 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 Sect. 5, could be developed in order to process real ARO 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 & 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. The authors would like to thank the anonymous referee for valuable suggestions and comments for improving the manuscript. A.B. is thankfull 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 Sect. 3, 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 section, a complete solution to Eqs. (1), (4), and (5). We start by determining the shape of the light trajectory, then we focus on the light time. We provide in Table A.1 a reminder about the different notations introduced in the main text and used thereafter. A.1:

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 are decomposed as

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

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 (A.3)

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 (A.4)

It is obvious from the presence of α in Eq. (A.3), that the light ray experiences the gravitational pull via the media in which the ray propagates through.

From Eqs. (A.2) and (A.3), it is seen that the absence of a transverse component (along ) in the gradient of the index of refraction implies that is a conserved quantity. From Eqs. (12) and (A.1) we immediately deduce (A.5)

Inserting this expression into the radial part of Eq. (A.2), we find (A.6)

We change the independent variable from s to θ adopting the convention that a prime denotes the differentiation with respect to θ. We also adopt u ≡ 1∕h as a convenient substitute for h and derive adifferential equation for it. Making all these substitutions into Eq. (A.6), we quickly arrive to (A.7)

with (A.8)

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

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

They are respectively known as the semi-latus rectum and the true anomaly. The last parameter, ω, is also a constant of integration and its role is better understood when h is minimal, that is to say when θ = ω. It is seen that ω 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 ω is completely determined by the initial conditions of the problem. For instance, if we call (h₁, θ₁) the components of the position vector of the transmission point, then from the solution in Eq. (A.9), we end up with (A.11)

where the sign is taken to be positive when − π∕2 ≤ θ₁ < π∕2 and negative elsewhere. The minimal separation is also determined from Eq. (A.9) evaluated at the closest approach, such that (A.12)

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 ≥ 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 (A.13)

which is a negative quantity.

Non-periodic solutions, which are hyperbolic trajectories, also imply the introduction of a new important angle, Δf_∞, describing the net change between the two asymptotes’ directions. Based on Eq. (A.9), it can be seen that h goes to infinity when cosκf tends to − 1∕e, which gives rise to two asymptotic solutions for κf. Let call the negative solution and the positive one. If we introduce , we immediately deduce (A.14)

Since, we focus on hyperbolic trajectories (e > 1), we will always have Δf_∞ < 2π. In ARO experiments, the geometry is such that Δf_∞ can be geometrically related to the refractive bending (indeed from Fig. 1, we deduce Δf_∞ = π + ɛ) which is itself directly in relation with the changes in the Doppler measurements, at first order (see Eq. (B.2) or Fjeldbo & Eshleman 1965, 1968). Thus, Eq. (A.14) 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. (A.1b)). We thus need to determine the expressions of the radial and the transverse components; the later depends on the angular velocity. Invoking Eqs. (A.9) and (A.5) and differentiating them making use of Eqs. (A.10), reveals that

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

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 for instance that the semi-latus rectum is linked to the conservation of the magnitude of the impact parameter (cf. Eq. (A.10a)). 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. (A.6) by n² ḣ. Then, one can see that each term is a total differential with respect to s. In addition, one can recognize the squared components of Eq. (A.1b) in the integrated equation, and we end up with (A.16)

where E is a constant of integration. That constant would be similar to the energy in classical mechanics. Indeed, by computing the square of from Eqs. (A.15) and (A.1b), we end up with (A.17)

where we have used Eq. (A.13). Then, insertion of Eq. (A.17) into (A.16), leads to (A.18)

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. (A.16) can be put under a more fundamental form. Indeed, by making use of the definition of n (cf. Eq. (17)), we deduce (A.19)

where (A.20)

is known as the Fresnel’s coefficient (η has been previously defined in Eq. (17)). This expression of E is more fundamentalthan Eq. (A.16) 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. (A.18), it is seen that the semi-major axis is linked to the conservation of the dimensionless energy parameter, E. At the sametime, the constancy of the eccentricity is assured, because of the definition (Eq. (A.13)). So, in that sense the eccentricity is linked to both the conservation of energy and the conservation of the magnitude of the impact parameter. Equations (A.18), (A.13), and (A.10a) can be used to infer the eccentricity (A.21)

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 κ is different from unity, the eccentricity vector does not show up as easily as in celestial mechanics. It is given by the following relationship (A.22)

where we have introduced the vector A which possesses the following components in the polar basis (A.23)

The vector A is expressed in terms of fundamental quantities after determining f from Eq. (A.9), and expressing and in terms of h and K. The eccentricity in Eq. (A.22), must be expressed in terms of E and K thanks to Eq. (A.21).

The constancy of e can be demonstrated by substituting Eqs. (A.15) into (A.1b), then inserting the result together with Eq. (A.5) into (A.22), which finally gives with Eq. (15) (A.24)

Since e and ω are two constants of integration, e is conserved during the propagation of light. Moreover, it is seen that the constancy of ω is linked to the conservation of the eccentricity vector. It always points toward 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, ω. Those three parameters can totally be expressed (cf. Eqs. (A.10a), (A.21), and (A.24)) 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.

Equation (A.15b) 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. (A.15b) as (A.25)

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. (A.25) into (A.26)

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

Equations (A.25) and (A.26) 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

and is summarized with the half-angles formula (A.27c)

We emphasis that determining the hyperbolic from the true anomaly using (A.27c) could lead to numerical errors. Indeed, the inverse function of the hyperbolic tangent is not well defined when the argument is close to ± 1, which generates numerical noise. To avoid this issue, the inverse function of the hyperbolic sine (cf. Eq. (A.27a)) can be preferred.

Equations (A.27) allow us to introduce the following relations which relate the differentials of f and F (A.28)

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

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

in which (A.30b)

is the geometrical length mean anomaly. Similarly, expressing Eq. (A.29c) in term of the light time instead of length path with the help of Eq. (5), let one to directly integrate Eq. (A.26), such that (A.31a)

in which (A.31b)

is the light time mean anomaly.

From these two Kepler-like equations, it is seen that the light time (cf. Eqs. (A.31)) is not just the geometrical length (cf. Eqs. (A.30)) 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, ψ. As seen in Fig. 1, ψ is the angle between the orthogonal direction to the line of nodes and the tangent to the light ray. In the context of ARO experiments, this angle plays an important role. Indeed, from Fig. 1, it might be seen that the net change in ψ between the transmitter (labeled 1) and the receiver (labeled 2) is exactly the refractive bending, ɛ ≡ ψ₂ − ψ₁. Moreover, as shown in Eq. (B.2) 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 ARO experiments.

The differential relationship between ψ and f can be inferred from Fig. 1 which allows us to deduce the following relationship (A.32)

We remind that the angle ϕ has been previously defined in Eq. (13) as the angle between and s. Hence, from Eq. (A.1b) we immediately deduce (A.33)

Then, differentiating Eqs. (13) and (A.32), and making use of Eqs. (A.33), we obtain the following relation (A.34)

which can be written as (A.35)

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

in which we have introduced (A.37)

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

Appendix B Mathematical details for applications

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

B.1 Non-linearity with U₀

We first simplify the problem by assuming that the generalized potential reduces to the spherical part of the Newtonian potential (Φ = U₀). Consequently, Eq. (61) and the perturbing gradient now read

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

Equations (B.2) can then be inserted into Eqs. (55) and the integrand is expanded in power of 1∕e (see Sect. 4.5). The change in the hyperbolic elements is given for any k (with k ≥ 2) by the following expressions

where we have kept the leading order in 1∕e, and where f₂ ≥ f₁.

Considering for instance 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

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. (B.4) are valid for all f ≥f₁.

B.2 Steady rotating atmosphere

The term Φ_C in Eq. (60) is the centrifugal potential due to proper rotation of the fluid body and is given by (Kaula 1966; Poisson & Will 2014) (B.5)

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

Once and are projected into the polar basis, we obtain the following components

Inserting Eqs. (B.7) 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

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. (B.8) are valid for any f ≥ f₁.

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 (B.9)

in which, we have assumed a steady rotating atmosphere such that v = w ×h. If we express the dragging contribution in the polar basis, we obtain the following non-null components

where N is the refractivity defined as N = n − η. We notice the absence of the -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

Once again these equations are defined within a constant which is determined by the transmitter’s position. Eqs. (B.11) are defined for f ≥ f₁.

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 & Dermott 2000) or the symmetric and trace free (STF) tensors formalism (Hartmann et al. 1994; Poisson & Will 2014). Using for instance the later one, it can be shown that the expansion of the non-spherical part of the gravitational potential is given by (B.12)

for l ≥ 2. We recall that G is the gravitational constant, I^⟨L⟩ is the multipole moment STF tensors of the mass distribution, and ∂_⟨L⟩ is the STF tensor made with the partial derivatives with respect to h. We refer to Hartmann et al. (1994) and Poisson & Will (2014) for notations and definitions about the use of the STF tensors.

We remind that the unit vector has been chosen in Sect. 5.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 -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 & Will 2014). This symmetry imposes the multipole moments STF tensors to be proportional to Ẑ^⟨L⟩ which is the STF tensor formed with the unit vector. Hence, we have the relation (B.13)

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 & Will 2014), such

where and P_l (χ) are the well-known Legendre polynomials of degree l, we can rewrite the perturbing gravitational potential in Eq. (B.12) in terms of the multipole moments J_l. Keeping the linear part in Φ in the development of the index of refraction (cf. Eq. (61)) with Φ = U₀ + ∑_lU_l, we determine an expression for n in term of the non-spherical part of the gravitational potential (B.15)

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 jth component of the gradient of the index of refraction which is given by (B.16)

Differentiating Eq. (B.12) and using the previous properties about STF tensors (cf. Eq. (B.14c)), we deduce (B.17)

In the polar basis the components of f_pert are

Usually, for modestly deformed bodies (e.g. the planets in the solar system), the most important term in the developments (B.15) and (B.17) is the one proportional to the quadrupolar moment J₂. 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

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

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

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 ≥ f₁.

B.4 Tidal potential

In the fluid rotating frame, the perturbing tidal potential can be expressed as (Poisson & Will 2014) (B.22)

where are the tidalmoments being time dependent STF tensors. They are defined by (B.23)

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

Following Poisson & Will (2014), U_¬P is given by (B.24)

in which the integration must be performed in the vicinity of the external bodies A where ρ_A (t, x) is supposed to be non-null. Thus, the dummy variable of integration will range over values close to r_A (t), and we can introduce x = r_A(t) + [x −r_A(t)]. Then, we are left with two vectors in the denominator, the first one, h + r_AP(t) with r_AP(t) ≡r_P(t) −r_A(t), representing the separation between the center-of-mass of A and a point h at which the external potential must be evaluated, and the second one, x − r_A(t), representing the separation between the center-of-mass of A and a point lyingin the vicinity of A. Therefore, because h + 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 (B.25)

We recall that means the multipole moment STF tensors of body A and that the derivative must be taken with respect to h and then be evaluated at P’s center-of-mass (i.e. h = 0). This notation is shorten noticing that , where the derivative is taken at the extremity of the vector r_AP(t). From Eqs. (B.25), (B.22), and (B.23), we can define the perturbing external tidal potential (B.26)

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 jth component of the gradient of the non-constant change in the index of refraction are explicitly written as

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 r(t) ≡−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≪ M, with M the mass of the central planet. Therefore we can safely consider only the monopole contribution of A, that is l′ =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, that is 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 position vector of the satellite is given at time t by (B.28)

where the mean motion of the satellite around the central body is given by Kepler third law of motion (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 (B.29)

in which we have defined λ(t) ≡ Nt − Ω 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

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

In those expressions, the components of the unit vector directed along r are noted and can be inferred from Eq. (B.29) in the polar basis.

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

Here, we have introduced the unitless vector and the unitless parameter to simplify the notations. They are explicitly given by

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

One can now insert Eqs. (B.32) 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, ι = 0 and the equation can be further simplified. However, as discussed in Sect. 4.3, when the inclination is null the longitude of the node and the argument of closest approach are undefined. Therefore, instead of considering ω and Ω separately, we consider the longitude of the closest approach, ϖ = ω + Ω, which is measured from the -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

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 ≥ f₁.

In order to account for the non-elastic response of the fluid following a tidal stress, it is shown in Poisson & Will (2014) that the first order effects can be assessed by replacing by where τ is called the time-lag and is function of the Energy dissipation. A computation of at the instant t − τ, changes the instantaneous position of the perturbing body (cf. Eq. (B.28)) into a delayed position, r(t − τ). Most of the time, τ 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 τ. In this condition, only internal process occurring on time scale faster than the orbital motion can significantly change . 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

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 unitvector is given by

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 τ by

As before, we have introduced the trigonometric coefficients and which are given by

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

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 (B.38)

where C represents the vector of the hyperbolic elements.

B.5 Horizontal gradients

Let us consider the Earth centered frame with the -axis directed through the north pole, the plane superimposed within the Earth equator, and the -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 . A station on Earth can be located thanks to itsspherical coordinates (φ₂, λ₂, R), with R the Earth radius, φ₂ the colatitude (measured from the north pole), and λ₂ the azimuth of the station. We consider here the case of a receiving station, but the same work can totally be applied for a transmitting station. Any point along the light ray trajectory can be located through its spherical coordinates (φ, λ, 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 & Herring (1997) or in Hulley & Pavlis (2007) (B.39)

where n₀(h) is the hyperbolic index of refraction, and n_NS and n_WE are respectively given by ∂n∕∂φ|_h,λ and ∂n∕∂λ|_h,φ. In other words, n_NS represents the variation of the horizontal refraction along the north-south direction and n_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_NS and n_WE.

The perturbative gradient is given by taking the gradient of δn, with δn(φ, λ) = n(h, φ, λ) − n₀(h). Expressing the gradient in spherical coordinates we end up with (B.40)

where we have assumed that the radial variations of n_NS and n_WE are of second order. The two unit vectors and are part of the usual spherical vectorial basis which is given by the triad of vectors . It is seen that the perturbing gradient does indeed not possess any radial component if we neglect the radial variations of n_NS and n_WE.

Now letus write the transverse and the normal components of the perturbing gradient by expressing and into the polar basis . For any point, we can change from the spherical to the Cartesian coordinates thanks to two angles φ and λ

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

We have introduced the trigonometric function which is defined in Table B.1. From Eqs. (B.41) and (19) and by making use of Eqs. (B.42), we find

where the coefficients are given in Table B.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 of the index of refraction in term of hyperbolic elements. That equation can be obtained from Eqs. (B.42). Indeed at first order in f − f₂ we deduce

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

Finally, the effect on the hyperbolic elements is given by inserting Eqs. (B.45) 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 (ι =π∕2). In addition, the azimuth of the source is the same as the azimuth of the receiving station, and because the inclination is π∕2, the longitude of the node is the same as the azimuth of both the source and the station (Ω = λ). Keeping the leading order in 1∕e, and integrating over the true anomaly leads to the following expressions

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 ≥ f₁.

Appendix C Doppler shift and refractive bending

In this section, 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

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

where s_vac is the light beam direction in a vacuum and δs ≡s −s_vac is the deviation of the actual direction of the ray with respect to s_vac. In vacuum, δs = 0, and n − 1 = 0, thus ν₂ ∕ν₁ reduces to .

For application within the solar system, Eq. (B.1) can be simplified considering the small velocity approximation () together with the low refractivity limit (n − 1 ≪ 1). The later condition imposes that the deviation with respect to the light ray direction in a vacuum is small, also |δs| ≪ 1, and substitutions into Eq. (B.1) leads to the following approximation

For an hyperbolic path, we can consider that since K is constant in direction overall the light path. The angle β represents the elevation of the photon path above the straight line connecting the transmitter and the receiver (see Fig. 1).

For spherical symmetry, β₂ is proportional to β₁ by means of a trigonometric factor: β₂ = −Λβ₁. This Λ-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 ofmass (Λ →{0;1} when h₂ →{+∞;h₁}). Moreover, from Fig. 1 it is seen that the total refractive bending can be expressed as ɛ ≡ β₂ − β₁, also we end up with (C.2)

This equation shows that in first approximation, the Doppler shift due to atmospheric effects is controlled by the total refractive bending, ɛ. As discussed in Appendix A, ɛ is defined as the net change in the argument of the refractive bending, ψ, between the transmitter and the receiver (ɛ ≡ ψ₂ − ψ₁), where the exact expression of ψ has been given in Eq. (A.36).

C.2 Simplifications

For one-way Earth ARO 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₁ = n₂ = 1), also the previous equation reduces to (C.3)

when .

For one-way planetary ARO 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). In addition, the receiver (DSN antenna on Earth) can safely be considered at infinity, also Λ → 0, and the simplified Doppler shift expression becomes (C.4)

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.

References

---

All Tables

All Figures