© ESO 2021

1 Introduction

Theoretical problems dealing with time and frequency transfers require knowledge of the function relating the (coordinate) time transfer to the coordinate time at reception and to the spatial coordinates of the reception and emission point-events. Such a function is called a reception time transfer function. An emission time transfer function can be introduced too. The formalism designed to determine the time transfer functions is called the time transfer functions formalism. It was first introduced by Linet & Teyssandier (2002) relying on the theory of the world function developed by Synge (1960). Later on a general post-Minkowskian expansion of the world and the time transfer functions was proposed by Le Poncin-Lafitte et al. (2004) and then the method was refined by Teyssandier & Le Poncin-Lafitte (2008) thanks to a simplified iterative procedure. The time transfer functions formalism confers a great advantage in that it spares the trouble of integrating the geodesic equation, which usually leads to heavy calculations, especially beyond the post-Minkowskian regime (Richter & Matzner 1983; Brumberg 1987). Until recently, the time transfer functions formalism was systematically applied to the physical spacetime considering only gravitational effects on a ray of light propagating in a vacuum. However, Bourgoin (2020) showed that theoretical problems dealing with optical rays propagating into flowing dielectric medium could also be solved making use of the time transfer functions formalism by means of a powerful theoretical tool known as the optical metric or the Gordon’s metric of spacetime.

When a ray of light is propagating into an optical medium its trajectory does not follow a null geodesic path of the physical spacetime because light and matter interact. Conventionally, the real trajectory is therefore determined by solving Maxwell’s equation within the framework of geometrical optics. However, another interesting possibility initially proposed by Gordon (1923) is to introduce an artificial optical metric which implicitly accounts for the interaction between light and matter such that optical rays follow null geodesics of that new optical metric. In other words, refractivity is treated as spacetime curvature in the optical metric. Therefore, the time and frequency transfers can still be determined with the time transfer functions formalism even when considering a ray of light crossing through an optical medium.

Occultation experiments are an example of an observing technique requiring careful treatment of refractivity while modeling the time and frequency transfers. The method consists in remotely measuring the physical properties of a planetary atmosphere while the source of an electromagnetic signal is being occulted by the atmosphere. When the source is a radio signal emitted by a spacecraft’s antenna the experiment is called an atmospheric radio occultation (Kliore et al. 1965; Fjeldbo & Eshleman 1965, 1968; Lindal et al. 1985, 1987; Lindal 1992; Schinder et al. 2012, 2015), whereas it is called an atmospheric stellar occultation when the source is made of visible or near-infrared light emitted by a distant star (Roques et al. 1994; Sicardy et al. 2006). In practice, two methods are usually employed for processing occultation data, namely the Abel inversion for spherical symmetry (Phinney & Anderson 1968; Steiner et al. 1999) and the numerical ray-tracing for any generic case (Schinder et al. 2015). The former is an exact expression providing the index of the refraction profile directly from the bending angle which is itself retrieved from the frequency transfer. The latter consists in a numerical integration of the equations for optical rays across a layered atmosphere. The refractivity in each layer and the initial pointing direction are iteratively determined such that the computed frequency coincides with the observed frequency. While the numerical ray-tracing method is the most general one, it does not provide a comprehensive description of the light path and requires significant computational time. On the other hand, if the Abel inversion method does not require numerically solving for the equations for optical rays, it can only be applied to spherically symmetric atmospheres and cannot account for the atmospheric time delay which eventually affects the determination of the emitter’s position (i.e., the spacecraft in a one-way downlink configuration during a radio occultation event). In addition, as in numerical ray-tracing, the Abel inversion does not provide a comprehensive description of the light path. An analytical expression describing the time and frequency transfers for radio occultation experiments would allow these issues to be overcome.

With this goal in mind, Bourgoin et al. (2019) proposed a new approach exploiting similarities between equations of geometrical optics and equations of celestial mechanics. It consists in expressing the equations of geometrical optics as a set of first-order perturbation equations (similar to Gauss equations of celestial mechanics) which are better suited for finding analytical solutions beyond the spherical symmetry assumption. However, even if the perturbation equations can be solved more easily than the equations of geometrical optics it is still a challenging task to solve for the second-order expressions or to incorporate the light-dragging effect caused by the motion of the optical medium.

In this paper, we make use of the time transfer functions formalism considering an optical metric in order to accurately model the time and frequency transfers for occultation experiments involving a flowing spherically symmetric atmosphere. The determination of the time transfer functions allows us to account for the atmospheric time delay. In addition, the formalism being fully covariant, it naturally accounts for the light-dragging effect due to the motion of the optical medium and constitutes in that sense a significant improvement with respect to the perturbation equations approach.

The paper is organized as follows. Section 2 lists the notations and assumptions we make throughout the paper. Section 3 recalls some basic information about relativistic geometrical optics and allows us to define the index of refraction, the refractivity, and the optical metric. The equations for optical rays propagating in an isotropic dispersive medium are derived in the same section. The reader who is already familiar with relativistic geometrical optics can skip this section. Section 4 recalls basic information about the time transfer functions formalism and introduces the refractive delay function and the post-Minkowskian parameter N₀. The integral form of the delay function is given at any post-Minkowskian order. Section 5 is an application to occultation experiments involving steadily rotating and spherically symmetric atmospheres. The refractivity profile is built in Sect. 6 considering an exponential pressure profile and a polynomial temperature profile of arbitrary degree. The refractive delay function is finally solved at first post-Minkowskian order in the limit where the angular velocity of the optical medium is small with respect to the speed of light in a vacuum. The expressions for the time and frequency transfers are given explicitly. Section 7 assesses the accuracy of the first-order solutions for the time and frequency transfers by comparing them to results of a numerical integration of the equations for optical rays propagating into a nondispersive isotropic medium. Finally, we give our conclusions in Sect. 8. Appendix A is a discussion about the Abel transform method for finding the refractivity from the frequency transfer while considering the light-dragging effect.

2 General assumptions and notations

