Analytical ray-tracing in planetary atmospheres
AAnalytical ray-tracing in planetary atmospheres
A. Bourgoin, M. Zannoni, P. Tortora
Dipartimento di Ingegneria Industriale, University of Bologna, via fontanelle 40, Forl`ı, Italy ∗ (Dated: January 25, 2019) Context.
Ground-based astro-geodetic observations and atmospheric occultations, are two examples of obser-vational techniques requiring a scrutiny analysis of atmospheric refraction. In both cases, the measured changesin observables (range, Doppler shift, or signal attenuation) are geometrically related to changes in the photonpath and the light time of the received electromagnetic signal. In the context of geometrical optics, the change inthe physical properties of the signal are related to the refractive profile of the crossed medium. Therefore, havinga clear knowledge of how the refractivity governs the photon path and the light time evolution is of prime impor-tance to clearly understand observational features. Analytical studies usually focused on spherically symmetricatmospheres and only few aimed at exploring the e ff ect of the non-spherical symmetry on the observables. Aims.
In this paper, we analytically perform the integration of the photon path and the light time of raystraveling across a planetary atmosphere. We do not restrict our attention to spherically symmetric atmospheresand introduce a comprehensive mathematical framework which allows to handle any kind of analytical studiesin the context of geometrical optics.
Methods.
Assuming that the index of refraction of the medium is a linear function of the Newtonian potential,we derive an exact solution to the equations of geometrical optics. The solution’s arbitrary constants of inte-gration are parametrized by the refractive response of the medium to the gravitational potential, and by the firstintegrals of the problem. Varying the constants, we are able to reformulate the equation of geometrical opticsinto a new set of osculating equations describing the constants’ evolution for any arbitrary changes in the indexof refraction profile.
Results.
The osculating equations are identical to the equation of geometrical optics, however, they o ff er acomprehensive framework to handle analytical studies which aim at exploring non-radial dependencies in therefractive profile. To highlight the capabilities of this new formalism, we carry out five realistic applications forwhich we derive analytical solutions. The accuracy of the method of integration is assessed by comparing ourresults to a numerical integration of the equations of geometrical optics in the presence of a quadrupolar moment J . This shows that the analytical solution leads to the determination of the light time and the refractive bendingwith relative errors at the level of one part in 10 and one part in 10 , for typical values of the refractivity andthe J parameter at levels of 10 − and 10 − , respectively. I. INTRODUCTION
In the context of Maxwell’s theory of electromagnetism, the refraction phenomenon is understood as the superposition ofan incident electromagnetic wave with wavelets generated by some electric particles being accelerated in response to a localelectromagnetic field. If the medium containing the particles is su ffi ciently tenuous (like for a gas), the local field is proportionalto the incident field itself. The waves’ superposition generates a resulting signal which exhibits a phase shift with respect tothe incident one. This apparent change in the primary signal’s phase is parametrized thanks to a dimensionless parameter calledthe index of refraction of the medium. It encloses the physical properties related to the interaction between the electromagneticsignal and the material filling the medium. Therefore, any variation of the physical properties of the medium will change theindex of refraction which in turn will modify the features of the incident signal. In Astronomy, those changes in the signal’scharacteristics are the observables, thus, it is of prime importance to understand the relationship between refraction and lightpropagation.Ground-based (GB) astro-geodetic observations and atmospheric occultations (AO) are two examples of observational tech-niques requiring a careful analysis of refraction. For GB observations, refraction is an accuracy-limiting factor and it has to bemodeled in the analysis software; for AO, it is thoroughly modeled to infer the physical properties of the refractive medium. Inboth cases, the refraction usually occurs at characteristic lengths which are much larger than the wavelength of the electromag-netic signal and the problem can be studied in the approximation of geometrical optics. In that context, the di ff erent observables(range, Doppler, or attenuation) can be related to the geometry of the photon path and the light time which are both controlled bythe index of refraction. Thus, a clear knowledge of how the geometry of the light propagation is related to the index of refractionis important to understand the changes in observables.For GB observations, especially for techniques operating for the realization of the International Earth Rotation and ReferenceSystems Service (IERS), the accuracy of the analysis is largely a ff ected by errors in modeling the group delay during propagation ∗ [email protected] a r X i v : . [ phy s i c s . c l a ss - ph ] J a n of the signal through the atmosphere [1]. Usually, the atmospheric delay is evaluated at zenith and is computed with an analyticalmodel of [2] or with the more recent model (also analytic) of [3]. Then, the projection to a given elevation angle is operated withthe help of the mapping function developed by [4]. Those models are formulated under the assumption that the atmosphere isspherically symmetric and assume a radial dependency for the index of refraction. Consequently, they leave aside some possiblee ff ects due to the presence of horizontal gradients in the atmosphere , which might be somehow inaccurate. This issue has onlybeen partially overcome by [5] for Very Long Baseline Interferometry (VLBI), and by [6] for Satellite Laser Ranging (SLR),by performing numerical integration (numerical ray-tracing) across multi-layered spherical atmosphere, assuming constant re-fractivity inside each layer. Even if in those cases the index of refraction may possess non-radial dependencies, the shape of theatmosphere is still considered for being spheric. So, it is worth mentioning that beyond the radial dependency of the refractiveindex, only numerical integrations are usually used for such GB observations.In AO experiments, the basic idea is to establish radio links between a transmitter and a receiver when the latter is beingocculted by a planetary atmosphere as seen from the transmitter. Conversely to GB observations, AO experiments use the e ff ectof the refraction by the atmosphere to determine the physical properties of the medium. This is often referred to as the inverseproblem. In the literature, there exists di ff erent approaches devoted to the resolution of the inverse problem. For instance, theinitial data processing of the Mariner IV occultations was based primarily on model fitting [7]. The model used for determiningthe index of refraction profile was approximate and only valid for a planet with a thin spherically symmetric atmosphere [8, 9].Later, it as been shown [10] that determining directly the refractivity profile from the Doppler frequencies (or from the bendingangle) is a special case of Abelian integral inversion when the index of refraction is assumed to be driven by spherical symmetry.In the past, this method has been successfully applied to the dense atmosphere of Venus [11] using the data of the occultation ofMariner V. More recently, Abelian integral inversion was used to process radio data of Cassini occultation by Titan [12, 13] orGPS / MET (Global Positioning System Meteorology Experiment) occultation by the Earth atmosphere [14, 15]. Later, with theVoyager 2 occultation by Neptune in 1989 and the more recent Cassini mission orbiting Saturn, the problem of studying non-spherically symmetric atmosphere was raised. A first analytical model was proposed by [16]. This solution based on a Taylorseries expansion is suited for a numerical evaluation, but does not provide a clear understanding of the e ff ect of non-sphericalsymmetry on the photon path and the time-delay. More recently, in order to account for additional e ff ects such the light dragginge ff ect by the winds, [17] proposed a purely numerical ray-tracing technique to solve the inverse problem. This method is ofcourse 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 thegravitational pull of the planet. This parameter is either an input for GB observations or an output for AO experiments. Thegeometric 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 usefulmathematical framework which allows to perform analytical studies beyond the usual spherically symmetric case. The validityof the derived solution is assessed by comparing it with numerical ray-tracing performed on a simulated atmosphere. We try tokeep the discussion as general as possible in order to let the possibility of applying the incoming description to any situationinvolving light propagation in planetary atmospheres e.g. GB observations or AO experiments. II. DEFINITIONS AND HYPOTHESIS
In this paper, we use the geometrical optics approximation which is characterized by neglecting the wavelength of electro-magnetic waves [18]. In geometrical optics, the concept of rays is introduced as curves whose tangents at each point coincidewith the direction of propagation of the electromagnetic waves. In the following, we will mainly focus on the evolution of theseparation vector, locating a point along the light ray, and the unit vector tangent to the ray at that point. In this chapter, weintroduce the basic principles for determining the evolution of these quantities.
A. 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 s ≡ d h d s . (1) Those horizontal gradients can arise at two di ff erent levels. Firstly, they can be generated by an atmosphere with a non-spherical shape, and secondly, theycan be due to a local horizontal variation of the physical quantities entering in the computation of the index of refraction e.g. the temperature. d h is an elemental displacement vector along the ray and d s represents its length, which is defined as d s = d h · d h . Because ofthis quadratic relation, one can infer that s is a unit vector, and this reduces the number of independent constants of integrationfrom six to five. Together, h and s provide the description of the evolution of the separation vector as well as the direction of theray.For convenience, we introduce the new following quantity ∇ S ≡ n s (2)where S ( h ) is a real scalar function of the position and is usually defined as the optical path [19]. Here, the gradient is computedwith respect to h , and n ( h ) is an arbitrary scalar field called the index of refraction of the medium. Remembering that s is a unitvector, the square of the gradient of the optical path leads to the well-known eikonal equation ∇ S · ∇ S = n , (3)which is the fundamental equation of geometrical optics. The function S is also referred to as the eikonal .The height of the ray, its direction of propagation or the eikonal, changes in space according to the spatial variation of theindex of refraction of the medium in which the light ray propagates. Thus, it is possible to directly specify the ray with the indexof refraction by di ff erentiating Eq. (2) with respect to s , which leads to [19]dd s (cid:0) ∇ S (cid:1) = ∇ n . (4)This expression is given in a reference frame at rest with respect to the refractive medium. Eq. (4) can only be integrated afterhaving introduced some hypotheses about the dependencies of n . In the following, we will refer to Eq. (4) as the equation ofgeometrical optics.In order to assess a time description along the path, we introduce an additional well-known formula describing the rate ofchange of the light time with respect to the geometrical lengthd t d s = nc . (5)The introduction of this additional equation, leads the total number of independent constants of integration from five to six.Eqs. (1), (4), and (5) are the equations that we have to deal with in order to solve simultaneously the photon path and the lighttime which in turn can be used to determine the observables, as shown in the next section. B. 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 ∇ S and t are known. Therefore, this discussion highlights the great benefit of having analytical expressions describing theevolution 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-delaydue to atmospheric e ff ects τ atm = c (cid:90) s s n d s − t vac . (6)Here, the integral represents the total light time accounting for the atmospheric e ff ects and t vac is the light time in vacuum. Theintegration 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 transmittedfrequencies can be expressed as [20, 21] ν ν = ( u α k α ) ( u α k α ) = ( u k ) ( u k ) (1 + ˆ u i ˆ k i ) (1 + ˆ u i ˆ k i ) , (7)where ν and ν are the transmitted and the received frequencies respectively, u α denotes the 4-velocity vector, k α is the covariantform of the 4-wave vector which is tangent to light ray, and ˆ u i and ˆ k i are respectively given by ˆ u i = u i / u , and ˆ k i = k i / k (see[21, 22] for definitions and notations). In this expression, the subscripts 1 and 2 specify that the terms between parenthesis mustbe evaluated at the position of the transmitter and the receiver respectively.We can specify the frequency ratio in term of ∇ S by following e.g. [23]. Indeed, for a static medium, it can be shown thatˆ k i = ∂ i S . (8)Using the flat Minkowski metric with the signature ( + , − , − , − ), we have k i = − k ˆ k i and u = Γ , where Γ is the Lorentz factorof special relativity which is defined as Γ ≡ (1 − ˆ u · ˆ u ) − / . In addition, it turns out that k is constant along null geodesic for thatspace-time metric, then Eq. (7) becomes ν ν = Γ Γ (cid:16) − ˆ u · ∇ S (cid:17) (cid:16) − ˆ u · ∇ S (cid:17) . (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 [24] as A def =
10 log (cid:32) − S d (cid:15) d K (cid:33) , (10)where S is the distance from the transmission till the closest approach, K is the impact parameter of the light beam trajectory,and (cid:15) is the refractive bending experienced by the ray between the transmission and the reception .Those three quantities [cf. Eqs. (6), (9), and (10)] are three examples of observables that can directly be determined knowingthe evolution of ∇ S and t between the transmission and the reception points. However, as seen in Eqs. (1)–(5), they can onlybe determined after specifying the index of refraction’s dependencies. C. Refractive profile dependency
Except in Sec. IV where we remain as general as possible, in this work we will consider that the index of refraction coincideswith the surface of constant Φ , where Φ is a generalized potential which will be defined later in Sec. V.This hypothesis is motivated by the hydrostatic equilibrium assumption applied to a body which is made up with fluid ele-ments. Indeed, assuming a non-rotating fluid body, it can be demonstrated (cf. p. 91–92 of [25]) that the surfaces of constantdensity, pressure, and gravitational potential all coincide. This result can even be extended to a steady rotating body as shownin [25]. However in that case, if surfaces of constant density and pressure still coincide together they do not coincide anymorewith 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 contextof elementary theory of dispersion [19]), it becomes obvious that the surface of constant n also coincide with the surfaces ofconstant Φ . Such a dependence between n and Φ has also been used recently in [17] to analyze the radioscience data of Cassiniin the context of occultations by Saturn’s atmosphere.In order to simplify future notations, we introduce α ≡ − d n d Φ (11)which possesses the same dimensions as the inverse of a potential ( L − T ).In the context of GB observations, α must be seen as an input parameter describing the variation of the index of refractionin di ff erent layers of the atmosphere. Conversely, in the context of AO experiments it is seen as an output parameter whichmust be determined following an accurate modeling of the observable. In general, α is not a constant parameter and severalmeasurements are needed through the atmosphere in order to integrate Eq. (11), as it is discussed in [17]. D. Presentation of the work
In this work, we introduce a mathematical framework which provides a comprehensive description of the light ray trajectoryinside planetary atmospheres using geometrical optics.In Chap. III, we derive a simple solution to Eq. (4) assuming that the index of refraction of a planetary atmosphere is a linearfunction 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 Without loss of generality for the discussion, we neglect general relativity e ff ects and focus our attention on the refraction only. If needed, a more generictreatment of the Doppler shift in general relativity can be found in [20]. If (cid:15) can be derived from the cross product between ∇ S and ∇ S as (cid:15) = arcsin( | ∇ S × ∇ S | / n n ), in the following, we will use an alternative way forcomputing it. referred to as the reference solution . It is given in terms of constant parameters called hyperbolic elements (in reference to ellipticelements of celestial mechanics) which describe together the shape and the spatial orientation of the photon path as well as thetime-delay.In Chap. IV, we derive a comprehensive mathematical framework which allows to study generic refractive profiles. It isobtained using the reference solution derived in Sec. III as an osculating solution in order to turn Eqs. (1), (4), and (5) intoa set of oscultating equations describing the rate of change of the hyperbolic elements following a change in the index ofrefraction profile. The osculating equations (also called the perturbation equations) are perfectly equivalent to the equation ofoptics, therefore, in that sense, they can used in place of Eqs. (1), (4), and (5). The main interest is that they are well suited foranalytical study of light propagation inside planetary atmospheres.In Chap. V, we find approximate solutions to the perturbation equations for a large class of additional contributions in therefractive profile of the central planet’s atmosphere. First, we study the e ff ects on the hyperbolic elements due to non-lineardependencies of the refractive profile to the generalized potential. We focus on the limiting case where the generalized potentialreduces to the monopole term of the gravitational potential of the central planet (this reduces to the spherically symmetricassumption). 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, andfinally an external gravitational potential raised by a perturbing body. We also spotlights the versatility o ff ered by the formalismby studying the light dragging e ff ect caused by a rotating medium. All those perturbations are mainly related to the shape of theatmosphere. The last application concerns the e ff ect due to horizontal gradients inside a spherically symmetric atmosphere fora specific application to GB observations.Finally in Chap. VI, we compare the analytical results obtained for the axisymmetric part of the gravitational potential with asolution obtained from a numerical integration of the equations of geometrical optics (numerical ray-tracing). III. REFERENCE SOLUTION
The purpose of this chapter is to build a reference solution, for an idealized problem. As we will see in the next chapter, themain advantage of this reference solution lays in the fact that it is exact and simple, providing a good understanding of the e ff ectof 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 leaststringent is the monopole approximation which implies spherical symmetry. The most stringent is the constancy of α whichimposes a particular type of spherical symmetry. Each are discussed in turn. A. Spherical symmetry
In this section, we remind some general results which arise assuming a spherically symmetric refractive profile. The sphericalsymmetry implies that the index of refraction is only function of the magnitude of the separation vector, i.e. n = n ( h ) with h = | h | . Therefore, it is shown from Eqs. (1) and (4) (see also [19]) that the quantity which is defined as K = h × ∇ S (12)is constant all along the path of the light ray. Regarding the magnitude, we deduce the following relation which is known as theBouguer’s rule K = nh sin φ , (13)in which φ is the angle between the separation vector h and the unit tangent vector to the ray s (see Fig. 1). K = | K | is a constantvalue 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 inmagnitude and the light ray propagates into a plane which remains orthogonal to K . In such a case, it is helpful to describe thelight path with the help of the polar coordinates h and θ which are defined such that the components of the separation vector aregiven by h = h cos θ ˆ x + h sin θ ˆ y . (14) Physically, such gradients may be generated e.g. by horizontal changes in temperature. ˆ e ˆ ρ ˆ τ h θ φψ π − θω β β f P x ˆ y Line of nodes s -directionFIG. 1. Trajectory of the light ray as seen in the propagation plane (ˆ x , ˆ y ). P is the position of the central planet, 1 and 2 are respectively thepositions of a transmitter and a receiver along the light trajectory. The unit vector ˆ e points toward the direction of the closest approach. We have imposed that the (ˆ x , ˆ y ) plane coincides with the propagation plane. The ˆ z -direction is collinear to K , the ˆ x -axis isdirected along the line of nodes (defined in Sec. III C), and the ˆ y -axis is orthogonal to ˆ x , and ˆ z . We take the triad of vector(ˆ x , ˆ y , ˆ z ) to form a right-handed vectorial basis that we will refer to as the propagation basis since (ˆ x , ˆ y ) coincides with thepropagation 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 ˆ ρ ≡ h / h . The third member is ˆ σ ≡ K / K , whichis necessarily orthogonal to ˆ ρ and collinear to ˆ z . We notice that in the specific case where the index of refraction is a function ofthe height only, ˆ σ is a constant vector. The second member of the basis is ˆ τ which is defined to be orthogonal to ˆ ρ and ˆ σ . Wetake the triad of vectors ( ˆ ρ , ˆ τ , ˆ σ ) to form a right-handed vectorial basis. In the propagation plane, we have the relationsˆ ρ = cos θ ˆ x + sin θ ˆ y , ˆ τ = − sin θ ˆ x + cos θ ˆ y . (15) B. Monopole approximation and constancy of α In Sec. II C, we have discussed that assuming hydrostatic equilibrium together with elementary theory of dispersion, one caninfer that n = n ( Φ ) where Φ is a generalized potential containing the gravitational potential U of the central planet. The firstapproximation of the gravitational potential is the so-called monopole term (point-mass limit), which is given by the well-knownNewtonian potential U = − µ h , (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 refractivityprofile becomes spherically symmetric (see Sec. III A 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 refractiveprofile evolves linearly with the monopole term of the gravitational potential. This is the case when α is constant. Therefore, theexpression of the index of refraction is explicitly inferred from Eqs. (11) and (16) assuming α constant n ( h ) = η + αµ h . (17)In this expression, η is a constant representing the value of the index of refraction at infinity . 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 isspherically symmetric and evolves as 1 / h . 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 the relation η = n H − αµ/ H . We can now determine the quantities related to the description of the geometrical path of photons and the light time of the rayin 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 theray, both expressed in the propagation plane [making use of Eqs. (15)] h = p + e cos κ f ˆ ρ , (18a) ∇ S = √ ηκ (cid:114) αµ p (cid:2) e κ sin κ f ˆ ρ + (1 + e cos κ f ) ˆ τ (cid:3) . (18b)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. (A8)] 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. (A10a) and (A21)] are both linked to the conservation of the magnitude of theimpact parameter. Simultaneously, e is also linked to the conservation of E which is the dimensionless energy [cf. Eq. (A19)].In Eqs. (18), f is the true anomaly and is defined in Eq. (A10b) by means of an other constant of integration, ω , known asthe argument of the closest approach. The constancy of ω is linked to the conservation of the eccentricity vector defined inEq. (A22).The planar description is completed by using the fact that the index of refraction does not depend explicitly on the length orthe time. This gives rise to an additional constant of integration, S , being the geometrical length till the closest approach. Withthat constant we now assess a length description of the trajectory given the true anomaly by means of a Kepler-like equation[cf. Eq. (A30)]. In order to directly assess a time description instead of a geometrical length, we derived a second Kepler-likeequation for the time [cf. Eq. (A31)]. For that last case, the arbitrary constant of integration is, T , the elapsed time till the closestapproach.The solution is consequently described in the propagation plane by means of five constants of integration ( p , e , ω, S , T ). Itsvalidity has been assessed by comparison with a numerical integration of Eqs. (1), (4), and (5) for the index of refraction givenby Eq. (17). The hyperbolic elements remain constants at the level of the numerical noise. C. Hyperbolic path in space
The solution which is described in Sec. III B and appendix A achieves a remarkable degree of simplicity. By assuming thatthe index of refraction is a linear function of the Newtonian potential of the planet, we have shown that the light path is totallycontained inside a plane which remains orthogonal to the parameter vector. We have taken advantage of this particularity bychoosing for convenience the polar coordinates ( h , θ ) attached to the propagation plane. Then, we have demonstrated that theequation of light propagation possesses an exact solution which is an hyperbola contained inside that plane. However, a moregeneral description involving a more generic frame can be preferred. Indeed, when the reference solution is perturbed, thepropagation plane is no longer fixed in space, and the introduction of a fixed frame is required to describe the changes. Theprojection 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 ( ˆ X , ˆ Y , ˆ Z ). We adopt the ( ˆ X , ˆ Y ) plane as a reference plane and theˆ Z -direction as a reference direction. As before, the (ˆ x , ˆ y ) plane is the propagation plane while the ˆ z -direction is the orthogonaldirection to the propagation plane. The reference plane could be superimposed with the inertial equator of the central planet ata given time and the reference direction would be the direction of the pole. In that case, we impose that the two frames sharethe same origin. To achieve a complete description of the light path in the reference frame, we have to introduce two additionalangles as it can be seen from Fig. 2. The first one, labeled Ω , is called the longitude of the node and is the angle betweenthe ˆ X -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 ˆ z -direction orthogonal to the propagation plane and the ˆ Z -direction orthogonalto the reference plane. The line forming the intersection between the two planes is called the line of nodes . From Fig. 2, we caneasily determine the components of the unit vectors ˆ ρ , ˆ τ , and ˆ σ into the reference frame. The radial, tangent, and normal unitvectors are given byˆ ρ = + (cid:2) cos Ω cos( f + ω ) − cos ι sin Ω sin( f + ω ) (cid:3) ˆ X + (cid:2) sin Ω cos( f + ω ) + cos ι cos Ω sin( f + ω ) (cid:3) ˆ Y + sin ι sin( f + ω ) ˆ Z , (19a)ˆ τ = − (cid:2) cos Ω sin( f + ω ) + cos ι sin Ω cos( f + ω ) (cid:3) ˆ X − (cid:2) sin Ω sin( f + ω ) − cos ι cos Ω cos( f + ω ) (cid:3) ˆ Y + sin ι cos( f + ω ) ˆ Z , (19b)ˆ σ = + sin ι (sin Ω ˆ X − cos Ω ˆ Y ) + cos ι ˆ Z . (19c)Regarding the direction of the closest approach, we haveˆ e = + (cid:2) cos Ω cos ω − cos ι sin Ω sin ω (cid:3) ˆ X + (cid:2) sin Ω cos ω + cos ι cos Ω sin ω (cid:3) ˆ Y + sin ι sin ω ˆ Z . (20) ˆ X ˆ Y ˆ Z ˆ e ˆ σ ˆ ρ ˆ τ f Ω ω ιι β h P p r op a g a ti onp l a n e FIG. 2. Trajectory of the light ray as viewed from the fundamental reference frame ( ˆ X , ˆ Y , ˆ Z ). Thanks to these expressions, the height and the direction of the ray [cf. Eqs. (18)] can now be expressed into the referenceframe. We can also provide some simple definitions for the angles of the problem using the solution in the reference frame. Forinstance, from the previous equations, we deduce the following useful relationships K cos ι ≡ K · ˆ Z , (21a) − K sin ι cos Ω ≡ K · ˆ Y , (21b) K sin ι sin Ω ≡ K · ˆ X , (21c) K sin ι cos ω ≡ ( K × ˆ e ) · ˆ Z , (21d)sin ι sin ω ≡ ˆ e · ˆ Z . (21e)We have used the previous expressions of ˆ ρ , ˆ τ , ˆ σ , and ˆ e . These definitions are fundamental since they are expressed withvectors 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 .The complete solution described here and in the previous section is the analogous to the solution of the two body problem incelestial mechanics [25–28]. 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 ) forstudying much more complex problems, as discussed in Sec. IV. D. 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 befurther 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 onpurpose, and the reason will become obvious in Sec. IV B. Here, we are about to further simplify the previous expressionsderived in appendix A by using d s = (d h · d h ) / , which implies that s is a unit vector. Thus, all formulas which are presentedin this section are determined under this assumption and should only be applied when the index of refraction is exactly given byEq. (17). Although they should not be used in general, they provide a simplified and compact version of the formulas derived inappendix A. If we do not specify that s is a unit vector (like in appendix A), we count seven hyperbolic elements. Otherwise, one of the hyperbolic elements is notindependent anymore and we end up with six hyperbolic elements. This point is further developed in Sec. III D. 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)], andthe total number of independent constants of integration becomes six. After inserting the eikonal equation into Eq. (A19), wededuce E = η /
2, and thus, from Eq. (A21), one infers K = αµ e , (22)which shows that the eccentricity reduces to the inverse of the refractivity at K . This expression allows to express κ [cf. Eq. (A8)]as a unique function of the eccentricity κ = e − e . (23)From these two last equations, it is straightforward to show that the semi-latus rectum [cf. (A10a)] is expressed as p = αµη ( e − 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 thetangent to the ray h = κ K η (cid:18) αµ K + cos κ f (cid:19) − ˆ ρ , (25a) ∇ S = ηκ sin κ f ˆ ρ + ηκ (cid:18) αµ K + cos κ f (cid:19) ˆ τ . (25b)From Eqs. (22), (24), and (A13), it is easily shown that the semi-major axis, a , is now given by − a = αµη . (26)Insertion of Eqs. (22) and (26) into the two Kepler-like equations [cf. Eqs. (A30) and (A31)], let one to deduce s ( f ) = S + K η sinh F ( f ), (27)as well as t ( f ) = T + Kc sinh F ( f ) + αµ c F ( f ) + ( αµ ) cK f . (28) F is the hyperbolic anomaly [cf. Eqs. (A27)]. From these two relationships, we deduce that the two events ( ct , s ) and ( cT , S ) arelinked by the expression c ( t − T ) − η ( s − S ) = αµ F ( f ) + ( αµ ) K f . (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 ¯ e [cf. Eq. (A37)], we obtain ¯ e =
0, and we deducethe evolution of the argument of the bending angle ψ ( f ) = ω + f − (cid:114) e − e + (cid:32) κ f (cid:33) . (30)We can now examine the vacuum limiting case. As already mentioned, the parameter α characterizes the coupling betweenthe geometrical propagation of the light ray and the gravitational attraction of the central planet through the refractive propertiesof the medium. Thus, it is understood that if α →
0, the light ray should become insensitive to the presence of the mediumand should therefore propagate along a straight line trajectory. This can be easily checked by inserting the expression of p [cf.Eq. (24)] into Eq. (A9), and by taking the limit when α →
0. In this case, it is seen from Eqs. (22) and (26) thatlim α → e = + ∞ , lim α → a = − , (31)0also, we end up with lim α → h = K η cos f , (32)which represents, indeed, the equation of a straight line in polar coordinates (let us note that the vacuum limit is obtained for η = α → c ( t − T ) − η ( s − S ) =
0, which just reflects the fact that the twoevents ( 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 η = α → ψ = ω , (33)which shows that ψ remains constant whatever is the value of the true anomaly, meaning that the bending angle is null and thusthat the trajectory is a straight line. E. Applications
At this level, the formulas presented in Sec. III D can already be used either to model atmospheric delays or to process AOdata. In this section, we describe succinctly the procedure to simulate an atmospheric time-delay or to analyze Doppler datausing the hyperbolic solution described so far.The solutions in Eqs. (25), (28), and (30) allow us to compute explicit expressions of the observables [cf. Eqs. (6), (9), and(10)] using directly the reference solution.The expression of the atmospheric time-delay is immediate, since we only have to subtract the light time in a vacuum fromthe 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 α → τ atm = Kc (sinh F − sinh F ) − Kc (tan f − tan f ) + αµ c ( F − F ) + ( αµ ) cK ( f − f ), (34)where F = F ( f ) and F = F ( f ).The expression for the Doppler shift [cf. Eq. (9)] is also immediate since it is given by the scalar product between Eq. (25b)and ˆ u , where ˆ u is the coordinate velocity vector of the transmitter (subscript 1) and the receiver (subscript 2), divided by thespeed of light in vacuum. The substitution leads to ν ν = Γ Γ (cid:104) − ηκ ( ˆ ρ · ˆ u ) sin κ f − ηκ ( αµ K + cos κ f )( ˆ τ · ˆ u ) (cid:105)(cid:104) − ηκ ( ˆ ρ · ˆ u ) sin κ f − ηκ ( αµ K + cos κ f )( ˆ τ · ˆ u ) (cid:105) , (35)where ˆ ρ = ˆ ρ ( f ) and ˆ ρ = ˆ ρ ( f ), likewise for ˆ τ . If Eq. (35) is the most general form that can be obtained for the hyperbolicsolution, it is also interesting to express the Doppler frequency shift (see appendix C for the complete derivation) in term of thetotal refractive bending angle (cid:15) ≡ ψ ( f ) − ψ ( f ), which is easily inferred from Eq. (30).Finally, we can also find an expression for the refractive defocussing e ff ect [cf. Eq. (10)] by performing directly the di ff er-entiation of Eq. (30) with respect to K , remembering that the constants of integration, e.g. the eccentricity, are function of theimpact parameter. Such an expression may be tedious to derive, thus a numerical evaluation of the derivative d (cid:15)/ d K from finitedi ff erences may be preferred.From the evolution of the observables, we can now think of an easy procedure to apply to simulate delays and to analyzerange-rate measurements. However, we remind once again that all the analytical computations made in Secs. III D and III E havebeen carried out assuming that the index of refraction is exactly given by Eq. (17). Furthermore, we have seen in Sec. III B thatthis refractive profile corresponds to a very peculiar type of spherical symmetry ( α = cst). So, in order to extend the referencesolution to any kind of spherical symmetries ( α (cid:44) cst), we are to assume (only for this discussion) that the atmosphere understudy is made with m concentric layers, where each layer possesses its own value of α k with k = , . . . , m . In addition, insideeach layer the index of refraction follows Eq. (17) where the value of η k must be chosen in such a way that it connects the valuesof the index of refraction at the interfaces between two layers.For GB observations, values of α k can be assumed or deduced from di ff erent sources. For instance, we can use the profileof temperature and pressure of the International Standard Atmosphere [29], or we can apply the same procedure as the onedescribed in [2] or [3]. For a given geometry, we can thus determine the value of the hyperbolic elements inside each layer. Theatmospheric time-delay inside each layer is simply given by Eq. (34), and finally the total delay is the sum of all the contributions.Conversely, the purpose of AO experiments is to deduce the value of α k at di ff erent heights inside the atmosphere. Theobservables are the Doppler frequency shift measurements, for which the analytical counterpart is given in Eq. (35). It dependson the evolution of ∇ S which relies on α k . Therefore, the di ff erence between the observed and the computed Doppler shiftcan be minimized by varying the value of α k . At the end, each layer provides its own determination of α k . Therefore, n can beretrieved overall the vertical profile by integrating Eq. (11).1 IV. PERTURBATION EQUATIONS
In the case of a spherically symmetric refractivity, we saw that the direction of K [see Eq. (12)] remains the same while thebeam ray propagates through the atmosphere. However, if the gradient of the index of refraction have non-null transverse ornormal components, the vector K does not remain constant anymore and the previous light path which used to occur in a fixedplane is no longer a solution.In this chapter, we address the problem of finding a new solution to a more complex situation where the index of refraction isas much general as possible. A. 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 thefollowing notation n = n + δ n (36)with n being the index of refraction given in Eq. (17) and δ n being a new non-constant contribution . Therefore, this newcontribution will modify the second member of Eq. (4) which is now given by a much more general expression than in Eq. (A3) ∇ n = − α g ˆ ρ + f pert . (37)Here, g still refers to the magnitude of the local spherical acceleration [cf. Eq. (A4)], and the first term on the right-handside is the gradient of the index of refraction which produces the hyperbolic solution for the light path. The additional term, f pert , represents a supplementary spatial variation of the index of refraction and is given by f pert = ∇ δ n . We do not make anyassumption about its magnitude with respect to the hyperbolic part, and we do not specify its dependency. In addition, we releasethe isopotential surfaces hypothesis. We only assume that the perturbing gradient can be decomposed in a certain rotating basissuch that its radial, transverse, and normal components are given by f pert = R ˆ ρ + T ˆ τ + S ˆ σ . (38)Here it must be pointed out that the basis ( ˆ ρ , ˆ τ , ˆ σ ) has a di ff erent meaning with respect to the one introduced in Chap. III.Indeed, as discussed above, if f pert = , 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 ˆ σ ≡ K / K does notremain fixed in space as before. Good insight can be obtained on the e ff ects of the perturbing gradient on the impact parametervector by di ff erentiating Eq. (12) with respect to s . Making use of Eqs. (4), (1), and (37), we quickly arrive to d K / d s = h × f pert ,which can be expressed in terms of the perturbing gradient’s componentsd K d s = h ( T ˆ σ − S ˆ τ ). (39)Then, a comparison between Eq. (39) and the derivative of ˆ σ with respect to s reveals that the magnitude of K changes accordingto d K d s = h T . (40)From these two last equations, we predict that a normal perturbation will only a ff ect the direction of the impact parametervector while a transverse one will also change its magnitude. These equations are perfectly equivalents to the ones of celestialmechanics where only the non-radial components of the perturbing acceleration change the angular momentum vector [25, 30].We can also demonstrate that K is not the only first integral being a ff ected by the perturbing gradient. Indeed, as seen fromthe form of the relationships in Eqs. (A16) and (A22), we can expect that changes in the index of refraction are going to impactthe dimensionless energy parameter and the eccentricity vector. Then, since all the hyperbolic elements ( p , e , ι, Ω , ω, S , T ) aredetermined from K , E , and e (see appendix. A), they are also expected for not being constant anymore. The determination oftheir evolution is the subject of the next two sections. If δ n is a pure constant, it must be treated as a new limit value at infinity, such η (cid:48) = η + δ n . B. Method of variation of arbitrary constants
In this section, we will refer to the general solution of the unperturbed problem as h = h (cid:0) s , C (cid:1) , ∇ S = ∇ S (cid:0) s , C (cid:1) , t = t (cid:0) s , C (cid:1) , (41)where h , ∇ S , and t are the solutions respectively for the separation vector, the direction, and the light time along the rayfor an hyperbolic path, i.e. for n . The C -vector represents the hyperbolic elements which are constants for the unperturbedtrajectory. These solutions have been fully described in Chap. III.We are now interested in studying the very general case where the perturbing gradient is non-null. In this context, the methodof variation of arbitrary constants is very well adapted (see e.g. [25] for application of the method of variation of arbitraryconstants in the context of celestial mechanics). The core of the method is to consider that h and ∇ S [cf. Eqs. (18)] arestill solutions of the perturbed problem, avoiding the apparent contradiction by allowing the constants of integration (hyperbolicelements) to change as s evolves. The physical meaning is that at any length s along the ray, the trajectory is taken to behyperbolic with the parameters C ( s ). The general solution must be expressed as h = h (cid:0) s , C ( s ) (cid:1) , ∇ S = ∇ S (cid:0) s , C ( s ) (cid:1) . (42)A subtle consequence of this choice, is that s and s cannot be two unit vectors at the same time. Indeed, from the definitionof ∇ S in Eq. (2) and from (42), we notice that s · s = ( n / n ) s · s . Moreover, considering that the eikonal equation (3) musthold true for the index of refraction under study (in this case n ), we infer that s is a unit vector, thus s is not. Consequently, fromnow, we are not allowed to use the hyperbolic expressions derived in Secs. III D and III E, which are based on the assumptionthat s is a unit vector. However, we remind that the solution which is derived in appendix A is free of this assumption, then allthe expressions in there can be used.Another important consequence is that, saving the expression of ∇ S [cf. Eq. (18b)], allows us to fix the expression of theDoppler frequency shift once and for all. Indeed, Eq. (9) depends only of ∇ S , and if Eq. (18b) must hold by virtue of themethod of variation of arbitrary constant, then the expression of the Doppler frequency shift [as inferred using Eq. (18b)] mustalso hold true for any arbitrary index of refraction. However, since the hyperbolic elements are not constants anymore, theirvalues at the reception point might be di ff erent from the ones at the transmission point.Applying the proposition (42), we can determine an explicit expression for the variation of the dimensionless hyperbolicenergy E . Di ff erentiating Eq. (A16) and making use of Eqs. (38), and (18), we get the following expressiond E d s = √ η (cid:114) αµ p (cid:104) e R sin κ f + κ (1 + e cos κ f ) T − e αµ p (1 + e cos κ f ) N sin κ f (cid:105) , (43)in which N = δ nn p − (44)is a new perturbing component. The term δ n has been introduced in Eq. (36) and refers to the di ff erence between the actualand the hyperbolic index of refraction. We recall that δ n must contain only the variable part of the di ff erence as discussed infootnote 7. We see from Eq. (43) that only the components of the perturbing gradient which are contained inside the propagationplane 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 infertheir length rate of change as functions of the perturbing gradient’s components using Eqs. (39), (40), and (43). However, wefirst need to express the di ff erentials of the hyperbolic elements as functions of d E / d s , d K / d s , and d e / d s .Let us start with the semi-major axis. Di ff erentiation of Eq. (A18) yieldsd a d s = αµη E d E d s . (45)The equation for the rate of change of the eccentricity is given by di ff erentiating Eq. (A21). After some little algebra we findd e d s = K Ee ( αµη ) (cid:32) K d K d s + κ E d E d s (cid:33) . (46)The variation of the semi-latus rectum is determined either from Eqs. (45) and (46) making use of Eq. (A13), or directly bydi ff erentiating Eq. (A10a). Both ways lead to d p d s = K αµη d K d s . (47)3 ˆ X ˆ Y ˆ Z Ω δ Ω θ θ + δθι ι + δι h FIG. 3. Illustration of the e ff ect on the θ angle due a tilt δ Ω along the longitude of the node. The darker plane corresponds to the propagationplane before the tilt while the lighter one is the same plane after the tilt. From that sketch, we infer the relation ( θ + δθ ) cos δι + δ Ω cos ι = θ ,which reduces to δθ = − δ Ω cos ι for infinitesimal tilts. The variation of the inclination is given after di ff erentiating Eq. (21a)d ι d s = − K sin ι (cid:32) d K d s · ˆ Z − cos ι d K d s (cid:33) . (48)Eqs. (21b) and (21c) let to determine the tangent of the longitude of the node. After di ff erentiation, we arrive tod Ω d s = K sin ι (cid:32) cos Ω d K d s · ˆ X + sin Ω d K d s · ˆ Y (cid:33) . (49)For the derivation of the rate of change of the argument of the closest approach, we do not di ff erentiate the eccentricity vectoras it is usually done in celestial mechanics [25]. Instead, we are to simplify the algebra by using a geometrical relation [30]. Thecomputation starts by di ff erentiating Eq. (A9) and by making use of Eqs. (A10), (A8), (47), and (46), which leads tod ω d s = − cos ι d Ω d s − κ e K ( αµη ) cot κ f d E d s + e κ K αµη csc κ f (cid:18) h − + e η κ ( αµ ) K f sin κ f − Ee αµη cos κ f (cid:19) d K d s . (50)On the right-hand side we have substituted the following relationshipd θ d s − ˆ τ · ∇ S nh = − cos ι d Ω d s (51)which has been determined geometrically from Fig. 3. The left-hand side is obtained by making use of the method of variationof arbitrary constants [cf. Eqs. (42)]. The right-hand side contribution is determined remembering that for a certain locationalong the ray, h remains unchanged, so we conclude that a tilt along Ω will decrease θ by an amount equal to the projection ofthe 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 [30]. With this geometrical relation, the expression of thechange 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 photonpath. The relation is obtained by di ff erentiating Eq. (A10b) and by substituting it into Eq. (51)d f d s − ˆ τ · ∇ S nh = − (cid:32) d ω d s + cos ι d Ω d s (cid:33) . (52)The left-hand side represents the hyperbolic contribution to the change in the true anomaly. The right-hand side represents thechange in the direction of the closest approach relative to the ˆ X -direction and measured inside the propagation plane (see Fig. 3).The same is also true in celestial mechanics as mentioned by [25]. Substituting Eq. (50) into Eq. (52) provides the followingexpression4d f d s = √ ηκ n (cid:114) αµ p (1 + e cos κ f ) + κ e K ( αµη ) cot κ f d E d s − e κ K αµη csc κ f (cid:18) h − + e η κ ( αµ ) K f sin κ f − Ee αµη cos κ f (cid:19) d K d s . (53)The right-hand side reduces to ˆ τ · ∇ S / nh when the perturbation is null (i.e. when d E / d s and d K / d s vanish). C. Osculating equations
We are now ready to derive the final expressions for the evolution of the hyperbolic elements as functions of the componentsof the perturbing gradient of the index of refraction. However, for a matter of compactness, we do not provide the di ff erentialsof ( S , T ), but rather the di ff erentials 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 elementswith d E / d s , d K / d s , and d K / d s , it is possible to infer the following relationshipsd p d s = κ √ η (cid:115) p αµ + e cos κ f T , (54a)d e d s = √ η (cid:114) p αµ (cid:40) R sin κ f + κ (cid:34) κ f + e (1 + cos κ f )1 + e cos κ f (cid:35) T − αµ p (1 + e cos κ f ) N sin κ f (cid:41) , (54b)d ι d s = κ √ η (cid:114) p αµ cos( f + ω )1 + e cos κ f S , (54c)d Ω d s = κ √ η (cid:114) p αµ sin( f + ω )1 + e cos κ f S csc ι , (54d)d ω d s = − e κ √ η (cid:114) p αµ (cid:40) R cos κ f − κ (cid:34) + e cos κ f + e cos κ f sin κ f + e κ (1 − κ ) f + e cos κ f (cid:35) T + e κ sin( f + ω )1 + e cos κ f S cot ι − αµ p (1 + e cos κ f ) N cos κ f (cid:41) , (54e)d f d s = √ ηκ n (cid:114) αµ p (1 + e cos κ f ) + e κ √ η (cid:114) p αµ (cid:40) R cos κ f − κ (cid:34) + e cos κ f + e cos κ f sin κ f + e κ (1 − κ ) f + e cos κ f (cid:35) T − αµ p (1 + e cos κ f ) N cos κ f (cid:41) , (54f)d t d s = nc . (54g)These equations provide a good understanding of how the components of the perturbing gradient govern the length rate of changeof the hyperbolic elements. For instance, it can be seen that p is only a ff ected by the transverse component of the perturbinggradient, where ι and Ω only change because of the normal component. We also see that e and f are a ff ected by the radial andthe transverse components, therefore any perturbation which is contained inside the propagation plane will induce variationsof these two elements. In addition, they are also expected to vary because of a term in δ n (through N ), which represents anon-constant change in the expression of the index of refraction. Finally, we see that ω is the only one being a ff ected by all thecomponents 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 di ff erence stands in the terms in N which does not possess any equivalent in celestial mechanics. Thereason is that the motion of planets is not directly sensitive to the gravitational potential itself but rather to its gradient. In thecase of geometrical optics, the equations of light propagation not only contain the gradient of the index of refraction, but alsothe value of that index, see e.g. Eq. (5). This can also be inferred once Eq. (4) is written as d s / d s = n − [ ∇ n − s ( s · ∇ ) n ], whichrepresents the rate of change of the curvature of the ray [19].As in celestial mechanics, when the inclination is null, ω and Ω are not defined since the direction of the line of nodes is notdetermined as seen in Fig 2. These singularities are related to the choice we have made in the definition of angles. They can beremoved [25, 27, 28, 31] by introducing new angles, such the longitude of the periapsis which is defined as (cid:36) = ω + Ω . In thenext, 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 thatsense, they are generic and no approximations have been introduced in the transcription. Their validity has been assessed by5comparing numerical integration of Eqs. (54) with numerical integration of Eqs. (1), (4), and (5) showing that the di ff erencesremain at the level of the numerical noise. The main advantage of the perturbation equations is to provide an easy geometricalpicture to tackle the e ff ects of a perturbing gradient with respect to an hyperbolic path.It is worth noticing that any exact solution to the fundamental equations of optics could have been used as a reference solution,or osculating solution. For instance, we could have chosen the straight line as the osculating solution instead of the hyperbolicpath. Actually, if we insert Eqs. (22)–(24) into Eqs. (54), and if we take the limit α → e → ∞ ), the perturbationequations 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 consideredsmall with respect to the spherical contribution of the gravitational potential of the central planet. D. First-order approximation
The great benefit of having written the equation of fundamental optics in terms of the components of the perturbing gradientwill become obvious in this section.Indeed, in the case where the perturbing gradient is small with respect to the hyperbolic contribution (i.e. α g (cid:29) | f pert | ), thechanges in the hyperbolic elements are expected to be small, too. Thus, a good understanding can be reached by inserting theconstant values of the hyperbolic elements in the right-hand side of Eqs. (54) in order to get the first-order approximation. Forthis, it is convenient to use the true anomaly as independent variable instead of the geometrical path length. This requires toinvert Eq. (54f) keeping only the linear order in the components of the perturbing gradient. We find the following relationsd p d f (cid:39) n η p αµ + e cos κ f ) T , (55a)d e d f (cid:39) κ n η p αµ (cid:40) sin κ f (1 + e cos κ f ) R + κ (cid:34) κ f + e (1 + cos κ f )(1 + e cos κ f ) (cid:35) T − αµ p N sin κ f (cid:41) , (55b)d ι d f (cid:39) κ n η p αµ cos( f + ω )(1 + e cos κ f ) S , (55c)d Ω d f (cid:39) κ n η p αµ sin( f + ω )(1 + e cos κ f ) S csc ι , (55d)d ω d f (cid:39) − n η e p αµ (cid:40) cos κ f (1 + e cos κ f ) R − κ (cid:34) + e cos κ f (1 + e cos κ f ) sin κ f + e κ (1 − κ ) f (1 + e cos κ f ) (cid:35) T + e κ sin( f + ω )(1 + e cos κ f ) S cot ι − αµ p N cos κ f (cid:41) , (55e)d δ s d f (cid:39) − κ e n (cid:112) η (cid:115) p ( αµ ) (cid:40) cos κ f (1 + e cos κ f ) R − κ (cid:34) + e cos κ f (1 + e cos κ f ) sin κ f + e κ (1 − κ ) f (1 + e cos κ f ) (cid:35) T− αµ p ( e η/ n + cos κ f )(1 + e cos κ f ) N (cid:41) , (55f)d δ t d f (cid:39) − κ ec n (cid:112) η (cid:115) p ( αµ ) (cid:40) cos κ f (1 + e cos κ f ) R − κ (cid:34) + e cos κ f (1 + e cos κ f ) sin κ f + e κ (1 − κ ) f (1 + e cos κ f ) (cid:35) T− αµ p (2 e η/ n + cos κ f )(1 + e cos κ f ) N (cid:41) , (55g)In the two last equations, we have introduced the non-hyperbolic contributions to the geometrical length and the light time, suchthat δ s = s − s and δ t = t − t , where s and t are the hyperbolic contributions which are given explicitly in Eqs. (A30) and(A31), respectively.We can also add an expression for the evolution of the argument of the refractive bending which is needed to determine thedefocussing and which can also be used to approximate the Doppler frequency shift (cf. discussion in Sec. III E). To do so, onecan di ff erentiate Eqs. (13) and (A32), and then makes use of Eqs. (A33) to deduced ψ = − hn d n d h (cid:32) + d ω d f − h d K d n d s d f (cid:33) d f ,6which is more general than Eq. (A34) since it includes, besides the hyperbolic term, the variations of K and ω . Substituting thedi ff erentials into that equation and keeping only the first order in the components of the perturbation, we end up withd δψ d f (cid:39) − p η e (cid:40) ( e η/ n + cos κ f )(1 + e cos κ f ) R − κ (cid:34) + e cos κ f (1 + e cos κ f ) sin κ f + e κ (1 − κ ) f (1 + e cos κ f ) + e κ η n sin κ f (1 + e cos κ f ) (cid:35) T + e κ sin( f + ω )(1 + e cos κ f ) S cot ι + αµ p ( e η/ n − cos κ f )(1 + e cos κ f ) N (cid:41) , (55h)where we have removed the hyperbolic contribution, ψ , which is given exactly in Eq. (A36). We have defined δψ = ψ − ψ asthe non-hyperbolic contribution to the argument of the refractive bending. Once integrated, that last equation will give us the netchange in the refractive bending due to a perturbation and then, Eq. (C2) gives us directly the first order e ff ect on the Dopplerfrequency shift. E. Method of integration
Eqs. (55) are general enough to be applied to di ff erent geometries and di ff erent kind of problems as AO experiments or GBobservations. In both cases, the integration can be performed along the hyperbolic light path to get the first order e ff ects on thehyperbolic elements. Calling C the vector of the following elements, C = ( p , e , ι, Ω , ω, δ s , δ t , δψ ), we have ∆ C ≡ C ( f ) − C ( f ) = (cid:90) f f d C d f (cid:48) d f (cid:48) , (56)where f and f are respectively the true anomaly at the transmission and at the reception point, with f ≥ f . The term d C / d f (cid:48) 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), aTaylor expansion in term of these quantities can first of all let to simplify the integration and secondly lead to a more user-friendlysolution.For instance, in the context of AO experiments and GB observations, we can think of the inverse of the eccentricity as a smallquantity. Indeed, in most experiments, we have to deal with a small refractivity which implies α (cid:28) e (cid:29) K (cid:29) αµ ). The meaning is that the shape of the hyperbola departs slightly from a straight-linetrajectory. However, if this approximation is largely relevant for AO, it can become somehow inaccurate for GB experiments,in particular at very high elevations (close to the zenith-direction of the site). Indeed, K might become the same order as αµ ,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 (cid:46)
10, we mustconsider only trajectories with K (cid:38)
30 km, which corresponds to elevations angles (cid:46) . ◦ .For GB observations at very high elevations, the light ray trajectory is close to be radial, that is to say that the change in thetrue anomaly between the transmitter and the receiver ( f − f ) is a small quantity. For that particular case, we do not need tointegrate Eqs. (56), instead an evaluation of Eqs. (55) at the level of the transmitter is su ffi cient to determine the change in thehyperbolic elements at first order in ( f − f ) ∆ C = ( f − f ) d C d f ( f ). (57)Therefore, for all the following applications (see Chap. V), we consider the leading order in 1 / e since: i) it encompasses almostall experimental cases (usually, at zenith angles less than few degrees, source tracking can be di ffi cult), ii) the transposition toGB observations at very high elevations (zenith angles less than few degrees) reduces to an evaluation of Eqs. (55) at the levelof the transmitter. V. APPLICATIONS
In this chapter, di ff erent 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 theobservables can be easily inferred as discussed in Sec. II B. Following [32], for an average parcel of air at sea level, the refractivity is approximately (3 ± × − , which represents a value of α = (5 ± × − km − · s .So, the region around the zenith-direction, in which the eccentricity is e (cid:46)
10, is enclosed inside a cone with a half top angle which remains (cid:46) . ◦ . All thelight path trajectories with K (cid:46)
30 km, are enclosed inside that region. ff ect on the light propagation due to themoving medium. We also apply the perturbation equations considering the non-spherical part of the self-gravitational potentialof the central planet. Then, we study the perturbation on the light beam trajectory caused by a perturbing body raising tides onthe central planet’s atmosphere. Finally, we close the section by studying the e ff ects due to the presence of horizontal gradientsin a spherically shaped atmosphere. A. 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 ageneric refractive profile assuming the influence of di ff erent sources of stress.Lately, we have made use of three reference frames, all centered at the planet’s center of mass and referred to as the propagationframe (ˆ x , ˆ y , ˆ z ), the polar basis ( ˆ ρ , ˆ τ , ˆ σ ), and the reference frame ( ˆ X , ˆ Y , ˆ Z ). The first one was helpful to derive the solution ofreference (hyperbolic path) by making use of the second one to introduce polar coordinates. The last one was introduced toinfer the perturbation equations in an arbitrary oriented frame. It was assumed that the medium was at rest simultaneously inthe propagation frame and in the frame of reference, also the equations of geometrical optics [cf. Eq. (4)] were consequentlyavailable 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 themedium. 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 [33]. Following [23] and [34], the equation of geometrical opticsexpressed in the medium comoving frame is dd s (cid:0) ∇ S (cid:1) = ∇ n + f drag , (58a)where n still refers to the value of the index of refraction as measured by an observer at rest in the refractive medium. Theadditional term refers to a dragging contribution which is given by f drag = c (cid:104) ( ∇ S · ∇ n ) v − ( v · ∇ S ) ∇ n (cid:105) + γ nc ∇ S × ( ∇ × v ) + O (cid:16) v / c (cid:17) , (58b)where v is the velocity vector of the medium, with v = | v | . γ is the Fresnel’s dragging coe ffi cient introduced in Eq. (A20). Itmight be seen from this expression, that the corrections to Eq. (4) are of the order of v / c . For gas giant planets in the SolarSystem (e.g. Jupiter), the typical upper bound value of the velocity for the solid rotation at the equator is 5 × km / h withzonal winds asymmetric velocity ranging between ±
550 km / h [35], also v / c (cid:46) − . Therefore, we will consider the e ff ect ofthe moving medium as a perturbation to the hyperbolic path, f pert = f drag . (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 ( ˆ X (cid:48) , ˆ Y (cid:48) , ˆ Z (cid:48) ) centered at the planet’s center of mass androtating with it (we call P the central planet). The ˆ Z (cid:48) -axis is chosen to be orthogonal to the equatorial plane of P . The angularvelocity vector is considered to be constant and given by w = w ˆ Z (cid:48) . Since w is independent of time, the ˆ Z (cid:48) -axis is spatially fixed.Therefore, we choose the ˆ Z -axis of the reference frame to be collinear with the ˆ Z (cid:48) -axis. The frame ( ˆ X (cid:48) , ˆ Y (cid:48) , ˆ Z (cid:48) ) is well suited forthe study of internal motions of P and will be referred to as the fluid rotating frame.It is shown in [25] (cf. p. 106) that Euler’s equation describing the time evolution of the velocity of a fluid element makesappear the following generalized potential Φ = U + (cid:88) l ≥ U l − Φ C + U tidal ( t ), (60)once written in the fluid rotating frame. In the following, we consider that the dominant term of that expression is the monopoleterm U of the gravitational potential of the fluid planet. The other terms will be considered as perturbations before U . The U l In the Solar System, the internal problem (study of the internal structure of a self-gravitating bodies) can be mainly decoupled from the external one (inter-bodydynamics). This is true as long as the inter-body distance is much more higher than a characteristic length scale within the extended bodies. Consequently,an external stress can be considered as a perturbation in the internal problem and conversely, an internal stress may be seen as a perturbation in the externalproblem. We took the convention that the Newtonian gravitational potential is a negative quantity, also the signs in that expression di ff ers from [25]. Φ C is the centrifugal potential due to planet’sproper rotation, and U tidal ( t ) is the external potential due to the presence of other massive bodies in the surrounding of the centralplanet. This last term is dynamic since it evolves with the positions of the perturbing bodies. Each contribution will be studiedin turn.We now have to define an expression for the index of refraction. If we assume that the refractive index is directly related tothe density of the fluid, we can consider that Eq. (60) constitutes the generalized potential in Eq. (11). In addition, successiveintegrations by part of Eq. (11) leads to the following infinite series n = η − α Φ + + ∞ (cid:88) k = ( − k k ! α ( k ) Φ k , (61)where Φ is given by Eq. (60), and α ( k ) ≡ d k n / d Φ k for k ≥
2. We keep the same notation as in previous chapters, i.e. α ≡ α (1) isstill the linear term with the generalized potential. This series can represent any function of the gravitational potential, and thechoice is governed by the numerical values which are assigned to the α k coe ffi cients.Eqs. (58)–(61) are all we need to formulate the internal problem and study first order e ff ect on the light propagation due tothe di ff erent pieces in Eqs. (60) and (61). B. Non-linearity
Previously in Chap. III, we have assumed d n / d Φ = cst [cf. Eq. (11)] in order to simplify the integration of Eq. (4). However,this could be too restrictive in some applications where it is not possible to build a multi-layers modeling of the atmosphere. Ifa multi-layers modeling can be achieved as in [17] or in [6], a constant value of α can be associated to each layer and the totalprofile can be recovered by integrating Eq. (11).In this section, we propose to explore the first order e ff ect of non-linearity with isopotentials in the refractive profile [seeEq. (61)]. We focus on the simple case where the generalized potential is given by the monopole term of the Newtonianpotential, Φ = U . All the mathematical details are exposed in Sec. B 1 and the change in hyperbolic elements for any degree k is given at leading order in 1 / e in Eq. (B3). The application to the quadratic contribution ( k =
2) is given in Eq. (B4).After integration, it is shown that whatever function of the Newtonian potential the index of refraction is, p , ι , and Ω remainalways constants regardless degree k . This is due to the fact that the index of refraction acts like a central field with no transverseor 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 = ∆ δ t = η α (2) µ cK ( f − f ), ∆ δ(cid:15) = η α (2) µ K ( f − f + cos f sin f − cos f sin f ) .Obviously, from an observational point of view, these expressions are of a great interest. The first one provide directly theadditional time-delay with respect to the hyperbolic path, and the second one provides the supplementary contribution to thehyperbolic refractive bending. Both are function of the geometry of the problem only. In addition, from Eq. (C2), it can be seenthat the later provides directly the change in the Doppler measurement due to α (2) . C. Steady rotating atmosphere
The second application concerns the rotation of the atmosphere. All bodies in the Solar System are rotating. The rotationmotion 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 andconsiders a steady rotating planet.In the frame rotating with the fluid, the media experiences the e ff ect of the proper rotation via the simplified followinggeneralized potential Φ = U − Φ C , where Φ C represents the centrifugal potential. In addition, the dragging e ff ect due to therotational motion of the medium is described by the additional contribution f drag in Eqs. (58). The solutions describing the changein hyperbolic elements due to a steady rotating atmosphere are given in Eqs. (B8) and (B11). All the mathematical details arepresented in Sec. B 2.9On one hand, Eqs. (B8) describe the change in hyperbolic elements due to a centrifugal potential induced by the properrotation of the planet. The expressions of the non-hyperbolic contribution to the light time and the refractive bending are ∆ δ t = − α w K c η sin ι (cid:104) ω (cid:16) tan f − tan f (cid:17) + ω (cid:16) sec f (cid:2) + f (cid:3) − sec f (cid:2) + f (cid:3)(cid:17) (cid:105) , ∆ δ(cid:15) = α w K η (cid:104) sin 2 ω sin ι (cid:16) sec f − sec f (cid:17) − (cid:0) + cos 2 ι + ω sin ι (cid:1)(cid:0) tan f − tan f (cid:1)(cid:105) ,with f ≥ f .On the other hand, Eqs. (B11) describe the change in hyperbolic elements induced by the light dragging e ff ect due to therotation of the fluid material which composes the atmosphere of the central planet. The expressions of the non-hyperboliccontribution to the light time and the refractive bending are ∆ δ t = − γ K wc (cid:104) (2 + cos 2 f ) sec f tan f − (2 + cos 2 f ) sec f tan f (cid:105) cos ι , ∆ δ(cid:15) = − γ Kwc (cid:0) tan f − tan f (cid:1) cos ι ,with f ≥ f .These relationships show that the first set is quadratic with the magnitude of the angular velocity vector while the latter evolveslinearly with it. Therefore, the change in hyperbolic elements due to the centrifugal potential is independent of the orientationof w and remains the same for either a direct or an indirect rotation. On the contrary, the change in hyperbolic elements due tothe dragging of light is dependent of the direction of the medium’s rotation. Because of a cosine of the inclination, we can alsosee that the e ff ect is maximum in the planetary equator, where the distance with respect to the axis of rotation is maximum.Another important di ff erence between the e ff ects due to the centrifugal potential and those due to the dragging of light canbe inferred by comparing the expressions of the non-hyperbolic contribution to the light time and the path length. In the caseof light dragging, it is seen in Eq. (B11g) that the non-hyperbolic contribution to the light time reduces to the non-hyperboliccontribution to the geometrical length divided by the speed of light in vacuum. This comes from the fact that the N perturbingcomponent is null for the light dragging e ff ect while it is not for the centrifugal potential contribution. D. Axisymmetric gravitational potential
Because of their proper rotation, non-rigid bodies tend to be flattened due to centrifugal forces. Then, the mass repartitionbecomes slightly di ff erent from what would be expected in spherical symmetry, and consequently, the gravitational potentialalso changes. This mass redistribution due to centrifugal forces is usually the most important departure from the monopolecontribution to the total self-gravitational potential of the planet.The non-spheric components of the gravitational potential of the central planet can be modeled using the simplified followinggeneralized potential Φ = U + (cid:80) l U l with l ≥
2. The solutions for the change in the hyperbolic elements are derived in Sec. B 3and the results for l = f ≥ f ∆ δ t = η J αµ R cK (cid:104) sin ι (cid:16) ω (cid:2) cos f − cos f (cid:3) + (cid:2) sin(3 f + ω ) − sin(3 f + ω ) (cid:3)(cid:17) + (cid:16) − sin ι (cid:2) − ω (cid:3)(cid:17)(cid:0) sin f − sin f (cid:1)(cid:105) , ∆ δ(cid:15) = η J αµ R K (cid:104)(cid:0) + ι (cid:1)(cid:0) sin 3 f − sin 3 f (cid:1) + (cid:16) − sin ι (cid:2) − ω (cid:3)(cid:17)(cid:0) sin f − sin f (cid:1) + ι (cid:16)(cid:2) + f (cid:3)(cid:2) f + ω ) + sin(3 f + ω ) (cid:3) − (cid:2) + f (cid:3)(cid:2) f + ω ) + sin(3 f + ω ) (cid:3)(cid:17)(cid:105) .Considering that the centrifugal potential and the quadrupole moment of the gravitational potential of the central planet areboth axisymmetric fields, one might expect similar signatures in the changes induced on the hyperbolic elements. However,because they evolve di ff erently with h (the centrifugal e ff ect tends to grow with h , where the J e ff ect decreases with h ), acomparison shows that the signatures produced on the hyperbolic elements di ff er. However, one can see that the ratio betweenEqs. (B21) and (B8) is always proportional to C J C C ∝ µ J R η w K ,0where C represents the set of hyperbolic elements. The subscript J refers to the e ff ect due to the quadrupolar moment of thegravitational 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 ratiois 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 themost important. However, it must be noticed that this ratio can grow while the light ray goes deeper inside the atmosphere, sincethe impact parameter decreases. This can make the J e ff ect being more important than the centrifugal one, especially for Titanfor 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. E. Tidal potential
In systems containing close orbiting bodies, e.g. a planet and its satellites, the question of knowing whether the tidal e ff ectdue to a perturbing body on the planetary atmosphere can a ff ect 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 rotationaldeformation and the non-spherical gravitational potential contribution which were treated previously. Hence, we consider thesimplified following generalized potential Φ = U + U tidal ( t ). In addition, we will assume that the characteristic time of variationof the external tidal field is so slow that it never takes the central planet’s atmosphere out of hydrostatic equilibrium. Thisassumption is known as the static tides approximation.All the computations are detailed in Sec. B 4, and the evolution of the hyperbolic elements are given in Eqs. (B35). We remindthat these equations are derived considering a tide raising body moving on an equatorial circular orbit around the central planet,and a photon path lying inside the equatorial plane (i.e. ι = f ≥ f ∆ δ t = − η αµ K cr µ A µ (cid:110) sec f (cid:2)
15 sin(5 f + Nt − (cid:36) ) +
15 sin(5 f − Nt + (cid:36) ) + f + Nt − (cid:36) ) + f − Nt + (cid:36) ) +
75 sin(3 f + Nt − (cid:36) ) +
75 sin(3 f − Nt + (cid:36) ) +
35 sin(3 f + Nt − (cid:36) ) +
35 sin(3 f − Nt + (cid:36) ) +
150 sin( f + Nt − (cid:36) ) −
150 sin( f − Nt + (cid:36) ) +
70 sin( f + Nt − (cid:36) ) +
10 sin( f − Nt + (cid:36) ) (cid:3) − sec f (cid:2)
15 sin(5 f + Nt − (cid:36) ) +
15 sin(5 f − Nt + (cid:36) ) + f + Nt − (cid:36) ) + f − Nt + (cid:36) ) +
75 sin(3 f + Nt − (cid:36) ) +
75 sin(3 f − Nt + (cid:36) ) +
35 sin(3 f + Nt − (cid:36) ) +
35 sin(3 f − Nt + (cid:36) ) +
150 sin( f + Nt − (cid:36) ) −
150 sin( f − Nt + (cid:36) ) +
70 sin( f + Nt − (cid:36) ) +
10 sin( f − Nt + (cid:36) ) (cid:3)(cid:111) , ∆ δ(cid:15) = − η αµ K r µ A µ (cid:110) sec f (cid:2) f + Nt − (cid:36) ) + f − Nt + (cid:36) ) + f + Nt − (cid:36) ) + f + Nt − (cid:36) ) + f − Nt + (cid:36) ) +
15 sin( f + Nt − (cid:36) ) −
15 sin( f − Nt + (cid:36) ) (cid:3) − sec f (cid:2) f + Nt − (cid:36) ) + f − Nt + (cid:36) ) + f + Nt − (cid:36) ) + f + Nt − (cid:36) ) + f − Nt + (cid:36) ) +
15 sin( f + Nt − (cid:36) ) −
15 sin( f − Nt + (cid:36) ) (cid:3)(cid:111) .Comparison between the magnitude of the change in hyperbolic elements in Eqs. (B35) and (B21) reveals C tide C J ∝ µ A K µ J r R η ,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 e ff ects are at least 5 ordersof magnitude smaller than the one of the oblateness. However, they can become important when the light ray passes throughthe atmosphere of a satellite orbiting a central massive planet. In such a configuration the tidal e ff ects are expected to be moreimportant. For instance, let us consider the case of a light ray crossing Titan’s atmosphere; substitution of numerical values (wetook the Titan’s J value from [36]) reveals that the changes due to Saturn’s tides are expected to be the same order of magnitudethan changes due to Titan’s oblateness.All the results in Eqs. (B35) are derived under the assumption that the main planet is made with a perfect fluid mediumhaving a non-unity index of refraction. However, such a perfect fluid model might be somewhat inaccurate in the context oftidal dynamics. The main reason is linked to the fact that perfect fluid model does not admit a mechanism to dissipate energy1and the fluid’s response to a tidal stress is purely elastic. In the context of celestial mechanics, energy dissipation mechanismsare of prime importance in a large number of applications [28]. In Sec. B 4, we compute the first order e ff ect on the hyperbolicelements caused by the non-elastic response of the atmosphere following a tidal stress. Referring to Eq. (B38), we see that thenon-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 [37], so the ratio is of theorder 2 × − . So, if the elastic contribution is important in some circumstances, the viscoelastic response of the atmosphererepresents a perturbation of a few percent of the amplitude of the elastic one. In the case of Saturn and Jupiter, we can use thefact that w τ ≈ (2 Q ) − where Q is the specific dissipation function [28]. The numerical substitution is done by taking value of k / Q (with k the gravitational Love number of degree 2) from [38] together with the value of k provided by [39] for Jupiterand the one provided by [40] for Saturn. We find the ratios to be respectively of the order of 8 × − and 7 × − , which is wellbeyond the elastic response. F. Horizontal gradients
The last study is a direct application for GB observations operating for the realization of IERS reference systems e.g. SLR,VLBI, or the Global Positioning System (GPS). In [5] or in [6] it is shown that the group delay due to the horizontal gradientsin the Earth atmosphere can overpass the centimetric level while observational techniques like SLR or VLBI currently operatewith sub-centimeter accuracy measurements. Therefore, not considering horizontal gradients may lead to systematic errors inthe estimation of the station coordinates, which can have repercussions on the determination of the International Terrestrial Ref-erence Frame (ITRF). Consequently, delays due to horizontal gradients must be taken into account, especially at low elevationswhere the e ff ect is maximum.The computations are carried out in Sec. B 5, and the equations describing the change in the hyperbolic elements followinghorizontal gradients are given in Eqs. (B46). 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. Thissituation is encountered when the azimuth of the transmitter equals the azimuth of the observer, that is to say when the sourceis at the meridian of the receiving site. In such a case, the propagation plane is aligned with the meridian and intersects Earth’sequator for ι = π/ λ = Ω . For this type of inclination, it is seen from Eqs. (B46) that only a West-East horizontal variationof the refractivity can change the inclination or the longitude of the node. Conversely, only a North-South horizontal gradient isexpected to change the other hyperbolic elements.We copy here the results for the non-hyperbolic contribution to the light time δ t ( f ) = sign( f − f ) Kn NS η c (cid:104) sec f − sec f + | cos f | − | cos f | − f − f ) tan f + f − f ) tan f (cid:105) . (62)For applications, the true anomaly is not a common way of expressing the results. Instead, the colatitude (referred from thedirection of the North pole) may be preferred. In order to pass from true anomalies to colatitudes, we make use of hyperbolicrelations defined in appendix A, in particular Eqs. (A10b) and (A11). Because, this transformation is only a matter of rewritingthe boundary conditions of Eq. (62), we can safely approximate the photon path assuming straight line between the sourceand the receiver. In this case, the hyperbolic relations in appendix A can be further simplified taking the limit α →
0, whichcorresponds to light propagation in vacuum.As an example, we give the transformation rule which allows us to pass from true anomaly to colatitude. For the photon path1 → (cid:48) , which is depicted on the left panel of Fig. 4, we get f = ¯ ϕ − ( ω (cid:48) − ω ), ¯ ϕ = ϕ − ϕ e . (63a)Here, for a matter of convenience for the next, we have introduced the colatitude ¯ ϕ which is referred from ˆ e ’s direction insteadof North’s (see Fig. 4). The di ff erence ω (cid:48) − ω appearing on the right-hand side of Eq. (63a), is function of the location of thereceiving antenna ( ¯ ϕ (cid:48) ) and is also function of the ratio between the magnitude of the separation vectors of the source ( h ) andthe receiver ( R ) as q ≡ h / R . It is given by the following expression ω (cid:48) − ω = arctan − cos ¯ ϕ (cid:48) − sign( ¯ ϕ (cid:48) ) (cid:112) q − + sin ¯ ϕ (cid:48) . (63b)Finally, insertion of Eqs. (63) into (62) allows one to infer the atmospheric time-delay along the photon path [ δ t ( ¯ ϕ )] and for thefollowing boundary conditions ( q , ¯ ϕ , ¯ ϕ ).In order to emphasize the behavior of δ t ( ¯ ϕ ; q , ¯ ϕ , ¯ ϕ ), we have represented in Fig. 4, the evolution of the North-South horizontalgradient contribution to the atmospheric time-delay. The computation has been carried out assuming three stations (labeled 2,2 • •• • North ˆ e ˆ e (cid:48) ˆ e (cid:48)(cid:48) ϕ e ωω (cid:48)(cid:48) ω (cid:48) f (cid:48) f (cid:48)(cid:48) f (cid:48)(cid:48) f (cid:48) f E2 12 (cid:48) (cid:48)(cid:48) L i n e o f n o d e s δ t ( ϕ ) Rn NS η c ¯ ϕ [rad] q = . , ¯ ϕ = q = . , ¯ ϕ = q = . , ¯ ϕ (cid:48) = − . q = . , ¯ ϕ (cid:48) = − . q = . , ¯ ϕ (cid:48)(cid:48) = − . q = . , ¯ ϕ (cid:48)(cid:48) = − . FIG. 4. Evolution of the North-South horizontal gradient contribution to the atmospheric time-delay for observations at meridian. The leftpanel is a sketch representing a transmitting source (labeled 1) being ground tracked simultaneously by three stations on Earth (labeled 2, 2 (cid:48) ,and 2 (cid:48)(cid:48) ). The source is passing in the meridian of the three stations, i.e. ι = π/ λ = Ω , thus the direction of the line of nodes is theintersection between the propagation plane of photons and the Earth’s equator. The right panel is a graph of Eq. (62) [making use of (63)], forthe three di ff erent paths which are depicted on the left panel . The plain and dotted curves are computed for q = . q = . ϕ = ϕ (cid:48) = − . ϕ (cid:48)(cid:48) = − . ϕ = ¯ ϕ (cid:48) = ¯ ϕ (cid:48)(cid:48) = − arccos[1 / q ]. (cid:48) , and 2 (cid:48)(cid:48) ) which are simultaneously ground tracking a unique source (labeled 1). We have assumed that the stations are locatedat di ff erent colatitudes along the same meridian.For a given ratio q , it is seen that the largest e ff ect is reached for station 2, which is ground tracking the source at the minimalelevation. This is due to three cumulative e ff ects. Firstly, the projection of the horizontal gradient along the direction of theray is maximum for the path 1 →
2. Secondly, the geometrical length inside the atmosphere is the most important for the path1 →
2. Finally, the horizontal contribution in ( ϕ − ϕ ) n NS into the expression of the index of refraction [see Eq. (B39)] is themost important for ϕ = ϕ and for the path 1 → ϕ − ϕ . Consequently, for thispath, the phase velocity is minimum and consequently the time-delay is maximum. VI. COMPARISON WITH NUMERICAL INTEGRATION
The validity of the approximated solutions which are derived in previous section can be assessed by a direct comparison withthe output of a numerical integration of the light path.As discussed so far, the solutions can indi ff erently be applied to AO experiments or astro-geodetic GB observations. Thedi ff erences remain in the role of α , which is either an output or an input, and also in the geometry through the values of thetrue anomaly at the transmission and at the reception (or at the entrance and the exit of the region of refractive influence).However, we have seen that the method of computation for AO experiments may di ff er for GB observations, especially at veryhigh elevations for which an expansion in 1 / e might not be accurate enough.In this chapter, we test the validity of the derived solutions in the context of AO experiments. We focus on Eqs. (B21)which solve, in a non-ubiquitous way, the problem of determining analytically the light path in the presence of the atmosphere’soblateness. As discussed in Chap. I, this problem has never been solved analytically in a complete and satisfactory way.In order to assess the validity of Eqs. (B21), we have performed a numerical integration of Eqs. (1), (4), and (5) across therefractive profile given in Eqs. (B15) and (B17) for l =
2. The numerical integration has been carried out in double precision,with a numerical error tolerance of 10 − , for di ff erent values of α and J . [Actually, instead of α , we work in the followingwith the dimensionless parameter αµ/ R , which represents the hyperbolic refractivity evaluated at the level of the radius of thecentral planet, N ( R ) = n ( R ) − η ]. The tested values of the refractivity and J , range from 10 − to 10 − , and from 10 − to10 − respectively. For each numerical integration, we have compared the total change in the hyperbolic elements, between thetransmission and the reception, with the analytical predictions given in Eqs. (B21). Let us emphasize that δ s , δ t , and δψ cannotbe determined easily from the numerical integration of Eqs. (1), (4), and (5), thus, instead of working with the non-hyperboliccontributions alone, we are considering the total change in s , t , and ψ , which is easily inferred from results of the numericalintegration. From an analytical point of view, the total change in s , t , and ψ , is simply given by the sum of hyperbolic [see3 err( p ) err( e )err( ι ) err( Ω )err( ω ) err( ψ )err( s ) err( t ) J J − − − − − FIG. 5. Evolution of the relative error on the change in elements C = ( p , e , ι, Ω , ω, s , t , ψ ), for di ff erent values of the refractivity at the surface(colored curves) and for di ff erent values of J ( x -axis). The relative error on the change in elements is computed as err( C ) = | ∆ C num − ∆ C ana | / | ∆ C num | , with ∆ C being the total change in the value of the elements between the transmission and the reception of the signal. Thesubscripts “num” refers to the numerical predictions while the subscript “ana” refers to the analytical solutions [cf. Eqs. (B21)]. The changein s , t , and ψ is the total change including the hyperbolic and the non-hyperbolic contributions. Eqs. (A30), (A31), and (A36)] and non-hyperbolic contributions [see Eqs. (B21f)–(B21h)]. In Fig. 5, we show the evolution ofthe relative error on the total change in the hyperbolic elements and s , t , and ψ .Firstly, it is seen that the solutions, for low refractivity [ N ( R ) below 10 − ] and for small values of J , seem to su ff er from alack of accuracy (c.f. errors on p , e , ι , Ω , and ω in Fig. 5). However, this loss of accuracy is only due to numerical noise, sincefor small values of J , the perturbing e ff ect is so small that the hyperbolic elements tend to be constants. The di ff erence betweentheir values at the reception and at the transmission remains at the level of the numerical noise. For instance, for N ( R ) = − ,we have p ( f ) ∼ , moreover the analytical solution predicts ∆ p ana ∼ − for J = − , which represents a change in ∆ p ana / p ( f ) ∼ − , i.e. one order of magnitude beyond the numerical double precision.Secondly, we also notice an important feature linked to the accuracy of the solutions. Indeed, by changing the order ofmagnitude of the refractivity we also change the order of magnitude of the relative error on the total change in the value ofall parameters. This feature reveals the approximation at leading order in 1 / e , which has been assumed in order to simplifythe integration of the perturbation equations. In fact, while changing the order of magnitude of the refractivity value, whichcorresponds to a change in α , we also modify 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 . Forvery high refractivity (e.g. 10 − which corresponds here to a value of eccentricity around 10), it is seen (cf. purple curves inFig. 5) that the change in the hyperbolic elements, as given by the analytical expressions, is accurate at the level of ∼ ∼ − ), the analytical solutions in Eqs. (B21a)–(B21e) are found to be accurate at thelevel of ∼ . J in thelog-log plot. This behavior is due to the fact that the hyperbolic contribution is included within the analytical computation of s , t ,and ψ . Indeed, for hyperbolic elements, as discussed above, when J tends to be null, the computation of ∆ C num generates a lotof 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 remainswhen the non-hyperbolic one vanishes. In order to prove that the non-hyperbolic evolution of s , t , and ψ is well described by4Eqs. (B21f)–(B21h), one can compute the relative error considering only the hyperbolic contribution in ∆ C ana . The computationreveals that, the relative error computed without the non-hyperbolic contributions is between one and two orders of magnitudelarger 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 refractivebending can be assessed by looking at Fig. 5. For very high refractivity (e.g. 10 − ), the maximum relative error is reached for themaximum value of the oblateness ( J = − ). For instance, considering the hyperbolic and the non-hyperbolic contributions,we are able to provide analytical solutions which are accurate at the level of ∼ .
01% for both s and t , and at the level of ∼ ψ . For a typical value of the refractivity ( ∼ − ), and for a value of Jupiter or Saturn’s J ( ∼ − ), we might expect errorsat the level of ∼ − % for both s and t , and at the level of ∼ . ψ . For the same values of the refractivity and J , if weonly consider the hyperbolic contribution, the relative error grows to ∼ − % for both s and t , and to 0 .
1% for ψ . VII. CONCLUSION AND PERSPECTIVES
In this paper we have demonstrated, in Sec. III, that the equations of geometrical optics possess an exact solution, referred toas the reference solution, when one assumes the hydrostatic equilibrium together with the constancy of the variation of the indexof refraction with the Newtonian potential. Not surprisingly, since the problem reduces to a central field problem, the solutionis found to be a conic-section similar to the solution of the well-known two-body problem in celestial mechanics. For thecases of interest, the reference solution is an hyperbola which is completely described with some constants of integration calledhyperbolic elements. These elements describe the geometry of the light beam, namely the shape and the spatial orientation of thehyperbola. We have shown that the hyperbolic elements are related to the refractive response of the medium to the gravitationalpotential which is parametrized by the derivative of the index of refraction with respect to the potential.The reference solution is found to be applicable in situation involving light propagation in planetary atmospheres, and thus,to be well suited for future applications to GB observations or AO experiments. In Sec. III D, we provided short indications onhow applying it in the context of real observations, but its practical application to real space missions’ data will be discussedin a successive paper. However, the main feature of this solution is not its direct application, but the fact that it provides acomprehensive framework that can be extended to carry out analytical studies in the case where the index of refraction hasgeneric dependencies.Indeed, based on the method of variation of arbitrary constants, in Chap. IV we have converted the equation of geometricaloptics as length rate of change in the hyperbolic elements. These new equations are shown to be the optical counterpart ofthe perturbation equations in celestial mechanics. When no assumption is made on the order of magnitude of the change inthe refractive profile, they are perfectly equivalent to equations of geometrical optics. If they are less compact, they presentthe advantage of providing a comprehensive alternative description to the path of light rays. They describe quantitatively thedeparture from the hyperbolic path due to additional terms besides the hyperbolic contribution into the index of refraction.In the case where the perturbing gradient can be assumed small with respect to the hyperbolic one, we have shown in Sec. IV E,how to approach analytically the solutions of the perturbation equations. This procedure seems to be generic enough to easilyhandle the first order e ff ects due to any kind of perturbation besides the hyperbolic contribution.To highlight the capabilities of this formalism, in Chap. V we have analyzed di ff erent sources of perturbing gradients. Forinstance, we have studied the e ff ects on the light path due to, the non-linear dependencies with the Newtonian potential, thecentrifugal potential and the light dragging e ff ect due to the rigid rotation of the atmosphere, the non-spherical gravitationalpotential of the planet, the tidal atmospheric bulge raised by the presence of a massive perturbing body, and the horizontalgradients of refractivity inside the atmosphere. These examples of utilization are just a non-exhaustive list showing the possibilityo ff ered by the perturbation equations formalism applied to geometrical optics.Finally, in Sec. VI, in order to assess the validity of, i) the perturbation equations, and ii) the method of integration, wehave compared the analytical solution derived in the context of a quadrupolar axisymmetry in the gravitational potential of theplanet with its numerical counterpart determined from a numerical integration of the equations of geometrical optics. We carriedout several numerical integrations for di ff erent values of the refractivity and the J parameter. For each one of them, we havecompared the changes in the hyperbolic elements as provided by the numerical integration to those obtained from the analyticalsolution. By doing this, we have been able to assess the accuracy of the analytical solution. We have shown that for standardrefractivity values ( ∼ − ), the relative error in the total change in the hyperbolic elements is of the order of 0 .
01% and isindependent 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 ( ∼ − ), the relative error isof the order of 10 − % for both the length and the time, and 0 . ff ering in this way the possibility of analyticallystudying other symmetries involving small refractivity (e.g. for application to occultations by the Io plasma torus [41–43]). Thisis also something which can can be easily achieved from Eqs. (54) by taking the limit α → ffi cient tool for improving existing models of the Earth tropospheric delay. In particular, an accuratemodeling of horizontal gradients is of prime importance for the determination of station coordinates which in turn define thescale and the origin of the ITRF. However, existing models of the tropospheric delay always stick to the spherical symmetryassumption, and only few of them consider horizontal gradients ([5] for VLBI, [6] for SLR). In addition, even for those whichconsider horizontal gradients, the integration across the atmosphere is always done numerically assuming spherical layers madeof constant index of refraction. The perturbation equations could easily tackle these issues by providing e.g. the contribution tothe time-delay due to the atmosphere’s oblateness and the one due to horizontal gradients, as discussed in Secs. V D and V F.In the context of AO experiments, the perturbation equations allow to assess for the first time a clear description of the photonpath, the light time or the refractive bending into an oblate atmosphere. As mentioned in the introduction, this question isusually tackled using numerical ray-tracing techniques [17], and only a few studies [16] aimed at exploring analytically thisproblem. Yet, the main advantage of analytical formulations is twofold. First, they help to acquire a physical understanding ofthe phenomenon, and secondly, they do not require high computation time while analyzing observational data. If the solutionproposed by [16] is a first analytical approach, it fails to provide a clear physical understanding of the photon path or theatmospheric time-delay due to the atmosphere’s oblateness. The reason is due to the fact that it consists in a Taylor seriesexpansion around a local point along the photon path, and in addition it is expressed in terms of Cartesian coordinates whichmight be somehow abstract. Because of this Taylor expansion, the solution cannot be employed for high refractivity and isonly defined in the surroundings of the expansion point. With the perturbation equations, we have solved all these di ffi cultiessince the solution is available all along the light path trajectory and can also be applied for high refractivity. In addition, thesolution provides a simple geometrical interpretation by fitting at any length along the ray an hyperbolic trajectory. Anotheradvantage which is worth mentioning is that the perturbation equations allows to express directly the light time as well as therefractive bending expressions which are needed for computing range and range-rate observables. Finally, let us mention thatif analytical studies are not as precise as purely numerical ones, the formalism of the perturbation equations has been shown tobe really accurate (errors of one part in 10 and 10 between the numerical and the analytical predictions on the light time andthe refractive bending, in the presence of an oblate atmosphere characterized by J = − and N ( R ) = − , respectively). Inaddition, we have demonstrated that the formalism presented in this paper o ff ers a large flexibility since it is able to handle a widerange of perturbations including light dragging e ff ects caused by atmosphere’s winds. Thus, a complete and purely analyticalmethod, only based on analytical solutions derived in Chap. V, could be developed in order to process real AO data in the contextof past, current, and future space missions. Our next opportunity, which actually motivated this work, is ESA’s L-class missionJUICE (JUpiter Icy moons Explorer) [44, 45]; here a radioscience experiment named 3GM (Gravity and Geophysics of Jupiterand the Galilean Moons) will take advantage of a careful Jupiter system tour design [46], o ff ering both frequent Jupiter moon’sflybys [47, 48] and radio occultation opportunities by Jupiter’s oblate atmosphere. ACKNOWLEDGMENTS
The Authors are grateful to the University of Bologna and to Italian Space Agency (ASI) for financial support through theAgreement 2013-056-RO in the context of ESA’s JUICE mission. A.B. is grateful to F. Mignard from Observatoire de la Cˆoted’Azur for interesting discussions about atmospheric and astronomic refraction. A.B. is also grateful to A. Hees from ParisObservatory for valuable comments about a preliminary version of this manuscript. [1] G. Petit, B. Luzum, and et al., IERS Technical Note (2010).[2] J. W. Marini and C. W. Murray, NASA Technical Memorandum NASA-TM-X-70555 (1973).[3] V. B. Mendes and E. C. Pavlis, Geophys. Res. Lett. , L14602 (2004).[4] V. B. Mendes, G. Prates, E. C. Pavlis, D. E. Pavlis, and R. B. Langley, Geophys. Res. Lett. , 1414 (2002).[5] G. Chen and T. A. Herring, J. Geophys. Res. , 20 (1997).[6] G. C. Hulley and E. C. Pavlis, Journal of Geophysical Research: Solid Earth (2007), 10.1029 / , 1243 (1965).[8] G. Fjeldbo and V. R. Eshleman, J. Geophys. Res. , 3217 (1965).[9] G. Fjeldbo and V. R. Eshleman, Planet. Space Sci. , 1035 (1968).[10] R. A. Phinney and D. L. Anderson, J. Geophys. Res. , 1819 (1968).[11] G. Fjeldbo, A. J. Kliore, and V. R. Eshleman, AJ , 123 (1971). [12] P. J. Schinder, F. M. Flasar, E. A. Marouf, R. G. French, C. A. McGhee, A. J. Kliore, N. J. Rappaport, E. Barbinis, D. Fleischman, andA. Anabtawi, Icarus , 460 (2011).[13] P. J. Schinder, F. M. Flasar, E. A. Marouf, R. G. French, C. A. McGhee, A. J. Kliore, N. J. Rappaport, E. Barbinis, D. Fleischman, andA. Anabtawi, Icarus , 1020 (2012).[14] E. R. Kursinski, G. A. Hajj, J. T. Schofield, R. P. Linfield, and K. R. Hardy, J. Geophys. Res. , 23429 (1997).[15] A. K. Steiner, G. Kirchengast, and H. P. Ladreiter, Annales Geophysicae , 122 (1999).[16] G. F. Lindal, AJ , 967 (1992).[17] P. J. Schinder, F. M. Flasar, E. A. Marouf, R. G. French, A. Anabtawi, E. Barbinis, and A. J. Kliore, Radio Science , 712 (2015).[18] L. D. Landau and E. M. Lifshitz, The classical theory of fields (Ellipse, 1975).[19] M. Born and E. Wolf,
Principles of Optics (Cambridge University Press, 1999).[20] J. L. Synge,
Relativity: The General Theory (North-Holland Publ. Co., Amsterdam, 1960) russian translation: IL (Foreign Literature),Moscow, 1963.[21] P. Teyssandier, C. Le Poncin-Lafitte, and B. Linet, in
Lasers, Clocks and Drag-Free Control: Exploration of Relativistic Gravity in Space ,Astrophysics and Space Science Library, Vol. 349, edited by H. Dittus, C. Lammerzahl, & S. G. Turyshev (2008) p. 153, arXiv:0711.0034[gr-qc].[22] A. Hees, S. Bertone, and C. Le Poncin-Lafitte, Phys. Rev. D , 064045 (2014), arXiv:1401.7622 [gr-qc].[23] U. Leonhardt and P. Piwnicki, Phys. Rev. A , 4301 (1999), physics / , 1 (2004).[25] E. Poisson and C. M. Will, Gravity (Cambridge University Press, 2014).[26] H. Goldstein,
Addison-Wesley World Student Series, Reading, Mass.: Addison-Wesley, 1950 (1950).[27] D. Brouwer and G. M. Clemence,
Mechanics of Composite Materials (1961).[28] C. D. Murray and S. F. Dermott,
Solar System Dynamics (Cambridge University Press, 2000).[29] W. W. Vaughan,
American Institute of Aeronautics and Astronautics , American Institute of Aeronautics and Astronautics
AIAA G-003C-2010 (2010).[30] J. A. Burns, American Journal of Physics , 944 (1976).[31] A. Morbidelli, Modern celestial mechanics : aspects of solar system dynamics (2002).[32] J. M. Aparicio and S. Laroche, Journal of Geophysical Research: Atmospheres (2011), 10.1029 / (1818).[34] N. N. Rozanov and G. B. Sochilin, Optics and Spectroscopy , 441 (2005).[35] Y. Kaspi, E. Galanti, W. B. Hubbard, D. J. Stevenson, S. J. Bolton, L. Iess, T. Guillot, J. Bloxham, J. E. P. Connerney, H. Cao, D. Durante,W. M. Folkner, R. Helled, A. P. Ingersoll, S. M. Levin, J. I. Lunine, Y. Miguel, B. Militzer, M. Parisi, and S. M. Wahl, Nature , 223(2018).[36] L. Iess, R. A. Jacobson, M. Ducci, D. J. Stevenson, J. I. Lunine, J. W. Armstrong, S. W. Asmar, P. Racioppa, N. J. Rappaport, andP. Tortora, Science , 457 (2012).[37] W. M. Folkner, J. G. Williams, D. H. Boggs, R. S. Park, and P. Kuchynka, Interplanetary Network Progress Report , 1 (2014).[38] V. Lainey, J.-E. Arlot, ¨O. Karatekin, and T. van Hoolst, Nature , 957 (2009).[39] L. Iess, W. M. Folkner, D. Durante, M. Parisi, Y. Kaspi, E. Galanti, T. Guillot, W. B. Hubbard, D. J. Stevenson, J. D. Anderson, D. R.Buccino, L. G. Casajus, A. Milani, R. Park, P. Racioppa, D. Serra, P. Tortora, M. Zannoni, H. Cao, R. Helled, J. I. Lunine, Y. Miguel,B. Militzer, S. Wahl, J. E. P. Connerney, S. M. Levin, and S. J. Bolton, Nature , 220 (2018).[40] V. Lainey, R. A. Jacobson, R. Tajeddine, N. J. Cooper, C. Murray, V. Robert, G. Tobie, T. Guillot, S. Mathis, F. Remus, J. Desmars, J.-E.Arlot, J.-P. De Cuyper, V. Dehant, D. Pascu, W. Thuillot, C. Le Poncin-Lafitte, and J.-P. Zahn, Icarus , 286 (2017), arXiv:1510.05870[astro-ph.EP].[41] F. Bagenal and J. D. Sullivan, Journal of Geophysical Research , 8447 (1981).[42] P. H. Phipps and P. Withers, Journal of Geophysical Research (Space Physics) , 1731 (2017), arXiv:1701.07435 [astro-ph.EP].[43] P. H. Phipps, P. Withers, D. R. Buccino, and Y.-M. Yang, Journal of Geophysical Research: Space Physics , 6207 (2018),https: // agupubs.onlinelibrary.wiley.com / doi / pdf / / , 1(2013), cited By 149.[45] O. Witasse, Proceedings of the International Astronautical Congress, IAC (2016).[46] A. Boutonnet and J. Schoenmaekers, Advances in the Astronautical Sciences , 1465 (2016).[47] G. Lari, Celestial Mechanics and Dynamical Astronomy , 50 (2018).[48] D. Dirkx, L. Gurvits, V. Lainey, G. Lari, A. Milani, G. Cim, T. Bocanegra-Bahamon, and P. Visser, Planetary and Space Science , 14(2017).[49] V. R. Eshleman, Planet. Space Sci. , 1521 (1973).[50] W. M. Kaula, Theory of satellite geodesy. Applications of satellites to geodesy (Waltham, Mass.: Blaisdell, 1966).[51] T. Hartmann, M. H. So ff el, and T. Kioustelidis, Celestial Mechanics and Dynamical Astronomy , 139 (1994). Appendix A: Hyperbolic path
As discussed in Chap. III, when the index of refraction is only function of the height of the ray (spherically symmetricassumption), all the light path is contained in the propagation plane. Using this fact, we derive, in this chapter, a completesolution to Eqs. (1), (4), and (5). We start by determining the shape of the light trajectory (i.e. the light path), then we focus onthe time-delay.
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 thepolar basis. The vectors h and ∇ S are decomposed as h = h ˆ ρ , (A1a) ∇ S = n ˙ h ˆ ρ + nh ˙ θ ˆ τ , (A1b)where a dot denotes the di ff erentiation with respect to s . Moreover, the left-hand side of Eq. (4) is given bydd s (cid:0) ∇ S (cid:1) = ( n ¨ h + ˙ n ˙ h − nh ˙ θ ) ˆ ρ + h dd s ( nh ˙ θ ) ˆ τ . (A2)The gradient of the index of refraction can also be expressed in the same coordinate system. With the help of Eq. (11) and themonopole approximation [cf. Eq. (16)], the right-hand side of Eq. (4) is given by ∇ n = − α g ˆ ρ . (A3)The term g represents the magnitude of the local acceleration (monopole contribution) experienced by the media which makesup the atmosphere of the central planet. It is explicitly given by g = d U d h = µ h . (A4)It is obvious from the presence of α in Eq. (A3), that the light ray experiences the gravitational pull via the media in which theray propagates through.From Eqs. (A2) and (A3), it is seen that the absence of a transverse component (along ˆ τ ) in the gradient of the index ofrefraction implies that nh ˙ θ is a conserved quantity. From Eqs. (12) and (A1) we immediately deduce K = nh ˙ θ . (A5)Inserting this expression into the radial part of Eq. (A2), we find¨ h + ˙ n ˙ hn − K n h = − α gn . (A6)We change the independent variable from s to θ adopting the convention that a prime denotes the di ff erentiation with respect to θ . We also adopt u ≡ / h as a convenient substitute for h and derive a di ff erential equation for it. Making all these substitutionsinto Eq. (A6), we quickly arrive to u (cid:48)(cid:48) + κ u = αµη K , (A7)with κ = − ( αµ ) K (A8)a constant parameter. The general solution to the simple Eq. (A7) is a conic section with origin at the focus. Returning to theoriginal radial variable the spatial solution is h = p + e cos κ f , (A9)8in which e is an arbitrary constant of integration called the eccentricity of the conic. We will adopt the convention that theeccentricity is a positive quantity. The two other parameters, p and f , are defined such that p ≡ κ K αµη , (A10a) f ≡ θ − ω . (A10b)They are respectively known as the semi-latus rectum and the true anomaly . The last parameter, ω , is also a constant ofintegration and its role is better understood when h is minimal, that is to say when θ = ω . It is seen that ω defines the argumentof the closest approach which is the minimal separation between the path of the ray and the center of the reference frame. Thisangle ω is completely determined by the initial conditions of the problem. For instance, if we call ( h , θ ) the components of theposition vector of the transmission point, then from the solution in Eq. (A9), we end up with ω = θ ± κ − arccos (cid:32) p − h eh (cid:33) , (A11)where the sign is taken to be positive when − π/ ≤ θ < π/ h C = p + e . (A12)In the following, we will consider that the signal was emitted or received (or both) from infinity, in a space region wellbeyond 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 ≥ e >
1, which corresponds to hyperbolic trajectories for the lightpath. Therefore, we can introduce a new constant known as the semi-major axisa ≡ − pe − ∆ f ∞ , describ-ing the net change between the two asymptotes’ directions. Based on Eq. (A9), it can be seen that h goes to infinity when cos κ f tends to − / e , which gives rise to two asymptotic solutions for κ f . Let call f ∞ in the negative solution and f ∞ out the positive one. Ifwe introduce ∆ f ∞ = f ∞ out − f ∞ in , we immediately deduce ∆ f ∞ ≡ κ arccos (cid:0) − / e (cid:1) , (A14)Since, we focus on hyperbolic trajectories ( e > ∆ f ∞ < π .In AO experiments, the geometry is such that ∆ f ∞ can be geometrically related to the refractive bending (indeed from Fig. 1,we deduce ∆ f ∞ = π + (cid:15) ) which is itself directly in relation with the changes in the Doppler measurements, at first order (seeEq. (C2) or [8, 9]). Thus, Eq. (A14) 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. (A1b)].We thus need to determine the expressions of the radial and the transverse components; the later depends on the angular velocity.Invoking Eqs. (A9) and (A5) and di ff erentiating them making use of Eqs. (A10), reveals that n ˙ h = e √ η (cid:114) αµ p sin κ f , (A15a) n ˙ θ = √ ηκ (cid:114) αµ p (1 + e cos κ f ) . (A15b)The evolution of the tangent to the ray is provided after inserting those relationships together with Eq. (A9) into Eq. (A1b).
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 fundamentalconstants of the problem, such as K . We saw e.g. that the semi-latus rectum is linked to the conservation of the magnitude of9the impact parameter [cf. Eq. (A10a)]. In this section, we look for additional first integrals of the light path, which then couldlet to express the constant of integration in a fundamental way.The first fundamental constant can be obtained by multiplying both side of Eq. (A6) by n ˙ h . Then, one can see that eachterm is a total di ff erential with respect to s . In addition, one can recognized the squared components of (A1b) in the integratedequation, and we end up with E = ( ∇ S ) − αµ h (cid:18) η + αµ h (cid:19) , (A16)where E is a constant of integration. That constant would be similar to the energy in classical mechanics. Indeed, by computingthe square of ∇ S from Eqs. (A15) and (A1b), we end up with( ∇ S ) = αµ (cid:32) η h − η a + αµ h (cid:33) , (A17)where we have used Eq. (A13). Then, insertion of Eq. (A17) into (A16), leads to E = − αµη a , (A18)which possesses the same form as the energy in the Kepler problem. From the definition of a , the right-hand side of the equationis 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. (A16) can be put under a more fundamental form. Indeed, by making use of thedefinition of n [cf. Eq. (17)], we deduce E = (cid:16) ∇ S · ∇ S − γ n (cid:17) , (A19)where γ ≡ (cid:32) − η n (cid:33) (A20)is known as the Fresnel coe ffi cient [ η has been previously defined in Eq. (17)]. This expression of E is more fundamental thanEq. (A16) since it does not require a prior knowledge of the index of refraction. Indeed, the constancy of E can be inferred usingonly Eqs. (1) and (4).From Eq. (A18), it is seen that the semi-major axis is linked to the conservation of the dimensionless energy parameter, E .At the same time, the constancy of the eccentricity is assured, because of the definition (A13). So, in that sense the eccentricityis linked to both the conservation of energy and the conservation of the magnitude of the impact parameter. Eqs. (A18), (A13),and (A10a) can be used to infer the eccentricity e ≡ (cid:115) + κ EK ( αµη ) . (A21)In addition to E , and K , we can determine an other fundamental constant of the problem. This last one, is equivalent to theeccentricity vector of celestial mechanics and comes from the very peculiar form of the index of refraction, which is linear withthe gravitational potential making it evolving as 1 / h . However, because κ is di ff erent from unity, the eccentricity vector does notshow up as easily as in celestial mechanics. It is given by the following relationship e = pK ∇ S × K − ˆ ρ + e A , (A22)where we have introduced the vector A which possesses the following components in the polar basis A = (cos f − cos κ f ) ˆ ρ + ( κ sin κ f − sin f ) ˆ τ . (A23)The vector A is expressed in terms of fundamental quantities after determining f from Eq. (A9), and expressing ˆ ρ and ˆ τ interms of h and K . The eccentricity in Eq. (A22), must be expressed in terms of E and K thanks to Eq. (A21).The constancy of e can be demonstrated by substituting Eqs. (A15) into (A1b), then inserting the result together with Eq. (A5)into (A22), which finally gives with Eq. (15) e = e (ˆ x cos ω + ˆ y sin ω ). (A24)0Since e and ω are two constants of integration, e is conserved during the propagation of light. Moreover, it is seen that theconstancy of ω is linked to the conservation of the eccentricity vector. It always points towards the closest approach and itsmagnitude 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 whichremains orthogonal to the impact parameter vector when the refractive profile is purely radial. In addition, we found that theshape of the trajectory is described by an hyperbola when the refractive profile is given by Eq. (17). That shape is parametrizedthanks 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. (A10a), (A21), and (A24)] in terms of fundamental quantities of theproblem, such that K , E , and e .
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 termsof constant parameters. In order to completely solve the problem in the plane, we need to determine a single location along thelight path given a geometrical length or a light time.The Eq. (A15b) 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. (A15b) as s ( f ) − S = κ √ η (cid:115) p αµ (cid:90) f n ( f (cid:48) )d f (cid:48) (1 + e cos κ f (cid:48) ) , (A25)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. (A25) into t ( f ) − T = κ c √ η (cid:115) p αµ (cid:90) f n ( f (cid:48) ) d f (cid:48) (1 + e cos κ f (cid:48) ) , (A26)where T is an arbitrary constant of integration representing the light time till the closest approach .Eqs. (A25) and (A26) are the analogous of the Kepler equation of celestial mechanics. They can be exactly solved byintroducing a new variable, F , known as the hyperbolic anomaly . The change from the true to the hyperbolic anomaly, andconversely, is given by cosh F = e + cos κ f + e cos κ f , sinh F = √ e − κ f + e cos κ f , (A27a)cos κ f = e − cosh Fe cosh F −
1, sin κ f = √ e − Fe cosh F − tan (cid:32) κ f (cid:33) = (cid:114) e + e − (cid:18) F (cid:19) . (A27c)Eqs. (A27) allow us to introduce the following relations which relate the di ff erentials of f and F d f d F = κ − √ e − e cosh F −
1, d F d f = κ √ e − + e cos κ f . (A28)The hyperbolic anomaly can be employed instead of the true anomaly as a convenient substitute. Making all these substitutionswe quickly arrive to h = − a ( e cosh F − n ˙ h = e √ η (cid:114) αµ − a sinh Fe cosh F − n ˙ F = √ η (cid:114) αµ − a e cosh F − Attention must be paid when one want to determine the hyperbolic from the true anomaly using (A27c). Indeed, the inverse function of the hyperbolic tangentis not well defined when the argument is close to ±
1, which might generates numerical noise. To avoid this issue, the inverse function of the hyperbolic sine[cf. Eq. (A27a)] can be preferred. F , we find the following analogousto the well-known Kepler equation √ η (cid:114) αµ − a ( s − S ) = η e sinh F − η F − αµ a F , (A30)in which the left-hand side is the geometrical length mean anomaly .Similarly, expressing Eq. (A29c) in term of the light time instead of length path with the help of Eq. (5), let one to directlyintegrate Eq. (A26), such that c √ η (cid:114) αµ − a ( t − T ) = η e sinh F − η F − η αµ a F + √ e − αµ ) a arctan (cid:114) e + e − (cid:18) F (cid:19) , (A31)in which the left-hand side is the light time mean anomaly From these two Kepler-like equations, it is seen that the light time [cf. Eqs. (A31)] is not just the geometrical length [cf.Eqs. (A30)] multiplied by the inverse of the speed of light in a vacuum. Indeed, it also involves some additional terms whichaccount for the apparent change in the speed of propagation of the electromagnetic waves. Consequently, the geometrical lengthmean anomaly considers a ray traveling at the speed of light along the bent trajectory while the time mean anomaly is related tothe phase path and takes into account the apparent change in the speed of the waves along the actual path of photons.
4. Argument of the refractive bending
In this section, we take advantage of the introduction of the hyperbolic anomaly to define a very important geometricalrelationship between the true anomaly and the argument of the refractive bending, ψ . As seen in Fig. 1, ψ is the angle betweenthe orthogonal direction to the line of nodes and the tangent to the light ray.In the context of AO experiments, this angle plays an important role. Indeed, from Fig. 1, it might be seen that the net changein ψ between the transmitter (labeled 1) and the receiver (labeled 2) is exactly the refractive bending, (cid:15) ≡ ψ − ψ . Moreover, asshown in Eq. (C2) and in [11] or [49], the refractive bending can be use to infer the Doppler frequency shift which is the mainobservable for AO experiments.The di ff erential relationship between ψ and f can be inferred from Fig. 1 which allows us deduce the following relationship ψ = θ + φ − π φ has been previously defined in Eq. (13) as the angle between ˆ ρ and s . Hence, from Eq. (A1b) weimmediately deduce cos φ = ˙ h , sin φ = h ˙ θ . (A33)Then, di ff erentiating Eqs. (13) and (A32), and making use of Eqs. (A33), we obtain the following relationd ψ = − hn d n d h d f , (A34)which can be written as ψ ( f ) − ω = αµ p (cid:90) f + e cos κ f (cid:48) n ( f (cid:48) ) d f (cid:48) . (A35)Making use of the relationships in Eqs. (A29), we change the variable of integration from the true to the hyperbolic anomaly,and after some algebra we integrate the previous expression as ψ − ω = f ( F ) − e κ (cid:114) e − − ¯ e arctan (cid:114) + ¯ e − ¯ e tanh (cid:18) F (cid:19) , (A36)in which we have introduced ¯ e = κ + ( κ − e − e κ . (A37)Since Eq. (A27c) provides a unique relationship between f and F , Eq. (A36) can be used to find the argument of the bendingangle knowing the true anomaly.2 Appendix B: Mathematical details for applications
In this chapter, we derive the change in hyperbolic elements for di ff erent contributions in the index of refraction.
1. Non-linearity with U We first simplify the problem by assuming that the generalized potential reduces to the spherical part of the Newtonianpotential (
Φ = U ). Consequently, Eq. (61) and the perturbing gradient now read n = n + + ∞ (cid:88) k = α ( k ) k ! (cid:18) µ h (cid:19) k , (B1a) f pert = − h + ∞ (cid:88) k = α ( k ) ( k − (cid:18) µ h (cid:19) k ˆ ρ . (B1b)We can directly see that the non-null components of the perturbing gradient are the radial and the change in the index of refraction R = − h + ∞ (cid:88) k = α ( k ) ( k − (cid:18) µ h (cid:19) k , (B2a) N = np + ∞ (cid:88) k = α ( k ) k ! (cid:18) µ h (cid:19) k . (B2b)Eqs. (B2) can then be inserted into Eqs. (55) and the integrand is expanded in power of 1 / e (see Sec. IV E). The change in thehyperbolic elements is given for any k (with k ≥
2) by the following expressions ∆ e ( k ) = − η k − ( k − e k − α ( k ) α k (cid:90) f f cos k − f sin f d f , (B3a) ∆ ω ( k ) = η k − ( k − e k α ( k ) α k (cid:90) f f cos k f d f , (B3b) ∆ s ( k ) = ( k + η k − k ! e k − α ( k ) α k − µ (cid:90) f f cos k − f d f , (B3c) ∆ t ( k ) = ( k + η k − k ! e k − α ( k ) α k − µ c (cid:90) f f cos k − f d f , (B3d) ∆ ψ ( k ) = ∆ ω ( k ) , (B3e)where we have kept the leading order in 1 / e , and where f ≥ f .Considering e.g. the quadratic dependence with the gravitational potential in the index of refraction profile ( k = e ( f ) = η e α (2) α cos f , (B4a) ω ( f ) = η e α (2) α ( f + cos f sin f ), (B4b) δ s ( f ) = e α (2) α µ f , (B4c) δ t ( f ) = η e α (2) α µ c f , (B4d) δψ ( f ) = ω ( f ). (B4e)All those equations are defined within a constant which is obviously given by the initial condition of the problem at the level ofthe transmitter, so Eqs. (B4) are valid for all f ≥ f .3
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 [25, 50] Φ C = w h (cid:104) ( ˆ ρ · ˆ X (cid:48) ) + ( ˆ ρ · ˆ Y (cid:48) ) (cid:105) . (B5)At linear order in the generalized potential, the index of refraction and its gradient are expressed in the frame comoving with thefluid as n = n + α w h (cid:104) ( ˆ ρ · ˆ X (cid:48) ) + ( ˆ ρ · ˆ Y (cid:48) ) (cid:105) , (B6a) f pert = α w h (cid:104) ( ˆ ρ · ˆ X (cid:48) ) ˆ X (cid:48) + ( ˆ ρ · ˆ Y (cid:48) ) ˆ Y (cid:48) (cid:105) . (B6b)Once ˆ X (cid:48) and ˆ Y (cid:48) are projected into the polar basis, we obtain the following components R = α w p (cid:2) + cos 2 ι + ι cos 2( f + ω ) (cid:3) + e cos κ f , (B7a) T = − α w p f + ω )1 + e cos κ f sin ι , (B7b) S = − α w p f + ω )1 + e cos κ f sin 2 ι , (B7c) N = α w p n (cid:2) + cos 2 ι + ι cos 2( f + ω ) (cid:3) (1 + e cos κ f ) . (B7d)Inserting Eqs. (B7) into Eqs. (55), expending them in power of 1 / e and performing the integration, allows one to find, at firstorder in 1 / e , that the centrifugal contribution impacts all the hyperbolic elements, such p ( f ) = − e η α µ w sin ι (cid:0) ω + sin 2 ω sec f sin 3 f (cid:1) sec f , (B8a) e ( f ) = e η α µ w (cid:16) + ι − sin ι (cid:2)
18 cos 2 ω + sin 2 ω (3 sin f + f ) sec f (cid:3)(cid:17) sec f , (B8b) ι ( f ) = − e η α µ w sin 2 ι (cid:0) ω + sin 2 ω sec f sin 3 f (cid:1) sec f , (B8c) Ω ( f ) = − e η α µ w cos ι (cid:0) f + ω cos f + [1 − ω ] sin 3 f (cid:1) sec f , (B8d) ω ( f ) = e η α µ w cos ω (cid:16) ω (cid:2) + cos 2 ι (cid:3) − cos ω (cid:2) + f − ι sin f (cid:3) tan f (cid:17) sec f , (B8e) δ s ( f ) = − e η α µ w (cid:16)(cid:2) + cos 2 ι (cid:3) sec f (cid:2) f + sin 3 f (cid:3) + ι sin 2 ω cos 2 f sec f +
12 sin ι cos 2 ω tan f (cid:3)(cid:17) sec f , (B8f) δ t ( f ) = − e c η α µ w sin ι (cid:16) ω sec f (cid:2) + f (cid:3) + ω tan f (cid:17) , (B8g) δψ ( f ) = e η α µ w (cid:16) sin 2 ω sin ι sec f − (cid:2) + cos 2 ι + ω sin ι (cid:3) tan f (cid:17) . (B8h)All these solutions are defined within a constant which is given by the initial condition of the problem at the level of thetransmitter, so Eqs. (B8) 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 thedragging e ff ect due to the moving medium. To quantify it, let us express Eq. (58b) as f drag = whc (cid:104) ( ˆ ρ · ˆ X ) ˆ Y − ( ˆ ρ · ˆ Y ) ˆ X (cid:105) ( ∇ S · ∇ n ) + whc (cid:104) ( ˆ ρ · ˆ Y )( ∇ S · ˆ X ) − ( ˆ ρ · ˆ X )( ∇ S · ˆ Y ) (cid:105) ∇ n + γ n wc (cid:104) ( ∇ S · ˆ Y ) ˆ X − ( ∇ S · ˆ X ) ˆ Y (cid:105) , (B9)4in which, we have assumed a steady rotating atmosphere such that v = w × h . If we express the dragging contribution in thepolar basis, we obtain the following non-null components R = √ ηκ wc (cid:114) αµ p ( n γ + N ) cos ι (1 + e cos κ f ), (B10a) T = − e √ η wc (cid:114) αµ p ( n γ + N ) cos ι sin κ f , (B10b) S = − √ ηκ wc (cid:114) αµ p sin ι (cid:104) n γ (1 + e cos κ f ) sin( f + ω ) − e κ ( n γ + N ) cos( f + ω ) sin κ f (cid:105) , (B10c)where N is the refractivity defined as N = n − η . We notice the absence of the N -component of the perturbation.The e ff ect of the fluid rotation on the light propagation is assessed by inserting those components into Eqs. (55), by keepingthe first order in 1 / e , and then by performing the integration over the true anomaly p ( f ) = − γ e η ( αµ ) wc cos ι sec f , (B11a) e ( f ) = − γ e αµ wc cos ι sec f , (B11b) ι ( f ) = γ e αµ wc sin ι sin 2 ω (cid:0) tan ω sec f − f (cid:1) , (B11c) Ω ( f ) = − γ e αµ wc sin 2 ω (cid:0) sec f + ω tan f (cid:1) , (B11d) ω ( f ) = γ e αµ wc cos ι cos ω (cid:0) tan ω sec f − f (cid:1) , (B11e) δ s ( f ) = − γ e η ( αµ ) wc cos ι (2 + cos 2 f ) sec f tan f , (B11f) δ t ( f ) = η c δ s ( f ), (B11g) δψ ( f ) = − γ e αµ wc cos ι tan f . (B11h)Once again let us precise that these equations are defined within a constant which is determine with the transmitter’s position.Eqs. (B11) are defined for f ≥ f .
3. Axisymmetric gravitational potential
A convenient way to handle non-spherical gravitational potential is to deal with either the spherical harmonic expansion[28, 50] or the symmetric and trace free (STF) tensors formalism [25, 51]. Using e.g. the later one, it can be shown that theexpansion of the non-spheric part of the gravitational potential is given by U l = − G ( − l l ! I < L > ∂ < L > h − , (B12)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 [51] and [25] for notations and definitionsabout the use of the STF tensors.We remind that the unit vector ˆ Z has been chosen in Sec. V A for being aligned with the rotation axis of the extended body. Wefocus on the special case where the spatial changes in the rotation axis are supposed to be negligible during the light propagationevent. Thus, we consider that the ˆ Z -direction is spatially fixed. In addition, we consider a fluid body in hydrostatic equilibriumrotating at a constant rate. Then, the resulting gravitational field is considered for being axisymmetric and independent of thelongitude [25]. This symmetry imposes the multipole moments STF tensors to be proportional to ˆ Z < L > which is the STF tensorformed with the ˆ Z unit vector. Hence, we have the relation I < L > ≡ − MR l J l ˆ Z < L > , (B13)where M is the total mass of the extended body, R its equatorial radius, and J l are its unitless multipole moments.5Making use of basic properties about STF tensors [25, 51], such ∂ < L > h − ≡ ( − l (2 l − ρ < L > h − ( l + , (B14a)ˆ Z < L > ˆ ρ < L > ≡ l !(2 l − P l ( χ ), (B14b)ˆ Z < L > ˆ ρ < jL > ≡ l !(2 l + (cid:32) d P l + d χ ˆ ρ j − d P l d χ ˆ Z j (cid:33) , (B14c)where χ = ˆ ρ · ˆ Z and P l ( χ ) are the well-known Legendre polynomials of degree l , we can rewrite the perturbing gravitationalpotential in Eq. (B12) in terms of the multipole moments J l . Keeping the linear part in Φ in the development of the indexof refraction [cf. Eq. (61)] with Φ = U + (cid:80) l U l , we determine an expression for n in term of the non-spherical part of thegravitational potential n = n − αµ h + ∞ (cid:88) l = J l (cid:18) Rh (cid:19) l P l ( χ ). (B15)This profile di ff ers 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 ofmagnitude than the equatorial radius. We can now, determine the j th component of the gradient of the index of refraction whichis given by ∂ j δ n = − α + ∞ (cid:88) l = ∂ j U l . (B16)Di ff erentiating Eq. (B12) and using the previous properties about STF tensors [cf. Eq. (B14c)], we deduce f pert = αµ h + ∞ (cid:88) l = J l (cid:18) Rh (cid:19) l (cid:34) d P l + d χ ˆ ρ − d P l d χ ˆ Z (cid:35) . (B17)In the polar basis the components of f pert are R = αµ h + ∞ (cid:88) l = J l (cid:18) Rh (cid:19) l (cid:34) d P l + d χ − ( ˆ Z · ˆ ρ ) d P l d χ (cid:35) , (B18a) T = − αµ h + ∞ (cid:88) l = J l (cid:18) Rh (cid:19) l ( ˆ Z · ˆ τ ) d P l d χ , (B18b) S = − αµ h + ∞ (cid:88) l = J l (cid:18) Rh (cid:19) l ( ˆ Z · ˆ σ ) d P l d χ . (B18c)Usually, for modestly deformed bodies (e.g. the planets in the Solar System), the most important term in the developments (B15)and (B17) is the one proportional to the quadrupolar moment J . This parameter measures the rotational flattening due to theproper rotation of the extended body. From now, we restrict our attention to that single parameter.Using the definition of Legendre polynomials P ( χ ) =
12 (3 χ − P ( χ ) =
12 (5 χ − χ ), (B19b)with Eqs. (19), we can determine the change in the index of refraction as well as the components of the perturbing gradient R = J αµ h R h (cid:2) ι sin ( f + ω ) − (cid:3) , (B20a) T = − J αµ h R h sin ι sin 2( f + ω ), (B20b) S = − J αµ h R h sin 2 ι sin( f + ω ), (B20c) N = − J αµ np R h (cid:2) ι sin ( f + ω ) − (cid:3) . (B20d)6Inserting those components into Eqs. (55) and applying Eq. (56), allows one to find the following estimations which are givenat leading order in 1 / ep ( f ) = η J e R αµ sin ι (cid:2) f + ω ) + cos(3 f + ω ) (cid:3) , (B21a) e ( f ) = η J e R ( αµ ) (cid:16) (cid:2) + ι cos 2 ω (cid:0) + f (cid:1) + ι (cid:3) cos f − sin ι sin 2 ω (cid:2)
30 sin f +
17 sin 3 f + f (cid:3)(cid:17) , (B21b) ι ( f ) = η J e R ( αµ ) sin 2 ι (cid:2) f + ω ) + cos(3 f + ω ) (cid:3) , (B21c) Ω ( f ) = η J e R ( αµ ) cos ι (cid:2) f + ω ) + sin(3 f + ω ) − f (cid:3) , (B21d) ω ( f ) = η J e R ( αµ ) (cid:16) (cid:2) + cos 2 ι (cid:0) − ω (cid:1) + f cos 2( f + ω ) (cid:3) sin f + (cid:0) + ι (cid:1) sin 3 f − cos 2 ι (cid:2)
42 sin 2 ω cos f +
19 sin(3 f + ω ) + f + ω ) (cid:3)(cid:17) , (B21e) δ s ( f ) = η J e R αµ (cid:16)(cid:2) − ι (cid:0) − cos 2 ω (cid:1)(cid:3) sin f + ι (cid:2) sin 2 ω cos f + sin(3 f + ω ) (cid:3)(cid:17) , (B21f) δ t ( f ) = η J e R αµ c (cid:16)(cid:2) − sin ι (cid:0) − ω (cid:1)(cid:3) sin f + sin ι (cid:2) ω cos f + f + ω ) (cid:3)(cid:17) , (B21g) δψ ( f ) = η J e R ( αµ ) (cid:16)(cid:2) + ι (cid:3) sin 3 f + (cid:2) − sin ι (cid:0) − ω (cid:1)(cid:3) sin f + ι (cid:0) + f (cid:1)(cid:2) f + ω ) + sin(3 f + ω ) (cid:3)(cid:17) . (B21h)These equations are defined within a constant which is determined from the initial conditions at the level of the transmitter. Theyare valid for all f ≥ f .
4. Tidal potential
In the fluid rotating frame, the perturbing tidal potential can be expressed as [25] U tidal = − ∞ (cid:88) l = l ! (cid:34) h L E L ( t ) − I < L > M h j E jL ( t ) (cid:35) , (B22)where E L ( t ) are the tidal moments being time dependent STF tensors. They are defined by E L ( t ) = − ∂ L U ¬ P ( t , ), (B23)where the gravitational potential created by the external bodies U ¬ P is di ff erentiated l times with respect to h and the result isevaluated at the central planet’s center-of-mass (i.e. h = ).Following [25], U ¬ P is given by U ¬ P ( t , h ) = − G (cid:88) A (cid:44) P (cid:90) A ρ ( t , x ) | h + r P ( t ) − x | d x , (B24)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 ranges 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 separationbetween 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 lying in 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 Taylorseries expansion. After some algebra, we arrive to U ¬ P ( t , h ) = − G (cid:88) A (cid:44) P ∞ (cid:88) l (cid:48) = ( − l (cid:48) l (cid:48) ! I < L (cid:48) > A ∂ L (cid:48) | r AP + h | − , (B25)7We recall that I < L (cid:48) > A 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 = ). This notation is shorten noticing that ∂ L (cid:48) | r AP + h | − = ∂ L (cid:48) r − AP , where thederivative is taken at the extremity of the vector r AP ( t ).From Eq. (B25), (B22), and (B23), we can define the perturbing external tidal potential U tidal = − G (cid:88) A (cid:44) P ∞ (cid:88) l = ∞ (cid:88) l (cid:48) = ( − l (cid:48) l ! l (cid:48) ! I < L (cid:48) > A (cid:34) h L ∂ LL (cid:48) r − AP − I < L > M h j ∂ jLL (cid:48) r − AP (cid:35) . (B26)That expression is very general since it takes into account the tidal e ff ects due to the total gravitational potential (all multipolemoments included) of all the external bodies acting on the total gravitational potential of the central planet.At linear order in the generalized potential, the refractive profile and the j th component of the gradient of the non-constantchange in the index of refraction are explicitly written as n = n + α G (cid:88) A (cid:44) P ∞ (cid:88) l = ∞ (cid:88) l (cid:48) = ( − l (cid:48) l ! l (cid:48) ! I < L (cid:48) > A (cid:34) h L ∂ LL (cid:48) r − AP − I < L > M h j ∂ jLL (cid:48) r − AP (cid:35) , (B27a) ∂ j δ n = α G (cid:88) A (cid:44) P ∞ (cid:88) l = ∞ (cid:88) l (cid:48) = ( − l (cid:48) l ! l (cid:48) ! I < L (cid:48) > A (cid:34) h L ∂ jLL (cid:48) r − AP − I < L > M h k ∂ jkLL (cid:48) r − AP (cid:35) . (B27b)For the application that we are to consider now, we admit only one perturbing body, A , and we define the separation vectorbetween P and A as r ( t ) ≡ − r AP ( t ). We consider that the central planet is close to spherical symmetry, also the couplingbetween its higher multipole moments I < L > and the other STF tensors can safely be neglected. In addition, we consider that theperturbing body is a satellite of the main planet P , and that M A (cid:28) M , with M the mass of the central planet. Therefore we cansafely consider only the monopole contribution of A , i.e. l (cid:48) =
0. In addition, assuming that inter-body distance is much morelarger than the typical scale of the ray closest approach to P (we are interested in the light propagation in the atmosphere of P ), itis su ffi cient to keep only the leading quadrupole term in the development, i.e. l =
2. For a matter of simplicity, we also considerthat the satellite follows a circular orbit of radius r which lies in the equatorial plane of the central planet. Then, in the referenceframe, the direction of the position vector of the satellite is given at any time t byˆ r ( t ) ≡ r ( t ) / r = cos Nt ˆ X + sin Nt ˆ Y , (B28)where the mean motion of the satellite around the central body is given by Kepler third law of motion N (cid:39) ( µ/ r ) / (For thisapplication, N is not the refractivity anymore). Making use of Eqs. (19a)-(19c), the orbital motion of the satellite is given, in thepolar frame, by ˆ r ( t ) = (cid:2) cos λ ( t ) cos( f + ω ) + cos ι sin λ ( t ) sin( f + ω ) (cid:3) ˆ ρ − (cid:2) cos λ ( t ) sin( f + ω ) − cos ι sin λ ( t ) cos( f + ω ) (cid:3) ˆ τ − sin ι sin λ ( t ) ˆ σ , (B29)in which we have defined λ ( t ) ≡ Nt − Ω for being the longitude of the perturbing body as measured from the intersection betweenthe propagation plane and the equatorial plane (longitude of the node); this angle remains within the equatorial plane of P sincethe perturbing body describes an equatorial circular orbit.Following all the simplifications, the change in the index of refraction and its gradient now read δ n = αµ A h kl ∂ kl r − AP , (B30a) ∂ j δ n = αµ A h kl ∂ jkl r − AP . (B30b)Making use of properties about the STF tensors [cf. Eqs. (B14)], the multi-index derivatives now read ∂ i j r − AP = − r ( δ i j − r i j ), (B31a) ∂ i jk r − AP = − r ( δ i j ˆ r k + δ ik ˆ r j + δ jk ˆ r i − r i jk ). (B31b)In those expressions, the components of the unit vector directed along r are noted ˆ r j ≡ r j / r and can be inferred from Eq. (B29)in the polar basis.8We can now express the components of the perturbing gradient and the non-constant change in the index of refraction f pert = − αµ A r h r A , (B32a) N = − αµ A np h r A N . (B32b)Here, we have introduced the unitless vector A and the unitless parameter A N to simplify the notations. They are explicitlygiven by A = r · ˆ ρ ) ˆ ρ + ˆ r − r · ˆ ρ ) ˆ r , (B33a) A N = − r · ˆ ρ ) , (B33b)and once expressed in term of the fundamental angles of the problem they read A R = A · ˆ ρ = (cid:2) cos( f + ω ) cos λ + cos ι sin( f + ω ) sin λ (cid:3)(cid:16) − (cid:2) cos( f + ω ) cos λ + cos ι sin( f + ω ) sin λ (cid:3) (cid:17) , (B34a) A T = A · ˆ τ = (cid:2) cos ι cos( f + ω ) sin λ − sin( f + ω ) cos λ (cid:3)(cid:16) − (cid:2) cos( f + ω ) cos λ + cos ι sin( f + ω ) sin λ (cid:3) (cid:17) , (B34b) A S = A · ˆ σ = − sin ι sin λ (cid:16) − (cid:2) cos( f + ω ) cos λ + cos ι sin( f + ω ) sin λ (cid:3) (cid:17) , (B34c) A N = − (cid:2) cos( f + ω ) cos λ + cos ι sin( f + ω ) sin λ (cid:3) . (B34d)One can now insert Eqs. (B32) into Eqs. (55). Before integrating, we are to further simplify the problem by assuming a raytraveling in the equator of the central planet and passing through the bulge raised by the perturbing body. In such a case, ι = ω and Ω separately, we consider thelongitude of the closest approach, (cid:36) = ω + Ω , which is measured from the ˆ X -direction. Keeping the leading order in 1 / e , andintegrating over the true anomaly leads to the following set of equations describing the change in hyperbolic elements followinga gravitational tidal stress p ( f ) = − e η ( αµ ) r µ A µ sec f (cid:2)
15 cos(4 f + Nt − (cid:36) ) −
15 cos(4 f − Nt + (cid:36) ) + cos(4 f + Nt − (cid:36) ) − cos(4 f − Nt + (cid:36) ) +
60 cos(2 f + Nt − (cid:36) ) +
60 cos(2 f − Nt + (cid:36) ) + f + Nt − (cid:36) ) − f − Nt + (cid:36) ) +
90 cos(3 Nt − (cid:36) ) + Nt − (cid:36) ) (cid:3) , (B35a) e ( f ) = − e η ( αµ ) r µ A µ sec f (cid:2)
25 cos(4 f + Nt − (cid:36) ) −
25 cos(4 f − Nt + (cid:36) ) + cos(4 f + Nt − (cid:36) ) − cos(4 f − Nt + (cid:36) ) +
100 cos(2 f + Nt − (cid:36) ) +
80 cos(2 f − Nt + (cid:36) ) + f + Nt − (cid:36) ) −
16 cos(2 f − Nt + (cid:36) ) +
150 cos(3 Nt − (cid:36) ) − Nt − (cid:36) ) (cid:3) , (B35b) (cid:36) ( f ) = − e η ( αµ ) r µ A µ sec f (cid:2) f + Nt − (cid:36) ) + f − Nt + (cid:36) ) + f + Nt − (cid:36) ) + f + Nt − (cid:36) ) + f − Nt + (cid:36) ) +
15 sin( f + Nt − (cid:36) ) −
15 sin( f − Nt + (cid:36) ) (cid:3) , (B35c) δ s ( f ) = − e η ( αµ ) r µ A µ sec f (cid:2)
15 sin(5 f + Nt − (cid:36) ) +
15 sin(5 f − Nt + (cid:36) ) + f + Nt − (cid:36) ) + f − Nt + (cid:36) ) +
75 sin(3 f + Nt − (cid:36) ) +
75 sin(3 f − Nt + (cid:36) ) +
35 sin(3 f + Nt − (cid:36) ) +
35 sin(3 f − Nt + (cid:36) ) +
150 sin( f + Nt − (cid:36) ) −
150 sin( f − Nt + (cid:36) ) +
70 sin( f + Nt − (cid:36) ) +
10 sin( f − Nt + (cid:36) ) (cid:3) , (B35d) δ t ( f ) = η c δ s ( f ), (B35e) δψ ( f ) = (cid:36) ( f ). (B35f)These equations are defined within a constant which is determined from the initial conditions at the level of the transmitter. Theyare 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 [25] that the first ordere ff ects can be assessed by replacing E L ( t ) by E L ( t − τ ) where τ is called the time-lag and is function of the Energy dissipation. Acomputation of E L at the instant t − τ , changes the instantaneous position of the perturbing body [cf. Eq. (B28)] into a delayed9position, r ( t − τ ). Most of the time, τ is a small quantity compared to orbital characteristic time scale. Consequently, the orbitalvelocity of the perturbing body is usually too small to change significantly it’s apparent position during the time interval τ . Inthis condition, only internal process occurring on time scale faster than the orbital motion can significantly change E L . Amongthese process, the rigid rotation of the atmosphere appears to be a good candidate for fast rotators. The delayed position cantherefore be approximated by r ( t − τ ) = r ( t ) + w τ r (cid:0) cos Nt ˆ Y − sin Nt ˆ X (cid:1) ,where w is the magnitude of the angular velocity which has been defined to be aligned with the normal to the equator of reference.The delayed unit vector is given by ˆ r ( t − τ ) = ˆ r ( t ) + w τ ˆ r ∗ ( t ),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 inthe index of refraction are given at first order in τ by f ∗ pert = − w τ αµ A r h r A ∗ , (B36a) N ∗ = − w τ αµ A np h r A ∗N . (B36b)As before, we have introduced the trigonometric coe ffi cients A ∗ and A ∗N which are given by A ∗ = r ∗ · ˆ ρ ) ˆ ρ + ˆ r ∗ − r · ˆ ρ ) (cid:104) (ˆ r · ˆ ρ )ˆ r ∗ + r ∗ · ˆ ρ )ˆ r (cid:105) , (B37a) A ∗N = − r · ˆ ρ )(ˆ r ∗ · ˆ ρ ). (B37b)One can notice from the form of Eqs. (B36), that( f ∗ pert · ˆ ρ )( f pert · ˆ ρ ) = w τ A ∗R A R , (idem along ˆ τ and ˆ σ ), N ∗ N = w τ A ∗N A N .Then, it is immediate to show that the change in hyperbolic elements due to the non-elastic response of the atmosphere to tidalstress is C ∗ C ∝ w τ , (B38)where C represents the vector of the hyperbolic elements.
5. Horizontal gradients
Let us consider the Earth centered frame ( ˆ X , ˆ Y , ˆ Z ) with the ˆ Z -axis directed through the North pole, the ( ˆ X , ˆ Y ) plane superim-posed within the Earth equator, and the ˆ X -axis pointed through the intersection between the equator and the Greenwich meridian.We assume that the atmosphere is filled with a stationary medium with respect to the frame attached to the Earth, whence Eq. (4)is valid in ( ˆ X , ˆ Y , ˆ Z ). A station on Earth can be located thanks to its spherical coordinates ( ϕ , λ , R ), with R the Earth radius, ϕ the colatitude (measured from the North pole), and λ the azimuth of the station . Any point along the light ray trajectory canbe 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 tothe 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 [5] or [6] n ( h , ϕ, λ ) = n ( h ) + ( ϕ − ϕ ) n NS + ( λ − λ ) n WE , (B39) Here, we consider the case of a receiving station but the same work can be applied for a transmitting station. n ( h ) is the hyperbolic index of refraction, and n NS and n WE are respectively given by ∂ n /∂ϕ | h ,λ and ∂ n /∂λ | h ,ϕ . In otherwords, n NS represents the variation of the horizontal refraction along the North-South direction and n WE represents the variationof the horizontal refraction along West-East direction. Here, we aim at computing the e ff ects 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 inspherical coordinates we end up with f pert = n NS h ˆ u ϕ + n WE h sin ϕ ˆ u λ , (B40)where we have assumed that the radial variations of n NS and n WE are of second order. The two unit vectors ˆ u ϕ and ˆ u λ are part ofthe usual spherical vectorial basis which is given by the triad of vectors ( ˆ u ϕ , ˆ u λ , ˆ ρ ). It is seen that the perturbing gradient doesindeed not possess any radial component if we neglect the radial variations of n NS and n WE .Now let us write the transverse and the normal components of the perturbing gradient by expressing ˆ u ϕ and ˆ u λ into the polarbasis ( ˆ ρ , ˆ τ , ˆ σ ).For any arbitrary point, we can change from the spherical to the Cartesian coordinates thanks to two angles ϕ and λ ˆ ρ = sin ϕ cos λ ˆ X + sin ϕ sin λ ˆ Y + cos ϕ ˆ Z , (B41a)ˆ u ϕ = cos ϕ cos λ ˆ X + cos ϕ sin λ ˆ Y − sin ϕ ˆ Z , (B41b)ˆ u λ = − sin λ ˆ X + cos λ ˆ Y . (B41c)The radial direction can be compared to Eq. (19a) which reveals the following relationshipscos ϕ = sin ι sin( f + ω ), (B42a)cos λ sin ϕ = cos Ω cos( f + ω ) − cos ι sin Ω sin( f + ω ), (B42b)sin λ sin ϕ = sin Ω cos( f + ω ) + cos ι cos Ω sin( f + ω ), (B42c)sin ϕ = (cid:112) A ( f ) cos( f + ω ). (B42d)We have introduced the trigonometric function A ( f ) which is defined in Tab. I. From Eqs. (B41) and (19) and by making useof Eqs. (B42), we find ˆ u ϕ = − sin ι (cid:112) A ( f ) ˆ τ − ι (cid:112) A ( f ) ˆ σ , (B43a)ˆ u λ = ι (cid:112) A ( f ) ˆ τ − sin ι (cid:112) A ( f ) ˆ σ , (B43b)where the A l ( f ) coe ffi cients are given in Tab. I. Those equations allow us to express the perturbing gradient into the polar basisand allow us to derive its transverse and normal components.At this level, we also need an equation for the change in the index of refraction in term of hyperbolic elements. That equationcan be obtained from Eqs. (B42). Indeed at first order in f − f we deduce ϕ − ϕ = − sin ι cos( f + ω ) (cid:112) A ( f ) ( f − f ), (B44a) λ − λ = cos ι A ( f ) cos ( f + ω ) ( f − f ), (B44b)which both can be inserted into Eq. (B39). We sum up the results into the following relations (the radial component is null) T = (cid:104) n WE cos ι − n NS (cid:112) A ( f ) sin ι cos( f + ω ) (cid:105) (1 + e cos κ f ) p A ( f ) , (B45a) S = − (cid:20) n NS cos ι (cid:112) A ( f ) + n WE sin ι A ( f ) cos( f + ω ) (cid:21) (1 + e cos κ f ) p , (B45b) N = (cid:20) n WE cos ι A ( f ) cos ( f + ω ) − n NS sin ι cos( f + ω ) (cid:112) A ( f ) (cid:21) ( f − f ) np . (B45c)Finally, the e ff ect on the hyperbolic elements is given by inserting Eqs. (B45) into Eqs. (55). Before integrating we are tosimplify the problem by assuming source tracking at the meridian of the site of observation. For this geometry, the propagation1 l A l ( f )1 1 + cos ι tan ( f + ω )2 1 − sin ι sin ( f + ω )3 3 + cos 2 ι + ι cos 2( f + ω )4 3 + cos 2( f + ω ) + ι sin ( f + ω )TABLE I. Definition of the trigonometric functions A l ( f ) which are used in this section. We do not specify the dependencies in ι and ω , sincethey are not expected to vary as quickly as f at first order in the components of the perturbation. plane of photons is orthogonal to the plane attached to Earth’s equator ( ι = π/ π/
2, the longitude of the node is the same as theazimuth of both the source and the station (
Ω = λ ). Keeping the leading order in 1 / e , and integrating over the true anomaly leadsto the following expressions p ( f ) = − e η αµ n NS tan f , (B46a) e ( f ) = − e η n NS ( f + tan f ), (B46b) ι ( f ) = − n WE η tan f , (B46c) Ω ( f ) = n WE η (cid:104) cot ω tan f + (cid:16) ln | cos( f + ω ) | − ln | cos f | (cid:17) csc ω (cid:105) , (B46d) ω ( f ) = n NS η ln | cos f | , (B46e) δ s ( f ) = − e η αµ n NS (cid:104) − ( f − f ) sin 2 f + (1 + cos 2 f ) ln | cos f | (cid:105) sec f , (B46f) δ t ( f ) = − e η c αµ n NS (cid:104) sec f + (cid:16) ln | cos f | − ( f − f ) tan f (cid:17)(cid:105) , (B46g) δψ ( f ) = ω ( f ). (B46h)These equations are defined within a constant which is determined from the initial conditions at the level of the transmitter. Theyare valid for all f ≥ f . Appendix C: Doppler shift and refractive bending
In this chapter, we show how the Doppler frequency shift can be related to the refractive bending.
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 (cid:34) ν ν (cid:35) vac = Γ Γ (1 − ˆ u · s vac ) (1 − ˆ u · s vac ) .Inserting the definition (2) into Eq. (9), we arrive to ν ν = (cid:34) ν ν (cid:35) vac − n ˆ u · δ s − ˆ u · s vac − ( n −
1) ˆ u · s vac − ˆ u · s vac − n ˆ u · δ s − ˆ u · s vac − ( n −
1) ˆ u · s vac − ˆ u · s vac , (C1)2where 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 respectto s vac . In vacuum, δ s = , and n − =
0, thus ν /ν reduces to [ ν /ν ] vac .For application within the Solar System, Eq. (C1) can be simplified considering the small velocity approximation ( | ˆ u | (cid:28) n − (cid:28) | δ s | (cid:28)
1, and substitutions into Eq. (C1) leads to the following approximation ν ν (cid:30) (cid:34) ν ν (cid:35) vac − ≈ n β ( ˆ u × s vac ) · ˆ σ + n β ( ˆ u × s vac ) · ˆ σ + ( n −
1) ˆ u · s vac − ( n −
1) ˆ u · s vac .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 from 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 berespectively zero or unity in the limiting cases where i) the receiver is at infinity or ii) the receiver and the transmitter are at equaldistance from the center of mass ( Λ → {
0; 1 } when h → { + ∞ ; h } ). Moreover, from Fig. 1 it is seen that the total refractivebending can be expressed as (cid:15) ≡ β − β , also we end up with ν ν (cid:30) (cid:34) ν ν (cid:35) vac − ≈ − (cid:15) n + Λ ( ˆ u × s vac ) · ˆ σ + Λ (cid:15) n + Λ ( ˆ u × s vac ) · ˆ σ + ( n −
1) ˆ u · s vac − ( n −
1) ˆ u · s vac . (C2)This equation shows that in first approximation, the Doppler shift due to atmospheric e ff ects is controlled by the total refractivebending, (cid:15) . As discussed in appendix A, (cid:15) is defined as the net change in the argument of the refractive bending, ψ , between thetransmitter and the receiver ( (cid:15) ≡ ψ − ψ ), where the exact expression of ψ has been given in Eq. (A36).
2. Simplifications
For one-way Earth AO experiments [14, 15], the index of refraction at the level of the two spacecrafts orbiting the Earth canbe considered for being unity ( n = n = ν ν (cid:30) (cid:34) ν ν (cid:35) vac − ≈ − (cid:15) + Λ (cid:16) [ ˆ u − Λ ˆ u ] × s vac (cid:17) · ˆ σ , (C3)when ˆ σ = ˆ σ = ˆ σ .For one-way planetary AO experiments (see e.g. [11, 49]), the index of refraction at the level of the spacecraft can be takenfor being unity ( n = Λ →
0, andthe simplified Doppler shift expression becomes ν ν (cid:30) (cid:34) ν ν (cid:35) vac − ≈ − (cid:15) ( ˆ u × s vac ) · ˆ σ − ( n −
1) ˆ u · s vac . (C4)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 constanto ffff