Accurate light-time correction due to a gravitating mass

# Accurate light-time correction due to a gravitating mass

Neil Ashby
Department of Physics
University of Colorado, Boulder, Co. (USA)

Bruno Bertotti
Dipartimento di Fisica Nucleare e Teorica
Università di Pavia (Italy)
###### Abstract

This technical paper of mathematical physics arose as an aftermath of Cassini’s 2002 experiment [6], in which the PPN parameter was measured with an accuracy and found consistent with the prediction of general relativity. The Orbit Determination Program (ODP) of NASA’s Jet Propulsion Laboratory, which was used in the data analysis, is based on an expression (8) for the gravitational delay which differs from the standard formula (2); this difference is of second order in powers of – the gravitational radius of the Sun – but in Cassini’s case it was much larger than the expected order of magnitude , where is the distance of closest approach of the ray. Since the ODP does not take into account any other second-order terms, it is necessary, also in view of future more accurate experiments, to revisit the whole problem, to systematically evaluate higher order corrections and to determine which terms, and why, are larger than the expected value. We note that light propagation in a static spacetime is equivalent to a problem in ordinary geometrical optics; Fermat’s action functional at its minimum is just the light-time between the two end points A and B. A new and powerful formulation is thus obtained. This method is closely connected with the much more general approach of [18], which is based on Synge’s world function. Asymptotic power series are necessary to provide a safe and automatic way of selecting which terms to keep at each order. Higher order approximations to the required quantities, in particular the delay and the deflection, are easily obtained. We also show that in a close superior conjunction, when is much smaller than the distances of A and B from the Sun, of order , say, the second-order correction has an enhanced part of order , which corresponds just to the second-order terms introduced in the ODP. Gravitational deflection of the image of a far away source when observed from a finite distance from the mass is obtained up to .

## 1 Introduction

In the framework of metric theories of gravity and the PPN formalism, the main violations of general relativity – those linear in the masses – are described by a single dimensionless parameter . The question, at what level and how general relativity is violated, in particular how much differs from unity, Einstein’s value, is still moot. No definite and consistent prediction about it are available, except for the inequality , which must be fulfilled in a scalar-tensor theory, in particular those arising as the low-energy limit of certain string theories. To date, the best measurement of has been obtained with Cassini’s experiment, which has provided the fit (at 1-)

 γ−1=(2.1±2.3)×10−5. (1)

Einstein’s prediction is still acceptable, but more accurate experiment are needed and planned.

While controls also other relativistic effects, in particular those related to gravito-magnetism, it mainly affects electromagnetic propagation. The differential displacement of the stellar images near the Sun historically was the first experimental effect to be investigated and is now of great importance in accurate astrometry. The bending of a light ray also increases the light-time between two points, an important effect usually named after its discoverer I. I. Shapiro [27]. Several experiments to measure this delay have been successfully carried out, using wide-band microwave signals passing near the Sun and transponded back, either passively by planets, or actively, by space probes (see [31], [24]).

Cassini’s 2002 experiment has implemented a third way to measure [4], in which coherent microwave trains sent from the ground station to the spacecraft (at that time about 7 AU far away) were transponded back continuously. The use of high-frequency carriers (in K band, 34 and 32 GHz) and the combination with standard X-band carriers (about 8 GHz) allowed successful elimination of the main hindrance, dispersive effects due to the solar corona traversed by the beam. The tracking was carried out around the 2002 superior conjunction; the minimum value of the impact parameter of the beam was , but in effect only 18 passages have been used, with a minimum impact parameter of . The two-way total amount of phase between the time of emission and the time of arrival has been continuously measured in each passage. In effect, however, NASA’s Deep Space Network provides the phase count in a given integration time . Mathematically, in the limit this would give the received frequency, in which Doppler effects and gravitational frequency shift are mixed up (Sec. 4). Cassini’s observable, therefore, can also be assessed in terms of the predicted change in frequency, as in [4]; but in practice, taking small would introduce unacceptable high-frequency noise. The change in light-time in a given integration time is the correct, theoretically available observable.

In the standard formulation for a superior conjunction, and taking the Sun at rest, the (one-way) light-time from an event A to an event B is:

 tB−tA=rAB+Δt=rAB+(1+γ)mlnrA+rB+rABrA+rB−rAB, (2)