In this paper, the influence of gravity on the propagation of light is regarded as negligible, and so the physical metric g of spacetime is assumed to be a Minkowski metric. Greek indices run from 0 to 3, and Latin indices run from 1 to 3. We systematically make use of an orthonormal Cartesian coordinate system (x^μ) = (x⁰, xⁱ), and so the components of the physical metric may be written as $$g_{\mu \nu }={\eta }_{\mu \nu }\text{,}$$(1)

where $${\eta }_{\mu \nu }=\text{diag}(+1,-1,-1,-1)\text{.}$$(2)

We put x⁰ = ct, with c being the speed of light in a vacuum and t a time coordinate, and we use x to denote the triplet of spatial coordinates (x¹, x², x³). More generally, we use the notation a = (aⁱ) = (a¹, a², a³) for a triplet constituted by the spatial components of a 4-vector, and $\underset{\_}{b}=(b_i)=(b_1,b_2,b_3)$ for a triplet built with the spatial components of a covariant 4-vector.

Given the triplet a, b, and $\underset{\_}{c}$, the usual Euclidean scalar product a ⋅b is denoted aⁱbⁱ = δ_ikaⁱb^k, where δ_ik is the Kronecker delta. Similarly, $a\cdot \underset{\_}{c}$ denotes the quantity aⁱc_i = δ_ikaⁱc_k. In each case, Einstein’s summation convention on repeated indices is used. Furthermore, ∥a∥ denotes the Euclidean norm of a: $\Vert a\Vert =\sqrt{a\cdot a}$. Similarly, $\Vert \underset{\_}{c}\Vert$ denotes the quantity $\Vert \underset{\_}{c}\Vert =\sqrt{\underset{\_}{c}\cdot \underset{\_}{c}}$.

For the sake of legibility, we employ ${(f)}_x$ or ${[f]}_x$ instead of f(x) whenever necessary. When a quantity f(x) is to be evaluated at two point-events x_A and x_B, we employ ${(f)}_{A/B}$ to denote f(x_A) and f(x_B), respectively. The partial differentiation of f with respect to x^μ is denoted ∂_μf or f_,μ.

3 Relativistic geometrical optics

This paper is devoted to the propagation of light rays through a linear, isotropic, and nondispersive medium filling a spatially bounded region $\mathcal{D}$ in Minkowski spacetime. The region exterior to $\mathcal{D}$ is assumed to be empty of any matter. According to these assumptions, the basic relations of the present paper are inferred by substituting the metric components g_μν from Eq. (1) into the equations supplied in Bourgoin (2020).

The electromagnetic properties of the medium are characterized by two scalar functions, the permittivity ϵ(x) and the permeability μ(x). The index of refraction of the medium is the scalar function defined by the well-known relationship $$n(x)=c\sqrt{ϵ(x)\mu (x)}\text{.}$$(3)

Moreover, it is assumed that the medium is made of a fluid schematized by a flow of particles which are not colliding. The unit 4-velocity vector of a particle of the fluid at a point-event x of its world-line is denoted w^μ(x). Outside $\mathcal{D}$, the permittivity and the permeability reduce to their vacuum values ϵ(x) = ϵ₀ and μ(x) = μ₀, respectively. As $c={({ϵ}_0{\mu }_0)}^{-1/2}$, the index of refraction then reduces to n(x) = 1. The expression for the refractivity is obtained by subtracting its vacuum value from the index of refraction, namely $$N(x)=n(x)-1.$$(4)

In the context of the geometrical optics approximation, the light rays propagating through our medium are the bicharacteristic curves of the so-called eikonal equation, which reads $${\overline{g}}^{\mu \nu }{\partial }_{\mu }\mathcal{S}{\partial }_{\nu }\mathcal{S}=0\text{,}$$(5)

where $\mathcal{S}(x)$ is the eikonal function and ḡ^μν is the contravariant tensor defined by $${\overline{g}}^{\mu \nu }=g^{\mu \nu }+{\kappa }^{\mu \nu }\text{,}{\kappa }^{\mu \nu }=(n^2-1)w^{\mu }{w}^{\nu }\text{.}$$(6)

Let us use ḡ_μν to denote the quantities such that $${\overline{g}}_{\mu \alpha }{\overline{g}}^{\alpha \nu }={\delta }_{\mu }^{\nu }\text{.}$$(7)

An elementary calculation leads to $${\overline{g}}_{\mu \nu }=g_{\mu \nu }+{\gamma }_{\mu \nu }\text{,}{\gamma }_{\mu \nu }=-\left(1-\frac{1}{n^2}\right)w_{\mu }{w}_{\nu }\text{.}$$(8)

The quantities ḡ_μν can be regarded as the components of a Lorentzian metric ḡ defined on the region $\mathcal{D}$. This new metric is called the optical metric associated with the optical medium, terminology justified by the following considerations, which were derived by Gordon (1923), Quan (1957), and Ehlers (1967).

Let us define the 4-wave covector field k_μ as $$k_{\mu }(x)={\partial }_{\mu }\mathcal{S}(x)\text{.}$$(9)

A contravariant vector field k^μ can be associated with this covector by putting $$k^{\mu }=g^{\mu \nu }{k}_{\nu }\text{.}$$(10)

As recalled by Bourgoin (2020), it is convenient to define the contravariant vector field ${\overline{k}}^{\mu }$ as $${\overline{k}}^{\mu }={\overline{g}}^{\mu \nu }{k}_{\nu }\text{.}$$(11)

Indeed, the light rays x^μ = x^μ(ζ) associated with a solution $\mathcal{S}$ of the eikonal Eq. (5) are the integral curves of the vector field ${\overline{k}}^{\mu }$, namely, they are solutions of the differential system (Perlick 2000) $$\frac{\text{d}{x}^{\mu }}{\text{d}\zeta }={\overline{k}}^{\mu }(x(\zeta ))={\overline{g}}^{\mu \nu }{\partial }_{\nu }\mathcal{S}(x(\zeta ))\text{.}$$(12)

A classical calculation shows that the solutions of Eq. (12) are null geodesics of the optical metric ḡ (Synge 1960; Perlick 2000), ζ being an affine parameter. The null character of the rays is directly inferred from Eqs. (7), (12), and (5). We have indeed: $${\overline{g}}_{\mu \nu }\frac{\text{d}{x}^{\mu }}{\text{d}\zeta }\frac{\text{d}{x}^{\nu }}{\text{d}\zeta }={\overline{g}}^{\alpha \beta }{\partial }_{\alpha }\mathcal{S}{\partial }_{\beta }\mathcal{S}=0.$$(13)

