Gravity Probe B data analysis: I. Coordinate frames and analysis models

Gravity Probe B (GP-B) was a cryogenic, space-based experiment testing the geodetic and frame-dragging predictions of Einstein's theory of general relativity (GR) by means of gyroscopes in Earth orbit. This first of three data analysis papers reviews the GR predictions and details the models that provide the framework for the relativity analysis. In the second paper we describe the flight data and their preprocessing. The third paper covers the algorithms and software tools that fit the preprocessed flight data to the models to give the experimental results published in Everitt et al (2011 Phys. Rev. Lett. 106 221101–4).


Introduction
In 1960 Schiff [2] showed that the spin axis direction of a free gyroscope orbiting the Earth would exhibit two drift rates due to general relativity (GR): (1) geodetic, due to the gyroscope's motion through the curved spacetime about the massive Earth, (2) frame-dragging, due to the dragging of spacetime by the rotating Earth. Gravity Probe B (GP-B) was designed to measure these two effects. An overview in this issue by Everitt et al [3] gives the concept and implementation of the experiment, critical on-orbit gyroscope observations that motivate the data analysis, and a summary of results. The early history of the GP-B development can be found in [4]; a detailed description of on-orbit operations is presented in [5]; the data analysis status at the end of 2008 is summarized in [6]. In this, the first of three data analysis papers, abbreviated DA I, we derive the required theoretical models supporting the experiment. The second part, DA II [7], discusses the experimental data and the preparation of these data for the final analysis. Estimation algorithms and their computational realization, as well as the resulting relativity estimates and associated uncertainties are found in the third part, DA III [8].
To perform the GP-B experiment, a gyroscope, located within a spacecraft, was placed in a circular polar orbit about the Earth. The two GR effects are perpendicular to one another in this orbit allowing them to be independently measured (figure 1). The spacecraft tracked a nearly fixed reference star (guide star, GS) with an on-board telescope, so that the gyroscope drift rate could be referenced to inertial space. A superconducting quantum interference device (SQUID) measured the drift using the gyroscope's London moment (LM), a magnetic marker that is intrinsically aligned with the spin axis. The spacecraft rolled about the line-ofsight to the star, thereby modulating the drift data at roll frequency; measurement of roll phase allowed demodulation. This paper has six further sections, as follows: II. The GR predictions for the orbiting gyroscope are discussed and the reference frames needed for the data analysis are defined, along with the nearly fixed vectors of gyro spin, spacecraft roll, and their difference the spin-to-roll misalignment.
III. We introduce the time-varying influence of the rotorʼs trapped magnetic flux on the LM signal. This influence is analyzed by using the SQUID signal after it has been lightly filtered to a 1 kHz ∼ bandwidth; a signal henceforth referred to as the high frequency (HF) signal. The model for its analysis, called trapped flux mapping (TFM), is derived in appendix A. TFM allows the trapped flux (TF) signal to be converted from volts to arcseconds; the conversion quantity is referred to as the scale factor. TFM, described in detail in DA III, also provides the evolving polhode, critical to analyzing scale factor variations and the torques discussed later in this paper (sections 4 and 5).
IV. The model for the SQUID signal, heavily filtered onboard, is derived in this section. This low frequency (LF) signal is described by the main measurement equation for the experiment. It contains a gyroscope scale factor term comprised of three components: the LM, the ∼1% contribution from TFM and a ∼1% adjustment from SQUID electronics variation. Other most important submodels involved in the measurement equation are gyro spin motion caused by patch effect (PE) torques (section 5), and the model of telescope pointing (section 6). The measurement equation is mainly used for that part of the orbit when the telescope tracks the guide star (guide star valid or GSV); a generalization of this central equation is found in appendix B for guide star invalid (GSI) periods during which the Earth occults the line-of-sight to the star. V. Contrary to the pre-launch expectations [9], disturbance torques had a significant impact on gyroscope drift rate. This section gives the exact solution to the equations of motion for the trajectory of the spin axis in the presence of torques due to patch potentials on the rotor and housing. In appendix C we derive these torques and examine implications.
VI. Here we consider telescope pointing and give the model for one phenomenon that is essential to the experiment, the orbital and annual aberration of starlight. This effect causes the telescope to point towards an apparent rather than true guide star position. The apparent position is precisely known so that aberration may be used to accurately calibrate the telescope and LM readouts. Next, we develop the model for apparent position variation due to parallax and the bending of the guide star light by the Sun. Pointing variation also arises from telescope and spacecraft limitations. We describe the model for this mispointing and give the scale factor matching method that removes its effect upon the relativity result.
VII. This section summarizes all of the GP-B analysis models, which in combination provide the foundation of the estimation algorithms and software that fit the flight data to give the GP-B relativity result.