where km is the gravitational radius of the Sun, are, in Euclidian geometry (See Fig. 1 left), the distances of A and B from the Sun and their distance. The velocity of light is unity. , the increase of the light-time over , is the gravitational delay.

In a close superior conjunction A and B are on the opposite sides of the mass and the Euclidian distance of the straight line AB from the mass fulfils , say. In this approximation eq. (2) reduces to

 tB−tA=rAB+Δt=rAB+(1+γ)mln(4rArBb20), (3)

with a logarithmic enhancement over the formal order of magnitude .111As stated in the supplementary material, in eq. (2) of [6] the two terms in the right-hand side should obviously be multiplied by a factor 2. This error, of course, had no consequence on the computer fit. Taking the logarithm equal to 10, this provides an estimate of the timing accuracy in terms of the error in :

 σΔt=1.43σγ×106cm, (4)

corresponding, in Cassini’s case, to 30 cm. (3) embodies also the one-way frequency change induced by gravity between A and B. Their motion makes (and the distances) change with time, so that, for a one-way experiment,

 Δνν=dΔtdt=−2(γ+1)mb0db0dt. (5)

The basic geometric setup is straightforward: a point mass at rest at the origin in an asymptotically flat space generates a line element with rotational symmetry. An invariant Killing time is defined; events on each constant surface are ‘simultaneous’ and the metric components are constant. The proper time of a static observer differs from by the red-shift factor . A null geodesic runs from the event A (with radial coordinate and time ) to the event B (with radial coordinate and time ); it stays on a plane, taken here as the equatorial plane . The (invariant) longitude difference completes the setup. In the PPN formalism and isotropic coordinates the metric reads:

 ds2 = A(r)dt2−B(r)dℓ2= (6) = (1−2mr+2βm2r2+…)dt2−(1+γ2mr+3ϵ2m2r2+…)dℓ2,

where

 dℓ2=dr2+r2(dθ2+sin2θdϕ2)=dr2+r2dΩ2

is the Euclidian line element. The parameters , and are equal to 1 in general relativity; while and are accurately known, currently no information is available about .

In our case the best mathematical tool to deal with electromagnetic propagation is not null geodesics, but the theory of eikonal. It is known (e.g., [19]) that in this problem Fermat’s Principle holds, corresponding to the refractive index

 N(r)=√B(r)A(r); (7)

we develop ab initio the eikonal and solve for it by separation of variables (Sec. 4). The radial part provides Fermat’s action as a radial integral containing and the impact parameter ; when computed at the true value , such action is just the required light-time. The solution can be obtained recursively, using appropriate expansions in powers of : the expansion for begins with , the distance of the straight line AB from the origin. In this way the variational nature of the problem brings about a great conceptual and algebraic simplification. At the linear approximation in one would expect that the light-time contains , the correction in the impact parameter linear in the mass; as one can see from (2), this is not the case. This property is generally true: the correction to the light-time does not contain (Sec. 6).

Cassini’s and many other space experiments have been analyzed using NASA’s Orbit Determination Program (ODP), developed by NASA at Jet Propulsion Laboratory in the 60’s and steadily improved since; a new version called MONTE is under development. The ODP, whose theoretical formulation is due to T. D. Moyer [21], integrates the equations of motion of the relevant bodies and provides their trajectories in the ephemeris time. This task is carried out in a reference system – called BCRS (Barycentric Coordinate Reference System) – in which the centre of gravity of the solar system is at rest and the Sun moves around with a velocity . As discussed in [2], the light-time in this this frame differs from the rest frame of the Sun essentially due to Lorentz time dilatation; being of order , this difference is quite below the sensitivity of Cassini’s experiment. We do not discuss this point any more; is just Killing time.

The ODP uses a fictitious Euclidian space , which corresponds to the isotropic coordinates of (6). This space is just a computational convenience and should not be considered as a physical background in which gravity acts. For example, replacing , the Euclidian distance from the origin, with , where is an arbitrary constant, is fully legitimate in a covariant theory, but it destroys the conformal flatness of space, introduces a gravitational potential and adds a second-order term to the delay . Strictly speaking, the word ‘delay’ is inappropriate: we just have a light-time and there is nothing with respect to which a delay can be reckoned. The object of the measurement is the time change of the delay. The arbitrariness of the radial coordinate also affects gravitational bending: its second-order approximation up to , depends on which radial coordinate is used (see [13], [10] and [25]) [7]).