The eikonal Eq. (5) is the Jacobi equation associated with the Hamiltonian $$H(x^{\mu },k_{\nu })=\frac{1}{2}{\overline{g}}^{\alpha \beta }(x)k_{\alpha }{k}_{\beta }\text{,}$$(14)

where k_ν must be regarded as conjugate canonical variables of x^μ. As a consequence, the light rays in the region $\mathcal{D}$ can be considered as solutions of the set of canonical equations (Quan 1957) $$\begin{array}{ll}\frac{\text{d}{x}^{\mu }}{\text{d}\zeta }\hfill & =\frac{\partial H}{\partial k_{\mu }}=k^{\mu }+(n^2-1)(w^{\nu }{k}_{\nu })w^{\mu }\text{,}\hfill \\ \frac{\text{d}{k}_{\mu }}{\text{d}\zeta }\hfill & =-\frac{\partial H}{\partial x^{\mu }}=-n{n}_{,\mu }{(w^{\nu }{k}_{\nu })}^2-(n^2-1)(w^{\nu }{k}_{\nu })w_{,\mu }^{\alpha }{k}_{\alpha }\text{.}\hfill \end{array}$$

It must be noted that the eikonal function is constant along a light ray. Indeed, it follows from Eqs. (12) and (13) that $$\frac{\text{d}\mathcal{S}}{\text{d}\zeta }=\frac{\text{d}{x}^{\mu }}{\text{d}\zeta }{\partial }_{\mu }\mathcal{S}={\overline{g}}^{\mu \nu }{\partial }_{\mu }\mathcal{S}{\partial }_{\nu }\mathcal{S}=0.$$(16)

This property is at the core of the procedure developed in the following section.

For numerical integration, it is relevant to separate space and time components in Eq. (15). Let l_i be the components defined by $$l_i=\frac{k_i}{k_0}\text{,}$$(17)

and let ℓ be a new parametrization of the curves for optical rays (Schinder et al. 2015) such that $$\text{d}\mathcal{l}=n{k}_0\text{d}\zeta \text{.}$$(18)

After inserting these new quantities into the set of canonical Eq. (15), we find for the time components: $$\begin{array}{ll}\frac{\text{d}{x}^0}{\text{d}\mathcal{l}}=\hfill & \frac{1}{n}\left[1+(n^2-1)\left(w^0+w^k{l}_k\right)w^0\right]\text{,}\hfill \\ \frac{\text{d}\ln \Vert k_0\Vert }{\text{d}\mathcal{l}}=\hfill & -n_{,0}{\left(w^0+w^k{l}_k\right)}^2\hfill \\ \hfill & -\frac{(n^2-1)}{n}\left(w^0+w^k{l}_k\right)\left(w_{,0}^0+w_{,0}^k{l}_k\right)\text{,}\hfill \end{array}$$

and, for the space components: $$\begin{array}{ll}\frac{\text{d}{x}^i}{\text{d}\mathcal{l}}=\hfill & -\frac{1}{n}\left[l_i-(n^2-1)\left(w^0+w^k{l}_k\right)w^i\right]\text{,}\hfill \\ \frac{\text{d}{l}_i}{\text{d}\mathcal{l}}=\hfill & -n_{,i}{\left(w^0+w^k{l}_k\right)}^2\hfill \\ \hfill & -\frac{(n^2-1)}{n}\left(w^0+w^k{l}_k\right)\left(w_{,i}^0+w_{,i}^k{l}_k\right)-\frac{\text{d}\ln \Vert k_0\Vert }{\text{d}\mathcal{l}}{l}_i\text{.}\hfill \end{array}$$

Equations (19) may be particularly interesting in the case where the index of refraction n and the unit 4-velocity w^μ do not depend on time, and we use them in our attempt at numerical integration (see Sect. 7). It is worthy of note that in this case the component k₀ is constant along each light ray.

4 Refractive delay function

In this section, we introduce the formalism of time transfer functions that we apply to optical metric in order to describe refraction due to a linear, isotropic, and nondispersive medium.

Fig. 1 Illustration of an occultation event in spacetime. The past light cone $\mathcal{C}(x_B)$ of the reception point-event xB intersects ${\mathcal{C}}_A$, the world-line of the emitter, at the point-event xA. The zeroth-order null light path (red line) joining the emission point-event xA to the reception point-event xB lies on the surface of $\mathcal{C}(x_B)$. The domain $\mathcal{D}$ represents the limit of the refractive region while x− and x+ are the intersection point-events between $\mathcal{D}$ and the zeroth-order null light path. The path of integration (dashed red line) is limited to the portion of the zeroth-order null light path crossing through $\mathcal{D}$.

4.1 Timetransfer functions formalism

Let us consider a light ray Γ_AB starting from an emission point-event $(x_A^0,x_A)$ and arriving at a reception point-event $(x_B^0,x_B)$. We consider that a part of Γ_AB travels through the domain $\mathcal{D}$, while the other part travels through a vacuum, that is a medium such that n = 1 (see Fig. 1).

It is immediately inferred from Eq. (16) that the phase $\mathcal{S}$ is constant along Γ_AB. As a consequence, the coordinates of the emission point x_A and of the reception point x_B are such that the following relation is satisfied: $$\mathcal{S}(x_B^0,x_B)-\mathcal{S}(x_A^0,x_A)=0.$$(20)

Equation (20) shows that $x_A^0$ may be considered as an implicit function of x_A and the coordinates of the reception point, namely $x_B^0$ and x_B. Hence, it is appropriate to introduce $\mathcal{T}$ the reception time transfer function associated with Γ_AB as $$x_B^0-x_A^0=c\mathcal{T}(x_A,x_B^0,x_B)\text{.}$$(21)

An emission time transfer function can be introduced too. Hereafter, we merely consider the case at reception which is better suited for a downlink one-way transfer during an occultation event with a radio signal that is recorded when received.