Coordinate frames. Relativistic drift predictions. Spin-to-roll misalignment
Practically all the models we used to analyze the GP-B data are either expressed in terms of certain coordinate frames, or are derived using them. So it is natural to start with the definitions of those frames and the quantities, mostly vectors, that are needed for our analysis. We also discuss here the GR predictions for the gyro drift, to clarify exactly what GP-B was measuring.
The  for the Cartesian coordinates of the inertial frame JE2000 [10], with the unit vectors x x y x z xˆ,ˆˆ,ˆ1 2 3 = = = along the corresponding axes. We need a Cartesian inertial frame more appropriate for the GP-B data reduction. It is natural to ( 1 ) gs δ α δ α δ = + + here δ is the declination and α is the right ascension of the GS (for IM Pegasi-used as the GS in the GP-B mission-343.25949461 α =°and 16.84126106 δ =° [11]). Now, the z axis of the JE2000 frame was exactly parallel to the Earth rotation axis on noon GMT, 1 January 2000, and stayed within a few arc-seconds (as; 1 as 4.848 10 rad The index WE stands for the West-East direction perpendicular to the ideal orbit plane, in which the gyroscope drifts due to the relativistic Lense-Thirring effect (see the next section). We define the third axis in the usual way (see figure 2): ns we gs

≡ ×
The unit vector ê ns lies in the ideal orbit plane and is orthogonal to ê gs ; as shown below, the geodetic relativistic drift goes in this NS direction. From the definitions (2) and (3) using expression (1) we find ê ns and ê we in the basic frame JE2000: The GS-NS-WE right-handed Cartesian inertial frame with the unit vectors e e e {ˆ,ˆ,ˆ} gs ns we is the most used in the GP-B data analysis.
2.2. GR drift rate predictions and the gyro spin axis in inertial space 2.2.1. Geodetic and frame-dragging effects predicted by GR. In 1960, L I Schiff [2], using the linearized GR theory (weak field approximation, e.g. [12], ch 18) showed that an ideal gyroscope in orbit around the Earth would undergo a relativistic precession according to the equations Here M is the Earthʼs mass, I is Earthʼs moment of inertia about its rotation axis, e ω⃗ is the Earth angular velocity vector (whose direction was close to zˆduring the mission), G is the gravitational constant, and r ⃗ and v⃗ are the S/C orbital radius and velocity, r r = | ⃗ |. Notation g Ω ⃗ stands for the geodetic (or gravito-electric, or de Sitter [13]) precession caused by the gyro motion in the spacetime curved by the massive Earth, and fd Ω ⃗ is the frame-dragging (or gravito-magnetic, or Lense-Thirring [14]) precession caused by the Earthʼs rotation.
The non-precession terms in the first equation (6) average out over a closed free-fall orbit [2,19], and the change in sˆis entirely negligible over such short periods of time; thus the orbit averaged dynamics becomes purely precessional: orb Ω = 〈 ⃗ 〉 × Since g Ω ⃗ is proportional to the orbital angular momentum, both the instant vector and its orbit average are perpendicular to the orbit plane. The average magnitude of the precession is GMv c r v v 3 2 ( ) 2 2 = | ⃗ | for a perfectly circular orbit by Schiffʼs formula (6). If the orbit is also precisely polar, and the GS is exactly in the orbit plane, then the geodetic precession is totally in the WE direction defined by (2). For a circular polar orbit, with zê e ω ω ⃗ = , the x-and ycomponents of the Lense-Thirring precession fd Ω ⃗ average to zero, leaving just the zcomponent equal to GI c r 2 e 2 3 ω , so the WE precession due to frame-dragging vanishes for a perfect orbit. Thus the proper polar orbit provides the best spatial separation of the two relativistic effects: they are orthogonal in this case. The geodetic drift is in the NS direction, and the frame-dragging drift in the WE one. In reality the GP-B orbit differed but slightly from the ideal one: the GS was off the orbit plane by 2 10 4 ∼ × − on average, and the Earthʼs rotation axis deviated from the z axis by less than 10 4 − . So the frame-dragging effect contribution to the NS drift rate ( 0.04 mas yr ≲ ) and geodetic effect contribution to the WE drift ( 1 mas yr ≲ ) were both negligible.
Formulas (6) are rigorously valid for the spherically symmetric gravitational field of a rotating mass M, which is not exactly the case with the Earth. The effect of the Earth quadrupole moment, J 1.083 10 2 3 ≈ × − , on the geodetic precession was found independently by D C Wilkins and R F O'Connel in 1968 (see papers [15] and [16]). The correction to the orbit averaged geodetic precession due to the Earth oblateness and small orbit eccentricity was calculated in [15,17] using rather elaborate analytical techniques; its simple derivation was presented later by J V Breakwell [18].
The general treatment of both the geodetic and frame-dragging precessions in the gravitational field that contains all the multipoles was given in [19]. In particular, it was shown that the desired accuracy 1 mas yr ∼ for measuring the geodetic effect requires including the Earth quadrupole moment only; all other terms are two to four orders of magnitude smaller [20]. The proper expressions from [19] based on Earthʼs gravitational potential (r e is the Earth equatorial radius), were used to calculate the instantaneous values of g Ω ⃗ and fd Ω ⃗ (in fact, the J 2 correction to the frame-dragging was negligible). These were computed for every orbit from the precise GPS orbit data [21] using each available valid data point. The averages found differed slightly from orbit to orbit due its small classical variation, so we Ω 〈 〉and ns Ω 〈 〉 were additionally averaged over all orbits, and relativistic drift rate predictions were determined from the result given by the formulas (11). These resulting numbers for r ( 6, 606 mas yr) ns − and r ( 39 mas yr) we − are shown in figure 1 (precise definitions of these quanitites are found in section 2.2.3); note that the contribution of the geodetic precession in the WE drift rate, due to the orbit imperfection, was not completely insignificant, 1 mas yr ∼ . As an important cross-check, the averages of the non-precession terms in the first of the equations (6) were also computed and found to be small as expected: 0.1 mas yr at most, and typically much less.
2.2.2. Gyro spin axis in inertial space. Starting here and throughout all the three papers on data analysis, we use the term '(gyro) spin axis' in place of 'instant (gyro) rotation axis', for brevity and out of the GP-B tradition. As explained below, for GP-B this axis pratically coincides with the axis of the angular momentum (usually called 'spin axis' in physics, and meant in the Schiffʼs equations (6)), so we will not make a difference between the two.
After the on-orbit spin-up of all the GP-B gyroscopes [1,3,5], their spin axes were torqued to point within less than 10 as of the direction to the GS. During the mission the axes remained within less than 20 as 10 4 ≈ − of the GS, with a few arc-seconds average over the whole mission. Consequently, we have (writing sˆfor the unit vector of a gyro spin axis looking towards the GS): The vector s t ( ) ⃗ to lowest order (l.o.) is two-dimensional and lies in the plane perpendicular to the GS direction. It describes the drift of the gyroscope, i.e., its change carries the information on the relativistic drift rates; thus we can call s ⃗ the drift vector.
More precisely, GP-B had to trace the change in the vector of angular momentum, rather than in the spin vector; however, the difference proves to be insignificant, as mentioned in the beginning of this section. Indeed, in inertial space the instantaneous rotation axis of a free rotor precesses about its angular momentum vector, and the sine of the angle between them is proportional to the relative difference of the corresponding moments of inertia (e.g. [22]). For the GP-B gyroscopes I I 10 6 Δ ∼ − [3,23,24] for any two body axes, so the angle between the spin and angular momentum is of the same order: very tiny. In addition, sˆaveraged over each spin period is exactly aligned with the angular momentum. Since the sampling rates of all the GP-B LF science signals ( Hz 0.5 ) are much smaller than the spin speed (∼60 or 80 Hz), the data are effectively spin-averaged, so we can use s t ( ) as the unit vector of angular momentum as well. Gyroscopes #1 and #3 were rotating counterclockwise, and #2 and #4 clockwise around the s t ( ) direction.
ns ns 0 ns 0 we we we 0 being the initial spin axis position at t t 0 = . However, if non-negligible classical torques are present, their contribution must be added to the rhs of equations (12), and it appears then in the gyro trajectory s t ( ) ⃗ as well. As stated in the introduction, this is exactly what happened with the GP-B gyroscopes because of the PE torques.   I  I  I  I  0 1, We call p γ and p ϕ the polhode angle and phase; given s ω , the free motion of a gyroscope is described by these two quantities (see appendix A). These two angles play a very significant role in the whole data analysis; if Q 0 = , both change with the time and depend on the value of the asymmetry parameter.   τ τ are small. Indeed, during GSV the roll axis nearly coincided with the S/C axis, i.e., the boresight axis of the telescope, the pointing sensor. The S/C axis, and hence also the roll axis, deviates from ê gs mainly by the amount of the optical aberration of stralight (see section 6.1), within 20 as 10 4 ∼ ≈ − . During GSI the roll axis deviated from the S/C axis by at most 1 as ∼ . The S/C attitude controlled by rate gyroscopes during GSI was still within 2 as ∼ of that during GSV, and significnatly less on the average. Therefore the deviation ofτ from ê gs was 25 as 1. 25 We call x y {ˆ,ˆ,ˆ} τ τ τ the roll axis frame; since it wanders in inertial space no more than 10 4 − , it is called quasi-inertial. This is clearly seen also from the zeroth order expressions of its unit vectors, e x e y êˆ,ˆˆ,ˆˆ, rendering them to the GS inertial frame of section 2.1. The roll axis frame is used to derive the main measurement equation for the SQUID signal (section 4.5 and appendix B). When deriving it, we are interested in exactly those small deviations, to establish the linearity of the GP-B readout. For this reason, we carefully kept quadratic terms in all the formulas (18) and (19).
Like the vector s ⃗ , the vector τ ⃗ is two-dimensional to l.o.; it is called the pointing vector, and gives the deviation of the roll axis from the true direction to the GS. We now introduce the vector of the spin-to-roll axis misalignment, μ⃗ , as the difference betweenτ and sˆ(figure 4). Using formulas (9) and (20) where μ is taken to l.o.,  To gain more insight in the meaning of μ⃗ , let t ( ) ψ ψ = be the angle between the spin and roll, andψ the unit vector of the corresponding rotation axis, s ŝˆˆ. ψ τ τ = × | × | Using sˆandτ first by equations (8), (17), and then by equations (9) and (20), one calculates: The above expressions imply This equality is rather obvious: the arc (ψ) between two close vectors of the same length ( sˆˆ1 τ | | = | | = ) is equal, to l.o., to the length (μ) of the chord connecting the ends of these vectors (see figure 4). This clarifies the meaning of μ: to l.o. it is the angle between the spin and roll axes. Results (24) and (26) The last expression using representation (23) remains finite and meaningful even when 0 μ → ; the corresponding limit value of θ is also well defined, since μ⃗ in the limit position is perpendicular to the aligned vectors sˆandτ.

Frame fixed in the S/C body and the pick-up loop orientation
Here we define the Cartesian frame x y z {ˆ,ˆ,ˆ} b b b attached to the spacecraft body. We choose the unit vector zˆb along the S/C axis, i.e., the boresight axis of the science telescope. We can take any other two orthogonal body-fixed unit vectors x b and yˆb such that the three form a right-handed Cartesian frame. It is convenient to take the vector x b in the plane formed by zˆb and the normal, n, to the SQUID pick-up loop. Since the frame we are building is only an intermediate, it might be different for each GP-B gyroscope; the GP-B instrument set-up allows, however, for one common frame for gyroscopes #1 and #2, and another one for gyroscopes #3 and #4.
Due to unavoidable assembly imperfections, the pick-up loop plane is at a constant angle to the axis zˆb, const 55 as 2.75 10 .
According to the frame choice, the body-fixed pick-up loop normal is thus We want now to convert this to the inertial frame for any GSV period. During GSV, the S/C roll axis coincides with the telescope axis, zˆˆ.
The body-fixed vectors x b and yˆb rotate about zˆτ = of the roll-axis frame in the previous section. Hence they are expressed through the vectors x τ , equation (19), and yˆτ, equation (19), by the rotation R ( ) r r 3 ϕ δ + with the satellite roll phase, r ϕ : These formulas are obtained by substituting expressions x τ and yˆτ from (19), with quadratic terms in them. We wrote const r δ = for the roll phase offset, i.e., the angle between x τ and x b at a moment when 0 (mod 2 ) We can now introduce expressions (31) and (32) into equation (29) to obtain n in the inertially fixed frame during GSV. After carefully gathering all similar terms, n acquires a relatively compact form: data allowing us to determine the gyro spin and polhode frequencies and the spin-down rate. The role of this HF signal for the main science analysis was significantly re-evaluated onorbit, after finding the changing polhode period, and especially after discovering the anomalously large PE torques whose imprint on the gyro motion strongly depends on polhoding (see section 5). We developed and implemented a special procedure to analyze the HF signal called TFM; it is described in detail in our paper [24] and briefly below. In appendix A we give the complete derivation of all the model equations for TFM for the first time. The ultimate goal of TFM was to estimate accurate polhode parameters and LF SQUID scale factor variations throughout the mission, which were critically important for the science analysis.

Origin and structure of HF SQUID signal
The dipole LM field is the main but not the only source of the GP-B magnetic readout. When a thin film superconductor is cooled below the critical temperature, residual magnetic field, which existed even under the unique low-field GP-B conditions [3], is trapped inside the rotor. Magnetic quanta 0 Φ ± comprise point sources on the rotor surface (fluxons). The multipole field of these point sources contributes to the total flux through the SQUID pick-up loop, which rolls with the S/C in inertial space; the roll period was T s 77.5 r = . The time signature of the TF data emerges from the following argument. Fluxons are frozen in the rotor surface and spin with it; the function that converts a fluxonʼs position into its magnetic flux through the pick-up loop is strongly nonlinear (almost a step-function, due to the relatively small gap between the rotor and the loop [25]). Therefore the TF signal consists of multiple harmonics of spin, and, since the pick-up loop is rolling, they are actually the harmonics of ± spin + roll frequency (recall that two gyroscopes were rotating counterclockwise, and the other two clockwise). Finally, since the spin axis moves in the rotor body (polhoding), the spin harmonics are modulated by the polhode frequency. All this combined leads to the following expression for TF through the pick-up loop (this formula is derived in detail in appendix A resulting in (A.18)-(A.20)): (by section 2.2, 0 s ϕ > for gyroscopes #1 and #3, and 0 s ϕ < for #2 and #4). Here t ( ) β is the angle between the gyro spin axis and the pick-up loop plane shown to be 10 4 ≲ − in section 4.3 and appendix B. The even harmonics of spin are multiplied by this small angle for a purely geometric reason: the spin axis resided almost in the pick-up loop plane in the GP-B experiment. Note the n = 0 contribution of TF to the LF signal: we will return to it in section 4.

Sources of the GP-B HF data
During the GP-B flight there were two telemetry sources for the HF SQUID signal: (1) the first six harmonics of spin obtained by the FFT, and (2) the so called snapshots, which were 2 s ∼ stretches of original SQUID signal sampled at 2200 Hz. Both were received only during the part of the orbit when the GS was occulted by the Earth, due to the limited capacity of the on-board CPU. When available, a snapshot would come about every 40 s, but gaps between the snapshot arrays might be up to 2 d.
HF FFT harmonics were analyzed during the mission as soon as the data were available. Snapshots, being the most accurate source of HF information, were thoroughly treated after the mission. All 976 478 snapshots from the mission science period were processed after the flight.

Changing polhode period and path: kinetic energy dissipation
As mentioned, one of the two on-orbit discoveries most strongly affecting GP-B data analysis was the changing polhode period. The measured time history of the polhode period for each of the GP-B gyroscopes is shown in figure 5, with very good agreement of the results obtained using the two HF signals, FFT harmonics (red) and HF snapshots (blue). Notably, the T t ( ) p plot for gyros 1, 2 has a peak where formally T p = ∞. Nevertheless, throughout the science mission (started 13 September 2004) polhode periods of all the four gyros decreased monotonically tending to specific asymptotic values T pa .
Such behavior contradicts the Euler solution for free gyro motion, which gives constant polhode frequency and period. It can only be explained by energy dissipation that reduces the kinetic energy, E, of the rotor without changing its angular momentum, L. Dissipation moves the spin axis in the body to the only stable position, the maximum inertia axis, at which energy is minimum when L const = . So dissipation also changes the polhode path controlled by the parameter L E 2 2 ; in contrast, the spin-down torque changes both quantities L and E keeping this ratio constant, and the polhode path unchanged.
The two types of behavior seen in the plots of figure 5 are also explained: by pure chance, gyroscopes 1 and 2 start their evolution with the spin vector precessing about Iˆ1 (see section 2.3 for notations), so T p tends to infinity when s ω⃗ crosses the separatrix [22,26]. Contrary to this, the evolution of gyrosopes 3 and 4 starts with s ω⃗ already precessing about Iˆ3, so T t ( ) p just decreases monotonically. From angular momentum conservation one can readily determine the relative energy loss and spin speed reduction when the spin axis moves in the rotor body all the way from the minimum, Iˆ1, to the maximum, Iˆ3, inertia axis. They both are equal to the same ratio The gyro kinetic energy is J 1 ∼ , so the total energy loss when s ω moves from Iˆ1 to Iˆ3 is less than J 4 μ ; in one year, the average dissipated power required for this is just 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 PE torques [27,28]. Independent of the origin, we have found a general dissipation model in the form of an additional term in the Euler equations unique up to a scalar factor (its brief derivation can be found in [29]). Fitting the polhode period time history obtained from the model to the measurements allowed us to determine the rotor asymmetry parameter, the asymptotic period T 1 2 pa ∼ − hours, and the characteristic time of dissipation 1 2 dis τ ∼ − months, depending on the gyroscope (see table 2 in [24] for details). The gyroscopes started the science period with T p ranging between 1.5 and 8 hours; thus the dissipation is slow (T p dis τ ≪ ), so that the polhode motion of the GP-B gyros is quasi-adiabatic.

TFM and its products
TFM is a procedure for finding the distribution of trapped magnetic field and characteristics of gyro motion from odd harmonics of the HF SQUID signal by fitting them to their theoretical model. It is based on the solution for the scalar magnetostatic potential outside the spherical rotor with fluxons on its surface in the rotor-fixed frame. The complete derivation of formula (34)  which is presented in our paper [24]. Here we briefly review some points of TFM relevant to science analysis.
To perform TFM, one first of all needs to process the measured HF SQUID signal, , observed in almost one million snapshots, to extract from them the amplitudes of multiple harmonics of spin, H t ( ) n , according to the formula (34): Because of the small factor β involved in the even harmonics, their amplitudes are orders of magnitude smaller than those of the odd ones. So one uses expressions (A. 19) for H t n ( ), n = odd, through the coeffcients A lm , polhode phase and angle to determine them by fiting to the measured H t ( ) n . 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 data processing provides estimates of the polhode phase and angle, t ( ) p ϕ and t ( ) p γ , the spin phase, t ( ) s ϕ , and the asymmetry parameter, Q; the estimation at levels A and B is essentially nonlinear. Once the three Euler angles are known at all times along with the rotor asymmetry, coefficients A lm can be found in level C by a linear fit. The details of processing are given in [24] and [8].
For the entire duration of the mission, TFM produced spin speed to 10 nHz and the spindown rate to 1 pHz s 1 − [30]; spin phases to 0.05 rad; polhode phases to 0.02 rad, and the polhode angles to 0.01 0.1 rad − [24]. With these results we computed the LF SQUID scale factor variations C t

Charactersic frequences and their values
It is already clear, and will be more evident from the rest of this paper and papers [7,8], that our data analysis was pretty much about proper handling of fundamental characteristic frequencies involved in the GP-B signals. They were numerous, of both natural and artificial origin, and ranged ten orders of magnitude, from 10 8 ∼ − to 10 Hz 1 ∼ . There were two groups of frequencies, one not related to gyroscopes, the other gyro-dependent. For convenience, we list all of them here along with their values.
The frequencies of the first group are: (1)  , to move f dy away from the second harmonics of the roll frequency.
Every gyroscope had its own frequency of rotation (spin), f s , and the polhode frequency, f p , growing through the Mission; their values are shown in table 1.

LF SQUID signal: model for science readout (main measurement equation)
Within the linear range of the SQUID, the LF gyroscope readout signal is proportional to the total magnetic flux, Φ, through the pick-up loop: Here b is the bias, and C el is the scale factor of the whole chain of the electronic onboard signal processing, starting with the SQUID to pick-up loop scale factor. The output signal is a changing voltage, so the physical dimension of both z t ( ) and b is volts, and the dimension of Unfortunately, both C el and b may be not perfectly constant; this issue is addressed later. Meanwhile, to be able to estimate relativistic drift rates, we need first to express the total flux through other variables, in particular, the gyro orientation s t ( ).
between the gyro spin axis and the pick-up loop plane changes due to the relativistic drift of the gyroscope spin axis (as well as the classical torques, if any), and thus contains the quantities we want to measure. Otherwise N * 4 ≈ is the effective number of turns of the pick-up loop, M L is the LM expressed through physical parameters, C Lm is the LM scale factor, r and R are the radii of the rotor and the pick-up loop, respectively, m e and e are the mass and charge of the electron, and f s is the rotor spin frequency. As we explicitly show in section 4.5, 3 β β β β = + = with more than sufficient accuracy of 10 12 ∼ − during GSV. For this reason, we permanently drop the sine in the formula (38) and write the linear LM flux through the pick-up loop as Lm Lm

Φ β =
The LM scale factor C Lm , whose dimension is magnetic flux per angle, with the standard unit as 0 Φ , depends only on the spin speed, which is extremely stable: the spin-down rate of the four gyros is shown to be within 0.3 1.3 Hz hr 1 μ − − [30,32]. Hence C Lm can be considered constant for all the short-and mid-term analyses. However, for 10 months of observations the relative change due to the spin-down can be up to 5 10 5 × − , so we write: According to equation (35), the LF part of the TF contributes to the LF readout, adding to the LM flux. These two contributions provide the total magnetic flux through the SQUID pick-up loop: so the TF adds, in fact, to the SQUID scale factor. At the start of the science mission, the added factor ranges between 0.1 and 5% of C Lm , depending on the gyroscope, and eventually approaches a constant.
, of a given gyroscope with its specific asymmetry parameter Q 0 > . As explained in appendix A, this representation is not convenient for the analysis, so it is converted into the final formula (A.29) using the polhode phase and angle ϵ tend to zero as the spin axis asymptotically approaches its stable position along the maximum axis of inertia. This is seen from the explicit model for these coefficients (formula (A.31), appendix A): Introducing expressions (41), (43) and (44) to the formula (42), we arrive at the final model of the total SQUID scale factor: Of course, in numerical calculations all the infinite series here and elsewhere are replaced with finite sums with the number of terms providing sufficient accuracy. A universal and unambiguously defined procedure for choosing these cutoffs is discussed in the DA III paper.

Long-term variations of electronics scale factor
Next we study the electronics scale factor, C el , from (37), which was supposed to be constant. Accurate investigation of the LF SQUID data has shown that at the level 10 10 3 5 − − − some long term variations were in it, with characterstic times of several weeks or months. Most probably they reflect the changing environment (thermal, electromagnetic, material aging, etc) of the onboard circuitry. There is no way to detect, trace, and model all of them, for which reason the model for C t ( ) el is the only one in the whole data analysis which is not derived from the underlying physics, but is motivated by computational convenience and efficiency; the same model is used also for the misalignment torque coefficient (see section 5.2).
Namely, C t ( ) el , as any function on a given interval of the length T sm ('sm' stands for the science misson, about one year), can be represented by its Fourier series with the basic frequency T 1 sm . Thus the analysis model is , are the constant dimensionless parameters to estimate, and a subharmonics (m 1 2 = ) was added to facilitate secular changes.

Complete scale factor model
Everything is ready now for the complete model of LF SQUID signal scale factor. Introducing (45) and (46) to the formula (37), we find: Since the corrections are small, their products are neglected, and the final model is:  (41), (45), and (46), respectively. When analyzing the real data numerically, the series in the formulas (45), and (46), as well as the series in all other models, are truncated at some finite number of terms, as is discussed in the DA III paper.

Main measurement equation: SQUID signal in terms of spin-to-roll misalignment
Representation (47) is not helpful for the analysis and thus for the determination of relativistic drift rates, until angle t ( ) β is expressed in terms of the experimentʼs observables. By its definition and the formula (39) it is determined to be s n sinˆ·ˆ, where n is the unit normal to the pick-up loop. Using its GS frame expression (33) for GSV, and expression (9) for the unit vector of gyro spin sˆin the same inertial frame, one calculates the dot product in (48)  The last estimate of the remainder here is due to the equality O ( ) β μ α = + implied by the leading term. It is amazing that all the quadratic terms on the right of equation (49) cancel out completely, leaving just the leading linear term and its cubic corrections. For any GSV outside the S/C anomalies the cubic corrections should be not larger than about one part in 10 11 (worst case); this is definitely negligible as compared to 10 4 ≲ − of the main linear term. In the majority of anomalies where pointing was not lost, the spin-to-roll misalignment was still within 500 as ∼ . Even then the corrections remained 10 10 ∼ − , which is negligible.

Class. Quantum Grav. 32 (2015) 224018 A S Silbergleit et al
After introducing expression (49) in the readout equation (47), we obtain + may be treated as the new time-varying bias. This is the main measurement equation of the whole GP-B data analysis. Here, as before, z t ( ) is the LF SQUID output (in volts), and C t ( ) (V as ) g 1 − is the scale factor whose model was given in the previous section; b is the constant bias, r δ is the constant roll phase offset, and α is the constant assembly angle. These are estimated in the signal analysis, the latter as a part of the total bias. Next, t ( ) r ϕ is the phase of the S/C roll measured independently by S/C star trackers, and t ( ) τ ⃗ , entering through its NS and WE projections, is also known for any GSV (see section 6.1). Relativistic drift rates r ns and r we are hidden in the gyro orientation described by the projections of the vector s t ( ) ⃗ from equation (9)  . If all the classical torques were negligible, these would be linear functions of time (13) with the slope coefficients r ns and r we ; the real dependences, much more complicated due to PE torques, will be derived in the next section. Note that the SQUID signal model for GSI is derived in appendix B; it was not used for estimating relativistic drift rates, but only for determining the S/C pointing, t ( ) τ ⃗ , during every GSI period. It has just one more term than equation (50). The bias variability can be modelled by some additional states, B j , However, the bias is totally removed before the final relativistic drift estimation. The analyzed signal was free of both the bias and pointing error removed using the telescope signal (see section 6.3).

Linear readout and nonlinear estimation problem
We have thus proved that the GP-B readout is linear in all the small angles , , μ α β, etc, with only cubic corrections; those are several orders smaller than required by the desired precision. Nevertheless, the problem of estimating the relativistic drift rates and other parameters by the equation (50)  ( ) ⃗ , are mutilplied by a scale factor that also needs to be estimated-yet another nonlinearity! However, under some natural simplifying conditions the measurement equation (50) can be reduced to the linear form (52) by just a change of variables [34][35][36]. This is important for fully understanding the difficulties of the GP-B data analysis, so we discuss it in more detail.
Let us assume for a moment that the scale factor variations, bias and classical torques are negligible, so that    can be determined uniquely through them. The covariance matrix of the latter can also be recalculated from that gotten for x n , giving the full solution to the problem.
Next it occurs that dropping the first of the simplifying assumptions (53), i.e., allowing for the scale factor variation, does not spoil the picture much, because those variations are relatively small and can be treated as perturbations. Indeed, the full model (47) for the scale factor can be written as . The next iteration starts with the LSQ fit for x n using the found values c n (1) in equation (57), and thus determines x n (2) ; inserting them in equation (57), one concludes the second iteration by finding c n (2) , and so on. In this way, only a linear estimation problem needs to be solved at each of the two steps of any iteration.
Our results show that it is only the combination of relatively small but significant variations of the gyro scale factor with a nonlinear time change of the gyro orientation s t ( ) ⃗ that makes the main GP-B estimation problem essentially and strongly nonlinear. It can be no longer reduced to iterations of linear estimation problems. Contributions of the PE torques to the gyro motion were not only non-negligible, but could not be even treated as small perturbations of the relativistic drift.

PE torques on GP-B rotors
The electrostatic PE [37][38][39] is a non-uniform distribution of electrical potential on metal surface. It occurs when a non-uniform dipole layer forms on the surface due to impurities or microcrystal structure, resulting in a varying potential and electric field not necessarily perpendicular to the surface. In the idealized case of an isolated conductor, the net force and torque on the body are still equal to zero. But for two metallic surfaces at a finite distance, the net force and torque on each of them do not in general vanish. Naturally, the effect is larger, the closer the surfaces, as first confirmed by the calculation of the PE force for two parallel conducting planes [40].
In GP-B, the two metallic surfaces close to each other are the rotor and its housing. Various manifestations of PE in our experiment were studied in paper [28], mostly by fitting a PE model to experimental data of different physical origin. PE torques acting on GP-B rotors are discussed in [41]: their experimental discovery is described, and some key theoretical results are presented (see also [3,5]). A theory of PE torques is developed in detail in appendix C. Its results effectively simplify due to two small parameters: (1) the ratio, d a ( ) 10 3 ∼ − , of the rotor-to-housing gap, d, to the rotor radius, a, and (2) the spin-to-roll misalignment, 10 4 μ ≲ − (see section 2.4, in particular equations (21) and (22)). Based on these results, we examine the motion of GP-B gyroscopes with PE torques acting on them and Two distinct types of PE torques were found during and after the flight mission, which were confirmed by the PE torque theoretical model derived in appendix C. They are the misalignment torque and the roll-polhode resonance torque.
The misalignment torque was discovered in the post-flight calibration period of the experiment when the S/C was deliberately manouvered up to 7 degrees away from the GS. The torque increased with the misalignment, and, within a 1 ≲°range, was proportional to it; the direction of the torque was perpendicular to the misalignment vector μ⃗ . This is exactly as one of the PE torques behaves, as shown in appendix C.
The roll-polhode resonance torque was discovered in the data as sizeable 'jumps' in the estimated gyro orientation: from time to time one of the gyroscopes would change its spin direction by as much as 20-100 mas in one day or less. The jumps occured when some high harmonic of changing polhode frequency became equal to the roll frequency. Before this, one additional PE torque was found theoretically; it had zero order in the misalignment (i.e., was independent of it), thus potentially it could be very large, without the misalignment factor 10 4 ≲ − . But it was proportional to the sine or cosine of roll and thus usually averaged out during each roll period; so it was neglected. However, since those roll harmonics were modulated by the polhode ones, the roll averaging broke down at the times when the resonance condition, n r p ω ω = , holds. This led to a significant gyro drift in the vicinity of the resonance, when the frequency detuning was small. The contributions of the roll-polhode resonance torque are given by the third term on the right of equations (58), with t n t t ( ) ( ) ( ) n r 0 Δϕ ϕ ϕ = − , r ϕ and 0 ϕ being the known roll and polhode phases (of the corresponding symmetric gyroscope, see below). The resonance occurs when the roll frequency coincides with some harmonics of the polhode frequency, n r p ω ω = . Note that the value of n necessary for a resonance is large, from 40 up to about 300, depending on the gyro. Therefore it was difficult to believe that polhode harmonics of such high order could cause a significant effect until it was seen in the data.
Unlike the general formulas (C.26), only the terms coming into resonance during the flight were included in equation (58)  for gyroscope #2. Its value is determined by the span between the start and end values of the polhode frequency, and by how they compare to the roll frequency.
The time distribution of resonances depends also on the rate of the polhode change, the largest in the beginning and monotonically decreasing to the asymptotic zero. Due to this, the resonances are dense at the start and become more and more scarce towards the end of the mission. The time gaps between the neighbor resonances range from a few days for gyroscope #2 to several weeks and even a month for gyroscope #3. The picture of resonances for gyroscope #2 may be seen in figures 13, 14 of paper [7]; only the resonances whose contribution to the drift was visible were included in these plots. Naturally, the number, n, of the occurring resonances decreases with the time, since the polhode frequency increases: higher order resonances come earlier.
As shown in appendix C, the misalignment torque coefficient k k t ( ) = is modulated by the polhode motion exactly like the TF scale factor. The only difference is that the modulation is caused not by the motion of magnetic fluxons on the rotor relative to the pick-up loop, but by the relative motion of voltage patches on the rotor and housing. Therefore the model for k t ( ) (formula (C.24)) is exactly the same as the model (44)  Just as with the electronic scale factor C el (section 4.3), small long-term variations in k were found by signal analysis. They could come from very high spherical harmonics of the patch distribution whose significant presence in the field is known. This translates, in particular, in a very large number of terms needed in the expansion (59) for k ( ) c 0 0 ϵ to get it accurately, which is difficult to implement.
To capture these long-term variations, a Fourier series in harmonics of the period T sm , the duration of the science mission, was added to the expression (59) for k t ( ). The addition has exactly the same structure as the variable part el Δ of the electronic scale factor in the formula (46). It is important to note that relativistic drift rates and the misalignment torque coefficient were found simultaneously from the science signal. For this reason, any model of the torque coefficient that adequately describes its variation is sufficient to determine the relativistic drift.
The resonance torque contribution to equations (58) already includes modulation by polhode harmonics, so the coefficients a t b t ( ), ( ) This expresses the uniqueness of a trajectory of the dynamical system (61) speicified by any initial point on it, an implication of the arbitrary time shift that does not change the system. In our case the system is linear, which is why it is straightforward to derive the property (67) from the definition (66) by simply multiplying appropriate matrices. The two properties of  imply Using expression (68) in the solution (65), we write the latter in a different way: A glance at our main measurement equation (50) shows that we have already provided models for all the quantities involved there, except for two components , ns we τ τ of the pointing vector τ ⃗ . The equation is used for relativity estimates only during GSV periods, when the the science telescope serves as the S/C attitude reference. Accordingly, we need the telescope pointing during GSV, given by The first three terms on the right are optical effects causing the telescope to deviate from the true direction to the GS, namely, the optical aberration (consisting of the annual and orbital motion contributions), the parallax, and the bending of the guide star light by the Sun. The first two are classical (although we use the first special relativity correction in the aberrations, see below), while the third one, the bending of light, is the famous GR effect absent in Newtonian gravity. The last term in equation (70), t ( ) in θ ⃗ , is the telescope inertial pointing Class. Quantum Grav. 32 (2015) 224018 A S Silbergleit et al error whose value is provided by the telescope signal as described in the next two sections. The magnitude of the annual and orbital aberration is 20 as ∼ and 5 as ∼ , respectively. The parallax, due to the Earth motion around the Sun, is about three orders of magnitude smaller, 10 mas ∼ , and the bending of light is comparable, 20 mas ∼ . The pointing error during GSV was oscillatory with 100 mas ∼ amplitude, for the majority of orbits. Let R E ⃗ and V ⃗ be the position and velocity of the Earth on its orbit, so that R E ⃗ is the vector from the Sun to the Earth. For any inertial velocity vector w⃗ let us introduce a vector operation gs ns ns we we Then the annual aberration including the first relativity correction [44] is: The orbital aberration differs only by the replacement of V ⃗ with v⃗ ,

( )
v⃗ being the S/C orbital velocity. However, the total aberration A ⃗ is not simply the sum of the annual and orbital ones: The inequality here is due to the relativistic term making the operation w ( )  ⃗ ⃗ nonlinear; to l. o., i.e., in the non-relativistic approximation, the total aberration is, of course, just the sum of its two constituents. By the way, the special relativistic correction to the velocity addition is too small to matter, vV c 10 2 3 13 ∼ − at most. The time signature of the annual aberration components is a sinusoid with annual frequency and the proper amplitude and phase. The same is true for the parallax, = × is the distance to the GS [11]. However, not only the amplitude, but also the phase of its sinusoid is different from that of the annual aberration: the parallax depends on the Earth position, rather than its velocity.
The time signature of the orbital aberration projected is also a sinusoid with orbital frequency. For the ideal GP-B orbit v 0 we = , and the entire orbital aberration is in the NS direction. In reality, the WE orbital aberration was small, with the amplitude 10 3 < − of the total 5 as.
It remains now to produce the model for the bending of star light. Although the derivation of this effect is found in a large number of books, the form required by our data analysis is hard to find anywhere. So we took the result closest to the one needed from paper [45], and developed it to the proper answer, namely: is the geometrical mass (gravitational radius) of the Sun. This result is valid in the weak field approximation of GR to l.o. in m R E , also taking into account that r R L E ≪ ≪ , with L being the distance to the GS, and r being the S/C orbital radius. The time signature of the star light deflection signal is a complicated function containing all harmonics of the annual frequency.
All three effects (72)-(74) were precalculated for the entire science period at the same moments with the s 2 intervals as all other signals (SQUID signal, telescope signal, etc). The Earthʼs position and velocity were taken from JPL Ephemerides, and the S/C velocity was provided by the most precise post-processed GPS data. Finally, the pointing error came from the telescope signal, and was eventually removed from the SQUID signal (see the next two sections). In principle, any other optical effect can be added to those included, if known; one such effect is the orbital motion of the guide star IM Pegassi, but its amplitude is way too small, 1 mas ∼ . This would be the end of the story, but for the fact that the misalignment torque was acting continulously, so t ( ) τ ⃗ was permanently involved in the drift rate, see formulas (64) and (69). For this reason, we needed S/C pointing during GSI as well, even though these periods were not used for relativity estimation: to calculate the misalignment torque contribution at a given moment of time, we had to find it for all previous moments of time, starting with the inital time of the data stretch. During GSI the S/C attitude reference was not the telescope (see paper [42]), so expression (70) was not valid there. In fact, t ( ) τ ⃗ during GSI was found using gyro SQUID signals in the following way.
The gyro orientation does not change much during one GSV: even with the PE torques acting, the change is within 10 mas 5 10 8 = × − at most. The SQUID scale factor also does not change much: its relative variation is just 10 3 − or less. So the average GSV values C g av and s av ⃗ are meaningful and accurate enough. Thus the SQUID signal of every science gyroscope was analyzed at each GSV period (more than 5000 throughout the mission) to determine C g av and s av ⃗ using measurement equation (50) without either torque or scale factor variation modeling (the estimation problem is then converted to a linear one using the two-step scheme, see section 4.6). After this, for any GSI one specifies C g * , s * ⃗ as the average of the corresponding values from the two neighboring GSV stretches. Then, by fitting the model (50) to the SQUID in a window of the length of the roll period or larger, one finds the estimates of the weighted components of the spin-to-roll misalignment, m C g ns ns μ = and m C g we we μ = . The approximate pointing direction at the middle of the window is calculated by the definition (21) of the misalignment: In practice, we used a moving window of the length twice roll period with a 2 s time step through each GSI, getting the values of pointing at the proper time instants. This procedure was carried out for each of the four gyroscopes, and the final numbers for pointing were found by averaging the four appropriate values. In this way t ( ) τ ⃗ was obtained for the whole mission, both at GSV and GSI; the GSI pointing was refined later by a more sophisticated modeling taking into account some features of the GSI spacecraft motion, as described in [42]. The telescope signal model is derived from the telescope system design and from on-orbit analysis of its pointing signal as a function of the pointing measured by the four gyroscopes, as described above. The scheme and construction of star-tracking telescope are given in detail in [3,5]. The telescope observes the guide star when it is not occulted by the Earth (GSV) serving as the S/C attitude sensor. The telescope pointing signals are produced as follows.
A beam splitter in front of the focal plane reflects half of the light and transmits the rest. This arrangement provides two focal planes. Roof prisms with very sharp and straight edges are placed at the focal positions: one roof prism is oriented to measure the telescope position relative to the S/C x-axis and the other about the y-axis. When the telescope is exactly pointed toward the GS, half of the light is reflected from each side of the roof prism, but when the telescope is pointed slightly away from the GS more light is reflected from one side than the other. The light reflected from each side of the roof prism is directed toward one photodiode of a photodiode pair. The measured photocurrents from the photodetector pair, i + and i − , whose values depend thus on the deviation from the GS, are used to derive the telescope pointing signal. For redundancy there is a beam splitter in front of the primary photodetector pair (side A) so that the second photodetector pair (side B) is activated. This photodiode set-up thus generates the following four pairs of signals: , , , , , , . If the ± channels of each photodiode pair were identical, then the normalized telescope signal could be taken as the ratio of the currents' difference to their sum. In reality, these channels are not perfectly symmetric, so a weighting factor compensating for the difference in the light transport from each side of the roof prism, as well as the difference in the photodiode efficiencies, is introduced. Using it, the normalized telescope signal, z T , for a particular axis (X or Y ) and side (A or B) is defined as (here and below we ignore the x y , and A B , subscripts). We estimated the weighting factor w for each photodetector pair by assuming that the weighted sum of the photocurrents is a constant and by using the natural variation in the pointing resulting from disturbances and the response of the S/C attitude control system, We only used data for which the normalized pointing signal was small, z t ( ) 0.2 T | | < , to avoid larger pointing angles where light can be lost in optical transport.
To l.o. in the small deviation from the direction to the GS, the normalized telescope signal (76) is just proportional to the pointing error θ, where C T is the telescope scale factor. However, even for a perfect telescope the normalized telescope signal is inherently a nonlinear measure of the telescope pointing angle, so the question is just how large the nonlinear corrections are. On-orbit measurements of the GP-B telescope pointing demonstrated that the linear approximation is not accurate enough. The correction was represented by a defocus and by a wavefront error using the axially symmetric 4th order Zernike polynomial. In this way we established the conversion from the normalized telescope signal to the telescope pointing angle θ for both axes in the form . With these numbers the nonlinearity correction was accurate to 4% for pointing angles less than 400 mas. The value of the scale factor was refined during data analysis by matching the telescope and gyroscope signals discussed in the following section.

SQUID to telsecope signal matching
As just pointed out, the estimate of the telescope scale factors C Tx and C Ty relating the bodyfixed pointing error to the normalized telescope signal needs to be refined using on-orbit science data. The body-fixed vector θ ⃗ and the inertial vector in θ ⃗ , appearing in the formula (70) for the GSV pointing, are connected by a rotation due to S/C roll, Using this and denoting O ⃗ the sum of the optical effects in the formula (70), we rewrite the SQUID measurement equation (50) with the body-fixed pointing error entering explicitly: cos , sin . (80) Tx Tx x x Ty Ty y y θ ν θ ν = + = + with ν ξ denoting the noise in the telescope signal channel ξ; the nonlinear correction included in the equation (78) was properly treated in the GP-B estimation process. Apparently, pointing errors can be eliminated from equation (80) by using equation (81) if the ratio of the SQUID and telescope scale factors is known. The procedure for finding this ratio and eliminating , x y θ θ is called SQUID-to-telescope scale factor matching. To implement it, a dither signal, a sinusoid at a known stable frequency, was injected into the attitude control system which forced the S/C to slowly oscillate about the line of sight to the GS. The two dither signals in the x-and y-channels of the attitude control system were at different frequencies denoted dx ω and dy ω . Cross-coupling due the S/C inertia and the attitude control system response made each of the dither frequencies appear in both channels, so the 'true' pointing errors t ( )  with the noise dropped from all the formulas, to save space. The subtraction method of matching (see [36]) consists of two stages. At the first stage we analyze the SQUID and telescope signals separately, finding the amplitudes at the cosines and sines of the dither frequencies in each signal. Then, using the dither component as a benchmark for the whole pointing error, we find such a linear combination of the SQUID and telescope signals which is free of it.
To complete the first stage, we use the model (83) to analyze the SQUID signal without modeling pointing errors themselves except for their dither part. We find thus the estimates of all the dither frequency amplitudes,  We then move on to the second stage of the subtraction method by looking for such a linear combination of the three signals, g g T x T y match 1 2 = + + that does not contain the dither terms, and, presumably, the pointing errors as well. Using equations (83), we require the coefficients at the two cosines and two sines of the dither frequencies to vanish. We thus obtain the following four linear equations for the two unknowns Q 1,2 (in fact, they are simply the negatives of the scale factor ratios): The described matching procedure was actually used at the preliminary stage of GP-B data analysis carried out by means of the 2-floor filter ( [7], section 5.3). In it, the model for finding the combination of signals free of pointing errors was augmented by several other frequencies known to be present in the pointing signal. Those are the harmonics of the roll, roll plus or minus twice orbital frequency, and some others. In this way both the injected and natural dithers were exploited to achieve the most accurate estimate of the scale factor ratios.
The final GP-B results were obtained using the 2 s estimator ( [8], section III) where the scale factor matching was not exploited: the algorithm used the ready pointing signal prepared by the recipes of section 6.1 (see also [7], section 5). However, the matching procedure not only worked nicely at preliminary analysis stages, but was also fundamentally important for validating the whole data analysis process.

Summary of science signal models
The description of all the models needed in the nonlinear multiparameter estimation problem of the GP-B data analysis is complete; here we summarize.
The reference frames needed for the model construction were defined in section 2, where we also discussed the GR predicitions of the orbiting gyroscope precessions, see sections 2.2.1 and 2.2.3. In sections 2.2.2 and 2.4 two-dimensional inertially fixed vectors in the plane perpendicular to the direction to the GS, which play the key role in our modeling, were defined. Those are: s t ( ) ⃗ , representing the gyro drift, t ( ) τ ⃗ , giving the S/C deviation from the direction to the GS, and their difference, the spin-to-roll misalignment vector t ( ) μ⃗ . The main measurement equation, i.e., the model of the LF SQUID signal, is given by the formula (50). Its submodels are: (a) SQUID scale factor model, equations (47), (41), (45) and (46), and bias model, equation (51) during GSV, are calculated by equations (72)-(74), respectively. Pointing errors, being also a part of t ( ) τ ⃗ , were removed from the SQUID signal using the normalized telescope signal whose model is given by the equation (78). The removal is possible due to the SQUID-to-telescope scale factor matching described in section 6.3; the resulting SQUID signal model for GSV free of pointing errors is given by equation (80)  Note that all the models, with two minor exceptions, were carefully derived from the underlying physics, and their complete derivations are given in the proper sections and the three appendices. In just two exceptional cases (see This expression is in the proper frame and is ready for the TF calculation.
A.2. TF through the pick-up loop. HF SQUID signal Let R be the pick-up loop radius. To get the total trapped magnetic flux, Tf Φ , through the loop we need to integrate the magnetic field normal to it over the circle r R ″ < in the plane z 0 ″ = .
Technically it is more convenient to find the equal flux through the upper hemisphere S r R { : , 0 Using (A.11) and integrating over the azimuthal angle The values of the complex harmonic amplitudes H t H t ( ) ( ) n n * = − are clearly seen from (A.14). However, we need to verify the second representation of the TF signal from (34) for our case of small spin-to-pick-up loop misalignment β. To do this, we use formula (4.30) from [46] which gives: Keeping only linear and quadratic terms in all the small angles involved, we obtain: We do not give yˆb Δ because it is not involved in n Δ , which we only need for the answer. Using formula (29) for normal in the body-fixed frame, we obtain where the last representation is true because, by (B.2), x b Δ itself is proportional to ν. In addition, zˆb Δ is needed only to linear order, due to the small factor α in front of it in the utmost right of the previous formula.
With that in mind, we use formulas (18) The inertial expression for n Δ follows from introducing these results to the equation (B.3). In the process we drop a lot of terms because in the correction s n ·Δβ Δ = we multiply n Δ by the spin whose GS component is of order unity, and the other two are small. Therefore the combination of (B.3) and (B.4) is remarkably simple: n e ê cosˆquadratic terms in , , orthogonal toˆ; (B.5) r gs gs Δ ν δ α τ ν = − + the dropped terms make at most a cubic contribution to the angle β. This is the needed GSI correction to the inertial expression (33) for the pick-up loop normal. To avoid any misunderstanding, let us emphasize that the vectors x ẑ ,b b , and n are fixed in the spacecraft body, so they are, of course, the same for both GSV and GSI. However, they are different for GSV and GSI in the inertial space due to the different ATC performance in keeping the S/C attitude. Exactly this difference is seen in the formulas (B.4), (B.5) (as usual, ATC stands for the attitude and translation control system).
By (B.5) the GSI correction to the GSV value of angle β is: The coefficients a b , nj nj ± ± are completely determined by the patch patterns on the rotor and its housing and the rotor asymmetry parameter Q.
We now introduce expressions (C.25) to the drift rate formulas (C.22), move the cosine and sine of roll inside the sums over n, and convert the products of trigonometric functions into sums. This leads to the following expressions:

C.5. Gyro motion through a roll-polhode resonance
Let us finally clarify the picture of gyro motion through a roll-polhode resonance assuming, first of all, that other drifts are negligible as compared to the one generated by the resonance. Then all terms but the resonance term (number m, say,) can be dropped from the gyro motion equations (58), which become: