Frequency Shifts Induced by Orbit Changes

9 Frequency Shifts Induced by Orbit Changes

Improvements in GPS motivate attention to other small relativistic effects that have previously been too small to be explicitly considered. For SV clocks, these include frequency changes due to orbit adjustments, and effects due to earth’s oblateness. For example, between July 25 and October 10, 2000, SV43 occupied a transfer orbit while it was moved from slot 5 to slot 3 in orbit plane F. I will show here that the fractional frequency change associated with a change $da$ in the semi-major axis $a$ (in meters) can be estimated as $9.429 × 10− 18da$ . In the case of SV43, this yields a prediction of $− 1.77 × 10−13$ for the fractional frequency change of the SV43 clock which occurred July 25, 2000. This relativistic effect was measured very carefully [12 ]. Another orbit adjustment on October 10, 2000 should have resulted in another fractional frequency change of $+1.75 × 10− 13$ , which was not measured carefully. Also, earth’s oblateness causes a periodic fractional frequency shift with period of almost 6 hours and amplitude $0.695 × 10−14$ . This means that quadrupole effects on SV clock frequencies are of the same general order of magnitude as the frequency breaks induced by orbit changes. Thus, some approximate expressions for the frequency effects on SV clock frequencies due to earth’s oblateness are needed. These effects will be discussed with the help of Lagrange’s planetary perturbation equations.

Five distinct relativistic effects, discussed in Section 5, are incorporated into the System Specification Document, ICD-GPS-200 [2 ]. These are:

the effect of earth’s mass on gravitational frequency shifts of atomic reference clocks fixed on the earth’s surface relative to clocks at infinity;
the effect of earth’s oblate mass distribution on gravitational frequency shifts of atomic clocks fixed on earth’s surface;
second-order Doppler shifts of clocks fixed on earth’s surface due to earth rotation;
gravitational frequency shifts of clocks in GPS satellites due to earth’s mass;
and second-order Doppler shifts of clocks in GPS satellites due to their motion through an Earth-Centered Inertial (ECI) Frame.

The combination of second-order Doppler and gravitational frequency shifts given in Eq. (27 ) for a clock in a GPS satellite leads directly to the following expression for the fractional frequency shift of a satellite clock relative to a reference clock fixed on earth’s geoid:

where $v$ is the satellite speed in a local ECI reference frame, $GME$ is the product of the Newtonian gravitational constant $G$ and earth’s mass $M$ , $c$ is the defined speed of light, and $Φ0$ is the effective gravitational potential on the earth’s rotating geoid. The term $Φ0$ includes contributions from both monopole and quadrupole moments of earth’s mass distribution, and the effective centripetal potential in an earth-fixed reference frame such as the WGS-84 (873) frame, due to earth’s rotation. The value for $Φ0$ is given in Eq. (18

), and depends on earth’s equatorial radius $a 1$ , earth’s quadrupole moment coefficient $J2$ , and earth’s angular rotational speed $ωE$ .Update

If the GPS satellite orbit can be approximated by a Keplerian orbit of semi-major axis $a$ , then at an instant when the distance of the clock from earth’s center of mass is $r$ , this leads to the following expression for the fraction frequency shift of Eq. (53 ):

Eq. (54

) is derived by making use of the conservation of total energy (per unit mass) of the satellite, Eq. (33

), which leads to an expression for $v2$ in terms of $GME ∕r$ and $GME ∕a$ that can be substituted into Eq. (53

). The first two terms in Eq. (54

) give rise to the “factory frequency offset”, which is applied to GPS clocks before launch in order to make them beat at a rate equal to that of reference clocks on earth’s surface. The last term in Eq. (54

) is very small when the orbit eccentricity $e$ is small; when integrated over time these terms give rise to the so-called “ $e sin E$ ” effect or “eccentricity effect”. In most of the following discussion we shall assume that eccentricity is very small.

Clearly, from Eq. (54 ), if the semi-major axis should change by an amount $δa$ due to an orbit adjustment, the satellite clock will experience a fractional frequency change

The factor $3∕2$ in this expression arises from the combined effect of second-order Doppler and gravitational frequency shifts. If the semi-major axis increases, the satellite will be higher in earth’s gravitational potential and will be gravitationally blue-shifted more, while at the same time the satellite velocity will be reduced, reducing the size of the second-order Doppler shift (which is generally a red shift). The net effect would make a positive contribution to the fractional frequency shift.