A relevant theorem about the components of the 4-wave covector can be directly inferred from Eqs. (20) and (21); see Le Poncin-Lafitte et al. (2004). Indeed after writing Eq. (21) as $$x_A^0=x_B^0-c\mathcal{T}(x_A,x_B^0,x_B)\text{,}$$(22)

and then inserting this relation into Eq. (20), we get the following expression: $$\mathcal{S}(x_B^0,x_B)-\mathcal{S}(x_B^0-c\mathcal{T}(x_A,x_B^0,x_B),x_A)=0.$$(23)

Equation (23) is in fact an identity. Therefore, straightforward differentiations of this last relation with respect to x_A, $x_B^0$, and x_B lead to the following set of equations: $$\begin{array}{l}\frac{\partial \mathcal{S}}{\partial x_A^i}-\frac{\partial \mathcal{S}}{\partial x_A^0}\frac{\partial \mathcal{R}}{\partial x_A^i}=0\text{,}\\ \frac{\partial \mathcal{S}}{\partial x_B^0}-\frac{\partial \mathcal{S}}{\partial x_A^0}\left(1-\frac{\partial \mathcal{R}}{\partial x_B^0}\right)=0\text{,}\\ \frac{\partial \mathcal{S}}{\partial x_B^i}+\frac{\partial \mathcal{S}}{\partial x_A^0}\frac{\partial \mathcal{R}}{\partial x_B^i}=0\text{,}\end{array}$$

where we introduced $\mathcal{R}$ a reception range transfer function as $$\mathcal{R}(x_A,x_B)=c\mathcal{T}(x_A,x_B^0,x_B)\text{.}$$(25)

Using Eqs. (9) and (17), it is easily seen that Eq. (24) imply the following relationships: $$\begin{array}{ll}{(l_i)}_A\hfill & =\frac{\partial \mathcal{R}}{\partial x_A^i}\text{,}\hfill \\ {(l_i)}_B\hfill & =-\frac{\partial \mathcal{R}}{\partial x_B^i}{\left(1-\frac{\partial \mathcal{R}}{\partial x_B^0}\right)}^{-1}\text{,}\hfill \end{array}$$

and $$\frac{{(k_0)}_B}{{(k_0)}_A}=1-\frac{\partial \mathcal{R}}{\partial x_B^0}\text{.}$$(26c)

Unsurprisingly, these relations are similar to Eqs. (40)–(42) of Le Poncin-Lafitte et al. (2004) established for optical rays in a vacuum. We actually see that they are still valid for optical rays propagating into a linear, isotropic, and nondispersive medium within the framework of geometrical optics.

The main interest of Eq. (26) lies in the fact that it enables us to calculate the Doppler effect between an emitter and a receiver when the explicit expression of the function $\mathcal{R}$ is known. Indeed, it is well known (Synge 1960; Blanchet et al. 2001) that the Doppler frequency shift measured between an emitter and a receiver can be expressed as $$\frac{{\nu }_B}{{\nu }_A}=\frac{{(u^{\mu }{k}_{\mu })}_B}{{(u^{\mu }{k}_{\mu })}_A}=\frac{{(u^0{k}_0)}_B}{{(u^0{k}_0)}_A}\frac{{(1+{\beta }^i{l}_i)}_B}{{(1+{\beta }^i{l}_i)}_A}\text{,}$$(27)

where ${(u^{\mu })}_{A/B}$ are the unit 4-velocity vectors of the emitter and receiver defined by $${(u^{\mu })}_{A/B}={\left(\frac{\text{d}{x}^{\mu }}{\text{d}s}\right)}_{A/B}\text{,}$$(28)

with $$\text{d}{s}^2=g_{\mu \nu }\text{d}{x}^{\mu }\text{d}{x}^{\nu }\text{.}$$(29)

The unit 4-velocity is by definition a unit vector for the physical metric of spacetime Eq. (1), hence $${(u^0)}_{A/B}={\left(\frac{1}{\sqrt{1-\Vert \beta {\Vert }^2}}\right)}_{A/B}\text{,}$$(30)

where ${({\beta }^i)}_{A/B}$ denote the coordinate 3-velocity vectors of the emitter and receiver, namely $${({\beta }^i)}_{A/B}={\left(\frac{u^i}{u^0}\right)}_{A/B}=\frac{1}{c}{\left(\frac{\text{d}{x}^i}{\text{d}t}\right)}_{A/B}\text{.}$$(31)

As shown by Linet & Teyssandier (2002) and Hees et al. (2012, 2014), the frequency transfer may be expressed in terms of the time (orsimilarly the range) transfer function after inserting Eqs. (26) and (30) into (27), namely $$\frac{{\nu }_B}{{\nu }_A}=\frac{{(u^0)}_B}{{(u^0)}_A}\frac{q_B}{q_A}\text{,}$$(32)

with $$\begin{array}{ll}{q}_A\hfill & =1+{\beta }_A^i\frac{\partial \mathcal{R}}{\partial x_A^i}\text{,}\hfill \\ q_B\hfill & =1-\frac{\partial \mathcal{R}}{\partial x_B^0}-{\beta }_B^i\frac{\partial \mathcal{R}}{\partial x_B^i}\text{.}\hfill \end{array}$$

Therefore, the time and frequency transfers in Eqs. (21) and (33) can be computed once the explicit form of the time (or similarly the range) transfer function is known.

The time transfer function may be determined from the eikonal equation. Indeed, after making use of Eq. (9) and (17) while evaluating Eq. (5) at x_A, we end up with $${\left({\overline{g}}^{00}+2{\overline{g}}^{0i}{l}_i+{\overline{g}}^{ij}{l}_i{l}_j\right)}_A=0.$$(34)

