Science Resultsâ NASA Final Report
December 2008
ii
Gravity Probe B Science ResultsâNASA Final Report, December 2008
Signatures & Approvals
Tom Langenstein: ITAR Assessment performed. ITAR Control required: ____Yes _
X
_No
Prepared by
Robert Kahn
Public Affairs Coordinator,
Stanford University
December 31, 2008
Date
Approved by
Francis Everitt
Principal Investigator,
Stanford University
December 31, 2008
Date
Approved by
Barry Muhlfelder
Program Manager
Stanford University
December 31, 2008
Date
Approved by
Tom Langenstein
Deputy Program Manager
Stanford University
December 31, 2008
Date
Contents
iii
Table of Contents
iv
Gravity Probe B Science ResultsâNASA Final Report, December 2008
Preface
1
Preface
This Final Scientific Report to NASA on the Gravity Probe B (GP-B) Relativity Mission is written to
accompany the Post Flight Analysis Final Report (March 2007, submitted September 2007). It covers
results of data analysis from 28 August 2004 when science data gathering began to 30 September 2008 when
NASA support under Contract NAS8-39225 ended. During the last period, 15 January 2008 through 30
September 2008, funding was received from three distinct sources:
1. "NASA ($500K to Stanford University, $133K to NASA MSFC)
2. "Mr. Richard Fairbank through a donation of $512K to Stanford University
3. "Stanford University ($350K directly to the program, ~$220K indirectly through waiving normal
University overhead on the NASA contract)
The objective of GP-B was to design, develop, conduct, and analyze the data of a flight experiment to test
two predictions of Albert Einstein's General Relativity Theory: 1)the space-time curvature in the vicinity of
the Earth caused by the Earth's mass (geodetic effect), and 2) the dragging of space-time caused by the
rotation of the Earth (frame dragging effect). Both have now been clearly measured. The NS drift
component yields a consistent determination of the geodetic effect to 0.5%. The frame-dragging effect is
plainly visible in the processed data with a present statistical uncertainty of 15%. See the reconstructed
relativistic East-West orientations of Figure 6 in the fourth chapter of this Report by M. Heifetz et al.
The data analysis leading up to this important result proved more subtle than expected. 'Patch-effect'
1
anomalies on the gyro rotor and housing have complicated the gyro behavior in two ways, a changing
polhode path affecting the determination of the gyro scale factor, and two larger than expected Newtonian
torques. Put simply, while mechanically both rotor and housing are exceedingly spherical, electrically they
are not. Steadily advancing progress, reported to NASA directly and via successive meetings of the NASA
appointed GP-B Science Advisory Committee (SAC) chaired by Professor Clifford Will, has brought a
rather complete understanding of these effects. A turning point came in August 2008 with the
incorporation of an elegant approach for computing the detailed history of the remarkable "roll-polhode
resonance" torques discovered a year earlier by J. Kolodziejczak of NASA MSFC. The result, reported here
2
,
was a large reduction in previously unexplained discrepancies between the four gyroscopes.
Much further work remains to bring the analysis to completion. To date, limits in computational power
have bounded the processing to essentially one point per 97-minute GP-B satellite orbit. The driving
period of the roll-polhode resonance torques is at the difference between the 77.5 sec roll period of the
spacecraft and a harmonic of the gyroscope polhode period. High-speed computing techniques now in
process will lead to more detailed analyses, and allow GP-B to approach the intrinsic limit of the gyro
readout.
Funding for this final effort has been provided from Saudi Arabia through the King Abdulaziz City for
Science and Technology (KACST) as part of a recently established Stanford University/KACST
collaboration with Professor Charbel Farhat of the Stanford Aero-Astro Department as Co-PI. It is
anticipated that completion of this more detailed study, together with publication of a number of
1. 'Patch effect' is the name given to contact potential differences between different crystalline or contamination regions on a metallic surface
2. The present report is based upon the texts of papers submitted for publication in the international, refereed journal Space Science Reviews, and reprinted as hardcover
book in the Space Sciences Series of ISSI (International Space Science Institute), both published by Springer. These papers follow from a set of invited talks given October
6-10, 2008, by five GP-B team members at the Workshop "The Nature of Gravity: Confronting Theory and Experiment in Space" held at ISSI in Bern, Switzerland.
2
Gravity Probe B Science ResultsâNASA Final Report, December 2008
technology papers, will be approximately June 2010. To maximize the benefit to the scientific and
engineering community, and NASA, we plan to make the capstone of the program a conference on
Fundamental Physics and Innovative Engineering in Space, in honor of William M. Fairbank who was one
of the three vital founders of Gravity Probe B.
We thank NASA for forty-four years of continued support since issuing the first research Grant NSG-582 to
the program in March 1964. The March 2007 Report contained an extensive history of GP-B and the NASA
personnel who guided it. It is appropriate here to express further special thanks to three individuals, the
MSFC Manager Mr. Anthony T. Lyons, the HQ Program Scientist for Physics of the Cosmos Dr. Michael H.
Salamon, and the HQ Program Executive Dr. Alan P. Smale. The SAC will continue to meet with us
regularly; we are glad that Dr. Salamon has agreed to participate in its meetings.
C.W. Francis Everitt
Bradford W. Parkinson
Principal Investigator
Co-Principal Investigator
December 2008
3
1
G
ravity
P
robe
b D
ata
a
nalysis
â
s
tatus
anD
P
otential
for
i
mProveD
a
ccuracy
of
s
cientific
r
esults
Francis Everitt et al.
4
Gravity Probe B Science ResultsâNASA Final Report, December 2008
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results
5
Gravity Probe B Data AnalysisâStatus and Potential
for Improved Accuracy of Scientific Results
C.W.F. Everitt, M. Adams, W. Bencze, S. Buchman, B. Clarke, J. Conklin, D.B. DeBra, M. Dolphin,
M. Heifetz, D. Hipkins, T. Holmes, G. M. Keiser, J. Kolodziejczak*, J. Li, J. Lipa, J.M. Lockhart,
B. Muhlfelder, Y. Ohshima, B.W. Parkinson, M. Salomon, A. Silbergleit, V. Solomonik, K. Stahl,
M. Taber, J.P. Turneaure, S. Wang, P.W. Worden Jr.
W.W. Hansen Experimental Physics Laboratory, Stanford University, Stanford, CA 94305-4085
*
NASA Marshall Space Flight Center, Huntsville, AL 35812
Abstract
This is the first of five connected papers giving details of progress on the Gravity Probe B (GP-B)
Relativity Mission. GP-B, launched on 20 April 2004, is a landmark physics experiment in space to test
two fundamental predictions of Einsteinâs theory of general relativity, the geodetic and frame-dragging
effects, by means of ultra-precise cryogenic gyroscopes in Earth orbit. Data collection began 28 August
2004 and science operations were completed 29 September 2005. The data analysis has proven deeper
than expected. Patch effect anomalies on the gyro rotor and housing have given rise to two mutually
reinforcing complications:
1)
A changing polhode path affecting the calibration of the gyroscope scale factor
C
g
against the
aberration of starlight
2)
Two larger than expected manifestations of a Newtonian gyro torque due to patch potentials on
the rotor and housing.
In earlier papers, we reported two distinct methods, âgeometricâ and âalgebraicâ, for identifying and
removing the first Newtonian effect (âmisalignment torqueâ), and also a preliminary method of treating
the second (âroll-polhode resonance torqueâ). Central to the progress in both torque modeling and
C
g
determination has been an extended effort from January 2007 on âTrapped Flux Mappingâ. A new
turning point came in August 2008 when it became possible to include a detailed history of the roll-
polhode resonance torques into the computation. The frame-dragging effect is now plainly visible in the
processed data. The current statistical uncertainty from an analysis of 154 days of data by the algebraic
method is 6 marcsec/yr (~15% of the frame-dragging effect). The systematic error will be added to this
statistical uncertainty using the methods discussed in an accompanying paper by Muhlfelder et al. A
covariance analysis, incorporating the impact of patch effect anomalies, indicates that a 3 to 5%
determination of frame-dragging is possible with a more complete, but computationally intensive data
analysis approach.
1 The NASA-Stanford Gravity Probe B Mission (GP-B)
In 1960, L. I. Schiff (Schiff (1960)) showed that an ideal gyroscope in orbit around the Earth or other
massive body would undergo two relativistic precessions with respect to a distant undisturbed inertial
frame: 1) a geodetic precession in the orbit plane due to the local curvature of space time; 2) a frame-
dragging precession due to the Earthâs rotation.
6
Gravity Probe B Science ResultsâNASA Final Report, December 2008
Figure 1.
Predicted precessions of the GP-B gyroscope (North and West are positive in this coordinate system).
The geodetic effect is fundamentally equivalent to the curvature precession of the Earth-Moon system
around the Sun first given by W. De Sitter in 1916. The Schiff frame-dragging effect is related to, but not
identical with, the dragging of the orbit plane of a satellite around a rotating planet computed in 1918 by
J. Lense and H. Thirring (Lense and Thirring 1918). It has, as Thirring observed, connections to the
possible relation between General Relativity and Machâs principle.
Schiffâs formula for the combined gyroscope precession
âŚ
in Einsteinâs theory is:
(
)
(
)
âĽâŚ
â¤
â˘âŁ
âĄ
â
â
+
Ă
=
Ď
R
Ď
R
v
R
2
3
2
3
2
3
2
3
R
R
c
GI
R
c
GM
âŚ
where the first term is the geodetic, and the second the frame-dragging effect, with
G
the gravitational
constant,
c
the velocity of light,
M
,
I
,
Ď
, the mass, moment of inertia, and angular velocity of the
Earth, and
R
and
v
the instantaneous distance and velocity of the gyroscope. In the polar orbit selected
for Gravity Probe B, the two effects are at right angles. Their magnitudes are found by integrating
Schiffâs formula around the orbit. For the geodetic effect, it is necessary to take into account the
oblateness of the Earth which modifies both the orbit and the term itself; a correction investigated first by
D. C. Wilkins, then B. M. Barker and R. F. OâConnell (Barker and O'Connell (1970)), and finally, most
elegantly, by J. V. Breakwell (Breakwell (1988)). In the 642 km orbit of GP-B, the predicted geodetic
precession is -6,606 marc-s/yr and the frame-dragging -39 marc-s/yr. GP-B was designed to measure
both to < 0.5 marc-s/yr.
The gyroscope is a spinning spherical body. Conceptually, therefore, Gravity Probe B is simple. All it
needs is a star, a telescope, and a spinning sphere. The difficulty lies in the numbers. To reach the
0.5 marc-s/yr experiment goal calls for:
1)
One or more exceedingly accurate gyroscopes with drift rates < 10
-11
deg/hr, i.e. 6 to 7 orders of
magnitude better than the best modeled inertial navigation gyroscopes
2)
A reference telescope ~3 orders of magnitude better than the best previous star trackers
3)
A sufficiently bright suitably located guide star (IM Pegasi was chosen) whose proper motion
with respect to remote inertial space is known to <0.5 marc-s/yr
4)
Sufficiently accurate orbit information to calibrate the science signal and calculate the two
predicted effects
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results
7
These four constraints led to 10 fundamental requirements on the design of GP-B, and then to 27 distinct
experiment subsystems for executing the mission.
Figure 2.
The Gravity Probe B Instrument.
The heart of the instrument is the integrated fused quartz structure, operating at cryogenic temperatures,
illustrated in Figure 2. It comprises a 0.14m aperture, 3.81m focal length folded Cassegrainian telescope
bonded to a quartz block containing four gyroscopes mounted in line to within 0.1mm of the boresight of
the telescope. Two gyroscopes spin clockwise and two counter-clockwise with their axes initially aligned
to the star to within a few arc-s. Each gyroscope measures both relativity effects. The spacecraft rolls
slowly (77.5s period) around the line to the star. Since the guide star is occulted for almost half of each
orbit (see Figure 1), the telescope must reacquire it every orbit.
The instrument is mounted in an evacuated chamber within a 2440l superfluid helium dewar (1.8 K)
fabricated by Lockheed Martin (Figure 3) (Parmley, Goodman et al. 1986).
Figure 3.
The Flight Dewar.
The boil-off gas provides thrust for attitude, translational and roll control of the spacecraft. With a
pointing of 150-250 marc-s, reaching the < 0.5 marc-s/yr experiment goal requires accurate subtraction of
the gyroscope signal
G
from the telescope signal
T
. Thus, in addition to near perfect telescope and
Quartz
Block
Gyro 1
Gyro 2
Reference
Telescope
Guide
Star
Gyro 3 Gyro 4
Mounting Flange
8
Gravity Probe B Science ResultsâNASA Final Report, December 2008
gyroscopes, we face the challenge of accurately matching the two scale factors to make valid
G
-
T
subtractions, and as will appear, an interesting separate process of calibration for
G
.
The gyroscope (Figure 4) comprises a near perfect niobium-coated fused quartz sphere (3.81 cm dia)
electrically suspended about its geometric center in a spherical quartz housing.
Figure 4.
Gravity Probe B gyroscope.
A simple argument reveals the benefit of space. Let the maximum allowed Newtonian drift rate be
0
Ί
(~2
Ă
10-17 rad/s for a < 0.5 marc-s/yr experiment). Picture a perfectly spherical, but not perfectly
homogeneous, rotor with its mass center a distance
δ
r
from the geometric center. An acceleration
f
at
right angles to the spin axis will cause a torque
T=
M(
δ
r
Ă
f
)
yielding the requirement
0
5
2
Ί
<
s
v
r
r
f
δ
,
where
s
v
is the peripheral velocity of the rotor (~10m/s). Inserting numbers, one finds for an experiment
on Earth (f = 1g) the ridiculous homogeneity requirement 6
Ă
10
-18
. For a satellite, the acceleration level
from air drag, solar radiation pressure, etc. might typically be f ~10
-8
g which yields the still exceedingly
difficult requirement ~6
Ă
10
-10
. This brings us to the âdrag-freeâ concept first suggested in 1959 by G. E.
Pugh (Pugh (1959)), and then extensively developed at Stanford University. The gyroscope, being
shielded from drag, tends to follow an ideal gravitational orbit. By sensing its position inside the satellite
and applying compensating thrust, the net acceleration is further reduced to f ~10
-11
g. The resulting
homogeneity requirement is 6
Ă
10
-7
; the homogeneity of the actual rotors was 3
Ă
10
-7
. A parallel
argument holds for suspension torques acting on the out-of-roundness of the gyro rotor, and is met
experimentally with an equally satisfactory margin.
Having a rotor spherical and mass-balanced to 10 nm poses a readout problem; for how is one to measure,
with marc-s accuracy, the direction of spin of a perfectly round, perfectly homogeneous unmarked
sphere? The key is cryogenics. Below 9K, the rotorâs niobium coating is superconducting; it generates a
small dipole magnetic moment, the London moment (
M
L
) (London (1960)) proportional to spin, aligned
with the instantaneous spin axis.
Detection is by a highly sensitive SQUID magnetometer coupled to a 4-turn thin film superconducting
coil patterned on the flat surface of the left-side housing half (Figure 4). Figure 5 indicates how the
amplitude and phase of the misalignment angle are determined at the spacecraftâs 77.5 s roll period
(Muhlfelder, Lockhart et al.(2003)).
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results
9
Figure 5.
London Moment Readout
Because the rotor is cooled initially in a finite field, superimposed on the London moment and
complicating the readout process is some trapped magnetic flux. See Section 3 below.
Cryogenics helps in three other ways: magnetic shielding, ultra-high vacuum, and mechanical stability.
â˘
Magnetic Shielding: The shielding has two formidable requirements, (
i
) a dc field level of 10
-7
gauss to reduce trapped flux in the rotor to ~1% of the London moment, and (
ii
) 240db ac
attenuation to isolate the gyroscopes from varying external fields. The 10
-7
gauss dc field was
obtained by an expanding lead-bag technique devised by B. Cabrera (Cabrera (1982))and applied
to GP-B by M. Taber, J. Mester, J. Lockhart, and colleagues (Mester, Lockhart et al. (2000)); ac
shielding was implemented with a succession of nested Cryoperm and superconducting shields
including transverse cylindrical shields around each gyroscope. Among the engineering
challenges, two were critical, inserting the warm probe into a cold dewar, and ensuring that the
probe itself was non-magnetic. Insertion was a severely disciplined 24-hour-long process
requiring a sophisticated airlock (or rather helium lock). As for the probe itself, this took an even
more disciplined process in which every component, tool, and manufacturing procedure was
subject to appropriate limits and testing. An engineer in good position to judge (R. Parmley of
Lockheed Martin) has estimated that 30% of the cost came from meeting the magnetic
requirements (Parmley, Bell et al. (2003)).
â˘
Ultra-high vacuum: To spin a gyroscope, a torque must be applied. For GP-B, the only viable
technique is by helium gas run at sonic velocity through the differentially pumped channel in the
right-side housing half (Figure 4). The spin-up gas and torque are then switched off. Let the spin
torque be
T
s
and its residual value after spin-up be
T
r
. With spin-up time
t
s
and the previously
stated Newtonian drift limit
0
Ί
,
T
r
/
T
s
has to be <
t
s
0
Ί
, i.e. <10
-13
, another formidable
requirement and one that makes cryogenic operation seem at first disadvantageous since at a
given pressure, gas torques increase as the temperature is lowered. The key is âlow-temperature
bakeoutâ (Turneaure, Cornell et al. (1982)). At 2.5K, everything except He
3
and He
4
is frozen
out. By raising the temperature from 2.5K to 6K for 10 hours, any adsorbed gas is driven off
from the exposed surfaces, first allowed to be exhausted to space and then taken up in a large
cryopump after the valve to space is closed. Tests on orbit showed that the pressure in the
chamber after bakeout was <10
-14
torr, perhaps the lowest ever attained in a scientific instrument.
â˘
Mechanical stability: Two obvious concerns have been: (
i
) possible thermal distortion of the
telescope/quartz block structure, and (
ii
) viscoelastic creep in the shape of the gyro rotor under
the surprisingly large centrifugal acceleration from spin. Both matter at room temperature; both
nearly vanish at 2.5K. At the 230K ambient satellite temperature, sunlight falling on the side of
the telescope would shift the star image many arc-s. The greatly reduced expansion coefficient,
and greatly improved thermal isolation in the GP-B dewar, together eliminate any significant
image motion. Likewise with creep. Metal rotors at room temperature have been known to
change shape significantly after 100 days. Fused quartz at 2.5K has a viscoelastic time well in
10 Gravity Probe B Science ResultsâNASA Final Report, December 2008
excess of 1000 years; no detectable change can occur. One further interesting aspect of cryogenic
operation is the reduced thermal relaxation time; the GP-B gyro rotor equilibrates thermally in
seconds, not minutes.
Issues of mechanical stability are, as H. A. Rowland remarked 120 years ago, the nightmare of precision
measurement. In GP-B, not only do cryogenics help; so does space. On Earth, the GP-B telescope
cantilevered under its own weight sags causing a 130 marc-s displacement of the image. In space, such
distortions vanish.
A focused 0.14m aperture diffraction-limited telescope has a 1.8 arc-s Airy disk corresponding at 3.81m
focal length to an image diameter of ~30
Îź
m. Measurement to < 0.5 marc-s requires image division at the
nm level, yielding two main design constraints, one on pointing, the other on the sharpness of the image
divider. First, two star images are formed by means of a beam splitter. These then separately fall on the
edges of two roof prisms to give reference signals in two orthogonal directions from the intensities of the
beams reflected from the prisms. The criteria determining performance are photon noise and the
sharpness of the prism edge. For IM Pegasi (5.9 magnitude), the effective pointing accuracy determined
by photon noise was ~50 marc-s. With that, as already remarked, an accurate gyro-telescope subtraction
(
G
-
T
) becomes necessary. Two separate methods have been used. In the first, dither signals of known
periods (~30s & ~40s) and amplitudes ~ 60 marc-s were injected into the pointing servo, and then
demodulated. In the second, correlation functions were drawn on the natural pointing variations. The
agreement between the two methods on orbit was excellent.
The output of the London moment readout is a voltage which, though proportional to angle, depends on
several experimental parameters and hence is not directly in angular measure. Nature itself provides the
calibration signal. The apparent position of the guide star varies by known amounts through the Earthâs
motion around the Sun (annual aberration) and the spacecraftâs motion around the Earth (orbital
aberration). Because the telescope is held pointing at the guide star, these aberration signals appear in the
output of the gyroscope. They provide exact cross-checking calibrations of the gyro scale factor
G
.
Recalling that the line of sight to the star is occulted by the Earth for nearly half of each orbit (Figure 1),
we define two parts for each orbit, Guide Star Valid (GSV) and Guide Star Invalid (GSI). Figure 6
illustrates data from one GSV period. The 5 arc-s orbital aberration envelope is clearly visible
superimposed on the output of the gyroscope modulated at the 77.5s roll period.
Figure 6.
Gyroscope London Moment Data Over One Guide Star Valid Period.
Peak to peak ~ 24 arc-sec
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results 11
Gravity Probe B was launched from Vandenberg Air Force Base, California on 20 April 2004. The dewar
stayed cold on orbit for 17 months and 9 days, surpassing the 16½ month design lifetime. Initial Orbital
Checkout (IOC) took 128 days during which time extensive calibration tests re-verified numerous pre-
launch performance measures. Spin-up of the four gyroscopes was completed on 14 July 2004; science
data collection began 28 August 2004. A final Post-Science Calibration Phase, starting 14 August 2005
and ending with the depletion of the dewar on 29 September 2005, involved tests in which certain
potential systematic effects were deliberately enhanced. These included, for example, pointing the
spacecraft for limited periods at other real or virtual guide stars to investigate the magnitude of
misalignment torques.
2 The Actual vs. the Ideal Gravity Probe B
Our aim for GP-B was an apparatus where all disturbances would be vanishingly small so that all four
gyroscopes would agree to within the <0.5 marc-s/yr design requirement. Since the gyroscopes have
different mechanical and electrical parameters, and function in different magnetic, electric, and
acceleration environments, that four-way agreement, if achieved, would be exceedingly powerful. Other
fundamentals included rolling the spacecraft about the line to the star for torque averaging and removing
1/f noise in the readouts, and the expectation of <10
-12
g cross-axis average of the drag-free control, and
<5 marc-s pointing accuracy at roll period. With a series of ground-based and IOC tests, along with
known natural variations during the Science Phase and calibrated increases of potential systematic effects
during the Calibration Phase, GP-B would become the most rigorously validated of all tests of Einsteinâs
theory. That aim, though attained by different means in the actual as against the ideal GP-B, is one we
still hold to.
Ideally, also there would be no data interruptions other than the once per orbit occultations of the guide
star. With relativity signals increasing as time
t
, and readout noise averaging as
t
-½
, the measurement
accuracy should then advance as
t
-3/2
, a result slightly modified by the process of calibrating the gyroscope
scale factor against the aberration of starlight. We knew, however, that interruptions of data would
probably occur, and had a provision for dealing with them (Duhamel (1984)).
Many of the ideal conditions were met. The performances of the telescope and dewar met or beat
expectation. The launch vehicle delivered GP-B into, as close as could be imagined, an ideal orbit. The
gyroscopes were nearly 10
6
times better than the best inertial navigation gyroscopes. The UV discharge
of electric charge on the rotors proved easier than expected. The low field and trapped flux levels
remained constant through the year. The SQUID noise in the gyro readouts was actually less than the
pre-launch expectation, reaching levels corresponding to ideal measurements ranging from 0.14 marc-s/yr
for Gyro 3 to 0.35 marc-s/yr for Gyro 4. Nevertheless, some deviations from the ideal did occur, of
which the most important were two mutually reinforcing complications in gyro behavior, a changing
polhode period affecting the determination of the scale factor, and two larger than expected forms of
Newtonian torque, known respectively as the âmisalignmentâ and âroll-polhode resonanceâ torques.
Figure 7 (April 2007) shows what is effectively raw data from the four gyroscopes in the North-South
(orbital) plane. The geodetic effect is at once obvious, but so are a number of unexplained wiggles. The
mean 1
Ď
result is
â
6,638 Âą 97 marc-s/yr (Everitt (2007)). After subtracting corrections for the solar
geodetic effect (+7 marc-s/yr) and the proper motion of the guide star (+28 Âą 1 marc-s/yr), the result is
â
6,673 Âą 97 marc-s/yr, to be compared with the predicted
â
6,606 marc-s/yr of General Relativity. The
East-West results required for frame-dragging were at this stage inconclusive. Remarkably, it has proved
possible to provide detailed understanding of all the more important disturbing effects, and provide
Newtonian methods for rigorously removing them.
12 Gravity Probe B Science ResultsâNASA Final Report, December 2008
10/30
12/19
02/07
03/29
05/18
07/07
-2
0
2
4
ar
cs
ec
Gyro 2. NS Inertial Orientation
10/30
12/19
02/07
03/29
05/18
07/07
-2
0
2
4
ar
cs
ec
Gyro 2. NS Inertial Orientation
Gyro 2: NS Inertial Orientation
10/30
12/19
02/07
03/29
05/18
07/07
0
2
4
ar
cs
ec
Gyro 4. NS Inertial Orientation
10/30
12/19
02/07
03/29
05/18
07/07
0
2
4
ar
cs
ec
Gyro 4. NS Inertial Orientation
Gyro 4: NS Inertial Orientation
10/30
12/19
02/07
03/29
05/18
07/07
-10
-8
-6
ar
cs
ec
Gyro 1. NS Inertial Orientation
10/30
12/19
02/07
03/29
05/18
07/07
-10
-8
-6
ar
cs
ec
Gyro 1. NS Inertial Orientation
Gyro 1: NS Inertial Orientation
10/30
12/19
02/07
03/29
05/18
07/07
0
2
4
6
ar
cs
ec
Gyro 3. NS Inertial Orientation
10/30
12/19
02/07
03/29
05/18
07/07
0
2
4
6
ar
cs
ec
Gyro 3. NS Inertial Orientation
Gyro 3: NS Inertial Orientation
Figure 7.
Direct measure of the North-South (NS) geodetic precession.
The original ideal <10
-11
deg/hr GP-B gyroscope depended on making seven distinct quantities
simultaneously ânear zeroâ. Three, homogeneity, sphericity, and electric dipole moment were properties
of the rotor. Three others, drag-free acceleration, pressure, and magnetic field level, were properties of
the environment. The seventh concerned maintaining the electric charge on the rotor (due either to lift-off
or cosmic ray impacts) at an acceptably low level. Table 1 gives the requirements set before launch.
Detailed on-orbit investigations showed that six of the seven were met, or even surpassed, but with regard
to the electric dipole ânear zeroâ, there was an issue, or to be more exact, a double issue due to the patch
effect.
Rotor Properties
â˘
Density homogeneity
<6
Ă
10
-7
met
â˘
Sphericity
< 10 nm
met
â˘
Electric dipole moment
< 0.1 V-m
issue
Environment
â˘
Cross track acceleration
< 10
-11
g
met
â˘
Gas pressure
< 10
-12
torr
met
â˘
Magnetic field
< 10
-6
gauss
met
Mixed
â˘
Rotor electric charge
< 10
8
electrons
met
Table 1
The Original Seven Near Zeros
âPatch effectâ is the name given to contact potential differences between different crystalline or
contamination regions on a metal surface (Darling (1989)). The concern prior to launch was of patches
on the rotor creating a net overall dipole moment which would interact with the gyro suspension voltages
to cause a torque. Kelvin probe measurements on flat samples indicated crystals so minute that any such
effect would vanish. This conclusion was, as later UV scanning measurements on one flight-quality rotor
revealed, incorrect, but gradually, more important, it became evident that individual patches on the rotor
could interact with nearby patches on the housing causing significant disturbing forces. Put simply, while
mechanically both rotor and housing are exceedingly spherical, electrically they are not. These patch
effect terms are now known to explain the two classes of anomalous Newtonian torques, and quite
probably also the changing polhode period just referred to.
Pre-launch calculations had given strong reason that the polhode periods of the gyroscopes would remain
fixed throughout the mission. The discovery, a few days into the Science Phase, illustrated in Figure 8,
that they were changing was a considerable surprise. With Gyros 1 and 2, the initial 1.8 and 6 hour
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results 13
periods increase rapidly, reaching separatrices with extremely high periods after 20 to 30 days, and then
decrease progressively to ~1 and ~3 hour asymptotic values over the rest of the year. They began by
spinning at or near their minimum moments of inertia, pass through their intermediate moments at the
separatrix, and finally approach spinning about their maximum moments. With Gyros 3 and 4 initially
spinning more nearly around their maximum moments, there is a steady decrease to asymptotic periods of
1.55 and 4.2 hours. Observe that while the body axes are going through these large motions, the angular
momentum of the gyroscope is unaffected.
Figure 8.
Polhode period history of GP-B gyros: red from HF FFT, blue from snapshots.
Why does a changing polhode period matter? The answer comes in the process of calibrating the gyro
scale factor against the aberration of starlight. To perform a <0.5 marc-s/yr measurement in the presence
of a 6.6 arc-s/yr geodetic effect, 5 arc-s orbital aberration, and 20 arc-s annual aberration, one needs a
readout scaled correctly to 1-2 parts in 10
5
. Now the signal used to provide that calibration is the SQUID
readout data itself, which is taken during the half-orbit GSV segments. Reaching the ~10
-5
accuracy
necessarily depends on making exact connections between the data from successive half-orbits. The
changing polhode period complicates this process owing to the presence of the trapped flux in the rotor.
Details of the two methods of dealing with this problem are given in the accompanying paper, âPolhode
Motion, Trapped Flux, and the GP-B Science Data Analysisâ by A. Silbergleit et al. (Silbergleit et al.
(2009)).
The accompanying paper by G. M. Keiser et al. âMisalignment and Resonance Torques and Their
Treatment in the GP-B Data Analysisâ (G. M. Keiser et al. (2009)) details these two unfolding surprises.
Briefly, the story is as follows. In pointing the spacecraft off to a series of real and virtual stars at much
larger angles to the gyro spins during the Calibration Phase, we observed much larger than expected gyro
drift rates (Figure 9). The rate proved to be proportional to the mean misalignment when <1°, and always
in the direction perpendicular to this misalignment. The data provided a useful initial calibration of these
torques, but not one sufficiently accurate for the final analysis. In Section 3, we describe the methods by
which more advanced calibrations are achieved.
An interesting aspect of GP-B is that whereas the guide star is only visible for about half the orbit, gyro-
to-gyro comparison signals can be processed continuously even during the GSI phase. The roll-polhode
resonance torque was discovered in April 2007 by J. Kolodziejczak through just such comparisons. From
time to time, the spin direction of one of the gyroscopes would shift over the course of a few days by
angles of 20 or more marc-s with no such effect in the other three.
14 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Figure 9.
Mean Drift Rate vs. Misalignment for Gyroscope 3 During Calibration Maneuvers
Further investigation revealed that these shifts took place at times when there was a high-order resonance
between the slowly changing polhode period and the 77.5s satellite roll period. As Keiser et al. show, the
gyro displacement broadly follows the Cornu spiral pattern of Figure 10.
Figure 10.
Expected Change in Orientation due to Resonance Torques. The units in this figure are dimensionless but
are equal in the North-South and West-East directions.
3 The Evolving Data Analysis Process
The geodetic effect, being in Einsteinâs theory 170 times larger, is naturally easier to separate from
Newtonian torques and other disturbances than the frame-dragging effect. As early as September 2005, it
was visible in the raw data, or more exactly, in the data after an initial processing in which the aberration
1000
2000
3000
4000
30
210
60
240
90
270
120
300
150
330
180
0
M
ea
n N
ort
h-S
ou
th
M
isa
lig
nm
en
t
Mean West-East Misalignment
0
500
1000
1500
2000
2500
3000
3500
4000
0
0.5
1
1.5
2
2.5
3
3.5
4
Gyro 3, Drift Rate Magnitude vs. Mean Misalignment
D
rif
t R
at
e M
ag
ni
tu
de
(a
rc
s
ec
/d
ay
)
Mean Misalignment (arc sec)
k = 2.5 arc sec/day/degree
Gyroscope 3, Mean Rate (mas/day) vs. Mean Misalignment (as)
-1
0
1
2
3
4
5
6
7
-3
-2
-1
0
1
2
3
Orientation of Gyroscope Spin Axis
N
or
th
-S
ou
th
O
rie
nt
at
io
n
West-East Orientation
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results 15
signals were used for approximate calibration, and then removed from the plotted curve (compare Figure
7).
Two complementary data processing methods, called for convenience âgeometricâ and âalgebraicâ,
emerged as the analysis proceeded. The algebraic method, partially set up before launch, depends on a
sophisticated filtering process with explicitly modeled torques. The geometric method, begun in August
2006, has the immediate benefit of illuminating why a clean separation of relativity from the
misalignment drift is possible. The relativity signal is in a fixed direction in inertial space; the
misalignment drift, depending as it does on spacecraft pointing, is varied over the course of a year in an
exactly known way by the annual aberration of starlight. Define two directions, varying in inertial space,
one
radial
, parallel to the misalignment, and the other
axial
, at right angles to it. The relativity signal can
be resolved into a radial term free of classical drift, and an axial one, where the two are superimposed.
The geometric separation of relativity from the misalignment drifts proved strikingly effective, and has
gone through a number of refinements. An essential feature in it and the corresponding âalgebraicâ
treatment is the need for a continuous history of the misalignment through both GSV and GSI periods.
During GSI, accurate pointing information is obtained from the science gyros themselves.
The misalignment torque is, as already noted, only one of the three principal data analysis issues, the
others being accurate calibration of the gyro scale factor
C
g
and treatment of the roll-polhode resonance
torques. Both depend on a detailed process of Trapped Flux Mapping, begun in November 2006 and now
complete.
Trapped fields
London field at 80 Hz: 57.2 ÎźG
Gyro 1 3.0 ÎźG
Gyro 2 1.3 ÎźG
Gyro 3 0.8 ÎźG
Gyro 4 0.2 ÎźG
Trapped Flux
Moment
M
L
M
T
Polhode
Path
Figure 11.
Ideal vs. Actual London Moment Readout.
Figure 11 illustrates the trapped flux in the rotor, and its effect on the gyro scale factor
C
g
. Superimposed
on the London moment is a flux ranging from 0.3% - 5%
M
L
depending on which gyro, with components
parallel and perpendicular to
M
L
. The perpendicular component provides a very accurate measurement of
the gyro spin speed. The parallel component adds to
M
L
, evolving with time as the polhode period varies.
The SQUID readout signal,
z
i
, is the dot product of the normal,
n
, to the gyro readout loop with the
difference between
Ď
, the satellite roll axis direction and
s
, the inertial orientation of the gyroscope,
multiplied by
C
g
. The
n
direction changes uniformly with the vehicle roll angle
Ď
r
. The angle
Ď
,
measured to high accuracy by the cryogenic telescope, includes both the annual (20.4958 arc-s peak) and
orbital (5.1856 arc-s peak) aberrations computed respectively from JPL ephemerides data and precise
orbit measurements from onboard GPS receivers and laser reflectors. Since the measurement noise,
ν
SQUID
, near roll frequency is ~200 marc-s/Hz½ and
C
g
needs to be known to ~1 part in 10
5
, data from
16 Gravity Probe B Science ResultsâNASA Final Report, December 2008
many successive half-orbits must be connected unambiguously. To do so requires a continuous accurate
determination of the polhode phase
Ď
P
.
Though the readout was not originally intended to measure
Ď
P
, data in the form of 2s long 2.2 kHz bursts
was available from engineering telemetry channels. From it we obtained âsnapshotsâ of the trapped flux
signal over a number of gyroscope spin cycles, and from the evolving flux pattern thus seen were able to
track the motion of the spin axis in the body frame over time with remarkable precision.
Figure 12
. Gyro 1 Trapped Flux Map (left) and resulting
C
g
variations (right).
This Trapped Flux Mapping (TFM) technique produces two critical results: 1) a history of the polhode
phase good to 1° over the entire 353 days of science data, and 2) a direct estimate of polhode-driven
variations in the
C
g
. The mapping is accomplished by fitting a constant-coefficient rotor-fixed flux map
to a dynamic model of the polhode motion. Figure 12 shows (left) the flux map for Gyro 1 which
remained perfectly constant within the measurement limit throughout the mission, and (right) the periodic
variation in scale factor over a 4-hour time period as the polhoding modulates the direction of the flux
pattern. The results for Gyros 2, 3 and 4 had a similar form with successively lower flux levels. With
this, data from multiple orbits could be chained together in order to calibrate
C
g
with enough precision to
resolve marc-s changes in the gyro orientation. Absolute polhode phase is now known to ~1° over the
entire mission for all four gyroscopes, spin speeds to ~1 nHz and spin-down rates to ~1 pHz/sec. One of
the dramatic moments of GP-B came in August 2007 when the first reasonably complete TFM data was
inserted into the âalgebraicâ method, and scatter in the processed data abruptly dropped from 100
Ď
to 2
Ď
.
The roll-polhode resonance torques have already been described. Coming to one gyro at a time, the
resonances were in general unambiguously identifiable. Our first approach was to exclude the resonant
period data. A second, slightly more advanced, was to bridge across with a simple ramp function. With
the recognition that the dynamics follow or nearly follow a Cornu spiral, a more complete treatment
became possible by incorporating an additional term into the basic data processing model with
continuously varying functions of the polhode phase and angle. The result was a remarkable reduction in
spread between the processed data from the four gyroscopes. Figure 13 shows the results from individual
gyroscopes, and the four gyros combined. The four separate error ellipses overlap. Details are in the
accompanying paper by Heifetz et al. (Heifetz et al. (2009)); the current limits on systematic errors are
reviewed in the accompanying paper by Muhlfelder et al. (Muhlfelder et al (2009)).
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results 17
Figure 13
. Preliminary Relativity Estimates. Gyroscopes 1, 2, 4: Segments 5, 6, and 9, (Gyro 3: Segments 5, 6), not including
systematic error or a model sensitivity analysis. Note that the âGR predictionâ in the Figure is the net expected value of the
attached table allowing for solar geodetic and guide star proper motion terms.
An important result of a different kind has emerged from continuing work on the geometric approach.
The guide star (IM Pegasi) is known to be a binary with period 24.6 days and amplitude 1-2 marc-s. The
motion has been determined with considerable accuracy by our colleagues at the Smithsonian
Astrophysical Observatory and York University, Canada (I. I. Shapiro et al. (2007)). Careful
investigation has identified an ellipse of the same period and approximate amplitude in the processed GP-
B data. When data processing has reached a more advanced stage, the two independent measurements
will be made the basis of a suitable blind comparison test.
The results of Figure 13 were limited by available computational speed. To bring the processing within
the available constraints, we had to average the data once per orbit. We are now actively pursuing high
speed computer techniques which will overcome this limitation and allow 2 second processing to more
completely investigate the details of the roll-polhode resonance torque.
4 Conclusion
The underlying physics of the major disturbances to the science measurement is understood. Methods to
remove them credibly from the measurement have a solid physics foundation and result in greatly
improved quality of fit to the data model. Precise knowledge of polhode phase, now known to 1 deg
throughout the mission, is the key synchronizing element needed to tie together the entire data set. The
TFM machinery can accurately model polhode-modulated mission-length variations of the readout scale
factor, thus allowing separation of them from estimates of gyro drift.
Four Gyros Combined
NS: -6550 Âą 14.0 marc-s/yr
EW: -69.1 Âą 5.8 marc-s/yr
Earth
Solar
Geodetic
Proper
Motion
Net
Expected
NS
-6606
+7
+28 Âą 1
-6571 Âą 1
EW
-39
-16
-20 Âą 1
-75 Âą 1
18 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Analysis through the Fall of 2007 has provided 20% determinations of the ~ 1 marc-s orbital motion of
the guide star and of the bending of starlight due to the sun. The results in Section 3 for the 155 day
processing period from December 10, 2004 through May 27, 2005 (155 effective days of data) are
encouraging.
A non-zero WE gyroscope drift establishes the frame-dragging effect with a present
statistical uncertainty of 15%. The NS drift component provides better than a 0.5% measurement
of the geodetic effect.
The systematic uncertainties for the GR effects need to be added to these
statistical uncertainties using the methods outlined by Muhlfelder et al. The current results look very
encouraging, and with the extension to 279 days or longer science data and the use of more advanced
processing techniques, we hope to approach the 3% - 5% covariance limit.
Reaching the final limit will include completing the work on the geometric and algebraic science analyses
and performing an even more extensive investigation of the systematics than is given in the
accompanying paper by Muhlfelder et al. The use of high speed computing to decrease the analysis time
step from once per orbit (5859 s) to once every two seconds will significantly increase the fidelity of the
analysis by providing multiple data points per roll cycle for the modeling of the roll-polhode resonance
torque. No fundamental âbreakthroughsâ appear to be needed to extend the analysis to these other
segments. The necessary techniques are being developed in collaboration with Prof. Charbel Farhat of
the Stanford University Department of Aeronautics and Astronautics and his colleagues.
Acknowledgements
This work was supported by NASA Contract NAS8-39225 through the George C. Marshall Space Flight
Center (MSFC). From January 15, 2008 through September 30, 2008 critical bridge funding was
provided through a generous personal gift by Mr. Richard Fairbank with matching funds from Stanford
University and NASA. Acknowledgement is also made for significant funding from KACST in the
period since September 30, 2008.
References
Barker, B. M. and R. F. O'Connell (1970). Phys. Rev. D 2: 1428.
Breakwell, J. V. (1988). The Stanford Relativity Gyroscope Experiment (F): Correction to the Predicted
Geodetic Precession of the Gyroscope Resulting from the Earth's Oblateness. Near Zero: New
Frontiers of Physics. J. D. Fairband, J. B. S. Deaver, C. W. F. Everitt and P. F. Michelson. New
York, W. H. Freeman and Company: 685-690.
Cabrera, B. (1982). Near Zero Magnetic Fields with Superconducting Shields. Near Zero: New Frontiers
in Physics, Stanford, CA, W. H. Freeman and Co., New York.
Darling, T. W. (1989). Electric Fields on Metal Surfaces at Low Temperatures. School of Physics.
Parkville, Victoria 3051 Australia, University of Melbourne: 88.
Duhamel, T.G. (1984). "Contributions to the Error Analysis in the Relativity Gyroscope Experiment.
Dept. of Aeronautics and Astronautics and Dept. of Physics", Stanford University: 112 pages.
Everitt, C. W. F. (2007). "Results from Gravity Probe B." Bulletin of the American Physical Society
52(3).
Heifetz, M., W. Bencze et al. (2009). "The Gravity Probe B Data Analysis Approach", Space Science
Reviews (Chapter 4 of this report).
Keiser. G. M., J. Kolodziejczak, A. S. Silbergleit (2009). "Misalignment and Resonance Torques and
Their Treatement in the GP-B Data Analysis", Space Science Reviews (Chapter 3 of this report).
Lense, J. and H. Thirring (1918). "On the Influence of the Proper Rotations of Central Bodies on the
Motions of Planets and Moons According to Einstein's Theory of Gravitation." Zeitschrift fur
Physik 19: 156. English translation in Mashhoon, B., F. W. Hehl, et al. (1984). "On the
Everitt et al: GP-B Data AnalysisâStatus & Potential for Improved Accuracy of Scientific Results 19
Gravitational Effects of Rotating Masses: The Lense-Thirring Papers." General Relativity and
Gravitation 16: 711-750.
London, F. (1960). Superfluids, Macroscopic Theory of Superconductivity. New York, Dover
Publications.
Mester, J. C., J. M. Lockhart, et al. (2000). "Ultralow Magnetic Fields and Gravity Probe B Gyroscope
Readout." Advances in Space Research 25(6): 1185-8.
Muhlfelder, B., J. M. Lockhart, et al. (2003). "The Gravity Probe B Gyroscope Readout System."
Advances in Space Research 32(7): 1397-1400.
Muhlfelder, B. et al. (2009). "GP-B Systematic Error", Space Science Reviews (Chapter 5 of this report).
Parmley, R. T., G. A. Bell, et al. (2003). "Performance of the Relativity Mission Superfluid Helium Flight
Dewar." Advances in Space Research 32(7): 1407-1416.
Parmley, R. T., J. Goodman, et al. (1986). "Gravity Probe B Dewar/Probe Concept." Proceedings of SPIE
619: 126-133.
Pugh, G. E. (1959). Proposal for a Satellite Test of the Coriolis Prediction of General Relativity.
Washington, D. C., Weapons Systems Evaluation Group, The Pentagon.
Shapiro, I. I. et al. (2007). Bulletin of the American Physical Society (3).
Silbergleit, A. S., J. Conklin, et al. (2009). "Polhode Motion, Trapped Flux, and the GP-B Science Data
Analysis", Space Science Reviews (Chapter 2 of this report).
Schiff, L. I. (1960). "Possible New Experimental Test of General Relativity Theory." Phyical Review
Letters 4(5): 215-217.
Turneaure, J. P., E. A. Cornell, et al. (1982). The Stanford Relativity Gyroscope Experiment (D);
Ultrahigh Vacuum Techniques for the Experiment. Near Zero: New Frontiers of Physics,
Stanford, CA, W. H. Freeman and Co., New York.
20 Gravity Probe B Science ResultsâNASA Final Report, December 2008
21
2
P
olhoDe
m
otion
, t
raPPeD
f
lux
,
anD
the
GP-b s
cience
D
ata
a
nalysis
Alex Silbergleit et al.
22 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Silbergleit et al: Polhode Motion, Trapped Flux, and the GP-B Science Data Analysis 23
Noname manuscript No.
(will be inserted by the editor)
Polhode Motion, Trapped Flux, and the GP-B Science
Data Analysis
A. Silbergleit
¡
J. Conklin
¡
D. DeBra
¡
M. Dolphin
¡
G. Keiser
¡
J. Kozaczuk
¡
D. Santiago
¡
M. Salomon
¡
P. Worden
Received: / Accepted:
Abstract
Magnetic ďŹeld trapped in the Gravity Probe B (GP-B) gyroscope rotors
contributes to the scale factor of the science readout signal. This contribution is mod-
ulated by the rotorâs polhode motion. In orbit, polhode period was observed to change
due to a small energy dissipation, which signiďŹcantly complicates data analysis. We
present precise values of spin phase, spin down rate, polhode phase and angle, and
scale factor variations obtained from the data by Trapped Flux Mapping. This method
ďŹnds the (unique) trapped ďŹeld distribution and rotor motion by ďŹtting a theoreti-
cal model to the harmonics of high (gyroscope spin) frequency signal. The results are
crucial for accurately determining the gyroscope relativistic drift rate from the science
signal.
Keywords
Gravity Probe B
¡
Polhode motion
¡
Dissipation
¡
Trapped Flux Mapping
PACS
04
¡
04.20-q
¡
45.40Cc
¡
74.25-q
¡
02.70.-c
1 Introduction
The GP-B relativity science mission to measure the geodetic and frameâdragging eďŹects
on the four spinning superconducting gyroscopes [see Adler, Silbergleit (2000) and the
refrences therein] is described in Everitt et al. (1980); Turneaure et al. (2003); Everitt et
al. (2009). It uses the low frequency (LF) London Moment magnetic readout provided
by a highâprecision lowânoise SQUID. In addition, a high frequency (HF) signal from
residual magnetic ďŹeld trapped in the type II superconductor (niobium) is present,
as a superposition of multiple harmonics of the spin frequency. This paper deals with
the analysis of the HF signal, and presents the analysis results, which are crucial for
obtaining the desired accuracy of the measurements.
In sect. 2 the relation between the gyro polhode motion, trapped ďŹeld, and science
readout is explained, and the sources of HF signal are described. The on-orbit discovery
Silbergleit (corresponding author)
Gravity Probe B, HEPL, Stanford University, Stanford, CA 94305-4085, USA
Tel.: +1-(650)-723-1641
Fax: +1-(650)-725-8312
E-mail: gleit@stanford.edu
24 Gravity Probe B Science ResultsâNASA Final Report, December 2008
2
of the time variation of the gyro polhode period, and the explanation of this eďŹect as
kinetic energy dissipation, are given in sect. 3, along with some results obtained from
the analysis of the measured polhode period time history using a general dissipation
model. Sect. 4 explains the concept of Trapped Flux Mapping (TFM; determination of
trapped magnetic ďŹeld and characteristics of gyro motion from the HF signal analysis),
its relevance and importance for the ďŹnal result of GP-B experiment, and the key points
of the procedure. The actual approach to TFM in the form of a three level computer
analysis scheme is detailed in sect. 5, while the results of TFM are given in sect. 6.
Sect. 7 contains some conclusions.
2 Gyroscope Polhode Motion, Trapped Flux, and GP-B Readout
Polhode motion is important in the context of GP-B data analysis: the trapped mag-
netic ďŹux couples to the SQUID pick-up loop, which contributes to the overall readout
scale factor. Below we describe the model of this phenomenon.
2.1 Free Gyroscope Motion: Polhoding
Let us choose the principal inertia axes of a GP-B gyroscope, Ë
I
1
,
Ë
I
2
,
Ë
I
3
, as a Cartesian
bodyâďŹxed frame, and number the moments of inertia in their non-decreasing order,
0
< I
1
â¤
I
2
â¤
I
3
. The GP-B ďŹight rotors are the best spheres ever made, as certiďŹed by
the Guinness World Records,
1
, with
I
i
â
I
j
/I
k
âź
10
â
6
. Still, the value of the inertial
asymmetry parameter
, deďŹned as
0
â¤
Q
âĄ
I
2
â
I
1
I
3
â
I
1
â¤
1
,
(1)
can be anywhere between zero and one. The case
Q
= 0 corresponds to a rotor sym-
metric about the maximum inertia axis, Ë
I
3
.
The instanteneous position of the spin axis,
Ď
s
, in the rotor body can be described
by the angle,
Îł
p
, between it and the axis Ë
I
3
, and the second angle,
Ď
p
, between the
projection of
Ď
s
on the plane
{
Ë
I
1
,
Ë
I
2
}
and the axis Ë
I
1
(see Fig. 1). The motion of
the spin axis in the body,
Ď
p
(
t
)
, Îł
p
(
t
), obeys the Euler equations. If the torques are
negligible and the motion is free, it is described by the exact Euler solution [e.g. Landau,
Lifshitz (1959); MacMillan (1960)] corresponding to the precession of the spin axis
Ď
s
about the inertia axis ( Ë
I
3
, in the case of GP-B gyros). This precession is called the
polhode motion, or polhoding; we also call
Îł
p
and
Ď
p
the polhode angle and phase,
respectively. The period of the polhode motion,
T
p
=
const
, as well as the path of the
vector
Ď
s
in the body, is determined by the values of the conserved angular momentum,
L
=
const
, and energy,
E
=
const
. The polhode angular rate,
âŚ
p
âĄ
2
Ď/T
p
â
î¨î°
(
I
3
â
I
1
)(
I
3
â
I
2
)
/I
3
îŠ
Ď
s
â
(1
â
4)
Ă
10
â
6
Ď
s
(or smaller)
,
for GP-B gyros. It is worthwhile noting that in the case of a symmetric rotor (
Q
= 0),
the polhode path is a circular cone and the motion is uniform, i.e., with
Îł
p
=
const
and
Ë
Ď
p
=
âŚ
p
=
const
, so the polhode phase is a linear function of time. In contrast, for an
1
See http://einstein.stanford.edu/MISSION/mission6.html#awards
Silbergleit et al: Polhode Motion, Trapped Flux, and the GP-B Science Data Analysis 25
3
Fig. 1
BodyâďŹxed frame and the instanteneous position of the spin axis.
asymmetric rotor (
Q >
0) the polhode cone is
not
circular, and the motion of
Ď
s
in the
body is no longer uniform. In particular,
Îł
p
(
t
) oscillates with the period
T
p
(
t
)
/
2,
Ď
p
(
t
)
is not a linear function of time (modulated with the same period), and the angular
velocity Ë
Ď
p
(
t
)
îś
=
const
does not coincide identically with
âŚ
p
(
t
).
2.2 GP-B Readout: London Moment and Trapped Flux
The main source of the GP-B magnetic readout is the dipole ďŹeld of the London
Moment (LM) aligned with the gyroscope spin axis [London (1961), sect. 12]. However,
this is not the only source: when a type II superconductor is cooled below the critical
temperature, residual magnetic ďŹeld (even under the GP-B conditions!) is trapped
inside the rotor. Magnetic quanta
Âą
ÎŚ
0
form (macroscopically) point sources on the
rotor surface (ďŹuxons). The multi-pole ďŹeld of these point sources contributes to the
total ďŹux through the SQUID pick-up loop, which rolls with the spacecraft in inertial
space with the roll period
T
r
â
77
.
5
s
.
The SQUID signal is proportional to the total magnetic ďŹux through the pick-up
loop, which is the sum of the LM and trapped ďŹux (TF) contributions,
ÎŚ
(
t
) =
ÎŚ
LM
(
t
) +
ÎŚ
T F
(
t
)
.
(2)
It is straightforward to see that the LM ďŹux is proportional to the small angle
β
(
âź
10
â
4
or smaller) between the LM and pick-up loop plane,
ÎŚ
LM
(
t
) =
C
LM
g
β
(
t
). The LM scale
factor
C
LM
g
, proportional to the spin speed, is constant, up to a very small spinâdown
rate of GP-B gyros [spinâdown time
âź
10
4
yr
, see Everitt et al. (2009)]. The angle
β
(
t
)
carries the
relativity signal
at the low roll frequency
f
r
â
0
.
01
Hz
.
The time signature of the TF signal emerges from the following argument. Fluxons
are frozen in the rotor surface and spin with it; the function that converts a ďŹuxonâs
position into its magnetic ďŹux through the pick-up loop is strongly nonlinear (almost a
stepâfunction, due to a relatviely small gap between the rotor and the loop). Therefore
the TF signal consists of multiple harmonics of spin, and, since the pick-up loop is
rolling, those are actually the harmonics of spin
Âą
roll frequency. Finally, since the
spin axis moves in the rotor body (polhoding), the spin harmonics are modulated by
26 Gravity Probe B Science ResultsâNASA Final Report, December 2008
4
the polhode frequency, so that
ÎŚ
T F
(
t
) =
î
n
H
n
(
t
)
e
in
(
Ď
s
Âą
Ď
r
)
=
î
|
n
|
=
odd
H
n
(
t
)
e
in
(
Ď
s
Âą
Ď
r
)
+
β
(
t
)
î
|
n
|
=
even
h
n
(
t
)
e
in
(
Ď
s
Âą
Ď
r
)
,
(3)
where
H
n
(
t
) =
H
n
(
Ď
p
(
t
)
, Îł
p
(
t
))
, h
n
(
t
) =
h
n
(
Ď
p
(
t
)
, Îł
p
(
t
)). The accurate derivation of
the formula (3) is given in Keiser, Silbergleit (1995); Nemenmann, Silbergleit (1999),
and for the case of non-rolling spacecraft outlined in Keiser, Cabrera (1983), where an
alternative approach was used [see Sect. 5.1].
The fact that even harmonics of spin are multiplied by the (same as above) small
angle
β
is of purely geometric nature: the spin axis is set to reside almost in the pick-up
loop plane in the GP-B experiment.
According to Eq. (2), the LM ďŹux and d.c. part of the TF [
n
= 0 in the expression
(3)] combine to provide the GP-B
low frequency
(LF)
science readout
(sampled after
the additional lowpass ďŹlter with a 4
Hz
cutoďŹ) given by:
ÎŚ
LF
(
t
) =
ÎŚ
LM
(
t
) +
ÎŚ
T F
DC
(
t
) =
C
LM
g
β
(
t
) +
C
T F
g
(
t
)
β
(
t
)
,
C
T F
g
(
t
)
âĄ
h
0
(
t
)
.
(4)
Thus, polhode motion combined with the TF creates polhode variations of the readout
scale factor. These variations are not larger than 5% of the (constant) LM scale factor
for all the GP-B gyros.
2.3 GP-B High Frequency Data
There were two telemetry sources of the high frequency (HF) SQUID signal: 1) FFT
of ďŹrst six harmonics of spin, and 2)
snapshots
, i.e.,
âź
2
s
stretches of original SQUID
signal sampled at 2200
Hz
. Both were recieved only during the part of an orbit when
the Guide Star was occulted by the earth. At time when the snapshots were available,
a snapshot would come about every 40
s
, but the gaps between the snapshot arrays
were up to 2 days.
HF FFT harmonics were analyzed during the mission as soon as the telemetry was
in place. Snapshots, being the most accurate source of HF information, were thoroughly
treated after the mission. All 976,478 snapshots available from the mission science
period were processed after the ďŹight (see sect. 4).
3 Changing Polhode Period: OnâOrbit Discovery and Its Explanation
One of the two onâorbit discoveries which aďŹected GP-B data analysis most strongly
was the changing polhode period. The full time history of the polhode period for each
of the GP-B gyroscopes is shown in Fig. 2, with very good agreement of the results
obtained using the two HF signals. Moreover, an entirely diďŹerent gyro position signal
(polhode modulation of its the spin component) gives the same result [see Dolphin
(2007)]. Notably, the
T
p
(
t
) plot for gyros 1, 2 has a peak where formally
T
p
=
â
.
Nevertheless, throughout the science period (starting Sept. 13, 2004), polhode periods
of all the four gyros decrease monotonically and tend to speciďŹc asymptotic values
T
pa
.
Silbergleit et al: Polhode Motion, Trapped Flux, and the GP-B Science Data Analysis 27
5
Fig. 2
Polhode period history of GP-B gyros: redâfrom HF FFT, blueâfrom snapshots.
Such behavior can only be explained by
energy dissipation
, which reduces the kinetic
energy of the rotor without changing its angular momentum. This phenomenon was
observed, in particular, in rolling satellites starting with the âExplorerâ launched in
1958 [Modi (1973)]. Dissipation moves the spin axis in the body to the
maximum
inertia axis where the energy is easily seen to have a
minimum
, under the conserved
angular momentum constraint. Thus the polhode path, controlled by the parameter
L
2
/
2
E
, also changes (note that the spinâdown torque changes both quantities
L
and
E
, but
not
this ratio, hence
not the polhode path
). The two types of behavior seen in
the plots of Fig. 1 are also explained: by pure chance, gyros 1, 2 start their evolution
with the spin vector precessing about Ë
I
1
, so
T
p
tends to inďŹnity when
Ď
s
crosses
the separatrix [Landau, Lifshitz (1959); MacMillan (1960)]. In contrast with this, the
evolution of gyros 3, 4 starts when
Ď
s
is already precessing about Ë
I
3
, so
T
p
(
t
) just
decreases monotonically.
Writing
Ď
s
=
Ď
1
,
3
for the case when
Ď
s
is parallel to Ë
I
1
,
3
, respectively, and using
the angular momentum conservation [
I
1
Ď
1
=
I
3
Ď
3
, Ď
1
= (
I
3
/I
1
)
Ď
3
] one can readily
determine the relative energy loss (and spin speed reduction) from the minimum, Ë
I
1
,
to the maximum, Ë
I
3
, inertia axis:
E
1
â
E
3
E
1
=
I
1
Ď
2
1
â
I
3
Ď
2
3
I
1
Ď
2
1
=
Ď
1
â
Ď
3
Ď
1
=
I
3
â
I
1
I
3
<
4
Ă
10
â
6
(5)
for GP-B gyros. The total energy loss needed to move the spin axis all the way from
min to max inertia axis is thus less than 4
ÂľJ
(
E
âź
1
J
); in one year, the average
dissipated power required for this is just 10
â
13
W
!
The physical origin of this energy dissipation is not completely clear. It might be
related to inelastic deformations of the rotor, but probably is dominated by dissipative
patch eďŹect torques [see Buchman, Gill (2007)]. Independent of the origin, we have
found a
general dissipation model
in the form of an additional term in the Euler equa-
tions unique up to a scalar factor. The detailed description of these new equations
along with their solution will be published separately [a brief derivation of the model
can be found in Salomon (2008)]. Fitting the polhode period time history obtained
from the model to the measurements allowed us to determine the rotor asymmetry
parameter (also found, in some cases, from the gyro position signal), the asymptotic
period
T
pa
âź
1
â
2 hours, and the characteristic time of dissipation
Ď
dis
âź
1
â
2 months,
depending on the gyro (see Table 1). Thus the dissipation is
slow
(
T
p
î
Ď
dis
), so that
the polhode motion of the GP-B gyros is
quasi-adiabatic
.
28 Gravity Probe B Science ResultsâNASA Final Report, December 2008
6
Table 1
Asymptotic Polhode Period,
T
pa
, and Dissipation Time,
Ď
dis
Parameter
Gyro 1
Gyro 2
Gyro 3
Gyro 4
T
pa
(hours)
0.867
2.851
1.529
4.137
Ď
dis
(days)
31.9
74.6
30.7
61.2
4 Trapped Flux Mapping: Concept and Importance
4.1 Trapped Flux Mapping
Trapped Flux Mapping
(TFM) is the procedure of ďŹnding the distribution of trapped
magnetic ďŹeld and charcteristics of gyro motion from
odd
harmonics of HF SQUID
signal, by ďŹtting them to their theoretical model. The latter is based on the follow-
ing solution for the scalar magnetostatic potential outside the spherical rotor with
K
pairs of ďŹuxons and anti-ďŹuxons on its surface (in the bodyâďŹxed frame described in
sect. 2.1):
Ψ
(
r, θ, Ď
) =
ÎŚ
0
2
r
g
â
î
l
=1
î
r
g
r
î
l
+1
1
l
+ 1
l
î
m
=
â
l
A
lm
Y
lm
(
θ, Ď
)
,
A
lm
=
K
î
k
=1
î
Y
â
lm
(
θ
+
k
, Ď
+
k
)
â
Y
â
lm
(
θ
â
k
, Ď
â
k
)
î
.
(6)
Here
ÎŚ
0
is the magnetic ďŹux quantum, and
r
g
is the rotor radius (more on this model
in sect. 5.1).
If
the ďŹuxon number,
K
, and positions,
θ
Âą
k
, Ď
Âą
k
, were known,
then
the
coeďŹcients
A
lm
would be found uniquely (and easily!) by the above formula. However,
neither are known, in reality, so the coeďŹcients
A
lm
are to be estimated
by the TFM
procedure, along with the main functions describing the gyro motion, such as the spin
phase, the polhode phase, etc.. As shown below, TFM provides the asymmetry param-
eter
Q
and all coeďŹcients
A
lm
up to
l
max
= 21
,
25
,
21
,
21 (gyros 1â4 respectively)
for each rotor. Likewise, for the entire duration of the mission, TFM gives the spin
speeds to 10
nHz
; the spin down rates to 1
pHz/s
; the spin phases to 0
.
05
rad
; the
polhode phases to 0
.
02
rad
, and the polhode angles to 0
.
01
â
0
.
1
rad
. With these
results we compute the complete scale factor variations
C
T F
g
(
t
) for the whole mission
using formula (9).
The TFM products play a crucial role in the GP-B science analysis, which simply
cannot be accomplished without them with the needed accuracy. First of all, accurate
polhode phase and angle values at any time of the mission are required for modeling
both the scale factor variations and patch eďŹect torques [ Heifetz et al. (2009); Keiser et
al. (2009)], because all the torque coeďŹcients are modulated by the polhode harmonics
in the same fashion as the LF SQUID scale factor. In addition, TFM brings in an
independent determination of those scale factor variations, which, ďŹrst, allows for a
separate determination of the LM scale factor and slowly varying D.C. part of the TF
scale factor. Second,
C
T F
g
(
t
) found by TFM can be directly used in the LF science
analysis thus dramatically reducing the number of estimated parameters, and even
turning the (originally nonâlinear) estimation problem into a linear one.
Silbergleit et al: Polhode Motion, Trapped Flux, and the GP-B Science Data Analysis 29
7
4.2 Key Points of TFM
To carry out TFM, one ďŹrst of all needs to process the measured HF SQUID signal,
z
HF
(
t
)
â
ÎŚ
HF
(
t
), [those almost one million snapshots mentioned in sect. 2.3], to
extract the multiple harmonics of spin,
H
n
(
t
), according to the formula (3):
measured
z
HF
(
t
) =
î
n
H
n
(
t
)
e
in
(
Ď
s
Âą
Ď
r
)
.
(7)
Then, starting with expression (6) for the magnetic potential in the rotorâďŹxed frame,
one transforms it to the frame of the pickâup loop using the rotation matrices for
spherical harmonics [Rose (1995)], and computes the magnetic ďŹeld. The ďŹux through
the loop can then be found by integrating the normal ďŹeld component [for details of
this derivation see Kozaczuk (2007); Conklin (2009)]. As in Eq. (3), the result is a sum
of the spin harmonics, but with an explicit expression for each of them. In particular,
the odd harmonics relevant to TFM are given by
H
n
(
t
) =
ÎŚ
0
2
â
î
l
=
|
n
|
l
odd
î
r
g
b
î
l
l
î
m
=
â
l
A
lm
d
l
0
n
(
Ď/
2)
d
l
mn
(
Îł
p
)
e
imĎ
p
I
l
,
n
odd
.
(8)
Here
b
is the pick-up loop radius,
r
g
/b
âź
1,
Ď
p
(
t
) and
Îł
p
(
t
) are the polhode phase
and angle deďŹned in sect. 2.1 [and roughly known from the polhode period modelling,
sect. 3],
d
l
mn
(
Îą
) is the rotation matrix, essentially, a hypergeometric polynomial of the
argument tan
2
(
Îą/
2) [see Rose (1995), formula (4.14)], and
I
l
is the known coeďŹcient.
Formula (8) is exactly the expression to ďŹt to the measured harmonics
H
n
(
t
). The
Euler angles
Ď
p
(
t
),
Îł
p
(
t
), and
Ď
s
(
t
), which enter it in a
nonlinear
way, are originally
derived from measurements, but their accuracy is insuďŹcient for properly determining
A
lm
. Expressions for these functions based on the dissipation model and Euler solution
are used in a
nonlinear
ďŹtting process described in sect. 5, which was developed to ďŹnd
them. Then a
linear
ďŹt is carried out to estimate the coeďŹcients
A
lm
and to complete
TFM.
Concluding this section we note that the same derivation that yields formula (8)
for odd harmonics provides similar expressions for even ones as well. The latter are
proportional to the spin-to-pickâup loop misalignment,
β
, to lowest order [see Eq. (3)].
According to Eq. (4), the D.C. harmonics (
n
= 0) gives the scale factor of the trapped
ďŹux contribution to the science signal, namely:
C
T F
g
(
t
) =
h
0
(
t
) =
ÎŚ
0
2
â
î
l
=1
l
odd
î
r
g
b
î
l
l
î
m
=
â
l
A
lm
d
l
m
0
(
Îł
p
)
e
imĎ
p
I
l
l P
l
â
1
(0)
(9)
Therefore, as soon as accurate values of
Ď
p
(
t
)
, Îł
p
(
t
) and
A
lm
are determined from
TFM, the scale factor variations can be immediately computed.
5 Trapped Flux Mapping: Three Analysis Levels
The TFM estimation process is broken up into three parts called
levels
, in which
various groups of parameters are estimated somewhat independently. Level A and B
30 Gravity Probe B Science ResultsâNASA Final Report, December 2008
8
data processing leads to the estimates of the polhode phase and angle,
Ď
p
(
t
) and
Îł
(
t
),
the spin phase,
Ď
s
(
t
), and the asymmetry parameter,
Q
. Once we know the rotor
orientation (these three Euler angles) and its asymmetry, in Level C we perform a
linear ďŹt to ďŹnd the coeďŹcients
A
lm
.
5.1 Level A: Asymmetry Parameter, Polhode Phase and Angle
The polhode phase,
Ď
p
(
t
), can be written as the sum of the
symmetric polhode phase
,
Ď
p
(
t, Q
= 0), and the polhode modulation,
âĎ
p
(
t, Q
), which has a period
T
p
/
2, and
is due to the inertial asymmetry described by the parameter
Q
from Eq. (1). The
symmetric polhode phase is the integral of the polhode rate (we choose
t
0
to be any
moment of time when
Ď
s
lies in the
{
Ë
I
1
,
Ë
I
3
}
plane),
Ď
p
(
t, Q
= 0) =
î
t
t
0
âŚ
p
(
t
î°
)
dt
î°
=
î
t
t
0
2
Ď
T
p
(
t
î°
)
dt
î°
,
(10)
and it is just a straight line if no dissipation is present,
T
p
=
const, âŚ
p
=
const
. An
expression for the deviation from the straight line (a sum of exponents with the time
scale
Ď
dis
) is implied by the dissipation model, so the total model for the polhode phase
becomes:
Ď
p
(
t
) =
Ď
p
0
+
âŚ
pa
(
t
â
t
ref
)
â
M
î
m
=1
D
m
e
â
n
t
â
t
ref
Ď
dis
+
âĎ
p
(
t, Q
)
(11)
Here
Ď
p
0
is the value of the symmetric polhode phase at time
t
ref
with no dissi-
pation present,
âŚ
pa
= 2
Ď/T
pa
is the asymptotic polhode rate, and
D
m
are con-
stant coeďŹcients. Note that the total number of the parameters in the formula (11),
{
Ď
p
0
, âŚ
pa
, Q, D
1
, D
2
, . . . D
M
}
, that need to be estimated is 3 +
M
, with
M <
12
suďŹcient for a good enough accuracy for all the gyros. The polhode angle,
Îł
p
, is com-
puted from the Euler solution given the polhode rate, the spin rate
Ď
s
(
t
) (discussed in
the next section), and
Q
.
Estimating
Q
is the most computationally intensive part of the TFM data pro-
cessing. This is because
Q
estimation is sensitive to all ďŹt parameters, even though
Q
is entirely independent physically. As a result, only three separate batches of data
spanning 1-2 days each are used to estimate
Q
for every gyroscope. The batches were
chosen to be relatively long and early in the mission, when the polhode variations are
largest and the inďŹuence of
Q
is strongest. The Level A code is run many times on a
single batch of data while incrementing
Q
from 0 to 1 by 0.02 at the each next run.
The value of
Q
for which the RMS of the post-ďŹt residuals is minimum is taken as the
best estimate of
Q
, with its associated error. Table 2 shows the resulting values of
Q
from Level A processing. The accuracy of these results is
âź
20%, enough for estimation
of the scale factor variations due to trapped ďŹux to
âź
1 %. This is due to the relative
insensitivity of the signal to the asymmetry parameter.
Once
Q
is estimated, individual ďŹts are performed on one batch of data at a time
depending on how the data are grouped. Spectral analysis of the data provides initial
estimates of the polhode rate
âŚ
p
(
t
) accurate to
âź
100
nHz
[Salomon (2008)]. Integrat-
ing this preliminary value of
âŚ
p
(
t
) is suďŹcient for the Level A ďŹtting process, but the
initial polhode phase,
Ď
p
0
, is unknown. Therefore the initial polhode phase is estimated
for each batch of data using a standard Nelder-Mead simplex method [see Lagarias et
Silbergleit et al: Polhode Motion, Trapped Flux, and the GP-B Science Data Analysis 31
9
Table 2
Asymmerty Parameter,
Q
Gyro 1
Gyro 2
Gyro 3
Gyro 4
0.303
Âą
0.069
0.143
Âą
0.029
0.127
Âą
0.072
0.190
Âą
0.048
al. (1998)] implemented by the function
fminsearch.m
in the MATLAB software. This
gives the estimates of
Ď
p
0
at roughly 100 diferent times throughout the science mission.
Each estimate of the initial polhode phase has some error associated with it. To im-
prove the overall accuracy of the polhode phase determination, a single model, Eq. (11),
is ďŹt to all the preliminary polhode phase estimates. Not only does this smooth the
estimates, reducing the overall error, but it also produces a polhode phase time-history
that is continuous and self-consistent over the entire mission. The unknown parameters
Ď
p
0
,
âŚ
pa
and
D
m
all appear linearly in this model, allowing for a simple least squares
ďŹt. An estimate of the polhode phase during each batch of data is constructed from
the estimates of
Ď
p
0
for the same batch, and the integration of the estimated polhode
rate. To construct a self-consistent polhode phase, the estimated polhode rate accurate
to
âź
600
nrad/sec
is adequate for resolving
Ď
ambiguities between batches of data
that are typically not more that 4 days apart. The RMS of the polhode phase post-ďŹt
residuals for each gyroscope is on the order of 0.1 rad.
5.2 Level B: Spin Phase
The model used to describe the spin phase is based on the observed behavior of the
measured spin frequency, which is accurately described over short stretches by a straight
line plus a small contribution due to the changing polhode frequency. This results in a
spin phase model consisting of a quadratic polynomial plus two corrections,
Ď
s
(
t
) =
Ď
si
+
Ď
si
(
t
â
t
i
) +
1
2
dĎ
s
(
t
â
t
i
)
2
+
Ď
p
(
t, Q
= 0) +
âĎ
s
(
t, Q
)
(12)
The ďŹrst correction, the symmetric polhode phase
Ď
p
(
t, Q
= 0), is added because the
measured HF harmonics prove to be of the spin plus polhode frequency rather than
the spin frequency. The last term,
âĎ
s
(
t, Q
), is a small modulation of the spin phase
caused by the precession of the angular velocity,
Ď
s
, about the angular momentum
vector in inertial space. This modulation term can be explicitly computed from the
solution to Eulerâs equations transformed into the inertial frame [see Landau, Lifshitz
(1959); Salomon (2008)], so it requires no additional parameters to estimate. The model
(12) is independently ďŹt to every batch of data about 1 day duration. Since at this
step the parameters in the symmetric polhode phase model are considered known, the
three parameters to estimate for each batch are
Ď
si
,
dĎ
s
, and
Ď
si
.
The nonlinear search routine for these three parameters is a modiďŹed version of
the Nelder-Mead simplex method developed by Kozaczuk (2007) and implemented in
the properly modiďŹed version of the MATLAB function
fminsearch.m
. The nonlinear
search for the initial spin rate
Ď
si
and the decay rate
dĎ
s
is performed separately from
the search for the initial spin phase
Ď
s
0
, in order to simplify the nonlinear estimation
process. The Level B code iterates between a search for (
Ď
si
, dĎ
s
) and a search for
Ď
s
0
.
A total of four iterations is a good compromise between accuracy and computation
time for all the four GP-B gyros.
32 Gravity Probe B Science ResultsâNASA Final Report, December 2008
10
5.3 Level C: Trapped Field Distribution and Polhode Phase ReďŹning
With the parameters estimated in the Level A and Level B data processing, the time
histories of the three Euler angles,
Ď
p
(
t
),
θ
p
(
t
) and
Ď
s
(
t
), are computed for the whole
science period. Then a linear least squares ďŹt of the model (6) to the measured odd
harmonics of spin,
H
2
k
+1
(
t
), is performed to estimate the coeďŹcients
A
lm
using the
data from the entire mission. Fig. 3 shows the best-ďŹt values of the real and imaginary
parts of
A
lm
with the various maximum number of odd harmonics,
N
max
, and coef-
ďŹcients,
L
max
, used,
N
max
=
L
max
= 19
,
23
,
29. A good agreement in the common
values from diďŹerent ďŹts is seen, with a rather regular distribution, which turns out to
be normal with zero mean.
After this step, estimates exist for all necessary parameters in the model (8). An
iterative approach is used to reďŹne these estimates further.
One of the most fundamental quantities in the determination of the trapped mag-
netic potential and the prediction of the scale factor variations
C
T F
g
(
t
) is the polhode
phase,
Ď
p
(
t
). The estimation of
Ď
p
(
t
) from the Level A data processing is typically
accurate to 50
â
100
mrad
, based on the post-ďŹt residuals. A more accurate technique
for determining the polhode phase error comes in at Level C and involves the absolute
value of the measured ďŹrst harmonics,
|
H
1
(
t
)
|
. Using this quantity has two signiďŹ-
cant advantages for polhode phase determination. The ďŹrst is that
H
1
is the largest
signal present in the HF SQUID output, providing conďŹdence in its accuracy. The
second, more important advantage is that the absolute value of any measured harmon-
ics,
|
H
n
(
t
)
|
, is totally independent of errors in the spin phase, which is apparent from
Eq. (7). In fact,
|
H
1
(
t
)
|
eďŹectively depends only on the polhode phase and
A
lm
, since
the polhode angle is simply computed from the measured parameters.
Upon completion of Levels A, B and C, a best-ďŹt
|
H
1
(
t
)
|
can be constructed from
the estimated parameters and compared with the measured value. This comparison
gives a correction,
δĎ
p
, to the best-ďŹt polhode phase from Level A. The exponential
polhode phase model, Eq.(11), is then re-ďŹt to the corrected polhode phase. With the
new polhode phase, Level B processing is carried out again to reďŹne the spin phase
parameters, and then Level C is rerun to reďŹne the coeďŹcients
A
lm
. The polhode
phase reďŹnement, the Level B processing, and then the Level C processing constitutes
a single iteration of the process that converges to the optimal solution. We continue
these iterations until the RMS of the Level C post-ďŹt residuals reaches a stable value.
6 Trapped Flux Mapping: Results
The accuracy achieved in the Level C ďŹts and trapped ďŹux scale factor,
C
T F
g
(
t
), is
given in Table 3. The postâďŹt residuals are roughly a few percent for all gyroscopes,
which translates into
âź
10
â
4
error in the predicted trapped ďŹux scale factor relative
to the total scale factor
C
g
=
C
LM
g
+
C
T F
g
. The highest uncertainty in the estimated
C
T F
g
is for gyroscope 3, even though the post-ďŹt residuals are larger for gyroscope 4.
This is because the polhode damps out relatively quickly for gyroscope 3, causing poor
observability of the coeďŹcients
A
lm
. The calculated distribution of trapped magnetic
potential of the gyro 1 surface is shown in Fig. 4, along with the polhode paths and
scale factor variations,
C
T F
g
, on diďŹerent dates separated by more than 4 months. Due
to the dissipation, the polhode path shrinks over time, and the magnitude of scale
factor oscillations drops accordingly.
Silbergleit et al: Polhode Motion, Trapped Flux, and the GP-B Science Data Analysis 33
11
Fig. 3
Best-ďŹt Re (
A
lm
)
,
Im (
A
lm
), with error bars, for gyroscope 1 using up to
L
max
odd
harmonics,
L
max
= 19
,
23
,
29 - in black, blue, and red, respectively.
Table 3
Trapped Flux Mapping Results of Level C Fits
Gyro
Relative
Number of
Relative Size
Formal Error in
C
T F
g
No.
Residuals
Harmonics
of Variations
Relative to
C
g
1
1.1%
21
0.2% - 3.0%
1
.
5
Ă
10
â
4
â
7
.
0
Ă
10
â
5
2
1.5%
25
0.5% - 1.5%
6
.
0
Ă
10
â
5
â
3
.
0
Ă
10
â
5
3
2.6%
21
0.01% - 1.0%
2
.
0
Ă
10
â
4
â
1
.
6
Ă
10
â
4
4
2.8%
21
0.1% - 0.3%
8
.
5
Ă
10
â
5
â
6
.
5
Ă
10
â
5
Fig. 5 shows polhode variations of the scale factor,
C
g
(
t
)
T F
, for all four gyroscopes
during a few hours of the same day, 10 November 2004. An independent estimate of
the total gyroscope readout scale factor,
C
g
(
t
), comes from the analysis of the LF
SQUID signal. This determination is limited by the constantly changing trapped ďŹux
contribution, and it is not possible to separate the London Moment scale factor from
the d.c. part of the Trapped Flux. Nevertheless, this determination provides a useful
comparison of TFM. The LF results agree with the TFM results to
âź
10
â
3
or better,
relative to the total scale factor,
C
g
. The best agreement is for gyroscope 3, which has
the highest uncertainty in the
C
T F
g
estimated by TFM: for gyroscope 3 the relative
polhode variations,
C
T F
g
/C
g
, are by far smaller than for any other gyro.
7 Conclusions
The change of polhode period and path discovered on orbit is explained by a slow
rotation energy loss while conserving the angular momentum, and is properly analyzed
using a new general model of dissipative gyroscope motion. This lays the ground for
the developed procedure of mapping magnetic trapped ďŹeld. Trapped Flux Mapping
has been carried out successfully for each of the four GP-B gyroscopes using the odd
34 Gravity Probe B Science ResultsâNASA Final Report, December 2008
12
Fig. 4
Gyroscope 1: magnetic potential distribution, and polhode paths and scale factor
variations,
C
T F
g
(
t
), on Oct. 4, 2004, and Feb. 20, 2005. The path shrinks due to dissipation,
and the magnitude of variations goes down accordingly.
Fig. 5
Scale factor variations,
C
T F
g
(
t
), relative to the total scale factor,
C
g
(
t
), on Nov. 10,
2004, as estimated by both TFM (green) and the LF analysis (red).
Silbergleit et al: Polhode Motion, Trapped Flux, and the GP-B Science Data Analysis 35
13
harmonics of the HF SQUID signal. Its results are necessary to determine the LF
scale factor variations and modulation of patch eďŹect torque coeďŹcients, thus they are
crucial for the best measurement of relativistic drift rate.
Acknowledgements
This work was supported by NASA contract NAS 8-39225, by dona-
tions from Richard Fairbank and Stanford University, and recently through a collaborative
agreement with KACST. The authors are extremely grateful to W. Bencze, J. Berberian, S.
Buchman, B. Clarke, F. Everitt, M. Heifetz, T. Holmes, J. Li, B. Muhlfelder, I. Nemenman,
V. Solomonik and J. Turneaure for their numerous insights and help in completing this work.
References
Adler R.J., A.S. Silbergleit, Internat. J. of Theor. Phys.
39
(5), 1287 (2000).
Buchman S., D.K. Gill,
Evidence for Patch EďŹect Forces on the Gravity Probe B Gyroscopes
,
APS Meeting, April 2007. Abstract: Bulletin of the APS,
52
(3) (2007).
Conklin J.W.,
Estimation of the Mass Center and Dynamics of a Spherical Test Mass for
Gravitational Reference Sensors
, Ph.D. Thesis, Aero/Astro, Stanford University, (Stan-
ford, 2009).
Dolphin M.,
Polhode Dynamics and Gyroscope Asymmetry Analysis on Gravity Probe B Using
Gyroscope Position Data
, Ph.D. Thesis, Aero/Astro, Stanford University, (Stanford, 2007).
Everitt C.W.F., W.W. Fairbank, D.B. DeBra, et al.,
Report on a Program to Develop a
Gyro Test of General Relativity in a Satellite and Associated Control Technlogy
, HEPL,
Aero/Astro, Stanford University, (Stanford, 1980).
Everitt C.W.F., M. Adams, W. Bencze, et al., Space Science Reviews, this issue (2009) (Chap-
ter 1 of this report).
Heifetz M., W. Bencze, T. Holmes, et al., Space Science Reviews, this issue (2009) (Chapter
4 of this report).
Keiser G.M., B.Cabrera, in
Proc.of The National Aerospace Meeting
(The Institute of Navi-
gation, Washington, D.C., 1983).
Keiser G.M., J. Kolodziejczak, A.S. Silbergleit, Space Science Reviews, this issue (2009)
(Chapter 3 of this report).
Keiser G.M., A.S. Silbergleit,
Pick-up Loop Symmetry and Centering
, Gravity Probe B docu-
ment S0243 (Stanford University, Stanford, 1991).
Kozaczuk J.A.,
Precise Determination of the Spin Speed and Spin Down Rate of Gravity Probe
B Gyroscopes
, Physics Honor Thesis, Stanford University, (Stanford, 2007).
Lagarias J., et al., SIAM J. on Optimization
9
(1), 112 (1998).
Landau L.D., E.M. Lifshitz,
Mechanics
(Pergamon Press, Oxford, 1959).
London F.,
SuperďŹuids, v. 1
(Dover Publications, New York, 1961).
MacMillan W.D.,
Dynamics of Rigid Bodies
(Dover Publications, New York, 1960).
Modi V.J. , J. of Spacecraft and Rockets
11
, 743 (1973).
Nemenmann I.M., A.S. Silbergleit, J. of Appl. Phys.
86
(1), 614 (1999).
Rose M.E.,
Elementary Theory of Angular Momentum
, (Dover Publications, New York, 1995).
Salomon M.,
Properties of Gravity Probe B Gyroscopes Obtained from High Frequency SQUID
Signal
, Ph.D. Thesis, Aero/Astro, Stanford University, (Stanford, 2008).
Turneaure J.P., C.W.F. Everitt, B.W. Parkinson et al., Advances in Space Research
32
(7),
1387 (2003).
36 Gravity Probe B Science ResultsâNASA Final Report, December 2008
37
3
m
isaliGnment
anD
r
esonance
t
orques
anD
t
heir
t
reatment
in
the
GP-b D
ata
a
nalysis
George (Mac) Keiser et al.
38 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Keiser et al: Misalignment and Resonance Torques and Their Treatment in the GP-B Data Analysis 39
Misalignment and Resonance Torques and Their
Treatment in the GP-B Data Analysis
G. M. Keiser
1
, J. Kolodziejczak
2
, A. S. Silbergleit
Abstract
Classical torques acting on the GP-B gyroscopes
decrease the accuracy in the measurement of the
relativistic drift rate. Based on measurements made during the year-long science data collection, there are
two dominant classical torques acting on the gyroscopes. The first torque, known as the misalignment
torque, has a magnitude proportional to the misalignment between the gyroscope spin axis and the
satellite roll axis and is aligned perpendicular to the plane containing these two vectors. The second
torque, known as the resonance torque, mainly produces a permanent offset in the orientation of the
gyroscope spin axis when a harmonic of the gyroscope polhode frequency is in the vicinity of the satellite
roll frequency. These two torques have the same physical origin: an electrostatic interaction between the
patch effect fields on the surfaces of the rotor and the housing. In the post-mission data analysis, the
change in the gyroscope orientation due to both of these torques can be clearly separated from the
relativistic drift rate.
1 Introduction
As discussed in the accompanying papers (Everitt et al. (2009), Heifetz et al. (2009), Mulhfelder et al.
(2009), Silbergleit et al. (2009)), the Gravity Probe B (GP-B) satellite was designed to accurately measure
the geodetic and frame-dragging precessions predicted by the general theory of relativity (de Sitter
(1916), Lense and Thirring (1918)). Averaged over the orbital motion of the satellite (Adler and
Silbergleit (2000)), both of these predicted effects produce a linear drift in the orientation of the
gyroscope spin axis with time. A classical torque acting on a gyroscope can produce an error in an
accurate measurement of the relativistic drift rate.
This paper describes the two classical torques acting on the GP-B gyroscopes that had a measurable effect
on the orientation of the gyroscope spin axis. The misalignment torque is described in section 2 of this
paper, while the resonance torque is described in section 3 of this paper. In each of these sections, we
present the observational evidence for these torques, discuss the physical origin of each torque, and
demonstrate post-mission data analysis methods of separating the gyroscope drift rate caused by each
torque from the relativistic drift rate. Both torques are caused by an interaction of the patch effect
potential on the surface of the gyroscope rotor with a patch effect potential on the surface of the
gyroscope housing. Additional observations that can be explained by the patch effect will be described in
a forthcoming paper (Buchman et al. (2009)). Because the misalignment torque acts in a specific direction
and because the resonance torque produces a unique time signature in the orientation of the gyroscope
spin axis, the effects of these torques on the gyroscope spin axis orientation can clearly be separated from
the relativistic drift rate. The final section of this paper discusses the impact of the torques on the
measurement of the relativistic drift rate.
1
Keiser (corresponding author)
Gravity Probe B, HEPL, Stanford University, Stanford, CA 94305-4085 USA
Tel: 1-650-725-4116
FAX: 1-650-725-8312
E-mail: mac@relgyro.stanford.edu
2
NASA/Marshall Space Flight Center
40 Gravity Probe B Science ResultsâNASA Final Report, December 2008
2 Misalignment Torques
2.1
Observations
The Gravity Probe B mission timeline was divided into three phases. The initialization phase, which
began after launch on April 20, 2004, through the spin-up and final alignment of the gyroscopes on
August 29, 2004, included the initial calibration of the instrument, tests of each of the four gyroscopes at
slow spin speeds, adjustment of the spacecraftâs attitude and translation control system, spin-up of each of
the four gyroscopes, and a final alignment of the spin axes. The main science data collection phase began
on August 29, 2004, and continued until August 15, 2005. The final phase of the mission was called the
calibration phase, where potential disturbances were deliberately increased to determine their impact on
the measurement of the spin axis orientation with respect to the guide star and to measure other possible
classical torques acting on the gyroscope.
During the calibration phase, 17 spacecraft operations were performed to study the effect of increasing the
angle between the satellite roll axis and the four gyroscope spin axes. Throughout the science data
collection phase, while the guide star, IM Pegasi (HR 8703), was not occulted by the earth and was being
used to measure the satelliteâs attitude, the satelliteâs attitude control system maintained the axis of the
science telescope and the satellite roll axis, which coincided with the science telescopeâs optical axis,
within 200 mas (milli-arc-sec) of the apparent position of the guide star. Because of the 20 arc second
annual optical aberration and the relativistic linear drift in the gyroscopesâ spin axis orientation, the angle
between the roll axis and the spin axis slowly varied over the course of the year but was maintained
within 30 arc seconds during normal operations.
Each of the 17 spacecraft operations included (1) measurement of the gyroscope orientation for 24 hours,
(2) maneuver of the satellite roll axis to a nearby star or virtual star, (3) maintenance of the gyroscope
spin axis in this orientation for a period of 12 or 24 hours, (4) another spacecraft maneuver back to the
guide star, IM Pegasi, and, finally, (5) another period of 24 hours to measure the gyroscope spin axis
orientation after the spacecraft maneuvers. The gyroscope drift rate at the increased misalignment could
then be found from the change in each gyroscope spin axis orientation divided by the time at the increased
misalignment. The majority of these spacecraft maneuvers were to stars or virtual stars within one degree
of the IM Pegasi (HR Pegasi is 0.4
deg to the East of IM Pegasi, and HD 216635 is 1 deg to the North of
IM Pegasi), but maneuvers were also made to one star as far as 7 deg from IM Pegasi.
In addition to increasing the misalignment between the satellite roll axis and the gyroscope spin axis,
during some of the these operations the spacecraft was deliberately accelerated as much as 10
-7
m/s
2
using
the helium thrusters to measure the combined effect of increasing the misalignment and the spacecraft
acceleration. Furthermore, at enhanced misalignment angles, tests were done to measure the effect of
increasing the 20 Hz control voltage of the electrostatic suspension system from its nominal 200 mV
level. The operating mode of the electrostatic suspension system was also changed so that the mean
control voltages on the pair of electrodes on one axis were held at a constant d.c. value rather than being
modulated at 20 Hz.
The average gyroscope drift rate for gyroscope 3 and the orientation of the spacecraft relative to IM
Pegasi during the maneuver are shown in Figure 1 for those operations where the satellite roll axis lay
within 1 degree of the guide star, and where the control voltages applied by the electrostatic suspension
system were modulated at 20 Hz. In the polar plot, each vector represents the magnitude and direction of
the average drift rate during the maneuver, and the base of the vector is at the location of the
misalignment between the roll axis and spin axes during the maneuver. In this figure, note that the
gyroscope drift rate lies in a direction perpendicular to the misalignment and is approximately
proportional to the magnitude of the misalignment. The right hand panel of Figure 1 shows the measured
Keiser et al: Misalignment and Resonance Torques and Their Treatment in the GP-B Data Analysis 41
magnitude of the drift rate compared to the magnitude of the misalignment for gyroscope 3. The drift rate
is proportional to the misalignment for misalignments less than 0.4 degrees but begins to show evidence
of nonlinearity at a misalignment 1 degree. At angles larger than 1 degree, the relation between the
misalignment and the drift rate was clearly nonlinear although the direction of the drift rate of the
gyroscope spin axis continues to be in a direction perpendicular to the misalignment.
Figure 1
. Mean Drift Rate vs. Misalignment for Gyroscope 3 During Calibration Maneuvers
Changing the operating mode of the electrostatic suspension system provided clear evidence that the
source of these torques is electrostatic. While the voltages applied to the electrodes were modulated at 20
Hz, the proportionality factor between the misalignment angle and the gyroscope drift rate remained the
same regardless of the magnitude of the modulated 20 Hz control voltage applied to the electrodes and
regardless of the acceleration of the spacecraft. However, when steady voltages were instead applied to
the electrodes, there was a significant change in the magnitude of the drift rate for a given misalignment
angle, although the direction of the gyroscope spin axis drift rate continued to lie in a direction
perpendicular to the misalignment. In addition, the magnitude and direction of the drift rate change
depended on the magnitudes of the steady voltages applied to the six electrodes.
2.2
Explanation and Calculation of Torque
The electrostatic patch effect (Darling (1989)) produces a nonuniform potential on the surface of a metal.
Charges in a metal with a finite conductivity move until the electrostatic potential is uniform. However, if
there is a nonuniform dipole layer on the surface of the metal due to the crystalline structure or impurities
on the surface of the metal, then the electrostatic potential on the surface is no longer uniform. In this
case, the electric field and the force are no longer necessarily perpendicular to the surface of the metal.
For an isolated metal surface, the net force and torque on the body are zero, but with two metallic surfaces
in close proximity, the net force or torque on each surfaces will, in general, be nonzero. For two plane
parallel surfaces, the net force due to patch effect fields has been calculated by Speake (Speake (1996)).
To investigate the effect of a nonuniform potential on the surface of the rotor and the housing, the
electrostatic potential on the surface of the rotor and the housing was expanded in terms of spherical
harmonics:
1000
2000
3000
4000
30
210
60
240
90
270
120
300
150
330
180
0
M
ea
n N
ort
h-S
ou
th
M
isa
lig
nm
en
t
Mean West-East Misalignment
0
500
1000
1500
2000
2500
3000
3500
4000
0
0.5
1
1.5
2
2.5
3
3.5
4
Gyro 3, Drift Rate Magnitude vs. Mean Misalignment
D
rif
t R
at
e M
ag
ni
tu
de
(a
rc
s
ec
/d
ay
)
Mean Misalignment (arc sec)
k = 2.5 arc sec/day/degree
Gyroscope 3, Mean Rate (mas/day) vs. Mean Misalignment (as)
42 Gravity Probe B Science ResultsâNASA Final Report, December 2008
)
,
(
)
,
(
)
,
(
)
,
(
0
0
Ď
θ
Ď
θ
Ď
θ
Ď
θ
lm
l
l
l
m
lm
H
lm
l
l
l
m
lm
R
Y
H
V
Y
R
V
â â
â â
â
=
â
=
â
=
â
=
=
=
(1)
Where
R
lm
and
H
lm
are coefficients in the expansion of the potentials in terms of the orthonormal spherical
harmonics,
Y
lm
(
θ
,
Ď
)
. With these boundary conditions, Laplaceâs equation may be solved for the potential,
ÎŚ
, in the region between the two surfaces. With this solution and the approximation that the gap,
Î
,
between the rotor and the housing is much smaller than the radius,
a
, of the rotor, the gradient in the
potential at the surface of the rotor and the housing is given by
(
)
)
,
(
1
0
Ď
θ
lm
l
l
l
m
lm
lm
b
r
a
r
Y
R
H
r
r
â â
â
=
â
=
=
=
â
Î
â
â
ÎŚ
â
=
â
ÎŚ
â
(2)
The energy stored in the electric field may be calculated from the potential and the gradient in the
potential at the two surfaces:
( )
(
)
(
)(
)
.
2
2
2
2
*
0
2
0
2
0
0
2
0
lm
lm
l
lm
lm
a
r
a
b
b
r
a
r
gap
R
H
R
H
a
r
d
a
dA
r
dA
r
dV
W
â
â
Î
=
â
ÎŚ
â
ÎŚ
â
ÎŚ
Ί
â
ââ
â
â
ââ
â
â
â
ÎŚ
â
ÎŚ
â
â
ÎŚ
â
ÎŚ
=
ÎŚ
â
=
â
âŤ
âŤ
âŤ
âŤ
â
=
=
=
=
Îľ
Îľ
Îľ
Îľ
(3)
In this expression,
Îľ
0
is the permittivity of free space. Here and in the equations which follow, the asterisk
denotes the complex conjugate of the coefficient. The orthonormal properties of the spherical harmonics
have been used here: it is important to note that the result (3) is valid only if the spherical harmonic
expansions for both surfaces are written in the same reference frame.
The expansion of both surface potentials are originally done in the frame fixed to the corresponding
surface. Their coefficients do not change as long as the patch potential remains constant. To calculate the
energy stored in the electrostatic field and the gyroscope torques, rotation matrices may be used to find
the expansion of the surface potential in a common reference frame. In the principal axis reference frame,
with the z
â˛
-axis chosen as the principal axis corresponding to the maximum moment of inertia, the rotor
potential may be expanded in terms of spherical harmonics in that reference frame (r,
θâ˛
,
Ďâ˛
):
.)
'
,'
(
'
)'
,'
(
0
â â
â
=
â
=
=
l
l
l
m
lm
lm
R
Y
R
V
Ď
θ
Ď
θ
(4)
Since the rotor potential is real,
R
lm
â˛
= (-1)
m
R
l,-m
â˛
. We chose a common reference frame to be the guide
star reference frame, where the
z-
axis lies in the apparent direction of the guide star and the
x-
axis lies in
the plane of the
z-
axis and the Earthâs rotation axis. The spherical harmonics in the principal axis
reference frame may be rotated to the guide star reference frame through two successive sets of Euler
angle rotations. The first is a transformation from the principal axis reference frame to a reference frame
with the
z-
axis aligned with the instantaneous spin axis. This 3-2-3 transformation uses the three Euler
angles
Ď
P
,
Îł
,
and
Ď
S
which correspond to a rotation about the principal axis with the largest moment of
inertia, a rotation about the new
y-
axis by an angle
Îł
which is the angle between the principal axis and the
spin axis, and finally a rotation about the spin axis by the spin phase,
Ď
S.
The second set of Euler angle
transformations is a 2-3 transformation: the first angle rotates about the
y-
axis of the coordinate system by
an angle
-
β
so the
z-
axis is aligned with the apparent direction to the guide star, while the second angle
Keiser et al: Misalignment and Resonance Torques and Their Treatment in the GP-B Data Analysis 43
rotates about the direction to the guide star by and angle -
Îą
. These two angles,
Îą
and
β
, determine the
orientation of the gyroscope spin axis in the guide star reference frame. With these two sets of rotations,
the rotor potential in the guide star reference frame becomes
â â
â
â
â â
â
â
â
=
â
=
â
=
â
=
â
â
=
â
=
â
=
â
=
â
=
â
â
â
=
0
0
)
,
(
)
(
)
(
'
)
,
(
)
,
,
(
)
0
,
,
(
'
)
,
(
l
l
l
m
l
l
q
l
l
p
lq
im
l
pm
ip
l
qp
iq
lm
l
l
l
m
l
l
q
l
l
p
lq
P
S
l
pm
l
qp
lm
R
Y
e
d
e
d
e
R
Y
D
D
R
V
P
S
Ď
θ
Îł
β
Ď
θ
Ď
Îł
Ď
β
Îą
Ď
θ
Ď
Ď
Îą
(5)
The elements of the rotation matrices,
D
and
d
, are defined in Rose (1995). If the housing potential is
initially expanded in a reference frame where the
z-
axis coincides with both the satellite roll axis and the
apparent direction to the guide star, then transforming to the inertially fixed reference frame involves only
a rotation about the
z-
axis by the satellite roll phase,
Ď
R
. In the guide star reference frame, the housing
potential may be written as
)
,
(
)
,
(
0
'
Ď
θ
Ď
θ
Ď
kj
ij
k
k
k
j
H
Y
e
H
V
R
j
k
â â
â
=
â
=
=
(6)
With these two expressions for the rotor and housing potentials in a common reference frame, the variable
part of the energy in the electrostatic field averaged over the spin of the gyroscope is
â â â
â
=
â
=
â
=
+
â
â
Î
â
=
0
)
(
0
0
*
'
'
2
0
)
(
)
(
l
l
l
m
l
l
q
iq
im
l
m
l
q
lq
R
P
lm
e
e
d
d
H
R
a
W
Ď
Îą
Ď
Îł
β
Îľ
(7)
In addition to averaging over the spin of the gyroscope, the energy in the electrostatic field may be
averaged over the roll of the housing and the polhode period. In this case, the averaged energy is
)
(
)
(
)
(
)
(
'
0
'
2
0
00
00
'
0
'
2
0
0
0
0
0
Îł
β
Îľ
Îł
β
Îľ
â
Î
â
=
â
Î
â
=
â
â
â
=
â
=
l
l
l
l
l
l
P
P
H
R
a
d
d
H
R
a
W
l
l
l
l
(8)
Here
P
l
is the Legendre polynomial of order
l
.
The torque acting on the gyroscope may be found from the derivative of the average energy with respect
to the angles,
Îą
and
β
, which define the orientation of the gyroscope spin axis in the guide star reference
frame. Since the average energy given above only depends on the angle
β
, the torque is given by
β
β
β
β
Îł
Îľ
β
Ď
e
P
P
H
R
a
e
W
l
l
l
l
l
Ë
)
(
)
(
Ë
1
'
0
'
0
2
0
â
â
â
Î
=
â
â
â
=
â
â
=
r
(9)
The torque is the negative derivative of the energy with respect to the angle. The unit vector,
e
β
, lies in a
direction perpendicular to the plane of the misalignment angle,
β
, so this torque will always be
perpendicular to the misalignment. From this result, the torque is clearly a nonlinear function of the
misalignment angle, but for small angles this expression becomes
1
),
(
)
1
(
Ë
1
'
0
'
0
2
0
<<
â
+
Î
=
â
â
=
β
Îł
β
Îľ
Ď
β
l
l
l
l
P
H
R
l
l
e
a
r
(10)
44 Gravity Probe B Science ResultsâNASA Final Report, December 2008
This result agrees with the observations. The torque always lies in a direction perpendicular to the
misalignment, and, for small angles, the torque is proportional to the misalignment. In addition, averaged
over the polhode period, the magnitude of the torque for a given misalignment will slowly change if the
angle,
Îł
, between the rotorâs principal axis and the spin axis slowly changes even though the patch effect
potentials on the surface of the rotor and the housing remain fixed. Since the absolute value of the
Legendre polynomial is always less than or equal to unity, for a given value of
l,
the polhode motion
reduces the average torque compared to the torque in the absence of polhode motion (
Îł
= 0
). The
magnitude of the torque depends on the magnitude and distribution of the patch effect potentials on the
surface of the rotor and the housing, but it is important to note that this torque is due to the interaction of
the patch effect on the rotor with the patch effect on the housing. Neither surface alone will produce a
torque.
In addition, this result explains the change in the magnitude of the torque when the operating mode of the
electrostatic suspension system was changed. Since the torque is linearly proportional to the potential on
the surface of the housing, a modulated voltage applied to the electrodes will, on the average, produce no
net torque regardless of its magnitude. However, when a steady voltage is applied to the electrodes, the
electrode potentials are the sum of the electrode patch effect potential and the voltage applied. This
change in the electrode potential produces a corresponding change in the magnitude of the misalignment
torque. The magnitude of the torque depends on the magnitude of the applied steady voltage as well as the
orientation of each electrode with respect to the gyroscope spin axis.
2.3
Data Analysis in the Presence of Misalignment Torques
In an inertial reference frame (Silbergleit (2001)), the equations for the components of the gyroscope drift
rate in the North-South and West-East directions in the presence of a misalignment torque may be written
as
NS
WE
WE
WE
NS
NS
k
r
dt
ds
k
r
dt
ds
Îź
Îź
â
=
+
=
(11)
Where
s
NS
and
s
WE
are the components of the gyroscope spin axis in the North-South and West-East
directions,
r
NS
and
r
WE
are the uniform components of the gyroscope drift rate,
Îź
NS
and
Îź
WE
are the
components of the misalignment,
β
, between the spin and roll axes in the North-South and West-East
directions. From equation (10), this coefficient,
k
, is given by
).
(
)
1
(
1
1
'
0
'
0
2
0
Îł
Îľ
Ď
â
+
Î
=
â
â
=
l
l
l
l
s
P
H
R
l
l
a
I
k
(12)
This coefficient may be time dependent because the polhode angle,
Îł
, is changing due to polhode
damping even though the patch effect potentials are fixed on the surface of the rotor and the housing.
There are several methods of analyzing the data to determine simultaneously the uniform components of
the drift rate and the misalignment torque coefficient based on the time history of the gyroscope spin axis
orientation. One method is to estimate the uniform components of the drift rate as well as the
misalignment torque coefficient, using an appropriate model of the time variation of the latter. This
method along with preliminary results is described in the accompanying paper by Heifetz et al. (2009). A
second, complementary method takes advantage of the linear combination of the two drift rate equations
(11) which eliminates the torque coefficient. Defining the misalignment phase as
Keiser et al: Misalignment and Resonance Torques and Their Treatment in the GP-B Data Analysis 45
NS
WE
Îź
Îź
θ
1
tan
â
=
.
(13)
the linear combination of the drift rate equations
θ
θ
θ
θ
sin
cos
sin
cos
WE
NS
WE
NS
R
r
r
dt
ds
dt
ds
dt
ds
+
=
+
=
(14)
is independent of the misalignment torque coefficient. For any given misalignment phase, it is possible to
determine this radial component of the drift rate without any knowledge of the torque coefficient.
Furthermore, with sufficient variation in the misalignment phase, it is possible to independently determine
the two drift rate components,
r
NS
and
r
WE
.
Equation (14) describes the rate of change of the gyroscope spin axis orientation in the presence of a
continuously changing, but known, misalignment. In the analysis of the data, the rate of change of the
spin axis orientation must be determined from the measured orientation. This rate may be estimated by
dividing the data into batches of 4 or 5 days and finding the rate of change of the orientation over this
interval. Then, the radial rate is a sinusoidal function of the average misalignment phase. The amplitude
and phase of this sinusoidal curve determine the uniform components of the gyroscope drift rate, r
NS
and
r
WE
.
Instead of dividing the data into 4 or 5 day batches, a better method is to find the gyroscope spin
axis drift rate on an orbit-by-orbit basis. In this case, the gyroscope drift rates in the North-South and
West-East directions may be determined from the first differences in the spin axis orientation,
.
1
1
i
i
i
i
i
i
t
t
s
s
dt
ds
r
â
â
=
=
+
+
(15)
The noise in these gyroscope drift rates is sequentially correlated, and it is essential to explicitly include
this sequential correlation in the measurement noise matrix. [See, for example, Reasenberg (1972)]. The
first differences in the orientation estimates may then be used to construct a series of derived
measurements of the gyroscope drift rate. These drift rate measurements are a function only of the
uniform components and the misalignment phase. If the misalignment phase is known, then a least
squares fit may be used to determine the uniform components of the drift rate as long as the sequential
correlation and the dependence on the misalignment phase are explicitly included.
3 Resonance Torques
3.1
Observations
In the course of the data analysis, we found that significant changes in the orientation of the gyroscope
spin axis occurred during those intervals when a harmonic of the gyroscope polhode frequency was nearly
equal to the satellite roll frequency. Although measurable changes did not always occur at those times
when this resonance condition held, each measurable change occurred at one of these resonance
conditions.
Figure 2 shows the magnitude of the components of the gyroscope 2 readout signal at the cosine and sine
of the roll frequency for each orbit during one of these resonances. The cosine component is proportional
to the angle between the gyroscope spin axis and the apparent direction to the guide star in the plane of
the orbit, while the sine component is proportional to the same angle in the perpendicular direction.
Figure 3 is a plot of these components versus one another.
46 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Figure 2.
Amplitude of Components of Gyroscope Readout Signal at Cosine and Sine of Satellite Roll Frequency
Figure 3
. Components of Gyroscope Readout Signal at Satellite Roll Frequency. Each point represents the
estimated average orientation for one orbit of the satellite.
5660
5680
5700
5720
5740
5760
-2
-1
0
1
2
Cosine of Roll After Removing DC and Linear Drift, Gyro 2, Res 277
m
V
5660
5680
5700
5720
5740
5760
-2
-1
0
1
2
Sine of Roll
m
V
Orbit Number
Gyroscope 2
From: May 10, 2005, 2am
To: May 15, 2005, 5 pm
1 day
50 mas
-1.5
-1
-0.5
0
0.5
1
1.5
-1.5
-1
-0.5
0
0.5
1
1.5
Sine and Cosine at Roll Frequency, Gyro 2, Res 277
C
os
in
e o
f R
ol
l F
re
qu
en
cy
(m
V)
Sine of Roll Frequency (mV)
North
West
Approximate Scale:
50 mas
Keiser et al: Misalignment and Resonance Torques and Their Treatment in the GP-B Data Analysis 47
3.2
Explanation and Calculation of Torque
From equation (7), the spin averaged torque may be calculated during those intervals where a harmonic of
the gyroscope polhode frequency is equal to the satellite roll frequency. In this case, there are two
components of the torque because the spin averaged energy in the electrostatic field depends on both the
angles,
Îą
and
β
, which define the orientation of the gyroscope spin axis with respect to the guide star
reference frame:
(
)
â â â
â â â
â
=
â
=
â
=
â
â
â
=
â
=
â
=
â
â
â
â
â
â
â
â
â
â
â
â
â
Î
=
â
â
â
=
â
Î
â
=
â
â
â
=
0
0
0
*
'
'
2
0
0
0
0
*
'
'
2
0
)
(
)
(
Re
)
(
)
(
Re
l
l
l
m
l
l
q
iq
im
l
m
l
q
iq
lq
l
l
l
m
l
l
q
iq
im
l
m
l
q
iq
lq
R
P
lm
R
P
lm
e
e
d
d
e
H
R
a
W
e
e
d
d
iqe
H
R
a
W
Ď
Ď
Îą
β
Ď
Ď
Îą
Îą
Îł
β
β
Îľ
β
Ď
Îł
β
Îľ
Îą
Ď
(16)
A third component of the spin-averaged torque, which lies along the instantaneous spin axis, is zero.
These components of the torque are not orthogonal. The torque in the guide reference frame in the two
perpendicular directions which are also perpendicular to the direction to the guide star may be found from
the relations
Îą
β
Îą
β
Îą
Ď
Ď
Îą
Ď
β
Îą
Ď
Ď
Îą
Ď
β
Îą
Ď
Ď
=
+
â
=
â
â
=
z
y
x
cos
cot
sin
sin
cot
cos
(17)
Combining these last two equations and transforming to the inertial reference frame (Silbergleit, 2001),
the components of the spin axis drift rate in the North-South and West-East directions are
Ď
Ď
Ď
Ď
Î
+
Î
=
Î
â
Î
=
cos
sin
sin
cos
m
m
WE
m
m
NS
b
a
dt
ds
b
a
dt
ds
(18)
where
ÎĎ
is the phase difference between the m
th
harmonic of the polhode phase and the satellite roll
phase. For small misalignment angles,
β
<< 1
, the coefficients a
m
and b
m
are given by
{
}
{
}
1
,
0
0
2
0
1
,
0
0
2
0
Im
)
(
)
1
(
Re
)
(
)
1
(
â
â
=
â
â
=
â
â
+
Î
=
+
Î
=
l
lm
l
l
m
m
l
lm
l
l
m
m
H
R
d
l
l
a
b
H
R
d
l
l
a
a
Îł
Îľ
Îł
Îľ
(19)
48 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Figure 4
. Expected Change in Orientation due to Resonance Torques. The units in this figure are arbitrary but are
equal in the North-South and West-East directions.
Assuming the polhode frequency is changing approximately linearly with time at a rate r close to the
resonance, the phase difference changes quadratically with time, the time history of the orientation of the
gyroscope spin axis may be found by integrating the drift rate equations
2
cos
2
sin
)
(
2
sin
2
cos
)
(
2
2
0
2
2
0
rt
b
rt
a
s
t
s
rt
b
rt
a
s
t
s
t
m
t
m
WE
WE
t
m
t
m
NS
NS
âŤ
âŤ
âŤ
âŤ
â
â
â
â
â
â
â
â
+
â
=
+
+
=
(20)
The integrals in these equation are Fresnelâs integrals, and the path traced out by the components of the
gyroscope spin axis is the Cornu spiral as shown in Figure 4. This curve agrees well with the
measurements shown in Figure 3.
3.3
Data Analysis in the Presence of Resonance Torques
As long as the polhode phase is known well enough, the times at which the resonances occur and the
interval over which the phase difference is less than 2
Ď
radians may be accurately determined. Outside of
this interval, although these resonance torques have a constant amplitude but time dependent frequency,
the deviation of the gyroscopeâs spin axis orientation is less a small fraction of the magnitude of the total
offset produced by the resonance torque. The most straightforward approach to the analysis of the data in
the presence of these resonance torques is simply to exclude the data during these intervals and to assume
that the torques produce an unknown change in the orientation during this time. With such a data analysis
approach, it is, of course, necessary to determine the sensitivity of the overall result to the width of the
interval over which the data is excluded.
-1
0
1
2
3
4
5
6
7
-3
-2
-1
0
1
2
3
Orientation of Gyroscope Spin Axis
N
or
th
-S
ou
th
O
rie
nt
at
io
n
West-East Orientation
Keiser et al: Misalignment and Resonance Torques and Their Treatment in the GP-B Data Analysis 49
An alternative approach would be to use the data during each of the resonances to determine the two
coefficients,
a
m
and
b
m
, for each resonance. Although this approach increases the number of coefficients
which need to be determined from the data analysis, the advantage is that additional, unsegmented data
may be used. Covariance analyses have shown that there is a net increase in the precision of the results
when data during the resonance is included and the two coefficients are determined along with the other
parameters. The caveat here is that it is necessary to clearly demonstrate that the model used for the
orientation or rate during each of these resonances is accurate.
4 Summary and Conclusion
The misalignment and resonance torques produced significant disturbances to the orientation of the
gyroscope spin axes. These torques may be explained by the same physical effect â an interaction of patch
effect fields on the surface of the rotor with patch effect fields on the surface of the housing. The
misalignment torque is proportional to the misalignment between the gyroscope spin axis and the satellite
roll axis. All four gyroscopes exhibited a drift rate due to the misalignment torque, but the magnitude of
the torque varied considerably from one gyroscope to another and over the duration of the mission. The
resonance torques, on the other hand, are independent of the misalignment to lowest order, but their
frequency depends on the difference between a harmonic of the polhode frequency and the satellite roll
frequency. When this frequency difference is small, they can produce a permanent offset of the gyroscope
spin axis. Because of the significantly different polhode frequencies and rates of change of the polhode
frequencies, the frequency of the offsets due to the resonance torques was very different for the different
gyroscopes. As the polhode motion was damped out, the frequency of these resonance conditions
decreased.
In the analysis of the telemetry data from the satellite, the effects of these torques may be clearly
separated from the uniform drift rate. In the case of the misalignment torques, the drift rate due to these
torques may be separated from the uniform drift rate because of the known direction of the misalignment
drift rate which is perpendicular to the misalignment. Two alternative data analysis methods can be used:
one only uses the component of the drift rate in a direction parallel to the misalignment, the other uses the
drift rate information in both directions but includes the misalignment torque coefficient as one of the
parameters in the data analysis. The effects of the resonance torques may also be separated from the
uniform drift rate either by not using the data during those times when the effects of the resonance torques
are most pronounced or by explicitly included two additional parameters for each resonance in the data
analysis.
Acknowledgements
The authors would like to acknowledge the remarkable contributions of numerous people at Stanford
University, Lockheed-Martin, and NASA to the construction and operation of the GP-B Satellite. We
would also like to acknowledge S. Buchman, B. Clarke, J. Conklin, F. Everitt, M. Heifetz, D. Hipkins, J.
Li, B. Muhlfelder, Y. Ohshima, J. Turnearue, and P. Worden for their insights and contributions to this
paper.
50 Gravity Probe B Science ResultsâNASA Final Report, December 2008
References
Adler, R. J. and A. S. Silbergleit (2000). "General Treatment of Orbiting Gyroscope Precession."
International Journal of Theoretical Physics
39
(5): 1291-1316.
S. Buchman, J. P. Turneaure, et al. (2009). to be published.
Darling, T. W. (1989). Electric Fields on Metal Surfaces at Low Temperatures. School of Physics.
Parkville, Victoria 3051 Australia, University of Melbourne
:
88.
de Sitter, W. (1916). "On Einstein's theory of Gravitation and its Astronomical Consequences." Monthly
Notices Royal Astronomical Society
77
: 155-184.
Everitt, C. W. F. et al. (2009). Space Science Reviews (Chapter 1 of this report)
Heifetz, M., W. Bencze, et al. (2009). "The Gravity Probe B Data Analysis Filtering Approach", Space
Science Reviews (Chapter 4 of this report).
Lense, J. and H. Thirring (1918). "On the Influence of the Proper Rotations of Central Bodies on the
Motions of Planets and Moons According to Einstein's Theory of Gravitation." Zeitschrift fur
Physik
19
: 156. English Translation: Mashhoon, B., F. W. Hehl, et al. (1984). "On the
Gravitational Effects of Rotating Masses: The Lense-Thirring Papers." General Relativity and
Gravitation
16
: 711-750.
Muhlfelder, B. (2009), Space Science Reviews (Chapter 5 of this report).
Reasenberg, R. (1972). "Filtering with Perfectly Correlated Measurement Noise." AIAA Journal
10
(7):
942-3.
Rose, M. E. (1995). Elementary Theory of Angular Momentum. New York, Dover.
Silbergleit, A. S., M. I. Heifetz, A. S. Krechetov (2001). "Model of Starlight Deflection and Parallax for
the GP-B Data Reduction", Gravity Probe B Document S0393, Stanford University, Stanford,
CA.
Silbergleit, A. S., J. Conklin et al. (2009). "Polhode Motion, Trapped Flux, and the GP-B Science Data
Analysis", Space Science Reviews (Chapter 2 of this report).
Speake, C. C. (1996). "Forces and force gradients due to patch fields and contact-potential differences."
Classical and Quantum Gravity
13
: A291-A297.
51
4
t
he
G
ravity
P
robe
b D
ata
a
nalysis
f
ilterinG
a
PProach
Michael Heifetz et al.
52 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Heifetz et al.: The Gravity Probe B Data Analysis Filtering Approach 53
The Gravity Probe B Data Analysis Filtering Approach
M.Heifetz, W.Bencze, T.Holmes, A.Silbergleit, V.Solomonik
Abstract
A simple strategy of the Gravity Probe B (GP-B) data analysis has evolved in the elaborate multi-level
structure after the discovery of the complex polhode motion, and of the patch effect torques. We describe
a cascade of four estimators (filters) that reduce the science data (SQUID and telescope signals) to the
estimates of the relativistic drift rates. Those estimators, structured in two âfloorsâ, are based on the
polhode-related models for the readout scale factor and patch effect torque. Results of the 1
st
Floor
processing - gyro orientation profiles - manifest clearly the strong geodedic effect but also the presence of
classical torque. Modeling of the patch effect torque at the 2
nd
Floor provides a successful compensation
of the torque contributions, and leads to consistent estimates of the relativistic drift rates.
1 Introduction
Data analysis in the GP-B experiment - estimation of the relativistic drift rate for each of the four unique
GP-B gyroscopes - was supposed to be rather straightforward [Haupt (1996), Heifetz et al. (2003)]. One
should take the low frequency (LF) SQUID readout signal for the duration of the mission, calibrate the
scale factor based on orbital and annual aberration known from the GPS orbit data and earth ephemeris;
obtain the time-history of the gyroscope inertial orientation in the projections on the directions in the
orbital plane and perpendicular to it, and then estimate the slopes of those projections, which were
supposed to be almost straight lines. The slope of the line in the orbital plane gives the measurement of
the geodetic effect. The slope of the projection perpendicular to the orbital plane gives the measurement
of the frame-dragging effect [this is exactly true for the ideal polar orbit with the Guide Star (GS) in the
orbit plane].
In reality, an unexpected damped polhode motion of all the GP-B rotors, and some larger than expected
classical torques on them, were discovered on orbit. In addition, about a year of science data was cut into
10 segments by various spacecraft (S/C) anomalies. This has turned a âsimpleâ data analysis strategy into
a challenging adaptive estimation process involving a multi-level filtering machinery.
There are three cornerstones of the GP-B filtering (estimation) method. 1) SQUID readout signal
structure:
measurement models
â these models are based on the underlying physics and engineering of the
GP-B instrument. 2) Gyroscope motion:
torque models
â accumulated understanding of the underlying
physics. 3) Filter implementation:
numerical techniques
â algorithms from the estimation theory [Kailath
et al. (2000), Bierman (2006)] specifically adjusted to GP-B data.
To successfully estimate the relativity parameters, the filtering approach must meet several major
challenges: a) the complicated variation of the readout scale factor (see details in Sect. 3); b) a
continuously acting torque due to electrostatic patch effect (see Sect. 4); c) simultaneous estimation of the
relativity parameters, multiple torque and scale factor coefficients; d) noisy measurements that depend on
unknown parameters in a nonlinear fashion; e) segmented science data (see Sect. 3); f) the large volume
of data (~one terabyte): four gyroscopes x 4605 orbits of science data.
In Sect. 2 the main original idea of a âsimpleâ data analysis is explained and the science signal
measurement model is introduced. Sect. 3 presents the current two-floor structure of the filtering
approach, gives the description of the model of scale factor variations, and shows the results of the
54 Gravity Probe B Science ResultsâNASA Final Report, December 2008
science signal analysis. Those are the gyroscope orientation profiles revealing the presence of Newtonian
torques. A model of the patch effect torque is introduced in Sect. 4, and a concept of the âreconstructionâ
of the gyroscopeâs relativistic motion by means of modeling, estimation, and elimination of the torque
contributions from the orientation profiles is described. A further expansion of the data analysis, the âtwo-
second filterâ, is described in Sect. 5. Section 6 contains some conclusions.
2 âSimpleâ GP-B Data Analysis: Pre-Launch Concept
2.1
Inertial frame and the relevant vectors
The drift experienced by the GP-B gyroscopes is measured in the
inertial Cartesian frame {
NS
EW
GS
e
e
e
Ë
,
Ë
,
Ë
} introduced by Silbergleit et al.
(2001) (see Fig. 1). Its first unit vector
GS
e
Ë
points from the Earth
center to the true position of the GS (IM Pegassi). The second unit
vector is
z
e
z
e
e
GS
GS
EW
Ë
Ë
/
)
Ë
Ë
(
Ë
Ă
Ă
=
, where
z
Ë
is the unit vector of the
inertial frame JE2000 [see Seidelmann (1992), McCarthy (1996)]. The
last axis
NS
e
Ë
is naturally defined as
GS
EW
NS
e
e
e
Ë
Ë
Ë
Ă
=
. The
z
Ë
axis
coincides with the Earth rotation axis at noon, January 1, 2000, and
stays very close to it later.
Thus the geodetic drift is in the direction of
NS
e
Ë
, and the frame-
dragging drift is along
EW
e
Ë
, if the orbit is polar and contains
GS
e
Ë
in its
plane. The real GP-B orbit is close enough to this ideal one, so the
frame is a natural choice for GP-B data analysis.
The unit vector
s
Ë
represents the gyroscope spin axis, while
Ď
Ë
points along the S/C roll axis. The
misalignment angle
,
Ď
, between these vectors is only ~10
-4
or smaller, and they are both at about the same
small angle to
GS
e
Ë
. For this reason, the
misalignment vector,
Îź
Ë
, defined as
Ď
Îź
Ë
Ë
Ë
â
=
s
, lies, to lowest
order, in the NS â EW plane, and
Îź
Ď
â
~ 10
-4
.
When the GS is visible by the telescope (GS valid, GSV), the S/C attitude and translation control system
points the S/C roll axis in the apparent direction to the GS. The latter differs from true direction to the GS
by the sum of the orbital and annual aberrations (up to 25 arcs), the pointing error,
θ
,
(within ~50 marcs
on average), and the contribution of the stellar parallax and starlight bending by the Sun (up to 25 mas).
In the pre-launch data reduction concept, only GSV periods were planned to be analyzed.
2.2
Measurement model
The GP-B readout system provides measurements of the time-varying angle between the direction of the
gyroscope spin axis aligned with the London moment [Turneaure et al. (2003), Everitt et al. (2008)], and
the SQUID pick-up loop that rolls with the S/C. The corresponding measurement model is shown, to
lowest order, to be [Silbergleit (2006)]:
[
]
[
]
.
)
)
(
sin(
)
(
)
(
)
)
(
cos(
)
(
)
(
)
(
)
(
noise
bias
t
t
s
t
t
t
s
t
t
C
t
Z
r
EW
EW
r
NS
NS
g
+
+
â
âŹ
âŤ
âŠ
â¨
â§
+
ÎŚ
â
+
+
+
ÎŚ
â
=
δĎ
Ď
δĎ
Ď
Figure 1
. Inertial Reference Frame.
(1)
Heifetz et al.: The Gravity Probe B Data Analysis Filtering Approach 55
There are two groups of variables in this equation. The first group consists of the available data:
)
(
t
Z
is
the SQUID signal,
)
(
t
r
ÎŚ
is the roll phase signal, and
)
(
),
(
t
t
EW
NS
Ď
Ď
are known (see the previous section;
pointing error is determined by the telescope signal). The second group of the unknown variables that
need to be estimated includes the readout scale factor
)
(
t
C
g
, the gyroscope inertial
orientation
)
(
),
(
t
s
t
s
EW
NS
, the roll phase offset
δĎ
, and the bias components, such as the SQUID
calibration signal.
In the original pre-launch strategy of the data analysis, the idea of a âsimpleâ batch-based estimator was
implemented. The science signals were processed in a sequence of short batches; within each batch, the
scale factor and the roll phase offset can be calibrated using the known total aberration, and subsequently,
the estimates
)
(
~
t
s
NS
and
)
(
~
t
s
EW
of the gyro orientation time history could be obtained. The estimated
curve of each of the orientations would then be fitted to a line, and the estimated slopes would give the
measurement of the relativity parameters.
Two surprising on-orbit findings led to the development of a much more elaborate two-floor structure of
the data analysis described below.
3 Two-Floor Structure of Data Analysis
3.1
Modification of the pre-launch scheme
Pre-launch data analysis strategy was built on the assumption that no modeling of classical torques is
needed. In this case processing of the SQUID signal could be done in one batch, or, to check the
instrument stability, in a sequence of batches. During the mission, after the discovery of the polhode-
related behavior of the scale factor
)
(
t
C
g
[see Silbergleit et al. (2009)] and of the so called
misalignment
torque [see Keiser et al. (2009)], it became evident that more sophisticated modeling of both the scale
factor and torque is necessary. The magnitude of the misalignment torque is proportional to the (small)
misalignment,
Îź
, and its direction is perpendicular to the misalignment vector (Sect. 2.1). This torque
acted continuously, so it requires continuous modeling for the entire orbit, both its GSV and GSI parts
(GSI=GS Invalid, the period when the GS is occulted by the Earth).
To simplify the development of a cumbersome analysis structure, save the computational resources and
stabilize filtering algorithms, we decided to carry out data analysis in two steps called
floors
. The 1
st
Floor
is a modification of the pre-launch scheme using batch processing with no torque modeling within each
batch. At the 2
nd
Floor torque and more accurate scale factor modeling occurs, as a way of integrating
information from all the batches. For various reasons, we have chosen the batch length equal to one orbit
(~97 min).
The implementation of the two-floor analysis structure faced an additional difficulty in the form of
several different characteristic time scales present in the signals. The longest of them is 1 year, the period
of the annual aberration. Unfortunately, the year-long duration of the experiment was interrupted by
several spacecraft anomalies (e.g. the solar flare on January 20, 1995), so science data came in 10
segments separated by day long gaps (the list of segments is given in Table 1). Naturally, one of the
questions to be resolved was how to connect those segments in the data analysis.
56 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Segment
Science Data
Duration (days)
1
September 13, 2004 â September 23
11
2
September 25 â November 10
47
3
November 12 â December 04
23
4
December 05 â December 09
5
5
December 10 â January 20, 2005
42
6
January 21 â March 04
43
7
March 07 â March 15
9
8
March 16 â March 18
3
9
March 19 â May 27
70
10
May 31 â July 23, 2005
54
Table 1
. Segmented Science Data
The next important time scale is the orbital period, the period of orbital aberration calibrating the scale
factor within each batch. A short interval (~ 2 min) of the GS reacquisition by the telescope is also a part
of the data for each orbit.
The third meaningful time scale is the period of the GP-B satellite roll (77.5 sec): the science signal is at
the roll frequency, as seen from Eq. (1). Other time scales are related to the gyro polhode motion: it is the
polhode period (~ few hours), which itself is changing slowly with the characteristic time of kinetic
energy dissipation (~1 or 2 months, depending on the gyro) [Silbergleit (2009)].
All these time scales were taken into account in the design of the data analysis software architecture. The
orbital period was the main batch length for the 1
st
Floor of the data reduction (see Sect. 3.2 below). The
changing polhode period and the energy dissipation times were the key parameters in the development of
the trapped flux mapping (TFM) [Conklin (2009), Silbergleit et al. (2009)], and the outputs of the TFM,
the polhode phase and polhode angle, were key ingredients of the scale factor and torque modeling (see
Sect. 3.3 and Sect. 4 below). The annual aberration, which is defined by the annual period, is the main
long-term calibrator of the scale factor and also, the main contributor to the misalignment time signature.
The latter is essential for the observability and identification of the misalignment torque coefficients (see
Sect. 4). Lastly, the roll frequency is, naturally, the fundamental frequency of the SQUID science signal.
But in addition, it is also one of the key elements of the roll-resonance torque and the algorithms of its
compensation (see Sect. 4 and Sect. 5).
3.2
Floor 1: one-orbit batch data reduction
The idea of the 1
st
Floor is the compression of science data within each orbit. The inputs are the LF
SQUID and telescope signals, roll phase signal, and pre-computed aberration data. The filtering
(estimation) is carried out independently for each orbit based on the measurement model (1). The state
vector includes the two components of the gyro inertial orientation, the coefficients of the scale factor
polhode variations, the roll phase offset, the two telescope scale factors (of two directional channels), and
the coefficients in the bias variation model. The 1/f nature of the SQUID noise is addressed by applying a
band-pass digital filter to the SQUID and telescope signals. The resulting signal spectrum is within the
roll Âą orbit frequency band. No torque modeling is performed at the 1
st
Floor. The output of the 1
st
Floor
data reduction consists of the estimates of the state vector (1 point per orbit), full information (inverse
Heifetz et al.: The Gravity Probe B Data Analysis Filtering Approach 57
covariance) matrix, and the post-fit residuals. An additional output, the misalignment and pointing
signals, is used for the misalignment torque modeling at the 2
nd
Floor (see Sect. 4).
Figure 2.
First Floor structure block-diagram.
Among various modules in the 1st Floor structure (Fig.2) several should be highlighted:
1)
Scale factor model
2)
Pointing error compensation - gyroscope/telescope scale factor matching
3)
Data grading (using only high quality data)
4)
Bias model (including polhode variations, bias jumps, calibration signal components, etc.)
5)
Nonlinear least-squares estimator. The latter is implemented in the form of the square-root
information filter [Bierman (2006)] that provides numerical stability and computational fidelity.
3.3
Floors 1 and 2: readout scale factor modeling
A mathematical model of the scale factor variations due to trapped magnetic flux and gyro polhoding has
been shown to be (
Îł
Cg -model):
58 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Here
0
ÎŚ
and
0
Îł
are the polhode phase and polhode angle calculated for a symmetric gyroscope based
on the information from the TFM [see Silbergleit et al. (2009) for all the details]. The numbers
M, N
0
, N
n
are chosen so to make these approximations good enough for our analysis accuracy without including
unnecessary terms (in theory, all the expansions are, of course, infinite). The number of terms required
decreases towards the end of the mission, with the damping of polhode variations.
In this model the scale factor is represented as an expansion in harmonics of the polhode phase, but the
coefficients of this expansion turn out to be functions of the polhode angle. The polhode frequency and
angle are slowly changing due to energy dissipation, with the characteristic time of this process (1-2
months). At Floor 1 only
Eq. (2a) is used, and the coefficients
m
a
and
m
b
are estimated once per orbit. At the 2
nd
Floor these
estimates are fit to the model (2b), with the unknown parameters
a
mn
and
b
mn
to be estimated. If the fit is
successful, then the complete physical model (2a), (2b) is justified. The fit results for real data are plotted
in Fig.3; clearly, the model explains the complicated profiles of the coefficients
m
a
and
m
b
very well. An
identical model is used at the 2
nd
Floor for the torque coefficients (Sect. 4). The scale factor model (2a),
(2b) will also be used in a new two-second filter described in Sect. 5.
The successful use of the
Îł
C
g
âmodel to a significant extent is due to TFM that computes a highly
accurate polhode phase and angle for each gyroscope for the whole duration of the science mission. It also
provides an alternative to the estimation of scale factor variations by the model (2a), (2b). Namely, scale
factor variations were computed from the TFM results (again, for each gyro and the whole mission), and
could be used as a known function, simplifying significantly the 1
st
Floor processing.
12/19
03/29
07/07
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
date
d-
le
ss
Gyro 1. Fit Cga 1
12/19
03/29
07/07
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
date
d-
le
ss
Gyro 1. Fit Cgb 1
12/19
03/29
07/07
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
date
d-
le
ss
Gyro 1. Fit Cga 1
12/19
03/29
07/07
-0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
date
d-
le
ss
Gyro 1. Fit Cgb 1
10/30 12/19 02/07 03/29
-2
-1
0
1
2
3
x 10
-3
date
d-
le
ss
Gyro 1. Fit Cga 4
10/30 12/19 02/07 03/29
-2
-1
0
1
2
3
4
x 10
-3
date
d-
le
ss
Gyro 1. Fit Cgb 4
12/19
03/29
07/07
-1
0
1
2
3
4
5
6
x 10
-3
date
d-
le
ss
Gyro 1. Fit Cga 3
12/19
03/29
07/07
-6
-4
-2
0
2
4
x 10
-3
date
d-
le
ss
Gyro 1. Fit Cgb 3
a
1
b
1
a
3
b
3
b
4
b
2
a
2
a
4
Figure 3.
Gyro 1 Polhode Harmonics from Batch Analysis with Fits to
Îł
C
g
-model (dimensionless values).
.)
2
/
)
(
tan(
)
(
,...;
2
,
1
),
(
)
(
)
(
),
(
)
(
)
(
,
))
(
sin(
))
(
cos(
1
)
(
0
0
1
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
t
t
m
t
b
a
b
a
t
b
a
b
a
t
m
b
t
m
a
C
t
C
n
N
n
mn
mn
m
m
n
N
n
n
n
M
m
m
m
g
g
n
Îł
Îľ
Îľ
Îł
Îł
Îľ
Îł
Îł
=
=
ââ
â
â
ââ
â
â
=
ââ
â
â
ââ
â
â
ââ
â
â
ââ
â
â
=
ââ
â
â
ââ
â
â
â
âŹ
âŤ
âŠ
â¨
â§
ÎŚ
+
ÎŚ
+
=
+
=
=
=
â
â
â
blue
â measured
a
m
(t), b
m
(t)
red
- fit to
Îľ
(t) model
(2a)
(2b)
Heifetz et al.: The Gravity Probe B Data Analysis Filtering Approach 59
3.4
Gyroscope orientation profiles: moving window method
The main inputs of the 2
nd
Floor are the gyro orientation profiles from the 1
st
Floor, along with their
covariance matrix. However, the 1
st
Floor profiles are not smooth enough for the torque modeling.
Smoothing of the noisy 1-orbit batch orientation estimates is performed by means of a âmoving windowâ
concept.
A window of a length of several polhode periods moves along the time axis with one-orbit shifts. Within
each window the information obtained by the one-orbit batch processing is integrated: components of the
1
st
Floor state vector for each orbit are treated as the âmeasurementsâ of the window state vector
parameters, and the information matrices available for each orbit are summed up. To propagate the gyro
orientation within a window, we apply the misalignment torque model (see Sect. 4.1). If there is a jump in
one or both orientation components, it is modeled as a ârampâ. The output of the âmoving windowâ
estimator is a pair of smoothed profiles of gyro orientation (1 point per orbit) in the NS and EW
directions, along with their uncertainty at each point (see Fig. 4 and Fig.5). There is a clear linear trend in
the NS plot that should probably represent the geodetic effect, and about the right size, relative to the GR
prediction. However, the EW component is by far not a straight line: the small frame-dragging effect is
completely concealed by classical torque contributions.
Figure 4.
Gyroscope Orientation Time History (NS direction; units of the vertical axis are arc-seconds).
Figure 5.
Gyroscope Orientation Time History (EW direction). Red vertical lines correspond to the resonances
)
(
t
m
p
r
Ď
Ď
=
(units of the vertical axis are arc-seconds).
60 Gravity Probe B Science ResultsâNASA Final Report, December 2008
4 Floor 2: Misalignment and Roll-Resonance Torque
Modeling. Torque Compensation and Relativity Estimation
The âmoving windowâ smoothing allows one to observe the fine structure of the gyro orientation. It was
instrumental in understanding and modeling of the resonance phenomena, i.e. jumps in the orientation
profiles. Kolodziejczak (2007) identified many significant jumps with the so called
roll
-
resonance
torques which are not averaged out over the roll period around times,
t
m
,
when the constant roll frequency,
r
Ď
, coincides with a harmonic of the time-varying polhode frequency,
p
Ď
, i.e.,
)
(
m
p
r
t
m
Ď
Ď
=
[Everitt
et al. (2008), Keiser et al. (2009)]. A physical model of both the misalignment and roll-resonance torque
and its use in the 2
nd
Floor analysis is described in this section.
4.1
Misalignment and roll-resonance torque model
One of the assumptions in the early phases of the post-flight GP-B data analysis was that only the roll-
averaged patch effect torque should be considered. Patch effect torque is caused by the variation of
electrostatic potential on the rotor and housing surfaces and their relative motion. Discovery of the
resonance torques and subsequent derivation of the roll-resonance torque model [Keiser et al. (2009)]
showed that the non-roll-averaged torque should be taken into account, that is, should be modeled and
compensated.
As mentioned, the patch effect produces two torques that have to be taken into account in the data
reduction. The misalignment torque is proportional to the misalignment
Îź
, with the coefficient
determined by the patch patterns. There is also a contribution of zero order in
Îź
(not vanishing when
Îź
=0), which varies at the harmonics of roll; the coefficients involved in it are also determined by the
patches.
Due to the polhode motion of the rotor, all the torque coefficients are superpositions of multiple polhode
frequency harmonics. When one of the multiple frequencies coincides with the roll frequency, a non-
averaging d.c. contribution appears and affects the long-term gyroscope drift. In addition, there is some
evidence that the resonance torque does not completely average to zero over every roll period even
between the resonances. Therefore it should also be modeled continuously in the estimation.
The equation of motion of the GP-B gyroscope spin axis (
s
Ë
) in the presence of the misalignment and
roll-resonance torque is:
where
NS
r
and
EW
r
are the relativistic drift rates,
NS
Îź
and
EW
Îź
are the inertial components of the
misalignment,
Îź
r
, (see Fig.1), and
)
(
),
(
t
c
t
k
Âą
are the torque coefficients. Expressions in brackets
represent the non-roll-average, roll-resonance torque model, where
r
ÎŚ
is the roll phase, and
θ
is the
misalignment phase,
)
/
arctan(
NS
EW
Îź
Îź
θ
=
. The misalignment phase can be used as a pre-calculated
function of time (as done in the analyses of data presented below), or set to be a constant. There is some
question whether
θ
has the form used here, or the constant value
2
/
Ď
θ
=
as in the resonance torque
model [Keiser et al. (2009)].
)],
cos(
)
(
)
sin(
)
(
[
)
(
)]
sin(
)
(
)
cos(
)
(
[
)
(
r
r
NS
EW
EW
r
r
EW
NS
NS
t
c
t
c
t
k
r
dt
ds
t
c
t
c
t
k
r
dt
ds
ÎŚ
+
+
ÎŚ
+
+
â
=
ÎŚ
+
â
ÎŚ
+
+
+
=
+
â
+
â
θ
θ
Îź
θ
θ
Îź
(3)
Heifetz et al.: The Gravity Probe B Data Analysis Filtering Approach 61
The torque coefficients
)
(
),
(
t
c
t
k
Âą
related to the patch distribution can be represented in exactly the
same way as the scale factor
g
C
in Eqs. (2a), (2b), namely:
(4a)
(4b)
Here
0
ÎŚ
and
0
Îł
are the polhode phase and polhode angle calculated for a symmetric gyroscope, and
c(t)
stands for any of the torque coefficients. Thus
),
(
),
(
)
(
1
1
1
t
c
t
k
t
c
m
m
m
Âą
Âą
=
,
,
1
1
1
Âą
Âą
=
mn
mn
mn
c
k
c
and
),
(
),
(
)
(
2
2
2
t
c
t
k
t
c
m
m
m
Âą
Âą
=
,
,
2
2
2
Âą
Âą
=
mn
mn
mn
c
k
c
respectively. The only difference
between the formulas (2b) and (4b) is in the meaning of the states to be estimated,
}
,
{
mn
mn
b
a
in (2b), and
}
,
{
2
1
mn
mn
c
c
in (4b). The values of the first set are determined by the trapped field distribution on the
rotor surface, while in the second case they depend on the electrostatic patch patterns and the inertial
asymmetry of the rotor.
The torque model (3), (4a), (4b) exploits the combinations of phases
p
r
m
m
ÎŚ
ÎŚ
=
ÎŚ
m
m
, and produces
shifts in the orientations
NS
s
and
EW
s
at the resonance times when
)
(
t
m
â
ÎŚ
passes through zero, which is
equivalent to the above resonance condition
)
(
t
m
p
r
Ď
Ď
=
. Naturally, the following conceptual questions
arise: Can the torque model (3) - (4) explain the observed multiple jumps in the gyroscope orientation
profiles? Could all the significant torque contributions be estimated using this model and subtracted from
the gyro orientation profiles, providing thus a reconstruction of the gyroscope relativistic motion and
measurement of its rate? As shown below, both answers are affirmative.
4.2
Elimination of torque contributions from the 1st Floor
orientation profiles: Proof of concept
The output of the 1
st
Floor includes, in general, the multidimensional state vector (gyro orientation, scale
factor coefficients, roll phase offset, SQUID bias components, etc.) and the full information matrix. To
address the above questions, we consider only the orientation profiles and their uncertainties as inputs to
the 2
nd
Floor processing. The Kalman filter is a natural tool for solving this simplified estimation
problem: orientation profiles give the measurement of the orientations
NS
s
and
EW
s
at the discrete, 1 point
per orbit, known time points. Their behavior,
)
(
t
s
NS
and
)
(
t
s
EW
, is described by the dynamic model (3),
(4a), (4b). One complication should be immediately addressed, namely, science data segmentation,
according to Table 1. The segments are separated by the S/C anomalies. Some of them, such as, for
instance, solar flare between segments 5 and 6, could change the patch patterns, and thus the values of the
torque coefficients. To overcome this obstacle, the dynamic estimator was designed in the following
framework (Segments 5, 6, and 9 were analyzed for gyros 1, 2, 4; the results for gyro 3 are from
Segments 5 and 6 only).
)],
(
sin
)
(
)
(
cos
)
(
[
)
(
0
0
2
0
~
0
0
1
t
m
c
t
m
c
t
c
m
M
m
m
ÎŚ
+
ÎŚ
=
â
=
Îł
Îł
.)
2
/
)
(
(
tan
)
(
,...;
2
,1
),
(
)
(
)
(
);
(
)
(
)
(
0
0
1
0
~
1
2
1
0
2
0
1
0
~
0
20
10
0
20
0
10
0
t
t
m
t
c
c
c
c
t
c
c
c
c
n
N
n
mn
mn
m
m
n
N
n
n
n
n
Îł
Îľ
Îľ
Îł
Îł
Îľ
Îł
Îł
=
=
ââ
â
â
ââ
â
â
=
ââ
â
â
ââ
â
â
ââ
â
â
ââ
â
â
=
ââ
â
â
ââ
â
â
+
=
=
â
â
62 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Step1. Independent Kalman filter/smoother for each segment and each gyroscope.
Introduce the state vector
{
}
{
}
[
]
mn
mn
mn
mn
EW
EW
NS
k
k
c
c
r
r
s
s
x
NS
2
1
2
1
,
,
,
,
,
,
,
Âą
Âą
=
, where the first two components
are dynamically changing states according to the model (3), and all other components are constants. The
time-varying roll-resonance torque coefficients
)
(
t
k
,
)
(
),
(
t
c
t
c
+
â
are represented by the model (4.a),
(4.b); the parameters of this model,
mn
mn
k
k
c
c
mn
mn
2
1
,
,
,
2
1
Âą
Âą
, remain constant within a given segment.
A Kalman filter technique was implemented in the form of a smoother with the forward-backward
propagation-update structure [Bryson (2002), Bar-Shalom (1998)]. The observability analysis showed that
the roll-resonance torque coefficients in the vicinity of the corresponding resonances are strongly
observable, and their estimates are statistically significant. The outputs of the Step 1 estimators are,
specifically, the estimates of the torque coefficients and the modeled orientation profiles: the solution of
the equations (3) with the estimated torque coefficient values. An example of the modeled orientation
profiles, that has the same shape as the input, is presented in Fig.6.
Step 2. Reconstruction of ârelativistic trajectoryâ: subtracting torque contributions.
We can combine the segments by subtracting the torque contributions estimated for each of the segments
in Step 1 from the input orientation profiles. The resulting curve is the reconstructed ârelativistic
trajectoryâ; it consists of the segment pieces that have the same ârelativisticâ slope. The reconstructed
trajectory can then be fit to the same-slope line over each segment; thus the state vector includes a
common slope and individual y-intercepts.
Heifetz et al.: The Gravity Probe B Data Analysis Filtering Approach 63
Figure 6
. Gyroscopes 4 (top) and 2 (bottom): a) East-West orientation profile; b) modeled orientation profile; c)
reconstructed ârelativisticâ trajectory.
The estimate of the slope gives the value of the corresponding relativistic drift rate. The reconstructed
trajectory noise is larger than the orientation profile noise by the known amount defined by the
uncertainty in the subtracted torque contributions.
The described two-step approach can be used not only for one particular gyroscope, but also for any
combination of the gyroscope signals. Step 2 easily allows for fitting the pieces of the reconstructed
trajectories of different gyroscopes over different segments, since they all have the same relativistic drift
rate (slope). Reconstructed trajectories (corrected for the segment shifts) of the gyroscopes 4 (10
resonances) and 2 (73 resonances) are presented in Fig. 6. Both the NS and EW reconstructed trajectories
now contain a clear linear trend,
the EW curve for the first time ever
.
Estimates of the relativistic drift rates can also be obtained directly from the output of Step 1 estimators.
Each individual estimator results in the state vector and information matrix whose dimensions differ from
segment to segment owing to the different number of resonances and torque coefficients. Nevertheless,
they also have a common part: relativistic drift rates. The integration of the information from all one-
segment estimators can be done by creating the augmented state vector and information matrix. The
augmented state vector includes the relativistic parameters and all torque coefficients from all segments,
and the augmented information matrix consists of blocks of the individual information matrices. The
combined covariance matrix is the inverse of the augmented information matrix, and the combined state
vector estimate is a linear combination of the individual estimates with weighting coefficients, which
depend on the combined covariance matrix and the individual information matrices. These two described
methods of relativity parameters estimation are statistically equivalent.
64 Gravity Probe B Science ResultsâNASA Final Report, December 2008
4.3
Two-Floor Analysis Results
The estimates of the relativistic drift rates with 50% - probability error ellipses obtained by the second
method described above are shown in Fig.7. The âindividualâ gyro estimates are consistent with each
other (the ellipses overlap). The error bars of the estimates given in Fig. 7 are formal statistical 50%
probability errors.
The data analysis using the methods outlined above is still preliminary since a number of factors have not
been fully investigated:
1)
More data segments for gyroscope analysis
2)
Variation of the number of parameters used to characterize the torque coefficients;
3)
Unmodeled error in the measurement of the telescope pointing during GSI;
4)
Use of the maximum information matrix to determine the error ellipse;
5)
Inclusion of systematic errors, etc.. (The accompanying paper Muhlfelder (2009) discusses
approaches to establishing systematic errors bounds.)
Nevertheless, the preliminary results given in Fig. 7 are, of course, subject to the limitations just
described. There are two important observations that can be made in this data analysis example: (1) the
residual errors are substantially smaller than found with earlier data analysis methods, and (2) the
relativity estimates for the individual gyroscopes are statistically consistent with each other. The
combined result for all four gyroscopes is in agreement with the GR values. (GR predictions for the real
GP-B orbit, taking into account the GS proper motion, the solar geodetic effect, and the EW projection of
the Earth geodetic effect, are -6571 marcs/year in the NS direction, and -75 marcs/year in the EW
direction).
Figure7
. Preliminary relativity estimates. Gyroscopes 1, 2, 4: segments 5, 6, and 9, (gyro 3: segments 5,6). Do not
include systematic error or a model sensitivity analysis.
Heifetz et al.: The Gravity Probe B Data Analysis Filtering Approach 65
5 Data Analysis Expansion: Two-second Filter
The accuracy of the results of the above two-floor structure is badly limited by the orientation profile time
step, which is one data point per orbit (1
st
Floor data compression). Effective direct modeling of the roll-
resonance torques and their compensation in the gyro orientation profiles, implemented as a proof of
concept, emphasized the existence of torques at the roll frequency corresponding to ~1 min roll period,
apparently, much smaller than the duration of one orbit (~97 min). Consequently the current two-floor
estimation structure does not fully address the roll-resonance torque. To improve the relativity
measurements, it is desirable to replace it with a single-floor estimator/filter, with the data sampling rate
much higher than the roll frequency, so that continuous dynamical modeling of all the pertinent features
takes place.
We are currently working on the design and implementation of the single-floor âtwo-second filterâ where
the SQUID and telescope science signals, sampled at the rate of 2 seconds, are processed according to the
readout model (1), scale factor model (2a)-(2b)
simultaneously
with the patch effect torque model (3),
(4a), (4b). The state vector includes the relativity parameters, scale factor coefficients, and all the
parameters currently estimated at the 1
st
Floor, as well as the torque coefficients, which are so far
estimated separately at the 2
nd
Floor. For the first time, S/C pointing will be determined consistently with
all other data. During a GSI mode it will be estimated and iteratively updated using the available four
SQUID signals, and the estimates of all the needed parameters obtained from the analysis of the previous
GSV mode. Information about the GSI pointing is needed for the gyro orientation propagation due to
continuously acting torques.
This enhanced two-second filter will provide accurate modeling of the underlying physical phenomena
implemented in a flexible, module-based structure. It will permit fast incorporation and evaluation of
instrument models, and the check of sensitivity to variations in the source data and parameter sets. The
architecture of the two-second filter incorporates the âone systemâ approach: all the GP-B instruments
(gyroscopes, telescope, spacecraft, etc.) should be treated as a sigle GP-B system. It includes the model
identification and verification on the sensor-to-system level (gyroscope, SQUID, telescope, spacecraft,
etc.). It also encompasses the development, verification and testing of filtering algorithms based on
distributed computing and parallel processing. The enhanced direct torque model toolset brings in
important benefits to the overall GP-B data analysis picture and increases confidence in final results by
testing for consistency during the analysis.
6 Conclusion
A
â
simpleâ strategy of the GP-B data analysis has evolved in the complex two-floor structure after on-
orbit discoveries of the changes in the rotorâs polhode period and path, and of patch effect torques. Direct
modeling of the readout scale factor at the 1
st
Floor, and the misalignment and roll-resonance torque
modeling at the 2
nd
Floor, allowed us to separate the relativistic drift from the drift induced by classical
torques. A cascade of four interconnected estimators applied to the GP-B science data has demonstrated a
consistent determination of the geodetic and frame-dragging effects for all GP-B gyroscopes, as well as
fidelity of the physical models used.
There is a clear way to improve the accuracy of GP-B results. Existence of a roll component in the patch
effect torque suggests the development of a one-level estimator (â2-second filterâ), with a sampling data
rate much higher than the roll period. This filter will simultaneously exploit all the models currently used
in the estimators of the 1
st
and 2
nd
Floors. Though the number of estimated parameters in the new filter
will be much larger than the one in the reduced 1
st
and 2
nd
Floor filters, the new estimator will be free of
the shortcomings associated with the current two-floor data analysis structure.
66 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Acknowledgments.
The authors would like to thank D. Bartlett, B.Clarke, J.Conklin, D.DeBra, F.Everitt, G.Keiser,
A.Krechetov, J.Li, I.Mandel, B.Muhlfelder, B.Parkinson, R.Reasenberg, J.Ries, P.Saulson, J.Turneaure,
and P.Worden for numerous discussions, insights and support in our data analysis journey.
References
Y. Bar-Shalom, X-R. Li, T. Kirubarajan,
Estimation with Applications to Tracking and Navigation
(Wiley, New York, 2001).
G.J. Bierman,
Factorization Methods for Discrete Sequential Estimation
(Dover Publications, New York,
2006).
A.E. Bryson,
Applied linear optimal control: examples and algorithms
(Cambridge University Press,
New York, 2002).
J.W.Conklin,
Estimation of the Mass Center and Dynamics of a Spherical Test Mass for Gravitational
Reference Sensors,
Ph.D. Thesis, Aero/Astro, Stanford University, (Stanford, 2008).
C.W.F. Everitt, M. Adams, W. Bencze, et al., Class. Quantum Gravity, 25 (2008) 114002 (11p).
C.W.F. Everitt, M. Adams, W. Bencze, et al., Space Science Reviews, (2009, Chapter 1 of this report).
G.T.Haupt,
Development and Experimental verification of a nonlinear data reduction algorithm for the
Gravity Probe B Relativity Mission,
Ph.D. Thesis, Aero/Astro, Stanford University, (Stanford, 1996).
M.I. Heifetz, G.M. Keiser, A.S. Silbergleit, in
Nonlinear Gravitodynamics,
Eds. R.J.Ruffini, C.
Sigismondi
(World Scientific, New Jersey-London-Singapore-Hong Kong, 2003).
T. Kailath, A.H. Sayed, B. Hassibi,
Linear estimation
(Prentice Hall, Upper Saddle River, New Jersey,
2000).
G.M. Keiser, J. Kolodziejczak, A.S. Silbergleit, Space Science Reviews, (2009, Chapter 3 of this report).
J.Kolodziejczak,
Correlation between times of observed gyro orientation changes and high polhode
harmonic - roll frequency resonances,
Science Advisory Committee to GP-B, Meeting 16, Stanford
University, (March 23-24, 2007).
D.D. McCarthy,
IERS Conventions (1996),
(U.S. Naval Observatory, 1996).
B.Muhlfelder, M. Adams, B. Clarke, et al., Space Science Reviews, (2009, Chapter 5 of this report).
P.K. Seidelmann,
Explanatory Supplement to the Astronomical Almanac
, University Science Books,
Mill Valley, California, 1992.
A.S. Silbergleit, M.I. Heifetz, A.S. Krechetov,
Model of Starlight Deflection and Parallax for the GP-B
Data Reduction,
Gravity Probe B document S0393 (Stanford University, Stanford, 2001)
A.S. Silbergleit,
On the linearity of the Gravity Probe B Measurement Model,
Gravity Probe B document
S0632, Rev. C (Stanford University, Stanford, 2006).
A.S. Silbergleit,
GP-B Science Signal Models and Their Analysis,
Gravity Probe B document S0996
(Stanford University, Stanford, 2009).
A.S. Silbergleit, J. Conklin, D. DeBra, et al., Space Science Reviews, (2009, Chapter 2 of this report).
J.P. Turneaure, C.W.F. Everitt, B.W. Parkinson, et al., Advances In Space Research 32 (7), 1387 (2003).
67
5
GP-b s
ystematic
e
rror
Barry Muhlfelder et al.
68 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Muhlfelder: GP-B Systematic Error 69
GP-B Systematic Error
B. Muhlfelder
1
, M. Adams, B. Clarke, G.M. Keiser, J. Kolodziejczak
2
, J. Li, J.M. Lockhart, P. Worden
Abstract:
We have evaluated the systematic error in the GP-B experiment using five different approaches and
estimated the contribution of more than 100 specific sources of the error. The systematic effects we
consider include those due to gyroscope torques, gyroscope readout, telescope readout, and guide star
proper motion. Effects with an estimated impact on the experiment error larger than 1 mas/yr are
discussed in detail. Examples of analyses that bound other sources to less than 1 mas/yr are included to
show the range of techniques employed to perform this work. We describe the remaining tasks to
complete the systematic error analysis and estimate the total experiment uncertainty.
1 Introduction
The bounding of systematic error or uncertainty is an essential part of every high accuracy physics
experiment. An early discussion of GP-B systematic uncertainty may be found in a report that has come
to be known as âThe Green Bookâ [Everitt et al. (1980)]. Many GP-B ground-based tests and analyses
were described there, and a theoretical framework was presented, providing an early quantification of the
experimentâs systematic uncertainty.
The GP-B experiment was designed to provide an accurate measurement of non-Newtonian (general
relativistic) gyroscope precession without modeling systematic effects. A detailed error tree, containing
both stochastic and systematic sources for every gyro, allowed trade studies to be performed for various
design parameters. Prior to launch in April 2004, the experiment uncertainty was expected to be limited
by stochastic noise dominated by the noise in the SQUID-based London Moment gyroscope readout
system. Based upon test and analysis techniques, the expectation was that the impact of systematics on
the overall experiment accuracy would be small, and that without modeling systematic effects, the total
experiment uncertainty would be less than 0.5 mas/yr.
During the first two weeks of on-orbit GP-B operations, the SQUID noise calibration was successfully
performed. The results are shown in Figure 1. Using a systematics-free measurement model, the
experiment uncertainty based upon this noise level is ~ 0.2 mas/yr for each gyroscope, or ~ 0.1 mas/yr
when combining the results from all 4 gyroscopes. Although we expected that there would be little
impact from systematic effects, after completing the SQUID noise calibration, we began an exhaustive set
of calibrations designed to test the pre-launch expectation.
1
Muhlfelder (corresponding author)
Gravity Probe B, HEPL, Stanford University, Stanford, CA 94305-4085 USA
Tel: 1-650-725-4125
Fax: 1-650-725-8312
E-mail: barry@relgyro.stanford.edu
2
NASA/Marshall Space Flight Center
70 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Requirement
Figure 1.
Low frequency SQUID noise
The on-orbit mission plan comprised three phases: (1) Initialization, (2) Science, and (3) Calibration.
These phases were designed to perform three tasks: configure the hardware for science during
Initialization, collect the science data during Science, and calibrate systematic effects during all phases.
Keiser (2009) has described the numerous SQUID and telescope readout operations, acquisition of the
guide star, and gyroscope suspension and spin up. All the science hardware, including each of the four
gyroscopes and their readout systems, the telescope and its eight detector channels, and many
instrumentation channels functioned successfully. Initialization of the experiment was completed in
August 2004. By comparison, few operations were required during the Science phase. The goal of this
phase was to collect a year of science data with no changes to the configuration of the science instrument
and minimum interference of any kind. More than 1 Terabyte of science data was collected prior to the
completion of the science phase on August 15, 2005. The final on-orbit operations were performed
during the 6½ week long Calibration phase.
The full set of GP-B calibrations probes more than 100 physical sources including guide star proper
motion, telescope readout, and gyroscope readout and classical torques. The uncertainty in the proper
motion is known from published sources to better than ~ 1 mas/yr. To date, all GP-B analyses have been
performed with information available from these sources. Recent measurements [Shapiro et al. (2007)]
will provide even tighter limits on proper motion uncertainty, however, by agreement, the actual proper
motion values have been withheld from the GP-B data analysis team to allow for future comparison and
incorporation. In this paper we focus on readout issues and gyroscope torques discovered in the course of
calibration experiments and data analysis. Both of these effects are much larger in magnitude than the
known guide star proper motion uncertainty.
Analyses of data collected on-orbit demonstrate that all but three GP-B systematic effects are small.
These three behaviors are damping of the gyroscope polhode motion and two related Newtonian
gyroscope torques. In accompanying papers, Silbergleit et al. (2009) describe a technique to model the
effect of the gyroscopeâs changing polhode on the gyroscope scale factor and Keiser et al. (2009) describe
Muhlfelder: GP-B Systematic Error 71
the Newtonian torques generated by the interaction of patch potentials on the rotor and housing surfaces.
Assuming only gyroscope spin and roll averaging, two distinct manifestations of these patches were
found: (1) torques that are perpendicular to and proportional to the gyroscope misalignment and (2)
torques that occur at the difference frequency of harmonics of polhode and the spacecraft roll frequency.
Left unmodeled, these effects would cause experiment error up to 1 arcsec/yr. Contrary to the initial plan,
we had to model these three effects in the data analysis to achieve the best accuracy for the experiment.
These modeling efforts, begun in 2004 and continuing to the present, have resulted in a dramatically
reduced systematic error. In the discussion here we focus on the residual error associated with imperfect
modeling of these effects. Before describing these and other smaller effects, we give the methodology for
evaluating the total systematic experiment uncertainty.
2 Approaches to estimating experiment uncertainty
The study of individual systematic effects involves three steps: identification, modeling (if required), and
bounding. The identification of a systematic effect connects a physical mechanism disturbing the
measurement to the experiment. Modeling is only required for large effects, specifically those associated
with patch potentials. Physics-based models governing the patch effect are incorporated into the science
data analysis to separate the relativistic drift from that induced by patch effect torques. Systematic effects
that can be bounded below 1 mas/yr do not require modeling.
Five approaches are used to bound the total GP-B experiment uncertainty:
1.
Bottom-up analyses
2.
Sensitivity tests
3.
Gyro-to-gyro comparisons
4.
Data sampling tests
5.
Comparison of results of independent analysis methods
2.1
Bottom-up Analyses
These analyses give the uncertainty associated with individual sources and are based upon physical
models of the hardware with parameter values determined from ground and/or flight measurement. For
instance, a bottom-up analysis gives the gyroscope precession rate resulting from the support force acting
on the gyroscope mass unbalance vector. The physical model here is the mass unbalance torque causing a
change in the angular momentum of the spinning mass. This model, combined with the known
parameters, bounds the effect. The parameters are the support force (measured on-orbit), the mass
unbalance (measured on the ground and refined on-orbit), the angular velocity (measured on-orbit), and
the gyroscope moment of inertia. Using these known parameter values gives a Newtonian drift rate
below 0.08 mas/yr.
2.2
Sensitivity tests
The assessment of the impact of an individual systematic effect on experiment uncertainty is determined
by using the science analysis algorithms to determine the sensitivity to a particular effect. One sensitivity
test method is to propagate the numerical uncertainty in a system parameter through the modeling
software. The resulting variation in the relativity results gives the experiment uncertainty due to this
parameter.
72 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Example: Timing error. The gyroscope orientation evolves in time, and therefore the model to describe
this evolution is susceptible to timing error. To assess the impact of timing error on the experiment, we
find the science result using the time stamp provided by the onboard clock for each data point. The
uncertainty in our knowledge of the timing is determined by performing an independent analysis of on-
orbit timing data. The time stamp for each data point is modified by this uncertainty and the analysis is
rerun to give a new science result. The change in it is a measure of the timing systematic error. This
analysis bounds the impact of the timing error on the relativity determination to better than 0.5 mas/yr.
2.3
Gyro-to-gyro variation
The range of relativistic drift rate estimates for the four gyroscopes gives a measure of the experiment
uncertainty due to stochastic and systematic sources, but it does not include common-mode effects. The
dominant stochastic error results from noise in the gyro readout system. This noise directly leads to gyro-
to-gyro differences. Gyroscope systematic errors result from modeling deficiencies in the gyroscope
readout system and gyroscope torque modeling. Since the gyroscopes are independent sensors, we expect
gyro-to-gyro science differences.
2.4
Data selection tests.
Ideally, the science result should be independent of which data are selected for analysis. Difference in the
results for different data sets provides a measure of experiment uncertainty. For example, the yearlong
mission was interrupted nine times due to spacecraft anomalies resulting in 10 segments of science data.
Variation in the experiment result based on the analysis of various combinations of these segments
provides a measure of the experiment uncertainty.
2.5
Two independent analysis methods and teams
Two independent teams, using significantly different modeling approaches, have performed full science
analyses. Although these teams use the same science data, and address the same physical effects, their
results are based upon separate analysis algorithms. The most significant difference in their analysis
methodologies is in the treatment of the misalignment torque. One approach takes a torque-free
projection of the gyroscope motion whereas the other explicitly models this effect. The two teams also
model the resonance torque and the evolving gyroscope scale factor with differing techniques described
below. Agreement in science results using these two modeling approaches increases our confidence in
the result.
3 Initial models of scale factor, torques, and the associated
early science results
The modeling of the evolving gyroscope scale factor and the patch effect torques has dramatically
improved during the past 4½ years. Our progress has resulted from insights into the underlying physics
and from key advances in essential supporting analyses.
Prior to launch, we had expected the scale factor to be periodic at polhode frequency and patch effect
torques to be negligible. Early in the Science phase, it became clear that the gyroscope polhode path and
frequency were changing with time due to polhode damping, thereby resulting in a slowly evolving
gyroscope scale factor. The first and simplest technique provided an estimate of the scale factor profile
for multi-day-long batches of data. This approach treated the scale factor profile for each batch of data as
a Fourier expansion of polhode harmonics. Although over the past four years we have developed
Muhlfelder: GP-B Systematic Error 73
significantly more sophisticated techniques, this first method was used during the science and calibration
phases of the mission to track the orientation of the gyroscopes. The gyroscope orientation information
was essential to the discovery of the patch effect torques. This connection of one modeling advance
leading to another has been repeatedly seen in our work and underlies our progress over the past several
years.
Large misalignment torques were discovered during the calibration phase of the mission. The
geometrical nature of these torques was established through a series of tests in which the space vehicle
was pointed to neighboring stars, thereby dramatically increasing the misalignment of the gyroscope spin
axis relative to the spacecraft roll axis [Keiser et al., 2009]. The first misalignment torque models either
explicitly included a drift term perpendicular to the measured misalignment or took a torque-free
projection of the gyroscope drift.
The resonance torque was discovered after completing the science phase. It was observed that when a
harmonic of polhode frequency coincided with the spacecraft roll frequency, the gyroscope orientation
would often shift. In some cases this shift, typically occurring in one or two days, was as large as 0.15
arcsec. The original method to model this behavior was to excise the data near the shift. Although the
underlying physical explanation for these shifts was unknown, this approach significantly improved the
fit of the data to the model.
In April of 2007, using the techniques as described above, we reported an experiment uncertainty of 97
mas/yr. This result was based upon a series of sensitivity analyses. In these analyses we probed the
sensitivity of the result to scale factor and torque modeling deficiencies. Scale factor modeling
deficiencies were bounded by varying the number of polhode harmonic terms included in the model and
by changing the batch length used to determine these terms. Misalignment torque error due to
misalignment uncertainty was derived from two independent analyses. Figure 2 shows the result of these
two analyses for gyroscope 1, plotted here as gyroscope orientation (the orientation is equal to the vehicle
pointing minus the misalignment). The difference in results gives an uncertainty of 30 mas. Sensitivity
of the science result to this level of uncertainty was determined by modifying the misalignment by this
amount. This change to the misalignment caused only a small change (4 mas/yr) shift to the science
result, suggesting some other source for the leading cause of error.
Figure 2
. Gyroscope 1 WE orientation using two different methods from an analysis performed in early 2007.
Differences give a measure of the uncertainty.
74 Gravity Probe B Science ResultsâNASA Final Report, December 2008
4 Model Improvements
We found that the largest contribution to experiment uncertainty in 2007 was caused by resonance torque
mis-modeling, although we lacked a detailed understanding of the modeling deficiency. The science
analysis method used in 2007 excised data from the analysis for each resonance when the difference
between an harmonic of polhode phase and the roll phase was less than a specified value. For those orbits
with phase difference greater than the specified value, the resonance torque was ignored. Sensitivity tests
were performed by varying the resonance width as defined by this phase difference. This variation in
phase difference criteria led to larger than a 50 mas/yr variation in the science result. A subsequent
method used linear ramps to model the NS and WE orientations during resonances. Although tests using
this model demonstrated approximately three times lower sensitivity to the resonance torque, both this
approach and the previous one were simplifications of the true physical phenomenon. These methods
assume that the resonance torque vanishes before and after each resonance. A study of flight data now
gives insight into the potential limitations of these approaches [Keiser et al. (2009)]. Oscillations, as
predicted from the full torque model and as observed in the flight data, before and after each resonance
are ignored in the simplified models. Our most recent efforts make explicit use of the full model. Two
comments are in order here: First, the similarity between the data and the full model increases our
confidence in our understanding of the underlying physics, and second, this model improves the
numerical fit of the model to the data, thereby reducing the residuals. Initial science results making use of
the improved model are encouraging [Heifetz et al. (2009)].
Use of the full resonance torque model requires accurate knowledge of the polhode phase, only recently
available. Given that the observed resonances occurred for high harmonics of the polhode frequency, and
that the uncertainty in the polhode phase scales with harmonic number, the need for accurate phase
information is clear. Trapped Flux Mapping (TFM) now provides a continuous and accurate estimate of
polhode phase [Silbergleit et al. (2009)]. The uncertainty in these estimates has improved dramatically
over the past year with a current uncertainty of < 0.5.deg. Sensitivity tests are in progress to determine
the impact of this uncertainty on the overall experiment uncertainty.
5 Sensitivity Tests
The improvements we described above to the resonance torque modeling were motivated by initial high
sensitivity. These recent improvements suggest that we reexamine the other modeling sensitivities to
determine an updated total uncertainty. Table 1 gives a list of the most important tests. The GP-B guide
star, IM Pegasi, is known to be a binary star with a ~25 day orbital period. Adding a model of this motion
to the analysis caused less than a 5 mas/yr
change to the science result. This finding not only
demonstrated insensitivity to this motion, but also has allowed a determination of the starâs orbital
parameters. These parameters will be made public once the guide star results from the SAO are released.
A sensitivity test comparing harmonic expansion,
Îł
C
g
[Heifetz (2009)] and TFM analyses provide a
determination of our experiment uncertainty due to readout scale factor error. This is a particularly
interesting test, making use of the low frequency gyroscope signal in two cases (harmonic expansion and
Îł
C
g
) and the high frequency gyroscope signal in the other (TFM). Use of the TFM results also simplifies
the science analysis from a nonlinear to a linear fit and thus provides a check of sophisticated non-linear
estimation techniques. Tests were performed to probe the sensitivity of the science result to misalignment
and resonance torque modeling.
Muhlfelder: GP-B Systematic Error 75
Sensitivity Parameter
Parameter Range
Guide star orbital motion
0 to fitted
Resonance torque
Continuous fit vs. fixed width (0.75, 1.5, 2)
Scale factor modulation
Harmonic expansion,
Îł
C
g
, TFM
Polhode phase
-0.5° to 0.5°
Misalignment torque averaging window
15 â 25 orbits
Misalignment torque uncertainty
-0.1 to 0.1 radians
Table 1
. Sensitivity tests and the associated range of parameter values.
A preliminary estimate of the uncertainty in science result for gyroscope 1, based upon all of these
sensitivity analyses gives a preliminary uncertainty of 30 mas/yr. Gyro-to-gyro comparisons by Heifetz
et al. (2009) suggest a still lower bound. We continue to work to reduce these limits.
Following launch, analysis of the flight data has confirmed many potential systematic effects are less than
1 mas/yr. Currently these effects do not require modeling. In the future these resolved effects may
require re-examination as the total experiment error decreases.
Rather than provide a long list of resolved systematic effects, it is instructive to describe several
representative analyses.
5.1
Roll phase uncertainty coupling to science result
Roll phase uncertainty or error cause the science signal in the NS direction to be rotated into the WE
direction and the science signal in the WE direction to be rotated into the NS direction. To limit the
impact on the science result to 1 mas/yr requires 10 arcsec roll phase knowledge.
Figure. 3.
Spacecraft roll phase uncertainty vs. time.
Arc-sec
240 290 340 390 440 490 540
Days from 01-Jan-2004
6 arcsec RMS variation
76 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Figure 3 shows the roll phase uncertainty for the duration of the mission. The roll phase and roll phase
uncertainty is determined using a primary and backup CCD-based roll star telescopes. As the spacecraft
rolls about the line of sight to the guide star, other stars 30-40 degrees from the roll axis are tracked as
they rotate through each of the CCDâs field of view. The uncertainty in the roll phase is found by
comparing the primary and backup roll phase measurements.
5.2
Gyroscope position coupling to the science result
The gyroscope readout signal will be modified by the gyroscopeâs position change, if the position couples
to the readout signal. We know that a gravity gradient force acts on the gyroscopeâs position at roll +/-
twice orbit frequency. If the coupling exists, the gravity gradient force will cause a readout signal at roll
+/- twice orbit. This signal will mix with and modify the signal at near roll +/- orbit. The mixing is due
to the use of only the guide-star-visible portion of each aberration cycle to determine the gyroscope scale
factor. Figure 4 shows the results of two analyses for gyroscope 2. The signal that is out-of-phase with
orbital aberration should be very small. Confirming a small out-of-phase term increases confidence in the
scale factor as determined from the in-phase term. When the analysis is done without including the roll
+/- twice orbit term, a signal of 5 mas amplitude is observed out-of-phase at roll +/- orbit. The switch of
the sign in the middle of June occurred when the drag free sensor was switched from gyroscope 1 to
gyroscope 3. When the model includes roll +/- twice orbit, the estimated out-of-phase signal nearly
vanishes. This analysis allows us to place a 1 mas/yr limit on the experiment uncertainty associated with
the gyroscope 2 position. Analyses for the other gyroscopes demonstrated even weaker coupling than that
found for gyroscope 2.
Figure 4
. Out-of-phase orbital aberration term estimates. Cosine and Sine terms of out-of-phase orbital aberration.
2
Ă
orbital not modeled: Black/Red. 2
Ă
orbital modeled: Blue/Green.
Muhlfelder: GP-B Systematic Error 77
Figure 3 shows the roll phase uncertainty for the duration of the mission. The roll phase and roll phase
uncertainty is determined using a primary and backup CCD-based roll star telescopes. As the spacecraft
rolls about the line of sight to the guide star, other stars 30-40 degrees from the roll axis are tracked as
they rotate through each of the CCDâs field of view. The uncertainty in the roll phase is found by
comparing the primary and backup roll phase measurements.
5.2
Gyroscope position coupling to the science result
The gyroscope readout signal will be modified by the gyroscopeâs position change, if the position couples
to the readout signal. We know that a gravity gradient force acts on the gyroscopeâs position at roll +/-
twice orbit frequency. If the coupling exists, the gravity gradient force will cause a readout signal at roll
+/- twice orbit. This signal will mix with and modify the signal at near roll +/- orbit. The mixing is due
to the use of only the guide-star-visible portion of each aberration cycle to determine the gyroscope scale
factor. Figure 4 shows the results of two analyses for gyroscope 2. The signal that is out-of-phase with
orbital aberration should be very small. Confirming a small out-of-phase term increases confidence in the
scale factor as determined from the in-phase term. When the analysis is done without including the roll
+/- twice orbit term, a signal of 5 mas amplitude is observed out-of-phase at roll +/- orbit. The switch of
the sign in the middle of June occurred when the drag free sensor was switched from gyroscope 1 to
gyroscope 3. When the model includes roll +/- twice orbit, the estimated out-of-phase signal nearly
vanishes. This analysis allows us to place a 1 mas/yr limit on the experiment uncertainty associated with
the gyroscope 2 position. Analyses for the other gyroscopes demonstrated even weaker coupling than that
found for gyroscope 2.
Figure 4
. Out-of-phase orbital aberration term estimates. Cosine and Sine terms of out-of-phase orbital aberration.
2
Ă
orbital not modeled: Black/Red. 2
Ă
orbital modeled: Blue/Green.
5.3
Electromagnetic Compatibility (EMC)
The gyroscope readout system is based on the detection of magnetic fields of 10
-16
T near the satellite roll
frequency. Extreme care was taken in the design and operation of the science instrument to limit spurious
magnetic signals. Although these efforts were largely successful, one extraneous signal source, the
Experiment Control Unit (ECU) was identified during science data taking. Figure 5 shows approximately
1 minute of the science signal for gyroscope 1. The raw time series data profile is shown in black; the
processed data profile with the ECU signal removed is shown in red. Fortunately, as is the case in the
Figure, the science signal at 12.9 mHz is far removed from the ECU signal frequency more than 90% of
time. The readout data at these times are used in the science analysis. When the ECU signal frequency is
near the science signal frequency, we excise the data. Sensitivity tests show negligible variation in the
science result for varied amounts of excised data. To facilitate these and other tests the science database
incorporates a data-grading tool that allows consistent, convenient, unbiased testing.
Figure 5
. Gyroscope 1 low frequency science data
5.4
SQUID non-linearity
One area of on-going investigation is the small SQUID readout non-linearity. Thus far, analyses show
evidence for such a non-linearity although we do not yet have sufficient understanding to incorporate a
specific model into the science analysis. Prior to launch the SQUID was specified to have a non-linearity
of less than 0.00015 relative to 10 volts. Although this specification was successfully confirmed, tests on
the ground and later on-orbit revealed subtle behavior in the digitizing hardware of the SQUIDâs
electronics. A misconfiguration of the SQUIDâs A/D converter caused some of the 2
16
digital output
levels in the 16 bit A/D converter to be preferred over adjacent levels. This preferential-levels behavior,
which may cause system non-linearity, was analyzed and documented in a GP-B science document. A
recent analysis of on-orbit data highlighting possible non-linear behavior is shown in Figure 6.
Arc-sec
Time (seconds)
78 Gravity Probe B Science ResultsâNASA Final Report, December 2008
Figure 6
. Gyroscope 4 misalignment June-July, 2005. Apparent shifts in gyroscope misalignment correspond to
times when the gyroscope readout calibration signal was toggled on/off.
Shifts in the apparent misalignment occur when the readout calibration signal was toggled on/off and they
occur in the radial direction relative to zero misalignment. These shifts in gyroscope misalignment are
suggestive of offsets in the dc gyroscope scale factor and are likely due to non-linear behavior in the
gyroscope readout system.
6 Torques other than those due to patch effect potentials
Gyroscope torques cause Newtonian gyroscope drift. There are two approaches to limit the impact of this
drift on the experiment error. One approach that has been used with patch potential torques is to model
the underlying physics. Sensitivity tests and gyro-to-gyro comparisons provide a measure of residual
error. A bottom-up approach has been used for other Newtonian torques including support dependent and
support independent effects. Sources for these torques include gyroscope physical contamination,
magnetic contamination, differential gas pressure damping, and gyroscope mass unbalance and gyroscope
geometry imperfections. We have performed a preliminary analysis and find that the Newtonian drift
associated with these torques is less than 1 mas/yr.
7 Summary
There has been significant progress in bounding the total GP-B systematic uncertainty. Many effects
have been bounded to less than 0.1 mas/yr to 1 mas/yr. We have described five approaches to
determining the total uncertainty. The most recent models yield reduced scatter in sensitivity test results
and provide very encouraging gyro-to-gyro agreement in the relativity rate estimates, however more work
Gyro 4 Misalignment
Calibration signal switches
(Arc sec)
(Arc sec)
Muhlfelder: GP-B Systematic Error 79
remains. We need to extend the analysis of support dependent and support independent torques to include
the effects of dissipation on the experiment uncertainty [Buchman et al. (2009)]. Completion of these
analyses will provide the underpinnings for a report detailing the contributions to the uncertainty from
more than 100 sources. We continue to improve the science analysis, refining the modeling of the
physics governing the behavior of the GP-B hardware. A large effort, now in the early stages of
implementation, will track the gyroscope dynamics using a 2 second time step rather than the current once
per orbit time step. It is likely that reduced uncertainty due to systematic effects will be a vital by-
product. Sensitivity tests and gyro-to-gyro comparison of the science results using these refined models
will provide the final elements in determining the total experiment uncertainty.
Acknowledgments
The authors thank S. Buchman, J. Conklin, F. Everitt, M. Heifetz, D.Hipkins, Y.Ohshima, A.Silbergleit,
and J. Turneaure for their useful discussions and insights into this work.
References
Buchman, et al. To be published in Review of Scientific Instruments (2009)
C.W.F.Everitt et al., Report on a Program to Develop a Gyro Test of General Relativity in a Satellite and
Associated Control Technology (1980)
M. Heifetz et al., Space Science Reviews, (2009, Chapter 4 of this report).
G. M. Keiser, "Gravity Probe B" in Proceedings of the International School of Physics 'Enrico Fermi',
Course CLXVII - Atom Optics and Space Physics, to be published.
G.M. Keiser et al., Space Science Reviews, (2009, Chapter 3 of this report).
I. Shapiro et al., APS Conference, Jacksonville, FL (2007).
A. Silbergleit et al., Space Science Reviews, (2009, Chapter 2 of this report).
80 Gravity Probe B Science ResultsâNASA Final Report, December 2008