Although it has long been known that orbit adjustments are associated with satellite clock frequency shifts, nothing has been documented and up until 2000 no reliable measurements of such shifts had been made. On July 25, 2000, a trajectory change was applied to SV43 to shift the satellite from slot F5 to slot F3. A drift orbit extending from July 25, 2000 to October 10, 2000 was used to accomplish this move. A “frequency break” was observed but the cause of this frequency jump was not initially understood. Marvin Epstein, Joseph Fine, and Eric Stoll [12 ] of ITT evaluated the frequency shift of SV43 arising from this trajectory change. They reported that associated with the thruster firings on July 25, 2000 there was a frequency shift of the Rubidium clock on board SV43 of amount

Epstein et al. [12 ] suggested that the above frequency shift was relativistic in origin, and used precise ephemerides obtained from the National Imagery and Mapping Agency to estimate the frequency shift arising from second-order Doppler and gravitational potential differences. They calculated separately the second-order Doppler and gravitational frequency shifts due to the orbit change. The NIMA precise ephemerides are expressed in the WGS-84 coordinate frame, which is earth-fixed. If used without removing the underlying earth rotation, the velocity would be erroneous. They therefore transformed the NIMA precise ephemerides to an earth-centered inertial frame by accounting for a (uniform) earth rotation rate.

The semi-major axes before and after the orbit change were calculated by taking the average of the maximum and minimum radial distances. Speeds were calculated using a Keplerian orbit model. They arrived at the following numerical values for semi-major axis and velocity:

Since the semi-major axis decreased, the frequency shift should be negative. The prediction they made for the frequency shift, which was based on Eq. (53

), was then

which is to be compared with the measured value, Eq. (56

). This is fairly compelling evidence that the observed frequency shift is indeed a relativistic effect.

Lagrange perturbation theory. Perturbations of GPS orbits due to earth’s quadrupole mass distribution are a significant fraction of the change in semi-major axis associated with the orbit change discussed above. This raises the question whether it is sufficiently accurate to use a Keplerian orbit to describe GPS satellite orbits, and estimate the semi-major axis change as though the orbit were Keplerian. In this section, we estimate the effect of earth’s quadrupole moment on the orbital elements of a nominally circular orbit and thence on the change in frequency induced by an orbit change. Previously, such an effect on the SV clocks has been neglected, and indeed it does turn out to be small. However, the effect may be worth considering as GPS clock performance continues to improve.

To see how large such quadrupole effects may be, we use exact calculations for the perturbations of the Keplerian orbital elements available in the literature [13 ]. For the semi-major axis, if the eccentricity is very small, the dominant contribution has a period twice the orbital period and has amplitude $3J2a21sin2 i0∕(2a0) ≈ 1658 m$ . WGS-84 (837) values for the following additional constants are used in this section: $a = 6.3781370 × 106 m 1$ ; $ω = 7.291151467 × 10 −5 s− 1 E$ ; $a = 2.656175 × 107 m 0$ , where $a 1$ and $a0$ are earth’s equatorial radius and SV orbit semi-major axis, and $ωE$ is earth’s rotational angular velocity.

The oscillation in the semi-major axis would significantly affect calculations of the semi-major axis at any particular time. This suggests that Eq. (33 ) needs to be reexamined in light of the periodic perturbations on the semi-major axis. Therefore, in this section we develop an approximate description of a satellite orbit of small eccentricity, taking into account earth’s quadrupole moment to first order. Terms of order $J2 × e$ will be neglected. This problem is non-trivial because the perturbations themselves (see, for example, the equations for mean anomaly and altitude of perigee) have factors $1∕e$ which blow up as the eccentricity approaches zero. This problem is a mathematical one, not a physical one. It simply means that the observable quantities – such as coordinates and velocities – need to be calculated in such a way that finite values are obtained. Orbital elements that blow up are unobservable.

Conservation of energy. The gravitational potential of a satellite at position $(x,y,z )$ in equatorial ECI coordinates in the model under consideration here is

Since the force is conservative in this model (solar radiation pressure, thrust, etc. are not considered), the kinetic plus potential energy is conserved. Let $ε$ be the energy per unit mass of an orbiting mass point. Then

