% IMPORTANT NOTES:
% - figures avaiable on request from the authors

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  THIS IS THE REVISED VERSION OF THE LAGEOS SPIN-ORBIT PAPER BY
%  FARINELLA ET AL., SUBMITTED TO JGR - 28/11/1995
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\documentstyle[12pt,agums]{article}
 \documentstyle[agupp]{article}
%\documentstyle[jgrga]{article}

%\tighten

\received{Receipt date to be added later}
\revised{Revision date to be added later}
\accepted{Accept date to be added later}
\journalid{Vol. to be added later}{Journal date to be added later}
\articleid{Start page to be added later}{End page to be added later}
\paperid{Paper identification to be added later}
\cpright{AGU}{1995}
\ccc{Copyright code to be added later}
   % \cpright{PD}{1993}
   % \cpright{Crown}{1993}
   % (Crown copyrights have no "\ccc{}" information.)
\lefthead{FARINELLA ET AL.}
\righthead{LAGEOS' SPIN--ORBIT EVOLUTION}
\authoraddr{P. Farinella, Dipartimento di Matematica,
 Universit\`a di Pisa, Via Buonarroti 2, I--56127 Pisa, Italy
 (e-mail: paolof{@}dm.unipi.it).}
\authoraddr{D. Vokrouhlick\'y, Institute of Astronomy, Charles
 University, \v Sv\'edsk\'a 8, CZ-15000 Prague 5, Czech Republic
 (e-mail: davok{@}aci.cvut.cz, david{@}ocar01.obs-azur.fr).}
\authoraddr{F. Barlier, Observatoire de la C\^ote d'Azur, Dept. CER\-GA,
 Avenue N. Copernic, F-06130 Grasse, France (e-mail:
 barlier{@}ocar01.obs-azur.fr).}

\slugcomment{{\it Journal of Geophysical Research},
 1995, in press.}
% The \slugcomments appear only in manuscript and preprint styles.

\begin{document}

\title{ The rotation of LAGEOS and its long--term semimajor axis
 decay: A self--consistent solution}

\author{Paolo Farinella}
\affil{ Gruppo di Meccanica Spaziale,
Dipartimento di Matematica, Universit\`a di Pisa, Pisa, Italy}
\author{David Vokrouhlick\'y\altaffilmark{1} and Fran{\c c}ois Barlier}
\affil{ Observatoire de la C\^ote d'Azur, D\'epartement
 CERGA, Grasse, France}
\altaffiltext{1}{On leave from the Institute
 of Astronomy, Charles University, Prague, Czech Republic.}

\begin{abstract} We develop a self--consistent model for the evolution
of the spin axis of LAGEOS and the related long--term semimajor axis
perturbations, due to asymmetric emission/reflection of radiation from
the satellite's surface.  We show that the theory developed by
Bertotti and Iess (1990) for the evolution of LAGEOS' rotation under
magnetic and gravitational torques, which we have somewhat generalized
here, can lead to a successful fit of the observed semimajor axis
residuals, provided the correct initial conditions for the direction
of the spin axis are chosen.  The remaining residuals have a r.m.s.
dispersion of $0.50\times 10^{-12}$~m/s$^2$, comparable to that of
previous solutions, based on purely empirical fits of the spin axis
direction as a function of time.  The spectrum of the residuals
indicates that they are probably due to unmodelled radiation forces
(e.g., from Earth albedo and/or penumbra passages).  Our solution
allows us to predict the future evolution of LAGEOS' rotation for
about another decade in the future, until the spin rate will become so
slow that some basic assumptions of the theory will fail.  A similar
model can also be used to model/predict the coupled spin--orbit
evolution of the LAGEOS~II satellite, launched in late 1992, although
the available data still cover a span of time too short for reaching
significant quantitative conclusions concerning this satellite.

\end{abstract}

\section{Introduction} The diverse non-gravitational mechanisms which
perturb the orbits of artificial satellites have been the subject of
active research since the beginning of the space age.  In the last 15
years, the availability of very accurate laser ranging data for the
passive satellite LAGEOS (more recently joined by other satellites of
the same type) has given rise to the challenging task of understanding
and modelling quantitatively the long--periodic behaviour of the
orbital elements as inferred from the observations.  Standard
dynamical models, which include all the main gravitational and
non--gravitational perturbing forces, are not enough when
accelerations of the order of $10^{-12}$~m/s$^2$ are involved.  A
number of subtle mechanisms, mostly related to the interaction of the
satellite with particles present in space and with electromagnetic
radiation, have been discussed over the years as possible sources for
the residuals appearing in the evolution of the orbital elements, in
particular those in semimajor axis, which result into an accumulating
drift of the satellite's along--track position (for reviews, see
\markcite{{Milani et~al.,} [1987]}, \markcite{{Tapley et~al.,}
[1993]}).

The situation is complicated by the fact that some of the radiation
forces depend on the rotational state of the satellite.  As first
pointed out by \markcite{{\it Rubincam} [1987, 1988]} ``thermal
drag/thrust'' (also named ``Yarkovsky'' and ``Yarkovsky--Schach'')
effects arise as a consequence of thermal radiation anisotropically
emitted from the satellite's surface.  This results in a pertubing
force directed along the rotation axis of the satellite, which may
explain some $70\%$ of the secular part (about $-3.4\times
10^{-12}$~m/s$^2$) of the observed along--track residual acceleration
of LAGEOS [\markcite{{\it Rubincam,} 1990}].  When variations in the
thermal state of the satellite due to the occurrence of eclipses are
also taken into account [\markcite{{\it Afonso et~al.,} 1989}], the
observed periodic acceleration spikes can also be explained by this
mechanism.

However, there are indications that other unmodelled or poorly
modelled effects, such as neutral and charged particle drag
[\markcite{{\it Afonso et~al.,} 1985}; \markcite{{\it Rubincam,}
1990}], optical anisotropy of the satellite's surface [\markcite{{\it
Scharroo et~al.,} 1991}], Earth albedo radiation pressure
[\markcite{{\it Anselmo et~al.,} 1983}; \markcite{{\it Barlier
et~al.,} 1986}; \markcite{{\it Vokrouhlick\'y et~al.,} 1993};
\markcite{{\it Vokrouhlick\'y and Farinella,} 1995}], and solar
radiation pressure during penumbra passages [\markcite{{\it
Vokrouhlick\'y et~al.,} 1994}] do provide minor but significant
contributions to the observed residuals.  Assessing the importance of
these effects is difficult, unless a ``nominal'' model for the main
thermal thrust effects is available. However, for this purpose a
prerequisite is the development of a ``nominal'' scenario for the
evolution of LAGEOS' rotational state.  The same requirement applies
to the interpretation of the observed acceleration residuals for the
recently launched LAGEOS~II satellite.  In principle, this would allow
one to test the dynamical models developed for LAGEOS with a twin
satellite having a different orbital configuration [\markcite{{\it
Anselmo and Pardini,} 1992}].

The dependence of thermal thrust forces upon the orientation of the
satellite's rotational axis has allowed \markcite{{\it Rubincam}
[1990]} to solve inversely for the celestial coordinates of the spin
axis, by matching the predictions of his model for the long--term
along--track perturbations to the observed residuals.  He concluded
that an almost--constant orientation of the spin axis could explain
well the data, with the best--fitting solution for the polar
coordinates (colatitude $\theta=22\deg$ and longitude
$\lambda=313\deg$ in an equatorial reference system) corresponding
well to the telemetry data reported by the flight control center at
the time of injection of the satellite in its current orbit
[\markcite{{\it Rubincam,} 1987}].

However, in his 1990 paper {\it Rubincam} did not mention a second
{\it equally possible} solution for the initial orientation of LAGEOS'
spin axis: $\theta=158\deg$ and $\lambda=133\deg$.  This
``Southern'' solution is related to the previous one by a simple
spatial inversion, with the satellite spinning about the same axis but
in the opposite sense.  This ambiguity results from the fact that the
thermal thrust effects are invariant under a spatial inversion,
because they just depend on the (assumedly axisymmetric) temperature
distribution on the satellite's surface [see formul\ae\ given in
\markcite{{\it Rubincam,} 1987, 1988}; \markcite{{\it Farinella
et~al.,} 1990}].  It is also worth mentioning that the ``Southern''
solution cannot be excluded on the basis of the available information
on LAGEOS' in--orbit injection --- on the contrary, \markcite{{\it
Anselmo and Pardini} [1992]} argued that it is more consistent than
the ``Northern'' one with the engineering data on the apogee kick
motor separation.

In our opinion, the fact that the ambiguity in the spin axis direction
derived by {\it Rubincam} was overlooked in most subsequent studies
--- so that the ``Northern'' solution only was used in modelling the
thermal thrust forces --- has been the source of some
misunderstandings.  In the 1990s, the availability of a longer time
series for the acceleration residuals has led to the conclusion that
the assumption of a fixed spin axis is no longer sufficient to match
the data [\markcite{{\it Ries et~al.,} 1993}]. \markcite{{\it Rubincam}
[1987]} himself had
suggested that magnetic torques associated with eddy currents in the
satellite body should gradually change the direction of LAGEOS' spin
axis, and that this should be reflected in the resulting thermal
thrust force.  In 1991, {\it Bertotti and Iess} published a
theoretical model for the evolution of LAGEOS' rotational state
(hereinafter referred to as the BI model), according to which the spin
rate would rapidly decrease under the dissipative effect of induced
currents, and the direction of the spin axis would undergo a complex
evolution pattern due to the interplay of magnetic and gravitational
torques (the latter being due to the nonspherical mass distribution of
the satellite).  A crucial feature of the BI theory was that the
evolution of both the spin rate and the spin axis direction is {\it
not} invariant with respect to the inversion of the initial sense of
rotation --- actually, {\it Bertotti and Iess} suggested that some of
the available measurements of LAGEOS' spin rate after launch allow one
to discriminate between the two sets of initial conditions, favouring
the ``Southern'' solution.

A different approach to the problem was that of developing
\markcite{{\it Rubincam's} [1990]} method and deriving a best--fit
solution for the
spin axis direction from the observed along--track residual
accelerations, by solving for the corresponding parameters not just
once, but over a sequence of successive time intervals.  Of course,
this method does not allow one to discriminate between the two senses
of rotation, due to the above--mentioned invariance of the thermal
thrust effects.  The successful application of this technique in
recent years [\markcite{{\it Ries et~al.,} 1993}; \markcite{{\it Robbins
et~al.,} 1994}; \markcite{{\it Kolenkiewicz et~al.,} 1995}]
and the poor agreement of the best--fitting
spin axis evolution with the predictions of the BI theory --- {\it
based on the Northern initial conditions} --- led to the conclusion
that the theory was either incorrect or incomplete, possibly due to
some (vaguely defined) ``chaotic'' character of the solutions of the
dynamical equations.  Even a recent re-analysis of LAGEOS' rotational
dynamics, where some assumptions of the BI theory were relaxed
[\markcite{{\it Habib et~al.,} 1994}], assumed the ``Northern''
initial conditions to predict the future evolution.

Yet, a self--consistent solution for LAGEOS' orbital along--track
perturbations and spin axis evolution based on a theoretical model for
the latter would be highly desirable, even when compared to the
``empirical'' best--fit solutions.  There several reasons for this
preference:
\begin{description}

\item{(i)} the future evolution  could be predicted with some
confidence, including in particular the qualitatively important
transitions to a slow--rotation regime characterized by a
non--axisymmetric surface temperature distribution and to a (possibly
chaotic) spin--orbit resonant state;

\item{(ii)} better constraints could be imposed on some poorly known
physical parameters appearing in the expression of the thermal
thrust perturbations;

\item{(iii)} once a fit of the theoretical predictions with the
observed residuals had been performed, the post--fit residuals could
be interpreted in terms of perturbation mechanisms not included in the
dynamical model (variable particle drag effects, Earth albedo
radiation pressure, etc.).
\end{description}

In this paper, we show that the BI model is indeed fully suitable
for recovering the basic features of the observed along--track orbit
residuals of LAGEOS, provided only the ``Southern'' initial
conditions are selected instead of the ``Northern'' ones.  Our
approach consists in modelling only the perturbation forces which are
related to the orientation of the satellite's spin axis, in such a way
that a small number of {\it a priori} unknown free parameters is
required and the post--fit residuals can be attributed to the
remaining non--gravitational forces mentioned earlier.

The resulting ``nominal solution'' for the evolution of LAGEOS'
rotation may be also valuable for the analysis of the observed
residuals in orbital elements other than the semimajor axis, which are
also affected by thermal thrust effects [see e.g. \markcite{{\it
Farinella et~al.,} 1990}].  In a forthcoming paper, we plan to analyse in
particular the eccentricity residuals [or equivalently their
``excitations''; see \markcite{{\it Ries et~al.,} 1993}], which show
a puzzling
behaviour.  Subtracting the thermal thrust effects estimated from the
semimajor axis data, the contribution of additional mechanisms
affecting the eccentricity perturbations might become easier to detect
and to interpret.

The remainder of this paper is organized as follows.  In the next
Section, we summarize the BI theory and derive a simple generalization
of it allowing for a misalignment of the rotational angular velocity
vector with respect to the satellite's axis of symmetry.  This does
not modify the qualitative features of the theory, but just affects
the value of an important dynamical parameter appearing in it.  Then,
we proceed to present our best--fit simultaneous solution for LAGEOS'
long--term along--track perturbations and its spin axis evolution, and
discuss the origin of the post--fit residuals; the same approach is
applied to the LAGEOS~II data as well.  In the final Section we
summarize the main conclusions of this work and discuss some possible
future developments.

\section{ Spin axis evolution model}
In this section, we describe our version of the theoretical
model for LAGEOS' spin axis evolution developed by \markcite{{\it
Bertotti and Iess} [1991]} (BI model), which we have implemented in
a computer
program after correcting some errors in the published formulae.  We
have also introduced a generalization of the BI model to the case when
the satellite's spin axis is not aligned with its axis of symmetry.
In what follows, we summarize the main ideas and assumptions
behind this theory, and refer to the original paper of {\it Bertotti
and Iess} for more details.

Whenever external torques are applied, the satellite's rotational
angular momentum ${\bf L}$ undergoes an evolution described by Euler's
kinematical and dynamical equations.  However, it is not convenient to
consider the full set of six dynamical equations for the satellite's
rotational state, as the individual parameters describing it change
significantly on different time scales.  For instance, the initial
rotation period of LAGEOS was about $0.55$~s, whereas the evolution of
${\bf L}$ is characterized by a time scale of the order of one
year.  If one neglects the occurrence of resonances, which do not take
place in the case of LAGEOS for at least three decades after launch,
one can apply the usual technique of averaging over the fast variables
describing the proper rotation and the precession of the angular
velocity vector, in analogy with the procedure adopted for the orbital
(translational) motion.  For more details on this method, we refer to
\markcite{{\it Beletskii} [1966]} and \markcite{{\it Fitzpatrick} [1970]}.
As a result, one
is left with four dynamical equations: three for the components of
${\bf L}$ and one for the angle $\vartheta$ giving the tilt of the angular
velocity ${\bf \omega}$ from the symmetry axis of the satellite (which is
assumed to have an axisymmetric mass distribution, with principal
moments of inertia $C$ and $A<C$).

In agreement with the BI model, we shall assume that two different
kinds of external torques affect LAGEOS' rotation:
\begin{enumerate}

\item {\it Magnetic torques}, due to eddy currents induced in the
satellite body by LAGEOS' fast rotation within the Earth's magnetic
field;

\item {\it Gravitational torques}, similar to the lunisolar
tor\-ques on the Earth, due to the (small) nonsphericity of LAGEOS'
body. \end{enumerate}

\noindent Magnetic torques can be expressed as
\begin{equation}
 {\bf M} = {\bf m} \times {\bf B} \; ,
\end{equation}
where ${\bf m}$ is the satellite's magnetic moment and ${\bf B}$ the
intensity of the geomagnetic field.  In the following, we shall take
into account only the main, inclined dipole part of the Earth's
internal magnetic field, having a tilt angle $D\approx 11\deg$ to
the
Earth's rotation axis, and neglect both the higher--pole terms of the
internal field and the magnetic field due to ionospheric currents
(more on this at the end of this Section).  The precise evaluation of
the induced magnetic moment ${\bf m}$ for an axially symmetric conductor
rotating in an external magnetic field ${\bf B}$ is a difficult problem.
However, a good approximation is
\begin{equation}
 {\bf m} = V\alpha' {\hat {\bf \omega}} \left({\hat {\bf \omega}}
 \cdot {\bf B}\right) - V\alpha'' {\bf B}\times{\hat {\bf \omega}}
\end{equation}
[\markcite{{\it Beletskii and Khentov,} 1985}; \markcite{{\it Bertotti
and Iess,} 1991}].  Here ${\hat {\bf \omega}}$ is the unit vector
directed along  is the satellite's rotation rate
vector (${\bf \omega}$), $\alpha'$ and $\alpha''$ are the average
polarizability
parameters per unit volume of the conductor, and $V$ is its volume.
Let $\delta$ be the $e$---fold depth for the eddy current penetration
into the conductor:
\begin{equation}
 \delta = {c \over \sqrt{2\pi \omega\sigma}} \; ,
\end{equation}
where $\sigma$ is the electrical conductivity and $c$ is the speed of
light.
Then, if $\rho$ is the satellite's radius, in the low--frequency limit
$\rho \ll \delta$ the polarizability parameters can be approximated by
\begin{mathletters}
\begin{eqnarray}
 \alpha' & = & -{1 \over 105\pi} \left({\rho \over
 \delta}\right)^4 \beta'\left[\left({\rho \over \delta}\right)^2
 \right] \nonumber \\
 &=& -{1 \over 105\pi} \left({\rho \over \delta}\right)^4
 \beta'\left(0\right) + {\cal O}\left[\left({\rho \over \delta}\right)^6
 \right] \; , \\
 \alpha'' &= & {1 \over 20\pi} \left({\rho \over
 \delta}\right)^2 \beta''\left[\left({\rho \over \delta}\right)^2
 \right] \nonumber \\
 &=& {1 \over 20\pi} \left({\rho \over \delta}\right)^2
 \beta''\left(0\right) + {\cal O}\left[\left({\rho \over \delta}\right)^4
 \right] \;
\end{eqnarray}
\end{mathletters} \noindent
[\markcite{{\it Landau and Lifschitz,} 1960}; \markcite{{\it Bertotti
and Iess,} 1991}].  Here $\beta'$ and $\beta''$ are dimensionless
functions of the
$\rho/\delta$ ratio specifying how the electrical properties of the
satellite differ from those of a homogeneous sphere.  In principle,
the parameters $\beta'(0)$ and $\beta''(0)$ are not known {\it a
priori}, and should be adjusted {\it a posteriori} from suitable
measurements of the spin rate evolution.

A possible future improvement of the BI magnetic torque model might
consist
in taking more terms in the series development of the magnetization
parameters $\beta'$ and $\beta''$ appearing in Eqs.~(4).  The reason
is that the $\rho/\delta$ ratio changes significantly during the
mission, due to the decreasing spin rate.  We suggest this extension
of the model rather than the application of the original formul\ae\
given by \markcite{{\it Landau and Lifschitz} [1960]}.
The latter are valid for
any value of the $\rho/\delta$ parameter, but apply only to the case
of a homogeneous satellite with isotropic polarizability properties.
Of course, these conditions are not satisfied for the real LAGEOS
satellite.

Eqs.~(1) to (4) provide the instantaneous values of the magnetic
torques generated in the satellite.  Applying the averaging procedure
mentioned earlier, the two terms appearing in Eq.~(2) can be shown to
produce the following contributions to the overall torque (1):
\begin{eqnarray}
 M^{(1)}_i &=& V {\alpha'(\omega) \over \omega^2} \epsilon_{ijk}
 \omega_j \langle B_k B_l\rangle \omega_l \nonumber \\
 & \equiv& V {\alpha'(\omega) \over \omega^2} {d^2 \over a^6}
 \epsilon_{ijk} \omega_j \alpha_{kl} \omega_l \; , \\
 M^{(2)}_i &=& -V {\alpha''(\omega) \over \omega}
 \omega_j \langle B^2\delta_{ij} - B_i B_j\rangle \nonumber \\
 &\equiv&  -V {\alpha''(\omega) \over \omega} {d^2 \over a^6}
 \omega_j \beta_{ij} \; ,
\end{eqnarray}
where $d$ is the Earth's magnetic dipole, $a$ is the orbital
radius, $\epsilon_{ijk}$ denotes the fully antisymmetric Levi-Civita
tensor, and $\langle ... \rangle$ means that an average has been taken
over (i) the proper rotation and precession of the satellite, (ii) the
satellite's revolution around the Earth, and (iii) the Earth's
rotation period (which also gives the periodicity of the the spinning
magnetic dipole in inertial space).  A straightforward calculation
yields
\begin{mathletters}
\begin{eqnarray}
 \alpha_{ij} &= & {9\over 8} \left[\sin^2I + {1 \over 2}
 \sin^2D\left(20-27\sin^2I\right)\right] \delta_{ij} \nonumber \\
 & & - {1 \over 8}\left(1-3\cos^2D\right) E_i E_j \nonumber \\
 & & +{3 \over 8}\cos I\left(1-3\cos^2D\right) \left(E_i c_j + E_j
  c_i\right) \nonumber \\
 & & - {3 \over 8} c_i c_j \biggl[3\left(1-3\cos^2I \right)
  \nonumber \\
 & & \qquad\quad\;\;\;\; -{1\over 2}\sin^2D\left(5-27\cos^2I\right)
  \biggr] \; , \\
 \beta_{ij} &= & \biggl\{1+{1 \over 8}\Bigl[3\sin^2I \nonumber \\
 & & \qquad\quad + {1 \over 2} \sin^2D\left(4-9\sin^2I\right)\Bigr]
 \biggr\}\delta_{ij} \nonumber \\
 & & +  {1 \over 8} \left(1-3\cos^2D\right) E_i E_j \nonumber \\
 & & -{3 \over 8}\cos I\left(1-3\cos^2D\right) \left(E_i c_j + E_j
  c_i\right)\nonumber \\
 & & + {3 \over 8} c_i c_j \biggl[3\left(1-3\cos^2I \right)
  \nonumber \\
 & & \qquad\quad\;\;\;\; -{1\over 2}\sin^2D\left(5-9\cos^2I\right)
  \biggr] \; .
\end{eqnarray}
\end{mathletters}
Here ${\bf E}$ is the $(0,0,1)$ unit vector, $\delta_{ij}$
is Kro\-ne\-cker's delta and ${\bf c} =
(\sin I\sin\Omega,- \sin I\cos\Omega,\cos I)$ is the unit vector
perpendicular to the satellite's orbit, whose inclination and nodal
longitude are denoted by $I$ and $\Omega$. Note that formula
(7b) corrects an error in Eq.~(A4) of \markcite{{\it Bertotti and Iess}
[1991]}.

According to {\it Bertotti and Iess}, the first term in Eq.~(2) causes
only a small additional precession of the satellite.  It is merely the
second term which is responsible for the major effect resulting in the
evolution of LAGEOS' rotation, namely the fast despinning process.  As
a consequence, only the value of the $\beta''(0)$ parameter can be
constrained effectively by the observed spin rate.  Hereinafter, we
shall thus include in the dynamical model only the ${\bf M}^{(2)}$
torque given by Eq.~(6).  A simple analysis suggests that the
characteristic time scale of the magnetic despinning process is given
by
\begin{equation}
 \nu_m = V {\alpha''(\omega) \over \omega} {d^2 \over a^6} =
  V {\alpha''(\omega) \over \omega} B_0^2 \left({R_\oplus \over a}
  \right)^6 \; ,
\end{equation}
where $B_0\approx 0.307$~G is the intensity of the
geomagnetic field at the Earth's surface and $R_\oplus$ is the Earth's
mean radius.  [Note that in the BI
Eq.~(21) for $\nu_m$, one should replace $20\pi$ with $10$ in the
denominator, while the ratio $a^2/d^6$ should be replaced by
$d^2/a^6$.]  For LAGEOS, $\nu_m^{-1} \approx 0.87\, {\rm yr}$.

Finally, in order to set up the differential equations of the
averaged model one has to choose a suitable reference frame [see {\it
Beletskii}, 1966 for details].  We shall use an inertial equatorial
system ($z$--axis directed to the North pole, $x$--axis pointing to
the vernal equinox at a specified date).  The evolution of the angular
momentum vector in this reference frame is given by
\begin{eqnarray}
 {d {\bf L} \over dt} &=& {\bf M} + {3n^2 \over 4L^2} \left(C-A\right)
 \left(3\cos^2\vartheta - 1\right) \nonumber \\
 & & \qquad\quad\quad\quad \left({\bf c} \cdot {\bf L}\right)
  \left({\bf c} \times {\bf L}\right)\; ,
\end{eqnarray}
where ${\bf M}$ now stands for the averaged magnetic torque
acting on the satellite and $n$ is its orbital mean motion.  The
second term corresponds to the gravitational torque of the Earth's
(monopole) gravity field, averaged over the proper rotation and
precession of the satellite and the orbital revolution around the
Earth [see e.g. \markcite{{\it Beletskii,} 1966}; \markcite{{\it
Fitzpatrick,} 1970}].  Note
that this term takes the same form as the corresponding term in the BI
Eq.~(48), provided that $\vartheta \equiv 0$.

The system of equations (9) must be complemented by the differential
equation for the evolution of the $\vartheta$ tilt angle.  Let us
introduce an auxiliary reference frame ($R'$) whose
$z$--axis is directed along the satellite's angular momentum vector,
while the $xz$--plane coincides with the plane defined by the angular
momentum vector and the vector perpendicular to the osculating
satellite orbit [see \markcite{{\it Beletskii,} 1966}].  Then, one obtains
\begin{equation}
 {d(\cos\vartheta) \over dt} = \nu_m \Delta {\hat \beta}_{12}
 \cos\vartheta \sin^2\vartheta + {\cal O}(\nu_m \Delta^2) \; ,
\end{equation}
where ${\hat \beta}_{12}$ is the $(1,2)$ component of the
$\beta$ matrix transformed to the $R'$ frame and $\Delta = (C-A)/C$ is
the usual nondimensional ``flattening parameter'' defined by the
principal moments of inertia $A$ and $C$.  Two comments are
appropriate for Eq.~(10).  First, it admits the trivial solution
$\vartheta =0$, as assumed in the BI treatment.  Second, note that
$\vartheta$ changes only on a very long time scale, of order $(\nu_m
\Delta)^{-1}$.  Therefore, on the time span for which our LAGEOS
evolution model can be applied, roughly three decades (see later), the
$\vartheta$ parameter to a good approximation can be taken as
constant.  In the solutions, it will then play the role of rescaling
the satellite's flattening parameter $\Delta$ to the ``effective
value'' $\Delta' \equiv (3 \cos^2\vartheta - 1)\, \Delta /2\, \, \leq
\Delta$.  A determination of $\vartheta$ might be important for the
future evolution of LAGEOS' rotation, in particular because the timing
at which the spin--orbit resonance state will be met depends on this
parameter.  Moreover, the thickness of the chaotic zone, associated
with the separatrix of this resonance, may depend sensitively on the
dynamical flattening parameter $\Delta'$.

The mathematical setting of the theory is completed by the
relationship between the angular momentum vector ${\bf L}$ and the
rotational angular velocity vector ${\bf \omega}$.  Averaging over a
rotational and a precessional cycle, one obtains
\begin{equation}
 {\bf L} \cong C \left(1-\Delta \sin^2\vartheta\right) {\bf \omega}
 \; ,
\end{equation}
which generalizes the BI relationship ${\bf L}=C{\bf \omega}$, consistent
with the $\vartheta=0$ assumption.  Note that Eq.~(11) holds only
neglecting terms of order ${\cal O}(\Delta^2)$.

The general behavior of LAGEOS' spin axis evolution has been already
discussed by {\it Bertotti and Iess}.  At the beginning, due to the
fast rotation, magnetic torques play the dominant role, causing a
tendency to align the rotation axis with the magnetic pole.  The
precessional angular velocity $\omega_p$ due to the gravitational
torque is:
\begin{equation}
 \omega_p = {3n^2 \over 2 L^2} (C-A) \left({\bf c} \cdot {\bf L}
 \right) \; .
\end{equation}
The inequality $\omega_p \leq \nu_m\Delta \approx {\dot
\Omega}\Delta$ holds at the beginning of the LAGEOS mission.  Later
on, when the rotation has been slowed down by dissipation of energy
via eddy currents, a complicated intermediate phase takes place for a
few years, with a comparable influence of the magnetic and
gravitational torques.  At the end, the rotation of the satellite
becomes so slow that the magnetic processes are negligible, and a
longitudinal precession is triggered by the gravitational torque
alone.  We will see some more details on the specific evolutions of
LAGEOS and LAGEOS~II in the next Section.

Let us now briefly discuss some of the simplifications involved in
this theory and the limits of their applicability.  First, the
averaging of $\beta_{ij}$ has been performed over one orbital
revolution, assuming the orbit as circular.  This corresponds to
neglecting the small eccentricity ($e=0.0044$) of LAGEOS' orbit.
However, it can be shown that the neglected terms in $\beta_{ij}$ are
only of the second order in the eccentricity, so this approximation is
justified.

A potentially more serious problem is due to the neglected quadrupole
part of the Earth's magnetic field, whose amplitude reaches $\approx
10\%$ of the main dipole field [\markcite{{\it Peddie and Fabiano,}
1982}; \markcite{{\it Bertotti and Farinella,} 1990}].  In a more
complete theory, at lowest order the Earth's magnetic quadrupole
should ap\-pear in the ex\-pres\-sion of the mag\-netic torque (which
is quadratic in the magnetic field ${\bf B}$) owing to a coupling with
the contribution of the magnetic dipole to the $\beta_{ij}$ quantity.
However, after some algebra, it can be seen that these coupling terms
vanish, when averaged over one (assumedly circular) orbit.  Thus the
largest non--vanishing terms due to the quadrupole in the torque are
caused by the coupling with the quadrupole terms in $\beta_{ij}$.
Hence they are smaller than $1\%$ than the contribution of the main
dipole field, and can be safely neglected.  One can also neglect the
induced magnetic field due to the ionospheric currents, which is
typically less than $1\%$ of the main field.

Finally, it does not appear likely that other types of torque play an
important role in the evolution of LAGEOS' rotational state.  As for
torques due to atmospheric skin forces and radiation pressure, an
estimation of their magnitude can be based on the formul\ae\ given by
\markcite{{\it Beletskii} [1966]} and \markcite{{\it Beletskii et~al.}
[1975]}.  The radiation pressure torque associated to the (assumed)
anisotropic reflectivity of LAGEOS' surface [\markcite{{\it Scharroo
et~al.,} 1991}] does not exceed a few percent of the gravitational
torque.  Another additional radiation pressure torque arises if
LAGEOS' centre of mass is shifted by some small amount with respect to
the geometric centre of the satellite ({\it V.~Slabinski}, personal
communication).  This might reach at most $\approx 10\%$ of the
gravitational torque (for a shift of several mm), so it is not likely
to affect drastically the spin evolution.  However, this effect is not
completely negligible, and although no independent data exist on the
size of the shift, we plan to further study this problem.

\section{ Along--track perturbations from thermal thrust forces}
A great deal of work has been devoted in the last 15 years
to modelling a variety of perturbing forces on LAGEOS which may yield
detectable along--track (i.e., semimajor axis) effects.  We refer to
the literature mentioned in the Introduction for more details.  For
the reasons
dicussed there, in this paper we are going to limit ourselves to model
the effects of those forces which depend on the satellite's rotational
state.  They are:
\begin{enumerate}

\item the so--called Yarkovsky ther\-mal thrust (or drag) effect,
 due to the recoil force asso\-ci\-ated with an\-iso\-tro\-pic thermal
 radiation
 emitted from the satellite's surface, asymmetrically heated by the
 incoming flux of infrared radiation from the Earth [\markcite{{\it
 Rubincam,} 1987, 1988}];
\item the so--called Yarkovsky--Schach effect, similar to the
 previous one but due to asymmetric heating of the satellite by the
 (visible) solar radiation flux, which undergoes interruptions when the
 orbit crosses the Earth's shadow, so that the surface temperature
 distribution and the resulting force undergo strong variations during
 the eclipse periods [\markcite{{\it Afonso et~al.,} 1989};
 \markcite{{\it Farinella et~al.,} 1990}];
\item the force resulting from a difference in the reflectivity of
 the two hemispheres of LAGEOS, as defined by its rotation axis, which
is also modulated by eclipses
 [\markcite{{\it Scharroo et~al.,} 1991}];
\end{enumerate}

Also, we added to these perturbations a constant deceleration term,
simulating the secular effect of drag by both charged and neutral
particles present in space at LAGEOS' altitude [\markcite{{\it
Afonso et~al.,} 1985}; \markcite{{\it Rubincam,} 1990}].
Apart from this constant term, all the other effects produce a
net force directed along the spin axis; the two
thermal thrust effects are independent of the sense of rotation,
whereas the force associated with the hemispheric reflectivity
asymmetry
changes sign when the North and South poles of the satellite are
interchanged; but since there is no {\it a priori} knowledge (e.g.,
from pre--launch laboratory measurements) on the real reflectivity
asymmetry for LAGEOS, we just changed the sign in the formul\ae\ given
by \markcite{{\it Scharroo et~al.} [1991]} when taking the opposite
sense of rotation.

The LAGEOS and LAGEOS II data sets of along--track acceleration
residuals which we have used in our analysis have been kindly provided
to us by the Center of Space Research of the University of Texas
(CSR/UT).  For the former satellite, the data (a sequence of average
values over bimonthly arcs) cover $17.9$ years since the first
accurate orbit determination following launch on May 4, 1976.  This is
a time span long enough to include the transition from the initial
period when the magnetic torque dominates the rotational evolution to
the later phase when the gravitational torque becomes more and more
important.  On the other hand, in the case of LAGEOS~II we have at our
disposal data covering only the first $2.2$ years of the mission.  We
will show that this period is still too short to provide accurate
information on the spin axis evolution, and in particular to constrain
effectively some physical parameters appearing in the theory by
adjustment to the observed data.

We have cross--checked our results derived from the CSR/UT data by
applying the same analysis to the independent series of acceleration
residuals provided by the NASA Goddard Space Flight Center (GSFC) at
Greenbelt, Maryland [\markcite{{\it Robbins et~al.,} 1994};
\markcite{{\it Kolen\-kie\-wicz et~al.,} 1995}].  The only significant
difference we have found has been a somewhat larger r.m.s. of the
post--fit residuals with the GSFC data, as we will see below.

We shall now discuss separately our best--fit solutions for the cases
of LAGEOS and LAGEOS~II.

\subsection{ Solution for LAGEOS}
Our solution for the along--track acceleration of LAGEOS vs. time
with the model described above is shown in Fig.~1.
The empty squares correspond to the first three years of data, which
have been considered with a lower statistical weight ($1/3$) due to
the lower accuracy of the laser tracking technique at the end of the
1970's [\markcite{{\it Tapley et al.,} 1993}].  Later on, the typical
error bars to be associated with the data points are of about
$0.3\times 10^{-12}~{\rm m/s}^2$.  Fig.~1 clearly shows that the fit
is fairly good in reproducing all the main observed features.  The
r.m.s. dispersion of the post--fit residuals is $0.50\times
10^{-12}~{\rm m/s}^2$, a value quite similar to those obtained in
previous fits, resulting from a purely empirical adjustment of the
spin axis direction vs. time relationship.  Since only a small number
of free parameters were adjusted in our fit, this result implies that
the BI theory is a powerful tool to model the rotation--dependent
radiation forces acting on LAGEOS.

Besides the initial coordinates of the spin axis, we solved for some
other physical parameters appearing in the force model.  One was the
reflectivity difference between the two LAGEOS hemispheres, for which
we derived a value of $0.014\pm 0.002$, in good agreement with
previous estimates [\markcite{{\it Scharroo et~al.,} 1991};
\markcite{{\it Robbins et~al.,} 1994}].  Other parameters are related
to the Yarkovsky and Yarkovsky--Schach effects.  While we kept the
same values adopted by \markcite{{\it Scharroo et~al.} [1991]} for
their normalizing coefficients (i.e., $3.08$ and $0.89\times 10^{-12}$
m/s$^2$, respectively), we adjusted the thermal lag time of LAGEOS'
retroreflectors ($\tau_{RR}$), which appears in both the force terms
and is often expressed through the lag angle $n\tau_{RR}$.  Our best
estimate for this lag angle was $72\deg$ ($\tau_{RR}=2740$~s), again
close to the values adopted by \markcite{{\it Rubincam} [1988]} and
\markcite{{\it Scharroo et~al.} [1991]} and by \markcite{{\it Afonso
et~al.} [1989]} ($83\deg$ and $87\deg$, respectively).  According
to our tests, the best--fit solution is not very sensitive
to this parameter.  Actually, we noted that the best--fit value of
$\tau_{RR}$ can be varied in a fairly wide range (say, $\pm 20\deg$)
provided one adjusts at the same time the normalizing coefficient of
the Yarkovsky--Schach effect.  In other words, these two parameters
have a strong correlation, a property which may result useful when the
model will be fitted to the residuals of other orbital elements, such
as the eccentricity and the argument of perigee.

For the initial spin period of the satellite and the $\beta''(0)$
parameter appearing in Eq.~($4b$), we kept the $0.55$~s and $0.213$
values determined by \markcite{{\it Bertotti and Iess} [1991]} from a
fit of the predicted spin rate evolution to the measured spin periods
of LAGEOS; although we did not adjust these parameters, we checked
that moderate changes in these values would not affect much our
results.  Finally, we chose the value $-0.95\times 10^{-12}~{\rm
m/s}^2$ for the constant deceleration term assumed to be due to
charged and neutral drag, that is about $25\%$ of the observed secular
effect, in agreement with the findings of \markcite{{\it Rubincam}
[1990]}.

The post--fit residuals are plotted in Fig.~2, while Fig.~3 shows
their power spectrum.  Fig.~2 displays somewhat larger residuals near
the beginning and the end of the time span elapsed since launch.  It
appears likely that at the beginning the poorer observational accuracy
is the main responsible, while at the end the faster movement of
LAGEOS' spin axis (see Fig.~4 below) increases the sensitivity of the
results to small model errors.  It is worthwhile noting that the
dominant spectral lines appearing in Fig.~3 correspond to frequencies
which are very close to linear combinations of the rates of the solar
longitude $\lambda_{\odot}$ and the satellite's nodal longitude
$\Omega$, namely $m({\dot \lambda}_\odot - {\dot \Omega})$, with $m=2$
(period of $280$ days), $m=3$ ($187$ days) and $m=4$ ($140$ days);
also, there is a $221$--days line corresponding to the $(2{\dot
\lambda}_\odot - {\dot \Omega})$ argument.  The amplitude of all these
lines is of the order of $0.2\times 10^{-12}~{\rm m/s}^2$.  We do not
find in the post--fit residuals a line corresponding to the main
($m=1$) $560$--days harmonic, but the corresponding effects have been
probably ``absorbed'' by the fitting procedure, because the
perturbation mechanisms depending on the spin axis direction do
include this periodicity (see e.g. \markcite{{\it Afonso et~al.,}
[1989]}).  The presence of the lines listed above supports the idea
that the post--fit residuals are mainly due to radiation forces not
included in our dynamical model: for instance, radiation pressure
associated to the specular part of Earth--reflected sunlight produces
harmonics with just these frequencies [\markcite{{\it Barlier et~al.,}
1986}; \markcite{{\it Vokrouhlick\'y et~al.,} 1993}].  Recent work by
\markcite{{\it Vokrouhlick\'y and Farinella} [1995]} has shown that,
with realistic models of the wavy geometry of the ocean surfaces, the
amplitude of the along--track perturbations due to Earth--reflected
radiation pressure is quite comparable to that of the post--fit
residuals shown in Fig.~2.  Similar considerations can be done
concerning the ``residual'' force terms due to solar radiation
pressure which arise during the penumbra transitions preceding and
following eclipses, provided the refraction and absorption of sunlight
by the atmosphere are realistically modelled [\markcite{{\it
Vokrouhlick\'y et~al.,} 1994}].

Fig.~4 shows the evolution of the equatorial projection $(s_x,s_y)$ of
LAGEOS' spin axis ($s_x$ and $s_y$ being the corresponding components
of the $\hat{\bf \omega}$ unit spin rate vector).  The adjusted values for the
initial conditions, expressed in terms of polar coordinates, are
$\theta=155\deg$ for the colatitude and $\lambda=133{\deg}.5$ for
the longitude.  The corresponding formal $\sigma$'s
resulting from the fitting procedure are $2{\deg}.5$ and $5{\deg}.7$,
respectively.  In Fig.~4 we can note that the first $11$ years of
spin axis evolution are characterized by a slow latitude drift toward
the South pole.  Then a few wide loops about the South pole follow for
a few years, caused by the interplay of the magnetic and gravitational
torques (which in this intermediate phase have comparable strengths),
and finally there is a phase of rapid longitudinal precession.  In the
future (dashed line in Fig.~4), this precession will approach a
limiting cone of $\approx 20\deg$ half--aperture.  During the whole
evolution, a fast quasi--exponential decay of the satellite's spin
rate always takes place (see Fig.~4 in {\it Bertotti and Iess}, 1991).
We recall that the averaging assumptions implicit in the BI theory
fail once the $1:1$ resonance between the spin period and
the orbital period is reached [\markcite{{\it Habib et~al.,} 1994}].
According to our solution, this will occur for LAGEOS approximately
$33$ years after launch, namely about the year 2009.

Even before that, however, i.e. approximately 10 years from present,
the rotational period of LAGEOS will reach $\approx 3000$~s, that is a
value close to the thermal lag time $\tau_{RR}$.  In this situation,
longitudinal temperature gradients on the satellite's surface will no
longer be negligible, and the thermal thrust force will be neither
aligned with the spin axis nor indenpendent of the sense of rotation.
More complex models will then be needed to interpret the residual
along--track accelerations.  It will be interesting to test the BI
theory in these new conditions.

Our successful application of the BI theory, where we kept the nominal
value of the LAGEOS flattening parameter ($\Delta=0.0335$) reported in
the pre--launch technical documentation [\markcite{{\it Johnson et~al.,}
1976}] and already adopted by \markcite{{\it Bertotti and Iess} [1991]},
suggests that this value is fairly accurate, contrary to the conclusions
of \markcite{{\it
Scharroo et~al.} [1991]}.  A significantly smaller value of $\Delta$
would have delayed the onset of the longitudinal drift of the spin
axis in the early 1990's, which is instead clearly reflected in the
observed residual accelerations.  This also limits the possible values
of the tilt angle $\vartheta$, which was so far unconstrained and
results into a reduction of the ``effective value'' of $\Delta'$.  We
have carried out some tests with nonzero values of $\vartheta$, and
the results suggest that $\vartheta \geq 30\deg$ would not be
consistent with the observed acceleration residuals.  However, the
current data set does not allow us to determine a tighter constraint
on $\vartheta$ at the present time.  In the forthcoming years, such
constraints might become available from new measurements of LAGEOS'
spin period, combined with more residual acceleration data.  On the
other hand, a small value of $\vartheta$ is consistent with the
specifications of the in--orbit injection system, which normally
require $\vartheta$ values smaller than a few degrees ({\it
L.~Anselmo}, personal communication).

Fig.~5 shows some contour lines of constant r.m.s. post--fit residuals
in the plane corresponding to the initial spin axis coordinates.  This
plot has been obtained by keeping all the other free parameters of the
model fixed at their ``nominal'' values mentioned earlier (adjusting
other parameters at the same time as the spin axis position would not
change much the results).  Fig.~5 clearly shows that the goodness of
our fit is much more sensitive (by about a factor 3) to the initial
(co)latitude of the spin axis than to the initial longitude, in
agreement with the formal uncertainties for these parameters as stated
earlier.  Mainly, this is due to simple geometrical reasons: since the
coordinates are not very far from the South pole, a longitude shift
results into a smaller displacement of the spin axis than a similar
latitude shift.

As we already mentioned, all the main results summarized above,
including the initial conditions for the spin axis, have been
confirmed using the GSFC along--track acceleration residuals.  The
discrepancies are always within the stated uncertainties of the
solve-for parameters.  The most significant difference is a slightly
greater value of the post-fit r.m.s. dispersion of the residuals,
about $0.68\times 10^{-12}~{\rm m/s}^2$.  The most likely reason for
this is just a somewhat larger scatter of the GSFC data when compared
to the CSR/UT ones, probably related to the different orbit
determination procedures used by the two groups.

Our best--fitting procedure was based on the idea of adjusting the
constants appearing in our spin axis model to directly fit the time
series of residual along--track accelerations.  An alternative
two--stage methodology might be that of matching instead the time
series of LAGEOS' spin axis directions derived by \markcite{{\it Ries
et~al.} [1993]}, who determined from the orbital data a set of
best--fitting spin axis directions over a sequence of short time
intervals, and fitted these successive positions with a regular curve.
As shown in Fig.~6, the qualitative features of the two best--fitting
evolutions are similar, but there are significant quantitative
differences (i.e., the latter method yields initial polar coordinates
$\theta=162\deg$ and $\lambda=152\deg$).  There are three
possible reasons for this discrepancy: (i) LAGEOS' semimajor axis and
orbital longitude are obtained from the along--track accelerations
through one and two integrations with respect to time, respectively,
so matching the accelerations does not necessarily mean obtaining a
good fit for the osculating elements or the instantaneous position of
the satellite (this was independently pointed to us by {\it
P.~Exertier} and {\it A.~Milani}); (ii) on the other hand, the series
of spin axis positions derived by {\it Ries et~al.} may be biased by
the fact that not all the non--gravitational perturbations
contributing to the observed residuals (e.g., Earth--reflected
radiation pressure, penumbra effects, variable part of charged
particle drag) depend on the satellite's rotational state, as
implicitly assumed by these authors; (iii) finally, the BI theory for
the spin axis evolution is itself based on a number of simplifying
assumptions and approximations, which might affect the results of the
fitting procedure.

\subsection{ Solution for LAGEOS~II}
In principle, an analysis of the same type as that carried out for
LAGEOS could be applied to its twin satellite LAGEOS~II, launched on
October 22, 1992, into an orbit similar to that of LAGEOS, but with an
inclination of $52\deg$ (instead of $110\deg$). Since the thermal
thrust effects depend upon the orientation of the orbit,
LAGEOS~II data may provide a valuable independent check for the
corresponding models. However, in
this case the short duration of the available data set prevents such a
detailed treatment for the time being.

Our best--fitting model for the LAGEOS~II along--track residuals is
plotted in Fig.~7.  The initial conditions for the spin axis direction
were adjusted to $\theta=169\deg \pm 2{\deg}.6$ and
$\lambda=260\deg \pm 15{\deg}.3$.  The longitude uncertainty is
larger for geometrical reasons, like in the case of LAGEOS.  These
values are not far from the polar coordinates ($\theta=173{\deg}.1$
and $\lambda=219{\deg}.1$) determined from both pre-- and post--launch
data by the flight control center [{\it L.~Anselmo}, personal
communication; see also \markcite{{\it Robbins et~al.,} 1994}].  Note
that in the LAGEOS~II case the launching phase data do not allow any
ambiguity in the spin axis direction (i.e., concerning the sense of
rotation).  According to our solution (see Fig.~8), during the
available data span of about $2.2$ years the spin axis has performed
only a small loop about the initial direction, under the influence of
the magnetic torque.  As a consequence, a model assuming a constant
direction would lead to equally good results.

The post--fit r.m.s. dispersions of the residuals is $0.71\times
10^{-12}~{\rm m/s}^2$, somewhat greater than in the LAGEOS case.  This
result --- together with the significant discrepancy in longitude with
respect to the launch data --- may point to some problem in the
dynamical model in this case, possibly due to a stronger contribution
by force mechanisms not included in it.  For instance, note that the
lower inclination and higher nodal rate of LAGEOS~II result into more
frequent and shorter eclipse periods, during which both the albedo and
the penumbra effects are most effective.  The variable component of
neutral and charged particle drag may also depend upon inclination.

As for the ``free'' parameters of the model, we have proceeded as in
the LAGEOS case.  For the constant negative acceleration due to
particle drag the value $-0.72\times 10^{-12}$ m/s$^2$ has been
selected.  The best--fitting value for the reflectivity difference
between the two LAGEOS~II hemispheres is about $0.023$, nearly twice
as great as for LAGEOS.  This is a puzzling finding because, according
to the technical staff of the Alenia Spazio Corporation which built
the satellite, none of the pre--launch tests had evidenced this
property [this information is in agreement with that reported by
\markcite{{\it Rubincam} [1990]}, based on measurements performed at
GSFC].  It has been recently pointed to us ({\it J.~Ries}, personal
communication) that a possible explanation may be related to some
darkening or blemishing of the satellite's surface in the final stage
of boosting.  A more detailed analysis of such effects is desirable in
the future.

For the rest, we have adopted the same value of the magnetization
parameter $\beta''(0)$ found by {\it Bertotti and Iess} from the
LAGEOS despinning data.  This assumption is reasonable, as the overall
structure of the two satellites is very similar.  On the other hand,
for LAGEOS~II accurate pre--launch measurements of the moments of
inertia are available, giving the value $\Delta=0.0379$ ({\it
L.~Anselmo}, personal communication).  There is an uncertainty of
a few percent in the initial rotation period, for which we
adopted the ``nominal'' value of $1$~s [\markcite{{\it Anselmo and
Pardini,} 1992}].  We have verified that the spin axis evolution is
not very sensitive to this parameter, provided the actual value is not
more than $50\%$ different from the assumed one.

To give an idea on the possible future evolution of the spin axis of
LAGEOS~II, we have plotted in Fig.~8 a numerical simulation covering a
time span of $15$ years starting at launch.  Of course, these results
provide only qualitative information, and more reliable predictions
will become possible once the initial conditions and other parameters
will be better constrained by the data.  However, Fig.~8 shows that
the spin axis evolution of LAGEOS~II will probably be be more
complicated than in the case of LAGEOS (see Fig.~4).  The main reason
for the difference is the higher nodal rate of LAGEOS~II and its
opposite sign with respect to LAGEOS.  Thus, for long--term,
high--accuracy orbit determination and prediction, reliable spin
evolution and thermal thrust models will be even more important for
LAGEOS~II than for LAGEOS.

\section{ Conclusions}
The main results obtained in this paper can be summarized as
follows:
\begin{enumerate}

\item The BI theoretical model for the evolution of LAGEOS'
rotational state, generalized here to account for a misalignment
between the satellite's spin axis and axis of symmetry, leads to a
successful fit of the observed along--track orbital residuals, once
the force mechanisms depending on the direction of the spin axis are
properly modelled.
\item The best--fit initial conditions for LAGEOS' spin axis
direction are $\theta=155\deg \pm 2{\deg}.5$ (colatitude) and
$\lambda=
133{\deg}.5 \pm 5{\deg}.7$ (longitude).  This Southward orientation of
the spin axis is consistent with both the observed despinning of the
satellite [\markcite{{\it Bertotti and Iess,} 1991}] and the
in--orbit injection information [\markcite{{\it Anselmo and Pardini,}
1992}], and must replace the
Northern orientation corresponding to the opposite sense of rotation,
which has been used in a number of recent studies.
\item The flattening parameter $\Delta=(C-A)/C$ for LAGEOS appears
to be close to the value indicated by pre--launch measurements, and
the tilt angle of the spin axis with respect to the axis of symmetry
is bound to be $<30\deg$.
\item In the forthcoming years, the rotation axis of LAGEOS will
undergo a rapid precession in longitude, approaching a limit cone of
$\approx 20\deg$ half--aperture around the South pole; at the same
time, the spin period will continue to increase in a
quasi--exponential fashion, until longitudinal temperature variations
across LAGEOS' surface will qualitatively change the thermal thrust
effects, and subsequently the $1:1$ resonance with the orbital period
will be reached.
\item The post--fit residuals are not randomly dis\-tri\-bu\-ted, but
include several harmonics related to the orientation of the orbital
plane with respect to the Sun, and therefore are probably mostly due
to radiation force mechanisms not included in our dynamical model
(e.g., Earth--reflected radiation pressure and penumbra effects).
\item The LAGEOS~II orbital data set also leads to a best--fit
solution with the BI theory corresponding to plausible initial
conditions and physical parameters for the satellite's properties.
However, the available data cover a too limited span of time to derive
from the best--fitting procedure reliable quantitative predictions on
the future evolution of the satellite's rotational state;
qualitatively, the spin axis of LAGEOS~II appears likely to follow a
more complex evolution pattern than in the case of LAGEOS.
\end{enumerate}

These results suggest some possible areas for future work.  As
remarked above, in about 10 years from now the ongoing magnetic
despinning of LAGEOS will force us to modify the basic physical
assumptions we used in this paper, by accounting for a more complex
behaviour of the thermal thrust forces and taking into account the
spin--orbit resonance in the rotational dynamics theory.  Although it
is possible that LAGEOS' rotation will undergo a chaotic behaviour
during resonance crossing, the small value of the $\Delta$ flattening
parameter  suggests that the ``thickness'' of the chaotic region in
the phase space may be limited. This issue deserves further studies,
based on non--averaged dynamical models such as that developed by
\markcite{{\it Habib et~al.} [1994]}.

The availability of a ``nominal'' model for the cou\-pled spin/orbital
evolution of LAGEOS provides us a more solid ground to compare with
observations the models for complex non--gravitational forces, such as
the albedo and penumbra effects and the variable part of neutral and
charged particle drag.  The ``nominal'' model can also be
ex\-ploit\-ed to re\-cons\-truct the per\-tur\-ba\-tions due to
ro\-ta\-tion--de\-pen\-dent force mechanisms on orbital elements other
than the semimajor axis.  Of course, all of this will depend on
further improvements in the knowledge of some free parameters
appearing in both the BI theory (including the initial conditions) and
the rotation--dependent force models, as well as on the use of high
accuracy tracking data from other satellites.  Further measurements of
LAGEOS' spin rate will also be extremely useful.

Finally, in the next decade a detailed comparison between the
rotational/orbital evolutions of LAGEOS and LAGEOS~II will become
possible.  A self--consistent model in the case of LAGEOS~II is
possibly more difficult than for LAGEOS, because of a more frequent
occurrence of the shadow--crossing periods and of the more complex
pattern of spin rate evolution.  On the other hand, the availability
of simultaneous data from two very similar satellites will probably
result into better observational constraints for the dynamical models.

\acknowledgements We are grateful to J. Ries (CSR, Univ. of Texas),
D.P. Rubincam and J. Robbins (NASA/ GSFC) for useful discussions and
for providing us with the time series of the LAGEOS along--track
residual accelerations derived at the Center for Space Research of the
University of Texas and at the NASA Goddard Space Flight Center,
respectively.  We also thank E.C. Pavlis and an anonymous referee for
a number of helpful remarks, L.~Anselmo, B.~Bertotti, L.~Iess and V.
Slabinski for many valuable suggestions and comments, and P.~Exertier
and G. M\'etris for help in testing our model.  P.F.
acknowledges support from the Italian Space Agency (ASI) and from the
Italian Ministry of University and Scientific Research (MURST).  D.V.
worked to this project while staying at OCA/CERGA, Grasse (France),
thanks to the ``H.~Poin\-ca\-r\'e'' research fellowship.

\newpage

\begin{references}
\reference
 Afonso, G., F. Barlier, C. Berger, F. Mignard, and J.J. Walch,
 Reassessment of the charge and neutral drag of LAGEOS and its
 geophysical implications, {\it J. Geo\-phys. Res., 90,} 9381-9398,
1985.

\reference
Afonso, G., F. Barlier, M. Carpino, et~al., Orbital effects of LAGEOS
 seasons and eclipses, {\it Annales Geophysic\ae,} {\it 7,} 501-514,
 1989.

\reference
Anselmo, L., and C. Pardini, Probing the Rubincam effect with
 LAGEOS--2, {\it Adv. Astr.  Sciences,} {\it 79,} 917-926, 1992.

\reference
Anselmo, L., P. Farinella, A. Milani, and A.M. Nobili, Effects of the
 Earth--reflected sunlight on the orbit of the LAGEOS satellite, {\it
 Astron.  Astrophys.,} {\it 117,} 3-8, 1983.

\reference
Barlier, F., M. Carpino, P. Farinella, et~al.,
Non-gra\-vi\-ta\-tion\-al
 perturbations on the semimajor axis of LAGEOS, {\it Annales
 Geophysic\ae,} {\it 4A,} 193-210, 1986.

\reference
Beletskii, V.V., {\it Motion of an Artificial Satellite about its
 Center of Mass}, Israel Program for Scientific Translation, Jerusalem,
 1966.

\reference
Beletskii, V.V., and A.A. Khentov, {\it Rotational Motion of a
 Magnetized Satellite}, Nauka, Moscow, 1985 (in Russian).

\reference
Beletskii, V.V., V.M. Grigorevsky, and S.Ya.  Ko\-les\-nik, The rotational
 motion of a satellite and atmospheric density determination, in {\it
 Satellite Dynamics}, G.E.O.  Giacaglia ed., pp. 233-279,
 Springer-Verlag, Berlin, 1975.

\reference
Bertotti, B., and P. Farinella, {\it Physics of the Earth and the
 Solar System}, Kluwer Academic Publ., Dordrecht, 1990.

\reference
Bertotti, B., and L. Iess, The rotation of LAGEOS, \jgr {\it , 96,}
 2431-2440, 1991.

\reference
Farinella, P., A.M. Nobili, F. Barlier, and F. Mignard, Effects of
 thermal thrust on the node and inclination of LAGEOS, {\it Astron.
 Astrophys.,} {\it 234,} 546-554, 1990.

\reference
Fitzpatrick, P.M., {\it Principles of Celestial Mechanics}, Academic
 Press, New York, 1970.

\reference
Habib, S., D.E. Holz, A. Kheyfetz, et~al., Spin dynamics of the LAGEOS
 satellite in support of a measurement of the Earth's gravitomagnetism,
 {\it Phys.~Rev., D50,} 6068-6079, 1994.

\reference
Johnson, C.W., C.A. Lundquist, and J.L. Zurasky, The LAGEOS satellite,
 paper presented at International Astronautical Federation XXVII
 Congress, Anaheim, California, Oct. 10-16, 1976.

\reference
Kolenkiewicz, R., J.W. Robbins, and M. Torrence, Non--conservative
 forces on LAGEOS I and II, {\it Adv. Space Res., 16,} 12(25)--12(28),
1995.

\reference
Landau, L.D., and E.M. Lifschitz, {\it Electrodynamics of Continuous
 Media}, Pergamon, New York, 1960.

\reference
Milani, A., A.M. Nobili, and P. Farinella, {\it
 Non-gravi\-ta\-ti\-onal perturbations and satellite geodesy},
 A. Hil\-ger, Bristol, 1987.

\reference
Peddie, N.W., and E.B. Fabiano, A proposed international geomagnetic
 reference field for 1965-1985, {\it J. Geomag.  Geoelectr.,} {\it
 34,} 357-364, 1982.

\reference
Ries, J.C., R.J. Eanes, and M.M. Watkins, Spin vector influence
 on LAGEOS ephemeris, presented at the Second Meeting of IAG
 Special Study Group 2.130, Baltimore, 1993.

\reference
Robbins, J.W., R.G. Williamson, and D.P. Rubincam, LAGEOS I and II
 spin axis evolution, {\it EOS Trans. AGU} {\it V.75 No.16,} 109,
 1994.

\reference
Rubincam, D.P, LAGEOS orbit decay due to infrared radiation from Earth,
 \jgr {\it , 92,} 1287-1294, 1987.

\reference
Rubincam, D.P., Yarkovsky thermal drag on LAGEOS, \jgr
 {\it , 93,} 13805-13810, 1988.

\reference
Rubincam, D.P., Drag on the LAGEOS satellite, \jgr {\it , 95,}
 4881-4886, 1990.

\reference
Scharroo, R., K.F. Wakker, B.A.C. Ambrosius, and R. Noo\-men, On the
 along-track acceleration of the LAGEOS satellite, \jgr
 {\it , 96,} 729-740, 1991.

\reference
Tapley, B.D., B.E. Schutz, R.J. Eanes, J.C. Ries, and M.M. Watkins,
 Lageos laser ranging contributions to geodynamics, geodesy, and
 orbital dynamics, in {\it Contributions of Space Geodesy to
 Geodynamics: Earth Dynamics}, AGU, pp. 147-173, 1993.

\reference
Vokrouhlick\'y, D., and P. Farinella, Specular reflection of sunlight
from wavy ocean surfaces and the albedo effect on satellite orbits,
II. LAGEOS long--term albedo perturbations revisited, {\it Astron.
Astrophys. 301}, 282-289, 1995.

\reference
Vokrouhlick\'y, D., P. Farinella, and D. Lucchesi, Long--perio\-dic
albedo perturbations on LAGEOS, {\it Astron.  Astrophys.,} {\it 280,}
282-294, 1993.

\reference
Vokrouhlick\'y, D., P. Farinella, and F. Mignard, Solar radiation
pressure perturbations for Earth satellites, II. An approxi\-mate
method to model penumbra transitions and their long--term orbital
effects on LAGEOS, {\it Astron.  Astrophys.,} {\it 285,} 333-334,
1994.

\end{references}

\newpage

\begin{figure}  % Fig. 1
\figurewidth{35pc}
\caption{ Time series of the along--track
 orbit--averaged residual accelerations of the LAGEOS satellite
 (starting time: May 29, 1976), as derived from laser
 tracking data at the Center of Space Research of the University of
 Texas, and compared with the predictions of the model for the
 spin--orbit evolution described in the text.  Empty and full squares
 correspond to the first three years and to later data, respectively
 (the former were considered with a $1/3$ statistical weight).  The
 solid curve shows the best--fit solution, obtained by adjusting the
 initial direction of LAGEOS' spin axis, the coefficient accounting for
 the optical asymmetry between the two hemispheres of the satellite and
 the thermal lag time of LAGEOS' retroreflectors.  Moreover, a constant
 term $-0.95\times 10^{-12}~{\rm m/s}^2$ was added to simulate the
 secular effect of charged and neutral particle drag.}
\end{figure}

\begin{figure}  % Fig. 2
\figurewidth{35pc}
\caption{ Post--fit residuals of the solution
 presented in Fig.~1.  The r.m.s. dispersion is $0.50\times
 10^{-12}~{\rm m/s}^2$.  Note the somewhat larger scatter of the data
 obtained during the first three years following launch (empty squares)
 and at the end of the treated time span.}
\end{figure}

\begin{figure} % Fig. 3
\figurewidth{35pc}
\caption{ Power spectrum of the post--fit
 residuals shown in Fig.~2.  Spectral density is plotted vs. frequency
 (in day$^{-1}$).  The main lines correspond to the frequencies
 $m({\dot \lambda}_\odot - {\dot \Omega})$, where $m=2,3,4$ and ${\dot
 \lambda}_\odot - {\dot \Omega}$ is the difference between the rate of
 the solar longitude and of the nodal longitude of the satellite's
 orbit. These frequencies are typically associated to Earth reflected
 and rediffused radiation pressure perturbations [see \markcite{{\it
 Barlier et~al.,} 1986}; \markcite{{\it Vokrouhlick\'y et~al.,} 1993}].
No $m=1$ line show up at the position indicated by an arrow, probably
because the corresponding harmonic has been ``absorbed'' in the fit
with the force model described in the text.} \end{figure}

\begin{figure}  % Fig. 4
\figurewidth{20pc}
\caption{ The LAGEOS spin axis evolution
 corresponding to the along--track perturbations shown by the solid
 curve in Fig.~1.  Here the horizonthal components $(s_x, s_y)$ of the
 spin axis unit vector in an equatorial reference frame are plotted at
 different times.  The initial conditions leading to the best fit with
 the orbital data (corresponding to polar coordinates
 $\theta=155\deg$ and $\lambda=133{\deg}.5$) were adopted here.
 The solid curve corresponds to the first $17.9$ years after
 launch
 (the time span covered by the data shown in Fig.~1), while the dashed
 curve shows the future evolution of LAGEOS' spin axis direction
 inferred from the BI model.  Square marks along the curve correspond
 to 1--year intervals.  The limiting circle of radius $\approx 0.34$ is
 reached after about $25$ years since launch.  The averaging
 assumptions implicit in the BI model are no longer valid after about
 $33$ years, when the 1:1 spin--orbit resonance occurs.}
\end{figure}

\begin{figure}  % Fig. 5
\figurewidth{20pc}
\caption{ Contour lines of constant r.m.s.
 dispersion for LAGEOS' residual along--track accelerations from the
 fit with the BI theory, in the plane of the initial colatitude and
 longitude of the spin axis direction.  The fit was carried out
 by keeping the other free parameters of the model fixed at
 their ``nominal values''. We have checked that adjusting at the same
 time the values of the amplitude of the
 Yarkovsky-Schach thermal effect, the thermal lag time $\tau_{RR}$ and
 the reflectivity asymmetry parameters leads to very similar
 results.}
\end{figure}

\begin{figure}   % Fig. 6
\figurewidth{20pc}
\caption{ A comparison between the LAGEOS
 spin axis evolution obtained by fitting the residual along--track
 accelerations with the BI model (the same as in Fig.~4; solid curve 1),
 that derived by \markcite{{\it Ries et~al.} [1993]} (dashed curve 3),
 and the best match to the Ries et~al. curve with the BI model
 (short--dashed
 curve 2). The latter method yields initial polar coordinates
 $\theta=162\deg$ and $\lambda=152\deg$.}
\end{figure}

\begin{figure} % Fig. 7
\figurewidth{35pc}
\caption{ Fit of the observed along--track
 residual accelerations for LAGEOS~II (the horizontal axis gixes
 the time elapsed since launch on October 22, 1992).  The dynamical
 model includes thermal
 thrust effects with the same values for the corresponding parameters
 as in the LAGEOS solution, while a reflectivity asymmetry of $0.023$
 was selected.  The constant secular term due to charged and neutral
 particle drag is assumed to be $-0.72\times 10^{-12}~{\rm m/s}^2$.
 The r.m.s. dispersion of the post--fit residuals is $0.71\times
 10^{-12}~{\rm m/s}^2$.}
\end{figure}

\begin{figure}
\figurewidth{20pc}
\caption{ Probable scenario for the evolution
 of the LAGEOS~II spin axis direction over a time span of $15$ years
 since launch.  The same magnetization parameters as in the LAGEOS
 case have been adopted.  The initial conditions were adjusted as
 discussed in the text, while the value $\Delta = 0.0379$ for the
 flattening parameter was derived by pre--launch
 measurements.  Note that at the beginning the spin axis latitude
 drifts away from the pole, opposite to the LAGEOS case.  The faster
 nodal rate of the orbit results into a more complicated evolution
 pattern with respect to LAGEOS.}
\end{figure}

\end{document}