Then, by invoking Eqs. (26), (6), and assumption (1), we obtain the Hamilton-Jacobi equation that is satisfied by $\mathcal{R}$. By replacing x_A by a variable x while considering $x_B^0$ and x_B as fixed parameters, the Hamilton-Jacobi equation eventually reads $$\begin{array}{ll}{[{\partial }_i\mathcal{R}{\partial }_i\mathcal{R}]}_{(x,x_B)}=\hfill & 1+{\kappa }^{00}(x_B^0-\mathcal{R}(x,x_B),x)\hfill \\ \hfill & +2{\kappa }^{0i}(x_B^0-\mathcal{R}(x,x_B),x){[{\partial }_i\mathcal{R}]}_{(x,x_B)}\hfill \\ \hfill & +{\kappa }^{ij}(x_B^0-\mathcal{R}(x,x_B),x){[{\partial }_i\mathcal{R}{\partial }_j\mathcal{R}]}_{(x,x_B)}\text{.}\hfill \end{array}$$(35)

The form of the optical metric in Eqs. (6) and (8) implies that the range transfer function can be looked for as $$\mathcal{R}(x,x_B)=\Vert x_B-x\Vert +\Delta (x,x_B)\text{,}$$(36)

where Δ is a refractive delay function. In the present context, the delay function depends on κ^μν and is due to refraction when the light ray is crossing through $\mathcal{D}$. After substituting $\mathcal{R}$ from Eqs. (36) into (35), we find the following expression: $$-2N^i{[{\partial }_i\Delta ]}_{(x,x_B)}=W(x,x_B)\text{,}$$(37)

where we introduced $$\begin{array}{ll}W(x,x_B)=\hfill & {({\kappa }^{00}-2{\kappa }^{0i}{N}^i+{\kappa }^{ij}{N}^i{N}^j)}_x\hfill \\ \hfill & +2{({\kappa }^{0i}-{\kappa }^{ij}{N}^j)}_x{[{\partial }_i\Delta ]}_{(x,x_B)}\hfill \\ \hfill & +{({\kappa }^{ij}-{\delta }^{ij})}_x{[{\partial }_i\Delta {\partial }_j\Delta ]}_{(x,x_B)}\text{,}\hfill \end{array}$$(38)

with N given by $$N=\frac{x_B-x}{\Vert x_B-x\Vert }\text{,}$$(39)

and where the point-event x is defined as $$x(x)=(x_B^0-\Vert x_B-x\Vert -\Delta (x,x_B),x)\text{.}$$(40)

As x is a free variable, we choose for convenience to consider the case where x is varying along the straight line segment joining x_A and x_B, that is x = z(σ) with $$z(\sigma )=x_B-\sigma R_{AB}{N}_{AB}\text{,}0⩽\sigma ⩽1\text{,}$$(41)

where R_AB = ∥x_B −x_A∥. In that respect, we have $$N=N_{AB}\text{.}$$(42)

A straightforward calculation shows that the total differentiation of Δ(z(σ), x_B) with respect to σ is given by $$\frac{\text{d}}{\text{d}\sigma }\Delta (z(\sigma ),x_B)=-R_{AB}{N}_{AB}^i{[{\partial }_i\Delta ]}_{(z(\sigma ),x_B)}\text{,}$$(43)

where ${[{\partial }_i\Delta ]}_{(z(\sigma ),x_B)}$ denotes the partial derivative of Δ(x, x_B) with respect to xⁱ taken at x = z(σ). Then, after inserting Eqs. (37) and (39) into (43) while accounting for Eq. (42), we infer $$\frac{\text{d}}{\text{d}\sigma }\Delta (z(\sigma ),x_B)=\frac{R_{AB}}{2}W(\tilde{z}(\sigma ),x_B)\text{,}$$(44)

where the components of the point-event $\tilde{z}(\sigma )$ are given by (40) with $\tilde{z}(\sigma )\equiv x(z(\sigma ))$, namely $$\tilde{z}(\sigma )=(x_B^0-\sigma R_{AB}-\Delta (z(\sigma ),x_B),z(\sigma ))\text{,}0⩽\sigma ⩽1.$$(45)

Equation (44) is the fundamental differential equation for the determination of the delay function, and can be integrated with the following boundary conditions $$\begin{array}{ll}\Delta (z(0),x_B)\hfill & =0\text{,}\hfill \\ \Delta (z(1),x_B)\hfill & =\Delta (x_A,x_B)\text{,}\hfill \end{array}$$

which follow from the requirement that Δ(x_B, x_B) = 0 and from z(0) = x_B. Therefore, Eq. (44) is such that $$\Delta (x_A,x_B)=\text{\&}\frac{R_{AB}}{2}{\displaystyle {\int }_{\mathcal{D}}\{}{({\kappa }^{00}-2{\kappa }^{0i}{N}_{AB}^i+{\kappa }^{ij}{N}_{AB}^i{N}_{AB}^j)}_{\tilde{z}(\sigma )}\text{nn \&}+2{({\kappa }^{0i}-{\kappa }^{ij}{N}_{AB}^j)}_{\tilde{z}(\sigma )}{\left[{\partial }_i\Delta \right]}_{(z(\sigma ),x_B)}.\text{nn \&}-{\left[{\partial }_i\Delta {\partial }_i\Delta \right]}_{(z(\sigma ),x_B)}\}\text{d}\sigma \text{,}$$(47)

where the integration goes from σ = 0 to 1 and is limited to the spacetime region within $\mathcal{D}$.

4.2 Post-Minkowskian expansion

As we plan to apply our general formalism to occultation measurements, we suppose that the refractivity of the medium satisfies the condition $$N(x)\ll 1\text{,}$$(48)

on the domain $\mathcal{D}$. This condition implies that the optical metric ḡ deviates only slightly from the Minkowski metric, namely $$\left|{\gamma }_{\mu \nu }\right|\ll 1\text{,}\left|{\kappa }^{\mu \nu }\right|\ll 1.$$(49)

It will be convenient to use the maximum value N₀ of the refractivity as a small parameter. For a rocky planet or satellite, N₀ can be defined as the refractivity at ground level. For gas giants, N₀ can be defined as the refractivity at a certain pressure level. In each case, we can put $$N(x)=N_0\mathcal{N}(x)\text{,}$$(50)

where $0<\mathcal{N}(x)⩽1$ at each point $x\in \mathcal{D}$. With this definition of $\mathcal{N}(x)$, the perturbation terms in the optical metric read $${\kappa }^{\mu \nu }(x,N_0)=N_0{\kappa }_{(1)}^{\mu \nu }(x)+{(N_0)}^2{\kappa }_{(2)}^{\mu \nu }(x)\text{,}$$(51a)