It should also be noted that the spacetime coordinates of the end events are not directly provided in the experimental setup and depend on the gravitational delay , the very quantity one sets out to measure. The trajectories and are given by the numerical code; the starting time is just a label of the ray, but the arrival time is greater than . The way out is to take for the end point

 rB(rB)=rB(tA+rAB)+ΔtuB(tA+rAB),

where For a typical velocity the correction is of order cm, and the a priori accuracy in is sufficient.

Since for electromagnetic propagation and in (6) are almost equal, (2) is the correct approximation to the delay to ; one would expect this to be the first term in an expansion in powers of , so that the next term should be

 ≈mmb0=mmR⊙R⊙b0=0.3R⊙b0cm,

quite below Cassini’s sensitivity. The present paper arose because the ODP (eq. (8-54) of [21]), in fact does not use (2), but, in our notation,

 tB−tA=rAB+Δt=rAB+(1+γ)mln(rA+rB+rAB+(1+γ)mrA+rB−rAB+(1+γ)m). (8)

We have not been able to fully reconstruct Moyer’s derivation of this expression. It introduces non linear corrections arising from non linear effects of linear metric terms, but no quadratic metric terms. However, the difference between the two expressions of the delay is much larger than the estimate above; this arises because in Cassini’s case, in (2) the denominator is much smaller than the numerator . Indeed,

 Δt−(Δt)ODP=−2(1+γ)2m2b20rArBrA+rB=−(1+γ)2m2Rb20, (9)

where we have introduced the harmonic mean of the distances

 2R=1rA+1rB=rA+rBrArB. (10)

If, as in Cassini’s experiment, , the correction is about

 1600mmR⊙(R⊙b0)2=500(R⊙b0)2cm.

Even at this correction is somewhat below the sensitivity (4) and it should not have affected the result. However, it cannot be excluded that neglected non linear terms relevant for Cassini’s experiment affect the fit (1). One could say, (8) is mendacious; a full clarification of the problem is needed.

Empirically dropping or keeping ‘small’ terms may lead to inconsistencies and does not work; the rigourous method of asymptotic perturbation theory (see, e. g., [11], [16]) must be used. We briefly sketch it now at a practical level. One begins with a wise choice of a dimensionless ‘smallness’ parameter, and expands every function in the corresponding power series. Our main choice will be , but convenience may suggest using other lengths, like in . An asymptotic series

 G=∑sGs(mb0)s

is a formal object assigned just by the sequence of its coefficients ; arithmetics and calculus follows the obvious rules for sum, multiplication and differentiation. Equality between two asymptotic series just means that the coefficients of the same order are equal. The value of as a function of does not play any role, and even the convergence of the series is irrelevant; what matters is only the truncated value at any order

 G(k)=k∑s=0(mb0)sGs+O(mb0)k+1. (11)

The parameter should not be understood as a fixed number, but as a variable which tends to zero. The symbol means order of infinitesimal; it states how fast the remainder tends to zero as the parameter diminishes. An asymptotic series can be constructed from an ordinary arbitrary function ; but a whole class of functions give rise to the same series; for example, if is the sequence generated by , the same sequence is also generated by

 [1+Pexp(−Qb0/m)]G(m)(Q>0).

In this way any recursive iteration then proceeds automatically and safely, even in the most complex situations.

In our case light-time will be provided as an asymptotic power series

 tB−tA=rAB+m∑s=1Δs(rAb0,rBb0)(mb0)s−1, (12)