where $′ V (x,y,z )$ is the perturbing potential due to the earth’s quadrupole potential. It is shown in textbooks [13

] that, with the help of Lagrange’s planetary perturbation theory, the conservation of energy condition can be put in the form

where $a$ is the perturbed (osculating) semi-major axis. In other words, for the perturbed orbit,

On the other hand, the net fractional frequency shift relative to a clock at rest at infinity is determined by the second-order Doppler shift (a red shift) and a gravitational redshift. The total relativistic fractional frequency shift is

The conservation of energy condition can be used to express the second-order Doppler shift in terms of the potential. Since in this paper we are interested in fractional frequency changes caused by changing the orbit, it will make no difference if the calculations use a clock at rest at infinity as a reference rather than a clock at rest on earth’s surface. The reference potential cancels out to the required order of accuracy. Therefore, from perturbation theory we need expressions for the square of the velocity, for the radius $r$ , and for the perturbing potential. We now proceed to derive these expressions. We refer to the literature [13

] for the perturbed osculating elements. These are exactly known, to all orders in the eccentricity, and to first order in $J2$ . We shall need only the leading terms in eccentricity $e$ for each element.

Perturbation equations. First we recall some facts about an unperturbed Keplerian orbit, which have already been introduced (see Section 5). The eccentric anomaly $E$ is to be calculated by solving the equation

where $M$ is the “mean anomaly” and $t0$ is the time of passage past perigee, and

Then, the perturbed radial distance $r$ and true anomaly $f$ of the satellite are obtained from

The observable $x,y,z$ -coordinates of the satellite are then calculated from the following equations:

x = r&#x0028;cos&#x03A9; cos&#x0028;f + &#x03C9; &#x0029; &#x2212; cos isin &#x03A9; sin &#x0028;f + &#x03C9;&#x0029;&#x0029;, &#x0028;67 &#x0029; y = r&#x0028;sin &#x03A9; cos&#x0028;f + &#x03C9;&#x0029; + cosicos &#x03A9; sin &#x0028;f + &#x03C9;&#x0029;&#x0029;, &#x0028;68 &#x0029; z = r&#x0028;sin isin &#x0028;f + &#x03C9;&#x0029;&#x0029;, &#x0028;69 &#x0029;

where $Ω$ is the angle of the ascending line of nodes, $i$ is the inclination, and $ω$ is the altitude of perigee. By differentiation with respect to time, or by using the conservation of energy equation, one obtains the following expression for the square of the velocity:

In these expressions $2 v$ and $−1 r$ are observable quantities. The combination $ecosE$ , where E is the eccentric anomaly, occurs in both of these expressions. To derive expressions for $v2$ and $r−1$ in the perturbed orbits, expressions for the perturbed elements $a,$ $e,$ $E$ are to be substituted into the right-hand sides of the Keplerian equations for $E$ , $r$ , and $v2$ . Therefore, we need the combination $e cosE$ in the limit of small eccentricity.

Perturbed eccentricity. To leading order, from the literature [13 ] we have for the perturbed eccentricity the following expression:

&#x005B;&#x0028; &#x0029; 3J2a21- 3- 2 1- 2 e = e0 + 2a2 1 &#x2212; 2 sin i0 cosf + 4 sin i0cos&#x0028;2&#x03C9;0 + f&#x0029; 0 7 &#x005D; + ---sin2i0cos&#x0028;2&#x03C9;0 + 3f &#x0029; , &#x0028;71 &#x0029; 12

where $e0$ is a constant of integration.

Perturbed eccentric anomaly. The eccentric anomaly is calculated from the equation

with perturbed values for $M$ and $e$ . Expanding to first order in $e$ gives the following expression for $cos E$ :

and multiplying by $e$ yields

We shall neglect higher order terms in $e$ . The perturbed expression for mean anomaly $M$ can be written as

where we indicate explicitly the terms in $e−1 0$ ; that is, the quantity $M 0$ contains all terms which do not blow up as $e → 0$ , and $ΔM ∕e0$ contains all the other terms. The perturbations of $M$ are known exactly but we shall only need the leading order terms, which are

3J2a21&#x005B;&#x0028; 3 2 &#x0029; 1 2 &#x0394;M &#x2215;e0 = &#x2212; -----2 1 &#x2212; -sin i0 sin f &#x2212; --sin i0sin&#x0028;2&#x03C9;0 + f&#x0029; 2e0a 0 2 4 &#x005D; + -7-sin2i sin&#x0028;2&#x03C9; + 3f&#x0029; , &#x0028;76 &#x0029; 12 0 0

and so for very small eccentricity,

Then after accounting for contributions from the perturbed eccentricity and the perturbed mean anomaly, after a few lines of algebra we obtain the following for $ecos E$ :

where the first term is the unperturbed part. The perturbation is a constant, plus a term with twice the orbital period.

Perturbation in semi-major axis. From the literature, the leading terms in the perturbation of the semi-major axis are

where $a0$ is a constant of integration. The amplitude of the periodic term is about 1658 meters.

Perturbation in radius. We are now in position to compute the perturbation in the radius. From the expression for $r$ , after combining terms we have

r = a0 &#x0028;1 &#x2212; e0 cosE0 &#x0029; + &#x0394;a &#x2212; &#x0394; &#x0028;ecos E&#x0029; 3J a2 &#x0028; 3 &#x0029; = a0 &#x0028;1 &#x2212; e0 cosE0 &#x0029; &#x2212; ---2-1 1 &#x2212; -sin2i0 2a0 2 J2a21 2 + -4a--sin i0cos 2&#x0028;&#x03C9;0 + f&#x0029;. &#x0028;80 &#x0029; 0

The amplitude of the periodic part of the perturbation in the observable radial distance is only 276 meters.

Perturbation in the velocity squared. The above results, after substituting into Eq. (70 ), yield the expression

v2 GME 3GMEJ2a21 &#x0028; 3 2 &#x0029; ---= -----&#x0028;1 + 2e0 cosE0 &#x0029; +------3---- 1 &#x2212; --sin i0 2 2a0 2 a0 2 GMEJ2a---1 2 + 2a30 sin i0 cos2&#x0028;&#x03C9;0 + f&#x0029;. &#x0028;81 &#x0029;

Perturbation in $GME ∕r$ . The above expression for the perturbed $r$ yields the following for the monopole contribution to the gravitational potential:

&#x0028; &#x0029; GME--- GME--- 3GMEJ2a21-- 3- 2 &#x2212; r = &#x2212; a0 &#x0028;1 + e0cos E0&#x0029; &#x2212; 2a30 1 &#x2212; 2 sin i0 2 2 + GMEJ2a--1sin--i0-cos2&#x0028;&#x03C9;0 + f &#x0029;. &#x0028;82 &#x0029; 4a30

Evaluation of the perturbing potential. Since the perturbing potential contains the small factor $J2$ , to leading order we may substitute unperturbed values for $r$ and $z$ into $V ′(x,y,z)$ , which yields the expression

Conservation of energy. It is now very easy to check conservation of energy. Adding kinetic energy per unit mass to two contributions to the potential energy gives

This verifies that the perturbation theory gives a constant energy. The extra term in the above equation, with $J 2$ in it, can be neglected. This is because the nominal inclination of GPS orbits is such that the factor $2 (1 − 3 sin i0∕2)$ is essentially zero. The near vanishing of this factor is pure coincidence in the GPS. There was no intent, in the original GPS design, that quadrupole effects would be simpler if the orbital inclination were close to $55∘$ . However, because this term is negligible, numerical calculations of the total energy per unit mass provide a means of evaluating the quantity $a 0$ .

Calculation of fractional frequency shift. The fractional frequency shift calculation is very similar to the calculation of the energy, except that the second-order Doppler term contributes with a negative sign. The result is

&#x0394;f v2 GME V &#x2032; ----= &#x2212; ---2 &#x2212; --2---+ -2- f 2c c r c 2 &#x0028; &#x0029; 3GME--- 2GME--- 7GMEJ2a--1- 3- 2 = &#x2212; 2a0c2 &#x2212; c2a0 e0 cosE0 &#x2212; 2a30c2 1 &#x2212; 2 sin i0 2 2 &#x2212; GMEJ2a---1sin-i0cos 2&#x0028;&#x03C9;0 + f&#x0029;. &#x0028;85 &#x0029; a30c2

The first term, when combined with the reference potential at earth’s geoid, gives rise to the “factory frequency offset”. The seond term gives rise to the eccentricity effect. The third term can be neglected, as pointed out above. The last term has an amplitude

which may be large enough to consider when calculating frequency shifts produced by orbit changes. Therefore, this contribution may have to be considered in the future in the determination of the semi-major axis, but for now we neglect it.