with $${\kappa }_{(1)}^{\mu \nu }(x)=2\mathcal{N}(x)w^{\mu }{w}^{\nu }\text{,}{\kappa }_{(2)}^{\mu \nu }(x)={\mathcal{N}}^2(x)w^{\mu }{w}^{\nu }\text{.}$$(51b)

Our modeling is based on the two following additional assumptions: (i) the function $\mathcal{N}(x)$ is independent of N₀, and (ii) the spatial positions of the emitter and of the receiver are such that the light ray Γ_AB, which can be observed in realistic situations, has a delay function Δ(x_A, x_B) given by the expansion $$\Delta (x_A,x_B,N_0)={\sum }_{m=1}^{\infty }{(N_0)}^m{\Delta }^{(m)}(x_A,x_B)\text{.}$$(52)

Inserting the expansions (51) and (52) into (47), one deduces the expressions for the quantities Δ^(m) as (Bourgoin 2020) $$\begin{array}{ll}{\Delta }^{(1)}(x_A,x_B)=\hfill & \frac{R_{AB}}{2}\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathcal{D}}{\left({\kappa }_{(1)}^{00}-2{\kappa }_{(1)}^{0i}{N}_{AB}^i+{\kappa }_{(1)}^{ij}{N}_{AB}^i{N}_{AB}^j\right)}_{z(\sigma )}}\text{d}\sigma \text{,}\hfill \\ {\Delta }^{(2)}(x_A,x_B)=\hfill & \frac{R_{AB}}{2}{\displaystyle {\int }_{\mathcal{D}}\{}{\left({\widehat{\kappa }}_{(2)}^{00}-2{\widehat{\kappa }}_{(2)}^{0i}{N}_{AB}^i+{\widehat{\kappa }}_{(2)}^{ij}{N}_{AB}^i{N}_{AB}^j\right)}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +2{\left({\kappa }_{(1)}^{0i}-{\kappa }_{(1)}^{ij}{N}_{AB}^j\right)}_{z(\sigma )}{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & -{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\frac{\partial {\Delta }^{(1)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\}\text{d}\sigma \text{,}\hfill \end{array}$$

and, for l ≥ 3 by $$\begin{array}{ll}{\Delta }^{(l)}(x_A,x_B)=\hfill & \frac{R_{AB}}{2}{\displaystyle {\int }_{\mathcal{D}}\{}{\left({\widehat{\kappa }}_{(l)}^{00}-2{\widehat{\kappa }}_{(l)}^{0i}{N}_{AB}^i+{\widehat{\kappa }}_{(l)}^{ij}{N}_{AB}^i{N}_{AB}^j\right)}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +2{\sum }_{m=1}^{l-1}{\left({\widehat{\kappa }}_{(m)}^{0i}-{\widehat{\kappa }}_{(m)}^{ij}{N}_{AB}^j\right)}_{(z(\sigma ),x_B)}{\left[\frac{\partial {\Delta }^{(l-m)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +{\sum }_{m=1}^{l-2}{\left({\widehat{\kappa }}_{(m)}^{ij}\right)}_{(z(\sigma ),x_B)}{\sum }_{n=1}^{l-m-1}{\left[\frac{\partial {\Delta }^{(n)}}{\partial x^i}\frac{\partial {\Delta }^{(l-m-n)}}{\partial x^j}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & -{\sum }_{m=1}^{l-1}{\left[\frac{\partial {\Delta }^{(m)}}{\partial x^i}\frac{\partial {\Delta }^{(l-m)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\}\text{d}\sigma \text{.}\hfill \end{array}$$(53c)

The integrations are limited to the refractive domain $\mathcal{D}$ (see Fig. 1) and the point-event z(σ) is defined by $$z(\sigma )=(x_B^0-\sigma R_{AB},z(\sigma ))\text{,}0⩽\sigma ⩽1.$$(54)

The quantities ${\widehat{\kappa }}_{(m)}^{\mu \nu }$ are given by $${\widehat{\kappa }}_{(1)}^{\mu \nu }(z(\sigma ),x_B)={\kappa }_{(1)}^{\mu \nu }(z(\sigma ))\text{,}$$(55a)

and, for l ≥ 2 by $$\begin{array}{ll}{\widehat{\kappa }}_{(l)}^{\mu \nu }(z(\sigma ),x_B)=\hfill & {\kappa }_{(l)}^{\mu \nu }(z(\sigma ))\hfill \\ \hfill & +{\sum }_{m=1}^{l-1}{\sum }_{n=1}^m{\Phi }^{(m,n)}(z(\sigma ),x_B){\left[\frac{{\partial }^n{\kappa }_{(l-m)}^{\mu \nu }}{{(\partial x^0)}^n}\right]}_{z(\sigma )}\text{.}\hfill \end{array}$$(55b)

The reception function Φ^(m,n)(x, x_B) is defined for m ≥ 1 and 1 ≤ n ≤ m (Teyssandier & Le Poncin-Lafitte 2008) and is given by $$\begin{array}{ll}{\Phi }^{(m,n)}(x,x_B)\hfill & =\frac{{(-1)}^n}{n!}{\sum }_{k_1+\cdots +k_n=m-n}[{\displaystyle \prod _{i=1}^n{\Delta }^{(k_i+1)}}(x,x_B)]\text{,}\hfill \end{array}$$(56)

with $k_1,\dots ,k_m\in {\mathbb{N}}_{⩾0}$. The summation in Eq. (56) is taken over all sequences of k₁ through k_n such that the sum of all k_n is m − n.

5 Application to radio occultations

The time transfer functions formalism is now successively applied to the cases of a static and a stationary optical metric. Subsequently, the discussion is specialized to spherical symmetry and to radio occultations considering the case of a downlink one-way transfer.

5.1 Static optical metric

Let us assume that the optical medium is still in the global coordinate system (x^μ), namely $$w^{\mu }=(1,0)\text{.}$$(57)

Therefore, the post-Minkowskian expansion of the optical metric can be derived straightforwardly from Eqs. (51b). The only non-null components are $${\kappa }_{(1)}^{00}=2\mathcal{N}\text{,}{\kappa }_{(2)}^{00}={\mathcal{N}}^2\text{.}$$(58)

Within such a coordinate system, the refractive properties of the medium are supposed to be independent of the component x⁰, that is to say $$n(x)-1=\{\begin{array}{ll}{N}_0\mathcal{N}(x)\hfill & \text{for}x\in \mathcal{D}\text{,}\hfill \\ 0\hfill & \text{for}x\notin \mathcal{D}\text{.}\hfill \end{array}$$(59)

Hence, the optical metric is constant and static, so that relations (55) reduce to $${\widehat{\kappa }}_{(l)}^{\mu \nu }(z(\sigma ),x_B)={\kappa }_{(l)}^{\mu \nu }(z(\sigma ))\text{for}l⩾1.$$(60)

By inserting Eqs. (58) into (53) while considering (60), we obtain the integral form of the delay function up to the lth post-Minkowskian order $$\begin{array}{ll}{\Delta }^{(1)}(x_A,x_B)=\hfill & R_{AB}{\displaystyle {\int }_{\mathcal{D}}(}\mathcal{N}{)}_{z(\sigma )}\text{d}\sigma \text{,}\hfill \\ {\Delta }^{(2)}(x_A,x_B)=\hfill & \frac{R_{AB}}{2}\hfill \\ \hfill & \times {\displaystyle {\int }_{\mathcal{D}}\{}{\left({\mathcal{N}}^2\right)}_{z(\sigma )}-{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\frac{\partial {\Delta }^{(1)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\}\text{d}\sigma \text{,}\hfill \end{array}$$

and, for l ≥ 3 $${\Delta }^{(l)}(x_A,x_B)=-\frac{R_{AB}}{2}{\displaystyle {\int }_{\mathcal{D}}}{\sum }_{m=1}^{l-1}{\left[\frac{\partial {\Delta }^{(m)}}{\partial x^i}\frac{\partial {\Delta }^{(l-m)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\text{d}\sigma \text{.}$$(61c)

As mentioned previously by Bourgoin (2020), these equations show that the first-order delay is the well-known excess path delay due to the change of the phase velocity when the signal is crossing through $\mathcal{D}$. The geometric delay shows up at the second post-Minkowskian order as well as the second-order correction to the excess path delay.

5.2 Stationary optical metric

Let us now assume that the fluid optical medium is at rest in a coordinate system rotating with respect to the global coordinate system. We use $(x^{\widehat{\alpha }})$ to denote the rotating coordinate system¹. The relationbetween the optical medium rest frame and the global coordinate system reads $$x^{\widehat{0}}=x^0\text{,}{x}^{\widehat{a}}={\Lambda }_i^{\widehat{a}}(t)x^i\text{,}$$(62a)

in which ${\Lambda }_i^{\widehat{a}}(t)\in \text{SO}(3)$ are the elements of a rotation matrix. The inverse relation is given by $$x^0=x^{\widehat{0}}\text{,}{x}^i={\Lambda }_{\widehat{a}}^i(t)x^{\widehat{a}}\text{,}$$(62b)

with ${\Lambda }_{\widehat{a}}^i(t)$ being the elements of the inverse matrix. The relationships $${\Lambda }_i^{\widehat{a}}{\Lambda }_{\widehat{b}}^i={\delta }_{\widehat{b}}^{\widehat{a}}\text{,}{\Lambda }_{\widehat{a}}^i{\Lambda }_j^{\widehat{a}}={\delta }_j^i\text{,}$$(63)

ensure that a transformation followed by its inverse is an identity transformation.

Differentiating Eqs. (62b) with respect to the global coordinate time returns the transformation law that relates the 3-velocities expressed in the two coordinate systems: $${\dot{x}}^i={\Lambda }_{\widehat{a}}^i({\dot{x}}^{\widehat{a}}+{\omega }_{\widehat{b}}^{\widehat{a}}{x}^{\widehat{b}})\text{.}$$(64)

Here, an overdot indicates differentiation with respect to t and ω_âĉ is the angular velocity tensor being defined by $${\omega }_{\widehat{c}}^{\widehat{a}}={\Lambda }_i^{\widehat{a}}{\dot{\Lambda }}_{\widehat{c}}^i\text{,}{\omega }_{\widehat{a}\widehat{c}}={\epsilon }_{\widehat{a}\widehat{b}\widehat{c}}{\omega }^{\widehat{b}}\text{.}$$(65)

The quantities ${\epsilon }_{\widehat{a}\widehat{b}\widehat{c}}$ represent the components of the permutation symbol and ω^â is the âth component of the angular velocity vector expressed in the rotating coordinate system, namely $${\omega }^{\widehat{a}}=\omega e^{\widehat{a}}\text{,}$$(66)

where ω is the magnitude of the angular velocity of rotation and e^â is the âth component of the unit direction of the spin axis.

Hereafter, we assume that the time-dependent rotation matrix corresponds to a rigid uniform rotation so that the angular velocity tensor is constant, namely $${\dot{\omega }}_{\widehat{a}\widehat{c}}=0.$$(67)

In addition, in the rotating frame, which is the rest frame of the medium, the 3-velocity vector of a particle of a fluid is null by definition, that is to say $${\dot{x}}^{\widehat{a}}=0.$$(68)

The expression for the velocity of a particle of a fluid expressed in the global coordinate system is derived after substituting ẋ^â from Eqs. (68) into (64) which leads to $$\frac{\text{d}{x}^i}{\text{d}t}={\Lambda }_{\widehat{a}}^i{\omega }_{\widehat{b}}^{\widehat{a}}{\Lambda }_j^{\widehat{b}}{x}^j\text{.}$$(69)

Let ξ(x) be a coordinate 3-velocity vector field at the point-event x belonging to a fluid element of the optical medium. It is defined in global coordinate notation by $${\xi }^i(x)=\frac{w^i}{w^0}=\frac{1}{c}\frac{\text{d}{x}^i}{\text{d}t}\text{.}$$(70)

Owing to Eqs. (65), (66), and (69), the coordinate 3-velocity vector of an element of the fluid dielectric medium is given in the global coordinate system by $$\xi (x)=\frac{\omega }{c}e\times x\text{.}$$(71)

Thus, the unit 4-velocity vector of the medium reads $$w^{\mu }(x)=w^0(1,\xi (x))\text{.}$$(72)

The expression for the time component w⁰ is straightforwardly inferred from the fact that the 4-velocity is a unit vector for the physical metric of spacetime (see assumption (1)), hence w⁰ = Γ, where Γ is defined as $$\Gamma (x)=\frac{1}{\sqrt{1-\xi (x)\cdot \xi (x)}}\text{.}$$(73)

We can now express the components of the optical metric from the relation in Eq. (6). The only non-null components are $$\begin{array}{ll}{\kappa }_{(1)}^{00}\hfill & =2{\Gamma }^2\mathcal{N}\text{,}{\kappa }_{(1)}^{0i}={\kappa }_{(1)}^{00}{\xi }^i\text{,}{\kappa }_{(1)}^{ij}={\kappa }_{(1)}^{00}{\xi }^i{\xi }^j\text{,}\hfill \\ {\kappa }_{(2)}^{00}\hfill & ={\Gamma }^2{\mathcal{N}}^2\text{,}{\kappa }_{(2)}^{0i}={\kappa }_{(2)}^{00}{\xi }^i\text{,}{\kappa }_{(2)}^{ij}={\kappa }_{(2)}^{00}{\xi }^i{\xi }^j\text{.}\hfill \end{array}$$

As seen in Sect. 3, the index of refraction is defined in the instantaneous rest frame of the medium, namely the rotating frame. Thus, the index of refraction is independent of the time component of the rotating coordinate system, that is to say $$n(x^{\widehat{a}})-1=\{\begin{array}{ll}{N}_0\mathcal{N}(x^{\widehat{a}})\hfill & \text{for}{x}^{\widehat{a}}\in \mathcal{D}\text{,}\hfill \\ 0\hfill & \text{for}{x}^{\widehat{a}}\notin \mathcal{D}\text{.}\hfill \end{array}$$(75)

This statement implies that ${\partial }_0\mathcal{N}=0$ according to the transformation rules in (62). Similarly, (71) also reveals that ∂₀ ξⁱ = 0, hence ∂₀ Γ = 0. These simplifications imply that the optical metric is stationary, and so Eqs. (55) eventually reduce to $${\widehat{\kappa }}_{(l)}^{\mu \nu }(z(\sigma ),x_B)={\kappa }_{(l)}^{\mu \nu }(z(\sigma ))\text{for}l⩾1.$$(76)

By inserting (74) into (53) while considering (76), we obtain the integral form of the delay function up to the lth post-Minkowskian order: $$\begin{array}{ll}{\Delta }^{(1)}(x_A,x_B)=\hfill & R_{AB}{\displaystyle {\int }_{\mathcal{D}}{\left({\Gamma }^2\mathcal{N}{C}^2\right)}_{z(\sigma )}}\text{d}\sigma \text{,}\hfill \\ {\Delta }^{(2)}(x_A,x_B)=\hfill & \frac{R_{AB}}{2}{\displaystyle {\int }_{\mathcal{D}}\{}{\left({\Gamma }^2{\mathcal{N}}^2{C}^2\right)}_{z(\sigma )}-{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\frac{\partial {\Delta }^{(1)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +4{\left({\Gamma }^2\mathcal{N}C{\xi }^i\right)}_{z(\sigma )}{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\}\text{d}\sigma \text{,}\hfill \\ {\Delta }^{(3)}(x_A,x_B)=\hfill & R_{AB}{\displaystyle {\int }_{\mathcal{D}}\{}{\left({\Gamma }^2\mathcal{N}{\xi }^i{\xi }^j\right)}_{z(\sigma )}{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\frac{\partial {\Delta }^{(1)}}{\partial x^j}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +{\left({\Gamma }^2{\mathcal{N}}^{2}C{\xi }^i\right)}_{z(\sigma )}{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +2{\left({\Gamma }^2\mathcal{N}C{\xi }^i\right)}_{z(\sigma )}{\left[\frac{\partial {\Delta }^{(2)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & -{\left[\frac{\partial {\Delta }^{(1)}}{\partial x^i}\frac{\partial {\Delta }^{(2)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\}\text{d}\sigma \text{,}\hfill \end{array}$$

and, for l ≥ 4, $$\begin{array}{ll}{\Delta }^{(l)}(x_A,x_B)=\hfill & \frac{R_{AB}}{2}{\displaystyle {\int }_{\mathcal{D}}\{}-{\sum }_{m=1}^{l-1}{\left[\frac{\partial {\Delta }^{(m)}}{\partial x^i}\frac{\partial {\Delta }^{(l-m)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +2{\left({\Gamma }^2\mathcal{N}{\xi }^i{\xi }^j\right)}_{z(\sigma )}{\sum }_{n=1}^{l-2}{\left[\frac{\partial {\Delta }^{(n)}}{\partial x^i}\frac{\partial {\Delta }^{(l-n-1)}}{\partial x^j}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +{\left({\Gamma }^2{\mathcal{N}}^2{\xi }^i{\xi }^j\right)}_{z(\sigma )}{\sum }_{n=1}^{l-3}{\left[\frac{\partial {\Delta }^{(n)}}{\partial x^i}\frac{\partial {\Delta }^{(l-n-2)}}{\partial x^j}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +4{\left({\Gamma }^2\mathcal{N}C{\xi }^i\right)}_{z(\sigma )}{\left[\frac{\partial {\Delta }^{(l-1)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\hfill \\ \hfill & +2{\left({\Gamma }^2{\mathcal{N}}^{2}C{\xi }^i\right)}_{z(\sigma )}{\left[\frac{\partial {\Delta }^{(l-2)}}{\partial x^i}\right]}_{(z(\sigma ),x_B)}\}\text{d}\sigma \text{.}\hfill \end{array}$$(77d)

We introduced C and D as $$C(x)=1-D(x)\text{,}D(x)=\xi (x)\cdot N_{AB}\text{.}$$(78)