---
abstract: 'When dealing with highly accurate modeling of time and frequency transfers into arbitrarily moving dielectrics medium, it may be convenient to work with Gordon’s optical spacetime metric rather than the usual physical spacetime metric. Additionally, an accurate modeling of the geodesic evolution of observable quantities (e.g. the range and the Doppler) requires us to know the reception or the emission time transfer functions. In the physical spacetime, these functions can be derived to any post-Minkowskian orders through a recursive procedure. In this work, we show that the time transfer functions can be determined to any order in Gordon’s optical spacetime as well. The exact integral forms of the gravitational, the refractive, and the coupling contributions are recursively derived. The expression of the time transfer function is given within the post-linear approximation assuming a stationary optical spacetime covered with GCRS coordinates. The light-dragging effect due to the steady rotation of the neutral atmosphere of the Earth is found to be at the threshold of visibility in many experiments involving accurate modeling of the time and frequency transfers.'
author:
- 'A. Bourgoin'
bibliography:
- 'TTF\_opt.bib'
title: General expansion of time transfer functions in optical spacetime
---

Introduction
============

In geometrical optics, the concept of light rays is introduced as curves whose tangents coincide with the direction of propagation of an electromagnetic wave [@1975ctf..book.....L]. In this approximation, refraction operates at two different levels. First, it causes the phase velocity of the electromagnetic wave to slow down or speed up while crossing a region of higher or lower refractivity, respectively. Secondly, light rays tend to bend towards region of higher refractivity. These outcomes produce an excess path delay and a geometric delay in the light time. Depending on the context, these two effects must be either thoroughly modeled or precisely measured while designing highly accurate experiments involving time and frequency transfers in presence of a refractive medium.

In many fields of astronomy such like planetary physics, astrometry, metrology, geodesy, fundamental physics, or even cosmology, we can think of situations were refractivity plays a significant role in the time and frequency transfers. For instance, we can mention that ground-based astro-geodetic techniques operating for the realization of the international terrestrial reference frame (ITRF) are currently limited by errors in modeling the group delay during the signal propagation through the Earth’s atmosphere [@2010ITN....36....1P; @2004GeoRL..3114602M; @2002GeoRL..29.1414M; @1997JGR...10220489C; @doi1010292006JB004834]. We can also mention the cases of atmospheric radio occultations [@1965Sci...149.1243K; @1985AJ.....90.1136L; @1987JGR....9214987L; @1992AJ....103..967L; @2012Icar..221.1020S; @2015RaSc...50..712S] and atmospheric stellar ocultation experiments . Indeed, both techniques aim at determining a refractivity profile towards an occulting atmosphere from precise measurements of an *a priori* known frequency (usually given in the frame at rest with the emitter) and from an accurate modeling of the frequency transfer in presence of the occulting refractive medium. To an even higher degree of accuracy, we can cite experiments involving frequency transfers between distant atomic clocks via ground-ground free-space optical (FSO) link [@Djerroud.10; @2013NaPho...7..434G; @doi.10.1063.1.4963130], space-ground FSO link [@Fujieda_2014; @Hachisu.14], and optical fiber links [@2010Sci...329.1630C; @Predehl441; @2015Metro..52..552G]. Finally, let us emphasize that in the context of cosmology, it has been shown that the accumulated effect of an artificial refractivity over the distance-redshift relation fits perfectly Hubble curve of type Ia supernovae data in the framework of a non-accelerating cosmological model [@PhysRevD.78.044040; @*PhysRevD.79.104007; @*PhysRevD.80.044019]. All these examples highlight how important refraction can be in highly accurate experiments involving time and frequency transfers.

In the past, two independent theoretical formalism have been introduced, namely Gordon’s optical metric and the time transfer functions. On one side, Gordon’s metric allows one to handle refraction in curved spacetime; on the other side, the time transfer functions formalism handles theoretical problems related to the time and frequency transfers in curved spacetime. In this work, we intend to combine the two formalisms which are discussed in turn in the next.

In the early 20’s, Gordon introduced [@doi101002andp19233772202] a useful theoretical tool to study light refraction caused by an arbitrarily moving fluid dielectric medium, namely Gordon’s optical metric. In this work, he showed that in the presence of a fluid whose electromagnetic properties are described by a permittivity $\epsilon(x)$ and a permeability $\mu(x)$, any solutions to the macroscopic Maxwell’s equations can be looked for indifferently either in the usual physical spacetime fitted with the metric tensor, or in an artificial optical spacetime fitted with Gordon’s metric. Conveniently, in the optical spacetime and within the geometric optics approximation, by means of a slightly different set of Maxwell’s equations, the electromagnetic properties of the fluid medium are reduced to their vacuum values, that is to say $\epsilon(x)=\epsilon_0$ and $\mu(x)=\mu_0$. In other words, in the physical spacetime, the interaction between the electromagnetic field and the dielectric fluid medium must be carefully modeled, whereas in the optical spacetime this interaction is implicitly accounted for in the vacuum limit of the macroscopic version of Maxwell’s equations. Consequently, within the geometric optics approximation, light rays propagate into dielectrics medium along null geodesic lines of the optical spacetime.

At the same time, theoretical problems dealing with the deflection of light rays or the frequency transfer, require us to know the function relating the (coordinate) time transfer to the coordinate time at the reception and to the spatial coordinates of the reception and the emission points. Such function is called a reception time transfer function. Obliviously, an emission time transfer function can be introduced as well. The formalism which aims at determining the time transfer functions was first introduced by [@2002PhRvD..66b4045L] relying on the theory of the world function developed by Synge [@SyngeBookGR]. General expansions of the world function and the time transfer functions was first proposed by [@2004CQGra..21.4463L], and then, a simplified recursive approach, based on the determination of time delay functions instead of Synge’s world function, was presented in [@2008CQGra..25n5020T]. The usefulness of the time transfer function formalism lies in the fact that it spares the trouble of solving explicitly the null geodesic equation which usually leads to heavy calculations beyond the post-Minkowskian regime (see e.g. for explicit resolution of the null geodesic equation in the linearized weak field limit, and see e.g. [@1983PhRvD..28.3007R; @1987KFNT....3....8B; @Linet_2013] for resolution in the post-post Minkowskian approximation). Indeed, assuming that the emission and reception points-events are linked by a null geodesic path (quasi-Minkowskian path approximation), the time transfer functions formalism achieves a complete resolution of the time and frequency transfers to any post-Minkowskian order by means of an algorithmic resolution method [@2008CQGra..25n5020T]. For this reason, this formalism is currently one of the most powerful theoretical tool to derive the time and frequency transfers along null geodesics of the curved physical spacetime.

The scope of this paper is to generalize the formalism of the time transfer functions to optical spacetime. The aim is to provide a recursive method allowing one to solve theoretical problems related to the propagation of light in presence of arbitrarily moving refractive medium.

The present work is organized as follows. In Sec. \[sec:not\], we present the notations and conventions used throughout this paper. Sec. \[sec:geO\] is a short reminder about the use of Gordon’s metric in relativistic geometrical optics. In this section, we derive the optical counterpart of the scalar Eikonal equation (fundamental equation of geometrical optics) which is at the basis of the demonstration which follows. Sec. \[sec:TF\] is a recall about the time transfer functions formalism. In Sec. \[sec:intTF\], by applying a method initially proposed by [@2008CQGra..25n5020T], we show that working in optical spacetime induces that the time transfer functions can be decomposed in three components, that we call the gravitational, the refractive, and the coupling time transfer functions. In Sec. \[sec:exp\], we present the general expansion of the three contributions. In Sec. \[sec:app\], we illustrate the method by computing the time transfer function of an optical spacetime describing Earth’s rotating atmosphere in the geocentric celestial reference system (GCRS) within the post-linear approximation. Finally, we discuss about the importance of taking into account the light-dragging effect in future generation of data reduction softwares.