The result suggests the following method of computing the fractional frequency shift: Averaging the shift over one orbit, the periodic term will average down to a negligible value. The third term is negligible. So if one has a good estimate for the nominal semi-major axis parameter, the term $2 − 3GME ∕2a0c$ gives the average fractional frequency shift. On the other hand, the average energy per unit mass is given by $ε = − GME ∕2a0$ . Therefore, the precise ephemerides, specified in an ECI frame, can be used to compute the average value for $ε$ ; then the average fractional frequency shift will be

The last periodic term in Eq. (85

) is of a form similar to that which gives rise to the eccentricity correction, which is applied by GPS receivers. Considering only the last periodic term, the additional time elapsed on the orbiting clock will be given by

where to a sufficient approximation we have replaced the quantity $f$ in the integrand by $∘ -------3- n = GME ∕a0$ ; $n$ is the approximate mean motion of GPS satellites. Integrating and dropping the constant of integration (assuming as usual that such constant time offsets are lumped with other contributions) gives the periodic relativistic effect on the elapsed time of the SV clock due to earth’s quadrupole moment:

The correction that should be applied by the receiver is the negative of this expression,

The phase of this correction is zero when the satellite passes through earth’s equatorial plane going northwards. If not accounted for, this effect on the SV clock time would give rise to a peak-to-peak periodic navigational error in position of approximately $2c × δtJ2 = 1.43 cm$ .

These effects were considered by Ashby and Spilker [9 ], pp. 685–686, but in that work the effect of earth’s quadrupole moment on the term $GME ∕r$ was not considered; the present calculations supercede that work.

Numerical calculations. Precise ephemerides were obtained for SV43 from the web site ftp://sideshow.jpl.nasa.gov/pub/gipsy_products/2000/orbits at the Jet Propulsion Laboratory. These are expressed in the J2000 ECI frame. Computer code was written to compute the average value of $ε$ for one day and thence the fractional frequency shift relative to infinity before and after each orbit change. The following results were obtained:

07&#x2215;22&#x2215;00 : a = &#x0028;2.65611575 &#x00D7; 107 &#x00B1; 69&#x0029; m, 3 07&#x2215;30&#x2215;00 : a = &#x0028;2.65423597 &#x00D7; 10 &#x00B1; 188&#x0029; m, 10&#x2215;07&#x2215;00 : a = &#x0028;2.65418742 &#x00D7; 107 &#x00B1; 95&#x0029; m, 7 10&#x2215;12&#x2215;00 : a = &#x0028;2.65606323 &#x00D7; 10 &#x00B1; 58&#x0029; m.

Therefore, the fractional frequency change produced by the orbit change of July 25 is calculated to be

which agrees with the measured value to within about 3.3%. The agreement is slightly better than that obtained in [12

], perhaps because they did not consider contributions to the energy from the quadrupole moment term.

A similar calculation shows that the fractional frequency shift of SV43 on October 10, 2001 should have been

No measurement of this shift is available.

On March 9, 2001, SV54’s orbit was changed by firing the thruster rockets. Using the above procedures, I can calculate the fractional frequency change produced in the onboard clocks. The result is

Using Eq. (55

) yields the following prediction for the fractional frequency change of SV54 on March 9, 2001:

The quoted uncertainty is due to the combined uncertainties from the determination of the energy per unit mass before and after the orbit change. These uncertainties are due to neglecting tidal forces of the sun and moon, radiation pressure, and other non-gravitational forces.

Summary. We note that the values of semi-major axis reported by Epstein et al. [12 ] differ from the values obtained by averaging as outlined above, by 200–300 m. This difference arises because of the different methods of calculation. In the present calculation, an attempt was made to account for the effect of earth’s quadrupole moment on the Keplerian orbit. It was not necessary to compute the orbit eccentricity. Agreement with measurement of the fractional frequency shift was only a few percent better than that obtained by differencing the maximum and minimum radii. This approximate treatment of the orbit makes no attempt to consider perturbations that are non-gravitational in nature, e.g., solar radiation pressure. The work was an investigation of the approximate effect of earth’s quadrupole moment on the GPS satellite orbits, for the purpose of (possibly) accurate calculations of the fractional frequency shifts that result from orbit changes.

As a general conclusion, the fractional frequency shift can be estimated to very good accuracy from the expression for the “factory frequency offset”.