with dimensionless coefficients . provides the lowest, standard approximation to (see (2). In principle, asymptotic analysis does not provide a numerical estimate of the remainder in a given situation; this is a physical, not a mathematical question. But when the problem, properly formulated, does not contain small dimensionless quantities other than the smallness parameter itself, one can expect the mathematical operations leading to the result to maintain the order of magnitude and to lead to expansions whose coefficients are numerically of the same order. This is the case of deflection, the angle between the asymptotes of the ray. There is only one length in the problem, the distance of the point of closest approach, or, equivalently, the impact parameter (see Fig. 4); hence in the expansion

 δ=∑sδs(mh)s (13)

the coefficients are dimensionless numbers, solely determined by the PPN parameters and, must be of order unity (see Sec. 9). But in the delay problem the coefficients depend on the geometrical configuration. They are of order unity in the generic (but scarcely interesting) case in which and are of the same order; but in a close superior conjunction – of crucial relevance in experimental gravitation – when , besides , there is another smallness parameter, namely, , and there is no reason to exclude that the increase with beyond the expected order of magnitude unity. This we call enhancement. We already saw in (3) that is enhanced, albeit only logarithmically; the ODP correction (9), formally of second order, is enhanced by . This could place serious limitations on the method and even invalidate the iteration itself. This would occur, for instance, when ; if , this corresponds to AU. The enhancement, which has never been discussed in the literature, has been fully understood and tamed in the present paper (Sec. 8). We have found, indeed, that the second-order terms embodied in the ODP expression (8) which was used in Cassini’s experiment are just the enhanced second-order terms; Cassini’s result (1) is still safe.

The problem can be reduced to one of ordinary optics; due to its variational nature, the eikonal function can be easily solved in an expansion in powers of . The second-order expression of the light-time for a static spacetime has been obtained; extension to third order is also easy. This approach should be compared with the much more general work of [18], who consider Synge’s world function in a generic spacetime for a generic geodesic (not necessarily null) between two events A and B. On the basis of Hamiltonian theory, they develop a method to solve for in a formal power series with respect to the gravitational constant and compute it up to the second order. In the null case the world function vanishes on the solution and becomes the eikonal function. Out method, limited of course to the spherically symmetric case, exploits directly the variational nature of the problem and leads to the second-order expression of the light-time, which agrees with the expression of [18]; extension to third order is also easy.

For a realistic observation of a distant source from a point B at a finite distance , (13) must be generalized to an expansion of the type

 δB=∑sδBs(rBh)(mh)s, (14)

where is the impact parameter. The linear term has been evaluated in [19], \$ 40.3; the quadratic correction will be obtained in Sec.9 .

## 2 Hyperbolic Newtonian dynamics

Newtonian dynamics of a test particle attracted by a point mass , an exactly soluble problem, illustrates these issues. We consider a motion in the equatorial plane , with radial coordinate and azimuthal longitude . The Lagrangian function

 LNew=12[(drdt)2+r2(dϕdt)2]+GMr (15)

keeps the total energy constant; , the ultimate speed of the particle at a large distance, plays a role analogous to the speed of light and will be taken equal to unity. Then

 (drdt)2+r2(dϕdt)2−2mr=1, (16)

where is the gravitational radius. is an ignorable coordinate, so that the angular momentum

 ∂LNew∂(dϕ/dt)=r2dϕdt=h (17)

is constant. Since the velocity at infinity is 1, is also the impact parameter. Eliminating we get:

 rdϕdr=±h√r(r+2m)−h2; (18)

hence

 h=√b(b+2m) (19)

determines , the distance of closest approach where . The sign depends upon whether the ray is ingoing or outgoing. Integrating we get the true anomaly

 f=arccos(b2+2mb−mrr(b+m)). (20)

Alternatively, the motion can be expressed in terms of the semi-major axis and the hyperbolic eccentricity :

 r=a(e2−1)1+ecosf. (21)

The acute angle between the asymptotes is given by

 sinδ=sin(2arccos(−1e))=2mb+m√1−m2(b+m)2 (22)

This angle has a regular expansion in powers of , with no enhancement.

Consider, however, the hyperbola determined by two points and on the opposite sides of the vertex (right side in Fig. 1). As in space navigation – in particular in the ODP – the end points are provided in terms of the initial and final position vectors, or equivalently, in terms of the initial and final distances and and the elongation angle ; the “unperturbed distance of closest approach” may then be calculated from elementary geometry:

 b0=rArBsinΦAB√r2A+r2B−2rArBcosΦAB. (23)

Choosing the angles , and positive, we can express in terms of with the condition

 ΦAB = αA+αB=arccosb0rA+arccosb0rB= (24) = arccos(b(b+2m)rA(b+m)−mb+m)+ + arccos(b(b+2m)rB(b+m)−mb+m).

The symmetric case is sufficient to exhibit the problem. The condition reads:

 bRb+2mb+m−mb+m−b0R=0, (25)

or

 b2+(2m−1)b−m(R+1)=0, (26)

with the solution

 2bb0=b0−mb0+√1+4mb0R+mb0. (27)

Expansion in powers of gives

 bb0=1+(Rb0−1)mb0−((Rb0)2−1)(mb0)2+O(mRb0)3. (28)

The enhancement is clear: when the truncation error at order is , with a coefficient of order unity, as naïvely expected; but when – as in a close superior conjunction – , the error is larger, . Formally this requires introducing another smallness parameter and expanding every coefficient of the primary -expansion in descending powers of . Of course, the condition

 mRb20≪1, (29)

must be fulfilled, lest the whole procedure breaks down. One could say, anchoring the trajectory at far away end points has a lever effect, so that an increase in the mass produces a large increase in closest approach.

The quantity (29) gives, in order of magnitude, the ratio between the deflection and the angle which separates the central mass and a distant star, as seen from a distance . Hence the limiting constraint above implies that the geometry of astronomical deflection is the same as in the classical case (see Fig. 4): sources in the sky near the Sun are displaced outward by an amount inversely proportional to the angular distance. The transition through the milestone marks the passage to the gravitational lensing regime, in which the image can appear on both sides.

In Sec. 8 the light-time enhancement is dealt with in the general case and it is shown that the dimensionless coefficients in (12) are

The metric of a spherical body at rest has the general form

 ds2=A(r)dt2−B(r)dr2−C(r)r2dΩ2, (30)

where are of power series the form:

 A(r)=∑sAs(mr)s. (31)

It is asymptotically flat, so that . The radial coordinate is otherwise arbitrary; this is the gauge freedom at our disposal. For consistency, however, any change must become an identity at infinity and have a similar expansion:

 g(r)=r+g1m+g2m2r+…; (32)

the coefficients are not gauge invariant. Two gauges are common. In the isotropic form – the canonical choice in space physics – , so that

 ds2=A(r)dt2−B(r)(dr2+r2dΩ2)=A(r)dt2−B(r)dℓ2; (33)

the space part of the metric is conformally flat. We define

 N(r)=√B(r)A(r)=∑sNs(mr)s=1+N1mr+N2(mr)2+O(mr)3. (34)

In the PPN scheme (e. g., [31])

 N1=γ+1,N2=6−4β+3ϵ+4γ−2γ24. (35)

In ‘Schwarzschild’ gauge and

 ds2=¯A(¯r)dt2−¯B(¯r)d¯r2−¯r2dΩ2;

the area of a sphere of radius is just the Euclidian expression , which defines in an invariant way. In the original Schwarzschild solution . To get the isotropic form one requires

 g2(r)=¯B(g(r))(dgdr)2r2; (36)

to first order

 ¯r=r+γm+…. (37)

In the present paper a third radial coordinate

 ρ=rN(r)=r√B(r)A(r)=r+mN1+m2N2r+… (38)

plays an important role. It is a monotonic function of and ensures . In the linear approximation it was introduced by Moyer in [21] (eq. (8-23)), and boils down to just adding to a constant term, equal to 2.95 km for the Sun.

## 4 Geometrical optics

It is convenient to reduce the problem to geometrical optics using the eikonal function . In a generic spacetime fulfils the eikonal equation

 gμν∂μS∂νS=0; (39)

its characteristics are the null rays (see, e. g., [1]). is the phase of the electromagnetic wave. Let be the trajectories of the end points, given as functions of their proper times ; let

 vμA=drμdsA,vμB=drμdsB

be the corresponding four-velocities. Clocks associated with them measure the proper frequencies

 ωA=−vμA∂μS=dSdsA,ωB=−vμB∂μS=dSdsB. (40)

In the simple case in which the end points are far away from the source, where the metric corrections can be neglected, the contribution to the frequency difference corresponds to the ordinary Doppler effect, and can be evaluated with a slow motion expansion; the change in between A and B is determined by the accumulated gravitational effect along the ray and mainly come from the region near the mass.

 gμν∂μS∂νS=0=N2(r)(∂tS)2−∇S⋅∇S, (41)

where is the Euclidian gradient operator. We are really interested only in the spherically symmetric case, but the reasoning of this Section holds also for an arbitrary .

is the phase; propagation occurs keeping it constant. Separating space and time variables with

 S=St(t)+¯¯¯¯¯Sr(r),

leads to the class of solutions

 S=ω0(¯¯¯¯¯S(r)−t), (42)

where is the spatial part of the phase and is a constant frequency. has the dimension of time and satisfies

 ∇¯¯¯¯¯S⋅∇¯¯¯¯¯S=N2(r). (43)

If a clock is at rest relative to the mass, , and the measured proper frequency includes the appropriate gravitational shift away from the asymptotic value . This is enough to reduce the problem to geometrical optics (see, e. g., [8], Ch. III). A ray , as function of the Euclidian arc length , is orthogonal to the eikonal surfaces and fulfils

 ddℓ(N(r)drdℓ)=∇N(r). (44)

The index of refraction is the rate of increase of the spatial phase along the ray:

 d¯¯¯¯¯Sdℓ=N(r).

Consider now Fermat’s action functional

where the trajectory, any path joining the end points, is expressed in terms of a generic parameter :

 r(λA)=rA,r(λB)=rB. (46)

Since the action is, in fact, independent of the choice of , no generality is lost if , the Euclidean line element. The Euler-Lagrange equation for the action (45) reduces to (44). The actual elapsed time

is just the value of computed at a local minimum – the actual ray (Fermat’s Principle). One should keep in mind the distinction between the action functional, with its argument in square brackets, and the action computed at the extremum, an ordinary function of the end points denoted with . In , but not in , it is allowed to replace the generic independent variable with a more convenient one related to the solution, like . For simplicity, the different functions denoted by the symbol are distinguished by their arguments; below, the quantity will be introduced to denote the action corresponding to a ray anchored at and , but with arbitrary (or ).

## 5 The solution

The eikonal function provides a deep simplification in the evaluation of the light-time. Having already separated out the time, the three-dimensional eikonal equation (43) in spherical symmetry and in the equatorial plane can be solved by separating out the longitude : setting

 ¯¯¯¯¯S(r,ϕ)=¯¯¯¯¯Sr(r)+¯¯¯¯¯Sϕ(ϕ).

It satisfies 222For a function of a single variable a prime indicates the derivative.

 r2(¯¯¯¯¯S′r)2+(¯¯¯¯¯S′ϕ)2=r2N2(r),

so that is a constant. Setting , the eikonal equation reduces to

 (¯¯¯¯¯S′r)2=1r2(r2N2(r)−h2),

with the primitive

 ¯¯¯¯¯S(r)=±∫rdrr√r2N2(r)−h2.

The and the signs correspond, respectively, to an outgoing and an incoming photon. The radial coordinate of closest approach , where , is the solution of

 bN(b)=h; (48)

since , is a real function. In Sec. 9 it will be shown that , just like in the Newtonian case, is the impact parameter (Fig. 4). The total phase is, therefore,

 S=ω0(hϕ±∫rdrr√r2N2(r)−h2−t). (49)

A wavefront propagates keeping constant, so that the time along the ray is

 t=±∫rdrr√r2N2(r)−h2+hϕ. (50)

In the usual case (see Fig. 1), in which the angle is obtuse, the ray has two branches, both taken with the positive sign: an incoming one from to and an outgoing one from to . In the acute case is never reached and we have just an outgoing ray from to . In both cases, in going from A to B the longitude increases by The quantity

 S(h)=∫rBbdrr√r2N2(r)−h2±∫rAbdrr√r2N2(r)−h2+hΦAB (51)

gives the phase change, hence the light-time, between the end points, but the quantity is still arbitrary. The upper (lower) sign corresponds to the case in which the angle is obtuse (acute); in the latter case the two integrals combine in a single one from to , and disappears as a lower limit. (51) is what Fermat’s action functional becomes when its variability is restricted to and the longitude constraint is not imposed; it shall be called reduced action. At the true value it satisfies

 S′(htrue)=0, (52)

keeping the end points fixed.

The present work aims at providing the theoretical foundation for the time delay in all configurations; the sign freedom allows dealing with both cases at the same time, but applications will be mainly given for a conjunction, with the . The origin of longitudes is arbitrary. This general approach is relevant, for ex ample, for a spacecraft on an almost parabolic orbit, as in the Solar Probe concept; with a perihelion as low as , it can have a strong enhancement of the light-time even in the acute configuration.

In the derivative there are no contributions from the lower limits; then (52) provides as an implicit function of the total total elongation :

 ΦAB+∫rBbdrr−h√(rN(r))2−h2±∫rAbdrr−h√(rN(r))2−h2=0. (53)

Hence (51) reads333In a slightly inconsistent notation, we often use to denote both an independent and variable quantity, and the fixed value determined by the elongation. The context should be sufficient to clear the ambiguity.

 S(h)=∫rBbdrN(r)rN(r)√(rN(r))2−h2±∫rAbdrN(r)rN(r)√(rN(r))2−h2. (54)

i

Both integrals are convergent (and in the acute case the singularity at is not even reached). (51) suggests the introduction of the function

 G(r,h)=∫rbdrr√(rN(r))2−h2, (55)

in terms of which

 S(h)=G(rB,h)±G(rA,h)+hΦAB. (56)

(53) reads444The suffix indicates partial derivative.

 Gh(rB,h)±Gh(rA,h)+ΦAB=0. (57)

While in (51) is an independent parameter, in (54) it is fixed by (53).

This expression for can also be derived directly from Fermat’s Principle, thus providing its significance. Fermat’s action (47), expressed as a function of , has the Lagrange functional

 LF[ϕ(r)]=N(r)√1+r2(dϕ/dr)2, (58)

with the (positive) constant of the motion

 ∂LF∂(dϕ/dr)=±r2N(r)√1+r2(dϕ/dr)2dϕdr=h. (59)

The upper (lower) holds for the outgoing (incoming) branch. Integrating

 rdϕdr=±h√r2N2(r)−h2=±h√ρ2−h2, (60)

(53) is recovered. Comparison with the Newtonian case (18) shows that the latter corresponds to the exact index of refraction

 NNew(r)=√1+2mr, (61)

corresponding, as expected, to , and , etc.

## 6 A variational argument

At this point one could proceed as follows: using power series, solve (53) for in terms of , a known quantity. The value of , inserted into (54), provides the required light-time. The stationary character of the action (52), however, brings about a deep and important simplification. This is already tacitly applied in the usual derivation of the gravitational delay (2). To first order, the integral of in (33) reads

the second integral can be carried out along the straight path from A to B, leading to the characteristic logarithmic term. In principle, however, the first integral should take into account the (first order) deflection; we should understand as the Euclidian length of the bent arc between A and B. But the length of the straight segment AB is a minimum in the set of all curves joining A and B, so that vanishes to .555A didactical remark is in order here. This minimum property, crucial to the argument, is often omitted in the usual derivation. See, e. g., [19] p. 1107, [9] p. 125; in equation (17.59) of [5], p. 581 the minimum is not mentioned and a factor 4 is missing in the argument of the logarithm. Ray bending is irrelevant here.

In order to exploit the variational nature of the problem it is convenient to apply power expansions before imposing the extremum condition (52). We just need the value of the reduced action (51) at the value which fulfils (52), namely, Setting and expanding, the solution to second order is obtained iteratively:

 S′0(h0) = 0, (62) h1S′′0(h0) + S′1(h0)=0, (63) h2S′′0(h0) + h212S′′′0(h0)+h1S′′1(h0)+S′2(h0)=0. (64)

In the expression

 S(h) = S0(h0)+m(h1S′0(h0)+S1(h0))+ (65) + m2(h2S′0(h0)+h212S′′0(h0)+h1S′1(h0)+S2(h0)),

the effect of the extremum property is clear: since , the first order term does not contain , and the second order term does not contain ; in general, the term in of order does not depend on . This important result is reflected in the general approach of [18]. Referring to the equation numbering of that paper, their world function fulfills the Hamilton-Jacobi equation (30). In the null case , (30) becomes the eikonal equation. Their Theorem 2 proves that the -order can be expressed in terms of integrals along the lowest order Minkowskian path. In our case this variational Lemma clarifies the matter and produces considerable simplifications. Using (63), the light-time to second order reads:

 S(h)=S0(h0)+mS1(h0)+m2(h212S′′0(h0)+h1S′1(h0)+S2(h0)), (66)

where is given by (63). The delay coefficients (12) read

 Δ1=S1(h0),Δ2=h212S′′0(h0)+h1S′1(h0)+S2(h0). (67)

The second-order correction in the impact parameter , given by (64), is needed only at third and higher orders. For the record, note the third-order contribution to the light-time:

 Δ3=h1h2S′′0+h316S′′′0+h2S′1+h212S′′1+h1S′2+S3, (68)

where, for simplicity, the arguments have been understood, and is provided by (64).

## 7 Power series

We now proceed to apply this simple and general Lemma to the light-time. To lowest order, in (55) we use and , so that

 S0(h0) = √r2B−h20±√r2A−h20+ (69) − h0(arccosh0rB±arccosh0rA−ΦAB).

The condition determines with the trigonometric relation (see Fig. 1)

 ΦAB=arccosh0rB±arccosh0rA, (70)

or

 h0=rArBrABsin