Notations and conventions {#sec:not}
=========================

In this work, the metric of spacetime is denoted $g$ and its signature is $(+,-,-,-)$. The optical metric (also called Gordon’s metric) is denoted $\bar g$. We suppose that spacetime is covered with some global coordinates system $\bm x=(x^0,\mathbf{x})$. We put $x^0=ct$ with $c$ being the speed of light in a vacuum and $t$ being the coordinate time.

Greek indices run from 0 to 3 and Latin indices run from 1 to 3. Straight bold letters denote ordered triples such as $\mathbf x=(x^i)$ and italic bold letters denote ordered quadruples such like $\bm x=(x^\mu)$. Einstein’s summation convention on repeated indices is used for expressions like $a^ib^i$ as well as for expressions like $A^\mu B_\mu$. The ordinary Euclidean norm of $\mathbf{x}$ is denoted $|\mathbf x|$ and is defined as $|\mathbf{x}|=(\delta_{ij}x^ix^j)^{1/2}$ where $\delta_{ij}$ is the Kronecker delta. The maximum absolute value of the component $A_{\mu\nu}$ is denoted $|A_{\mu\nu}|_{\mathrm{max}}$. The 3-dimensional antisymmetric Levi-Civita tensor is denoted $e^{ijk}$.

For the sake of legibility, we would 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 would employ $(f)_{A/B}$ to denote $f(x_A)$ and $f(x_B)$, respectively. The partial differentiation with respect to coordinates $x^\mu$ is denoted $\partial_\mu A_\nu$. The physical and the optical covariant differentiation with respect to $x^\mu$ are denoted $\nabla_\mu A_\nu$ and $\bar\nabla_\mu A_\nu$, respectively. Given a scalar function $f(\bm x)$, we have the relation $\bar\nabla_\mu f=\nabla_\mu f=\partial_\mu f$.

Throughout the paper, we assume the presence of an arbitrarily moving fluid dielectric medium filling a finite domain $\mathcal{D}$ of spacetime. We call $\bm w(x)$ the unit 4-velocity vector of a point-event $x$ belonging to a fluid element of the optical medium. The expression of $\bm w(x)$ is given by $$\bm w(x)\equiv\frac{{\mathrm{d}}\bm x}{{\mathrm{d}}s}
  \label{eq:4velel}$$ where the spacetime interval ${\mathrm{d}}s$ is defined by $${\mathrm{d}}s^2=g_{\mu\nu}(x){\mathrm{d}}x^{\mu} {\mathrm{d}}x^{\nu}\text{.}$$

We call $\xi^i(x)$ the coordinate 3-velocity vector of the point-event $x$ belonging to a fluid element of the optical medium. Its expression is given by $$\xi^i(x)\equiv\frac{w^i}{w^0}=\frac{1}{c}\frac{{\mathrm{d}}x^i}{{\mathrm{d}}t}\text{.}
  \label{eq:3velel}$$

Finally, $G$ is the Newtonian gravitational constant.

Relativistic Geometrical optics {#sec:geO}
===============================

We assume the presence of an optical medium filling $\mathcal{D}$. Additionally, we consider for simplicity that the fluid’s electromagnetic properties are linear, isotropic, non-dispersive and can be summarized by two scalar functions, namely the permittivity $\epsilon(x)$ and the permeability $\mu(x)$. These two quantities completely determine the refractive properties of the optical medium through the following relationship $$n(x)\equiv c\sqrt{\epsilon(x)\mu(x)}
  \label{eq:n}$$ where $n$ is the index of refraction of the medium.

When $x\notin\mathcal{D}$, the permittivity and the permeability reduce to their vacuum values $\epsilon(x)=\epsilon_0$ and $\mu(x)=\mu_0$, respectively. Thus, considering that $c\equiv(\epsilon_0\mu_0)^{-1/2}$, the index of refraction becomes $n(x)=1$. By subtracting its vacuum value from the index of refraction, we obtain the refractivity $$N(x)\equiv n(x)-1\text{,}
  \label{eq:N}$$ which is obviously null in a vacuum.

In the physical spacetime, the evolution of electromagnetic phenomenon occurring in the presence of an optical medium are usually described by the macroscopic version of Maxwell’s equations. These equations are separated into two distinct sets involving a covariant antisymmetric tensor $F_{\mu\nu}$ called the electromagnetic field tensor (or Faraday tensor), and a contravariant antisymmetric tensor $B^{\mu\nu}$ called the electromagnetic field excitation tensor (or Maxwell tensor), respectively. The macroscopic version of Maxwell’s equations are given by [@1960ecm..book.....L; @doi101002andp200810313]

\[eq:Max\] $$\begin{aligned}
  \partial_{[\sigma} F_{\mu\nu]}&=0\text{,}\label{eq:Max1}\\
  \nabla_\mu B^{\mu\nu}&=j^\nu\label{eq:Max2}\end{aligned}$$

where $\bm j(x)$ is a 4-vector denoting the free charge density-current. The square brackets denote the complete antisymmetrization of the enclosed indices.

The first set , allows one to postulate the existence of a covector field $A_\mu(x)$, such that the electromagnetic field tensor $F_{\mu\nu}$ can be locally written as the rotational of the covector field, that is $$F_{\mu\nu}=\mathrm{Re}\big\{\partial_{\mu}A_\nu-\partial_{\nu}A_\mu\big\}\text{.}
  \label{eq:F}$$

The second set cannot be used alone to fully determine the six independent components of the electromagnetic field excitation tensor $B^{\mu\nu}$. In addition, it does not provide a way to determine the components of the electromagnetic field tensor $F_{\mu\nu}$ which yet governs the motion of particles through the Lorentz force. Therefore, Maxwell’s equations must be supplemented with constitutive relations.

For an arbitrarily moving medium of permittivity $\epsilon(x)$ and permeability $\mu(x)$ the covariant constitutive relationships are given by [@1960ecm..book.....L]

\[eq:HFctv\] $$\begin{aligned}
  B^{\mu\nu}w_\nu&=\epsilon c^2F^{\mu\nu}w_\nu\text{,}\label{eq:Hctv}\\
  \mu B_{[\mu\nu}w_{\sigma]}&=F_{[\mu\nu}w_{\sigma]}\text{.}\label{eq:Fctv}\end{aligned}$$

Eqs. can be written as a single relationship involving $B^{\mu\nu}$, $F_{\mu\nu}$, and $\bm w(x)$. Indeed, as initially shown by Gordon [@doi101002andp19233772202], when dealing with problems of electromagnetic waves propagating into dielectrics, it is convenient to introduce an optical spacetime in which refractivity is considered as spacetime curvature. The Gordon’s metric (or optical metric) is defined by

\[eq:Gor\] $$\bar g_{\mu\nu}\equiv g_{\mu\nu}+\gamma_{\mu\nu}\text{,} \qquad \gamma_{\mu\nu}=-\left(1-\frac{1}{n^2}\right)w_\mu w_\nu\text{,}\label{eq:Gorcov}$$ with inverse $$\bar g^{\mu\nu}\equiv g^{\mu\nu}+\kappa^{\mu\nu}\text{,} \qquad \kappa^{\mu\nu}=(n^2-1)w^\mu w^\nu\text{.}\label{eq:Gorcon}$$

Making use of Eq. , one can see that Eqs. are summarized within the single following relation [@doi101002andp19233772202] $$\mu B^{\mu\nu}=\bar F^{\mu\nu}
  \label{eq:ctv}$$ where the optical metric has been used to raise covariant indices of $F_{\alpha\beta}$, that is $$\bar F^{\mu\nu}\equiv\bar g^{\mu\alpha}\bar g^{\nu\beta} F_{\alpha\beta}\text{.}
  \label{eq:bF}$$

It is now possible to express Maxwell’s equation in the optical spacetime. Because, the covariant components of the electromagnetic field tensor are equivalents in both spacetimes [^1], the first pair of Maxwell’s equations remains unchanged. The optical form of the second pair is obtained after substituting for $B^{\mu\nu}$ from while introducing the optical covariant derivative [@PhysRevD.78.044040; @PhysRevD.79.104007]. After little algebra, we find $$\bar{\nabla}_\mu\left(\sqrt{\frac{\epsilon}{\mu}}\,\bar F^{\mu\nu}\right)=\sqrt{\epsilon\mu}\,j^{\nu}\text{.}
  \label{eq:bMax}$$

Eq. is perfectly equivalent to Eq. equipped with the constitutive relations . While working in the optical spacetime, Eq. allows to find $\bar F^{\mu\nu}$ and Eq.  lets to express the components of the electromagnetic field tensor in the physical spacetime. Hereafter, we work in the optical spacetime where the light propagation into the dielectric medium is simply given by the vacuum limit of the macroscopic version of Maxwell’s equations (no free density current i.e. $j^\nu=0$).

In this work, we consider geometrical optics approximation, so we assume that the 4-potential covector $A_\mu(x)$ of a traveling quasi-monochromatic wave possesses an expansion of the form [@1975ctf..book.....L] $$A_\mu=\big[a_\mu+\mathcal{O}(\omega^{-1})\big]e^{i\omega\mathscr{S}}\text{.}
  \label{eq:pot}$$ Here, $\mathscr{S}(x)$ is the usual eikonal function which determines the surfaces of constant phase for the wave, $a_\mu(x)$ is the complex covector amplitude varying slowly in comparison to $\mathscr{S}(x)$, and $\omega$ is a book-keeping parameter that we take to be high during our manipulations [@gravitationBook].

Then, substituting for $A_\mu$ from Eq. into , allows one to infer $$F_{\mu\nu}=\mathrm{Re}\Big\{\big[i\omega f_{\mu\nu}+\mathcal{O}\left(\omega^0\right)\big]e^{i\omega\mathscr{S}}\Big\}\text{,}
  \label{eq:Fopt}$$ where $f_{\mu\nu}(x)$ represent the electromagnetic field tensor amplitude’s components, that is $$f_{\mu\nu}=k_\mu a_\nu-k_\nu a_\mu\text{,}
  \label{eq:A}$$ with $k_\mu$ being the wave covector defined by $$k_\mu\equiv\partial_\mu\mathscr{S}\text{.}
  \label{eq:kdef}$$

We can introduce the contravariant optical wave vector such that $$\bar k^\mu\equiv\bar g^{\mu\nu}k_\nu
  \label{eq:bkdef}$$ where the low index has been raised with the help of the optical spacetime metric. We can directly check from the inverse conditions $\bar g_{\mu\sigma}\bar g^{\sigma\nu}=\delta_\mu^\nu$, that the covariant components of the wave vector are identical in physical and optical spacetimes, that is to say $$\bar k_\mu =k_\mu\text{.}
  \label{eq:kbk}$$

Assuming that the 4-potential fulfills the Lorentz gauge in the optical spacetime, that is $$\bar\nabla_\mu \bar A^\mu=0
  \label{eq:bLG}$$ where we introduced $\bar A^\mu\equiv\bar g^{\mu\nu}A_\nu$, we find $$\bar g^{\mu\nu}k_\mu a_\nu=0
  \label{eq:bLGopt}$$ within the geometrical optics approximation. This relationship states the orthogonality between the optical wave vector $\bar{\bm k}$ and the wave covector amplitude $a_\mu$.

Finally, the fundamental equations of geometrical optics can be derived from the vacuum limit of the optical version of Maxwell’s equations. We first determine the optical electromagnetic field tensor by making use of Eqs.  and . Then, by taking the covariant derivative of $\bar F^{\mu\nu}$, we find $$\bar\nabla_\mu\bar F^{\mu\nu}=-\mathrm{Re}\Big\{\big[\omega^2k_\mu\bar f^{\mu\nu}+\mathcal{O}(\omega)\big]e^{i\omega\mathscr{S}}\Big\}
  \label{eq:CDFopt}$$ where we introduced $\bar f^{\mu\nu}\equiv\bar g^{\mu\alpha}\bar g^{\nu\beta}f_{\alpha\beta}$. By substituting this result into the vacuum limit of Eq. , and by keeping the geometrical optics order only, one deduces $$\mathrm{Re}\big\{k_\mu\bar f^{\mu\nu}\big\}=0\text{.}
  \label{eq:beikA}$$

Then, substituting for $\bar f^{\mu\nu}$ from Eq. into and considering , we finally deduce $$\bar g^{\mu\nu}k_\mu k_\nu=0\text{.}
  \label{eq:beik}$$ This is the fundamental equation of geometrical optics expressed in optical spacetime. After inserting Eq. , we infer that the phase $\mathscr{S}(x)$ satisfies the well-known scalar Eikonal equation $$\bar g^{\mu\nu}\partial_\mu\mathscr{S}\partial_\nu\mathscr{S}=0\text{.}
  \label{eq:beikS}$$

We close this section by showing that $\bar{\bm k}$ is a null vector satisfying the geodesic equation for the optical metric. From Eqs. and , we easily infer $$\bar g_{\mu\nu}\bar k^\mu \bar k^\nu=0\text{.}
  \label{eq:beikcont}$$ This relation shows that $\bar{\bm k}$ is indeed isotropic for the optical metric $\bar g_{\mu\nu}$. Then, we differentiate Eq. with respect to $x^\sigma$. Considering the symmetry of the components of the optical metric together with Eq. , it comes $$\bar k^\nu(\bar\nabla_\sigma k_\nu)=0\text{.}$$ Making use of the definition , we infer $\bar\nabla_\sigma k_\nu=\bar\nabla_\nu k_\sigma$. Finally, considering Eqs. and , we deduce $$(\bar k^\nu\bar\nabla_\nu) \bar k^\sigma=0\text{,}
  \label{eq:geod}$$ which states (together with Eq. ) that curves admitting $\bar{\bm k}$ as tangent vector are null geodesic lines of the optical metric. In that respect, a null line which is solution of Eq. can be interpreted as a ray of light whose tangent at any point $x$ is orthogonal to the surface of constant phase $\mathscr{S}(x)$ [@SyngeBookGR].

Time transfer formalism {#sec:TF}
=======================

Let us consider a light ray $\Gamma$ propagating in a region of spacetime covered with some coordinate system $(x^\mu)$. Let $\bm x_A=(ct_A,\mathbf{x}_A)$ be the coordinates of a point-event $x_A$. We introduce $\mathcal{C}_A$ as the curve of parametric equations $x=x_A(\tau)$ with $\tau$ being a parametrization of $\mathcal{C}_A$. Let us suppose that $(x^\mu)$ is chosen in such away that $\mathcal{C}_A$ is a time-like worldline for any $x_A$, which means that $\partial/\partial x^0$ is a time-like vector field, that is to say $g_{00}>0$ everywhere. Let $x_A$ be the point-event where $\Gamma$ is emitted and let $x_B$ be the point-event of coordinates $\bm x_B=(ct_B,\mathbf{x}_B)$ where it is observed. The quantity $t_B-t_A$ is the (coordinate) travel time of the light ray connecting the emission point-event $x_A$ and the reception point-event $x_B$. This quantity allows us to introduce the time transfer functions $\mathcal{T}_{r,\Gamma}$ and $\mathcal{T}_{e,\Gamma}$ [@2004CQGra..21.4463L] as $$t_B-t_A\equiv\mathcal{T}_{r,\Gamma}(\mathbf x_A,t_B,\mathbf x_B)\equiv\mathcal{T}_{e,\Gamma}(t_A,\mathbf x_A,\mathbf x_B)\text{.}
  \label{eq:TTFdefGen}$$ We call $\mathcal{T}_{r,\Gamma}$ the reception time transfer function and $\mathcal{T}_{e,\Gamma}$ the emission time transfer function associated to $\Gamma$.

As shown in [@PhysRevD.93.044028], given a point-event $x_B$ and a spatial position $\mathbf x_A$, $\Gamma$ is not unique in general. Let $\Gamma^{[\sigma]}(\mathbf x_A,x_B)$ be a light ray intersecting $x_B$ and flowing from the point-event $x_A^{[\sigma]}\in\mathcal{C}_A$ of coordinates $(ct_A^{[\sigma]},\mathbf x_A)$. For each $\Gamma^{[\sigma]}$, there exists a reception time transfer function, denoted $\mathcal{T}_{r,\Gamma^{[\sigma]}}(\mathbf x_A,t_B,\mathbf{x}_B)$, such that $$t_B-t_A^{[\sigma]}=\mathcal{T}_{r,\Gamma^{[\sigma]}}(\mathbf x_A,t_B,\mathbf{x}_B)\text{,}$$ with $\sigma\in\mathbb{N}_{>0}$ (the same reasoning works for the emission time transfer function as well).

This fact shows that, in general, we cannot expect to find a unique reception (or emission) time transfer function. However, for a very particular type of null geodesics, referred to as quasi-Minkowskians [@2012CQGra..29x5010T; @PhysRevD.93.044028], it has been shown that the reception (or the emission) time transfer function, if it exists, can be uniquely determined [@2008CQGra..25n5020T].

Henceforth, we assume that $\Gamma$ is a quasi-Minkowskian light ray so that the corresponding time transfer functions are indeed unique. In agreement with this assumptions, we suppose that the past null cone at $x_B$ intersects $\mathcal{C}_A$ at one and only one point $x_A$. Therefore, we can rewrite Eq.  as $$t_B-t_A\equiv\mathcal{T}_r(\mathbf x_A,t_B,\mathbf x_B)\equiv\mathcal{T}_e(t_A,\mathbf x_A,\mathbf x_B)\text{.}
  \label{eq:TTFdef}$$

Hereafter, in order to shorten future notations, we introduce the reception and the emission range transfer functions being defined by

\[eq:RTFdef\] $$\mathcal{R}_r(\mathbf x_A,x_B)\equiv c\mathcal{T}_r(\mathbf x_A,t_B,\mathbf x_B)\text{,}
  \label{eq:RTFrdef}$$ and $$\mathcal{R}_e(x_A,\mathbf x_B)\equiv c\mathcal{T}_e(t_A,\mathbf x_A,\mathbf x_B)\text{.}
  \label{eq:RTFedef}$$

An important theorem (cf. Theorem 1 of [@2004CQGra..21.4463L]) states that the covariant components of the tangent vector are totally known, as soon as one of the time transfer functions (or equivalently, one of the range transfer functions) is explicitly determined. Therefore, if we define $$(l_i)_{A/B}\equiv \left(\frac{k_i}{k_0}\right)_{A/B}\text{,}
  \label{eq:l}$$ we have the following relationships

\[eq:kTFdef\] $$\begin{aligned}
  (l_i)_A&\equiv\frac{\partial\mathcal{R}_r}{\partial x_A^i}=\frac{\partial\mathcal{R}_e}{\partial x_A^i}\left(1+\frac{\partial\mathcal{R}_e}{\partial x_A^0}\right)^{-1}\text{,}\label{eq:kTFA}\\
  (l_i)_B&\equiv-\frac{\partial\mathcal{R}_r}{\partial x_B^i}\left(1-\frac{\partial\mathcal{R}_r}{\partial x_B^0}\right)^{-1}=-\frac{\partial\mathcal{R}_e}{\partial x_B^i}\text{,}\label{eq:kTFB}\end{aligned}$$ and $$\frac{(k_0)_B}{(k_0)_A}\equiv1-\frac{\partial\mathcal{R}_r}{\partial x_B^0}=\left(1+\frac{\partial\mathcal{R}_e}{\partial x_A^0}\right)^{-1}\text{.}
  \label{eq:kTF0}$$

Consequently, Eqs. completely solve theoretical problems related to frequency transfer. Indeed, it is well-known that the instantaneous expression of the Doppler shift along the null-geodesic path between the emitter and the receiver can be expressed as [@SyngeBookGR] $$\frac{\nu_B}{\nu_A}\equiv\frac{(u^{\mu}k_\mu)_B}{(u^{\mu}k_\mu)_A}=\frac{(u^0k_0)_B}{(u^0k_0)_A}\frac{(1+\beta^il_i)_B}{(1+\beta^il_i)_A}
  \label{eq:dop}$$ where $(\bm u)_{A/B}$ is the emitter/receiver’s unit 4-velocity vectors being defined as $$(\bm u)_{A/B}\equiv\left(\frac{{\mathrm{d}}\bm x}{{\mathrm{d}}s}\right)_{A/B}\text{,}
  \label{eq:4velAB}$$\
with $({\mathrm{d}}s^2)_{A/B}=(g_{\mu\nu}{\mathrm{d}}x^{\mu} {\mathrm{d}}x^{\nu})_{A/B}$.

By definition, the 4-velocities satisfy the unity condition $(g_{\mu\nu}u^\mu u^\nu)_{A/B}=1$, which implies $$(u^0)_{A/B}=(g_{00}+2g_{0i}\beta^i+g_{ij}\beta^i\beta^j)^{-1/2}_{A/B}\text{.}
  \label{eq:4vel}$$ The quantities $(\beta^i)_{A/B}$ in Eq. are the components of the emitter/receiver’s coordinate 3-velocity vectors and are defined such that $$(\beta^i)_{A/B}\equiv\left(\frac{u^i}{u^0}\right)_{A/B}=\frac{1}{c}\left(\frac{{\mathrm{d}}x^i}{{\mathrm{d}}t}\right)_{A/B}\text{.}
  \label{eq:3velAB}$$

It is then straightforward to determine the exact expression of the instantaneous Doppler formulation in terms of the range transfer functions . Indeed, after inserting Eqs.  and within , we infer $$\frac{\nu_B}{\nu_A}=\frac{(u^0)_B}{(u^0)_A}\frac{q_B}{q_A}\text{,}
  \label{eq:dopTF}$$ with

\[eq:qTF\] $$\begin{aligned}
  q_A&=1+\beta^i_A\frac{\partial\mathcal{R}_r}{\partial x_A^i}=1+\frac{\partial\mathcal{R}_e}{\partial x_A^0}+\beta^i_A\frac{\partial\mathcal{R}_e}{\partial x_A^i}\text{,}\label{eq:qA}\\
  q_B&=1-\frac{\partial\mathcal{R}_r}{\partial x_B^0}-\beta^i_B\frac{\partial\mathcal{R}_r}{\partial x_B^i}=1-\beta^i_B\frac{\partial\mathcal{R}_e}{\partial x_B^i}\text{,}\label{eq:qB}\end{aligned}$$

and $$\frac{(u^0)_B}{(u^0)_A}=\frac{(g_{00}+2g_{0i}\beta^i+g_{ij}\beta^i\beta^j)^{1/2}_{A}}{(g_{00}+2g_{0i}\beta^i+g_{ij}\beta^i\beta^j)^{1/2}_{B}}\text{.}
  \label{eq:uBuA}$$

From the fundamental equation of geometrical optics (see Eq. ), we know that the covariant components of the 4-wave optical vector at point-events $x_A$ or $x_B$ satisfy a relation as follows $$(\bar g^{\mu\nu}k_\mu k_\nu)_{A/B}=0\text{.}
  \label{eq:beikAB}$$ Then, dividing by $[(k_0)_{A/B}]^2$ and making use of Eqs. -, we infer the following theorem which generalizes theorem 1 of [@2008CQGra..25n5020T] to optical spacetime.

Within geometrical optics approximation, the range transfer functions $\mathcal{R}_r$ and $\mathcal{R}_e$ satisfy the following Hamilton-Jacobi-like equations over the optical spacetime, namely

\[eq:beikHJ\] $$\bar g^{00}(x_B^0-\mathcal{R}_r,\mathbf{x}_A)+2\bar g^{0i}(x_B^0-\mathcal{R}_r,\mathbf{x}_A)\frac{\partial\mathcal{R}_r}{\partial x_A^i}+\bar g^{ij}(x_B^0-\mathcal{R}_r,\mathbf{x}_A)\frac{\partial\mathcal{R}_r}{\partial x_A^i}\frac{\partial\mathcal{R}_r}{\partial x_A^j}=0\text{,}\label{eq:beikHJA}$$ and $$\bar g^{00}(x_A^0+\mathcal{R}_e,\mathbf{x}_B)-2\bar g^{0i}(x_A^0+\mathcal{R}_e,\mathbf{x}_B)\frac{\partial\mathcal{R}_e}{\partial x_B^i}+\bar g^{ij}(x_A^0+\mathcal{R}_e,\mathbf{x}_B)\frac{\partial\mathcal{R}_e}{\partial x_B^i}\frac{\partial\mathcal{R}_e}{\partial x_B^j}=0\text{,}\label{eq:beikHJB}$$

respectively. \[th:theo1\]

This theorem is at the basis of the demonstration for deriving the integral form of the range and then the time transfer functions. Henceforth, in order to avoid repetitions, we pursue the demonstration giving details only for the reception time delay function. However, the same results can be derived for the emission time delay function by applying the exact same reasoning.

Integral form of the time delay functions {#sec:intTF}
=========================================

Now, let us assume that the physical spacetime metric takes the following form

\[eq:gm\] $$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}
  \label{eq:gmcov}$$ throughout spacetime, where $\eta_{\mu\nu}$ is the Minkowski metric and $h_{\mu\nu}$ is the gravitational perturbation. In Cartesian coordinates, $\eta_{\mu\nu}=\text{diag}(+1,-1,-1,-1)$. The contravariant components of the physical spacetime metric can be decomposed as $$g^{\mu\nu}=\eta^{\mu\nu}+k^{\mu\nu}
  \label{eq:gmcon}$$

where the components $k^{\mu\nu}$ satisfy $$k^{\mu\nu}=-\eta^{\mu\alpha}\eta^{\beta\nu}h_{\alpha\beta}-\eta^{\mu\alpha}h_{\alpha\beta}k^{\beta\nu}\text{.}
  \label{eq:kmcon}$$

Therefore, the optical spacetime metric can be expressed as

\[eq:Gordec\] $$\bar g_{\mu\nu}=\eta_{\mu\nu}+H_{\mu\nu}\text{,}
  \label{eq:Gordeccov}$$ with the contravariant components $$\begin{aligned}
  \bar g^{\mu\nu}&=\eta^{\mu\nu}+K^{\mu\nu}\text{.}
  \label{eq:Gordeccon}\end{aligned}$$

Thus, the optical metric reduces to the sum of the flat Minkowski metric plus a spacetime curvature contribution which is given by

\[eq:Cur\] $$H_{\mu\nu}=h_{\mu\nu}+\gamma_{\mu\nu}\text{,}\label{eq:Curcov}$$ with the contravariant components $$K^{\mu\nu}=k^{\mu\nu}+\kappa^{\mu\nu}\text{.}\label{eq:Curcon}$$

From here, we suppose that the curvature contribution is small so that spacetime is mainly flat, that is to say $|h_{\mu\nu}|_{\mathrm{max}}\ll|\eta_{\mu\nu}|_{\mathrm{max}}$ and $|\gamma_{\mu\nu}|_{\mathrm{max}}\ll|\eta_{\mu\nu}|_{\mathrm{max}}$. In other words, we focus on the post-Minkowskian approximation. Under this condition, we ensure that the null geodesic path is quasi-Minkowskian. The form of the optical metric in Eqs. implies that the reception and the emission range transfer functions can be looked for according to the following expressions

\[eq:RTFdec\] $$\mathcal{R}_r(\mathbf{x}_A,x_B)=|\mathbf{x}_B-\mathbf{x}_A|+\Delta(\mathbf{x}_A,x_B)\text{,}\label{eq:RTFdecR}$$ and $$\mathcal{R}_e(x_A,\mathbf{x}_B)=|\mathbf{x}_B-\mathbf{x}_A|+\Xi(x_A,\mathbf{x}_B)\text{,}\label{eq:RTFdecE}$$

respectively. Following [@2008CQGra..25n5020T], we will call $\Delta/c$ the reception time delay function and $\Xi/c$ the emission time delay function [^2].

Now, if we assume that the reception point-event $x_B$ is perfectly known, then, we can regard its components $t_B$ and $\mathbf x_B$ as fixed parameters. Hence, the reception time delay function becomes a function of the spatial components of the emission point-event $x_A$ [^3]. Thus, if we now substitute $\mathbf x$ to $\mathbf x_A$, the reception time delay function $\Delta(\mathbf x,x_B)/c$ uniquely defines the point-event $x_-(\mathbf x)$ for the given set of spatial components $\mathbf{x}$, that is $$\bm x_-(\mathbf x)=\big(x_B^0-|\mathbf x_B-\mathbf x|-\Delta(\mathbf{x},x_B),\mathbf{x}\big)\text{.}
  \label{eq:x-}$$ Furthermore, assuming that the point-event $x_-$ lies in the vicinity of $x_B$, we can determine the spatial variation of the reception time delay function. Indeed, after inserting Eqs.  and into taken at $x_-$ instead of $x_A$, we deduce the following relation [@2008CQGra..25n5020T] $$-2N^i\partial_i\Delta(\mathbf{x},x_B)=\Omega_-(x_-,x_B)
  \label{eq:beikDP}$$ where $\mathbf{N}=(\mathbf{x}_B-\mathbf{x})/|\mathbf{x}_B-\mathbf{x}|$, and $$\begin{aligned}
  \Omega_-(x_-&,x_B)=\big(K^{00}-2K^{0i}N^i+K^{ij}N^iN^j\big)_{x_-}\nonumber\\
  &+2\big(K^{0i}-K^{ij}N^j\big)_{x_-}\partial_i\Delta(\mathbf{x},x_B)\nonumber\\
  &+\big(\eta^{ij}+K^{ij}\big)_{x_-}\partial_i\Delta(\mathbf{x},x_B)\partial_j\Delta(\mathbf{x},x_B)\text{.}\label{eq:Om}\end{aligned}$$

Since $\mathbf x$ is a free variable, we follow [@2008CQGra..25n5020T] and choose for convenience to focus on the case where $\mathbf x$ is varying along the straight line segment connecting $\mathbf x_A$ to $\mathbf x_B$, that is to say $\mathbf{x}=\mathbf{z}_-(\lambda)$ where $$\begin{aligned}
  \mathbf{z}_-(\lambda)&=\mathbf{x}_B-\lambda R_{AB}\mathbf{N}_{AB}\text{,} & 0\leqslant\lambda\leqslant 1\text{,}\label{eq:z-i}\end{aligned}$$ with $R_{AB}=|\mathbf{x}_B-\mathbf{x}_A|$ and $\mathbf N_{AB}=(\mathbf{x}_{B}-\mathbf{x}_{A})/R_{AB}$. In that respect, we also have the relation $$\mathbf{N}=\mathbf N_{AB}\text{.}
  \label{eq:NAB}$$

We can now determine the integral form of the time delay function by differentiating $\Delta(\mathbf{z}_-(\lambda),x_B)$ with respect to $\lambda$. Using Eq. , we can always write $$\frac{{\mathrm{d}}}{{\mathrm{d}}\lambda}\Delta(\mathbf{z}_-(\lambda),x_B)=-R_{AB}N_{AB}^i\big[\partial_i\Delta\big]_{(\mathbf{z}_-(\lambda),x_B)}
  \label{eq:DTdel}$$ where $[\partial_i\Delta]_{(\mathbf{z}_-(\lambda),x_B)}$ denotes the partial derivative of $\Delta(\mathbf{x},x_B)$ with respect to $x^i$ taken at $\mathbf x=\mathbf z_-(\lambda)$. Then, after inserting Eqs. and into , we infer $$\frac{{\mathrm{d}}}{{\mathrm{d}}\lambda}\Delta(\mathbf{z}_-(\lambda),x_B)=\frac{R_{AB}}{2}\,\Omega_-(\widetilde{z}_-(\lambda),x_B)
  \label{eq:DTdeleik}$$ where the coordinates of the point-event $\widetilde{z}_-(\lambda)$ are obtained from Eq. , that is $\widetilde{\bm z}_-(\lambda)=\bm x_-(\mathbf z_-(\lambda))$. The components are given explicitly later on in Eq. .

By fixing the following boundary conditions

\[eq:bdR\] $$\begin{aligned}
  \Delta(\mathbf{z}_-(0),x_B)&=0\text{,}\label{eq:bdR0}\\
  \Delta(\mathbf{z}_-(1),x_B)&=\Delta(\mathbf{x}_A,x_B)\text{,}\label{eq:bdR1}\end{aligned}$$

which follow from the requirement that $\Delta(\mathbf{x}_B,x_B)=0$ when $\mathbf{z}_-(0)=\mathbf{x}_B$, we find $$\Delta(\mathbf{x}_A,x_B)=\frac{R_{AB}}{2}\int_0^1\Omega_-(\widetilde{z}_-(\lambda),x_B)\,{\mathrm{d}}\lambda\text{.}
  \label{eq:intdel}$$ Then, insertion of Eq. allows us to recover the theorem 2 of [@2008CQGra..25n5020T] which would be expressed here in terms of the contravariant components $K^{\mu\nu}$ instead of $k^{\mu\nu}$’s.

In principle, the all machinery developed in [@2008CQGra..25n5020T] for computing the delay functions could be applied using directly the components $K^{\mu\nu}$. However, such an approach possesses the inconvenience of hiding the role played by the different components $k^{\mu\nu}$ and $\kappa^{\mu\nu}$ during the determination of the total delay functions.

Indeed, according to Eqs. the curvature of the optical spacetime is described simultaneously with the help of the components $h_{\mu\nu}$ and $\gamma_{\mu\nu}$ which might act on different characteristic lengths (e.g. $\gamma_{\mu\nu}(x)=0$ for $x\notin\mathcal{D}$) and might possess completely different orders of magnitude *a priori*. Therefore, in order to disentangle the contribution of each perturbation into the determination of the total delay, we must perform a complete separation between the physical quantities in Eq. .

As it might be seen from Eqs. and such a separation can be achieved when the total time delay functions take the following forms

\[eq:del\] $$\begin{aligned}
  \Delta(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{g}}(\mathbf{x}_A,x_B)+\Delta_{\mathrm{r}}(\mathbf{x}_A,x_B)\nonumber\\
  &+\Delta_{\mathrm{gr}}(\mathbf{x}_A,x_B)\text{,}\label{eq:delR}\end{aligned}$$ and $$\begin{aligned}
  \Xi(x_A,\mathbf{x}_B)&=\Xi_{\mathrm{g}}(x_A,\mathbf{x}_B)+\Xi_{\mathrm{r}}(x_A,\mathbf{x}_B)\nonumber\\
  &+\Xi_{\mathrm{gr}}(x_A,\mathbf{x}_B)\text{.}
  \label{eq:delE}\end{aligned}$$

The subscripts $\mathrm{g}$, $\mathrm{r}$, and $\mathrm{gr}$ refer to the gravitational, the refractive, and the coupling contributions, respectively.

The gravitational and the refractive time delay functions are expected to be driven by gravitational and refractive perturbations, respectively. Instead, the coupling time delay functions is expected to be of the order of the product of both perturbations.

By Substituting for $\Delta(\mathbf{x}_A,x_B)$ from into and , and then by making use of the contravariant components of the optical and the physical spacetime metrics (see Eq. ), we deduce the following theorem.

In the optical spacetime, the function $\Delta$ introduced in Eq. , can be decomposed as shown in Eq.  where each term in the summation satisfy an integro-differential equation

\[eq:intdeltR\] $$\begin{aligned}
      \Delta_{\mathrm{g}}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int_0^1\bigg\{\big(k^{00}-2k^{0i}N^i_{AB}+k^{ij}N^i_{AB}N^j_{AB}\big)_{\widetilde{z}_-(\lambda)}\nonumber\\
      &+2\big(k^{0i}-k^{ij}N^j_{AB}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{g}}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}+\big(\eta^{ij}+k^{ij}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{g}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\bigg\}{\mathrm{d}}\lambda\text{,}\label{eq:intdeltRg}\\
      \Delta_{\mathrm{r}}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int_0^1\bigg\{\big(\kappa^{00}-2\kappa^{0i}N^i_{AB}+\kappa^{ij}N^i_{AB}N^j_{AB}\big)_{\widetilde{z}_-(\lambda)}\nonumber\\
      &+2\big(\kappa^{0i}-\kappa^{ij}N^j_{AB}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{r}}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}+\big(\eta^{ij}+\kappa^{ij}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{r}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{r}}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\bigg\}{\mathrm{d}}\lambda\text{,}\label{eq:intdeltRr}
    \end{aligned}$$ and $$\begin{aligned}
      \Delta_{\mathrm{gr}}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int_0^1\bigg\{\big(\eta^{ij}+k^{ij}+\kappa^{ij}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{gr}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+2\big(k^{0i}+\kappa^{0i}-(k^{ij}+\kappa^{ij})N^j_{AB}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{gr}}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\big(k^{ij}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{r}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{r}}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}+2\big(k^{0i}-k^{ij}N^j_{AB}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{r}}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\big(\kappa^{ij}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{g}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}+2\big(\kappa^{0i}-\kappa^{ij}N^j_{AB}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{g}}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+2\big(\eta^{ij}+k^{ij}+\kappa^{ij}\big)_{\widetilde{z}_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{g}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{r}}}{\partial x^j}+\frac{\partial\Delta_{\mathrm{g}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}}{\partial x^j}+\frac{\partial\Delta_{\mathrm{r}}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\bigg\}{\mathrm{d}}\lambda\text{.}\label{eq:intdeltRgr}
    \end{aligned}$$

The coordinates of the point-event $\widetilde{z}_-(\lambda)$ are given by $$\widetilde{\bm z}_-(\lambda)=\big(x_B^0-\lambda R_{AB}-\Delta(\mathbf z_-(\lambda),x_B),\mathbf z_-(\lambda)\big)\text{,}
      \label{eq:tz-m}$$ $\mathbf z_-(\lambda)$ being defined in Eq. . \[th:theo2\]

Following the exact same reasoning, we state a similar theorem for the emission time delay function. However, for the emission case, the straight line segment connecting the emitter $x_A$ to the receiver $x_B$ is defined by [^4] $$\begin{aligned}
  \mathbf z_+(\mu)&=\mathbf x_A+\mu R_{AB}\mathbf{N}_{AB}\text{,} & 0\leqslant\mu\leqslant 1\text{.}\label{eq:z+i}\end{aligned}$$ Then, from the requirement that $\Xi(x_A,\mathbf{x}_A)=0$ when $\mathbf{z}_+(0)=\mathbf{x}_A$, we can set the following boundary conditions

\[eq:bdE\] $$\begin{aligned}
  \Xi(x_A,\mathbf{z}_+(0))&=0\text{,}\label{eq:bdE0}\\
  \Xi(x_A,\mathbf{z}_+(1))&=\Xi(x_A,\mathbf{x}_B)\text{.}\label{eq:bdE1}\end{aligned}$$

Hence, the theorem for the emission time delay function $\Xi(x_A,\mathbf x_B)/c$ reads as follows.

In the optical spacetime, the function $\Xi$ introduced in Eq. , can be decomposed as shown in Eq.  where each term in the summation satisfy an integro-differential equation

\[eq:intdeltE\] $$\begin{aligned}
      \Xi_{\mathrm{g}}(x_A,\mathbf x_B)&=\frac{R_{AB}}{2}\int_0^1\bigg\{\big(k^{00}-2k^{0i}N^i_{AB}+k^{ij}N^i_{AB}N^j_{AB}\big)_{\widetilde{z}_+(\mu)}\nonumber\\
      &-2\big(k^{0i}-k^{ij}N^j_{AB}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{g}}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}+\big(\eta^{ij}+k^{ij}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{g}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\bigg\}{\mathrm{d}}\mu\text{,}\label{eq:intdeltEg}\\
      \Xi_{\mathrm{r}}(x_A,\mathbf x_B)&=\frac{R_{AB}}{2}\int_0^1\bigg\{\big(\kappa^{00}-2\kappa^{0i}N^i_{AB}+\kappa^{ij}N^i_{AB}N^j_{AB}\big)_{\widetilde{z}_+(\mu)}\nonumber\\
      &-2\big(\kappa^{0i}-\kappa^{ij}N^j_{AB}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{r}}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}+\big(\eta^{ij}+\kappa^{ij}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{r}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{r}}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\bigg\}{\mathrm{d}}\mu\text{,}\label{eq:intdeltEr}
    \end{aligned}$$ and $$\begin{aligned}
      \Xi_{\mathrm{gr}}(x_A,\mathbf x_B)&=\frac{R_{AB}}{2}\int_0^1\bigg\{\big(\eta^{ij}+k^{ij}+\kappa^{ij}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{gr}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &-2\big(k^{0i}+\kappa^{0i}-(k^{ij}+\kappa^{ij})N^j_{AB}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{gr}}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\lambda))}\nonumber\\
      &+\big(k^{ij}\big)_{\widetilde{z}_+(\lambda)}\left[\frac{\partial\Xi_{\mathrm{r}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{r}}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}-2\big(k^{0i}-k^{ij}N^j_{AB}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{r}}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\big(\kappa^{ij}\big)_{\widetilde{z}(\lambda)}\left[\frac{\partial\Xi_{\mathrm{g}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}-2\big(\kappa^{0i}-\kappa^{ij}N^j_{AB}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{g}}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+2\big(\eta^{ij}+k^{ij}+\kappa^{ij}\big)_{\widetilde{z}_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{g}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{r}}}{\partial x^j}+\frac{\partial\Xi_{\mathrm{g}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}}{\partial x^j}+\frac{\partial\Xi_{\mathrm{r}}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\bigg\}{\mathrm{d}}\mu\text{.}\label{eq:intdeltEgr}
    \end{aligned}$$

The coordinates of the point-event $\widetilde{z}_+(\mu)$ are given by $$\widetilde{\bm z}_+(\mu)=\big(x_A^0+\mu R_{AB}+\Xi(x_A,\mathbf z_+(\mu)),\mathbf z_+(\mu)\big)\text{,}
      \label{eq:tz+m}$$ $\mathbf z_+(\mu)$ being defined in Eq. . \[th:theo3\]

Theorems \[th:theo2\] and \[th:theo3\] generalize theorems 2 and 3 of [@2008CQGra..25n5020T] for the optical spacetime. Indeed, in the limit where refractivity vanishes, that is to say $\kappa^{\mu\nu}\rightarrow 0$, theorems 2 and 3 of [@2008CQGra..25n5020T] are recovered.

From Eqs. , we see that the choice does achieve the separation between the different physical quantities entering the computation of the total time delay. As a matter of fact, the right-hand sides of Eqs.  and contains purely gravitational and purely refractive quantities, respectively. The right-hand side of Eq. regroups all terms being a mixture of both.

However, as it may be observed from the presence of the total delay in Eq. , the expressions of the different contributions are not fully independents but remain linked via the path of integration. In the next section, we shall further discuss this point and shall present a recursive resolution method for determining the time delay functions at any order.

General expansions of the time delay functions {#sec:exp}
==============================================

Because the line integrals in Eqs.  are taken along the path $\widetilde{z}_-(\lambda)$ for $0\leqslant\lambda\leqslant1$, the time delay functions $\Delta_{\mathrm{g}}/c$, $\Delta_{\mathrm{r}}/c$, and $\Delta_{\mathrm{gr}}/c$ cannot be solved independently from each other. Indeed, the total delay appearing in Eq.  depends on the three functions as it can be seen from the decomposition . Therefore, a systematic and recursive resolution of $\Delta/c$ can only be achieved, once the relative contributions of $\Delta_{\mathrm{g}}/c$, $\Delta_{\mathrm{r}}/c$, and $\Delta_{\mathrm{gr}}/c$ to the total time delay are known.

In Sec. \[subsec:QM\], we first show how to determine the relative importance between the different contributions. Then, within the approximation of a quasi-Minkowskian path, we show that the interdependence between each functions $\Delta_{\mathrm{g}}/c$, $\Delta_{\mathrm{r}}/c$, and $\Delta_{\mathrm{gr}}/c$ can always be rejected to the following order during the resolution of $\Delta/c$. This fact allows one to sort out the occurrence of the different contributions within the determination of the total delay function (cf. theorems \[th:theo4\] and \[th:theo5\]). In Sec. \[sec:detref\], we assume that the refractive components of the optical metric admit a series expansion in term a parameter $N_0$. Then, we show that the refractive delay functions can be determined to any order through a recursive resolution method presented in theorems \[th:theo6\] and \[th:theo7\]. In Sec. \[sec:detgra\], we assume that the gravitational components of the spacetime metric admit a post-Minkowskian expansion (series expansion in ascending power of $G$). Then, the recursive method allowing one to determine the gravitational delay expressions to any order is presented in theorems \[th:theo8\] and \[th:theo9\]. Finally in Sec. \[sec:detcou\], we determine the coupling delay expressions up to any order within theorems \[th:theo10\] and \[th:theo11\].

Quasi-Minkowskian path regime {#subsec:QM}
-----------------------------

As shown in theorems \[th:theo2\] and \[th:theo3\], the relative magnitude between each contribution to the total time delay rely on the line integrals of the gravitational and the refractive perturbations. Generally speaking, if gravity acts all along the light path $\Gamma$ joining $x_A$ to $x_B$, the refractive domain $\mathcal{D}$ is localized in spacetime and it follows that the action of refractivity remains bounded to a certain portion of $\Gamma$. Therefore, in order to determine the relative contributions of each time delay functions, not only the relative magnitude between the gravitational and refractive perturbations must be known, but also the typical length scales over which each perturbation acts.

Let $\ell$ be the length of $\Gamma$ passing through $\mathcal{D}$. For a Minkowskian path, we always have $\ell\leqslant R_{AB}$ whatever is the size of $\mathcal{D}$.

From Eq.  which has been formulated under the assumption that the light path is quasi-Minkowskian, we deduce $\Delta/R_{AB}\ll 1$. This implies that $\Delta_{\mathrm{g}}/R_{AB}\ll 1$, $\Delta_{\mathrm{r}}/R_{AB}\ll 1$, and $\Delta_{\mathrm{gr}}/R_{AB}\ll 1$. Considering that $\Delta_{\mathrm{gr}}$ represents the coupling contributions, its magnitude is expected to be of the order $$\frac{\Delta_{\mathrm{gr}}}{R_{AB}}\sim\left(\frac{\Delta_{\mathrm{g}}}{R_{AB}}\right)\left(\frac{\Delta_{\mathrm{r}}}{R_{AB}}\right)\text{.}$$ Therefore, we can first focus on the relative importance between the gravitational and refractive contributions.

To do so, let us introduce the parameter $s$ defined by $$s=\left\lfloor\frac{\log_{10}(\Delta_{\mathrm{g}}/R_{AB})}{\log_{10}(\Delta_{\mathrm{r}}/R_{AB})}\right\rceil\text{,}
  \label{eq:s}$$ with $\lfloor i\rceil$ denoting the operation of rounding to the nearest integer of $i$. Hereafter, we intend to show that the expansion pattern of the delay functions can totally be determined once $s$ is known. Indeed, it allows one to sort out the occurrences of the gravitational and refractive terms in the determination of the total delay functions.

Because we are only focusing on the main integer value of $s$ in Eq. , it is sufficient to get the first-order expressions of $\Delta_{\mathrm{g}}$ and $\Delta_{\mathrm{r}}$. Therefore, in Eqs. , line integrals can be changed into line integrals along the Minkowskian path between $x_A$ and $x_B$ by performing a Taylor series expansion of $\kappa^{\mu\nu}(\widetilde{z}_-(\lambda))$ and $k^{\mu\nu}(\widetilde{z}_-(\lambda))$ about the point-event $z_-(\lambda)$ whose coordinates are given by $$\bm z_-(\lambda)=\big(x_B^0-\lambda R_{AB},\mathbf{z}_-(\lambda)\big)\text{.}
  \label{eq:z-m}$$ Thus, optical metric components become infinite series in ascending power of the total time delay

\[eq:Tay\] $$\begin{aligned}
  &k^{\mu\nu}\left(\widetilde{z}_-(\lambda),\frac{\Delta(\mathbf{z}_-(\lambda),x_B)}{R_{AB}}\right)=k^{\mu\nu}(z_-(\lambda))\nonumber\\
  &+\sum_{l=1}^{\infty}\frac{(-R_{AB})^l}{l!}\left(\frac{\Delta(\mathbf{z}_-(\lambda),x_B)}{R_{AB}}\right)^l\big[\partial_0^lk^{\mu\nu}\big]_{z_-(\lambda)}\text{,}
  \label{eq:kTay}\end{aligned}$$ and $$\begin{aligned}
  &\kappa^{\mu\nu}\left(\widetilde{z}_-(\lambda),\frac{\Delta(\mathbf{z}_-(\lambda),x_B)}{R_{AB}}\right)=\kappa^{\mu\nu}(z_-(\lambda))\nonumber\\
  &+\sum_{l=1}^{\infty}\frac{(-R_{AB})^l}{l!}\left(\frac{\Delta(\mathbf{z}_-(\lambda),x_B)}{R_{AB}}\right)^{l}\big[\partial_0^l\kappa^{\mu\nu}\big]_{z_-(\lambda)}\text{.}
  \label{eq:kapTay}\end{aligned}$$

After inserting these expressions into Eqs. and , we infer that the zeroth-order terms in Eqs. , namely $k^{\mu\nu}(z_-(\lambda))$ and $\kappa^{\mu\nu}(z_-(\lambda))$, correspond to the first-order determination of the gravitational and the refractive delays

\[eq:dgdrrelmag\] $$\frac{\Delta_{\mathrm{g}}^{(1)}}{R_{AB}}=\frac{1}{2}\int_0^1(k^{00}-2k^{0i}N^i_{AB}+k^{ij}N^i_{AB}N^j_{AB})_{z_-(\lambda)}{\mathrm{d}}\lambda\text{,}$$ and $$\frac{\Delta_{\mathrm{r}}^{(1)}}{R_{AB}}=\frac{1}{2}\int_0^1(\kappa^{00}-2\kappa^{0i}N^i_{AB}+\kappa^{ij}N^i_{AB}N^j_{AB})_{z_-(\lambda)}{\mathrm{d}}\lambda\text{.}$$

(We will see with Eqs. and , that in the context of a quasi-Minkowskian path, these equations can be further simplified. But for now let us pursue the discussion with Eqs. ). These equations can be inserted into in order to determine the value of $s$. We shall discuss now how the expansion pattern of the delay functions can be worked out form $s$.

Henceforth, we consider the case $s\in\mathbb{N}_{>1}$ (the result will be still valid for $s\in\mathbb{N}_{>0}$). In other words, we suppose that the refractive perturbation is dominant with respect to the gravitational one [^5].

In order to simplify the next discussion, and without loss of generality, we focus on orders of magnitude only. In addition, we consider that the light path occurs in a sufficiently small region of spacetime where the metric components do not vary significantly. Thus, we deduce from Eqs. that $$\frac{\Delta_{\mathrm{g}}^{(1)}}{R_{AB}}\sim |k^{\mu\nu}|_{\mathrm{max}}\text{,} \qquad \frac{\Delta_{\mathrm{r}}^{(1)}}{R_{AB}}\sim\frac{\ell}{R_{AB}}|\kappa^{\mu\nu}|_{\mathrm{max}}\text{.}
  \label{eq:del0ord}$$ In order to keep track of the relative magnitude between the gravitational and the refractive terms, we introduce a dimensionless parameter denoted $\varepsilon$ and being of the order of the dominant term, that is $$\varepsilon=\frac{\Delta_{\mathrm{r}}^{(1)}}{R_{AB}}\text{.}
  \label{eq:varepDr}$$ Thus, from Eqs. and , we immediately infer $$\frac{\Delta_{\mathrm{g}}^{(1)}}{R_{AB}}\sim\mathcal{O}(\varepsilon^s)\text{.}
  \label{eq:varepDg}$$ Therefore, the first-order expression of the total delay is driven by the refractive term only $$\Delta^{(1)}(\mathbf{x}_A,x_B)=\Delta_{\mathrm{r}}^{(1)}(\mathbf{x}_A,x_B)
  \label{eq:del1delr}$$ which means that $$\frac{\Delta^{(1)}}{R_{AB}}=\varepsilon\text{,}
  \label{eq:deltot1}$$ when $s>1$ in Eq. .

Let us take a look at the relation between metric components. Eqs. , , and , allow us to deduce $$|k^{\mu\nu}|_{\mathrm{max}}\sim\mathcal{O}(\varepsilon^s)\text{,} \qquad \frac{\ell}{R_{AB}}|\kappa^{\mu\nu}|_{\mathrm{max}}\sim\mathcal{O}(\varepsilon)\text{.}
  \label{eq:kkap0}$$ At the same time, it might be seen from Eqs. , , and that the second-order refractive delay is driven by terms such like $$\frac{\ell}{R_{AB}}|\kappa^{\mu\nu}|_{\mathrm{max}}\frac{\Delta^{(1)}}{R_{AB}}\text{,} \qquad \frac{\Delta^{(1)}}{R_{AB}}\frac{\Delta^{(1)}}{R_{AB}}\text{,}$$ which, according to Eqs. and , are of the order of $\varepsilon^2$. Therefore, we conclude that the series expansion of the total delay goes on like $$\Delta^{(l)}(\mathbf{x}_A,x_B)=\Delta_{\mathrm{r}}^{(l)}(\mathbf{x}_A,x_B)\text{.}$$ for $1\leqslant l<s$.

The first occurrence of the gravitational contribution to the total delay arises for $l=s$ as anticipated in Eq. . Therefore, the $s$th-order expression of the total delay is given by $$\Delta^{(s)}(\mathbf{x}_A,x_B)=\Delta_{\mathrm{r}}^{(s)}(\mathbf{x}_A,x_B)+\Delta_{\mathrm{g}}^{(1)}(\mathbf{x}_A,x_B)\text{.}$$ Then, by looking at the first-order term in Eq. , one might see that the second-order expression of the gravitational delay is proportional to $$\Delta_{\mathrm g}^{(2)}\sim|k^{\mu\nu}|_{\mathrm{max}}\frac{\Delta^{(1)}}{R_{AB}}
  \label{eq:delg2}$$ which, according to Eqs. and , is of the order of $\varepsilon^{s+1}$. Additionally, after inserting Eqs. into , we infer that the first-order expression of the coupling delay is driven by terms such like $$\frac{\ell}{R_{AB}}|\kappa^{\mu\nu}|_{\mathrm{max}}\frac{\Delta_{\mathrm g}^{(1)}}{R_{AB}}\text{,} \qquad |k^{\mu\nu}|_{\mathrm{max}}\frac{\Delta_{\mathrm{r}}^{(1)}}{R_{AB}}\text{,} \qquad \frac{\Delta_{\mathrm{r}}^{(1)}}{R_{AB}}\frac{\Delta_{\mathrm g}^{(1)}}{R_{AB}}\text{,}$$ which are of the order of $\varepsilon^{s+1}$ too. Therefore, the $(s+1)$th-order expression of the total delay is given by $$\begin{aligned}
  \Delta^{(s+1)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{r}}^{(s+1)}(\mathbf{x}_A,x_B)+\Delta_{\mathrm{g}}^{(2)}(\mathbf{x}_A,x_B)\nonumber\\
   &+\Delta_{\mathrm{gr}}^{(1)}(\mathbf{x}_A,x_B)\text{.}\end{aligned}$$

A quick look at the second-order expression of the coupling delay shows that it is driven by terms proportional to $\varepsilon^{s+2}$. Consequently, one deduces that $$\begin{aligned}
  \Delta^{(l)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{r}}^{(l)}(\mathbf{x}_A,x_B)+\Delta_{\mathrm{g}}^{(l-s+1)}(\mathbf{x}_A,x_B)\nonumber\\
   &+\Delta_{\mathrm{gr}}^{(l-s)}(\mathbf{x}_A,x_B)\text{,}\end{aligned}$$ for $l\geqslant s+1$.

Let us summarize the discussion. Within the quasi-Minkowskian regime, the total delay satisfies $\Delta/R_{AB}\ll 1$, so the line integrals in Eqs. are simplified into line integrals along the Minkowskian path by performing a Taylor series expansion about the point-event $z_-(\lambda)$. Then, by assuming that refractivity is the dominant effect all along the light path $\Gamma$, it results that, in general, the total time delay admits an expansion as follows $$\Delta(\mathbf{x}_A,x_B)=\sum_{l=1}^{\infty}\Delta^{(l)}(\mathbf{x}_A,x_B)
  \label{eq:delRPM}$$ where the terms $\Delta^{(l)}$ are proportional to $\varepsilon^lR_{AB}$. In that respect, the different contributions to the total delay, namely the refractive, the gravitational, and the coupling delays, all admit series expansion as follows

\[eq:delRrgrPM\] $$\begin{aligned}
  \Delta_{\mathrm{r}}(\mathbf{x}_A,x_B)&=\sum_{l=1}^{\infty}\Delta_{\mathrm{r}}^{(l)}(\mathbf{x}_A,x_B)\text{,}\label{eq:delRrPM}\\
  \Delta_{\mathrm{g}}(\mathbf{x}_A,x_B)&=\sum_{l=1}^{\infty}\Delta_{\mathrm{g}}^{(l)}(\mathbf{x}_A,x_B)\text{,}\label{eq:delRgPM}\end{aligned}$$ and $$\Delta_{\mathrm{gr}}(\mathbf{x}_A,x_B)=\sum_{l=1}^{\infty}\Delta_{\mathrm{gr}}^{(l)}(\mathbf{x}_A,x_B)
  \label{eq:delRgrPM}$$

where the terms $\Delta_{\mathrm{r}}^{(l)}$, $\Delta_{\mathrm{g}}^{(l)}$, and $\Delta_{\mathrm{gr}}^{(l)}$ are of the order of $$\frac{\Delta_{\mathrm{r}}^{(l)}}{R_{AB}}\sim\mathcal{O}(\varepsilon^{l})\text{,} \quad \frac{\Delta_{\mathrm{g}}^{(l)}}{R_{AB}}\sim\mathcal{O}(\varepsilon^{l+s-1})\text{,} \quad \frac{\Delta_{\mathrm{gr}}^{(l)}}{R_{AB}}\sim\mathcal{O}(\varepsilon^{l+s})\text{.}$$ By making use of the Heaviside step function $$\Theta(i)=\left\{
  \begin{array}{l l}
    1 & \text{for }i\geqslant 0\text{,}\\
    0 & \text{otherwise,}
  \end{array}
  \right.
  \label{eq:stp}$$ we can write the terms $\Delta^{(l)}$ in Eq. as $$\begin{aligned}
  \Delta^{(l)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{r}}^{(l)}(\mathbf{x}_A,x_B)+\Theta(l\!-\!s)\Delta_{\mathrm{g}}^{(l-s+1)}(\mathbf{x}_A,x_B)\nonumber\\
  &+\Theta(l\!-\!s\!-\!1)\Delta_{\mathrm{gr}}^{(l-s)}(\mathbf{x}_A,x_B)\text{.}\label{eq:delRPMl}\end{aligned}$$ Let us recall that $\varepsilon$ is of the order of $\ell/R_{AB}|\kappa^{\mu\nu}|_{\mathrm{max}}$ only for a light path occurring in a sufficiently small region of spacetime where the metric components do not vary significantly. In general, it is given by Eq. .

In the next two sections, according to the fact that the light ray follows a quasi-Minkowskian path, we will assume that the components $\kappa^{\mu\nu}$ and $k^{\mu\nu}$ admit series expansion in ascending power of parameters $N_0$ and $G$, respectively (see Eqs.  and ). If this fact does not change the pattern of the series expansions and , we should nevertheless, for completeness, specify that the quasi-Minkowskian path is parametrized by the expansion coefficients $N_0$ and $G$. Therefore, by making use of Eq. , we can state a theorem as follows.

Within the approximation that the light path is quasi-Minkowskian and is parametrized by $N_0$ and $G$, the function $\Delta$ admits a series expansion as follows $$\Delta(\mathbf{x}_A,x_B,N_0,G)=\sum_{l=1}^{\infty}\Delta^{(l)}(\mathbf{x}_A,x_B)\text{,}
      \label{eq:delRPMN0G}$$ with $$\Delta^{(l)}(\mathbf{x}_A,x_B)=\Delta_{\mathrm{r}}^{(l)}(\mathbf{x}_A,x_B)+\Theta(l-s)\Delta_{\mathrm{g}}^{(l-s+1)}(\mathbf{x}_A,x_B)+\Theta(l-s-1)\Delta_{\mathrm{gr}}^{(l-s)}(\mathbf{x}_A,x_B)\text{.}$$ The parameter $s\in\mathbb{N}_{>0}$ is determined from Eq. by making use of the first-order expressions and . \[th:theo4\]

A similar reasoning works for the emission time delay function as well. Indeed, the line integrals in Eqs. can be expanded in a Taylor series about the point-event $z_+(\mu)$ defined by $$\bm z_+(\mu)=\big(x_B^0+\mu R_{AB},\mathbf{z}_+(\mu)\big)\text{.}
  \label{eq:z+m}$$ Then, $\kappa^{\mu\nu}(\widetilde{z}_+(\mu))$ and $k^{\mu\nu}(\widetilde{z}_+(\mu))$ become infinite series in ascending power of $\Xi$ similarly to what have been done in Eqs. .

Therefore, we end up with similar expansion for $\Xi$ than for $\Delta$ and we state the following theorem.

Within the approximation that the light path is quasi-Minkowskian and is parametrized by $N_0$ and $G$, the function $\Xi$ admits a series expansion as follows $$\Xi(x_A,\mathbf{x}_B,N_0,G)=\sum_{l=1}^{\infty}\Xi^{(l)}(x_A,\mathbf{x}_B)\text{,}
      \label{eq:delEPMN0G}$$ with $$\Xi^{(l)}(x_A,\mathbf{x}_B)=\Xi_{\mathrm{r}}^{(l)}(x_A,\mathbf{x}_B)+\Theta(l-s)\Xi_{\mathrm{g}}^{(l-s+1)}(x_A,\mathbf{x}_B)+\Theta(l-s-1)\Xi_{\mathrm{gr}}^{(l-s)}(x_A,\mathbf{x}_B)\text{.}$$ The parameter $s\in\mathbb{N}_{>0}$ is determined from Eq. by making use of the first-order expressions and . \[th:theo5\]

Equipped with theorems \[th:theo2\] to \[th:theo5\], we can now recursively determine the integral form of each component in the time delay function expressions .

The refractive time delay functions {#sec:detref}
-----------------------------------

We saw in Sec. \[sec:intTF\] that a quasi-Minkowskian path implies the small refractivity approximation, that is to say $N(x)\ll1$. Let $N_0$ be the refractivity at a well-chosen point-event $x_0\in\mathcal{D}$ located on $\Gamma$. $N_0$ is a constant defined by $N_0=N(x_0)$. We can always write $$N(x)=N_0\left(\frac{N(x)}{N_0}\right)\text{.}$$ Hence, considering a quasi-Minkowskian light path, it follows that $N_0\ll1$. Therefore, we can always expand the components $\kappa^{\mu\nu}$ in ascending power of $N_0$ such like $$\kappa^{\mu\nu}(x,N_0)=\sum_{l=1}^{\infty}\kappa^{\mu\nu}_{(l)}(x)
  \label{eq:kapPM}$$ where $\kappa^{\mu\nu}_{(l)}\propto(N_0)^l$.

Considering that the wave-vector is by definition a covector (see Eq. ), the optical metric is intrinsically defined for its contravariant components as seen from Eq. . Therefore, the covariant components of the optical metric are not needed *a priori* to solve the time and frequency transfers. However, we provide their expressions in Sec. \[sec:genexp\] for completeness.

As discussed previously, line integrations in Eq.  are taken along the real light path $\widetilde{z}_-(\lambda)$ for $0\leqslant\lambda\leqslant 1$. Within the quasi-Minkowskian path regime, we saw in Eq. that the metric components $\kappa^{\mu\nu}(\widetilde{z}_-(\lambda))$ can be expanded in ascending power of the total delay. In that respect, the right-hand side of Eq.  involves terms such as $\kappa^{\mu\nu}(z_-(\lambda))$. By making use of Eq. , we immediately find $$\kappa^{\mu\nu}(z_-(\lambda),N_0)=\sum_{l=1}^{\infty}\kappa^{\mu\nu}_{(l)}(z_-(\lambda))\text{.}
  \label{eq:kapPMz-}$$

Therefore, the general expansion of $\kappa^{\mu\nu}(\widetilde{z}_-(\lambda))$ is obtained after substituting for $\Delta$ and $\kappa^{\mu\nu}$ from Eqs.  and into , respectively. After some algebra, we find a relation as follows $$\kappa^{\mu\nu}(\widetilde{z}_-(\lambda),N_0,G)=\sum_{l=1}^{\infty}\widehat{\kappa}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)
  \label{eq:kapRPM}$$ where the quantities $\widehat{\kappa}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)$ are given by

\[eq:kapRPMdec\] $$\widehat{\kappa}^{\mu\nu}_{-(1)}(z_-(\lambda),x_B)=\kappa^{\mu\nu}_{(1)}(z_-(\lambda))\text{,}$$ and $$\begin{aligned}
  &\widehat{\kappa}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)=\kappa^{\mu\nu}_{(l)}(z_-(\lambda))\nonumber\\
  &+\sum_{m=1}^{l-1}\sum_{n=1}^{m}\Phi_-^{(m,n)}(\mathbf z_-(\lambda),x_B)\left[\frac{\partial^n\kappa^{\mu\nu}_{(l-m)}}{(\partial x^0)^n}\right]_{z_-(\lambda)}\end{aligned}$$

for $l\geqslant 2$.

The function $\Phi_-^{(l,m)}(\mathbf x,x_B)$ is called a reception function and is defined [@2008CQGra..25n5020T] such that $$\begin{aligned}
  \Phi_-^{(l,m)}&(\mathbf x,x_B)=\frac{(-1)^m}{m!}\nonumber\\
  &\times\sum_{n_1+\cdots+n_m=l-m}\Bigg[\prod_{d=1}^m\Delta^{(n_d+1)}(\mathbf x,x_B)\Bigg]\text{,}\label{eq:recfun}\end{aligned}$$ with $l\geqslant 1$ and $1\leqslant m\leqslant l$, and with $n_1,\ldots,n_m\in\mathbb{N}_{\geqslant 0}$. The summation in is taken over all sequences of $n_1$ through $n_m$ such that the sum of all $n_m$ is $l-m$.

Finally, by substituting for $\kappa^{\mu\nu}(\widetilde{z}_-(\lambda))$ from Eq. into , we infer the theorem which follows concerning the refractive time delay function at reception.

In the optical spacetime, within the quasi-Minkowskian path regime, $\Delta$ admits the series expansion introduced in theorem \[th:theo4\], so the function $\Delta_{\mathrm{r}}$ is given by $$\Delta_{\mathrm{r}}(\mathbf{x}_A,x_B,N_0,G)=\sum_{l=1}^{\infty}\Delta_{\mathrm{r}}^{(l)}(\mathbf{x}_A,x_B)
      \label{eq:delRrPMN0G}$$ where

\[eq:intdelRrPM\] $$\begin{aligned}
      \Delta_{\mathrm{r}}^{(1)}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\big(\kappa^{00}_{(1)}-2\kappa^{0i}_{(1)}N^i_{AB}+\kappa^{ij}_{(1)}N^i_{AB}N^j_{AB}\big)_{z_-(\lambda)}{\mathrm{d}}\lambda\text{,}\label{eq:intdelRrPM1}\\
      \Delta_{\mathrm{r}}^{(2)}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{\kappa}^{00}_{-(2)}-2\widehat{\kappa}^{0i}_{-(2)}N^i_{AB}+\widehat{\kappa}^{ij}_{-(2)}N^i_{AB}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\nonumber\\
      &+2\big(\kappa^{0i}_{(1)}-\kappa^{ij}_{(1)}N_{AB}^j\big)_{z_-(\lambda)}\left[\frac{\partial\Delta_{\mathrm{r}}^{(1)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}+\eta^{ij}\left[\frac{\partial\Delta_{\mathrm{r}}^{(1)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{r}}^{(1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda\text{,}
    \end{aligned}$$ and $$\begin{aligned}
      \Delta_{\mathrm{r}}^{(l)}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{\kappa}^{00}_{-(l)}-2\widehat{\kappa}^{0i}_{-(l)}N^i_{AB}+\widehat{\kappa}^{ij}_{-(l)}N^i_{AB}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\nonumber\\
      &+2\sum_{m=1}^{l-1}\big(\widehat\kappa^{0i}_{-(m)}-\widehat\kappa^{ij}_{-(m)}N_{AB}^j\big)_{(z_-(\lambda),x_B)}\left[\frac{\partial\Delta_{\mathrm{r}}^{(l-m)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}+\eta^{ij}\sum_{m=1}^{l-1}\left[\frac{\partial\Delta_{\mathrm{r}}^{(m)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{r}}^{(l-m)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=1}^{l-2}\big(\widehat{\kappa}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-m-1}\left[\frac{\partial\Delta_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{r}}^{(l-m-n)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$

for $l \geqslant 3$. The quantities $\widehat{\kappa}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)$ are defined in Eqs. . \[th:theo6\]

Applying the exact same reasoning, the analogous theorem for the refractive time delay function at emission can be stated as well. Line integrations in Eq.  are taken along the light path $\widetilde{z}_+(\mu)$ for $0\leqslant\mu\leqslant 1$. After Taylor expanding the light path about the point-event $z_+(\mu)$, the right-hand side of Eq.  involves terms such as $\kappa^{\mu\nu}(z_+(\mu))$. By making use of Eq. , we find $$\kappa^{\mu\nu}(z_+(\mu),N_0)=\sum_{l=1}^{\infty}\kappa^{\mu\nu}_{(l)}(z_+(\mu))\text{.}
  \label{eq:kapPMz+}$$

The general expansion of $\kappa^{\mu\nu}(\widetilde{z}_+(\mu))$ is obtained after substituting for $\Xi$ from Eq.  and for $\kappa^{\mu\nu}$ from Eq. , into . After some algebra, we find $$\kappa^{\mu\nu}(\widetilde{z}_+(\mu),N_0,G)=\sum_{l=1}^{\infty}\widehat{\kappa}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))
  \label{eq:kapEPM}$$ where the quantities $\widehat{\kappa}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))$ are given by

\[eq:kapEPMdec\] $$\widehat{\kappa}^{\mu\nu}_{+(1)}(x_A,z_+(\mu))=\kappa^{\mu\nu}_{(1)}(z_+(\mu))\text{,}$$ and $$\begin{aligned}
  &\widehat{\kappa}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))=\kappa^{\mu\nu}_{(l)}(z_+(\mu))\nonumber\\
  &+\sum_{m=1}^{l-1}\sum_{n=1}^{m}\Phi_+^{(m,n)}(x_A,\mathbf z_+(\mu))\left[\frac{\partial^n\kappa^{\mu\nu}_{(l-m)}}{(\partial x^0)^n}\right]_{z_+(\mu)}\end{aligned}$$

for $l\geqslant 2$.

The function $\Phi_+^{(l,m)}(x_A,\mathbf x)$ is called an emission function and is defined such as $$\begin{aligned}
  \Phi_+^{(l,m)}&(x_A,\mathbf x)=\frac{1}{m!}\nonumber\\
  &\times\sum_{n_1+\cdots+n_m=l-m}\Bigg[\prod_{d=1}^m\Xi^{(n_d+1)}(x_A,\mathbf x)\Bigg]\text{,}\label{eq:emifun}\end{aligned}$$ with $l\geqslant 1$ and $1\leqslant m\leqslant l$, and with $n_1,\ldots,n_m\in\mathbb{N}_{\geqslant 0}$. The summation in is taken over all sequences of $n_1$ through $n_m$ such that the sum of all $n_m$ is $l-m$.

Finally, the theorem for the refractive time delay function at emission is obtained after substituting for $\kappa^{\mu\nu}(\widetilde{z}_+(\mu))$ from Eq. into .

In the optical spacetime, within the quasi-Minkowskian path regime, $\Xi$ admits the series expansion introduced in theorem \[th:theo5\], so the function $\Xi_{\mathrm{r}}$ is given by $$\Xi_{\mathrm{r}}(x_A,\mathbf{x}_B,N_0,G)=\sum_{l=1}^{\infty}\Xi_{\mathrm{r}}^{(l)}(x_A,\mathbf{x}_B)
      \label{eq:delErPM}$$ where

\[eq:intdelErPM\] $$\begin{aligned}
      \Xi_{\mathrm{r}}^{(1)}(x_A,\mathbf{x}_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\big(\kappa^{00}_{(1)}-2\kappa^{0i}_{(1)}N^i_{AB}+\kappa^{ij}_{(1)}N^i_{AB}N^j_{AB}\big)_{z_+(\mu)}{\mathrm{d}}\mu\text{,}\\
      \Xi_{\mathrm{r}}^{(2)}(x_A,\mathbf{x}_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{\kappa}^{00}_{+(2)}-2\widehat{\kappa}^{0i}_{+(2)}N^i_{AB}+\widehat{\kappa}^{ij}_{+(2)}N^i_{AB}N^j_{AB}\big)_{(x_A,z_+(\mu))}\nonumber\\
      &-2\big(\kappa^{0i}_{(1)}-\kappa^{ij}_{(1)}N_{AB}^j\big)_{z_+(\mu)}\left[\frac{\partial\Xi_{\mathrm{r}}^{(1)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}+\eta^{ij}\left[\frac{\partial\Xi_{\mathrm{r}}^{(1)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{r}}^{(1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu\text{,}
    \end{aligned}$$ and $$\begin{aligned}
      \Xi_{\mathrm{r}}^{(l)}(x_A,\mathbf{x}_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{\kappa}^{00}_{+(l)}-2\widehat{\kappa}^{0i}_{+(l)}N^i_{AB}+\widehat{\kappa}^{ij}_{+(l)}N^i_{AB}N^j_{AB}\big)_{(x_A,z_+(\mu))}\nonumber\\
      &-2\sum_{m=1}^{l-1}\big(\widehat\kappa^{0i}_{+(m)}-\widehat\kappa^{ij}_{+(m)}N_{AB}^j\big)_{(x_A,z_+(\mu))}\left[\frac{\partial\Xi_{\mathrm{r}}^{(l-m)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}+\eta^{ij}\sum_{m=1}^{l-1}\left[\frac{\partial\Xi_{\mathrm{r}}^{(m)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{r}}^{(l-m)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\sum_{m=1}^{l-2}\big(\widehat{\kappa}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-m-1}\left[\frac{\partial\Xi_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{r}}^{(l-m-n)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$

for $l \geqslant 3$. The quantities $\widehat{\kappa}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))$ are defined in Eqs. . \[th:theo7\]

In the case where the components $\kappa^{\mu\nu}$ represent the leading perturbation (see Eq. ), it may be seen that the expansion pattern in theorems \[th:theo6\] and \[th:theo7\] is almost the same as the one in theorems 4 and 5 of [@2008CQGra..25n5020T]. Actually, if one assumes that the gravitational components $k^{\mu\nu}$ are the leading perturbations, one obtains quasi similar theorems than theorems 4 and 5 of [@2008CQGra..25n5020T]. The difference would be in the definition of the quantities $\widehat{k}^{\mu\nu}_{+(l)}$. In our case, they would have involved the total time delay containing the gravitational, the refractive, and the coupling components, instead of the unique gravitational contribution.

The gravitational time delay functions {#sec:detgra}
--------------------------------------

Following [@2008CQGra..25n5020T], we suppose that the gravitational perturbation terms $h_{\mu\nu}$ can be expressed as a post-Minkowskian expansion, such like $$h_{\mu\nu}(x,G)=\sum_{l=1}^{\infty}h_{\mu\nu}^{(l)}(x)\text{.}
  \label{eq:hPM}$$ The contravariant components are given by $$k^{\mu\nu}(x,G)=\sum_{l=1}^{\infty}k^{\mu\nu}_{(l)}(x)
  \label{eq:kPM}$$ where the components $k^{\mu\nu}_{(l)}$ can be recursively determined using the following relationships

\[eq:kPMord\] $$\begin{aligned}
  k^{\mu\nu}_{(1)}&=-\eta^{\mu\alpha}\eta^{\beta\nu}h_{\alpha\beta}^{(1)}\text{,}\label{eq:kPMord1}\\
  k^{\mu\nu}_{(l)}&=-\eta^{\mu\alpha}\eta^{\beta\nu}h_{\alpha\beta}^{(l)}-\sum_{m=1}^{l-1}\eta^{\mu\alpha}h_{\alpha\beta}^{(m)}k^{\beta\nu}_{(l-m)}\label{eq:kPMord2}\end{aligned}$$

for $l\geqslant 2$.

The right-hand side of Eq.  involves terms such as $k^{\mu\nu}(z_-(\lambda))$. Thus, by making use of Eq. , we find $$k^{\mu\nu}(z_-(\lambda),G)=\sum_{l=1}^{\infty}k^{\mu\nu}_{(l)}(z_-(\lambda))
  \label{eq:kPMz-}$$ where $k^{\mu\nu}_{(l)}\propto G^l$. Therefore, the general expansion of $k^{\mu\nu}(\widetilde{z}_-(\lambda))$ is obtained after substituting for $\Delta$ and $k^{\mu\nu}$ from Eqs.  and into , respectively. After some algebra, we find the following expression $$k^{\mu\nu}(\widetilde{z}_-(\lambda),N_0,G)=\sum_{l=s}^{\infty}\widehat{k}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)
  \label{eq:kRPM}$$ where the quantities $\widehat{k}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)$ are given for $l\geqslant s$ by the following expressions

\[eq:kRPMdec\] $$\widehat{k}^{\mu\nu}_{-(s)}(z_-(\lambda),x_B)=k_{(1)}^{\mu\nu}(z_-(\lambda))\text{,}$$ and $$\begin{aligned}
  &\widehat{k}^{\mu\nu}_{-(ps)}(z_-(\lambda),x_B)=k^{\mu\nu}_{(p)}(z_-(\lambda))\nonumber\\
  &+\sum_{m=1}^{p-1}\sum_{n=1}^{ms}\Phi_-^{(ms,n)}(\mathbf z_-(\lambda),x_B)\left[\frac{\partial^nk_{(p-m)}^{\mu\nu}}{(\partial x^0)^n}\right]_{z_-(\lambda)}\label{eq:kRPMdecl}\end{aligned}$$ for $p\geqslant 2$, and $$\begin{aligned}
  &\widehat{k}^{\mu\nu}_{-(ps+q)}(z_-(\lambda),x_B)=\nonumber\\
  &+\sum_{n=1}^{q}\Phi_-^{(q,n)}(\mathbf z_-(\lambda),x_B)\left[\frac{\partial^nk_{(p)}^{\mu\nu}}{(\partial x^0)^n}\right]_{z_-(\lambda)}\nonumber\\
  &+\sum_{m=1}^{p-1}\sum_{n=1}^{ms+q}\Phi_-^{(ms+q,n)}(\mathbf z_-(\lambda),x_B)\left[\frac{\partial^nk_{(p-m)}^{\mu\nu}}{(\partial x^0)^n}\right]_{z_-(\lambda)}\label{eq:kRPMdecl}\end{aligned}$$

for $p\geqslant 1$ and $1\leqslant q\leqslant s-1$, where $p$ and $q$ are determined from $l$ using the following relationships $$\begin{aligned}
  p&=\lfloor l/s\rfloor\text{,} & q&=l-ps\text{,}\label{eq:pq}\end{aligned}$$ with $\lfloor i\rfloor$ denoting the integer part of $i$.

By substituting for $k^{\mu\nu}(\widetilde{z}_-(\lambda))$ from Eq. into , we infer the theorem which follows concerning the gravitational time delay function at reception.

In the optical spacetime, within the quasi-Minkowskian path regime, $\Delta$ admits the series expansion introduced in theorem \[th:theo4\], so the function $\Delta_{\mathrm{g}}$ is given by $$\Delta_{\mathrm{g}}(\mathbf{x}_A,x_B,N_0,G)=\sum_{l=1}^{\infty}\Delta_{\mathrm{g}}^{(l)}(\mathbf{x}_A,x_B)
      \label{eq:delRgPMN0G}$$ where

\[eq:intdelRgPM\] $$\begin{aligned}
      \Delta_{\mathrm{g}}^{(1)}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\big(k^{00}_{(1)}-2k^{0i}_{(1)}N^i_{AB}+k^{ij}_{(1)}N^i_{AB}N^j_{AB}\big)_{z_-(\lambda)}{\mathrm{d}}\lambda\text{,}\label{eq:intdelRgPM1}\\
      \Delta_{\mathrm{g}}^{(l)}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\big(\widehat k^{00}_{-(s+l-1)}-2\widehat k^{0i}_{-(s+l-1)}N^i_{AB}+\widehat k^{ij}_{-(s+l-1)}N^i_{AB}N^j_{AB}\big)_{(z_-(\lambda),x_B)}{\mathrm{d}}\lambda
    \end{aligned}$$ for $2\leqslant l\leqslant s$, and $$\begin{aligned}
      \Delta_{\mathrm{g}}^{(l)}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{k}^{00}_{-(s+l-1)}-2\widehat{k}^{0i}_{-(s+l-1)}N^i_{AB}+\widehat{k}^{ij}_{-(s+l-1)}N^i_{AB}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\nonumber\\
      &+2\sum_{m=s}^{l-1}\big(\widehat k^{0i}_{-(m)}-\widehat k^{ij}_{-(m)}N_{AB}^j\big)_{(z_-(\lambda),x_B)}\left[\frac{\partial\Delta_{\mathrm{g}}^{(l-m)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\eta^{ij}\sum_{m=1}^{l-s}\left[\frac{\partial\Delta_{\mathrm{g}}^{(m)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}^{(l-s-m+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $s+1\leqslant l\leqslant 2s$, and finally $$\begin{aligned}
      \Delta_{\mathrm{g}}^{(l)}(\mathbf{x}_A,x_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{k}^{00}_{-(s+l-1)}-2\widehat{k}^{0i}_{-(s+l-1)}N^i_{AB}+\widehat{k}^{ij}_{-(s+l-1)}N^i_{AB}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\nonumber\\
      &+2\sum_{m=s}^{l-1}\big(\widehat k^{0i}_{-(m)}-\widehat k^{ij}_{-(m)}N_{AB}^j\big)_{(z_-(\lambda),x_B)}\left[\frac{\partial\Delta_{\mathrm{g}}^{(l-m)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\eta^{ij}\sum_{m=1}^{l-s}\left[\frac{\partial\Delta_{\mathrm{g}}^{(m)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}^{(l-s-m+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=s}^{l-s-1}\big(\widehat{k}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-s-m}\left[\frac{\partial\Delta_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}^{(l-s-m-n+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$

for $l\geqslant 2s+1$. The quantities $\widehat{k}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)$ are defined in Eqs. . \[th:theo8\]

A similar reasoning allows us to state an analogous theorem for the gravitational time delay function at emission. Indeed, the right-hand side of Eq.  involves line integrals along the light path $\widetilde{z}_+(\mu)$ parametrized by $0\leqslant\mu\leqslant1$. After Taylor expanding the light path about the point-event $z_+(\mu)$, the right-hand side of Eq.  involves terms such as $k^{\mu\nu}(z_+(\mu))$. Thus, by making use of Eq. , we find $$k^{\mu\nu}(z_+(\mu),G)=\sum_{l=1}^{\infty}k^{\mu\nu}_{(l)}(z_+(\mu))\text{.}
  \label{eq:kPMz+}$$

Then, the general expansion of $k^{\mu\nu}(\widetilde{z}_+(\mu))$ is obtained after substituting for $\Xi$ from Eq. , and for $k^{\mu\nu}$ from Eq. , into . After some algebra, we find $$k^{\mu\nu}(\widetilde{z}_+(\mu),N_0,G)=\sum_{l=s}^{\infty}\widehat{k}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))
  \label{eq:kEPM}$$ where the quantities $\widehat{k}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))$ are given for $l\geqslant s$ by the following expressions

\[eq:kEPMdec\] $$\widehat{k}^{\mu\nu}_{+(s)}(x_A,z_+(\mu))=k_{(1)}^{\mu\nu}(z_+(\mu))\text{,}$$ and $$\begin{aligned}
  &\widehat{k}^{\mu\nu}_{+(ps)}(x_A,z_+(\mu))=k^{\mu\nu}_{(p)}(z_+(\mu))\nonumber\\
  &+\sum_{m=1}^{p-1}\sum_{n=1}^{ms}\Phi_+^{(ms,n)}(x_A,\mathbf z_+(\mu))\left[\frac{\partial^nk_{(p-m)}^{\mu\nu}}{(\partial x^0)^n}\right]_{z_+(\mu)}\label{eq:kEPMdecl}\end{aligned}$$ for $p\geqslant 2$, and $$\begin{aligned}
  &\widehat{k}^{\mu\nu}_{+(ps+q)}(x_A,z_+(\mu))=\nonumber\\
  &+\sum_{n=1}^{q}\Phi_+^{(q,n)}(x_A,\mathbf z_+(\mu))\left[\frac{\partial^nk_{(p)}^{\mu\nu}}{(\partial x^0)^n}\right]_{z_+(\mu)}\nonumber\\
  &+\sum_{m=1}^{p-1}\sum_{n=1}^{ms+q}\Phi_+^{(ms+q,n)}(x_A,\mathbf z_+(\mu))\left[\frac{\partial^nk_{(p-m)}^{\mu\nu}}{(\partial x^0)^n}\right]_{z_+(\mu)}\label{eq:kEPMdecl}\end{aligned}$$

\
for $p\geqslant 1$ and $1\leqslant q\leqslant s-1$ where $p$ and $q$ are determined from $l$ using the relationships in Eq. .

By substituting for $k^{\mu\nu}(\widetilde{z}_-(\lambda))$ from Eq. into , we infer the theorem which follows.

In the optical spacetime, within the quasi-Minkowskian path regime, $\Xi$ admits the series expansion introduced in theorem \[th:theo5\], so the function $\Xi_{\mathrm{g}}$ is given by $$\Xi_{\mathrm{g}}(x_A,\mathbf{x}_B,N_0,G)=\sum_{l=1}^{\infty}\Xi_{\mathrm{g}}^{(l)}(x_A,\mathbf{x}_B)
      \label{eq:delEgPM}$$ where

\[eq:intdelEgPM\] $$\begin{aligned}
      \Xi_{\mathrm{g}}^{(1)}(x_A,\mathbf{x}_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\big(k^{00}_{(1)}-2k^{0i}_{(1)}N^i_{AB}+k^{ij}_{(1)}N^i_{AB}N^j_{AB}\big)_{z_+(\mu)}{\mathrm{d}}\mu\text{,}\\
      \Xi_{\mathrm{g}}^{(l)}(x_A,\mathbf{x}_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\big(\widehat k^{00}_{+(s+l-1)}-2\widehat k^{0i}_{+(s+l-1)}N^i_{AB}+\widehat k^{ij}_{+(s+l-1)}N^i_{AB}N^j_{AB}\big)_{(x_A,z_+(\mu))}{\mathrm{d}}\mu
    \end{aligned}$$ for $2\leqslant l\leqslant s$, and $$\begin{aligned}
      \Xi_{\mathrm{g}}^{(l)}(x_A,\mathbf{x}_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{k}^{00}_{+(s+l-1)}-2\widehat{k}^{0i}_{+(s+l-1)}N^i_{AB}+\widehat{k}^{ij}_{+(s+l-1)}N^i_{AB}N^j_{AB}\big)_{(x_A,z_+(\mu))}\nonumber\\
      &-2\sum_{m=s}^{l-1}\big(\widehat k^{0i}_{+(m)}-\widehat k^{ij}_{+(m)}N_{AB}^j\big)_{(x_A,z_+(\mu))}\left[\frac{\partial\Xi_{\mathrm{g}}^{(l-m)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\eta^{ij}\sum_{m=1}^{l-s}\left[\frac{\partial\Xi_{\mathrm{g}}^{(m)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}^{(l-s-m+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$ for $s+1\leqslant l\leqslant 2s$, and finally $$\begin{aligned}
      \Xi_{\mathrm{g}}^{(l)}(x_A,\mathbf{x}_B)&=\frac{R_{AB}}{2}\int^{1}_{0}\Bigg\{\big(\widehat{k}^{00}_{+(s+l-1)}-2\widehat{k}^{0i}_{+(s+l-1)}N^i_{AB}+\widehat{k}^{ij}_{+(s+l-1)}N^i_{AB}N^j_{AB}\big)_{(x_A,z_+(\mu))}\nonumber\\
      &-2\sum_{m=s}^{l-1}\big(\widehat k^{0i}_{+(m)}-\widehat k^{ij}_{+(m)}N_{AB}^j\big)_{(x_A,z_+(\mu))}\left[\frac{\partial\Xi_{\mathrm{g}}^{(l-m)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\eta^{ij}\sum_{m=1}^{l-s}\left[\frac{\partial\Xi_{\mathrm{g}}^{(m)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}^{(l-s-m+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\sum_{m=s}^{l-s-1}\big(\widehat{k}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-s-m}\left[\frac{\partial\Xi_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}^{(l-s-m-n+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\lambda))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$

for $l\geqslant 2s+1$. The quantities $\widehat{k}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))$ are defined in Eqs. . \[th:theo9\]

As a final remark, let us emphasize that Eqs.  and are independent of the total delay function. Therefore, as mentioned previously, they can directly be used in Eq. for the determination of $s$.

The coupling time delay functions {#sec:detcou}
---------------------------------

All the basic ingredients needed for the establishment of the general expansion of the coupling time delay functions, have been introduced in Secs. \[sec:detref\] and \[sec:detgra\]. The general expansions of the refractive and gravitational spacetime perturbations are given in Eqs. and , respectively. Then, the expansions of the reception time delay functions can be found in Eqs. , , and .

Therefore, by substituting for $\kappa^{\mu\nu}(\widetilde{z}_-(\lambda))$ and $k^{\mu\nu}(\widetilde{z}_-(\lambda))$ from Eqs.  and into , respectively, we obtain the theorem which follows.

In the optical spacetime, within the quasi-Minkowskian path regime, $\Delta$ admits the series expansion introduced in theorem \[th:theo4\], so the function $\Delta_{\mathrm{gr}}$ is given by $$\Delta_{\mathrm{gr}}(\mathbf{x}_A,x_B,N_0,G)=\sum_{l=1}^{\infty}\Delta_{\mathrm{gr}}^{(l)}(\mathbf{x}_A,x_B)
      \label{eq:delRgrPMN0G}$$ where

\[eq:intdelRgrPM\] $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant 1)}(\mathbf{x}_A,x_B)&=R_{AB}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l}\left[\frac{\partial\Delta_{\mathrm{r}}^{(m)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}^{(l-m+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=1}^{l}\big(\widehat{\kappa}^{0i}_{-(m)}-\widehat{\kappa}^{ij}_{-(m)}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\left[\frac{\partial\Delta_{\mathrm{g}}^{(l-m+1)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=s}^{l+s-1}\big(\widehat{k}^{0i}_{-(m)}-\widehat{k}^{ij}_{-(m)}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\left[\frac{\partial\Delta_{\mathrm{r}}^{(l+s-m)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $l\geqslant 1$, and $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant 2)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{gr}}^{(l\geqslant 1)}(\mathbf{x}_A,x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l-1}\left[\frac{\partial\Delta_{\mathrm{r}}^{(m)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-m)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=1}^{l-1}\big(\widehat{\kappa}^{0i}_{-(m)}-\widehat{\kappa}^{ij}_{-(m)}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\left[\frac{\partial\Delta_{\mathrm{gr}}^{(l-m)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=1}^{l-1}\big(\widehat{\kappa}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-m}\left[\frac{\partial\Delta_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}^{(l-m-n+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\frac{1}{2}\sum_{m=s}^{l+s-2}\big(\widehat{k}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l+s-m-1}\left[\frac{\partial\Delta_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{r}}^{(l+s-m-n)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $l\geqslant 2$, and $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant 3)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{gr}}^{(l\geqslant 2)}(\mathbf{x}_A,x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\sum_{m=1}^{l-2}\big(\widehat{\kappa}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-m-1}\left[\frac{\partial\Delta_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-m-n)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $l\geqslant 3$, and $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant s+1)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{gr}}^{(l\geqslant s)}(\mathbf{x}_A,x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l-s}\left[\frac{\partial\Delta_{\mathrm{g}}^{(m)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-s-m+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=s}^{l-1}\big(\widehat{k}^{0i}_{-(m)}-\widehat{k}^{ij}_{-(m)}N^j_{AB}\big)_{(z_-(\lambda),x_B)}\left[\frac{\partial\Delta_{\mathrm{gr}}^{(l-m)}}{\partial x^i}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\sum_{m=s}^{l-1}\big(\widehat{k}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-m}\left[\frac{\partial\Delta_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}^{(l-m-n+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+\frac{1}{2}\sum_{m=1}^{l-s}\big(\widehat{\kappa}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-s-m+1}\left[\frac{\partial\Delta_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{g}}^{(l-s-m-n+2)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $l\geqslant s+1$, and $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant s+2)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{gr}}^{(l\geqslant s+1)}(\mathbf{x}_A,x_B)\nonumber\\
      &+\frac{R_{AB}}{2}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l-s-1}\left[\frac{\partial\Delta_{\mathrm{gr}}^{(m)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-s-m)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+2\sum_{m=s}^{l-2}\big(\widehat{k}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-m-1}\left[\frac{\partial\Delta_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-m-n)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\nonumber\\
      &+2\sum_{m=1}^{l-s-1}\big(\widehat{\kappa}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-s-m}\left[\frac{\partial\Delta_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-s-m-n+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $l\geqslant s+2$, and $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant s+3)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{gr}}^{(l\geqslant s+2)}(\mathbf{x}_A,x_B)\nonumber\\
      &+\frac{R_{AB}}{2}\int_0^1\Bigg\{\sum_{m=1}^{l-s-2}\big(\widehat{\kappa}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-s-m-1}\left[\frac{\partial\Delta_{\mathrm{gr}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-s-m-n)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $l\geqslant s+3$, and $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant 2s+1)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{gr}}^{(l\geqslant 2s)}(\mathbf{x}_A,x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\sum_{m=s}^{l-s-1}\big(\widehat{k}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-s-m}\left[\frac{\partial\Delta_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-s-m-n+1)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$ for $l\geqslant 2s+1$, and finally $$\begin{aligned}
      \Delta_{\mathrm{gr}}^{(l\geqslant 2s+2)}(\mathbf{x}_A,x_B)&=\Delta_{\mathrm{gr}}^{(l\geqslant 2s+1)}(\mathbf{x}_A,x_B)\nonumber\\
      &+\frac{R_{AB}}{2}\int_0^1\Bigg\{\sum_{m=s}^{l-s-2}\big(\widehat{k}^{ij}_{-(m)}\big)_{(z_-(\lambda),x_B)}\sum_{n=1}^{l-s-m-1}\left[\frac{\partial\Delta_{\mathrm{gr}}^{(n)}}{\partial x^i}\frac{\partial\Delta_{\mathrm{gr}}^{(l-s-m-n)}}{\partial x^j}\right]_{(\mathbf{z}_-(\lambda),x_B)}\Bigg\}{\mathrm{d}}\lambda
    \end{aligned}$$

for $l\geqslant 2s+2$. The quantities $\widehat{k}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)$ and $\widehat{\kappa}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)$ are defined in Eqs. and , respectively. \[th:theo10\]

Applying the exact same reasoning, the analogous theorem for the coupling time delay function at emission can be stated. Indeed, substituting for $\kappa^{\mu\nu}(\widetilde{z}_+(\mu))$ from Eq. and for $k^{\mu\nu}(\widetilde{z}_+(\mu))$ from Eq. into , we obtain the theorem which follows.

In the optical spacetime, within the quasi-Minkowskian path regime, $\Xi$ admits the series expansion introduced in theorem \[th:theo5\], so the function $\Xi_{\mathrm{gr}}$ is given by $$\Xi_{\mathrm{gr}}(x_A,\mathbf x_B,N_0,G)=\sum_{l=1}^{\infty}\Xi_{\mathrm{gr}}^{(l)}(x_A,\mathbf x_B)
      \label{eq:delEgrPM}$$ where

\[eq:intdelEgrPM\] $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant 1)}(x_A,\mathbf x_B)&=R_{AB}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l}\left[\frac{\partial\Xi_{\mathrm{r}}^{(m)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}^{(l-m+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &-\sum_{m=1}^{l}\big(\widehat{\kappa}^{0i}_{+(m)}-\widehat{\kappa}^{ij}_{+(m)}N^j_{AB}\big)_{(x_A,z_+(\mu))}\left[\frac{\partial\Xi_{\mathrm{g}}^{(l-m+1)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &-\sum_{m=s}^{l+s-1}\big(\widehat{k}^{0i}_{+(m)}-\widehat{k}^{ij}_{+(m)}N^j_{AB}\big)_{(x_A,z_+(\mu))}\left[\frac{\partial\Xi_{\mathrm{r}}^{(l+s-m)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$ for $l\geqslant 1$, and $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant 2)}(x_A,\mathbf x_B)&=\Xi_{\mathrm{gr}}^{(l\geqslant 1)}(x_A,\mathbf x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l-1}\left[\frac{\partial\Xi_{\mathrm{r}}^{(m)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-m)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &-\sum_{m=1}^{l-1}\big(\widehat{\kappa}^{0i}_{+(m)}-\widehat{\kappa}^{ij}_{+(m)}N^j_{AB}\big)_{(x_A,z_+(\mu))}\left[\frac{\partial\Xi_{\mathrm{gr}}^{(l-m)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\sum_{m=1}^{l-1}\big(\widehat{\kappa}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-m}\left[\frac{\partial\Xi_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}^{(l-m-n+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\frac{1}{2}\sum_{m=s}^{l+s-2}\big(\widehat{k}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l+s-m-1}\left[\frac{\partial\Xi_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{r}}^{(l+s-m-n)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$ for $l\geqslant 2$, and $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant 3)}(x_A,\mathbf x_B)&=\Xi_{\mathrm{gr}}^{(l\geqslant 2)}(x_A,\mathbf x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\sum_{m=1}^{l-2}\big(\widehat{\kappa}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-m-1}\left[\frac{\partial\Xi_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-m-n)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$ for $l\geqslant 3$, and $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant s+1)}(x_A,\mathbf x_B)&=\Xi_{\mathrm{gr}}^{(l\geqslant s)}(x_A,\mathbf x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l-s}\left[\frac{\partial\Xi_{\mathrm{g}}^{(m)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-s-m+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &-\sum_{m=s}^{l-1}\big(\widehat{k}^{0i}_{+(m)}-\widehat{k}^{ij}_{+(m)}N^j_{AB}\big)_{(x_A,z_+(\mu))}\left[\frac{\partial\Xi_{\mathrm{gr}}^{(l-m)}}{\partial x^i}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\sum_{m=s}^{l-1}\big(\widehat{k}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-m}\left[\frac{\partial\Xi_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}^{(l-m-n+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+\frac{1}{2}\sum_{m=1}^{l-s}\big(\widehat{\kappa}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-s-m+1}\left[\frac{\partial\Xi_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{g}}^{(l-s-m-n+2)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$ for $l\geqslant s+1$, and $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant s+2)}(x_A,\mathbf x_B)&=\Xi_{\mathrm{gr}}^{(l\geqslant s+1)}(x_A,\mathbf x_B)\nonumber\\
      &+\frac{R_{AB}}{2}\int_0^1\Bigg\{\eta^{ij}\sum_{m=1}^{l-s-1}\left[\frac{\partial\Xi_{\mathrm{gr}}^{(m)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-s-m)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+2\sum_{m=s}^{l-2}\big(\widehat{k}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-m-1}\left[\frac{\partial\Xi_{\mathrm{r}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-m-n)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\nonumber\\
      &+2\sum_{m=1}^{l-s-1}\big(\widehat{\kappa}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-s-m}\left[\frac{\partial\Xi_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-s-m-n+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$ for $l\geqslant s+2$, and $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant s+3)}(x_A,\mathbf x_B)&=\Xi_{\mathrm{gr}}^{(l\geqslant s+2)}(x_A,\mathbf x_B)\nonumber\\
      &+\frac{R_{AB}}{2}\int_0^1\Bigg\{\sum_{m=1}^{l-s-2}\big(\widehat{\kappa}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-s-m-1}\left[\frac{\partial\Xi_{\mathrm{gr}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-s-m-n)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu\text{,}
    \end{aligned}$$ for $l\geqslant s+3$, and $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant 2s+1)}(x_A,\mathbf x_B)&=\Xi_{\mathrm{gr}}^{(l\geqslant 2s)}(x_A,\mathbf x_B)\nonumber\\
      &+R_{AB}\int_0^1\Bigg\{\sum_{m=s}^{l-s-1}\big(\widehat{k}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-s-m}\left[\frac{\partial\Xi_{\mathrm{g}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-s-m-n+1)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$ for $l\geqslant 2s+1$, and finally $$\begin{aligned}
      \Xi_{\mathrm{gr}}^{(l\geqslant 2s+2)}(x_A,\mathbf x_B)&=\Xi_{\mathrm{gr}}^{(l\geqslant 2s+1)}(x_A,\mathbf x_B)\nonumber\\
      &+\frac{R_{AB}}{2}\int_0^1\Bigg\{\sum_{m=s}^{l-s-2}\big(\widehat{k}^{ij}_{+(m)}\big)_{(x_A,z_+(\mu))}\sum_{n=1}^{l-s-m-1}\left[\frac{\partial\Xi_{\mathrm{gr}}^{(n)}}{\partial x^i}\frac{\partial\Xi_{\mathrm{gr}}^{(l-s-m-n)}}{\partial x^j}\right]_{(x_A,\mathbf{z}_+(\mu))}\Bigg\}{\mathrm{d}}\mu
    \end{aligned}$$

for $l\geqslant 2s+2$. The quantities $\widehat{k}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))$ and $\widehat{\kappa}^{\mu\nu}_{+(l)}(x_A,z_+(\mu))$ are defined in Eqs. and , respectively. \[th:theo11\]

Finally, from the gravitational, the refractive, and the coupling components, we can now determine the time delay expression up to $l$th-order by applying theorem \[th:theo4\]. Then, the expressions for the range and the time transfer functions are determined from Eqs. and , respectively.

Let us emphasis that line integrals occurring in Eqs. , , and are now zeroth-order null geodesics with parametric equations $x=z_-(\lambda)$. Similarly, Eqs. , , and are integrated along the zeroth-order null geodesic path with parametric equations $x=z_+(\mu)$. This specificity of the time transfer functions formalism considerably simplifies the integrations and constitutes one of the most important advantage with respect to an explicit resolution of the null geodesic equation.

The usefulness of the decomposition which has been performed in Eq.  is really apparent in stationary optical spacetimes. Indeed, when the coordinates $(x^\mu)$ are chosen so that the optical spacetime metric does not depend on $x^0$, it is seen that the series expansion in Eqs.  and reduce to

\[eq:kapkRPMsta\] $$\begin{aligned}
  \kappa^{\mu\nu}(\widetilde{z}_-(\lambda),N_0)&=\sum_{l=1}^{\infty}\widehat{\kappa}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)\text{,}\label{eq:kapRPMsta}\\
  k^{\mu\nu}(\widetilde{z}_-(\lambda),G)&=\sum_{l=1}^{\infty}\widehat{k}^{\mu\nu}_{-(l)}(z_-(\lambda),x_B)\label{eq:kRPMsta}\end{aligned}$$

where

\[eq:kapkRPMdecsta\] $$\begin{aligned}
  \widehat{\kappa}^{\mu\nu}_{(l)}(z_-(\lambda),x_B)&=\kappa^{\mu\nu}_{(l)}(z_-(\lambda))\text{,}\label{eq:kapRPMdecsta}\\
  \widehat{k}^{\mu\nu}_{(ps+q)}(z_-(\lambda),x_B)&=\delta(q)\,k^{\mu\nu}_{(p)}(z_-(\lambda))\text{,}\label{eq:kRPMdecsta}\end{aligned}$$

respectively. We recall that $p$ and $q$ are determined from $l$ using .

Hence, Eq. does no longer involve the formal expansion parameter $G$, while is now independent of $N_0$. In other words, the different theorems can be solved separately, each one being completely independent from the others. As a matter of fact, theorems involving gravitational perturbation become independent of $\widehat{k}^{\mu\nu}_{(l)}$ for any $l$ which is not a multiple of $s$.

Application to stationary optical spacetime in geocentric celestial reference system {#sec:app}
====================================================================================

Let us now illustrate the method by determining the time transfer function up to the post-linear approximation. We investigate, the light-dragging effect experienced by a signal during its propagation inside a flowing media of non-null refractivity. In the GCRS, the effect shows up at the post-linear approximation. In the case where the motion of the Earth’s atmosphere is mainly a steady rotation (e.g. in GCRS), we show that the light-dragging effect reduce to a geometrical factor scaling the static atmospheric contribution. During the computation, we never make use of an index of refraction profile in order to keep the equations under a form which is as general as possible.

Notations and definitions
-------------------------

We consider that spacetime is covered with some global coordinates $(x^\mu)$. We choose the coordinate system such that the optical metric components are independent of $x^0$. In addition, the coordinate system shall be chosen in such a way that it is convenient to model the outcomes of an experiment taking place in the Earth’s close vicinity. Therefore, we consider that $(x^\mu)$ are the GCRS coordinates. We recall that the GCRS is centered in the Earth’s center of mass and is non-rotating with respect to distant stars. We suppose that the domain $\mathcal{D}$ defines the spacetime boundaries of the Earth’s neutral atmosphere. In that sense, $\mathcal{D}$ draws a time-like tube in spacetime. The Earth’s atmosphere is considered spherically symmetric and we suppose that it is filled with a fluid dielectric medium whose refractive properties are independent of the component $x^0$, that is to say $$n(\mathbf{x})= 1+N(\mathbf{x})\text{.}$$ We consider that the atmosphere is still in the reference system rotating with the Earth, thus we assume that the unit 4-velocity vector is given in GCRS by $$w^{\mu}=w^0(1,\xi^i)
  \label{eq:4velmed}$$ where $\xi^i$ is the coordinate 3-velocity vector of the fluid dielectric medium. Hereafter, we assume that the 3-velocity vector is given by the following expression $$\xi^i(\mathbf x)=\frac{\omega_\oplus}{c}\,e^{ijk}S_\oplus^jx^k  \label{eq:3velmed}$$ where $\omega_\oplus$ is the magnitude of the Earth’s angular velocity of rotation and $\mathbf S_\oplus$ is the direction of the spin axis.

Moreover, we consider the case of a one-way transfer, with the transmitter being right outside $\mathcal{D}$ and the receiver being comoving with the fluid dielectrics medium, that is to say at rest in the reference system rotating with the Earth. To fix ideas for future discussion, let us assume that the emitter is transmitting from the International Space Station (ISS) at an altitude of $h\simeq 400\,\mathrm{km}$. Furthermore, let us consider that the emitter is moving along the time-like worldline $\mathcal{C}_A$ with the unit 4-velocity vector $\bm u_{A}$ defined by $$u^{\mu}_A=u^0_A(1,\beta^i_A)
  \label{eq:4velA}$$ where $\beta^i_A$ is the coordinate 3-velocity vector expressed in GCRS coordinates. Similarly, we assume that the receiver moves along the time-like worldline $\mathcal{C}_B$ with the unit 4-velocity vector $\bm u_{B}$ defined by $$u^{\mu}_B=u^0_B(1,\beta^i_B)
  \label{eq:4velB}$$ where $\beta^i_B$ is the coordinate 3-velocity vector expressed in GCRS coordinates. For a receiver comoving with the medium, we have $$\beta^i_B=\xi^i(\mathbf x_B)\text{.}$$

Expansion of the delay functions
--------------------------------

The components of the physical spacetime metric expressed in GCRS coordinates are given in [@2003AJ....126.2687S]. Let us emphasize that i) our convention for the signature of spacetime is $(+,-,-,-)$, and ii) our $g_{\mu\nu}$ corresponds to $G_{\alpha\beta}$ in [@2003AJ....126.2687S]. By keeping terms in $1/c^2$ only, the first-order gravitational perturbation reads as follows

\[eq:gGCRS\] $$h_{00}^{(1)}=-\frac{2U}{c^2}\text{,} \qquad h_{0i}^{(1)}=0\text{,} \qquad h_{ij}^{(1)}=-\frac{2U}{c^2}\delta_{ij}
    \label{eq:gGCRScov}$$ where the contravariant components are determined from Eq. $$k^{00}_{(1)}=\frac{2U}{c^2}\text{,} \qquad k^{0i}_{(1)}=0\text{,} \qquad k^{ij}_{(1)}=\frac{2U}{c^2}\delta_{ij}\text{.}
    \label{eq:gGCRScon}$$

In these expressions, we restrict $U$ to the monopole term of the Newtonian gravitational potential of the Earth, that is $$U(\mathbf x)=\frac{Gm_\oplus}{|\mathbf x|}
  \label{eq:U}$$ where $m_\oplus$ is the mass of the Earth. In that respect, at the level of the surface of the Earth, we find $$|k^{\mu\nu}_{(1)}(R_\oplus)|_{\mathrm{max}}\propto \frac{U(R_\oplus)}{c^2}\sim10^{-10}
  \label{eq:ordk}$$ where $R_\oplus$ denotes Earth’s equatorial radii.

According to [@doi1010292010JD015214], at the sea level an average parcel of air possesses a refractivity $N(R_\oplus)\simeq 3\times10^{-4}$, so we consider $N_0\sim10^{-4}$. Additionally, at the Earth’s surface, the 3-velocity of the refractive medium expressed in GCRS coordinates is $|\xi^i(R_\oplus)|_{\mathrm{max}}\propto\omega_\oplus R_\oplus/c\sim 10^{-6}$. Consequently, we can expand the refractive perturbation in term of the refractivity at the Earth’s surface and in the approximation of small velocities. Therefore, it can be seen that the first-order term of the refractive perturbation is given by (see Eqs. ) $$|\kappa^{\mu\nu}_{(1)}(R_\oplus)|_{\mathrm{max}}\propto N_0\sim 10^{-4}
  \label{eq:ordkap}$$ once evaluated at the Earth’s surface.

At the same time, typical measurement profile for the neutral atmosphere using GPS/MET (Global Positioning System Meteorology) occultations data starts at $\ell\simeq100\,\mathrm{km}$ [@1999AnGeo..17..122S], so that $\ell/h\simeq 0.4$. For observations at lower elevation than the zenith direction, we can roughly take $\ell/R_{AB}\sim0.1$. Then, if we consider that the light path is sufficiently small so that the metric components vary slowly during the integration, we can get a rough estimation of $s$ by making use of Eqs. . We quickly infer that $s$ must satisfy the zeroth-order following relation $$\frac{\ell}{R_{AB}}|\kappa^{\mu\nu}|_{\mathrm{max}}\sim |k^{\mu\nu}|_{\mathrm{max}}^{1/s}\text{.}
  \label{eq:relmag}$$ Inserting numerical values, we deduce $s=2$. This results can be double-checked by inserting the first-order expressions of the gravitational and refractive delays (see Eqs.  and ) into Eq. .

In this application, we exclude third-order terms and beyond, that is to say, all terms of the order of $\varepsilon^3$ with $\varepsilon\sim\ell/R_{AB}|\kappa^{\mu\nu}|_{\mathrm{max}}\sim 10^{-5}$. This means that a post-linear expression of the range transfer function neglects terms of the order $\varepsilon^3R_{AB}$. Therefore, the coupling terms which are of third-order are neglected too.

A look at Eq.  and allows one to infer that the time component of the 4-velocity vector of the fluid dielectric medium must be known up to $10^{-5}$ in order to account for all second-order terms. Considering that the four-velocity of the medium must be a unit-vector for the spacetime metric $g_{\mu\nu}$, we have the relation $$w^0=\left(g_{00}+2g_{0i}\xi^i+g_{ij}\xi^i\xi^j\right)^{-1/2}\text{.}
  \label{eq:u0}$$ Therefore, to sufficient accuracy, we can safely consider for the rest of the application that $w^0=1$.

Hence, we end up with the following contravariant components for the refractive perturbation

\[eq:kapGCRS\] $$\kappa^{00}_{(1)}=2N\text{,} \qquad \kappa^{0i}_{(1)}=0\text{,} \qquad \kappa^{ij}_{(1)}=0\text{,}
    \label{eq:kapGCRS1}$$ with the second-order $$\kappa^{00}_{(2)}=N^2\text{,} \qquad \kappa_{(2)}^{0i}=2N\xi^i\text{,} \qquad \kappa_{(2)}^{ij}=0\text{.}
    \label{eq:kapGCRS2}$$

Let us note that the cross component is non-null at the post-linear approximation. It represents the light-dragging effect due to the motion of the fluid dielectric medium in GCRS coordinates.

Additionally, let us mention that the optical spacetime is stationary as seen from Eqs.  and . In that respect, the emission or the reception time transfer functions become identical. As a consequence, the distinction between emission and reception functions is not relevant anymore meaning that the time component at emission or reception is no longer an independent variable [@2008CQGra..25n5020T]. Hence, $\widehat \kappa_{(l)}^{\mu\nu}$ and $\widehat k_{(l)}^{\mu\nu}$ are now given by Eqs. and which are independent of the total time delay. Therefore, the refractive and the gravitational delays may be solved separately.

A straightforward application of theorem \[th:theo4\] assuming $s=2$, allows us to infer the following form of the expansion of the total time delay function

\[eq:ord2\] $$\begin{aligned}
    \Delta^{(1)}(\mathbf{x}_A,\mathbf x_B)&=\Delta^{(1)}_{\mathrm r}(\mathbf{x}_A,\mathbf x_B)\text{,}\\
    \Delta^{(2)}(\mathbf{x}_A,\mathbf x_B)&=\Delta^{(2)}_{\mathrm r}(\mathbf{x}_A,\mathbf x_B)+\Delta^{(1)}_{\mathrm g}(\mathbf{x}_A,\mathbf x_B)\text{.}
  \end{aligned}$$

Thus, we deduce that the contributions in Eq. are given by

$$\begin{aligned}
  \Delta_{\mathrm g}(\mathbf{x}_A,\mathbf x_B)&=\Delta^{(1)}_{\mathrm g}(\mathbf{x}_A,\mathbf x_B)\text{,}\\
  \Delta_{\mathrm r}(\mathbf{x}_A,\mathbf x_B)&=\Delta^{(1)}_{\mathrm r}(\mathbf{x}_A,\mathbf x_B)+\Delta^{(2)}_{\mathrm r}(\mathbf{x}_A,\mathbf x_B)\text{.}\end{aligned}$$

Then, theorems \[th:theo6\] and \[th:theo8\] together with Eqs. and allow us to determine the refractive and the gravitational contributions up to the appropriate order.

Time transfer function and Doppler
----------------------------------

Using the fact that spacetime is stationary, we first focus on the gravitational time delay. By making use of theorem \[th:theo8\], we soon arrive to the well-known formula $$\Delta_{\mathrm g}(\mathbf{x}_A,\mathbf x_B)=\frac{2R_{AB}}{c^2}\int_0^1U(\mathbf z_-(\lambda))\,{\mathrm{d}}\lambda\text{,}$$ which leads after integration to the Shapiro delay [@1964PhRvL..13..789S] $$\Delta_{\mathrm g}(\mathbf{x}_A,\mathbf x_B)=\frac{2Gm_\oplus}{c^2}\mathrm{ln}\left(\frac{r_A+r_B+R_{AB}}{r_A+r_B-R_{AB}}\right)\text{.}\label{eq:delShap}$$ We introduced the notations $r_{A/B}=|\mathbf{x}_{A/B}|$.

The first-order refractive contribution is derived from theorem \[th:theo6\] and is given by $$\Delta^{(1)}_{\mathrm{r}}(\mathbf{x}_A,\mathbf x_B)=R_{AB}\int_0^1N(\mathbf z_-(\lambda))\,{\mathrm{d}}\lambda\text{.}
  \label{eq:delr1GCRS}$$ We find almost similar expression for the atmospheric delay in [@inproceedingsMendes; @2002GeoRL..29.1414M; @2004GeoRL..3114602M; @2010ITN....36....1P]. When applied to the Earth’s neutral atmosphere, it is common to define the refractivity within a factor $10^6$ and to split it into a hydrostatic and a non-hydrostatic components [@2010ITN....36....1P].

The main differences stand, firstly, in the integration path which is performed along the Euclidean path between the emitter and the receiver in Eq. . This holds true even for non-zenithal observations. Instead, in the literature the atmospheric delay is usually computed at zenith, and then, mapping functions are used to convert the zenithal delay into a delay in the line-of-sight direction as discussed in [@2002GeoRL..29.1414M]. In our case, the first-order refractive delay is the well-known excess path delay due to the change of the phase velocity experienced by the signal during the crossing of the dielectric medium. The geometric delay due to the refractive bending of the ray arises at the post-linear order as we shall see in the next.

The other difference stands in the upper limit of integration path in Eq. . However, considering that the refractive region is bounded to a certain domain $\mathcal{D}$ of spacetime, the integration out of $\mathcal{D}$ results in integrating null refractivity which thus does not contribute to the final results. In that respect, the difference in the upper integration limit is only superficial.

According to theorem \[th:theo6\], the second-order is given by $$\begin{aligned}
  \Delta^{(2)}_{\mathrm r}(\mathbf{x}_A,\mathbf x_B)&=\frac{R_{AB}}{2}\int_0^1\Big\{\big(N^2-4N\xi^iN_{AB}^i\big)_{\mathbf z_-(\lambda)}\nonumber\\
  &-\big[\partial_i\Delta^{(1)}_{\mathrm r}\partial_i\Delta^{(1)}_{\mathrm r}\big]_{(\mathbf z_-(\lambda),\mathbf x_B)}\Big\}{\mathrm{d}}\lambda\label{eq:delr2GCRS}\end{aligned}$$ where $\partial_i\Delta^{(1)}_{\mathrm r}$ is computed by differentiating Eq. $$\begin{aligned}
  \big[\partial_i\Delta^{(1)}_{\mathrm r}\big]_{(\mathbf x,\mathbf x_B)}&=-\frac{(\mathbf x_B-\mathbf x)^i}{|\mathbf x_B-\mathbf x|}\int_0^1N(\mathbf y_-(\mu,\mathbf x))\,{\mathrm{d}}\mu\nonumber\\
  &+|\mathbf x_B-\mathbf x|\int_0^1\mu\big[\partial_i N\big]_{\mathbf y_-(\mu,\mathbf x)}{\mathrm{d}}\mu\text{.}\label{eq:DPdel}\end{aligned}$$ We have introduced $$\mathbf y_-(\mu,\mathbf x)=(1-\mu)\,\mathbf x_B+\mu\,\mathbf x\text{,}$$ which, in the case where $\mathbf x=\mathbf z_-(\lambda)$, reduces to $$\mathbf y_-(\mu,\mathbf z_-(\lambda))=\mathbf z_-(\mu\lambda)\text{.}
  \label{eq:ym}$$

We can rearrange Eq. by first noticing that the light-dragging contribution can be further simplified. Indeed, after making use of Eq. , it may be seen that $$\big(\xi^iN_{AB}^i\big)_{\mathbf{z}_-(\lambda)}\equiv\xi^i(\mathbf x_B)N_{AB}^i$$ which is obviously independent of $\lambda$. Then, substituting for $\partial_i\Delta^{(1)}_{\mathrm r}$ from Eq. into while accounting for and , one can apply the following change of variables $\mu'=\mu\lambda$. Finally, by integrating by parts the double integrals, we infer the post-linear refractive order $$\begin{aligned}
  \Delta^{(2)}_{\mathrm{r}}(\mathbf{x}_A,\mathbf x_B)&=\Delta^{(2)}_{\mathrm{r,exc}}(\mathbf{x}_A,\mathbf x_B)+\Delta^{(2)}_{\mathrm{r,geo}}(\mathbf{x}_A,\mathbf x_B)\nonumber\\
  &+\Delta^{(2)}_{\mathrm{r,drag}}(\mathbf{x}_A,\mathbf x_B)\label{eq:delr2}\end{aligned}$$ where we have introduced

\[eq:delr2GCRSbis\] $$\begin{aligned}
  \Delta^{(2)}_{\mathrm{r,exc}}&(\mathbf{x}_A,\mathbf x_B)=\frac{R_{AB}}{2}\int_0^1N^2(\mathbf z_-(\lambda))(1+\mathrm{ln}\lambda)\,{\mathrm{d}}\lambda\text{,}\label{eq:delr2GCRSexc}\\
  \Delta^{(2)}_{\mathrm{r,geo}}&(\mathbf{x}_A,\mathbf x_B)=\frac{R_{AB}}{2}\int_0^1\Big\{\lambda R_{AB}^2\big[\partial_i N\partial_i N\big]_{\mathbf z_-(\lambda)}\nonumber\\
  &-2R_{AB}N_{AB}^i \big[N\partial_i N\big]_{\mathbf z_-(\lambda)}\Big\}\lambda\mathrm{ln}\lambda\,{\mathrm{d}}\lambda\text{,}\label{eq:delr2GCRSgeo}\end{aligned}$$ and $$\Delta^{(2)}_{\mathrm{r,drag}}(\mathbf{x}_A,\mathbf x_B)=-2\xi^i(\mathbf x_B)N_{AB}^i\Delta^{(1)}_{\mathrm r}(\mathbf{x}_A,\mathbf x_B)\text{.}
  \label{eq:delr2GCRSdrag}$$

We see that the post-linear approximation of the refractive time delay function can be split into three components. The first one, namely Eq. , is the second-order correction to the excess path delay . The second component, that is Eq. , is the geometric delay which accounts for the bending of the ray. These two components together with Eq. , constitute the static part of the refractive time delay $$\begin{aligned}
  \Delta_{\mathrm{r,stat}}(\mathbf x_A,\mathbf x_B)&=\Delta_{\mathrm r}^{(1)}(\mathbf x_A,\mathbf x_B)+\Delta_{\mathrm{r,exc}}^{(2)}(\mathbf x_A,\mathbf x_B)\nonumber\\
  &+\Delta_{\mathrm{r,geo}}^{(2)}(\mathbf x_A,\mathbf x_B)\text{,}\label{eq:delrstat}\end{aligned}$$ namely the refractive part of delay that would be measured or modeled in a frame comoving with the media. Instead, the last term in Eq. , namely Eq. , is the delay due to the dragging of light caused by the motion of the dielectric medium.

Interestingly, when one considers that the main component of the motion is due to the steady rotation of the atmosphere in GCRS coordinates, we may emphasis that the light-dragging component shows up about the same order than the geometric delay. Even more interesting is the fact that the light-dragging contribution can be expressed as a geometric factor scaling the static refractive part of the previous order (see Eq. ). This fact is not a specificity of the post-linear approximation but must hold true for higher order terms too. Indeed, it results from the really specific form of the refractive components $\kappa^{\mu\nu}$. As a matter of fact, the components $\kappa^{0i}$ can always be written such as $$\kappa^{0i}=\kappa^{00}\xi^i\text{.}$$ Therefore, because the scalar product $\xi^iN_{AB}^i$ is independent of the path of integration for a steady rotating atmosphere, the integration of $\kappa^{0i}N_{AB}^i$ reduces to $$\xi^i(\mathbf x_B)N_{AB}^i\int_0^1(\kappa^{00})_{\mathbf z_-(\lambda)}{\mathrm{d}}\lambda$$ where the integrated term corresponds to the static part of the refraction.

Solving the line integrals in Eqs. , , and for a realistic index of refraction is not an easy task. Moreover it is beyond the scope of this paper which aims at introducing a recursive method allowing one to determine the integral form of the time transfer functions up to any order in optical spacetime. For this reason, we address the effective resolution of the line integrals to future work.

Hereafter, we derive the range transfer function at the post-linear approximation from Eqs.  by making use of Eqs. , , and $$\begin{aligned}
  \mathcal{R}(\mathbf x_A,\mathbf x_B)&=R_{AB}+\Delta_{\mathrm{g}}(\mathbf{x}_A,\mathbf x_B)\nonumber\\
  &+D(\mathbf x_A,\mathbf x_B)\,\Delta^{(1)}_{\mathrm{r}}(\mathbf{x}_A,\mathbf x_B)\nonumber\\
  &+\Delta^{(2)}_{\mathrm{r,exc}}(\mathbf{x}_A,\mathbf x_B)+\Delta^{(2)}_{\mathrm{r,geo}}(\mathbf{x}_A,\mathbf x_B)\label{eq:RTFpl}\end{aligned}$$ where we have introduced the light-dragging coefficient $D(\mathbf x_A,\mathbf x_B)$ being defined by $$D(\mathbf x_A,\mathbf x_B)=1-2\xi^i(\mathbf x_B)N_{AB}^i\text{.}
  \label{eq:Kdrag}$$ According to the previous discussion, we can rewrite Eq. , within the same accuracy, such as $$\begin{aligned}
  \mathcal{R}(\mathbf x_A,\mathbf x_B)&=R_{AB}+\Delta_{\mathrm{g}}(\mathbf{x}_A,\mathbf x_B)\nonumber\\
  &+D(\mathbf x_A,\mathbf x_B)\,\Delta_{\mathrm{r,stat}}(\mathbf{x}_A,\mathbf x_B)\text{.}\label{eq:RTFplD}\end{aligned}$$ The time transfer function can be directly obtained by making use of Eq. $$\begin{aligned}
  \mathcal{T}(\mathbf x_A,\mathbf x_B)&=\frac{1}{c}\Big[R_{AB}+\Delta_{\mathrm{g}}(\mathbf{x}_A,\mathbf x_B)\nonumber\\
  &+D(\mathbf x_A,\mathbf x_B)\,\Delta_{\mathrm{r,stat}}(\mathbf{x}_A,\mathbf x_B)\Big]\label{eq:TTFpl}\end{aligned}$$ where we recall that $\Delta_{\mathrm{r,stat}}$ is given in Eq. .

Let us emphasis how remarkably simple result is. As a matter of fact, the dragging coefficient is just a geometrical factor scaling the static part of the refractive delay. Indeed, in a covariant theory such like general relativity, let us recall that the light-dragging contribution is naturally handled through the cross components of the Gordon’s metric. Moreover, let us emphasis that the demonstration which ends up with Eq.  is perfectly independent of the refractive profile inside the steady rotating optical medium. In comparison, a derivation of the light-dragging effect using perturbation equations applied to geometrical optics [@2019AeA...624A..41B] requires heavier calculations (the integration must be performed along hyperbolic path) highlighting the advantage of using the covariant formalism developed so far.

From the range or the time transfer functions, we can derive the expression of the frequency transfer within the post-linear approximation as well. After inserting Eq.  within , we deduce

\[eq:qAqB\] $$\begin{aligned}
  q_A&=1-\beta^i_AN_{AB}^i+\beta^i_A\frac{\partial\Delta_{\mathrm{g}}}{\partial x_A^i}+\widehat\beta^i_A\frac{\partial\Delta_{\mathrm{r,stat}}}{\partial x_A^i}\nonumber\\
  &+\beta^i_A\frac{\partial D}{\partial x_A^i}\Delta_{\mathrm{r,stat}}\text{,}\end{aligned}$$ and $$\begin{aligned}
  q_B&=1-\beta^i_BN_{AB}^i-\beta^i_B\frac{\partial\Delta_{\mathrm{g}}}{\partial x_B^i}-\widehat\beta^i_B\frac{\partial\Delta_{\mathrm{r,stat}}}{\partial x_B^i}\nonumber\\
  &-\beta^i_B\frac{\partial D}{\partial x_B^i}\Delta_{\mathrm{r,stat}}\end{aligned}$$

where we have introduced two artificial dragging coordinate velocities defined by $$\widehat\beta^i_{A/B}=D(\mathbf x_A,\mathbf x_B)\,\beta^i_{A/B}\text{.}
  \label{eq:velD}$$

Most of the time, while modeling range and Doppler observables in GCRS coordinates, the dragging coefficient is arbitrarily fixed to $D=1$. In order to investigate the consequences, we close the section by discussing orders of magnitude and variabilities due to the light-dragging contribution in the expressions of the time and the frequency transfers.

Light-dragging magnitude and variability
----------------------------------------

In GCRS coordinates, the velocity of the fluid medium at $\mathbf x_B$ is given by Eq. , that is $$\xi^i(\mathbf x_B)=\frac{\omega_\oplus r_B}{c}\,e^{ijk}S_\oplus^jn_B^k
  \label{eq:3velmedN}$$ where $\mathbf n_B=\mathbf x_B/r_B$. For a ground-based receiver, we have $r_B=R_\oplus$ and the light-dragging coefficient becomes $$D(\mathbf x_A,\mathbf x_B)=1-\frac{2\omega_\oplus R_\oplus}{c}\,(\mathbf S_\oplus\times\mathbf n_B)\cdot\mathbf N_{AB}\text{.}
  \label{eq:KdragGCRS}$$ Thus, the maximum value of $D-1$ is about $$\frac{2\omega_\oplus R_\oplus}{c}\simeq 3.099\times 10^{-6}\text{.}
  \label{eq:Knum}$$

A typical value of the static refractive delay in the zenith direction is approximately $2.5\,\mathrm m$ and can reach $15\,\mathrm m$ for elevation angle of $10\degres$ [@1997JGR...10220489C; @doi1010292006JB004834]. Therefore, the light-dragging contribution to the time transfer is expected to remain lower that $0.05\,\mathrm{mm}$ in GCRS coordinates. However, for experiments whose data are mainly analyzed in the Barycentric Celestial Reference System (BCRS), the velocity of the media possesses an orbital components which is of the order of $30\,\mathrm{km}\cdot\mathrm{s}^{-1}$. Thus, the maximum value of $D-1$ becomes of the order of $2\times10^{-4}$, and the dragging contribution can reach $3\,\mathrm{mm}$ in BCRS coordinates.

Experiments such like Satellite or Lunar Laser Ranging (SLR or LLR) are currently operating at the millimeter and the centimeter level of precision on range measurements . Therefore, the light-dragging effect is just below the threshold of visibility on both experiments. However, as it may be inferred from Eq. , the effect is mainly suppressed in the case of a round-trip light path. In other words, it might play a role only for one-way and three-way observations.

From Eq. , considering a slowly varying refractivity, we can infer that $$\beta_A^i\frac{\partial\Delta_{\mathrm{r,stat}}}{\partial x_A^i}\sim\frac{\ell}{R_{AB}}N_0(\beta_A^iN_{AB}^i)\sim10^{-5}\,(\beta_A^iN_{AB}^i)\text{.}$$ Therefore, for a one-way frequency transfer experiment, the static atmospheric contribution relative to the classical effect $(\beta_A^iN_{AB}^i)$, represents one part in $10^{5}$. Then, the contribution due to the dragging velocity in Eqs. is approximately given by $$\widehat{\beta}_A^i\frac{\partial\Delta_{\mathrm{r,stat}}}{\partial x_A^i}\sim10^{-5}\,(\widehat{\beta}_A^iN_{AB}^i)\text{.}$$ Making use of Eqs. and , one infers that, in GCRS coordinates, the light-dragging contribution relative to the static atmospheric effect, represents one part in $10^6$ and one part in $10^{11}$ relative to the classical effect. On the opposite, if we take a look at orders of magnitude in BCRS coordinates, the light-dragging contribution relative to the the static atmospheric effect, reaches one part in $10^4$ and one part in $10^9$ relative to the classical effect. Therefore, for typical spacecraft’s velocities of $10^{-5}$ and $10^{-4}$ in GCRS and BCRS, respectively, one infers that the effect of the light-dragging contribution produces a fractional frequency change of the order of one part in $10^{16}$ in GCRS coordinates and one part in $10^{13}$ in BCRS coordinates. For one-way radio links, these fractional frequency changes translate into radio signal frequencies at the level of $1\,\mu\mathrm{Hz}$ for $\mathrm{X}/\mathrm{Ka}$-bands and $0.1\,\mu\mathrm{Hz}$ for $\mathrm{S}$-band in GCRS coordinates. In BCRS coordinates, the frequencies of the radio signal due to the dragging of light should arise at $1\mathrm{mHz}$ for $\mathrm{X}/\mathrm{Ka}$-bands and $0.1\,\mathrm{mHz}$ for $\mathrm{S}$-band. The correspondence in term of velocity precision in the Doppler is at the level of $0.01\,\mu\mathrm{m}\cdot\mathrm{s}^{-1}$ and $10\,\mu\mathrm{m}\cdot\mathrm{s}^{-1}$ in GCRS and BCRS coordinates, respectively.

Past and future space missions like Cassini [@2003Natur.425..374B; @Kliore2004; @2007IJMPD..16.2117I], BepiColombo [@2002PhRvD..66h2001M; @2009AcAau..65..666I], or JUICE [@Grasset20131] have reached or will reach the level of the $\mu\mathrm{m}\cdot\mathrm{s}^{-1}$ for the Doppler. Therefore, the light-dragging effect is clearly at the threshold of visibility in Doppler observables and should be modeled in data reduction softwares in a close future.

In order to understand what could be the signature of an unaccounted light-dragging effect, let us now focus on the computation of the time variability of $D(\mathbf x_A,\mathbf x_B)$. For a ground-based instrument, the spatial coordinates expressed in an Earth centered frame are given by $\mathbf x_B=(R_\oplus,\phi_B,\lambda_B)$, where $\phi_B$ is the latitude and $\lambda_B$ the longitude of the instrument on the surface of the Earth. The variable part in Eq. is better understood if we introduce $(a,e,\iota,\varOmega,\omega,f)$ denoting the set of Keplerian elements of the emitter. In GCRS coordinates the direction $\mathbf n_A$ of the emitter is given for instance in Eq. (3.42) of [@2014gravbookP]. Then, the expression of the light-dragging coefficient reads as follows $$\begin{aligned}
  &D-1=\frac{\omega_\oplus R_\oplus}{c}\frac{a(1-e^2)}{R_{AB}}\frac{\cos\phi_B}{(1+e\cos f)}\nonumber\\*
  &\times\Big\{\sin\varOmega\big[I_-\cos(F_++P_+)+I_+\cos(F_-+P_-)\big]\nonumber\\*
  &-\cos\varOmega\big[I_-\sin(F_++P_+)+I_+\sin(F_-+P_-)\big]\Big\}\end{aligned}$$ where we have set

$$I_\pm=1\pm\cos\iota\text{,}$$

and $$F_\pm=f\pm \omega_\oplus t\text{,} \qquad P_\pm=\omega\pm\lambda_B\text{.}$$

Considering a quasi-circular orbit ($e\ll 1$), we have $r_A=a+\mathcal{O}(e)$ and $$f=n (t-t_0)+\mathcal{O}(e)$$ where $t_0$ is the time of perigee passage and where $n$ is the mean motion being given by Kepler’s third law $$n=\sqrt{\frac{Gm_\oplus}{a^3}}\,\text{.}$$ Therefore, the magnitude of $D-1$ oscillates with frequencies $n\pm\omega_\oplus$ around zero and $10^{-4}$ (maximum amplitude of the orbital barycentric velocity) in GCRS and BCRS coordinates, respectively. The peak to peak amplitude is of the order of $10^{-6}$ in both reference systems. In the limit case where $\lim_{a\to\infty}n=0$, the same magnitudes oscillate at diurnal frequency.

Consequently, while modeling the time and frequency transfers using Eqs. and in GCRS or BCRS coordinates, the fact of imposing $D=1$ leads to an unaccounted contribution which may lead to systematic errors for instance in the estimations of the spacecraft velocity (considering Eq. ) or in the receiver coordinates (considering that diurnal signatures mainly concern ground-based stations). This last example could be particularly relevant for ground-based techniques operating within the International Earth Rotation and Reference System Service (IERS) for which an error in the estimation of the station coordinates can result in a bias in the determination of the ITRF.

Conclusion
==========

This paper generalizes the algorithmic approach introduced in [@2008CQGra..25n5020T] by making the time transfer functions formalism applicable in optical spacetime. The main results stand in the theorems \[th:theo4\]–\[th:theo11\] which allow one to determine the integral form of the time transfer functions up to any order. The great benefit of using the time transfer functions formalism rely on the fact that all integrals in theorems \[th:theo6\]–\[th:theo11\] are line integrals taken along the zeroth-order null geodesic path between the emitter and the receiver, independently of the order which is considered.

In optical spacetime, the method requires us to know the order of magnitudes of both the gravitational and the refractive perturbations. Then, one can deduce the integer parameter $s$ from Eq. and use theorems \[th:theo4\]–\[th:theo5\] in order to determine the general expansion of the total time delay functions. The different components are the gravitational, the refractive, and the coupling contributions. Each of them is determined recursively making use of theorems \[th:theo7\]–\[th:theo11\]. We emphasis that these theorems have been derived assuming i) a post-Minkowskian expansion and ii) a general expansion in term of an arbitrary refractivity $N_0$. Both choices are motivated by the quasi-Minkowskian path regime which is assumed throughout the paper.

We have illustrated the method by determining the integral form of the time transfer function up to the post-linear approximation. We have considered the case of a one-way transfer between a low orbit emitter and a receiving station on the Earth’s surface. We have shown that the time and frequency transfers are both impacted by the light-dragging effect due to the motion of the atmosphere, as seen from a frame which is not comoving with the flowing optical media. With respect to other methods [@2019AeA...624A..41B], we have highlighted the great advantage of the covariant formalism developed in this paper which naturally takes into account the effect of the dragging of light. In addition, we have shown that the light-dragging contribution is independent of the refractive profile which is considered. At the end of the day, the dragging component reduces to a geometrical factor which scales the static part of the atmospheric time delay (where the term static refers to the delay which would be measured in a frame comoving with the refractive medium). Concerning the frequency transfer, we have shown that the light-dragging contribution scales the coordinate velocities of both the emitter and the receiver resulting in the introduction of artificial dragging coordinate velocities. Finally, we have discussed the necessity, in a close future, for taking into account the dragging of light in data reduction softwares modeling the time and frequency transfers within inertial frames (e.g. GCRS and BCRS).

Acknowledgments {#acknowledgments .unnumbered}
===============

The author is grateful to the university of Bologna and to Italian Space Agency (ASI) for financial support through the Agreement 2013-056-RO in the context of EAS’s JUICE mission. The author is also thankful to A. Hees, P. Teyssandier, and C. Le Poncin-Lafitte from SYRTE in Observatoire de Paris for interesting discussions and valuable comments about a preliminary version of the manuscript.

General expansion of $\gamma_{\mu\nu}$ {#sec:genexp}
======================================

The covariant components of $\gamma_{\mu\nu}$ are determined from the inverse conditions which lead to the following implicit expression $$\gamma_{\mu\nu}=-g_{\mu\alpha}g_{\beta\nu}\kappa^{\alpha\beta}-g_{\mu\alpha}\kappa^{\alpha\beta}\gamma_{\beta\nu}\text{.}
  \label{eq:gmrec}$$ Usually, assuming that $\gamma_{\mu\nu}=f(n)w^{\mu}w^{\nu}$ with $f(n)$ being a sought function of the index of refraction and using Eq. , we infer Eq. . However, the situation slightly changes if we expand the contravariant components $\kappa^{\mu\nu}$ as it is done in Eq. .

At the same time, we have assumed that the physical spacetime metric, which is given in Eq. , satisfies a post-Minkowskian expansion (see Eq. ). Thus, considering that the refractive components are the dominant order according to Eq. for $s\in\mathbb{N}_{>0}$, we deduce that the covariant components $\gamma_{\mu\nu}$ satisfy the following expansion $$\gamma_{\mu\nu}(x,N_0,G)=\sum_{l=1}^{\infty}\gamma_{\mu\nu}^{(l)}(x)
  \label{eq:gamPM}$$ where the quantities $\gamma_{\mu\nu}^{(l)}$ can be recursively determined from Eq. , that is

\[eq:kapPMord\] $$\begin{aligned}
  \gamma_{\mu\nu}^{(1)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(1)}\text{,}\label{eq:kapPMord1}\\
  \gamma_{\mu\nu}^{(q)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(q)}-\eta_{\mu\alpha}\sum_{m=1}^{q-1}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(q-m)}\label{eq:kapPMord2}\end{aligned}$$ for $2\leqslant q\leqslant s$, and $$\begin{aligned}
  \gamma_{\mu\nu}^{(s+1)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(s+1)}-2\eta_{\mu\alpha}\kappa^{\alpha\beta}_{(1)}h_{\beta\nu}^{(1)}\nonumber\\
  &-\eta_{\mu\alpha}\sum_{m=1}^{s}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(s-m+1)}\text{,}\label{eq:kapPMord3}\\
  \gamma_{\mu\nu}^{(s+q)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(s+q)}-2\eta_{\mu\alpha}\kappa^{\alpha\beta}_{(q)}h_{\beta\nu}^{(1)}\nonumber\\
  &-\eta_{\mu\alpha}\sum_{m=1}^{s+q-1}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(s-m+q)}\nonumber\\
  &-h_{\mu\alpha}^{(1)}\sum_{m=1}^{q-1}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(q-m)}\label{eq:kapPMord4}\end{aligned}$$ for $2\leqslant q\leqslant s$, and $$\begin{aligned}
  \gamma_{\mu\nu}^{(2s+1)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(2s+1)}-2\eta_{\mu\alpha}\kappa^{\alpha\beta}_{(s+1)}h_{\beta\nu}^{(1)}\nonumber\\
  &-2\eta_{\mu\alpha}\kappa^{\alpha\beta}_{(1)}h_{\beta\nu}^{(2)}-h_{\mu\alpha}^{(1)}\kappa^{\alpha\beta}_{(1)}h^{(1)}_{\beta\nu}\nonumber\\
  &-\eta_{\mu\alpha}\sum_{m=1}^{2s}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(2s-m+1)}\nonumber\\
  &-h_{\mu\alpha}^{(1)}\sum_{m=1}^{s}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(s-m+1)}\text{,}\label{eq:kapPMord7}\end{aligned}$$ and $$\begin{aligned}
  \gamma_{\mu\nu}^{(2s+q)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(2s+q)}-2\eta_{\mu\alpha}\kappa^{\alpha\beta}_{(q)}h_{\beta\nu}^{(2)}\nonumber\\
  &-2\eta_{\mu\alpha}\kappa^{\alpha\beta}_{(s+q)}h_{\beta\nu}^{(1)}-h_{\mu\alpha}^{(1)}\kappa^{\alpha\beta}_{(q)}h^{(1)}_{\beta\nu}\nonumber\\
  &-\eta_{\mu\alpha}\sum_{m=1}^{2s+q-1}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(2s+q-m)}\nonumber\\
  &-h_{\mu\alpha}^{(1)}\sum_{m=1}^{s+q-1}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(s+q-m)}\nonumber\\
  &-h_{\mu\alpha}^{(2)}\sum_{m=1}^{q-1}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(q-m)}\label{eq:kapPMord8}\end{aligned}$$ for $2\leqslant q\leqslant s$, and $$\begin{aligned}
  \gamma_{\mu\nu}^{(ps+1)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(ps+1)}-2\eta_{\mu\alpha}\sum_{m=0}^{p-1}\kappa^{\alpha\beta}_{(ms+1)}h_{\beta\nu}^{(p-m)}\nonumber\\
  &-\sum_{m=1}^{p-1}h_{\mu\alpha}^{(p-m)}\sum_{n=0}^{m-1}\kappa^{\alpha\beta}_{(ns+1)}h^{(m-n)}_{\beta\nu}\nonumber\\
  &-\eta_{\mu\alpha}\sum_{m=1}^{ps}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(ps+1-m)}\nonumber\\
  &-\sum_{m=1}^{p-1}h_{\mu\alpha}^{(p-m)}\sum_{n=1}^{ms}\kappa^{\alpha\beta}_{(n)}\gamma_{\beta\nu}^{(ms+1-n)}\label{eq:kapPMord9}\end{aligned}$$ for $p\geqslant 3$, and finally $$\begin{aligned}
  \gamma_{\mu\nu}^{(ps+q)}&=-\eta_{\mu\alpha}\eta_{\beta\nu}\kappa^{\alpha\beta}_{(ps+q)}-2\eta_{\mu\alpha}\sum_{m=0}^{p-1}\kappa^{\alpha\beta}_{(ms+q)}h_{\beta\nu}^{(p-m)}\nonumber\\
  &-\sum_{m=1}^{p-1}h_{\mu\alpha}^{(p-m)}\sum_{n=0}^{m-1}\kappa^{\alpha\beta}_{(ns+q)}h^{(m-n)}_{\beta\nu}\nonumber\\
  &-\eta_{\mu\alpha}\sum_{m=1}^{ps+q-1}\kappa^{\alpha\beta}_{(m)}\gamma_{\beta\nu}^{(ps+q-m)}\nonumber\\
  &-\sum_{m=0}^{p-1}h_{\mu\alpha}^{(p-m)}\sum_{n=1}^{ms+q-1}\kappa^{\alpha\beta}_{(n)}\gamma_{\beta\nu}^{(ms+q-n)}\label{eq:kapPMord10}\end{aligned}$$

for $p\geqslant 3$ and $2\leqslant q\leqslant s$, where $p$ and $q$ are determined from $l$ using Eqs. .

[^1]: This can be easily shown from Eq. using the inverse conditions $\bar g_{\mu\sigma}\bar g^{\sigma\nu}=\delta_\mu^\nu$.

[^2]: The reception and the emission time delay functions in [@2008CQGra..25n5020T] are denoted $\Delta_r/c$ and $\Delta_e/c$, respectively. Instead, in this work, we use $\Delta/c$ and $\Xi/c$ in order to keep incoming indices notations as clear as possible.

[^3]: We recall that the time component of $x_A$ is constrained by the fact that the past null cone at $x_B$ intersects the world line $\mathbf x=\mathbf x_A$ at only one point of coordinates $\bm x_A=(ct_A,\mathbf{x}_A)$. This is true as long as the null geodesic is quasi-Minkowskian.

[^4]: Hereafter, $\mu$ is an affine parameter and does not represent the permeability of the dielectric medium anymore.

[^5]: Results assuming that the gravitational perturbation is the leading term can easily be derived by applying the exact same following approach. However, in that case, it must be noted that $s$ should be introduced switching the numerator and the denominator in the right-hand side of Eq. .
